diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m index 6a49bbcac2d153a1935fc92352fd3c225faf25b8..647be1a95f627d57b58eff2bd70c0517df537cfc 100644 --- a/matlab/+occbin/kalman_update_engine.m +++ b/matlab/+occbin/kalman_update_engine.m @@ -32,7 +32,6 @@ if nargin>26 is_multivariate = false; end -use_relaxation = false; if is_multivariate [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options); @@ -109,23 +108,79 @@ if info==0 oldstart = regimes_(1).regimestart(end); newstart = regx(1).regimestart(end); diffstart = newstart-oldstart; + regname = 'regimestart'; 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); + [diffstart, diffregime] = max([newstart1-oldstart1,newstart2-oldstart2]); + switch diffregime + case 1 + regname = 'regimestart1'; + case 2 + regname = 'regimestart2'; end end +end +if options_.occbin.filter.use_relaxation && diffstart>options_.occbin.filter.use_relaxation + guess_regime = [base_regime base_regime]; + options_.occbin.filter.guess_regime = true; + guess_regime(1) = regx(1); + regx2= regx; + while isequal(guess_regime(1),regx2(1)) + % we reduce length until the converged regime does not change -if options_.occbin.filter.use_relaxation && diffstart>2 - if info0==0 - % make sure we match criteria to enter further solution attempts - info1=1; + guess_regime(1).(regname)(end) = regx2(1).(regname)(end)-1; + if guess_regime(1).(regname)(end-1)==guess_regime(1).(regname)(end) + guess_regime(1).(regname)(end-1) = guess_regime(1).(regname)(end-1)-1; + end + if is_multivariate + [ax2, a1x2, Px2, P1x2, vx2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, etahat2, alphahat2, V2] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,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); + else + [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 + 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; + Px=Px2; + P1x=P1x2; + vx=vx2; + Tx=Tx2; + Rx=Rx2; + Cx=Cx2; + regx=regx2; + info=info2; + likx=likx2; + M_= M_2; + etahat=etahat2; + alphahat=alphahat2; + if is_multivariate + V=V2; + else + Fix = Fix2; + Kix = Kix2; + TTx = TTx2; + RRx = RRx2; + CCx = CCx2; + end + end end - options_.occbin.likelihood.brute_force_regime_guess = true; - options_.occbin.likelihood.loss_function_regime_guess = true; - use_relaxation = true; + options_.occbin.filter.guess_regime = false; + options_.occbin.likelihood.brute_force_regime_guess = false; + options_.occbin.likelihood.loss_function_regime_guess = false; end if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime)) @@ -145,9 +200,8 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| ( end if info2==0 use_index= 1; - if not(info==0 && isequal(regx2{1},regx)) && not(use_relaxation && likx2{1}>=likx) - % found a solution, different from previous or - % use_relaxation and likelihood is better + if not(info==0 && isequal(regx2{1},regx)) + % found a solution, different from previous break end end @@ -161,8 +215,6 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| ( end if info2==0 use_index = 2; - % if use_relaxation and we are here, previous guess did not - % improve solution, so we test for this one end if use_index @@ -224,18 +276,16 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| ( end if info2==0 use_index= gindex; - if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx) + if not(info==0 && isequal(regx2{gindex},regx)) % found a solution, different from previous one - % use_relaxation and likelihood improves break end end end % loop over other regime slack, binding in expectation or binding in current period if info2==0 - if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx) + if not(info==0 && isequal(regx2{gindex},regx)) % found a solution, different from previous one - % use_relaxation and likelihood improves break end end @@ -243,9 +293,8 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| ( end % loop over current regime binding in expectation vs binding in current period if info2==0 - if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx) + if not(info==0 && isequal(regx2{gindex},regx)) % found a solution, different from previous one - % use_relaxation and likelihood improves break end end diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m index b720160c5392e6bcbba83ba096921b5141cf4d5d..04e13d54df14d98f8a7b783848481936f19cddeb 100644 --- a/matlab/+occbin/set_default_options.m +++ b/matlab/+occbin/set_default_options.m @@ -46,6 +46,7 @@ if ismember(flag,{'filter','all'}) options_occbin_.filter.guess_regime = false; options_occbin_.filter.periodic_solution = true; options_occbin_.filter.use_relaxation = false; + options_occbin_.filter.use_relaxation_tol_period = 1; end if ismember(flag,{'forecast','all'})