diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m index bf25201a0..42705e12d 100644 --- a/matlab/+occbin/kalman_update_algo_1.m +++ b/matlab/+occbin/kalman_update_algo_1.m @@ -105,17 +105,17 @@ PZI = P1(:,:,t)*ZZ'*iF(di,di,t); % L(:,:,t) = T-K(:,di,t)*ZZ; L(:,:,t) = eye(mm)-PZI*ZZ; -if ~options_.occbin.filter.use_relaxation && ~isempty(fieldnames(regimes0)) +if ~isempty(fieldnames(regimes0)) if options_.occbin.filter.guess_regime [~,~,~,~,~,~, TTx, RRx, CCx] ... = occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, regimes0(1),'reduced_state_space',T0,R0); - if M_.occbin.constraint_nbr==1 - bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart); - bindx = bindx(2:end); - [regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx); - bindx = bindx(2:end); - [regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx); - else + if M_.occbin.constraint_nbr==1 + bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart); + bindx = bindx(2:end); + [regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx); + bindx = bindx(2:end); + [regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx); + else bindx1 = occbin.backward_map_regime(regimes0(1).regime1, regimes0(1).regimestart1); bindx2 = occbin.backward_map_regime(regimes0(1).regime2, regimes0(1).regimestart2); bindx1 = bindx1(2:end); @@ -126,7 +126,7 @@ if ~options_.occbin.filter.use_relaxation && ~isempty(fieldnames(regimes0)) bindx2 = bindx2(2:end); [regimes0(3).regime1, regimes0(3).regimestart1, error_flag]=occbin.map_regime(bindx1); [regimes0(3).regime2, regimes0(3).regimestart2, error_flag]=occbin.map_regime(bindx2); - end + end % regimes0=[regimes0 base_regime base_regime]; TT(:,:,2) = TTx(:,:,end); RR(:,:,2) = RRx(:,:,end); @@ -200,48 +200,13 @@ end lik_hist=lik; niter=1; is_periodic=0; -if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime) - nguess=1; -else - nguess=0; -end -newguess=0; if any(myregime) || ~isequal(regimes_(1),regimes0(1)) while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations niter=niter+1; - oldstart=1; - if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1 - oldstart = regimes0(1).regimestart(end); - end - newstart=1; - if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1 - newstart = regimes_(1).regimestart(end); - end - if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation - regimestart = max(oldstart+2,round(0.5*(newstart+oldstart))); - regimestart = min(regimestart,oldstart+4); - if regimestart<=regimes_(1).regimestart(end-1) - if length(regimes_(1).regimestart)<=3 - regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart)); - else - regimes_(1).regime = regimes_(1).regime(1:end-2); - regimes_(1).regimestart = regimes_(1).regimestart(1:end-2); - regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1); - end - end - regimes_(1).regimestart(end)=regimestart; - [~,~,~,~,~,~, TTx, RRx, CCx] ... - = occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space', T0, R0); - TT(:,:,2) = TTx(:,:,end); - RR(:,:,2) = RRx(:,:,end); - CC(:,2) = CCx(:,end); - elseif newguess==0 - TT(:,:,2)=ss.T(my_order_var,my_order_var,1); - RR(:,:,2)=ss.R(my_order_var,:,1); - CC(:,2)=ss.C(my_order_var,1); - end - newguess=0; + TT(:,:,2)=ss.T(my_order_var,my_order_var,1); + RR(:,:,2)=ss.R(my_order_var,:,1); + CC(:,2)=ss.C(my_order_var,1); regime_hist(niter) = {regimes_(1)}; if M_.occbin.constraint_nbr==1 regime_end(niter) = regimes_(1).regimestart(end); @@ -257,9 +222,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) else opts_simul.endo_init = alphahat(dr.inv_order_var,1); end - if not(options_.occbin.filter.use_relaxation) - opts_simul.init_regime=regimes_(1); - end + opts_simul.init_regime=regimes_(1); if M_.occbin.constraint_nbr==1 myregimestart = [regimes_.regimestart]; else @@ -287,11 +250,11 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) is_periodic_iter = find(is_periodic); is_periodic = any(is_periodic); if is_periodic - % re-set to previous regime - if options_.occbin.filter.periodic_solution - % force projection conditional on most likely regime - [m, im]=min(lik_hist(is_periodic_iter:end)); - opts_simul.init_regime=regime_hist{is_periodic_iter+im-1}; + % re-set to previous regime + if options_.occbin.filter.periodic_solution + % force projection conditional on most likely regime + [m, im]=min(lik_hist(is_periodic_iter:end)); + opts_simul.init_regime=regime_hist{is_periodic_iter+im-1}; if M_.occbin.constraint_nbr==1 myregimestart = [regimes0.regimestart]; else @@ -303,23 +266,22 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) [~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state); if out.error_flag error_flag = out.error_flag; - etahat=etahat(:,2); - lik=inf; - return; - else - regimes_ = out.regime_history; - TT(:,:,2)=ss.T(my_order_var,my_order_var,1); - RR(:,:,2)=ss.R(my_order_var,:,1); - CC(:,2)=ss.C(my_order_var,1); - [a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance); - end - else - error_flag = 330; etahat=etahat(:,2); lik=inf; return; + else + regimes_ = out.regime_history; + TT(:,:,2)=ss.T(my_order_var,my_order_var,1); + RR(:,:,2)=ss.R(my_order_var,:,1); + CC(:,2)=ss.C(my_order_var,1); + [a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance); end -% end + else + error_flag = 330; + etahat=etahat(:,2); + lik=inf; + return; + end end end end diff --git a/matlab/+occbin/kalman_update_algo_3.m b/matlab/+occbin/kalman_update_algo_3.m index 797bd62ac..77c4a8cfc 100644 --- a/matlab/+occbin/kalman_update_algo_3.m +++ b/matlab/+occbin/kalman_update_algo_3.m @@ -105,7 +105,7 @@ if options_.smoother_redux end mm=size(a,1); -if ~options_.occbin.filter.use_relaxation && ~isempty(fieldnames(regimes0)) +if ~isempty(fieldnames(regimes0)) if options_.occbin.filter.guess_regime [~,~,~,~,~,~, TTx, RRx, CCx] ... = occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, regimes0(1),myrestrict,T0,R0); @@ -198,52 +198,13 @@ end lik_hist=lik; niter=1; is_periodic=0; -if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime) - nguess=1; -else - nguess=0; -end -newguess=0; if any(myregime) || ~isequal(regimes_(1),regimes0(1)) while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations niter=niter+1; - newstart=1; - if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1 - newstart = regimes_(1).regimestart(end); - end - oldstart=1; - if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1 - oldstart = regimes0(1).regimestart(end); - end - if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation - regimestart = max(oldstart+2,round(0.5*(newstart+oldstart))); - regimestart = min(regimestart,oldstart+4); - if regimestart<=regimes_(1).regimestart(end-1) - if length(regimes_(1).regimestart)<=3 - regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart)); - else - regimes_(1).regime = regimes_(1).regime(1:end-2); - regimes_(1).regimestart = regimes_(1).regimestart(1:end-2); - regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1); - end - end - - % % if (newstart-oldstart)>3 - % % regimestart = regimes_(1).regimestart(end-1)+oldstart+2; - % % % regimestart = regimes_(1).regimestart(end-1)+round(0.5*(newstart+oldstart))-1; - regimes_(1).regimestart(end)=regimestart; - [~,~,~,~,~,~, TTx, RRx, CCx] ... - = occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0); - TT(:,:,2) = TTx(:,:,end); - RR(:,:,2) = RRx(:,:,end); - CC(:,2) = CCx(:,end); - elseif newguess==0 - TT(:,:,2)=ss.T(my_order_var,my_order_var,1); - RR(:,:,2)=ss.R(my_order_var,:,1); - CC(:,2)=ss.C(my_order_var,1); - end - newguess=0; + TT(:,:,2)=ss.T(my_order_var,my_order_var,1); + RR(:,:,2)=ss.R(my_order_var,:,1); + CC(:,2)=ss.C(my_order_var,1); regime_hist(niter) = {regimes_(1)}; if M_.occbin.constraint_nbr==1 regime_end(niter) = regimes_(1).regimestart(end); @@ -251,11 +212,6 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) [a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol); lik_hist(niter) = lik; opts_simul.SHOCKS(1,:) = etahat(:,2)'; - % if opts_simul.restrict_state_space - % tmp=zeros(M_.endo_nbr,1); - % tmp(dr.restrict_var_list,1)=alphahat(:,1); - % opts_simul.endo_init = tmp(dr.inv_order_var,1); - % else if opts_simul.restrict_state_space tmp=zeros(M_.endo_nbr,1); tmp(dr.restrict_var_list,1)=alphahat(:,1); @@ -263,10 +219,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) else opts_simul.endo_init = alphahat(dr.inv_order_var,1); end - % end - if not(options_.occbin.filter.use_relaxation) - opts_simul.init_regime=regimes_(1); - end + opts_simul.init_regime=regimes_(1); if M_.occbin.constraint_nbr==1 myregimestart = [regimes_.regimestart]; else @@ -293,7 +246,6 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1)) % re-set to previous regime if options_.occbin.filter.periodic_solution % force projection conditional on most likely regime -% [m, im]=max(lik_hist(is_periodic_iter:end)); [m, im]=min(lik_hist(is_periodic_iter:end)); opts_simul.init_regime=regime_hist{is_periodic_iter+im-1}; if M_.occbin.constraint_nbr==1 diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m index dfc7350cc..fe2346436 100644 --- a/matlab/+occbin/kalman_update_engine.m +++ b/matlab/+occbin/kalman_update_engine.m @@ -32,6 +32,8 @@ 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); else @@ -93,6 +95,31 @@ else end 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 +end + +if options_.occbin.filter.use_relaxation && diffstart>2 + if info0==0 + % make sure we match criteria to enter further solution attempts + info1=1; + end + options_.occbin.likelihood.brute_force_regime_guess = true; + options_.occbin.likelihood.loss_function_regime_guess = true; + use_relaxation = true; +end + if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime)) guess_regime = [base_regime base_regime]; @@ -110,8 +137,9 @@ 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)) - % found a solution, different from previous one + 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 break end end @@ -125,6 +153,8 @@ 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 @@ -186,16 +216,18 @@ 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)) + 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 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)) + 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 end end @@ -203,8 +235,9 @@ 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)) + 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 end end diff --git a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc index d006a785d..cb08e0773 100644 --- a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc +++ b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc @@ -89,6 +89,9 @@ occbin_graph(noconstant) c erra lambdak k i a k; occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200); occbin_setup(smoother_periods=200,smoother_maxit=200,smoother_curb_retrench,smoother_check_ahead_periods=200); calib_smoother(datafile=datasim); + oo0=oo_; + occbin_setup(filter_use_relaxation); + calib_smoother(datafile=datasim); oo_= occbin.unpack_simulations(M_,oo_,options_); titlelist = char('c','lambdak','k','i','a','erra');