diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 36ab1f0c917bacb5e9a8048f7b4b7f215692d8ff..3ed80498ca5c29e803c93ac8abced89d9213c628 100644
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
+function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
%function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,PKK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DSGE_smoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_, dataset_info)
% Runs a DSGE smoother with occasionally binding constraints
%
@@ -39,6 +39,8 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P
% - oo_ [structure] storing the results
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
+% - alphahat0 [double] (m*1) array, smoothed endogenous variables in period 0 (a_{0|T}) (decision-rule order)
+% - state_uncertainty0 [double] (K,K,1) array, storing the uncertainty in period 0
% Copyright © 2021 Dynare Team
%
@@ -118,7 +120,7 @@ occbin_options.first_period_occbin_update = options_.occbin.smoother.first_perio
occbin_options.opts_regime = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history;
-[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
+[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history;
regime_history = oo_.occbin.smoother.regime_history;
@@ -135,9 +137,9 @@ opts_regime.binding_indicator=[];
regime_history0 = regime_history;
fprintf('Occbin smoother iteration 1.\n')
-opts_simul.SHOCKS = [etahat(:,2:end)'; zeros(1,M_.exo_nbr)];
+opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
opts_simul.exo_pos = 1:M_.exo_nbr;
-opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
+opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
options_.occbin.simul=opts_simul;
options_.noprint = true;
@@ -145,15 +147,18 @@ options_.noprint = true;
regime_history = out.regime_history;
if options_.smoother_redux
occbin_options.opts_simul.restrict_state_space =1;
- oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
- oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1);
+ [T0,R0] = dynare_resolve(M_,options_,oo_);
+ oo_.occbin.linear_smoother.T0=T0;
+ oo_.occbin.linear_smoother.R0=R0;
+ % oo_.occbin.linear_smoother.T0=ss.T(oo_.dr.order_var,oo_.dr.order_var,1);
+ % oo_.occbin.linear_smoother.R0=ss.R(oo_.dr.order_var,:,1);
end
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:);
-TT = cat(3,TT(:,:,1),TT);
-RR = cat(3,RR(:,:,1),RR);
-CC = cat(2,CC(:,1),CC);
+% TT = cat(3,TT(:,:,1),TT);
+% RR = cat(3,RR(:,:,1),RR);
+% CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history;
opts_regime.binding_indicator = [];
@@ -184,7 +189,7 @@ while is_changed && maxiter>iter && ~is_periodic
iter=iter+1;
fprintf('Occbin smoother iteration %u.\n', iter)
occbin_options.opts_regime.regime_history=regime_history;
- [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
+ [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0, diffuse_steps] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history;
if occbin_smoother_debug
@@ -195,17 +200,17 @@ while is_changed && maxiter>iter && ~is_periodic
sto_RR = RR;
sto_TT = TT;
- opts_simul.SHOCKS = [etahat(:,2:end)'; zeros(1,M_.exo_nbr)];
- opts_simul.endo_init = alphahat(oo_.dr.inv_order_var,1);
+ opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
+ opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,oo_,options_);
regime_history = out.regime_history;
TT = ss.T(oo_.dr.order_var,oo_.dr.order_var,:);
RR = ss.R(oo_.dr.order_var,:,:);
CC = ss.C(oo_.dr.order_var,:);
- TT = cat(3,TT(:,:,1),TT);
- RR = cat(3,RR(:,:,1),RR);
- CC = cat(2,CC(:,1),CC);
+% TT = cat(3,TT(:,:,1),TT);
+% RR = cat(3,RR(:,:,1),RR);
+% CC = cat(2,CC(:,1),CC);
opts_regime.regime_history = regime_history;
[TT, RR, CC, regime_history] = occbin.check_regimes(TT, RR, CC, opts_regime, M_, oo_, options_);
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 7ea0cd462d5c220a625b129fa79d2961cc687af7..e184f34ca5eba22c594743ae27030fcde5bb9aef 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
+function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty,M_,oo_,bayestopt_,alphahat0,state_uncertainty0,d] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,varargin)
% Estimation of the smoothed variables and innovations.
%
% INPUTS
@@ -265,7 +265,7 @@ if kalman_algo == 1 || kalman_algo == 3
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=T*a_initial; %set state prediction for first Kalman step;
- [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ...
+ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ...
Z,R1,Q,H,Pinf,Pstar, ...
data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty,options_.filter_covariance,options_.smoother_redux);
@@ -321,7 +321,7 @@ if kalman_algo == 2 || kalman_algo == 4
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=ST*a_initial; %set state prediction for first Kalman step;
- [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
+ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, alphahat0, state_uncertainty0, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ...
@@ -538,17 +538,20 @@ else
static_var_list0(static_var_list) = ilagged;
static_var_list(static_var_list) = ~ilagged;
% reconstruct smoothed variables
- aaa=zeros(M_.endo_nbr,gend);
- aaa(oo_.dr.restrict_var_list,:)=alphahat;
+ aaa=zeros(M_.endo_nbr,gend+1);
+ aaa(oo_.dr.restrict_var_list,1)=alphahat0;
+ aaa(oo_.dr.restrict_var_list,2:end)=alphahat;
for k=1:gend
- aaa(static_var_list,k) = C(~ilagged,:)*alphahat(:,k)+D(~ilagged,:)*etahat(:,k);
+ aaa(static_var_list,k+1) = C(~ilagged,:)*alphahat(:,k)+D(~ilagged,:)*etahat(:,k);
end
if any(ilagged)
for k=2:gend
- aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k);
+ aaa(static_var_list0,k+1) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k);
end
+ aaa(static_var_list0,2) = Tstar(ilagged,:)*alphahat0+Rstar(ilagged,:)*etahat(:,1);
end
- alphahat=aaa;
+ alphahat0=aaa(:,1);
+ alphahat=aaa(:,2:end);
% reconstruct updated variables
aaa=zeros(M_.endo_nbr,gend);
@@ -619,6 +622,18 @@ else
state_uncertainty=sstate_uncertainty;
clear sstate_uncertainty
end
+ if ~isempty(state_uncertainty0)
+ mm=size(T,1);
+ ss=length(find(static_var_list));
+ sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr);
+ sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list)=state_uncertainty0(1:mm,1:mm);
+ sstate_uncertainty(static_var_list,static_var_list)=[C(~ilagged,:) D(~ilagged,:)]*state_uncertainty0*[C(~ilagged,:) D(~ilagged,:)]';
+ tmp = [C(~ilagged,:) D(~ilagged,:)]*state_uncertainty0;
+ sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list)=tmp(1:ss,1:mm);
+ sstate_uncertainty(oo_.dr.restrict_var_list,static_var_list)=transpose(sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list));
+ state_uncertainty0=sstate_uncertainty;
+ clear sstate_uncertainty
+ end
% reconstruct PK
if ~isempty(PK)
diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 8e0fad922e5192312baf1ad6ad658e857eb3db49..09d5672637c0a1c1639139f003057e05261af733 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -1,4 +1,4 @@
-function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux)
+function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d,alphahat0,V0] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux)
% function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V,aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH1_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag,filter_covariance_flag,smoother_redux)
% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix.
@@ -140,8 +140,21 @@ if state_uncertainty_flag
else
V=[];
end
+alphahat0=[];
+V0=[];
t = 0;
+if rank(Pinf(:,:,1),diffuse_kalman_tol)
+ % this is needed to get smoothed states in period 0 with diffuse steps
+ % i.e. period 0 is a filter step without observables
+ Pinf_init = Pinf(:,:,1);
+ Pstar_init = Pstar(:,:,1);
+ a_init = a(:,1);
+ a(:,t) = T*a(:,t);
+ % only non-stationary part is affected by following line,
+ % hence Pstar on EXIT from diffuse step will NOT change.
+ Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ;
+end
while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
t = t+1;
di = data_index{t};
@@ -162,7 +175,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
@@ -339,8 +352,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 +387,7 @@ if d %diffuse periods
if ~Finf_singular(1,t)
r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf
r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ...
- + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1)
+ + Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1)
if state_uncertainty_flag
L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1)
N(:,:,t)=Linf(:,:,t)'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.19, noting that L^(0) is named Linf
@@ -374,7 +404,7 @@ if d %diffuse periods
r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14)
if state_uncertainty_flag
N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)...
- +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14)
+ +Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14)
N_1(:,:,t)=T'*N_1(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14)
N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14)
end
@@ -395,8 +425,8 @@ if d %diffuse periods
V(:,:,t) = pstmp - pstmp*ntmp*pstmp...
-(pitmp*ntmp1*pstmp)'...
- pitmp*ntmp1*pstmp...
- - pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30
-
+ - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30
+
else
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
@@ -405,6 +435,33 @@ if d %diffuse periods
end
end
end
+ % compute states and covarinace in period t=0
+ r10 = T'*r1(:,1);
+ r00 = T'*r0(:,1);
+ alphahat0 = a_init + Pstar_init*r00 + Pinf_init*r10; % DK (2012), eq. 5.23
+ if state_uncertainty_flag
+ N0=T'*N(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf
+ N_10=T'*N_1(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf
+ N_20=T'*N_2(:,:,1)*T; % DK (2012), eq. 5.19, noting that L^(0) is named Linf
+ if smoother_redux
+ pstmp = [Pstar_init R*Q; (R*Q)' Q];
+ pitmp = [Pinf_init zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp = [N0 zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp1 = [N_10 zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp2 = [N_20 zeros(mm,rr); zeros(rr,mm+rr)];
+ V0 = pstmp - pstmp*ntmp*pstmp...
+ -(pitmp*ntmp1*pstmp)'...
+ - pitmp*ntmp1*pstmp...
+ - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30
+
+ else
+ V0 = Pstar_init-Pstar_init*N0*Pstar_init...
+ -(Pinf_init*N_10*Pstar_init)'...
+ - Pinf_init*N_10*Pstar_init...
+ - Pinf_init*N_20*Pinf_init; % DK (2012), eq. 5.30
+ end
+ end
+
end
if decomp_flag
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index 8f100f5286e1730b9b3a0674345fd334ea83cccf..b8fcd416c93c67900b53c674c96d15789cc7ac59 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -1,4 +1,4 @@
-function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d,varargout] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_)
+function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d,alphahat0,V0,varargout] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_)
% function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V, aalphahat,eetahat,d] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag, filter_covariance_flag, smoother_redux, occbin_)
% Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
% Univariate treatment of multivariate time series.
@@ -149,6 +149,8 @@ if state_uncertainty_flag
else
V=[];
end
+alphahat0=[];
+V0=[];
if ~occbin_.status
isoccbin = 0;
@@ -180,7 +182,7 @@ else
T0=occbin_.info{5};
R0=occbin_.info{6};
else
-
+
TT=occbin_.info{5};
RR=occbin_.info{6};
CC=occbin_.info{7};
@@ -201,7 +203,7 @@ else
CC=repmat(zeros(mm,1),1,smpl+1);
end
end
-
+
else
TT=repmat(T,1,1,smpl+1);
RR=repmat(R,1,1,smpl+1);
@@ -210,7 +212,7 @@ else
if ~smoother_redux
T0=T;
R0=R;
-
+
end
if ~isinf(occbin_options.first_period_occbin_update)
% initialize state matrices (otherwise they are set to 0 for
@@ -219,7 +221,7 @@ else
RRR=repmat(R0,1,1,smpl+1);
CCC=repmat(zeros(length(T0),1),1,smpl+1);
end
-
+
end
t = 0;
@@ -229,6 +231,14 @@ if ~isempty(Pinf(:,:,1))
else
newRank = rank(Pinf(:,:,1),diffuse_kalman_tol);
end
+if newRank
+ % add this to get smoothed states in period 0
+ Pinf_init = Pinf(:,:,1);
+ Pstar_init = Pstar(:,:,1);
+ Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ;
+ ainit = a1(:,1);
+end
+
while newRank && t < smpl
t = t+1;
a(:,t) = a1(:,t);
@@ -270,7 +280,7 @@ while newRank && t < smpl
else
oldRank = 0;
end
- if isoccbin,
+ if isoccbin
TT(:,:,t+1)= T;
RR(:,:,t+1)= R;
end
@@ -298,7 +308,11 @@ while newRank && t < smpl
end
if isoccbin
- first_period_occbin_update = max(t+2,occbin_options.first_period_occbin_update);
+ first_period_occbin_update = occbin_options.first_period_occbin_update;
+ if d>0
+ first_period_occbin_update = max(t+2,occbin_options.first_period_occbin_update);
+ % kalman update is not yet robust to accommodate diffuse steps
+ end
if occbin_options.opts_regime.waitbar
hh = dyn_waitbar(0,'Occbin: Piecewise Kalman Filter');
set(hh,'Name','Occbin: Piecewise Kalman Filter.');
@@ -322,6 +336,10 @@ Pinf1 = Pinf1(:,:,1:d);
notsteady = 1;
while notsteady && t<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 +347,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
@@ -372,6 +412,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)
@@ -474,13 +526,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 +626,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 +714,8 @@ if d
V(:,:,t) = pstmp - pstmp*ntmp0*pstmp...
-(pitmp*ntmp1*pstmp)'...
- pitmp*ntmp1*pstmp...
- - pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30
-
+ - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30
+
else
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
@@ -654,14 +727,41 @@ if d
r0(:,t-1) = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
r1(:,t-1) = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
if state_uncertainty_flag
- N_0(:,:,t-1)= T'*N_0(:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T
- N_1(:,:,t-1)= T'*N_1(:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T
- N_2(:,:,t-1)= T'*N_2(:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T
+ N_0(:,:,t-1)= T'*N_0(:,:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T
+ N_1(:,:,t-1)= T'*N_1(:,:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T
+ N_2(:,:,t-1)= T'*N_2(:,:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T
+ end
+ else
+ r00 = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
+ r10 = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
+ if state_uncertainty_flag
+ N_00= T'*N_0(:,:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T
+ N_10= T'*N_1(:,:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T
+ N_20= T'*N_2(:,:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T
end
end
end
-else
- alphahat0 = 0*a1(:,1) + P1(:,:,1)*ri;
+ % get smoothed states in t=0
+ alphahat0 = ainit + Pstar_init*r00 + Pinf_init*r10; % KD (2000), eq. (26)
+ if state_uncertainty_flag
+ if smoother_redux
+ pstmp = [Pstar_init R*Q; (R*Q)' Q];
+ pitmp = [Pinf_init zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp0 = [N_00 zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp1 = [N_10 zeros(mm,rr); zeros(rr,mm+rr)];
+ ntmp2 = [N_20 zeros(mm,rr); zeros(rr,mm+rr)];
+ V0 = pstmp - pstmp*ntmp0*pstmp...
+ -(pitmp*ntmp1*pstmp)'...
+ - pitmp*ntmp1*pstmp...
+ - pitmp*ntmp2*pitmp; % DK (2012), eq. 5.30
+
+ else
+ V0=Pstar_init-Pstar_init*N_00*Pstar_init...
+ -(Pinf_init*N_10*Pstar_init)'...
+ - Pinf_init*N_10*Pstar_init...
+ - Pinf_init*N_20*Pinf_init; % DK (2012), eq. 5.30
+ end
+ end
end
if decomp_flag
@@ -675,7 +775,7 @@ if decomp_flag
ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
end
end
-
+
% calculate eta_tm1t
if isoccbin
if isqvec