diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 6fa4a27fe13d33781d1da1a7d268e41f5a47d6d9..364b328862aaf3b86735f10fcbe5c26117419457 100644
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@@ -124,6 +124,8 @@ occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_bin
occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;
error_indicator=false;
+options_.noprint = true;
+
try
%blanket try-catch should be replaced be proper error handling, see https://git.dynare.org/Dynare/dynare/-/merge_requests/2226#note_20318
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
@@ -178,7 +180,6 @@ opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
opts_simul.periods = size(opts_simul.SHOCKS,1);
options_.occbin.simul=opts_simul;
-options_.noprint = true;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if out.error_flag
fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
diff --git a/matlab/+occbin/IVF_core.m b/matlab/+occbin/IVF_core.m
index 30666cd3220ee8bcd99793c54a9f9442514e4f45..f3d4c4650798a14d2e9dcc542b158fec5f60401d 100644
--- a/matlab/+occbin/IVF_core.m
+++ b/matlab/+occbin/IVF_core.m
@@ -1,4 +1,4 @@
-function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
+function [filtered_errs, resids, Emat, stateval, error_code, regime_history] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
% Calls the solver to get the shocks explaining the data for the inversion filter
%
@@ -92,6 +92,9 @@ for this_period=1:sample_length
filtered_errs=NaN;
error_code(1) = 304;
error_code(4) = 1000;
+ if this_period == 1
+ regime_history(this_period) = [];
+ end
if options_.occbin.likelihood.waitbar; dyn_waitbar_close(hh_fig); end
return
end
@@ -99,9 +102,10 @@ for this_period=1:sample_length
opts_simul.SHOCKS = err_vals_out;
- [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_] = occbin.match_function(...
+ [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_, out] = occbin.match_function(...
err_vals_out,obs_list,current_obs,opts_simul, M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_);
init_val = stateval(this_period,:); %update
+ regime_history(this_period) = out.regime_history(1);
if max(abs(err_vals_out))>1e8
error_code(1) = 306;
error_code(4) = max(abs(err_vals_out))/1000;
diff --git a/matlab/+occbin/IVF_posterior.m b/matlab/+occbin/IVF_posterior.m
index 69c6c110ddb34815d31c8df8701ae1e24c5fe29a..dbfcaa9c8eac2c42f8f8f96918f59228017e1780 100644
--- a/matlab/+occbin/IVF_posterior.m
+++ b/matlab/+occbin/IVF_posterior.m
@@ -1,4 +1,4 @@
-function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
+function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov, regime_history] = IVF_posterior(xparam1,...
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
% dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
@@ -58,6 +58,8 @@ trend_coeff = [];
obs = dataset_.data;
obs_list = options_.varobs(:);
exit_flag = 1;
+% initialize output argument in case of early return
+regime_history = [];
if size(xparam1,1)<size(xparam1,2)
@@ -107,7 +109,7 @@ end
sample_length = size(obs,1);
filtered_errs_init = zeros(sample_length,sum(err_index));
-[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
+[filtered_errs, resids, Emat, stateval, info, regime_history] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
if info(1)
fval = Inf;
exit_flag = 0;
diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m
index 11205e63e0c80159a2af1994e51d4ca4a5a5acd6..83cbd587ba8f3ae9e440bd46ce7e1743d54a5f1e 100644
--- a/matlab/+occbin/kalman_update_algo_1.m
+++ b/matlab/+occbin/kalman_update_algo_1.m
@@ -213,7 +213,12 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
- [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);
+ [a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = 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);
+ if error_flag
+ etahat=NaN(size(QQQ,1),1);
+ warning(orig_warning_state);
+ return;
+ end
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
@@ -278,7 +283,12 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
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);
+ [a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = 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);
+ if error_flag
+ etahat=NaN(size(QQQ,1),1);
+ warning(orig_warning_state);
+ return;
+ end
end
else
error_flag = 330;
@@ -294,7 +304,7 @@ end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1))
- error_flag = 1;
+ error_flag = 331;
end
if ~error_flag
@@ -357,7 +367,7 @@ else
F = ZZ*P(:,:,t)*ZZ' + H(di,di);
sig=sqrt(diag(F));
if any(any(isnan(F)))
- error_flag=1;
+ error_flag=325;
warning(orig_warning_state);
return;
end
diff --git a/matlab/+occbin/kalman_update_algo_3.m b/matlab/+occbin/kalman_update_algo_3.m
index 457c47480666f57757ceeca79d35c83fe6883b5e..efffc9badbb50593472ba8681a4f0d2133733efc 100644
--- a/matlab/+occbin/kalman_update_algo_3.m
+++ b/matlab/+occbin/kalman_update_algo_3.m
@@ -282,7 +282,7 @@ end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1)) %fixed point algorithm did not converge
- error_flag = 1;
+ error_flag = 331;
end
if ~error_flag
diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m
index fe2346436ad94fa82febe265bb9bbd70ac83fcad..3459a26712cdb504bbc9852ee227ee886914a9fd 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)
@@ -32,13 +32,14 @@ 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
[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 +49,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 +59,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 +79,7 @@ else
etahat=etahat1;
alphahat=alphahat1;
if is_multivariate
- V=V1;
+ V=V1;
else
Fix = Fix1;
Kix = Kix1;
@@ -85,8 +92,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,27 +104,83 @@ 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);
+ if M_.occbin.constraint_nbr==1
+ 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, 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))
@@ -137,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
@@ -153,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
@@ -216,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
+ 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
@@ -235,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
@@ -261,6 +318,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 +370,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 +409,17 @@ 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
+
+if info(1)==0
+ if isnan(likx)
+ info = 323;
+ elseif any(any(isnan(ax)))
+ info = 324;
+ end
+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'})
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 52a51a4cf5862903152be09b763afef77b5c07d7..a1a01039280d585d8009bea5f34d4244804ef580 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -189,7 +189,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
if options_.order==1 && ~options_.particle.status
if options_.smoother
if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
- [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov, oo_.occbin.smoother.regime_history] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ismember(info(1),[303,304,306])
fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
else
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 9680a675737846095a2609abd8d036494dc67189..7156e3c253c4e7d2dc1b60b35d92491e4d016227 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -197,6 +197,16 @@ switch info(1)
message = 'Occbin: there was a problem in running the smoother. Simulation within smoother failed.';
case 322
message = 'Occbin: smoother did not converge.';
+ case 323
+ message = 'Piecewise linear Kalman filter: the likelihood is NaN.';
+ case 324
+ message = 'Piecewise linear Kalman filter: updated state vector is NaN.';
+ case 325
+ message = 'Piecewise linear Kalman filter: filter covariance NaN.';
+ case 330
+ message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (periodic loop).';
+ case 331
+ message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (max number of iterations reached).';
case 401
message = 'Cycle reduction reached the iteration limit. Try increasing maxit.';
case 402
diff --git a/matlab/kalman/DsgeSmoother.m b/matlab/kalman/DsgeSmoother.m
index f4fb1f801f1a922b92b86e9e25cfbf3b934c3cf4..c84b800286e301ff5a65c02db894c4c4c5f723e9 100644
--- a/matlab/kalman/DsgeSmoother.m
+++ b/matlab/kalman/DsgeSmoother.m
@@ -506,6 +506,13 @@ else
aaa = zeros(nk,M_.endo_nbr,gend+nk);
aaa(:,oo_.dr.restrict_var_list,:)=aK;
+ if isoccbin
+ tstart = 1;
+ else
+ % we enter here in the first occbin smoother iteration
+ % occbin kalman update is not yet able to accommodate diffuse steps!
+ tstart=d+2;
+ end
for k=2:gend+1
opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
opts_simul.waitbar = options_.occbin.smoother.waitbar;
@@ -526,8 +533,24 @@ else
% period ahead (so if regimestart was [1 5] it should be [1 4]
% in out
% end
- for jnk=1:nk
- aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.piecewise(jnk,:) - out.ys';
+ if out.error_flag==0
+ for jnk=1:nk
+ aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.piecewise(jnk,:) - out.ys';
+ end
+ elseif k>tstart
+ % the issue only matters non-stationary models, with
+ % diffuse filter, and for the first occbin smoother iteration,
+ % where tstart>1
+ %
+ % if k>tstart, the same simulation should have been done
+ % already in occbin.kalman_update, so it should never give
+ % an error
+ %
+ % if k<=tstart, the simulation may crash, since we ignore OBC in the first (diffuse) steps
+ % and it may happen that, given the linear updated states,
+ % the occbin simulation does not converge
+ error('this error should not occur, please contact the developers!')
+
end
end
aK=aaa;
diff --git a/matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m
index 2b6f9046681e1916fd64a62465cf26df5539391c..366e415f4a8a51b728427a1d363dfde2abaa0cf7 100644
--- a/matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m
@@ -544,7 +544,7 @@ while notsteady && t<smpl
opts_simul.init_regime = []; %regimes_(t);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul;
- [~, out] = occbin.solver(M_,options_,dr,steady_state,exo_steady_state,exo_det_steady_state);
+ [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
for jnk=1:nk
if filter_covariance_flag
diff --git a/matlab/kalman/save_display_classical_smoother_results.m b/matlab/kalman/save_display_classical_smoother_results.m
index 48dcad0b5b31a0e7eae91bbe043acf70331dac18..f7f2f981e710722c51b9f8ff466d92773f743cac 100644
--- a/matlab/kalman/save_display_classical_smoother_results.m
+++ b/matlab/kalman/save_display_classical_smoother_results.m
@@ -36,7 +36,7 @@ gend= dataset_.nobs;
n_varobs = length(options_.varobs);
if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
- [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+ [~, info, ~, ~, ~, ~, ~, ~, ~, ~, oo_.dr, atT, innov, oo_.occbin.smoother.regime_history] = occbin.IVF_posterior(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ismember(info(1),[303,304,306])
fprintf('\nIVF: smoother did not succeed. No results will be written to oo_.\n')
else
diff --git a/matlab/shock_decomposition/annualized_shock_decomposition.m b/matlab/shock_decomposition/annualized_shock_decomposition.m
index 962a0910c13b6c261da19af17dc824114a6f0588..834e37c9f3acc964661d8ca0a1dca11e8008239e 100644
--- a/matlab/shock_decomposition/annualized_shock_decomposition.m
+++ b/matlab/shock_decomposition/annualized_shock_decomposition.m
@@ -200,6 +200,7 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
% make annualized shock decomp
z = annualiz(z,t0,q2a,aux,steady_state);
+ oo_.annualized_realtime_shock_decomposition.(['yr_' int2str(yr)]) = z;
end
oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr)]) = z(:,:,end-nfrcst:end);
diff --git a/tests/occbin/filter/NKM.mod b/tests/occbin/filter/NKM.mod
index 1b552290f4dce06214ae6430c3b87fdea665d224..a54142f94ad7f1c928af1183ac1c8b716bd12ba1 100644
--- a/tests/occbin/filter/NKM.mod
+++ b/tests/occbin/filter/NKM.mod
@@ -309,7 +309,7 @@ varobs yg inom pi;
datafile=dataobsfile2, mode_file=NKM_mh_mode_saved,
mode_compute=0, nobs=120, first_obs=1,
mh_replic=0, plot_priors=0, smoother,
- consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty);
+ consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1:8],smoothed_state_uncertainty);
// plot regimes
occbin.plot_regimes(oo_.occbin.smoother.regime_history,M_,options_)
@@ -350,7 +350,7 @@ varobs yg inom pi;
datafile=dataobsfile2, mode_file=NKM_mh_mode_saved,
mode_compute=0, nobs=120, first_obs=1,
mh_replic=0, plot_priors=0, smoother, smoother_redux,
- consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty);
+ consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1:8],smoothed_state_uncertainty);
// check consistency of smoother_redux
for k=1:M_.endo_nbr,
@@ -362,6 +362,33 @@ varobs yg inom pi;
disp('smoother redux successfully recovers full smoother results!')
end
+ for k=1:M_.endo_nbr,
+ mer(k)=max(abs(oo_.UpdatedVariables.(M_.endo_names{k})-oo0.UpdatedVariables.(M_.endo_names{k})));
+ end
+ if max(mer)>1.e-10
+ error('smoother redux does not recover full updated variables results!')
+ else
+ disp('smoother redux successfully recovers full updated variables results!')
+ end
+
+
+ for k=1:M_.endo_nbr,
+ mer(k)=max(abs(oo_.FilteredVariables.(M_.endo_names{k})-oo0.FilteredVariables.(M_.endo_names{k})));
+ end
+ if max(mer)>1.e-10
+ error('smoother redux does not recover full filtered variables results!')
+ else
+ disp('smoother redux successfully recovers full filtered variables results!')
+ end
+
+ for k=1:M_.endo_nbr,
+ mer(k)=max(max(max(abs(oo_.FilteredVariablesKStepAhead(:,k,:)-oo0.FilteredVariablesKStepAhead(:,k,:)))));
+ end
+ if max(mer)>1.e-10
+ error('smoother redux does not recover full k-step ahead variables results!')
+ else
+ disp('smoother redux successfully recovers full k-step ahead variables results!')
+ end
// use inversion filter (note that IF provides smoother together with likelihood)
occbin_setup(likelihood_inversion_filter,smoother_inversion_filter);
@@ -369,7 +396,7 @@ varobs yg inom pi;
datafile=dataobsfile2, mode_file=NKM_mh_mode_saved,
mode_compute=0, nobs=120, first_obs=1,
mh_replic=0, plot_priors=0, smoother,
- consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty);
+ consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1:8],smoothed_state_uncertainty);
// show initial condition effect of IF
figure('Name','OccBin: Smoothed shocks')