Skip to content
Snippets Groups Projects
Commit 17451344 authored by Marco Ratto's avatar Marco Ratto
Browse files

- when multiple solutions are found, compute likelihood as the sum of all solutions found

parent 7041e107
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/+occbin/kalman_update_engine.m
+ 53
− 17
View file @ 17451344
function [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = ...
kalman_update_engine(a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,...
dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
kalman_update_engine(a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,...
dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
% [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = kalman_update_engine(
% a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,
% dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
......@@ -39,6 +39,8 @@ if is_multivariate
else
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
likvec = likx;
regvec = regx(1);
info0=info;
if info
if ~isequal(regimes_(1:2),[base_regime base_regime])
......@@ -48,6 +50,8 @@ if info
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
end
likvec = likx;
regvec = regx(1);
info1=info;
else
if ~isequal(regimes_(1:2),[base_regime base_regime])
......@@ -56,6 +60,10 @@ else
else
[ax1, a1x1, Px1, P1x1, vx1, Fix1, Kix1, Tx1, Rx1, Cx1, regx1, info1, M_1, likx1, alphahat1, etahat1,TTx1,RRx1,CCx1] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info1==0 && not(isequal(regx1,regx))
likvec = [likvec likx1];
regvec = [regvec; regx1(1)];
end
if info1==0 && likx1<likx
ax=ax1;
a1x=a1x1;
......@@ -72,7 +80,7 @@ else
etahat=etahat1;
alphahat=alphahat1;
if is_multivariate
V=V1;
V=V1;
else
Fix = Fix1;
Kix = Kix1;
......@@ -85,8 +93,8 @@ else
if t>options_.occbin.likelihood.number_of_initial_periods_with_extra_regime_guess
info1=0;
else
% may help in first 2 periods to try some other guess regime, due to
% larger state uncertainty
% may help in first 2 periods to try some other guess regime, due to
% larger state uncertainty
info1=1;
options_.occbin.likelihood.brute_force_regime_guess = true;
options_.occbin.likelihood.loss_function_regime_guess = true;
......@@ -97,17 +105,17 @@ end
diffstart=0;
if info==0
if M_.occbin.constraint_nbr==1
oldstart = regimes_(1).regimestart(end);
newstart = regx(1).regimestart(end);
diffstart = newstart-oldstart;
else
newstart1 = regx(1).regimestart1(end);
newstart2 = regx(1).regimestart2(end);
oldstart1 = regimes_(1).regimestart1(end);
oldstart2 = regimes_(1).regimestart2(end);
diffstart = max(newstart1-oldstart1,newstart2-oldstart2);
end
if M_.occbin.constraint_nbr==1
oldstart = regimes_(1).regimestart(end);
newstart = regx(1).regimestart(end);
diffstart = newstart-oldstart;
else
newstart1 = regx(1).regimestart1(end);
newstart2 = regx(1).regimestart2(end);
oldstart1 = regimes_(1).regimestart1(end);
oldstart2 = regimes_(1).regimestart2(end);
diffstart = max(newstart1-oldstart1,newstart2-oldstart2);
end
end
if options_.occbin.filter.use_relaxation && diffstart>2
......@@ -219,7 +227,7 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
% found a solution, different from previous one
% use_relaxation and likelihood improves
break
break
end
end
end % loop over other regime slack, binding in expectation or binding in current period
......@@ -261,6 +269,18 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
% so that we DO NOT enter IVF step
info0=0;
info1=0;
isnew=true;
for kr=1:length(likvec)
% make sure likelihood does not differ by rounding issue
% but truly for different regimes
if isequal(regx2{use_index}(1),regvec(kr))
isnew = false;
end
end
if isnew
likvec = [likvec likx2{use_index}];
regvec = [regvec; regx2{use_index}(1)];
end
end
if info2==0 && likx2{use_index}<likx
ax=ax2{use_index};
......@@ -301,6 +321,18 @@ if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %||
[ax2, a1x2, Px2, P1x2, vx2, Fix2, Kix2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, alphahat2, etahat2,TTx2,RRx2,CCx2] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
options_.occbin.filter.guess_regime = false;
isnew=true;
for kr=1:length(likvec)
% make sure likelihood does not differ by rounding issue
% but due to different regimes
if isequal(regx2(1),regvec(kr))
isnew = false;
end
end
if isnew
likvec = [likvec likx2];
regvec = [regvec; regx2(1)];
end
if info2==0 && likx2<likx
ax=ax2;
a1x=a1x2;
......@@ -328,4 +360,8 @@ if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %||
end
end
end
if length(likvec)>1
% sum the likelihood of multiple solutions
likx = -2*log(sum(exp(-likvec./2)));
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment