stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'occbin_init_smo' into 'master' 

occbin enhancements and bug fixes

See merge request Dynare/dynare!2123
silicon
Sébastien Villemot 2023-02-07 14:23:14 +00:00
parent c3524d33d0 63a299f64e
commit d6370a6d28
16 changed files with 436 additions and 178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
matlab/+occbin/DSGE_smoother.m
View File

@ -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) %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 % 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 % - oo_ [structure] storing the results
% - options_ [structure] describing the options % - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors % - 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 % 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 = 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.binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history; 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; oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history;
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; regime_history0 = regime_history;
fprintf('Occbin smoother iteration 1.\n') 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.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! opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
options_.noprint = true; options_.noprint = true;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, 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; regime_history = out.regime_history;
if options_.smoother_redux if options_.smoother_redux
occbin_options.opts_simul.restrict_state_space =1; occbin_options.opts_simul.restrict_state_space =1;
oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1); [T0,R0] = dynare_resolve(M_,options_,oo_);
oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1); 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 end
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:); TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:); RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:); CC = ss.C(oo_.dr.order_var,:);
TT = cat(3,TT(:,:,1),TT); % TT = cat(3,TT(:,:,1),TT);
RR = cat(3,RR(:,:,1),RR); % RR = cat(3,RR(:,:,1),RR);
CC = cat(2,CC(:,1),CC); % CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history; opts_regime.regime_history = regime_history;
opts_regime.binding_indicator = []; opts_regime.binding_indicator = [];
@ -184,7 +195,7 @@ while is_changed && maxiter>iter && ~is_periodic
iter=iter+1; iter=iter+1;
fprintf('Occbin smoother iteration %u.\n', iter) fprintf('Occbin smoother iteration %u.\n', iter)
occbin_options.opts_regime.regime_history=regime_history; 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}; sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history; regime_history0(iter,:) = regime_history;
if occbin_smoother_debug if occbin_smoother_debug
@ -195,17 +206,23 @@ while is_changed && maxiter>iter && ~is_periodic
sto_RR = RR; sto_RR = RR;
sto_TT = TT; sto_TT = TT;
opts_simul.SHOCKS = [etahat(:,2:end)'; zeros(1,M_.exo_nbr)]; opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(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);
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, 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; regime_history = out.regime_history;
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:); TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:); RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:); CC = ss.C(oo_.dr.order_var,:);
TT = cat(3,TT(:,:,1),TT); % TT = cat(3,TT(:,:,1),TT);
RR = cat(3,RR(:,:,1),RR); % RR = cat(3,RR(:,:,1),RR);
CC = cat(2,CC(:,1),CC); % CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history; opts_regime.regime_history = regime_history;
[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_); [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 else
fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n') fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
oo_.occbin.smoother=[]; oo_.occbin.smoother=[];
oo_.occbin.smoother.error_flag=1; oo_.occbin.smoother.error_flag=322;
end end
else else
disp('occbin.DSGE_smoother: smoother converged.') disp('occbin.DSGE_smoother: smoother converged.')

4
matlab/+occbin/kalman_update_algo_1.m
View File

@ -132,6 +132,10 @@ else
end end
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_); [~, 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 if out.error_flag
error_flag = out.error_flag; error_flag = out.error_flag;
etahat=etahat(:,2); etahat=etahat(:,2);

4
matlab/+occbin/set_default_options.m
View File

@ -67,7 +67,7 @@ end
if ismember(flag,{'likelihood','all'}) if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.curb_retrench = false; 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.full_output = false;
options_occbin_.likelihood.IF_likelihood = false; options_occbin_.likelihood.IF_likelihood = false;
options_occbin_.likelihood.init_regime_history = []; options_occbin_.likelihood.init_regime_history = [];
@ -187,7 +187,7 @@ if ismember(flag,{'smoother','all'})
options_occbin_.smoother.curb_retrench = false; options_occbin_.smoother.curb_retrench = false;
options_occbin_.smoother.debug = false; options_occbin_.smoother.debug = false;
options_occbin_.smoother.fast = 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.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_mode = 1; % 0 = standard; 1 = unconditional frcsts zero shocks+smoothed states in each period
options_occbin_.smoother.init_regime_history = []; options_occbin_.smoother.init_regime_history = [];

7
matlab/+occbin/solve_one_constraint.m
View File

@ -105,7 +105,7 @@ end
SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); 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.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
SS_out.C = nan(DM.n_vars,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(); regime_history = struct();
guess_history = false; guess_history = false;
else else
@ -152,9 +152,7 @@ for shock_period = 1:n_shocks_periods
if length(binding_indicator)<(nperiods_0 + 1) if length(binding_indicator)<(nperiods_0 + 1)
binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)]; binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)];
end end
binding_indicator_history{iter}=binding_indicator;
if iter==1 && guess_history_it if iter==1 && guess_history_it
regime = regime_history_guess(shock_period).regime; regime = regime_history_guess(shock_period).regime;
regime_start = regime_history_guess(shock_period).regimestart; regime_start = regime_history_guess(shock_period).regimestart;
@ -164,6 +162,7 @@ for shock_period = 1:n_shocks_periods
end end
nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
end end
binding_indicator_history{iter}=binding_indicator;
% analyze when each regime starts based on current guess % analyze when each regime starts based on current guess
[regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug); [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
regime_history(shock_period).regime = regime; regime_history(shock_period).regime = regime;

5
matlab/+occbin/solve_two_constraints.m
View File

@ -115,7 +115,7 @@ end
SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); 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.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
SS_out.C = NaN(DM.n_vars,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(); regime_history = struct();
guess_history = false; guess_history = false;
else else
@ -164,7 +164,6 @@ for shock_period = 1:n_shocks_periods
% unconstrained regime (nperiods_0) % unconstrained regime (nperiods_0)
binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)]; binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)];
end end
binding_indicator_history{iter}=binding_indicator;
if iter==1 && guess_history_it if iter==1 && guess_history_it
regime_1 = regime_history_guess(shock_period).regime1; regime_1 = regime_history_guess(shock_period).regime1;
@ -181,6 +180,8 @@ for shock_period = 1:n_shocks_periods
end end
nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
end end
binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime. % 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_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
regime_history(shock_period).regime1 = regime_1; regime_history(shock_period).regime1 = regime_1;

3
matlab/+occbin/solver.m
View File

@ -54,6 +54,9 @@ if solve_dr
if isempty(options_.qz_criterium) if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6; options_.qz_criterium = 1+1e-6;
end end
if options_.order>1
options_.order = 1;
end
[dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_); [dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_);
out.error_flag=error_flag; out.error_flag=error_flag;

62
matlab/DsgeSmoother.m
View File

@ -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. % Estimation of the smoothed variables and innovations.
% %
% INPUTS % INPUTS
@ -265,7 +265,7 @@ if kalman_algo == 1 || kalman_algo == 3
a_initial = zeros(np,1); a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=T*a_initial; %set state prediction for first Kalman step; 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, ... Z,R1,Q,H,Pinf,Pstar, ...
data1,vobs,np,smpl,data_index, ... 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); 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 = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=ST*a_initial; %set state prediction for first Kalman step; 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), ... Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ... options_.nk,kalman_tol,diffuse_kalman_tol, ...
@ -391,8 +391,9 @@ else
opts_simul = options_.occbin.simul; opts_simul = options_.occbin.simul;
end end
% reconstruct smoothed variables % reconstruct smoothed variables
aaa=zeros(M_.endo_nbr,gend); aaa=zeros(M_.endo_nbr,gend+1);
aaa(oo_.dr.restrict_var_list,:)=alphahat; aaa(oo_.dr.restrict_var_list,1)=alphahat0;
aaa(oo_.dr.restrict_var_list,2:end)=alphahat;
iTx = zeros(size(TTx)); iTx = zeros(size(TTx));
for k=1:gend for k=1:gend
if isoccbin if isoccbin
@ -416,13 +417,19 @@ else
static_var_list0 = static_var_list; static_var_list0 = static_var_list;
static_var_list0(static_var_list) = ilagged; static_var_list0(static_var_list) = ilagged;
static_var_list(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); aaa(static_var_list,k+1) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged);
if any(ilagged) && k>1 if any(ilagged)
aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(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
end end
alphahat=aaa; alphahat0=aaa(:,1);
alphahat=aaa(:,2:end);
% reconstruct updated variables % reconstruct updated variables
bbb=zeros(M_.endo_nbr,gend); bbb=zeros(M_.endo_nbr,gend);
@ -488,6 +495,13 @@ else
state_uncertainty=sstate_uncertainty; state_uncertainty=sstate_uncertainty;
clear sstate_uncertainty clear sstate_uncertainty
end 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 = zeros(nk,M_.endo_nbr,gend+nk);
aaa(:,oo_.dr.restrict_var_list,:)=aK; aaa(:,oo_.dr.restrict_var_list,:)=aK;
@ -538,17 +552,22 @@ else
static_var_list0(static_var_list) = ilagged; static_var_list0(static_var_list) = ilagged;
static_var_list(static_var_list) = ~ilagged; static_var_list(static_var_list) = ~ilagged;
% reconstruct smoothed variables % reconstruct smoothed variables
aaa=zeros(M_.endo_nbr,gend); aaa=zeros(M_.endo_nbr,gend+1);
aaa(oo_.dr.restrict_var_list,:)=alphahat; 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 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 end
if any(ilagged) if any(ilagged)
for k=2:gend 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 end
aaa(static_var_list0,2) = Tstar(ilagged,:)*alphahat0+Rstar(ilagged,:)*etahat(:,1);
end end
alphahat=aaa; alphahat0=aaa(:,1);
alphahat=aaa(:,2:end);
% reconstruct updated variables % reconstruct updated variables
aaa=zeros(M_.endo_nbr,gend); aaa=zeros(M_.endo_nbr,gend);
@ -559,6 +578,9 @@ else
if any(ilagged) if any(ilagged)
% bbb=zeros(M_.endo_nbr,gend); % bbb=zeros(M_.endo_nbr,gend);
% bbb(oo_.dr.restrict_var_list,:)=aahat; % 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 for k=d+2:gend
aaa(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k); aaa(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k);
end end
@ -619,6 +641,18 @@ else
state_uncertainty=sstate_uncertainty; state_uncertainty=sstate_uncertainty;
clear sstate_uncertainty clear sstate_uncertainty
end 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 % reconstruct PK
if ~isempty(PK) if ~isempty(PK)

8
matlab/dynare_estimation_1.m
View File

@ -194,8 +194,10 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
else else
if options_.occbin.smoother.status 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); [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); [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 end
else 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_); [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 else
if options_.occbin.smoother.status 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); [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); [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 end
else 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_); [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_);

55
matlab/evaluate_likelihood.m
View File

@ -47,24 +47,24 @@ end
if ischar(parameters) if ischar(parameters)
switch parameters switch parameters
case 'posterior mode' case 'posterior mode'
parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_); parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
case 'posterior mean' case 'posterior mean'
parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
case 'posterior median' case 'posterior median'
parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_); parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
case 'prior mode' case 'prior mode'
parameters = bayestopt_.p5(:); parameters = bayestopt_.p5(:);
case 'prior mean' case 'prior mean'
parameters = bayestopt_.p1; parameters = bayestopt_.p1;
otherwise otherwise
disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:') disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:')
disp(' ''posterior mode'', ') disp(' ''posterior mode'', ')
disp(' ''posterior mean'', ') disp(' ''posterior mean'', ')
disp(' ''posterior median'', ') disp(' ''posterior median'', ')
disp(' ''prior mode'' or') disp(' ''prior mode'' or')
disp(' ''prior mean''.') disp(' ''prior mean''.')
error error
end end
end end
@ -73,10 +73,23 @@ if isempty(dataset)
end end
options_=select_qz_criterium_value(options_); 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 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_); llik = -occbin.IVF_posterior(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
else else
llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_); llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
end end
ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_); ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_);
llik = llik - ldens; llik = llik - ldens;

4
matlab/get_error_message.m
View File

@ -192,6 +192,10 @@ switch info(1)
message = 'Occbin: Simulation did not converge -- infinite loop of guess regimes'; message = 'Occbin: Simulation did not converge -- infinite loop of guess regimes';
case 320 case 320
message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.'; 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 case 401
message = 'Cycle reduction reached the iteration limit. Try increasing maxit.'; message = 'Cycle reduction reached the iteration limit. Try increasing maxit.';
case 402 case 402

126
matlab/kalman/likelihood/missing_observations_kalman_filter.m
View File

@ -105,7 +105,11 @@ if occbin_.status
occbin_options=occbin_.info{4}; occbin_options=occbin_.info{4};
opts_regime.regime_history = occbin_options.opts_simul.init_regime; opts_regime.regime_history = occbin_options.opts_simul.init_regime;
opts_regime.binding_indicator = occbin_options.opts_simul.init_binding_indicator; 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) if isempty(opts_regime.binding_indicator) && isempty(opts_regime.regime_history)
opts_regime.binding_indicator=zeros(last+2,M_.occbin.constraint_nbr); opts_regime.binding_indicator=zeros(last+2,M_.occbin.constraint_nbr);
end end
@ -123,8 +127,8 @@ if occbin_.status
TT=repmat(T,1,1,last+1); TT=repmat(T,1,1,last+1);
RR=repmat(R,1,1,last+1); RR=repmat(R,1,1,last+1);
CC=repmat(zeros(mm,1),1,last+1); CC=repmat(zeros(mm,1),1,last+1);
end end
end end
else else
first_period_occbin_update = inf; first_period_occbin_update = inf;
@ -141,6 +145,10 @@ while notsteady && t<=last
if ~(isqvec) if ~(isqvec)
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 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
end end
s = t-start+1; s = t-start+1;
d_index = data_index{t}; d_index = data_index{t};
@ -178,61 +186,81 @@ while notsteady && t<=last
% badly_conditioned_F = true; % badly_conditioned_F = true;
end end
end 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 badly_conditioned_F && (~occbin_.status || (occbin_.status && t<first_period_occbin_update)) if ~all(abs(F(:))<kalman_tol)
if ~all(abs(F(:))<kalman_tol) % Use univariate filter.
% Use univariate filter. return
return
else
% Pathological case, discard draw
return
end
else else
F_singular = false; % Pathological case, discard draw
if rescale_prediction_error_covariance return
log_dF = log(det(F./(sig*sig')))+2*sum(log(sig)); end
iF = inv(F./(sig*sig'))./(sig*sig'); else
rescale_prediction_error_covariance=rescale_prediction_error_covariance0; F_singular = false;
else end
log_dF = log(det(F)); if ~occbin_.status || (occbin_.status && (options_.occbin.likelihood.use_updated_regime==0 || t<first_period_occbin_update))
iF = inv(F); if rescale_prediction_error_covariance
end log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
lik(s) = log_dF + transpose(v)*iF*v + length(d_index)*log(2*pi); iF = inv(F./(sig*sig'))./(sig*sig');
if t<first_period_occbin_update rescale_prediction_error_covariance=rescale_prediction_error_covariance0;
if Zflag else
K = P*z'*iF; log_dF = log(det(F));
if occbin_.status iF = inv(F);
P0(:,:,t) = (P-K*z*P); end
end lik(s) = log_dF + transpose(v)*iF*v + length(d_index)*log(2*pi);
if t<first_period_occbin_update
P = T*(P-K*z*P)*transpose(T)+QQ; if Zflag
else K = P*z'*iF;
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 if occbin_.status
a0(:,t) = (a+K*v); P0(:,:,t) = (P-K*z*P);
vv(d_index,t) = v;
end end
a = T*(a+K*v)+C;
if t>=no_more_missing_observations && ~isqvec && ~occbin_.status P = T*(P-K*z*P)*transpose(T)+QQ;
notsteady = max(abs(K(:)-oldK))>riccati_tol; else
oldK = K(:); K = P(:,z)*iF;
if occbin_.status
P0(:,:,t) = (P-K*P(z,:));
end 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
end end
end end
if occbin_.status && t>=first_period_occbin_update if occbin_.status && t>=first_period_occbin_update
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
occbin_options.opts_simul.waitbar=0; 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 info
if options_.debug if options_.debug
fprintf('\nmissing_observations_kalman_filter:PKF failed with: %s\n', get_error_message(info,options_)); 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); P0(:,:,t) = Px(:,:,1);
P1(:,:,t) = P1x(:,:,2); P1(:,:,t) = P1x(:,:,2);
P = Px(:,:,2); P = Px(:,:,2);
end end
t = t+1; t = t+1;
end end

78
matlab/missing_DiffuseKalmanSmootherH1_Z.m
View File

@ -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) % 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. % 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 else
V=[]; V=[];
end end
alphahat0=[];
aalphahat0=[];
V0=[];
t = 0; 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 while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
t = t+1; t = t+1;
di = data_index{t}; 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) Finf = ZZ*Pinf(:,:,t)*ZZ'; % (5.7) in DK (2012)
if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0 if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 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; alphahat = Inf;
return return
else %rank of F_{\infty,t} is 0 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 eetahat(:,st) = QRt*ri; %DK (2012), eq. 4.63
end end
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 end
a(:,t+1) = T*atilde(:,t); a(:,t+1) = T*atilde(:,t);
if filter_covariance_flag if filter_covariance_flag
@ -339,8 +357,25 @@ while t>d+1
end end
end end
if d %diffuse periods if d==0 % get smoother in period t=0
% initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23 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 = zeros(mm,d+1);
r0(:,d+1) = r(:,d+1); %set r0_{d}, i.e. shifted by one period 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 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) if ~Finf_singular(1,t)
r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf 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)) ... 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 if state_uncertainty_flag
L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1) 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 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) r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14)
if state_uncertainty_flag if state_uncertainty_flag
N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)... 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_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) N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14)
end end
@ -395,8 +430,8 @@ if d %diffuse periods
V(:,:,t) = pstmp - pstmp*ntmp*pstmp... V(:,:,t) = pstmp - pstmp*ntmp*pstmp...
-(pitmp*ntmp1*pstmp)'... -(pitmp*ntmp1*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 else
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)... V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'... -(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
@ -405,6 +440,33 @@ if d %diffuse periods
end end
end 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 end
if decomp_flag if decomp_flag

152
matlab/missing_DiffuseKalmanSmootherH3_Z.m
View File

@ -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_) % 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. % Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
% Univariate treatment of multivariate time series. % Univariate treatment of multivariate time series.
@ -149,6 +149,9 @@ if state_uncertainty_flag
else else
V=[]; V=[];
end end
alphahat0=[];
aalphahat0=[];
V0=[];
if ~occbin_.status if ~occbin_.status
isoccbin = 0; isoccbin = 0;
@ -180,7 +183,7 @@ else
T0=occbin_.info{5}; T0=occbin_.info{5};
R0=occbin_.info{6}; R0=occbin_.info{6};
else else
TT=occbin_.info{5}; TT=occbin_.info{5};
RR=occbin_.info{6}; RR=occbin_.info{6};
CC=occbin_.info{7}; CC=occbin_.info{7};
@ -201,7 +204,7 @@ else
CC=repmat(zeros(mm,1),1,smpl+1); CC=repmat(zeros(mm,1),1,smpl+1);
end end
end end
else else
TT=repmat(T,1,1,smpl+1); TT=repmat(T,1,1,smpl+1);
RR=repmat(R,1,1,smpl+1); RR=repmat(R,1,1,smpl+1);
@ -210,7 +213,7 @@ else
if ~smoother_redux if ~smoother_redux
T0=T; T0=T;
R0=R; R0=R;
end end
if ~isinf(occbin_options.first_period_occbin_update) if ~isinf(occbin_options.first_period_occbin_update)
% initialize state matrices (otherwise they are set to 0 for % initialize state matrices (otherwise they are set to 0 for
@ -219,7 +222,7 @@ else
RRR=repmat(R0,1,1,smpl+1); RRR=repmat(R0,1,1,smpl+1);
CCC=repmat(zeros(length(T0),1),1,smpl+1); CCC=repmat(zeros(length(T0),1),1,smpl+1);
end end
end end
t = 0; t = 0;
@ -229,6 +232,14 @@ if ~isempty(Pinf(:,:,1))
else else
newRank = rank(Pinf(:,:,1),diffuse_kalman_tol); newRank = rank(Pinf(:,:,1),diffuse_kalman_tol);
end 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 while newRank && t < smpl
t = t+1; t = t+1;
a(:,t) = a1(:,t); a(:,t) = a1(:,t);
@ -270,7 +281,7 @@ while newRank && t < smpl
else else
oldRank = 0; oldRank = 0;
end end
if isoccbin, if isoccbin
TT(:,:,t+1)= T; TT(:,:,t+1)= T;
RR(:,:,t+1)= R; RR(:,:,t+1)= R;
end end
@ -298,7 +309,11 @@ while newRank && t < smpl
end end
if isoccbin 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 if occbin_options.opts_regime.waitbar
hh = dyn_waitbar(0,'Occbin: Piecewise Kalman Filter'); hh = dyn_waitbar(0,'Occbin: Piecewise Kalman Filter');
set(hh,'Name','Occbin: Piecewise Kalman Filter.'); set(hh,'Name','Occbin: Piecewise Kalman Filter.');
@ -322,6 +337,10 @@ Pinf1 = Pinf1(:,:,1:d);
notsteady = 1; notsteady = 1;
while notsteady && t<smpl while notsteady && t<smpl
t = t+1; t = t+1;
if t==1
Pinit = P(:,:,1);
ainit = a1(:,1);
end
a(:,t) = a1(:,t); a(:,t) = a1(:,t);
P1(:,:,t) = P(:,:,t); P1(:,:,t) = P(:,:,t);
di = data_index{t}'; di = data_index{t}';
@ -329,11 +348,33 @@ while notsteady && t<smpl
if waitbar_indicator if waitbar_indicator
dyn_waitbar(t/smpl, hh, sprintf('Period %u of %u', t,smpl)); dyn_waitbar(t/smpl, hh, sprintf('Period %u of %u', t,smpl));
end end
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
occbin_options.opts_regime.waitbar=0; 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 if ~error_flag
regimes_(t:t+2)=tmp; regimes_(t:t+2)=tmp;
else else
@ -347,10 +388,10 @@ while notsteady && t<smpl
return return
end end
if smoother_redux if smoother_redux && t>1
aalphahat(:,t-1) = aha(:,1); aalphahat(:,t-1) = aha(:,1);
eetahat(:,t) = etaha(:,2);
end end
eetahat(:,t) = etaha(:,2);
a(:,t) = ax(:,1); a(:,t) = ax(:,1);
a1(:,t) = a1x(:,2); a1(:,t) = a1x(:,2);
a1(:,t+1) = ax(:,2); a1(:,t+1) = ax(:,2);
@ -372,6 +413,18 @@ while notsteady && t<smpl
aK(jnk,:,t+jnk) = ax(:,1+jnk); aK(jnk,:,t+jnk) = ax(:,1+jnk);
end end
else 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 for i=di
Zi = Z(i,:); Zi = Z(i,:);
v(i,t) = Y(i,t) - Zi*a(:,t); % nu_{t,i} in 6.13 in DK (2012) 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 end
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}} ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
end end
if t==1
aalphahat0 = P1(:,:,st)*ri;
end
end end
if isoccbin if isoccbin
if isqvec if isqvec
@ -474,13 +530,13 @@ while notsteady && t<smpl
end end
end end
if waitbar_indicator if waitbar_indicator
dyn_waitbar_close(hh); dyn_waitbar_close(hh);
end end
P1(:,:,t+1) = P(:,:,t+1); P1(:,:,t+1) = P(:,:,t+1);
if ~isinf(first_period_occbin_update) && isoccbin if ~isinf(first_period_occbin_update) && isoccbin
regimes_ = regimes_(2:smpl+1); regimes_ = regimes_(1:smpl+1);
else else
regimes_ = struct(); regimes_ = struct();
TTT=TT; 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}} Ni = T'*Ni*T; % KD (2000), eq. (23), equation for N_{t-1,p_{t-1}}
end end
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 = zeros(mm,d);
r0(:,d) = ri; r0(:,d) = ri;
r1 = zeros(mm,d); r1 = zeros(mm,d);
@ -641,8 +718,8 @@ if d
V(:,:,t) = pstmp - pstmp*ntmp0*pstmp... V(:,:,t) = pstmp - pstmp*ntmp0*pstmp...
-(pitmp*ntmp1*pstmp)'... -(pitmp*ntmp1*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 else
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)... V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,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} 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} 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 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_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_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_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 end
end end
else % get smoothed states in t=0
alphahat0 = 0*a1(:,1) + P1(:,:,1)*ri; 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 end
if decomp_flag 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,:)'; ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
end end
end end
% calculate eta_tm1t % calculate eta_tm1t
if isoccbin if isoccbin
if isqvec if isqvec

5
matlab/prior_posterior_statistics_core.m
View File

@ -226,8 +226,13 @@ for b=fpar:B
if options_.occbin.smoother.status if options_.occbin.smoother.status
opts_local.occbin.simul.waitbar=0; opts_local.occbin.simul.waitbar=0;
opts_local.occbin.smoother.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_] = ... [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_); 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 else
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,~,~,P,~,~,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = ... [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_); DsgeSmoother(deep,gend,Y,data_index,missing_value,M_,oo_,opts_local,bayestopt_,estim_params_);

42
tests/kalman_filter_smoother/fs2000_smoother_redux.mod
View File

@ -112,28 +112,13 @@ for k=1:length(vlist1)
merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo1.FilteredVariables.(vlist1{k}))); merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo1.FilteredVariables.(vlist1{k})));
end end
if max(merr1)>1.e-12 if max(merr1)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=1 does not replicate original smoothed static variables!')
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
end end
if max(merr1U)>1.e-12 if max(merr1U)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=1 does not replicate original updated static variables!')
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
end end
if max(merr1F)>1.e-12 if max(merr1F)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=1 does not replicate original filtered static variables!')
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
end end
merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo1.FilteredVariablesKStepAhead)))); merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo1.FilteredVariablesKStepAhead))));
if max(merrK)>1.e-12 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}))); merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo2.FilteredVariables.(vlist1{k})));
end end
if max(merr1)>1.e-12 if max(merr1)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=2 does not replicate original smoothed static variables!')
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
end end
if max(merr1U)>1.e-12 if max(merr1U)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=2 does not replicate original updated static variables!')
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
end end
if max(merr1F)>1.e-12 if max(merr1F)>1.e-12
for k=1:length(vlist1) error('smoother_redux with kalman_algo=2 does not replicate original filtered static variables!')
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
end end
merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo2.FilteredVariablesKStepAhead)))); merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo2.FilteredVariablesKStepAhead))));
if max(merrK)>1.e-12 if max(merrK)>1.e-12

10
tests/occbin/filter/NKM.mod
View File

@ -319,6 +319,16 @@ varobs yg inom pi;
mh_replic=0, plot_priors=0, smoother, smoother_redux, mh_replic=0, plot_priors=0, smoother, smoother_redux,
nodisplay,consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty); 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) // use inversion filter (note that IF provides smoother together with likelihood)
occbin_setup(likelihood_inversion_filter,smoother_inversion_filter); occbin_setup(likelihood_inversion_filter,smoother_inversion_filter);