From b90dff1986b9af57cfe3db29d74e72cbf18e53e5 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 11:09:39 +0100 Subject: [PATCH 01/21] retrieve states in period 0 from smoother and deploy this to: - allow having binding regimes in period 1 with occbin - further improve recovering of smoothed variables in period 1 under smoother_redux option --- matlab/+occbin/DSGE_smoother.m | 35 ++--- matlab/DsgeSmoother.m | 31 +++-- matlab/missing_DiffuseKalmanSmootherH1_Z.m | 73 +++++++++-- matlab/missing_DiffuseKalmanSmootherH3_Z.m | 144 +++++++++++++++++---- 4 files changed, 230 insertions(+), 53 deletions(-) diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m index 36ab1f0c9..3ed80498c 100644 --- a/matlab/+occbin/DSGE_smoother.m +++ b/matlab/+occbin/DSGE_smoother.m @@ -1,4 +1,4 @@ -function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info) +function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info) %function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info) % Runs a DSGE smoother with occasionally binding constraints % @@ -39,6 +39,8 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P % - oo_ [structure] storing the results % - options_ [structure] describing the options % - bayestopt_ [structure] describing the priors +% - alphahat0 [double] (m*1) array, smoothed endogenous variables in period 0 (a_{0|T}) (decision-rule order) +% - state_uncertainty0 [double] (K,K,1) array, storing the uncertainty in period 0 % Copyright © 2021 Dynare Team % @@ -118,7 +120,7 @@ occbin_options.first_period_occbin_update = options_.occbin.smoother.first_perio occbin_options.opts_regime = opts_simul; % this builds the opts_simul options field needed by occbin.solver occbin_options.opts_regime.binding_indicator = options_.occbin.likelihood.init_binding_indicator; occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history; -[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT; +[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT; oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history; regime_history = oo_.occbin.smoother.regime_history; @@ -135,9 +137,9 @@ opts_regime.binding_indicator=[]; regime_history0 = regime_history; fprintf('Occbin smoother iteration 1.\n') -opts_simul.SHOCKS = [etahat(:,2:end)'; zeros(1,M_.exo_nbr)]; +opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)]; opts_simul.exo_pos = 1:M_.exo_nbr; -opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1); +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! options_.occbin.simul=opts_simul; options_.noprint = true; @@ -145,15 +147,18 @@ options_.noprint = true; regime_history = out.regime_history; if options_.smoother_redux occbin_options.opts_simul.restrict_state_space =1; - oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1); - oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1); + [T0,R0] = dynare_resolve(M_,options_,oo_); + oo_.occbin.linear_smoother.T0=T0; + oo_.occbin.linear_smoother.R0=R0; + % oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1); + % oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1); end TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:); RR = ss.R(oo_.dr.order_var,:,:); CC = ss.C(oo_.dr.order_var,:); -TT = cat(3,TT(:,:,1),TT); -RR = cat(3,RR(:,:,1),RR); -CC = cat(2,CC(:,1),CC); +% TT = cat(3,TT(:,:,1),TT); +% RR = cat(3,RR(:,:,1),RR); +% CC = cat(2,CC(:,1),CC); opts_regime.regime_history = regime_history; opts_regime.binding_indicator = []; @@ -184,7 +189,7 @@ while is_changed && maxiter>iter && ~is_periodic iter=iter+1; fprintf('Occbin smoother iteration %u.\n', iter) occbin_options.opts_regime.regime_history=regime_history; - [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT; + [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT; sto_etahat(iter)={etahat}; regime_history0(iter,:) = regime_history; if occbin_smoother_debug @@ -195,17 +200,17 @@ while is_changed && maxiter>iter && ~is_periodic sto_RR = RR; sto_TT = TT; - opts_simul.SHOCKS = [etahat(:,2:end)'; zeros(1,M_.exo_nbr)]; - opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1); + opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)]; + opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1); options_.occbin.simul=opts_simul; [~, out, ss] = occbin.solver(M_,oo_,options_); regime_history = out.regime_history; TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:); RR = ss.R(oo_.dr.order_var,:,:); CC = ss.C(oo_.dr.order_var,:); - TT = cat(3,TT(:,:,1),TT); - RR = cat(3,RR(:,:,1),RR); - CC = cat(2,CC(:,1),CC); +% TT = cat(3,TT(:,:,1),TT); +% RR = cat(3,RR(:,:,1),RR); +% CC = cat(2,CC(:,1),CC); opts_regime.regime_history = regime_history; [TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_); diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 7ea0cd462..e184f34ca 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -1,4 +1,4 @@ -function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin) +function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin) % Estimation of the smoothed variables and innovations. % % INPUTS @@ -265,7 +265,7 @@ if kalman_algo == 1 || kalman_algo == 3 a_initial = zeros(np,1); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=T*a_initial; %set state prediction for first Kalman step; - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ... Z,R1,Q,H,Pinf,Pstar, ... data1,vobs,np,smpl,data_index, ... options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty,options_.filter_covariance,options_.smoother_redux); @@ -321,7 +321,7 @@ if kalman_algo == 2 || kalman_algo == 4 a_initial = zeros(np,1); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=ST*a_initial; %set state prediction for first Kalman step; - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ... Z,R1,Q,diag(H), ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ... options_.nk,kalman_tol,diffuse_kalman_tol, ... @@ -538,17 +538,20 @@ else static_var_list0(static_var_list) = ilagged; static_var_list(static_var_list) = ~ilagged; % reconstruct smoothed variables - aaa=zeros(M_.endo_nbr,gend); - aaa(oo_.dr.restrict_var_list,:)=alphahat; + aaa=zeros(M_.endo_nbr,gend+1); + aaa(oo_.dr.restrict_var_list,1)=alphahat0; + aaa(oo_.dr.restrict_var_list,2:end)=alphahat; for k=1:gend - aaa(static_var_list,k) = C(~ilagged,:)*alphahat(:,k)+D(~ilagged,:)*etahat(:,k); + aaa(static_var_list,k+1) = C(~ilagged,:)*alphahat(:,k)+D(~ilagged,:)*etahat(:,k); end if any(ilagged) for k=2:gend - aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k); + aaa(static_var_list0,k+1) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k); end + aaa(static_var_list0,2) = Tstar(ilagged,:)*alphahat0+Rstar(ilagged,:)*etahat(:,1); end - alphahat=aaa; + alphahat0=aaa(:,1); + alphahat=aaa(:,2:end); % reconstruct updated variables aaa=zeros(M_.endo_nbr,gend); @@ -619,6 +622,18 @@ else state_uncertainty=sstate_uncertainty; clear sstate_uncertainty end + if ~isempty(state_uncertainty0) + mm=size(T,1); + ss=length(find(static_var_list)); + sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr); + sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list)=state_uncertainty0(1:mm,1:mm); + sstate_uncertainty(static_var_list,static_var_list)=[C(~ilagged,:) D(~ilagged,:)]*state_uncertainty0*[C(~ilagged,:) D(~ilagged,:)]'; + tmp = [C(~ilagged,:) D(~ilagged,:)]*state_uncertainty0; + sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list)=tmp(1:ss,1:mm); + sstate_uncertainty(oo_.dr.restrict_var_list,static_var_list)=transpose(sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list)); + state_uncertainty0=sstate_uncertainty; + clear sstate_uncertainty + end % reconstruct PK if ~isempty(PK) diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index 8e0fad922..09d567263 100644 --- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m @@ -1,4 +1,4 @@ -function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) +function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d,alphahat0,V0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) % function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix. @@ -140,8 +140,21 @@ if state_uncertainty_flag else V=[]; end +alphahat0=[]; +V0=[]; t = 0; +if rank(Pinf(:,:,1),diffuse_kalman_tol) + % this is needed to get smoothed states in period 0 with diffuse steps + % i.e. period 0 is a filter step without observables + Pinf_init = Pinf(:,:,1); + Pstar_init = Pstar(:,:,1); + a_init = a(:,1); + a(:,t) = T*a(:,t); + % only non-stationary part is affected by following line, + % hence Pstar on EXIT from diffuse step will NOT change. + Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ; +end while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && td+1 end end -if d %diffuse periods - % initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23 +if d==0 % get smoother in period t=0 + a0 = a(:,1); + P0 = P(:,:,1); + L0=T; + r0 = L0'*r(:,1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0 + alphahat0 = a0 + P0*r0; %DK (2012), eq. 4.35 + if state_uncertainty_flag + N0=L0'*N(:,:,1)*L0; %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0 + if smoother_redux + ptmp = [P0 R*Q; (R*Q)' Q]; + ntmp = [N0 zeros(mm,rr); zeros(rr,mm+rr)]; + V0 = ptmp - ptmp*ntmp*ptmp; + else + V0 = P0-P0*N0*P0; %DK (2012), eq. 4.43 + end + end + +else %diffuse periods + % initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23 r0 = zeros(mm,d+1); r0(:,d+1) = r(:,d+1); %set r0_{d}, i.e. shifted by one period r1 = zeros(mm,d+1); %set r1_{d}, i.e. shifted by one period @@ -357,7 +387,7 @@ if d %diffuse periods if ~Finf_singular(1,t) r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ... - + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1) + + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1) if state_uncertainty_flag L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1) N(:,:,t)=Linf(:,:,t)'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.19, noting that L^(0) is named Linf @@ -374,7 +404,7 @@ if d %diffuse periods r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14) if state_uncertainty_flag N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)... - +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) + +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) N_1(:,:,t)=T'*N_1(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14) N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14) end @@ -395,8 +425,8 @@ if d %diffuse periods V(:,:,t) = pstmp - pstmp*ntmp*pstmp... -(pitmp*ntmp1*pstmp)'... - pitmp*ntmp1*pstmp... - - pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30 - + - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30 + else V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)... -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... @@ -405,6 +435,33 @@ if d %diffuse periods end end end + % compute states and covarinace in period t=0 + r10 = T'*r1(:,1); + r00 = T'*r0(:,1); + alphahat0 = a_init + Pstar_init*r00 + Pinf_init*r10; % DK (2012), eq. 5.23 + if state_uncertainty_flag + N0=T'*N(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf + N_10=T'*N_1(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf + N_20=T'*N_2(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf + if smoother_redux + pstmp = [Pstar_init R*Q; (R*Q)' Q]; + pitmp = [Pinf_init zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp = [N0 zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp1 = [N_10 zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp2 = [N_20 zeros(mm,rr); zeros(rr,mm+rr)]; + V0 = pstmp - pstmp*ntmp*pstmp... + -(pitmp*ntmp1*pstmp)'... + - pitmp*ntmp1*pstmp... + - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30 + + else + V0 = Pstar_init-Pstar_init*N0*Pstar_init... + -(Pinf_init*N_10*Pstar_init)'... + - Pinf_init*N_10*Pstar_init... + - Pinf_init*N_20*Pinf_init; % DK (2012), eq. 5.30 + end + end + end if decomp_flag diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m index 8f100f528..b8fcd416c 100644 --- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m @@ -1,4 +1,4 @@ -function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d,varargout] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_) +function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d,alphahat0,V0,varargout] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_) % function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_) % Computes the diffuse kalman smoother in the case of a singular var-cov matrix. % Univariate treatment of multivariate time series. @@ -149,6 +149,8 @@ if state_uncertainty_flag else V=[]; end +alphahat0=[]; +V0=[]; if ~occbin_.status isoccbin = 0; @@ -180,7 +182,7 @@ else T0=occbin_.info{5}; R0=occbin_.info{6}; else - + TT=occbin_.info{5}; RR=occbin_.info{6}; CC=occbin_.info{7}; @@ -201,7 +203,7 @@ else CC=repmat(zeros(mm,1),1,smpl+1); end end - + else TT=repmat(T,1,1,smpl+1); RR=repmat(R,1,1,smpl+1); @@ -210,7 +212,7 @@ else if ~smoother_redux T0=T; R0=R; - + end if ~isinf(occbin_options.first_period_occbin_update) % initialize state matrices (otherwise they are set to 0 for @@ -219,7 +221,7 @@ else RRR=repmat(R0,1,1,smpl+1); CCC=repmat(zeros(length(T0),1),1,smpl+1); end - + end t = 0; @@ -229,6 +231,14 @@ if ~isempty(Pinf(:,:,1)) else newRank = rank(Pinf(:,:,1),diffuse_kalman_tol); end +if newRank + % add this to get smoothed states in period 0 + Pinf_init = Pinf(:,:,1); + Pstar_init = Pstar(:,:,1); + Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ; + ainit = a1(:,1); +end + while newRank && t < smpl t = t+1; a(:,t) = a1(:,t); @@ -270,7 +280,7 @@ while newRank && t < smpl else oldRank = 0; end - if isoccbin, + if isoccbin TT(:,:,t+1)= T; RR(:,:,t+1)= R; end @@ -298,7 +308,11 @@ while newRank && t < smpl end if isoccbin - first_period_occbin_update = max(t+2,occbin_options.first_period_occbin_update); + first_period_occbin_update = occbin_options.first_period_occbin_update; + if d>0 + first_period_occbin_update = max(t+2,occbin_options.first_period_occbin_update); + % kalman update is not yet robust to accommodate diffuse steps + end if occbin_options.opts_regime.waitbar hh = dyn_waitbar(0,'Occbin: Piecewise Kalman Filter'); set(hh,'Name','Occbin: Piecewise Kalman Filter.'); @@ -322,6 +336,10 @@ Pinf1 = Pinf1(:,:,1:d); notsteady = 1; while notsteady && t d+1 Ni = T'*Ni*T; % KD (2000), eq. (23), equation for N_{t-1,p_{t-1}} end end -if d + +if d==0 % recover states in period t=0 + a0 = ainit; + r0 = ri; + P0 = Pinit; + % if OCCBIN, P1 in t=1 must be consistent with the regime in 1 + alphahat0 = a0 + P0*r0; + + % we do NOT need eps(0)! + % alphahat is smoothed state in t=0, so that S(1)=T*s(0)+R*eps(1); + if state_uncertainty_flag + N0 = Ni; % DK (2012), below 6.15, N_{t-1}=N_{t,0} + if smoother_redux + ptmp = [P0 R*Q; (R*Q)' Q]; + ntmp = [N0 zeros(mm,rr); zeros(rr,mm+rr)]; + V0 = ptmp - ptmp*ntmp*ptmp; + else + V0 = P0-P0*N0*P0; % KD (2000), eq. (7) with N_{t-1} stored in N(:,:,t) + end + end + +else % diffuse filter r0 = zeros(mm,d); r0(:,d) = ri; r1 = zeros(mm,d); @@ -641,8 +714,8 @@ if d V(:,:,t) = pstmp - pstmp*ntmp0*pstmp... -(pitmp*ntmp1*pstmp)'... - pitmp*ntmp1*pstmp... - - pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30 - + - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30 + else V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)... -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... @@ -654,14 +727,41 @@ if d r0(:,t-1) = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} r1(:,t-1) = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} if state_uncertainty_flag - N_0(:,:,t-1)= T'*N_0(:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T - N_1(:,:,t-1)= T'*N_1(:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T - N_2(:,:,t-1)= T'*N_2(:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T + N_0(:,:,t-1)= T'*N_0(:,:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T + N_1(:,:,t-1)= T'*N_1(:,:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T + N_2(:,:,t-1)= T'*N_2(:,:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T + end + else + r00 = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} + r10 = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0} + if state_uncertainty_flag + N_00= T'*N_0(:,:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T + N_10= T'*N_1(:,:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T + N_20= T'*N_2(:,:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T end end end -else - alphahat0 = 0*a1(:,1) + P1(:,:,1)*ri; + % get smoothed states in t=0 + alphahat0 = ainit + Pstar_init*r00 + Pinf_init*r10; % KD (2000), eq. (26) + if state_uncertainty_flag + if smoother_redux + pstmp = [Pstar_init R*Q; (R*Q)' Q]; + pitmp = [Pinf_init zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp0 = [N_00 zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp1 = [N_10 zeros(mm,rr); zeros(rr,mm+rr)]; + ntmp2 = [N_20 zeros(mm,rr); zeros(rr,mm+rr)]; + V0 = pstmp - pstmp*ntmp0*pstmp... + -(pitmp*ntmp1*pstmp)'... + - pitmp*ntmp1*pstmp... + - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30 + + else + V0=Pstar_init-Pstar_init*N_00*Pstar_init... + -(Pinf_init*N_10*Pstar_init)'... + - Pinf_init*N_10*Pstar_init... + - Pinf_init*N_20*Pinf_init; % DK (2012), eq. 5.30 + end + end end if decomp_flag @@ -675,7 +775,7 @@ if decomp_flag ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)'; end end - + % calculate eta_tm1t if isoccbin if isqvec From 0622f01f4efacc39de3cc1715e42780b571f5334 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 11:56:42 +0100 Subject: [PATCH 02/21] filter out t=1 when storing aalphahat and eetahat --- matlab/missing_DiffuseKalmanSmootherH3_Z.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m index b8fcd416c..179b28377 100644 --- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m @@ -387,7 +387,7 @@ while notsteady && t1 aalphahat(:,t-1) = aha(:,1); eetahat(:,t) = etaha(:,2); end From 21c9d59e8ceabcaf3b80c4591e930032ebf35d0b Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 11:57:41 +0100 Subject: [PATCH 03/21] incorporate information about states in period 0 for occbin with smoother_redux --- matlab/DsgeSmoother.m | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index e184f34ca..59b3b4c6e 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -391,8 +391,9 @@ else opts_simul = options_.occbin.simul; end % reconstruct smoothed variables - aaa=zeros(M_.endo_nbr,gend); - aaa(oo_.dr.restrict_var_list,:)=alphahat; + aaa=zeros(M_.endo_nbr,gend+1); + aaa(oo_.dr.restrict_var_list,1)=alphahat0; + aaa(oo_.dr.restrict_var_list,2:end)=alphahat; iTx = zeros(size(TTx)); for k=1:gend if isoccbin @@ -416,13 +417,19 @@ else static_var_list0 = static_var_list; static_var_list0(static_var_list) = ilagged; static_var_list(static_var_list) = ~ilagged; - aaa(static_var_list,k) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged); - if any(ilagged) && k>1 - aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged); + aaa(static_var_list,k+1) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged); + if any(ilagged) + if k>1 + aaa(static_var_list0,k+1) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged); + else + aaa(static_var_list0,2) = Tstar(ilagged,:)*alphahat0+Rstar(ilagged,:)*etahat(:,1)+Cstar(ilagged); + end + end end - alphahat=aaa; + alphahat0=aaa(:,1); + alphahat=aaa(:,2:end); % reconstruct updated variables bbb=zeros(M_.endo_nbr,gend); @@ -488,6 +495,13 @@ else state_uncertainty=sstate_uncertainty; clear sstate_uncertainty end + if ~isempty(state_uncertainty0) + mm=size(T,1); + sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr); + sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list)=state_uncertainty0(1:mm,1:mm); + state_uncertainty0=sstate_uncertainty; + clear sstate_uncertainty + end aaa = zeros(nk,M_.endo_nbr,gend+nk); aaa(:,oo_.dr.restrict_var_list,:)=aK; From 21443e0988228cc77811e180ae6e9f0cf418e927 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 11:59:27 +0100 Subject: [PATCH 04/21] add consistency check of results using smoother_redux with occbin --- tests/occbin/filter/NKM.mod | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/occbin/filter/NKM.mod b/tests/occbin/filter/NKM.mod index 8c04a966e..81b25d1d7 100644 --- a/tests/occbin/filter/NKM.mod +++ b/tests/occbin/filter/NKM.mod @@ -319,6 +319,16 @@ varobs yg inom pi; mh_replic=0, plot_priors=0, smoother, smoother_redux, nodisplay,consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty); + // check consistency of smoother_redux + for k=1:M_.endo_nbr, + mer(k)=max(abs(oo_.SmoothedVariables.(M_.endo_names{k})-oo0.SmoothedVariables.(M_.endo_names{k}))); + end + if max(mer)>1.e-10 + error('smoother redux does not recover full smoother results!') + else + disp('smoother redux successfully recovers full smoother results!') + end + // use inversion filter (note that IF provides smoother together with likelihood) occbin_setup(likelihood_inversion_filter,smoother_inversion_filter); From 74acb2d1f52ca75a171e512d907ed555d1e527ab Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 12:44:37 +0100 Subject: [PATCH 05/21] bug fix using index t --- matlab/missing_DiffuseKalmanSmootherH1_Z.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index 09d567263..3cbc664d6 100644 --- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m @@ -150,7 +150,7 @@ if rank(Pinf(:,:,1),diffuse_kalman_tol) Pinf_init = Pinf(:,:,1); Pstar_init = Pstar(:,:,1); a_init = a(:,1); - a(:,t) = T*a(:,t); + a(:,1) = T*a(:,1); % only non-stationary part is affected by following line, % hence Pstar on EXIT from diffuse step will NOT change. Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ; From a91c36f920bf029949e84f6369888697de8e9136 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 17:33:38 +0100 Subject: [PATCH 06/21] first_period_occbin_update is an index, not a boolean --- matlab/+occbin/set_default_options.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m index 17e74a323..34f31b47f 100644 --- a/matlab/+occbin/set_default_options.m +++ b/matlab/+occbin/set_default_options.m @@ -67,7 +67,7 @@ end if ismember(flag,{'likelihood','all'}) options_occbin_.likelihood.curb_retrench = false; - options_occbin_.likelihood.first_period_occbin_update = true; + options_occbin_.likelihood.first_period_occbin_update = 1; options_occbin_.likelihood.full_output = false; options_occbin_.likelihood.IF_likelihood = false; options_occbin_.likelihood.init_regime_history = []; @@ -187,7 +187,7 @@ if ismember(flag,{'smoother','all'}) options_occbin_.smoother.curb_retrench = false; options_occbin_.smoother.debug = false; options_occbin_.smoother.fast = false; - options_occbin_.smoother.first_period_occbin_update = true; + options_occbin_.smoother.first_period_occbin_update = 1; options_occbin_.smoother.full_output = false; % options.occbin.smoother.init_mode = 1; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period options_occbin_.smoother.init_regime_history = []; From 6598d615d19fcb60cf4fba71551452f139227ead Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 17:34:55 +0100 Subject: [PATCH 07/21] in some cases, order is not yet set to 1 --- matlab/+occbin/solver.m | 3 +++ 1 file changed, 3 insertions(+) diff --git a/matlab/+occbin/solver.m b/matlab/+occbin/solver.m index 9dd44f88f..3dc465ba3 100644 --- a/matlab/+occbin/solver.m +++ b/matlab/+occbin/solver.m @@ -54,6 +54,9 @@ if solve_dr if isempty(options_.qz_criterium) options_.qz_criterium = 1+1e-6; end + if options_.order>1 + options_.order = 1; + end [dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_); out.error_flag=error_flag; From a346639020bf751350a794ea5a5e83977cd1aa42 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 18 Nov 2022 17:36:29 +0100 Subject: [PATCH 08/21] when filter is NOT diffuse, allow likelihood to handle violation of constraints in period 1. --- .../missing_observations_kalman_filter.m | 126 +++++++++++------- 1 file changed, 77 insertions(+), 49 deletions(-) diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m index ba5145fad..65024aada 100644 --- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m +++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m @@ -105,7 +105,11 @@ if occbin_.status occbin_options=occbin_.info{4}; opts_regime.regime_history = occbin_options.opts_simul.init_regime; opts_regime.binding_indicator = occbin_options.opts_simul.init_binding_indicator; - first_period_occbin_update = max(t+1,options_.occbin.likelihood.first_period_occbin_update); + if t>1 + first_period_occbin_update = max(t+1,options_.occbin.likelihood.first_period_occbin_update); + else + first_period_occbin_update = options_.occbin.likelihood.first_period_occbin_update; + end if isempty(opts_regime.binding_indicator) && isempty(opts_regime.regime_history) opts_regime.binding_indicator=zeros(last+2,M_.occbin.constraint_nbr); end @@ -123,8 +127,8 @@ if occbin_.status TT=repmat(T,1,1,last+1); RR=repmat(R,1,1,last+1); CC=repmat(zeros(mm,1),1,last+1); - end - + end + end else first_period_occbin_update = inf; @@ -142,6 +146,10 @@ while notsteady && t<=last QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. end end + if t==1 + Pinit = P1(:,:,1); + ainit = a1(:,1); + end s = t-start+1; d_index = data_index{t}; if isqvec @@ -178,61 +186,81 @@ while notsteady && t<=last % badly_conditioned_F = true; end end - if ~occbin_.status || (occbin_.status && (options_.occbin.likelihood.use_updated_regime==0 || t=no_more_missing_observations && ~isqvec && ~occbin_.status - notsteady = max(abs(K(:)-oldK))>riccati_tol; - oldK = K(:); + + P = T*(P-K*z*P)*transpose(T)+QQ; + else + K = P(:,z)*iF; + if occbin_.status + P0(:,:,t) = (P-K*P(z,:)); end + P = T*(P-K*P(z,:))*transpose(T)+QQ; + end + if occbin_.status + a0(:,t) = (a+K*v); + vv(d_index,t) = v; + end + a = T*(a+K*v)+C; + if t>=no_more_missing_observations && ~isqvec && ~occbin_.status + notsteady = max(abs(K(:)-oldK))>riccati_tol; + oldK = K(:); end end end end if occbin_.status && t>=first_period_occbin_update - - if isqvec - Qt = Qvec(:,:,t-1:t+1); - end + occbin_options.opts_simul.waitbar=0; - [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,oo_,options_,occbin_options); + if t==1 + if isqvec + Qt = cat(3,Q,Qvec(:,:,t:t+1)); + end + a00 = ainit; + a10 = [a00 a(:,1)]; + P00 = Pinit; + P10 = P1(:,:,[1 1]); + data_index0{1}=[]; + data_index0(2)=data_index(1); + v0(:,2)=vv(:,1); + Y0(:,2)=Y(:,1); + Y0(:,1)=nan; + TT01 = cat(3,T,TT(:,:,1)); + RR01 = cat(3,R,RR(:,:,1)); + CC01 = zeros(size(CC,1),2); + CC01(:,2) = CC(:,1); + [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a00, a10, P00, P10, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), M_, oo_, options_, occbin_options); + else + if isqvec + Qt = Qvec(:,:,t-1:t+1); + end + [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,oo_,options_,occbin_options); + end if info if options_.debug fprintf('\nmissing_observations_kalman_filter:PKF failed with: %s\n', get_error_message(info,options_)); @@ -252,7 +280,7 @@ while notsteady && t<=last P0(:,:,t) = Px(:,:,1); P1(:,:,t) = P1x(:,:,2); P = Px(:,:,2); - + end t = t+1; end From a4dc58db8ee12d5a2b9027aa299a4503b9e16521 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 2 Dec 2022 09:16:39 +0100 Subject: [PATCH 09/21] adjust smoother redux test, to consider that now ALL smoothed variables must be properly reconstructed also in period 1 --- .../fs2000_smoother_redux.mod | 28 +++---------------- 1 file changed, 4 insertions(+), 24 deletions(-) diff --git a/tests/kalman_filter_smoother/fs2000_smoother_redux.mod b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod index 72d261a59..057dd82ae 100644 --- a/tests/kalman_filter_smoother/fs2000_smoother_redux.mod +++ b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod @@ -112,12 +112,7 @@ for k=1:length(vlist1) merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo1.FilteredVariables.(vlist1{k}))); end if max(merr1)>1.e-12 - for k=1:length(vlist1) - merr2(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})(2:end)-oo1.SmoothedVariables.(vlist1{k})(2:end))); - end - if max(merr2)>1.e-12 - error('smoother_redux with kalman_algo=1 does not replicate original smoothed static variables!') - end + error('smoother_redux with kalman_algo=1 does not replicate original smoothed static variables!') end if max(merr1U)>1.e-12 for k=1:length(vlist1) @@ -128,12 +123,7 @@ if max(merr1U)>1.e-12 end end if max(merr1F)>1.e-12 - for k=1:length(vlist1) - merr2F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})(2:end)-oo1.FilteredVariables.(vlist1{k})(2:end))); - end - if max(merr2F)>1.e-12 - error('smoother_redux with kalman_algo=1 does not replicate original filtered static variables!') - end + error('smoother_redux with kalman_algo=1 does not replicate original filtered static variables!') end merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo1.FilteredVariablesKStepAhead)))); if max(merrK)>1.e-12 @@ -184,12 +174,7 @@ for k=1:length(vlist1) merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo2.FilteredVariables.(vlist1{k}))); end if max(merr1)>1.e-12 - for k=1:length(vlist1) - merr2(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})(2:end)-oo2.SmoothedVariables.(vlist1{k})(2:end))); - end - if max(merr2)>1.e-12 - error('smoother_redux with kalman_algo=2 does not replicate original smoothed static variables!') - end + error('smoother_redux with kalman_algo=2 does not replicate original smoothed static variables!') end if max(merr1U)>1.e-12 for k=1:length(vlist1) @@ -200,12 +185,7 @@ if max(merr1U)>1.e-12 end end if max(merr1F)>1.e-12 - for k=1:length(vlist1) - merr2F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})(2:end)-oo2.FilteredVariables.(vlist1{k})(2:end))); - end - if max(merr2F)>1.e-12 - error('smoother_redux with kalman_algo=2 does not replicate original filtered static variables!') - end + error('smoother_redux with kalman_algo=2 does not replicate original filtered static variables!') end merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo2.FilteredVariablesKStepAhead)))); if max(merrK)>1.e-12 From e51a30eaf5cb73c59792aac005187f567272a84c Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 2 Dec 2022 09:17:38 +0100 Subject: [PATCH 10/21] Bug fix for the case where no prior is defined for estim params --- matlab/evaluate_likelihood.m | 55 ++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 21 deletions(-) diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m index 652a258cc..315663333 100644 --- a/matlab/evaluate_likelihood.m +++ b/matlab/evaluate_likelihood.m @@ -47,24 +47,24 @@ end if ischar(parameters) switch parameters - case 'posterior mode' - parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_); - case 'posterior mean' - parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); - case 'posterior median' - parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_); - case 'prior mode' - parameters = bayestopt_.p5(:); - case 'prior mean' - parameters = bayestopt_.p1; - otherwise - disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:') - disp(' ''posterior mode'', ') - disp(' ''posterior mean'', ') - disp(' ''posterior median'', ') - disp(' ''prior mode'' or') - disp(' ''prior mean''.') - error + case 'posterior mode' + parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_); + case 'posterior mean' + parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); + case 'posterior median' + parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_); + case 'prior mode' + parameters = bayestopt_.p5(:); + case 'prior mean' + parameters = bayestopt_.p1; + otherwise + disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:') + disp(' ''posterior mode'', ') + disp(' ''posterior mean'', ') + disp(' ''posterior median'', ') + disp(' ''prior mode'' or') + disp(' ''prior mean''.') + error end end @@ -73,10 +73,23 @@ if isempty(dataset) end options_=select_qz_criterium_value(options_); +[~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); +if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0) + % Plot prior densities. + % Set prior bounds + bounds = prior_bounds(bayestopt_, options_.prior_trunc); +else % estimated parameters but no declared priors + % No priors are declared so Dynare will estimate the model by + % maximum likelihood with inequality constraints for the parameters. + bounds.lb = lb; + bounds.ub = ub; +end + + if options_.occbin.likelihood.status && options_.occbin.likelihood.inversion_filter - llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); -else - llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); + llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); +else + llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); end ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_); llik = llik - ldens; \ No newline at end of file From 50aa20c742595f46d1b5a5612dbb98d4803e00d9 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 2 Dec 2022 09:23:22 +0100 Subject: [PATCH 11/21] eliminate buggy and useless check for existence of guess regime variable. this may create bugs in running PKF update step! --- matlab/+occbin/solve_one_constraint.m | 2 +- matlab/+occbin/solve_two_constraints.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m index f5f0c1446..a55ebd0e0 100644 --- a/matlab/+occbin/solve_one_constraint.m +++ b/matlab/+occbin/solve_one_constraint.m @@ -105,7 +105,7 @@ end SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods); SS_out.C = nan(DM.n_vars,n_shocks_periods); -if ~exist('regime_history_','var') || isempty(regime_history_guess) +if isempty(regime_history_guess) regime_history = struct(); guess_history = false; else diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m index 80744ac50..9869f6339 100644 --- a/matlab/+occbin/solve_two_constraints.m +++ b/matlab/+occbin/solve_two_constraints.m @@ -115,7 +115,7 @@ end SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods); SS_out.C = NaN(DM.n_vars,n_shocks_periods); -if ~exist('regime_history_','var') || isempty(regime_history_guess) +if isempty(regime_history_guess) regime_history = struct(); guess_history = false; else From 7132cb593e17eee4c6783f8a6bf60fb7f7da94e0 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 2 Dec 2022 11:11:31 +0100 Subject: [PATCH 12/21] bug fix for init variables to be defined only with occbin PKF --- .../likelihood/missing_observations_kalman_filter.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m index 65024aada..e966ffaa3 100644 --- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m +++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m @@ -145,10 +145,10 @@ while notsteady && t<=last if ~(isqvec) QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. end - end - if t==1 - Pinit = P1(:,:,1); - ainit = a1(:,1); + if t==1 + Pinit = P1(:,:,1); + ainit = a1(:,1); + end end s = t-start+1; d_index = data_index{t}; From 22e6fb181492d52c8586cc0ec69248d8c3197ac3 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 2 Dec 2022 12:05:00 +0100 Subject: [PATCH 13/21] make sure also ALL updated variables in period 1 are reconstructed with smoother_redux. test model changed accordingly. --- matlab/DsgeSmoother.m | 9 +++++++-- matlab/missing_DiffuseKalmanSmootherH1_Z.m | 7 ++++++- matlab/missing_DiffuseKalmanSmootherH3_Z.m | 8 ++++++-- .../fs2000_smoother_redux.mod | 14 ++------------ 4 files changed, 21 insertions(+), 17 deletions(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 59b3b4c6e..def6aaf23 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -265,7 +265,7 @@ if kalman_algo == 1 || kalman_algo == 3 a_initial = zeros(np,1); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=T*a_initial; %set state prediction for first Kalman step; - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, aalphahat0, state_uncertainty0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ... Z,R1,Q,H,Pinf,Pstar, ... data1,vobs,np,smpl,data_index, ... options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty,options_.filter_covariance,options_.smoother_redux); @@ -321,7 +321,7 @@ if kalman_algo == 2 || kalman_algo == 4 a_initial = zeros(np,1); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=ST*a_initial; %set state prediction for first Kalman step; - [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ... + [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, aalphahat0, state_uncertainty0, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ... Z,R1,Q,diag(H), ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ... options_.nk,kalman_tol,diffuse_kalman_tol, ... @@ -553,7 +553,9 @@ else static_var_list(static_var_list) = ~ilagged; % reconstruct smoothed variables aaa=zeros(M_.endo_nbr,gend+1); + if ~isempty(alphahat0) aaa(oo_.dr.restrict_var_list,1)=alphahat0; + end aaa(oo_.dr.restrict_var_list,2:end)=alphahat; for k=1:gend aaa(static_var_list,k+1) = C(~ilagged,:)*alphahat(:,k)+D(~ilagged,:)*etahat(:,k); @@ -576,6 +578,9 @@ else if any(ilagged) % bbb=zeros(M_.endo_nbr,gend); % bbb(oo_.dr.restrict_var_list,:)=aahat; + if ~isempty(aalphahat0) + aaa(static_var_list0,d+1) = Tstar(ilagged,:)*aalphahat0+Rstar(ilagged,:)*eehat(:,d+1); + end for k=d+2:gend aaa(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k); end diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index 3cbc664d6..cf1430a6a 100644 --- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m @@ -1,4 +1,4 @@ -function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d,alphahat0,V0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) +function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d,alphahat0,aalphahat0,V0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) % function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix. @@ -141,6 +141,7 @@ else V=[]; end alphahat0=[]; +aalphahat0=[]; V0=[]; t = 0; @@ -276,6 +277,10 @@ while notsteady && t1 aalphahat(:,t-1) = aha(:,1); - eetahat(:,t) = etaha(:,2); end + eetahat(:,t) = etaha(:,2); a(:,t) = ax(:,1); a1(:,t) = a1x(:,2); a1(:,t+1) = ax(:,2); @@ -469,6 +470,9 @@ while notsteady && t1.e-12 error('smoother_redux with kalman_algo=1 does not replicate original smoothed static variables!') end if max(merr1U)>1.e-12 - for k=1:length(vlist1) - merr2U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})(2:end)-oo1.UpdatedVariables.(vlist1{k})(2:end))); - end - if max(merr2U)>1.e-12 - error('smoother_redux with kalman_algo=1 does not replicate original updated static variables!') - end + error('smoother_redux with kalman_algo=1 does not replicate original updated static variables!') end if max(merr1F)>1.e-12 error('smoother_redux with kalman_algo=1 does not replicate original filtered static variables!') @@ -177,12 +172,7 @@ if max(merr1)>1.e-12 error('smoother_redux with kalman_algo=2 does not replicate original smoothed static variables!') end if max(merr1U)>1.e-12 - for k=1:length(vlist1) - merr2U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})(2:end)-oo2.UpdatedVariables.(vlist1{k})(2:end))); - end - if max(merr2U)>1.e-12 - error('smoother_redux with kalman_algo=2 does not replicate original updated static variables!') - end + error('smoother_redux with kalman_algo=2 does not replicate original updated static variables!') end if max(merr1F)>1.e-12 error('smoother_redux with kalman_algo=2 does not replicate original filtered static variables!') From 77e87be48b1b9fad48efb11820768259726132bb Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 26 Jan 2023 09:23:38 +0100 Subject: [PATCH 14/21] try with base regime before returning error --- matlab/+occbin/kalman_update_algo_1.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m index cfe97cf70..8381e7011 100644 --- a/matlab/+occbin/kalman_update_algo_1.m +++ b/matlab/+occbin/kalman_update_algo_1.m @@ -132,6 +132,10 @@ else end options_.occbin.simul=opts_simul; [~, out, ss] = occbin.solver(M_,oo_,options_); +if out.error_flag + options_.occbin.simul.init_regime=regimes0; + [~, out, ss] = occbin.solver(M_,oo_,options_); +end if out.error_flag error_flag = out.error_flag; etahat=etahat(:,2); From 035db8b8c4fbed490586f7e5782fd5fb852bec75 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 26 Jan 2023 09:24:17 +0100 Subject: [PATCH 15/21] bug fixes when using guess regime --- matlab/+occbin/solve_two_constraints.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m index 9869f6339..eab9c269b 100644 --- a/matlab/+occbin/solve_two_constraints.m +++ b/matlab/+occbin/solve_two_constraints.m @@ -164,7 +164,6 @@ for shock_period = 1:n_shocks_periods % unconstrained regime (nperiods_0) binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)]; end - binding_indicator_history{iter}=binding_indicator; if iter==1 && guess_history_it regime_1 = regime_history_guess(shock_period).regime1; @@ -181,6 +180,8 @@ for shock_period = 1:n_shocks_periods end nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required end + binding_indicator_history{iter}=binding_indicator; + % analyse violvec and isolate contiguous periods in the other regime. [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug); regime_history(shock_period).regime1 = regime_1; From 783c237d1799c7047bbf2293c2c6cce4a5c215ca Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Thu, 26 Jan 2023 09:27:07 +0100 Subject: [PATCH 16/21] evaluate set_prior only when needed --- matlab/evaluate_likelihood.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m index 315663333..d25bfc927 100644 --- a/matlab/evaluate_likelihood.m +++ b/matlab/evaluate_likelihood.m @@ -73,7 +73,6 @@ if isempty(dataset) end options_=select_qz_criterium_value(options_); -[~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0) % Plot prior densities. % Set prior bounds @@ -81,6 +80,7 @@ if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0) else % estimated parameters but no declared priors % No priors are declared so Dynare will estimate the model by % maximum likelihood with inequality constraints for the parameters. + [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); bounds.lb = lb; bounds.ub = ub; end From cdd195576ef0019efe4a161ebbc540cc95058b0a Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 3 Feb 2023 14:24:58 +0100 Subject: [PATCH 17/21] store binding_indicator_history once it has been set --- matlab/+occbin/solve_one_constraint.m | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m index a55ebd0e0..2f7834df9 100644 --- a/matlab/+occbin/solve_one_constraint.m +++ b/matlab/+occbin/solve_one_constraint.m @@ -152,9 +152,7 @@ for shock_period = 1:n_shocks_periods if length(binding_indicator)<(nperiods_0 + 1) binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)]; end - - binding_indicator_history{iter}=binding_indicator; - + if iter==1 && guess_history_it regime = regime_history_guess(shock_period).regime; regime_start = regime_history_guess(shock_period).regimestart; @@ -164,6 +162,7 @@ for shock_period = 1:n_shocks_periods end nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required end + binding_indicator_history{iter}=binding_indicator; % analyze when each regime starts based on current guess [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug); regime_history(shock_period).regime = regime; From b886de92dc99098cce4b3d388650839cbb6b2365 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 3 Feb 2023 14:25:32 +0100 Subject: [PATCH 18/21] trap case when simulation does not converge within smoother --- matlab/+occbin/DSGE_smoother.m | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m index 3ed80498c..59175ccf1 100644 --- a/matlab/+occbin/DSGE_smoother.m +++ b/matlab/+occbin/DSGE_smoother.m @@ -144,6 +144,11 @@ opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid options_.occbin.simul=opts_simul; options_.noprint = true; [~, out, ss] = occbin.solver(M_,oo_,options_); +if out.error_flag + fprintf('Occbin smoother:: simulation within smoother did not converge.\n') + print_info(error_flag, false, options_) + return; +end regime_history = out.regime_history; if options_.smoother_redux occbin_options.opts_simul.restrict_state_space =1; @@ -204,6 +209,11 @@ while is_changed && maxiter>iter && ~is_periodic opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1); options_.occbin.simul=opts_simul; [~, out, ss] = occbin.solver(M_,oo_,options_); + if out.error_flag + fprintf('Occbin smoother:: simulation within smoother did not converge.\n') + print_info(error_flag, false, options_) + return; + end regime_history = out.regime_history; TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:); RR = ss.R(oo_.dr.order_var,:,:); From f0aa2fb86fd91728186cf8c650a91a62c644b6d0 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Fri, 3 Feb 2023 14:25:59 +0100 Subject: [PATCH 19/21] cosmetic change --- matlab/DsgeSmoother.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index def6aaf23..3075e7a0b 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -554,7 +554,7 @@ else % reconstruct smoothed variables aaa=zeros(M_.endo_nbr,gend+1); if ~isempty(alphahat0) - aaa(oo_.dr.restrict_var_list,1)=alphahat0; + aaa(oo_.dr.restrict_var_list,1)=alphahat0; end aaa(oo_.dr.restrict_var_list,2:end)=alphahat; for k=1:gend From f565a0a84a2b5ef451fe29f1b1a30cd9379870c5 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 6 Feb 2023 11:39:24 +0100 Subject: [PATCH 20/21] bug fix and avoid missing output arguments and when smoother exits after non convergence of simulation. --- matlab/+occbin/DSGE_smoother.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m index 59175ccf1..ead33b243 100644 --- a/matlab/+occbin/DSGE_smoother.m +++ b/matlab/+occbin/DSGE_smoother.m @@ -146,7 +146,8 @@ options_.noprint = true; [~, out, ss] = occbin.solver(M_,oo_,options_); if out.error_flag fprintf('Occbin smoother:: simulation within smoother did not converge.\n') - print_info(error_flag, false, options_) + print_info(out.error_flag, options_.noprint, options_) + oo_.occbin.smoother.error_flag=1; return; end regime_history = out.regime_history; @@ -211,7 +212,8 @@ while is_changed && maxiter>iter && ~is_periodic [~, out, ss] = occbin.solver(M_,oo_,options_); if out.error_flag fprintf('Occbin smoother:: simulation within smoother did not converge.\n') - print_info(error_flag, false, options_) + print_info(out.error_flag, false, options_) + oo_.occbin.smoother.error_flag=1; return; end regime_history = out.regime_history; From 63a299f64e658ddcbf61c025584308f9a370c6ba Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Mon, 6 Feb 2023 14:58:42 +0100 Subject: [PATCH 21/21] new error codes introduced and applied for occbin smoother. - 321 when simulation within occbin smoother fails - 322 when occbin smoother does not converge. --- matlab/+occbin/DSGE_smoother.m | 6 +++--- matlab/dynare_estimation_1.m | 8 ++++++-- matlab/get_error_message.m | 4 ++++ matlab/prior_posterior_statistics_core.m | 5 +++++ 4 files changed, 18 insertions(+), 5 deletions(-) diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m index ead33b243..656e84a3c 100644 --- a/matlab/+occbin/DSGE_smoother.m +++ b/matlab/+occbin/DSGE_smoother.m @@ -147,7 +147,7 @@ options_.noprint = true; if out.error_flag fprintf('Occbin smoother:: simulation within smoother did not converge.\n') print_info(out.error_flag, options_.noprint, options_) - oo_.occbin.smoother.error_flag=1; + oo_.occbin.smoother.error_flag=321; return; end regime_history = out.regime_history; @@ -213,7 +213,7 @@ while is_changed && maxiter>iter && ~is_periodic if out.error_flag fprintf('Occbin smoother:: simulation within smoother did not converge.\n') print_info(out.error_flag, false, options_) - oo_.occbin.smoother.error_flag=1; + oo_.occbin.smoother.error_flag=321; return; end regime_history = out.regime_history; @@ -363,7 +363,7 @@ if (maxiter==iter && is_changed) || is_periodic else fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n') oo_.occbin.smoother=[]; - oo_.occbin.smoother.error_flag=1; + oo_.occbin.smoother.error_flag=322; end else disp('occbin.DSGE_smoother: smoother converged.') diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 9de11f35b..8bb1337ae 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -194,8 +194,10 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_. else if options_.occbin.smoother.status [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info); - if oo_.occbin.smoother.error_flag==0 + if oo_.occbin.smoother.error_flag(1)==0 [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); + else + fprintf('\nOccbin: smoother did not succeed. No results will be written to oo_.\n') end else [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); @@ -619,8 +621,10 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha else if options_.occbin.smoother.status [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info); - if oo_.occbin.smoother.error_flag==0 + if oo_.occbin.smoother.error_flag(1)==0 [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); + else + fprintf('\nOccbin: smoother did not succeed. No results will be written to oo_.\n') end else [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_); diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m index 57063a7d4..4f57c869b 100644 --- a/matlab/get_error_message.m +++ b/matlab/get_error_message.m @@ -192,6 +192,10 @@ switch info(1) message = 'Occbin: Simulation did not converge -- infinite loop of guess regimes'; case 320 message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.'; + case 321 + message = 'Occbin: there was a problem in running the smoother. Simulation within smoother failed.'; + case 322 + message = 'Occbin: smoother did not converge.'; case 401 message = 'Cycle reduction reached the iteration limit. Try increasing maxit.'; case 402 diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m index 1125328e2..7d04d2845 100644 --- a/matlab/prior_posterior_statistics_core.m +++ b/matlab/prior_posterior_statistics_core.m @@ -226,8 +226,13 @@ for b=fpar:B if options_.occbin.smoother.status opts_local.occbin.simul.waitbar=0; opts_local.occbin.smoother.waitbar = 0; + opts_local.occbin.smoother.linear_smoother=false; % speed-up [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = ... occbin.DSGE_smoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_); + if oo_.occbin.smoother.error_flag(1) + message=get_error_message(oo_.occbin.smoother.error_flag,opts_local); + fprintf('\nprior_posterior_statistics: One of the draws failed with the error:\n%s\n',message) + end else [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = ... DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);