diff --git a/matlab/DiffuseKalmanSmoother1.m b/matlab/DiffuseKalmanSmoother1.m
index 8526bb096acdae132e6c8dcfd868437ba517ba48..7f54bb63cc4c57fb2a66daa2c0d105cd2ef0dfb6 100644
--- a/matlab/DiffuseKalmanSmoother1.m
+++ b/matlab/DiffuseKalmanSmoother1.m
@@ -17,10 +17,10 @@ function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar
% mf: observed variables index in the state vector
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% etahat: smoothed shocks
-% atilde: matrix of filtered variables (E_t Y_t)
-% aK: 3D array of k step ahead filtered state variables
+% atilde: matrix of updated variables (a_{t|t})
+% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
% SPECIAL REQUIREMENTS
% See "Filtering and Smoothing of State Vector for Diffuse State Space
diff --git a/matlab/DiffuseKalmanSmoother1_Z.m b/matlab/DiffuseKalmanSmoother1_Z.m
index cfb33eaac7b3760c851344a8df061f1f33c9da2e..754dbeb8f8c4c02ac7524bcdb75cc78cae6b239c 100644
--- a/matlab/DiffuseKalmanSmoother1_Z.m
+++ b/matlab/DiffuseKalmanSmoother1_Z.m
@@ -18,7 +18,7 @@ function [alphahat,etahat,atilde,P,aK,PK,d,decomp] = DiffuseKalmanSmoother1_Z(T,
% OUTPUTS
% alphahat: smoothed variables (a_{t|T})
% etahat: smoothed shocks
-% atilde: matrix of filtered variables (a_{t|t})
+% 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
diff --git a/matlab/DiffuseKalmanSmoother3.m b/matlab/DiffuseKalmanSmoother3.m
index f5d65cfcf903ac033c9af8feec13f69fca78ddbf..be2c22235b88821703b93dac2281028f6169e4f9 100644
--- a/matlab/DiffuseKalmanSmoother3.m
+++ b/matlab/DiffuseKalmanSmoother3.m
@@ -17,10 +17,10 @@ function [alphahat,etahat,a, aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf1,Pstar1,Y,t
% mf: observed variables index in the state vector
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% etahat: smoothed shocks
-% a: matrix of filtered variables
-% aK: 3D array of k step ahead filtered state variables
+% a: matrix of updated variables (a_{t|t})
+% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
%
% SPECIAL REQUIREMENTS
% See "Filtering and Smoothing of State Vector for Diffuse State Space
@@ -296,4 +296,4 @@ else
etahat(:,1) = QRt*r(:,1);
end
-a=a(:,1:end-1);
+
diff --git a/matlab/DiffuseKalmanSmoother3_Z.m b/matlab/DiffuseKalmanSmoother3_Z.m
index 170ab24467c9a8f1ec7057712d7926124e7efa63..72679dbdb7ac6358d88ce680ec061d11d11ce51f 100644
--- a/matlab/DiffuseKalmanSmoother3_Z.m
+++ b/matlab/DiffuseKalmanSmoother3_Z.m
@@ -16,10 +16,10 @@ function [alphahat,etahat,a,P,aK,PK,d,decomp] = DiffuseKalmanSmoother3_Z(T,Z,R,Q
% smpl: sample size
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% etahat: smoothed shocks
-% a: matrix of filtered variables
-% aK: 3D array of k step ahead filtered state variables
+% a: 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
@@ -300,8 +300,6 @@ else
etahat(:,1) = QRt*r(:,1);
end
-a=a(:,1:end-1);
-
if nargout > 7
decomp = zeros(nk,mm,rr,smpl+nk);
ZRQinv = inv(Z*QQ*Z');
diff --git a/matlab/DiffuseKalmanSmootherH1.m b/matlab/DiffuseKalmanSmootherH1.m
index 2150c22d9ffaae1fc767d3972a2f8784e607dc2f..86dcfbc01757a8e4c012c5eea81f309284dca9ee 100644
--- a/matlab/DiffuseKalmanSmootherH1.m
+++ b/matlab/DiffuseKalmanSmootherH1.m
@@ -1,11 +1,12 @@
-function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
-% function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
+gfunction [alphahat,epsilonhat,etahat,atilde, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
+% function [alphahat,epsilonhat,etahat,atilde, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
% Computes the diffuse kalman smoother with measurement error, in the case of a non-singular var-cov matrix
%
% INPUTS
% T: mm*mm matrix
% R: mm*rr matrix
% Q: rr*rr matrix
+% H: pp*pp matrix variance of measurement errors
% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
% Pstar1: mm*mm variance-covariance matrix with stationary variables
% Y: pp*1 vector
@@ -16,11 +17,11 @@ function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pi
% mf: observed variables index in the state vector
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% epsilonhat:smoothed measurement errors
-% etahat: smoothed shocks
-% a: matrix of one step ahead filtered state variables
-% aK: 3D array of k step ahead filtered state variables
+% 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)}
%
% SPECIAL REQUIREMENTS
% See "Filtering and Smoothing of State Vector for Diffuse State Space
@@ -57,6 +58,7 @@ spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
+atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+nk);
iF = zeros(pp,pp,smpl);
Fstar = zeros(pp,pp,smpl);
@@ -92,8 +94,10 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
return
end
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
- Kinf(:,:,t) = T*Pinf(:,mf,t)*iFinf(:,:,t);
- a(:,t+1) = T*a(:,t) + Kinf(:,:,t)*v(:,t);
+ PZI = Pinf(:,mf,t)*iFinf(:,:,t);
+ atilde(:,t) = a(:,t) + PZI*v(:,t);
+ Kinf(:,:,t) = T*PZI;
+ a(:,t+1) = T*atilde(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
@@ -120,29 +124,33 @@ while notsteady & t<smpl
return
end
iF(:,:,t) = inv(P(mf,mf,t) + H);
- K(:,:,t) = T*P(:,mf,t)*iF(:,:,t);
+ PZI = P(:,mf,t)*iF(:,:,t);
+ atilde(:,t) = a(:,t) + PZI*v(:,t);
+ K(:,:,t) = T*PZI;
L(:,:,t) = T-K(:,:,t)*Z;
- a(:,t+1) = T*a(:,t) + K(:,:,t)*v(:,t);
+ a(:,t+1) = T*atilde(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
end
-K_s = K(:,:,t);
-iF_s = iF(:,:,t);
-P_s = P(:,:,t+1);
if t<smpl
- t_steady = t+1;
- P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
- iF = cat(3,iF(:,:,1:t),repmat(inv(P_s(mf,mf)+H),[1 1 smpl-t_steady+1]));
- L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t_steady+1]));
- K = cat(3,K(:,:,1:t),repmat(T*P_s(:,mf)*iF_s,[1 1 smpl-t_steady+1]));
+ PZI_s = PZI;
+ K_s = K(:,:,t);
+ iF_s = iF(:,:,t);
+ P_s = P(:,:,t+1);
+ t_steady = t+1;
+ P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
+ iF = cat(3,iF(:,:,1:t),repmat(inv(P_s(mf,mf)+H),[1 1 smpl-t_steady+1]));
+ L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t_steady+1]));
+ K = cat(3,K(:,:,1:t),repmat(T*P_s(:,mf)*iF_s,[1 1 smpl-t_steady+1]));
end
while t<smpl
t=t+1;
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- a(:,t+1) = T*a(:,t) + K_s*v(:,t);
+ atilde(:,t) = a(:,t) + PZI_s*v(:,t);
+ a(:,t+1) = T*atilde(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
@@ -172,4 +180,4 @@ else
alphahat(:,1) = a(:,1) + P(:,:,1)*r0;
etahat(:,1) = QRt*r(:,1);
end
-epsilonhat = Y-alphahat(mf,:)-trend;
\ No newline at end of file
+epsilonhat = Y-alphahat(mf,:)-trend;
diff --git a/matlab/DiffuseKalmanSmootherH1_Z.m b/matlab/DiffuseKalmanSmootherH1_Z.m
new file mode 100644
index 0000000000000000000000000000000000000000..398f6a3824506dfdea0caa349158330c466b3265
--- /dev/null
+++ b/matlab/DiffuseKalmanSmootherH1_Z.m
@@ -0,0 +1,220 @@
+function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,d,decomp] = DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
+
+% function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
+% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix
+%
+% INPUTS
+% T: mm*mm matrix
+% Z: pp*mm matrix
+% R: mm*rr matrix
+% Q: rr*rr matrix
+% H: pp*pp matrix variance of measurement errors
+% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
+% Pstar1: mm*mm variance-covariance matrix with stationary variables
+% Y: pp*1 vector
+% pp: number of observed variables
+% mm: number of state variables
+% smpl: sample size
+%
+% 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)
+% d: number of periods where filter remains in diffuse part
+% (should be equal to the order of integration of the model)
+% decomp: decomposition of the effect of shocks on filtered values
+%
+% SPECIAL REQUIREMENTS
+% See "Filtering and Smoothing of State Vector for Diffuse State Space
+% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
+% Analysis, vol. 24(1), pp. 85-98).
+
+% Copyright (C) 2004-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+
+% modified by M. Ratto:
+% new output argument aK (1-step to k-step predictions)
+% new options_.nk: the max step ahed prediction in aK (default is 4)
+% new crit1 value for rank of Pinf
+% it is assured that P is symmetric
+
+
+global options_
+
+d = 0;
+decomp = [];
+nk = options_.nk;
+spinf = size(Pinf1);
+spstar = size(Pstar1);
+v = zeros(pp,smpl);
+a = zeros(mm,smpl+1);
+atilde = zeros(mm,smpl);
+aK = zeros(nk,mm,smpl+nk);
+PK = zeros(nk,mm,mm,smpl+nk);
+iF = zeros(pp,pp,smpl);
+Fstar = zeros(pp,pp,smpl);
+iFinf = zeros(pp,pp,smpl);
+K = zeros(mm,pp,smpl);
+L = zeros(mm,mm,smpl);
+Linf = zeros(mm,mm,smpl);
+Kstar = zeros(mm,pp,smpl);
+P = zeros(mm,mm,smpl+1);
+Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
+Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
+crit = options_.kalman_tol;
+crit1 = 1.e-8;
+steady = smpl;
+rr = size(Q,1);
+QQ = R*Q*transpose(R);
+QRt = Q*transpose(R);
+alphahat = zeros(mm,smpl);
+etahat = zeros(rr,smpl);
+epsilonhat = zeros(size(Y));
+r = zeros(mm,smpl);
+
+t = 0;
+while rank(Pinf(:,:,t+1),crit1) & t<smpl
+ t = t+1;
+ v(:,t)= Y(:,t) - Z*a(:,t);
+ F = Z*Pinf(:,:,t)*Z';
+ if rcond(F) < crit
+ return
+ end
+ iFinf(:,:,t) = inv(F);
+ PZI = Pinf(:,:,t)*Z'*iFinf(:,:,t);
+ atilde(:,t) = a(:,t) + PZI*v(:,t);
+ Kinf(:,:,t) = T*PZI;
+ a(:,t+1) = T*atilde(:,t);
+ aK(1,:,t+1) = a(:,t+1);
+ % isn't a meaningless as long as we are in the diffuse part? MJ
+ for jnk=2:nk,
+ aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
+ end
+ Linf(:,:,t) = T - Kinf(:,:,t)*Z;
+ Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
+ Kstar(:,:,t) = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
+ Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*F*Kstar(:,:,t)' + QQ;
+ Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
+end
+d = t;
+P(:,:,d+1) = Pstar(:,:,d+1);
+iFinf = iFinf(:,:,1:d);
+Linf = Linf(:,:,1:d);
+Fstar = Fstar(:,:,1:d);
+Kstar = Kstar(:,:,1:d);
+Pstar = Pstar(:,:,1:d);
+Pinf = Pinf(:,:,1:d);
+notsteady = 1;
+while notsteady & t<smpl
+ t = t+1;
+ v(:,t) = Y(:,t) - Z*a(:,t);
+ P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
+ F = Z*P(:,:,t)*Z' + H;
+ if rcond(F) < crit
+ return
+ end
+ iF(:,:,t) = inv(F);
+ PZI = P(:,:,t)*Z'*iF(:,:,t);
+ atilde(:,t) = a(:,t) + PZI*v(:,t);
+ K(:,:,t) = T*PZI;
+ L(:,:,t) = T-K(:,:,t)*Z;
+ a(:,t+1) = T*atilde(:,t);
+ Pf = P(:,:,t);
+ for jnk=1:nk,
+ Pf = T*Pf*T' + QQ;
+ aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
+ PK(jnk,:,:,t+jnk) = Pf;
+ end
+ P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*Z'*K(:,:,t)' + QQ;
+ notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
+end
+if t<smpl
+ PZI_s = PZI;
+ K_s = K(:,:,t);
+ iF_s = iF(:,:,t);
+ P_s = P(:,:,t+1);
+ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
+ iF = cat(3,iF(:,:,1:t),repmat(iF_s,[1 1 smpl-t]));
+ L = cat(3,L(:,:,1:t),repmat(T-K_s*Z,[1 1 smpl-t]));
+ K = cat(3,K(:,:,1:t),repmat(T*P_s*Z'*iF_s,[1 1 smpl-t]));
+end
+while t<smpl
+ t=t+1;
+ v(:,t) = Y(:,t) - Z*a(:,t);
+ atilde(:,t) = a(:,t) + PZI*v(:,t);
+ a(:,t+1) = T*atilde(:,t);
+ Pf = P(:,:,t);
+ for jnk=1:nk,
+ Pf = T*Pf*T' + QQ;
+ aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
+ PK(jnk,:,:,t+jnk) = Pf;
+ end
+end
+t = smpl+1;
+while t>d+1 & t>2
+ t = t-1;
+ r(:,t-1) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t);
+ alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t-1);
+ etahat(:,t) = QRt*r(:,t);
+end
+if d
+ r0 = zeros(mm,d);
+ r0(:,d) = r(:,d);
+ r1 = zeros(mm,d);
+ for t = d:-1:2
+ r0(:,t-1) = Linf(:,:,t)'*r0(:,t);
+ r1(:,t-1) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t)) + Linf(:,:,t)'*r1(:,t);
+ alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t-1) + Pinf(:,:,t)*r1(:,t-1);
+ etahat(:,t) = QRt*r0(:,t);
+ end
+ r0_0 = Linf(:,:,1)'*r0(:,1);
+ r1_0 = Z'*(iFinf(:,:,1)*v(:,1)-Kstar(:,:,1)'*r0(:,1)) + Linf(:,:,1)'*r1(:,1);
+ alphahat(:,1) = a(:,1) + Pstar(:,:,1)*r0_0 + Pinf(:,:,1)*r1_0;
+ etahat(:,1) = QRt*r0(:,1);
+else
+ r0 = Z'*iF(:,:,1)*v(:,1) + L(:,:,1)'*r(:,1);
+ alphahat(:,1) = a(:,1) + P(:,:,1)*r0;
+ etahat(:,1) = QRt*r(:,1);
+end
+
+if nargout > 7
+ decomp = zeros(nk,mm,rr,smpl+nk);
+ ZRQinv = inv(Z*QQ*Z');
+ for t = d:smpl
+ ri_d = Z'*iF(:,:,t)*v(:,t);
+
+ % calculate eta_tm1t
+ eta_tm1t = QRt*ri_d;
+ % calculate decomposition
+ Ttok = eye(mm,mm);
+ for h = 1:nk
+ for j=1:rr
+ eta=zeros(rr,1);
+ eta(j) = eta_tm1t(j);
+ decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta;
+ end
+ end
+ end
+end
+epsilonhat = Y-alphahat(mf,:)-trend;
diff --git a/matlab/DiffuseKalmanSmootherH3.m b/matlab/DiffuseKalmanSmootherH3.m
index 7a69f56fb3095514cb0fd3aee043452bc826fc1c..926b8f2ce659f01a91b5409ee6b4a2c59815e87e 100644
--- a/matlab/DiffuseKalmanSmootherH3.m
+++ b/matlab/DiffuseKalmanSmootherH3.m
@@ -1,4 +1,4 @@
-function [alphahat,epsilonhat,etahat,a1, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
+function [alphahat,epsilonhat,etahat,a, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
% function [alphahat,epsilonhat,etahat,a1, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
% Computes the diffuse kalman smoother with measurement error, in the case of a singular var-cov matrix.
% Univariate treatment of multivariate time series.
@@ -17,11 +17,11 @@ function [alphahat,epsilonhat,etahat,a1, aK] = DiffuseKalmanSmootherH3(T,R,Q,H,P
% mf: observed variables index in the state vector
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% epsilonhat:smoothed measurement error
% etahat: smoothed shocks
-% a1: matrix of one step ahead filtered state variables
-% aK: 3D array of k step ahead filtered state variables
+% a: matrix of updated variables (a_{t|t})
+% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
%
% SPECIAL REQUIREMENTS
% See "Filtering and Smoothing of State Vector for Diffuse State Space
@@ -63,8 +63,8 @@ nk = options_.nk;
spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
-a = zeros(mm,smpl+1);
-a1 = a;
+a = zeros(mm,smpl);
+a1 = zeros(mm,smpl+1);
aK = zeros(nk,mm,smpl+nk);
if isempty(options_.diffuse_d),
@@ -107,7 +107,7 @@ icc=0;
newRank = rank(Pinf(:,:,1),crit1);
while newRank & t < smpl
t = t+1;
- a1(:,t) = a(:,t);
+ a(:,t) = a1(:,t);
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
Pstar1(:,:,t) = Pstar(:,:,t);
@@ -163,7 +163,7 @@ while newRank & t < smpl
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
end
@@ -191,7 +191,7 @@ Pinf1 = Pinf1(:,:,1:d);
notsteady = 1;
while notsteady & t<smpl
t = t+1;
- a1(:,t) = a(:,t);
+ a(:,t) = a1(:,t);
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
P1(:,:,t) = P(:,:,t);
for i=1:pp
@@ -205,7 +205,7 @@ while notsteady & t<smpl
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
end
@@ -225,19 +225,18 @@ if t<smpl
end
while t<smpl
t=t+1;
- a1(:,t) = a(:,t);
+ a(:,t) = a1(:,t);
for i=1:pp
v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
if Fi_s(i) > crit
a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
end
end
-a1(:,t+1) = a(:,t+1);
ri=r;
t = smpl+1;
while t>d+1 & t>2,
@@ -300,4 +299,4 @@ else
end
epsilonhat = Y-alphahat(mf,:)-trend;
-a=a(:,1:end-1);
+
diff --git a/matlab/DiffuseKalmanSmootherH3_Z.m b/matlab/DiffuseKalmanSmootherH3_Z.m
new file mode 100644
index 0000000000000000000000000000000000000000..5775c8ebc63502eb0d81d766c78492acb1fdf682
--- /dev/null
+++ b/matlab/DiffuseKalmanSmootherH3_Z.m
@@ -0,0 +1,331 @@
+function [alphahat,epsilonhat,etahat,a,P,aK,PK,d,decomp] = DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
+% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp_filt] = DiffuseKalmanSmootherH3(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl)
+% Computes the diffuse kalman smoother without measurement error, in the case of a singular var-cov matrix.
+% Univariate treatment of multivariate time series.
+%
+% INPUTS
+% T: mm*mm matrix
+% Z: pp*mm matrix
+% R: mm*rr matrix
+% Q: rr*rr matrix
+% Pinf1: mm*mm diagonal matrix with with q ones and m-q zeros
+% Pstar1: mm*mm variance-covariance matrix with stationary variables
+% Y: pp*1 vector
+% pp: number of observed variables
+% mm: number of state variables
+% smpl: sample size
+%
+% OUTPUTS
+% alphahat: smoothed state variables (a_{t|T})
+% etahat: smoothed shocks
+% epsilonhat:smoothed measurement error
+% a: 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)
+% d: number of periods where filter remains in diffuse part
+% (should be equal to the order of integration of the model)
+% decomp: decomposition of the effect of shocks on filtered values
+%
+% SPECIAL REQUIREMENTS
+% See "Filtering and Smoothing of State Vector for Diffuse State Space
+% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
+% Analysis, vol. 24(1), pp. 85-98).
+
+% Copyright (C) 2004-2008 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+
+% Modified by M. Ratto
+% New output argument aK: 1-step to nk-stpe ahed predictions)
+% New input argument nk: max order of predictions in aK
+% New option options_.diffuse_d where the DKF stops (common with
+% diffuselikelihood3)
+% New icc variable to count number of iterations for Finf steps
+% Pstar % Pinf simmetric
+% New termination of DKF iterations based on options_.diffuse_d
+% Li also stored during DKF iterations !!
+% some bugs corrected in the DKF part of the smoother (Z matrix and
+% alphahat)
+
+global options_
+
+d = 0;
+decomp = [];
+nk = options_.nk;
+spinf = size(Pinf1);
+spstar = size(Pstar1);
+v = zeros(pp,smpl);
+a = zeros(mm,smpl);
+a1 = zeros(mm,smpl+1);
+aK = zeros(nk,mm,smpl+nk);
+
+if isempty(options_.diffuse_d),
+ smpl_diff = 1;
+else
+ smpl_diff=rank(Pinf1);
+end
+
+Fstar = zeros(pp,smpl_diff);
+Finf = zeros(pp,smpl_diff);
+Ki = zeros(mm,pp,smpl);
+Li = zeros(mm,mm,pp,smpl);
+Linf = zeros(mm,mm,pp,smpl_diff);
+L0 = zeros(mm,mm,pp,smpl_diff);
+Kstar = zeros(mm,pp,smpl_diff);
+P = zeros(mm,mm,smpl+1);
+P1 = P;
+aK = zeros(nk,mm,smpl+nk);
+PK = zeros(nk,mm,mm,smpl+nk);
+Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
+Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
+Pstar1 = Pstar;
+Pinf1 = Pinf;
+crit = options_.kalman_tol;
+crit1 = 1.e-6;
+steady = smpl;
+rr = size(Q,1); % number of structural shocks
+QQ = R*Q*transpose(R);
+QRt = Q*transpose(R);
+alphahat = zeros(mm,smpl);
+etahat = zeros(rr,smpl);
+epsilonhat = zeros(size(Y));
+r = zeros(mm,smpl);
+
+t = 0;
+icc=0;
+newRank = rank(Pinf(:,:,1),crit1);
+while newRank & t < smpl
+ t = t+1;
+ a(:,t) = a1(:,t);
+ Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
+ Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
+ Pstar1(:,:,t) = Pstar(:,:,t);
+ Pinf1(:,:,t) = Pinf(:,:,t);
+ for i=1:pp
+ Zi = Z(i,:);
+ v(i,t) = Y(i,t)-Zi*a(:,t);
+ Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i,i);
+ Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
+ Kstar(:,i,t) = Pstar(:,:,t)*Zi';
+ if Finf(i,t) > crit & newRank
+ icc=icc+1;
+ Kinf(:,i,t) = Pinf(:,:,t)*Zi';
+ Linf(:,:,i,t) = eye(mm) - Kinf(:,i,t)*Z(i,:)/Finf(i,t);
+ L0(:,:,i,t) = (Kinf(:,i,t)*Fstar(i,t)/Finf(i,t) - Kstar(:,i,t))*Zi/Finf(i,t);
+ a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
+ Pstar(:,:,t) = Pstar(:,:,t) + ...
+ Kinf(:,i,t)*Kinf(:,i,t)'*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
+ (Kstar(:,i,t)*Kinf(:,i,t)' +...
+ Kinf(:,i,t)*Kstar(:,i,t)')/Finf(i,t);
+ Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
+ Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
+ Pinf(:,:,t)=tril(Pinf(:,:,t))+tril(Pinf(:,:,t),-1)';
+ % new terminiation criteria by M. Ratto
+ P0=Pinf(:,:,t);
+ if ~isempty(options_.diffuse_d),
+ newRank = (icc<options_.diffuse_d);
+ if newRank & (any(diag(Z*P0*Z')>crit)==0 & rank(P0,crit1)==0);
+ disp('WARNING!! Change in OPTIONS_.DIFFUSE_D in univariate DKF')
+ options_.diffuse_d = icc;
+ newRank=0;
+ end
+ else
+ newRank = (any(diag(Z*P0*Z')>crit) | rank(P0,crit1));
+ if newRank==0,
+ options_.diffuse_d = icc;
+ end
+ end,
+ % end new terminiation criteria by M. Ratto
+ elseif Fstar(i,t) > crit
+ %% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
+ %% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
+ %% rank(Pinf)=0. [st�phane,11-03-2004].
+ Li(:,:,i,t) = eye(mm)-Kstar(:,i,t)*Z(i,:)/Fstar(i,t); % we need to store Li for DKF smoother
+ a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
+ Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
+ Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)';
+ end
+ end
+ a1(:,t+1) = T*a(:,t);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ end
+ Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
+ Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
+ P0=Pinf(:,:,t+1);
+ if newRank,
+ newRank = rank(P0,crit1);
+ end
+end
+
+
+d = t;
+P(:,:,d+1) = Pstar(:,:,d+1);
+Linf = Linf(:,:,:,1:d);
+L0 = L0(:,:,:,1:d);
+Fstar = Fstar(:,1:d);
+Finf = Finf(:,1:d);
+Kstar = Kstar(:,:,1:d);
+Pstar = Pstar(:,:,1:d);
+Pinf = Pinf(:,:,1:d);
+Pstar1 = Pstar1(:,:,1:d);
+Pinf1 = Pinf1(:,:,1:d);
+notsteady = 1;
+while notsteady & t<smpl
+ t = t+1;
+ a(:,t) = a1(:,t);
+ P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
+ P1(:,:,t) = P(:,:,t);
+ for i=1:pp
+ Zi = Z(i,:);
+ v(i,t) = Y(i,t) - Zi*a(:,t);
+ Fi(i,t) = Zi*P(:,:,t)*Zi' + H(i,i);
+ Ki(:,i,t) = P(:,:,t)*Zi';
+ if Fi(i,t) > crit
+ Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
+ a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
+ P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
+ P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
+ end
+ end
+ a1(:,t+1) = T*a(:,t);
+ Pf = P(:,:,t);
+ for jnk=1:nk,
+ Pf = T*Pf*T' + QQ;
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ PK(jnk,:,:,t+jnk) = Pf;
+ end
+ P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
+ notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
+end
+P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
+P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
+Fi_s = Fi(:,t);
+Ki_s = Ki(:,:,t);
+L_s =Li(:,:,:,t);
+if t<smpl
+ P = cat(3,P(:,:,1:t),repmat(P_s,[1 1 smpl-t]));
+ P1 = cat(3,P1(:,:,1:t),repmat(P1_s,[1 1 smpl-t]));
+ Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t]));
+ Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t]));
+ Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t]));
+end
+while t<smpl
+ t=t+1;
+ a(:,t) = a1(:,t);
+ for i=1:pp
+ Zi = Z(i,:);
+ v(i,t) = Y(i,t) - Zi*a(:,t);
+ if Fi_s(i) > crit
+ a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
+ end
+ end
+ a1(:,t+1) = T*a(:,t);
+ Pf = P(:,:,t);
+ for jnk=1:nk,
+ Pf = T*Pf*T' + QQ;
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ PK(jnk,:,:,t+jnk) = Pf;
+ end
+end
+ri=r;
+t = smpl+1;
+while t>d+1 & t>2,
+ t = t-1;
+ for i=pp:-1:1
+ if Fi(i,t) > crit
+ ri(:,t) = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri(:,t);
+ end
+ end
+ r(:,t-1) = ri(:,t);
+ alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t-1);
+ etahat(:,t) = QRt*r(:,t);
+ ri(:,t-1) = T'*ri(:,t);
+end
+if d
+ r0 = zeros(mm,d); r0(:,d) = ri(:,d);
+ r1 = zeros(mm,d);
+ for t = d:-1:2
+ for i=pp:-1:1
+ if Finf(i,t) > crit & ~(t==d & i>options_.diffuse_d), % use of options_.diffuse_d to be sure of DKF termination
+ r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
+ L0(:,:,i,t)'*r0(:,t) + Linf(:,:,i,t)'*r1(:,t);
+ r0(:,t) = Linf(:,:,i,t)'*r0(:,t);
+ elseif Fstar(i,t) > crit % step needed whe Finf == 0
+ r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+Li(:,:,i,t)'*r0(:,t);
+ end
+ end
+ alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
+ r(:,t-1) = r0(:,t);
+ etahat(:,t) = QRt*r(:,t);
+ r0(:,t-1) = T'*r0(:,t);
+ r1(:,t-1) = T'*r1(:,t);
+ end
+ r0_0 = r0(:,1);
+ r1_0 = r1(:,1);
+ for i=pp:-1:1
+ if Finf(i,1) > crit,
+ r1_0 = Z(i,:)'*v(i,1)/Finf(i,1) + ...
+ L0(:,:,i,1)'*r0_0 + Linf(:,:,i,1)'*r1_0;
+ r0_0 = Linf(:,:,i,1)'*r0_0;
+ elseif Fstar(i,1) > crit, % step needed when Finf=0
+ r0_0=Z(i,:)'/Fstar(i,1)*v(i,1)+Li(:,:,i,1)'*r0_0;
+ end
+ end
+ alphahat(:,1) = a1(:,1) + Pstar1(:,:,1)*r0_0 + Pinf1(:,:,1)*r1_0;
+ etahat(:,1) = QRt*r(:,1);
+else
+ r0 = ri(:,1);
+ for i=pp:-1:1
+ if Fi(i,1) > crit
+ r0 = Z(i,:)'/Fi(i,1)*v(i,1)+Li(:,:,i,1)'*r0;
+ end
+ end
+ alphahat(:,1) = a1(:,1) + P1(:,:,1)*r0;
+ etahat(:,1) = QRt*r(:,1);
+end
+
+if nargout > 7
+ decomp = zeros(nk,mm,rr,smpl+nk);
+ ZRQinv = inv(Z*QQ*Z');
+ for t = d:smpl
+ ri_d = zeros(mm,1);
+ for i=pp:-1:1
+ if Fi(i,t) > crit
+ ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+Li(:,:,i,t)'*ri_d;
+ end
+ end
+
+ % calculate eta_tm1t
+ eta_tm1t = QRt*ri_d;
+ % calculate decomposition
+ Ttok = eye(mm,mm);
+ for h = 1:nk
+ for j=1:rr
+ eta=zeros(rr,1);
+ eta(j) = eta_tm1t(j);
+ decomp(h,:,j,t+h) = Ttok*P1(:,:,t)*Z'*ZRQinv*Z*R*eta;
+ end
+ Ttok = T*Ttok;
+ end
+ end
+end
+
+epsilonhat = Y-alphahat(mf,:)-trend;
diff --git a/matlab/DiffuseKalmanSmootherH3corr.m b/matlab/DiffuseKalmanSmootherH3corr.m
index fb7b060cc06b932e06809418d089337a029572a8..438a0da9b362531601dc72ddf5644df07d32b3d2 100644
--- a/matlab/DiffuseKalmanSmootherH3corr.m
+++ b/matlab/DiffuseKalmanSmootherH3corr.m
@@ -1,4 +1,4 @@
-function [alphahat,epsilonhat,etahat,a1] = DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
+function [alphahat,epsilonhat,etahat,a,aK] = DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
% function [alphahat,epsilonhat,etahat,a1] = DiffuseKalmanSmootherH3corr(T,R,Q,H,Pinf1,Pstar1,Y,trend,pp,mm,smpl,mf)
% Computes the diffuse kalman smoother with measurement error, in the case of a singular var-cov matrix.
@@ -18,10 +18,11 @@ function [alphahat,epsilonhat,etahat,a1] = DiffuseKalmanSmootherH3corr(T,R,Q,H,P
% mf: observed variables index in the state vector
%
% OUTPUTS
-% alphahat: smoothed state variables
+% alphahat: smoothed state variables (a_{t|T})
% epsilonhat:smoothed measurement error
% etahat: smoothed shocks
-% a1: matrix of one step ahead filtered state variables
+% a: matrix of updated variables (a_{t|t})
+% aK: matrix of one step ahead filtered state variables (a_{t+k|t})
% SPECIAL REQUIREMENTS
% See "Fast Filtering and Smoothing for Multivariate State Space
@@ -47,6 +48,7 @@ function [alphahat,epsilonhat,etahat,a1] = DiffuseKalmanSmootherH3corr(T,R,Q,H,P
global options_;
+nk = options_.nk;
rr = size(Q,1);
T = cat(1,cat(2,T,zeros(mm,pp)),zeros(pp,mm+pp));
R = cat(1,cat(2,R,zeros(mm,pp)),cat(2,zeros(pp,rr),eye(pp)));
@@ -62,11 +64,12 @@ Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
v = zeros(pp,smpl);
-a = zeros(mm+pp,smpl+1);
-a1 = a;
+a = zeros(mm+pp,smpl);
+a1 = zeros(mm+pp,smpl+1);
+aK = zeros(nk,mm,smpl+nk);
Fstar = zeros(pp,smpl);
Finf = zeros(pp,smpl);
-Fi = zeros(pp,smpl);
+Fi = zeros(pp,smpl);
Ki = zeros(mm+pp,pp,smpl);
Li = zeros(mm+pp,mm+pp,pp,smpl);
Linf = zeros(mm+pp,mm+pp,pp,smpl);
@@ -93,7 +96,7 @@ t = 0;
newRank = rank(Pinf(:,:,1),crit);
while newRank & t < smpl
t = t+1;
- a1(:,t) = a(:,t);
+ a(:,t) = a1(:,t);
Pstar1(:,:,t) = Pstar(:,:,t);
Pinf1(:,:,t) = Pinf(:,:,t);
for i=1:pp
@@ -118,7 +121,10 @@ while newRank & t < smpl
Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t);
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
P0=Pinf(:,:,t+1);
@@ -138,7 +144,7 @@ Pinf1 = Pinf1(:,:,1:d);
notsteady = 1;
while notsteady & t<smpl
t = t+1;
- a1(:,t) = a(:,t);
+ a(:,t) = a1(:,t);
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
P1(:,:,t) = P(:,:,t);
for i=1:pp
@@ -152,7 +158,10 @@ while notsteady & t<smpl
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ end
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
end
@@ -169,16 +178,19 @@ if t<smpl
end
while t<smpl
t=t+1;
- a1(:,t) = a(:,t);
- for i=1:pp
+ a(:,t) = a1(:,t);
+ for i=1:pp
v(i,t) = Y(i,t) - a(mf(i),t) - a(mm+i,t) - trend(i,t);
- if Fi_s(i) > crit
+ if Fi_s(i) > crit
a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
end
end
- a(:,t+1) = T*a(:,t);
+ a1(:,t+1) = T*a(:,t);
+ for jnk=1:nk,
+ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
+ end
end
-a1(:,t+1) = a(:,t+1);
+
%% [2] Kalman smoother...
ri=r;
t = smpl+1;
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 4be555372ce55c0e6efe4ebb36adab70f4c041e9..025461b76ea2aac0afde13930f5288de520d060c 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -31,11 +31,11 @@ function dynare_estimation(var_list_)
global M_ options_ oo_ estim_params_ bayestopt_
var_list_ = check_list_of_variables(options_, M_, var_list_)
-if isempty(var_list_)
- return
-else
+%if isempty(var_list_)
+% return
+%else
options_.varlist = var_list_;
-end
+%end
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
for i = 1:size(M_.endo_names,1)
@@ -303,7 +303,7 @@ initial_estimation_checks(xparam1,gend,data);
if options_.mode_compute == 0 & length(options_.mode_file) == 0
if options_.smoother == 1
- [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
@@ -316,7 +316,8 @@ if options_.mode_compute == 0 & length(options_.mode_file) == 0
end
for i=1:M_.endo_nbr
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = atT(i,:)'';']);
- eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = filtered_state_vector(i,:)'';']);
+ eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = squeeze(aK(1,i,:))'';']);
+ eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = updated_variables(i,:)'';']);
end
for i=1:M_.exo_nbr
eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
@@ -956,7 +957,7 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
> 0) & options_.load_mh_file)) ...
| ~options_.smoother ) & M_.endo_nbr^2*gend < 1e7 % to be fixed
%% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
- [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
+ [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
@@ -975,7 +976,9 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
end
for i=1:M_.endo_nbr
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = atT(i,:)'';']);
- eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = filtered_state_vector(i,:)'';']);
+ eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = squeeze(aK(1,i,:))'';']);
+ eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ...
+ ' = updated variables(i,:)'';']);
end
[nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
if options_.TeX
diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index b72ecef2e7b9cd748c85bbe1ee6dfc1ed33add84..7ecb3d89c46d92f3b351c23aeb6f8e52ef47584b 100755
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ