Commit abe8a05b authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Merge branch 'smoother_redux' into 'master'

Implement new option smoother_redux, to allow fast smoother for very large...

See merge request !1859
parents c4e0cd05 5fa62659
Pipeline #5384 passed with stages
in 74 minutes and 8 seconds
......@@ -6330,6 +6330,26 @@ block decomposition of the model (see :opt:`block`).
(:math:`E_{t}{y_t}`). See below for a description of all these
variables.
.. option:: smoother_redux
Triggers a faster computation of the smoothed endogenous variables and shocks for large models.
It runs the smoother only for the state variables (i.e. with the same representation used for
likelihood computations) and computes the remaining variables ex-post.
Static unobserved objects (filtered, smoothed, updated, k-step ahead) are recovered, but there are
exceptions to a full recovery, depending on how static unobserved variables depend on the restricted
state space adopted. For example, lagged shocks which are ONLY used to recover NON-observed static
variables will not be recovered).
For such exceptions, only the following output is provided:
``FilteredVariablesKStepAhead``: will be fully recovered
``SmoothedVariables``, ``FilteredVariables``, ``UpdatedVariables``: recovered for all periods beyond period ``d+1``,
where ``d`` denotes the number of diffuse filtering steps.
``FilteredVariablesKStepAheadVariances``, ``Variance``, and ``State_uncertainty`` cannot be recovered, and ZERO is provided as output.
If you need variances for those variables, either do not set the option, or declare the variable as observed, using NaNs as data points.
.. option:: forecast = INTEGER
Computes the posterior distribution of a forecast on INTEGER
......@@ -8417,6 +8437,18 @@ Dynare can also run the smoother on a calibrated model:
See :opt:`filter_decomposition`.
.. option:: filter_covariance
See :opt:`filter_covariance`.
.. option:: smoother_redux
See :opt:`smoother_redux`.
.. option:: kalman_algo = INTEGER
See :opt:`kalman_algo <kalman_algo = INTEGER>`.
.. option:: diffuse_filter = INTEGER
See :opt:`diffuse_filter`.
......
......@@ -97,22 +97,30 @@ end
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
%store old setting of restricted var_list
oldoo.restrict_var_list = oo_.dr.restrict_var_list;
oldoo.restrict_columns = oo_.dr.restrict_columns;
oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if ~options_.smoother_redux
%store old setting of restricted var_list
oldoo.restrict_var_list = oo_.dr.restrict_var_list;
oldoo.restrict_columns = oo_.dr.restrict_columns;
oo_.dr.restrict_var_list = bayestopt_.smoother_var_list;
oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns;
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
%get location of observed variables and requested smoothed variables in
%decision rules
bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
else
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
bayestopt_.mf = bayestopt_.mf1;
end
if info~=0
print_info(info,options_.noprint, options_);
return
end
%get location of observed variables and requested smoothed variables in
%decision rules
bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
if options_.noconstant
constant = zeros(vobs,1);
else
......@@ -233,10 +241,10 @@ 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] = missing_DiffuseKalmanSmootherH1_Z(a_initial,ST, ...
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d] = 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_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty,options_.filter_covariance,options_.smoother_redux);
if isinf(alphahat)
if kalman_algo == 1
fprintf('\nDsgeSmoother: Switching to univariate filter. This may be a sign of stochastic singularity.\n')
......@@ -289,14 +297,14 @@ 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] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d] = 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, ...
options_.filter_decomposition,options_.smoothed_state_uncertainty);
options_.filter_decomposition,options_.smoothed_state_uncertainty,options_.filter_covariance,options_.smoother_redux);
end
if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_algo == 4)
% extracting measurement errors
% removing observed variables from the state vector
......@@ -320,10 +328,122 @@ if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_al
end
end
%reset old setting of restricted var_list
oo_.dr.restrict_var_list = oldoo.restrict_var_list;
oo_.dr.restrict_columns = oldoo.restrict_columns;
if ~options_.smoother_redux
%reset old setting of restricted var_list
oo_.dr.restrict_var_list = oldoo.restrict_var_list;
oo_.dr.restrict_columns = oldoo.restrict_columns;
else
if options_.block == 0
ic = [ M_.nstatic+(1:M_.nspred) M_.endo_nbr+(1:size(oo_.dr.ghx,2)-M_.nspred) ]';
else
ic = oo_.dr.restrict_columns;
end
[A,B] = kalman_transition_matrix(oo_.dr,(1:M_.endo_nbr)',ic,M_.exo_nbr);
iT = pinv(T);
Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
C = Tstar*iT;
D = Rstar-C*R;
static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
ilagged = any(abs(C*T-Tstar)'>1.e-12);
static_var_list0 = static_var_list;
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;
for k=1:gend
aaa(static_var_list,k) = 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);
end
end
alphahat=aaa;
% reconstruct updated variables
aaa=zeros(M_.endo_nbr,gend);
aaa(oo_.dr.restrict_var_list,:)=ahat;
for k=1:gend
aaa(static_var_list,k) = C(~ilagged,:)*ahat(:,k)+D(~ilagged,:)*eehat(:,k);
end
if any(ilagged)
for k=d+2:gend
aaa(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k);
end
end
ahat1=aaa;
% reconstruct aK
aaa = zeros(options_.nk,M_.endo_nbr,gend+options_.nk);
aaa(:,oo_.dr.restrict_var_list,:)=aK;
for k=1:gend
for jnk=1:options_.nk
aaa(jnk,static_var_list,k+jnk) = C(~ilagged,:)*dynare_squeeze(aK(jnk,:,k+jnk));
end
end
if any(ilagged)
for k=1:gend
aaa(1,static_var_list0,k+1) = Tstar(ilagged,:)*ahat(:,k);
for jnk=2:options_.nk
aaa(jnk,static_var_list0,k+jnk) = Tstar(ilagged,:)*dynare_squeeze(aK(jnk-1,:,k+jnk-1));
end
end
end
aK=aaa;
ahat=ahat1;
% reconstruct P
if ~isempty(P)
PP=zeros(M_.endo_nbr,M_.endo_nbr,gend+1);
PP(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:)=P;
DQD=D(~ilagged,:)*Q*transpose(D(~ilagged,:))+C(~ilagged,:)*R*Q*transpose(D(~ilagged,:))+D(~ilagged,:)*Q*transpose(C(~ilagged,:)*R);
DQR=D(~ilagged,:)*Q*transpose(R);
for k=1:gend+1
PP(static_var_list,static_var_list,k)=C(~ilagged,:)*P(:,:,k)*C(~ilagged,:)'+DQD;
PP(static_var_list,oo_.dr.restrict_var_list,k)=C(~ilagged,:)*P(:,:,k)+DQR;
PP(oo_.dr.restrict_var_list,static_var_list,k)=transpose(PP(static_var_list,oo_.dr.restrict_var_list,k));
end
P=PP;
clear('PP');
end
% reconstruct state_uncertainty
if ~isempty(state_uncertainty)
mm=size(T,1);
ss=length(find(static_var_list));
sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr,gend);
sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:)=state_uncertainty(1:mm,1:mm,:);
for k=1:gend
sstate_uncertainty(static_var_list,static_var_list,k)=[C(~ilagged,:) D(~ilagged,:)]*state_uncertainty(:,:,k)*[C(~ilagged,:) D(~ilagged,:)]';
tmp = [C(~ilagged,:) D(~ilagged,:)]*state_uncertainty(:,:,k);
sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list,k)=tmp(1:ss,1:mm);
sstate_uncertainty(oo_.dr.restrict_var_list,static_var_list,k)=transpose(sstate_uncertainty(static_var_list,oo_.dr.restrict_var_list,k));
end
state_uncertainty=sstate_uncertainty;
clear('sstate_uncertainty');
end
% reconstruct PK TO DO!!
if ~isempty(PK)
PP = zeros(options_.nk,M_.endo_nbr,M_.endo_nbr,gend+options_.nk);
PP(:,oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:) = PK;
DQD=D(~ilagged,:)*Q*transpose(D(~ilagged,:))+C(~ilagged,:)*R*Q*transpose(D(~ilagged,:))+D(~ilagged,:)*Q*transpose(C(~ilagged,:)*R);
DQR=D(~ilagged,:)*Q*transpose(R);
for f=1:options_.nk
for k=1:gend
PP(f,static_var_list,static_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))*C(~ilagged,:)'+DQD;
PP(f,static_var_list,oo_.dr.restrict_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))+DQR;
PP(f,oo_.dr.restrict_var_list,static_var_list,k+f)=transpose(squeeze(PP(f,static_var_list,oo_.dr.restrict_var_list,k+f)));
end
end
PK=PP;
clear('PP');
end
bayestopt_.mf = bayestopt_.smoother_var_list(bayestopt_.smoother_mf);
end
function a=set_Kalman_smoother_starting_values(a,M_,oo_,options_)
% function a=set_Kalman_smoother_starting_values(a,M_,oo_,options_)
......@@ -352,3 +472,4 @@ if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
end
end
end
......@@ -428,6 +428,7 @@ options_.prefilter = 0;
options_.presample = 0;
options_.prior_trunc = 1e-10;
options_.smoother = false;
options_.smoother_redux = false;
options_.posterior_max_subsample_draws = 1200;
options_.sub_draws = [];
options_.ME_plot_tol=1e-6;
......
function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = 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)
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,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
% 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.
%
% INPUTS
......@@ -16,27 +16,34 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = missing_DiffuseK
% pp: number of observed variables
% mm: number of state variables
% smpl: sample size
% data_index [cell] 1*smpl cell of column vectors of indices.
% nk number of forecasting periods
% kalman_tol tolerance for reciprocal condition number
% diffuse_kalman_tol tolerance for reciprocal condition number (for Finf) and the rank of Pinf
% decomp_flag if true, compute filter decomposition
% state_uncertainty_flag if true, compute uncertainty about smoothed
% data_index: [cell] 1*smpl cell of column vectors of indices.
% nk: number of forecasting periods
% kalman_tol: tolerance for reciprocal condition number
% diffuse_kalman_tol: tolerance for reciprocal condition number (for Finf) and the rank of Pinf
% decomp_flag: if true, compute filter decomposition
% state_uncertainty_flag: if true, compute uncertainty about smoothed
% state estimate
% decomp_flag: if true, compute filter decomposition
% filter_covariance_flag: if true, compute filter covariance
% smoother_redux: if true, compute smoother on restricted
% state space, recover static variables from this
%
% OUTPUTS
% alphahat: smoothed variables (a_{t|T})
% epsilonhat:smoothed measurement errors
% etahat: smoothed shocks
% atilde: matrix of updated variables (a_{t|t})
% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)
% (meaningless for periods 1:d)
% P: 3D array of one-step ahead forecast error variance
% matrices
% PK: 4D array of k-step ahead forecast error variance
% matrices (meaningless for periods 1:d)
% decomp: decomposition of the effect of shocks on filtered values
% V: 3D array of state uncertainty matrices
% alphahat: smoothed variables (a_{t|T})
% epsilonhat: smoothed measurement errors
% etahat: smoothed shocks
% atilde: matrix of updated variables (a_{t|t})
% P: 3D array of one-step ahead forecast error variance
% matrices
% aK: 3D array of k step ahead filtered state variables (a_{t+k|t)
% (meaningless for periods 1:d)
% PK: 4D array of k-step ahead forecast error variance
% matrices (meaningless for periods 1:d)
% decomp: decomposition of the effect of shocks on filtered values
% V: 3D array of state uncertainty matrices
% aalphahat: filtered states in t-1|t
% eetahat: updated shocks in t|t
% d: number of diffuse periods
%
% Notes:
% Outputs are stored in decision-rule order, i.e. to get variables in order of declaration
......@@ -50,7 +57,7 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = missing_DiffuseK
% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
% Second Edition, Ch. 5
% Copyright (C) 2004-2018 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
......@@ -82,7 +89,11 @@ a = zeros(mm,smpl+1);
a(:,1) = a_initial;
atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+nk);
PK = zeros(nk,mm,mm,smpl+nk);
if filter_covariance_flag
PK = zeros(nk,mm,mm,smpl+nk);
else
PK = [];
end
iF = zeros(pp,pp,smpl);
Fstar = zeros(pp,pp,smpl);
iFstar = zeros(pp,pp,smpl);
......@@ -109,11 +120,22 @@ QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
if smoother_redux
aalphahat = alphahat;
eetahat = etahat;
else
aalphahat = [];
eetahat = [];
end
epsilonhat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
Finf_singular = zeros(1,smpl);
if state_uncertainty_flag
V = zeros(mm,mm,smpl);
if smoother_redux
V = zeros(mm+rr,mm+rr,smpl);
else
V = zeros(mm,mm,smpl);
end
N = zeros(mm,mm,smpl+1);
else
V=[];
......@@ -222,20 +244,45 @@ while notsteady && t<smpl
L(:,:,t) = T-K(:,di,t)*ZZ;
P(:,:,t+1) = T*P(:,:,t)*L(:,:,t)' + QQ;
end
if smoother_redux
ri=zeros(mm,1);
for st=t:-1:max(d+1,t-1)
di = data_index{st};
if isempty(di)
% in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
ri = L(:,:,t)'*ri; %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0
else
ZZ = Z(di,:);
ri = ZZ'*iF(di,di,st)*v(di,st) + L(:,:,st)'*ri; %compute r_{t-1}, DK (2012), eq. 4.38
end
if st==t-1
% get states in t-1|t
aalphahat(:,st) = a(:,st) + P(:,:,st)*ri; %DK (2012), eq. 4.35
else
% get shocks in t|t
eetahat(:,st) = QRt*ri; %DK (2012), eq. 4.63
end
end
end
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t+1);
if filter_covariance_flag
Pf = P(:,:,t+1);
end
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
if jnk>1
Pf = T*Pf*T' + QQ;
if filter_covariance_flag
if jnk>1
Pf = T*Pf*T' + QQ;
end
PK(jnk,:,:,t+jnk) = Pf;
end
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
end
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
end
% $$$ if t<smpl
% $$$ PZI_s = PZI;
% $$$ K_s = K(:,:,t);
......@@ -282,9 +329,16 @@ while t>d+1
end
etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63
if state_uncertainty_flag
V(:,:,t) = P(:,:,t)-P(:,:,t)*N(:,:,t)*P(:,:,t); %DK (2012), eq. 4.43
if smoother_redux
ptmp = [P(:,:,t) R*Q; (R*Q)' Q];
ntmp = [N(:,:,t) zeros(mm,rr); zeros(rr,mm+rr)];
V(:,:,t) = ptmp - ptmp*ntmp*ptmp;
else
V(:,:,t) = P(:,:,t)-P(:,:,t)*N(:,:,t)*P(:,:,t); %DK (2012), eq. 4.43
end
end
end
if d %diffuse periods
% initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23
r0 = zeros(mm,d+1);
......@@ -332,10 +386,23 @@ if d %diffuse periods
end
etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135
if state_uncertainty_flag
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
- Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)...
- Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30
if smoother_redux
pstmp = [Pstar(:,:,t) R*Q; (R*Q)' Q];
pitmp = [Pinf(:,:,t) zeros(mm,rr); zeros(rr,mm+rr)];
ntmp = [N(:,:,t) zeros(mm,rr); zeros(rr,mm+rr)];
ntmp1 = [N_1(:,:,t) zeros(mm,rr); zeros(rr,mm+rr)];
ntmp2 = [N_2(:,:,t) zeros(mm,rr); zeros(rr,mm+rr)];
V(:,:,t) = pstmp - pstmp*ntmp*pstmp...
-(pitmp*ntmp1*pstmp)'...
- pitmp*ntmp1*pstmp...
- pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30
else
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
- Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)...
- Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30
end
end
end
end
......
function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = 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)
% function [alphahat,epsilonhat,etahat,a1,P1,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
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)
% 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)
% Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
% Univariate treatment of multivariate time series.
%
......@@ -16,13 +16,17 @@ function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalma
% pp: number of observed variables
% mm: number of state variables
% smpl: sample size
% data_index [cell] 1*smpl cell of column vectors of indices.
% nk number of forecasting periods
% kalman_tol tolerance for zero divider
% diffuse_kalman_tol tolerance for zero divider
% decomp_flag if true, compute filter decomposition
% state_uncertainty_flag if true, compute uncertainty about smoothed
% data_index: [cell] 1*smpl cell of column vectors of indices.
% nk: number of forecasting periods
% kalman_tol: tolerance for zero divider
% diffuse_kalman_tol: tolerance for zero divider
% decomp_flag: if true, compute filter decomposition
% state_uncertainty_flag: if true, compute uncertainty about smoothed
% state estimate
% decomp_flag: if true, compute filter decomposition
% filter_covariance_flag: if true, compute filter covariance
% smoother_redux: if true, compute smoother on restricted
% state space, recover static variables from this
%
% OUTPUTS
% alphahat: smoothed state variables (a_{t|T})
......@@ -37,6 +41,9 @@ function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalma
% matrices (meaningless for periods 1:d)
% decomp: decomposition of the effect of shocks on filtered values
% V: 3D array of state uncertainty matrices
% aalphahat: filtered states in t-1|t
% eetahat: updated shocks in t|t
% d: number of diffuse periods
%
% Notes:
% Outputs are stored in decision-rule order, i.e. to get variables in order of declaration
......@@ -78,7 +85,6 @@ function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalma
% New output argument aK: 1-step to nk-stpe ahed predictions)
% New input argument nk: max order of predictions in aK
if size(H,2)>1
error('missing_DiffuseKalmanSmootherH3_Z:: H is not a vector. This must not happens')
end
......@@ -102,7 +108,11 @@ Kstar = zeros(mm,pp,smpl);
Kinf = zeros(spstar(1),pp,smpl);
P = zeros(mm,mm,smpl+1);
P1 = P;
PK = zeros(nk,mm,mm,smpl+nk);
if filter_covariance_flag
PK = zeros(nk,mm,mm,smpl+nk);
else
PK = [];
end
Pstar = zeros(spstar(1),spstar(2),smpl);
Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl);
......@@ -120,10 +130,21 @@ QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
if smoother_redux
aalphahat = alphahat;
eetahat = etahat;
else
aalphahat = [];
eetahat = [];
end
epsilonhat = zeros(rr,smpl);
r = zeros(mm,smpl);
if state_uncertainty_flag
V = zeros(mm,mm,smpl);
if smoother_redux
V = zeros(mm+rr,mm+rr,smpl);
else
V = zeros(mm,mm,smpl);
end
N = zeros(mm,mm,smpl);
else
V=[];
......@@ -229,18 +250,40 @@ while notsteady && t<smpl
% p. 157, DK (2012)
end
end
if smoother_redux
ri=zeros(mm,1);
for st=t:-1:max(d+1,t-1)
di = flipud(data_index{st})';
for i = di
if Fi(i,st) > kalman_tol
Li = eye(mm)-Ki(:,i,st)*Z(i,:)/Fi(i,st);
ri = Z(i,:)'/Fi(i,st)*v(i,st)+Li'*ri; % DK (2012), 6.15, equation for r_{t,i-1}
end
end
if st==t-1
aalphahat(:,st) = a1(:,st) + P1(:,:,st)*ri;
else
eetahat(:,st) = QRt*ri;
end
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
end
end
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
a1(:,t+1) = T*a(:,t); %transition according to (6.14) in DK (2012)
P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %transition according to (6.14) in DK (2012)
Pf = P(:,:,t+1);
if filter_covariance_flag
Pf = P(:,:,t+1);
end
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk
if jnk>1
Pf = T*Pf*T' + QQ;
if filter_covariance_flag
if jnk>1
Pf = T*Pf*T' + QQ;
end
PK(jnk,:,:,t+jnk) = Pf;
end