diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m
index fe2346436ad94fa82febe265bb9bbd70ac83fcad..6a49bbcac2d153a1935fc92352fd3c225faf25b8 100644
--- a/matlab/+occbin/kalman_update_engine.m
+++ b/matlab/+occbin/kalman_update_engine.m
@@ -1,6 +1,6 @@
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