diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 970fc5b20e7a3ade78974fbdf1d02b0d95b7e80d..9bd4a1c63b548b4ad602bd915d53c6b64b0323e7 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -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`.
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 112c1b3f622b6465b85d3bdb1f9d081ca61e304d..8dbee6f4df9a416393ceb1b5e0a4bf09b5dc9a72 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -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
+
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 86ef05dd5166393d0cdd4fcf24fe916724844d43..195622de220271248ca8f03c77de4b4f4a994f71 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -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;
diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index cda2b912336748531dd740b509f18af5913be366..e6c9335ee35ec89b013ed949d9d92fb44aec452f 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -1,6 +1,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
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index b0b9def45cfe0c9ffddd2936ccdeef78a5509ad1..0765eacca74b6ecf080671cf19bf614642f3d2b4 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -1,5 +1,5 @@
-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
- PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
@@ -248,6 +291,9 @@ while notsteady && t<smpl
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
end
P1(:,:,t+1) = P(:,:,t+1);
+
+P1(:,:,t+1)=P(:,:,t+1);
+
% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
% $$$ Fi_s = Fi(:,t);
@@ -307,7 +353,13 @@ while t > d+1
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
if state_uncertainty_flag
N(:,:,t) = Ni; % DK (2012), below 6.15, N_{t-1}=N_{t,0}
- V(:,:,t) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); % KD (2000), eq. (7) with N_{t-1} stored in N(:,:,t)
+ if smoother_redux
+ ptmp = [P1(:,:,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) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); % KD (2000), eq. (7) with N_{t-1} stored in N(:,:,t)
+ end
Ni = T'*Ni*T; % KD (2000), eq. (23), equation for N_{t-1,p_{t-1}}
end
end
@@ -359,10 +411,23 @@ if d
end
etahat(:,t) = QRt*r(:,t); % KD (2000), eq. (27)
if state_uncertainty_flag
- 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)...
- - 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)];
+ ntmp0 = [N_0(:,:,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*ntmp0*pstmp...
+ -(pitmp*ntmp1*pstmp)'...
+ - pitmp*ntmp1*pstmp...
+ - pitmp*ntmp2*Pinf(:,:,t); % DK (2012), eq. 5.30
+
+ else
+ 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)...
+ - Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30
+ end
end
if t > 1
r0(:,t-1) = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
diff --git a/preprocessor b/preprocessor
index b12b254f2c90c92ba8218d07a9860762d41be0fd..1e9ea1acf4e7447f3d336b5c7d571c9a51699723 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit b12b254f2c90c92ba8218d07a9860762d41be0fd
+Subproject commit 1e9ea1acf4e7447f3d336b5c7d571c9a51699723
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b02809da4a63b81a55df087d599396728380d9e5..e76ce7fc1c9269fd9d371f4714b96572e5dd4bee 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -297,6 +297,7 @@ MODFILES = \
kalman_filter_smoother/fs2000a.mod \
kalman_filter_smoother/fs2000_smoother_only.mod \
kalman_filter_smoother/fs2000_smoother_only_ns.mod \
+ kalman_filter_smoother/fs2000_smoother_redux.mod \
kalman_filter_smoother/test_compute_Pinf_Pstar_data.mod \
kalman_filter_smoother/test_compute_Pinf_Pstar.mod \
kalman_filter_smoother/SOE.mod \
diff --git a/tests/kalman_filter_smoother/fs2000_smoother_redux.mod b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod
new file mode 100644
index 0000000000000000000000000000000000000000..72d261a59a2988fa62c0abeede02bd66edacde6f
--- /dev/null
+++ b/tests/kalman_filter_smoother/fs2000_smoother_redux.mod
@@ -0,0 +1,225 @@
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m e_b;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*exp(e_b)*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*exp(e_b)*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+write_latex_dynamic_model;
+
+shocks;
+var e_a; stderr 0.014;
+var e_b; stderr 0.1;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+
+varobs gp_obs gy_obs;
+
+calib_smoother(datafile=fsdat_simul, filtered_vars, filter_step_ahead = [3:4],filter_covariance,smoothed_state_uncertainty) m P c e W R k d n l y dA;
+oo0=oo_;
+
+calib_smoother(datafile=fsdat_simul, filtered_vars, filter_step_ahead = [3:4],filter_covariance,smoothed_state_uncertainty,smoother_redux) m P c e W R k d n l y dA;
+oo1=oo_;
+
+calib_smoother(datafile=fsdat_simul, filtered_vars, filter_step_ahead = [3:4],filter_covariance,smoothed_state_uncertainty,kalman_algo=2,smoother_redux) m P c e W R k d n l y dA;
+oo2=oo_;
+
+
+for k=1:M_.exo_nbr
+ mserr(k)=max(abs(oo0.SmoothedShocks.(M_.exo_names{k})-oo1.SmoothedShocks.(M_.exo_names{k})));
+end
+if max(mserr)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original smoother for shocks!')
+end
+
+vlist = M_.endo_names(oo_.dr.order_var(oo_.dr.restrict_var_list));
+for k=1:length(vlist)
+ merr(k)=max(abs(oo0.SmoothedVariables.(vlist{k})-oo1.SmoothedVariables.(vlist{k})));
+ merrU(k)=max(abs(oo0.UpdatedVariables.(vlist{k})-oo1.UpdatedVariables.(vlist{k})));
+ merrF(k)=max(abs(oo0.FilteredVariables.(vlist{k})-oo1.FilteredVariables.(vlist{k})));
+
+end
+if max(merr)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original smoothed restricted var list!')
+end
+if max(merrU)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original updated restricted var list!')
+end
+if max(merrF)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original filtered restricted var list!')
+end
+
+vlist1 = M_.endo_names(~ismember(M_.endo_names,vlist));
+for k=1:length(vlist1)
+ merr1(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})-oo1.SmoothedVariables.(vlist1{k})));
+ merr1U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})-oo1.UpdatedVariables.(vlist1{k})));
+ merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo1.FilteredVariables.(vlist1{k})));
+end
+if max(merr1)>1.e-12
+ for k=1:length(vlist1)
+ merr2(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})(2:end)-oo1.SmoothedVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original smoothed static variables!')
+ end
+end
+if max(merr1U)>1.e-12
+ for k=1:length(vlist1)
+ merr2U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})(2:end)-oo1.UpdatedVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2U)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original updated static variables!')
+ end
+end
+if max(merr1F)>1.e-12
+ for k=1:length(vlist1)
+ merr2F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})(2:end)-oo1.FilteredVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2F)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original filtered static variables!')
+ end
+end
+merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo1.FilteredVariablesKStepAhead))));
+if max(merrK)>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original k-step ahead forecasts!')
+end
+verrK = max(max(max(max(abs(oo0.FilteredVariablesKStepAheadVariances(:,[1:14 16],[1:14 16],:)-oo1.FilteredVariablesKStepAheadVariances(:,[1:14 16],[1:14 16],:))))));
+if verrK>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original k-step ahead forecast variances!')
+end
+verr=max(max(max(abs(oo0.Smoother.Variance([1:14 16],[1:14 16],:)-oo1.Smoother.Variance([1:14 16],[1:14 16],:)))));
+if verr>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original filter covariance!')
+end
+verrS=max(max(max(abs(oo0.Smoother.State_uncertainty([1:14 16],[1:14 16],:)-oo1.Smoother.State_uncertainty([1:14 16],[1:14 16],:)))));
+if verrS>1.e-12
+ error('smoother_redux with kalman_algo=1 does not replicate original state covariance!')
+end
+
+// now I check kalman_algo=2
+for k=1:M_.exo_nbr
+ mserr(k)=max(abs(oo0.SmoothedShocks.(M_.exo_names{k})-oo2.SmoothedShocks.(M_.exo_names{k})));
+end
+if max(mserr)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original smoother for shocks!')
+end
+
+vlist = M_.endo_names(oo_.dr.order_var(oo_.dr.restrict_var_list));
+for k=1:length(vlist)
+ merr(k)=max(abs(oo0.SmoothedVariables.(vlist{k})-oo2.SmoothedVariables.(vlist{k})));
+ merrU(k)=max(abs(oo0.UpdatedVariables.(vlist{k})-oo2.UpdatedVariables.(vlist{k})));
+ merrF(k)=max(abs(oo0.FilteredVariables.(vlist{k})-oo2.FilteredVariables.(vlist{k})));
+
+end
+if max(merr)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original smoothed restricted var list!')
+end
+if max(merrU)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original updated restricted var list!')
+end
+if max(merrF)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original filtered restricted var list!')
+end
+
+vlist1 = M_.endo_names(~ismember(M_.endo_names,vlist));
+for k=1:length(vlist1)
+ merr1(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})-oo2.SmoothedVariables.(vlist1{k})));
+ merr1U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})-oo2.UpdatedVariables.(vlist1{k})));
+ merr1F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})-oo2.FilteredVariables.(vlist1{k})));
+end
+if max(merr1)>1.e-12
+ for k=1:length(vlist1)
+ merr2(k)=max(abs(oo0.SmoothedVariables.(vlist1{k})(2:end)-oo2.SmoothedVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original smoothed static variables!')
+ end
+end
+if max(merr1U)>1.e-12
+ for k=1:length(vlist1)
+ merr2U(k)=max(abs(oo0.UpdatedVariables.(vlist1{k})(2:end)-oo2.UpdatedVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2U)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original updated static variables!')
+ end
+end
+if max(merr1F)>1.e-12
+ for k=1:length(vlist1)
+ merr2F(k)=max(abs(oo0.FilteredVariables.(vlist1{k})(2:end)-oo2.FilteredVariables.(vlist1{k})(2:end)));
+ end
+ if max(merr2F)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original filtered static variables!')
+ end
+end
+merrK = max(max(max(abs(oo0.FilteredVariablesKStepAhead-oo2.FilteredVariablesKStepAhead))));
+if max(merrK)>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original k-step ahead forecasts!')
+end
+verrK = max(max(max(max(abs(oo0.FilteredVariablesKStepAheadVariances(:,[1:14 16],[1:14 16],:)-oo2.FilteredVariablesKStepAheadVariances(:,[1:14 16],[1:14 16],:))))));
+if verrK>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original k-step ahead forecast variances!')
+end
+verr=max(max(max(abs(oo0.Smoother.Variance([1:14 16],[1:14 16],:)-oo2.Smoother.Variance([1:14 16],[1:14 16],:)))));
+if verr>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original filter covariance!')
+end
+verrS=max(max(max(abs(oo0.Smoother.State_uncertainty([1:14 16],[1:14 16],:)-oo2.Smoother.State_uncertainty([1:14 16],[1:14 16],:)))));
+if verrS>1.e-12
+ error('smoother_redux with kalman_algo=2 does not replicate original state covariance!')
+end