diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 36ab1f0c917bacb5e9a8048f7b4b7f215692d8ff..656e84a3ce6fbb9e455208997f26e836d9700971 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,25 +137,34 @@ 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;
[~, out, ss] = occbin.solver(M_,oo_,options_);
+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=321;
+ return;
+end
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 +195,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 +206,23 @@ 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_);
+ 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=321;
+ 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,:,:);
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_);
@@ -346,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/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m
index cfe97cf70a4926176e688be291d504983be0e29b..8381e701110bbedcd76120ca777f495e19012f1b 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);
diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index 17e74a323448935fcdb9a57d6bb23a3f2c59798d..34f31b47f485515713ad1e107dde9f4c1d0185f3 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 = [];
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index f5f0c1446c79ff6edf9b1c829e8b98c90f2216ca..2f7834df9d589ea14b7c7f497f30411d432eaad2 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
@@ -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;
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 80744ac50943fa24b2e3ee285cb7ad3d3c587152..eab9c269b6614321f7e9cbd40a36567957834c52 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
@@ -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;
diff --git a/matlab/+occbin/solver.m b/matlab/+occbin/solver.m
index 9dd44f88fa68a76464de38df9d1c838a29f1166c..3dc465ba3cfeb102ed811a6e78cd4fade92d275e 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;
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 7ea0cd462d5c220a625b129fa79d2961cc687af7..3075e7a0b32e1a3c3635b5a3b36ef00a729afcbb 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, 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, 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, ...
@@ -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;
@@ -538,17 +552,22 @@ 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);
+ 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) = 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);
@@ -559,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
@@ -619,6 +641,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/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 9de11f35bdeec6339a639078f4501eb7599b83df..8bb1337aeeac801bcbc98683ef86fcac38a202a9 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/evaluate_likelihood.m b/matlab/evaluate_likelihood.m
index 652a258cc9e8bc1ebed3d3731cd138cb63f36c38..d25bfc927884c3b69e215971b9f0a7b45bddcc5c 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_);
+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.
+ [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_);
+ 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
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 57063a7d4bb2dc005b92d14b2fb05318945250c7..4f57c869b1c42cac7be76b2c08a57ab0636151ca 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/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
index ba5145fadc2f03b0df534f6f295d455ebc23d1c3..e966ffaa3b8b137e20c72d7042a9830eb568fe8f 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;
@@ -141,6 +145,10 @@ while notsteady && t<=last
if ~(isqvec)
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
end
+ if t==1
+ Pinit = P1(:,:,1);
+ ainit = a1(:,1);
+ end
end
s = t-start+1;
d_index = data_index{t};
@@ -178,61 +186,81 @@ while notsteady && t<=last
% badly_conditioned_F = true;
end
end
+ if badly_conditioned_F && (~occbin_.status || (occbin_.status && t<first_period_occbin_update))
+ if ~all(abs(F(:))<kalman_tol)
+ % Use univariate filter.
+ return
+ else
+ % Pathological case, discard draw
+ return
+ end
+ else
+ F_singular = false;
+ end
if ~occbin_.status || (occbin_.status && (options_.occbin.likelihood.use_updated_regime==0 || t<first_period_occbin_update))
- if badly_conditioned_F && (~occbin_.status || (occbin_.status && t<first_period_occbin_update))
- if ~all(abs(F(:))<kalman_tol)
- % Use univariate filter.
- return
- else
- % Pathological case, discard draw
- return
- end
+ if rescale_prediction_error_covariance
+ log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
+ iF = inv(F./(sig*sig'))./(sig*sig');
+ rescale_prediction_error_covariance=rescale_prediction_error_covariance0;
else
- F_singular = false;
- if rescale_prediction_error_covariance
- log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
- iF = inv(F./(sig*sig'))./(sig*sig');
- rescale_prediction_error_covariance=rescale_prediction_error_covariance0;
- else
- log_dF = log(det(F));
- iF = inv(F);
- end
- lik(s) = log_dF + transpose(v)*iF*v + length(d_index)*log(2*pi);
- if t<first_period_occbin_update
- if Zflag
- K = P*z'*iF;
- if occbin_.status
- P0(:,:,t) = (P-K*z*P);
- end
-
- 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
+ log_dF = log(det(F));
+ iF = inv(F);
+ end
+ lik(s) = log_dF + transpose(v)*iF*v + length(d_index)*log(2*pi);
+ if t<first_period_occbin_update
+ if Zflag
+ K = P*z'*iF;
if occbin_.status
- a0(:,t) = (a+K*v);
- vv(d_index,t) = v;
+ P0(:,:,t) = (P-K*z*P);
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(:);
+
+ 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
diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 8e0fad922e5192312baf1ad6ad658e857eb3db49..cf1430a6abe9bd6cf8f2a9668a4d9e1f917b5cf8 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,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.
@@ -140,8 +140,22 @@ if state_uncertainty_flag
else
V=[];
end
+alphahat0=[];
+aalphahat0=[];
+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(:,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;
+end
while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
t = t+1;
di = data_index{t};
@@ -162,7 +176,7 @@ while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
Finf = ZZ*Pinf(:,:,t)*ZZ'; % (5.7) in DK (2012)
if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0
- % The univariate diffuse kalman filter should be used.
+ % The univariate diffuse kalman filter should be used.
alphahat = Inf;
return
else %rank of F_{\infty,t} is 0
@@ -263,6 +277,10 @@ while notsteady && t<smpl
eetahat(:,st) = QRt*ri; %DK (2012), eq. 4.63
end
end
+ if t==1
+ ri = T'*ri; %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0
+ aalphahat0 = P(:,:,1)*ri; %DK (2012), eq. 4.35
+ end
end
a(:,t+1) = T*atilde(:,t);
if filter_covariance_flag
@@ -339,8 +357,25 @@ while t>d+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 +392,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 +409,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 +430,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 +440,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 8f100f5286e1730b9b3a0674345fd334ea83cccf..a5894922aa5815144f5e6bdc7767a99d1611b16b 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,aalphahat0,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,9 @@ if state_uncertainty_flag
else
V=[];
end
+alphahat0=[];
+aalphahat0=[];
+V0=[];
if ~occbin_.status
isoccbin = 0;
@@ -180,7 +183,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 +204,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 +213,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 +222,7 @@ else
RRR=repmat(R0,1,1,smpl+1);
CCC=repmat(zeros(length(T0),1),1,smpl+1);
end
-
+
end
t = 0;
@@ -229,6 +232,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 +281,7 @@ while newRank && t < smpl
else
oldRank = 0;
end
- if isoccbin,
+ if isoccbin
TT(:,:,t+1)= T;
RR(:,:,t+1)= R;
end
@@ -298,7 +309,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 +337,10 @@ Pinf1 = Pinf1(:,:,1:d);
notsteady = 1;
while notsteady && t<smpl
t = t+1;
+ if t==1
+ Pinit = P(:,:,1);
+ ainit = a1(:,1);
+ end
a(:,t) = a1(:,t);
P1(:,:,t) = P(:,:,t);
di = data_index{t}';
@@ -329,11 +348,33 @@ while notsteady && t<smpl
if waitbar_indicator
dyn_waitbar(t/smpl, hh, sprintf('Period %u of %u', t,smpl));
end
- if isqvec
- Qt = Qvec(:,:,t-1:t+1);
- end
occbin_options.opts_regime.waitbar=0;
- [ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),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,kalman_tol,nk);
+ if t==1
+ if isqvec
+ Qt = cat(3,Q,Qvec(:,:,t:t+1));
+ end
+ a0 = a(:,1);
+ a10 = [a0 a(:,1)];
+ P0 = P(:,:,1);
+ P10 = P1(:,:,[1 1]);
+ data_index0{1}=[];
+ data_index0(2)=data_index(1);
+ v0(:,2)=v(:,1);
+ Y0(:,2)=Y(:,1);
+ Y0(:,1)=nan;
+ Fi0 = Fi(:,1);
+ Ki0 = Ki(:,:,1);
+ 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, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a10,P0,P10,data_index0,Z,v0,Fi0,Ki0,Y0,H,Qt,T0,R0,TT01,RR01,CC01,regimes_(t:t+1),M_,oo_,options_,occbin_options,kalman_tol,nk);
+ else
+ if isqvec
+ Qt = Qvec(:,:,t-1:t+1);
+ end
+ [ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),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,kalman_tol,nk);
+ end
if ~error_flag
regimes_(t:t+2)=tmp;
else
@@ -347,10 +388,10 @@ while notsteady && t<smpl
return
end
- if smoother_redux
+ if smoother_redux && t>1
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);
@@ -372,6 +413,18 @@ while notsteady && t<smpl
aK(jnk,:,t+jnk) = ax(:,1+jnk);
end
else
+ if isoccbin && t==1
+ if isqvec
+ QQ = RR(:,:,t)*Qvec(:,:,t)*transpose(RR(:,:,t));
+ else
+ QQ = RR(:,:,t)*Q*transpose(RR(:,:,t));
+ end
+ T = TT(:,:,t);
+ C = CC(:,t);
+ a1(:,t) = T*a(:,t)+C; %transition according to (6.14) in DK (2012)
+ P(:,:,t) = T*P(:,:,t)*T' + QQ; %transition according to (6.14) in DK (2012)
+ P1(:,:,t) = P(:,:,t);
+ end
for i=di
Zi = Z(i,:);
v(i,t) = Y(i,t) - Zi*a(:,t); % nu_{t,i} in 6.13 in DK (2012)
@@ -417,6 +470,9 @@ while notsteady && t<smpl
end
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
end
+ if t==1
+ aalphahat0 = P1(:,:,st)*ri;
+ end
end
if isoccbin
if isqvec
@@ -474,13 +530,13 @@ while notsteady && t<smpl
end
end
if waitbar_indicator
- dyn_waitbar_close(hh);
+ dyn_waitbar_close(hh);
end
P1(:,:,t+1) = P(:,:,t+1);
if ~isinf(first_period_occbin_update) && isoccbin
- regimes_ = regimes_(2:smpl+1);
+ regimes_ = regimes_(1:smpl+1);
else
regimes_ = struct();
TTT=TT;
@@ -574,7 +630,28 @@ while 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 +718,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 +731,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 +779,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
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 1125328e2c5a0fb6542106a9ed182771b9fefe2d..7d04d28450e4df46fbf949220aa6eeaf5116c63f 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_);
diff --git a/tests/kalman_filter_smoother/fs2000_smoother_redux.mod b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod
index 72d261a59a2988fa62c0abeede02bd66edacde6f..ee13d7a2c4401d8edc725c81363fd01c77db13c6 100644
--- a/tests/kalman_filter_smoother/fs2000_smoother_redux.mod
+++ b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod
@@ -112,28 +112,13 @@ 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)
- 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
- 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,28 +169,13 @@ 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)
- 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
- 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
diff --git a/tests/occbin/filter/NKM.mod b/tests/occbin/filter/NKM.mod
index 8c04a966e9aecbbf09de9cfabf43eedee06f4268..81b25d1d7e54dc2d7cdcaf38c2f129bad377618b 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);