From 5fa6265944426ee5c50531eb9d59c2c11fbfe707 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Tue, 4 May 2021 09:01:27 +0200
Subject: [PATCH] Implement new option smoother_redux, to allow fast smoother
for very large models. It runs smoother only for state variables and computes
the remaining ones ex-post
Contains improvements, in order to recover as much as possible static unobserved (filtered, smoothed, updated, k-step ahead), Variance, State_uncertainty, k-step ahead variances trying to map lagged states onto current ones using pinv(T). This has exceptions (namely lagged shocks which are ONLY used to recover static NON observed variables). this exception is also trapped.
For such extensions we can only recover smoothed variables starting from d+1. Variances CANNOT be recovered for such variables (the smoother gives ZERO.)
---
doc/manual/source/the-model-file.rst | 32 +++
matlab/DsgeSmoother.m | 159 +++++++++++--
matlab/default_option_values.m | 1 +
matlab/missing_DiffuseKalmanSmootherH1_Z.m | 131 +++++++---
matlab/missing_DiffuseKalmanSmootherH3_Z.m | 105 ++++++--
preprocessor | 2 +-
tests/Makefile.am | 1 +
.../fs2000_smoother_redux.mod | 225 ++++++++++++++++++
8 files changed, 584 insertions(+), 72 deletions(-)
create mode 100644 tests/kalman_filter_smoother/fs2000_smoother_redux.mod
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 970fc5b20e..9bd4a1c63b 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 112c1b3f62..8dbee6f4df 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 86ef05dd51..195622de22 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 cda2b91233..e6c9335ee3 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 b0b9def45c..0765eacca7 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 b12b254f2c..1e9ea1acf4 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit b12b254f2c90c92ba8218d07a9860762d41be0fd
+Subproject commit 1e9ea1acf4e7447f3d336b5c7d571c9a51699723
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b02809da4a..e76ce7fc1c 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 0000000000..72d261a59a
--- /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
--
GitLab