2 bugfixes for missing_DiffuseKalmanSmootherH1_Z.m that led to wrong results
1. In case of missing observations, a_t was not propagated forward to updated a_t+1 2. In the rank-deficient Finf case, Kstar was defined as T^(-1)*K^(0), while in the full rank it was defined as Kstar=K^(0), leading to wrong results when switches between the two clauses occurred. Moreover, the later backwards pass relied on Kstar=K^(0), leading to wrong results when the rank-deficient Finf case was triggered. The implementation now consistently follows the one in kalman_filter_d.mtime-shift
Johannes Pfeifer
parent
57d600301b
commit
ff1522a571
|
|
@ -43,6 +43,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal
|
|||
% 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).
|
||||
% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
|
||||
% Second Edition, Ch. 5
|
||||
|
||||
% Copyright (C) 2004-2016 Dynare Team
|
||||
%
|
||||
|
|
@ -83,16 +85,16 @@ K = zeros(mm,pp,smpl);
|
|||
L = zeros(mm,mm,smpl);
|
||||
Linf = zeros(mm,mm,smpl);
|
||||
Kstar = zeros(mm,pp,smpl);
|
||||
Kinf = 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;
|
||||
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(rr,smpl);
|
||||
epsilonhat = zeros(rr,smpl);
|
||||
r = zeros(mm,smpl+1);
|
||||
|
||||
t = 0;
|
||||
|
|
@ -100,54 +102,57 @@ while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
|
|||
t = t+1;
|
||||
di = data_index{t};
|
||||
if isempty(di)
|
||||
atilde(:,t) = a(:,t);
|
||||
%no observations, propagate forward without updating based on
|
||||
%observations
|
||||
atilde(:,t) = a(:,t);
|
||||
a(:,t+1) = T*atilde(:,t);
|
||||
Linf(:,:,t) = T;
|
||||
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T' + QQ;
|
||||
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
|
||||
else
|
||||
ZZ = Z(di,:);
|
||||
v(di,t)= Y(di,t) - ZZ*a(:,t);
|
||||
Finf = ZZ*Pinf(:,:,t)*ZZ';
|
||||
if rcond(Finf) < diffuse_kalman_tol
|
||||
if ~all(abs(Finf(:)) < diffuse_kalman_tol)
|
||||
ZZ = Z(di,:); %span selector matrix
|
||||
v(di,t)= Y(di,t) - ZZ*a(:,t); %get prediction error v^(0) in (5.13) DK (2012)
|
||||
Finf = ZZ*Pinf(:,:,t)*ZZ'; % (5.7) in DK (2012)
|
||||
if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
|
||||
if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0
|
||||
% The univariate diffuse kalman filter should be used.
|
||||
alphahat = Inf;
|
||||
return
|
||||
else
|
||||
Fstar(:,:,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);
|
||||
if rcond(Fstar(:,:,t)) < kalman_tol
|
||||
else %rank of F_{\infty,t} is 0
|
||||
Fstar(:,:,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); % (5.7) in DK (2012)
|
||||
if rcond(Fstar(:,:,t)) < kalman_tol %F_{*} is singular
|
||||
if ~all(abs(Fstar(:,:,t))<kalman_tol)
|
||||
% The univariate diffuse kalman filter should be used.
|
||||
alphahat = Inf;
|
||||
return
|
||||
else
|
||||
else %rank 0
|
||||
a(:,t+1) = T*a(:,t);
|
||||
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+QQ;
|
||||
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
|
||||
end
|
||||
else
|
||||
iFstar = inv(Fstar(:,:,t));
|
||||
Kstar(:,:,t) = Pstar(:,:,t)*ZZ'*iFstar(:,:,t);
|
||||
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
|
||||
Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kstar(:,:,t)')*T'+QQ;
|
||||
a(:,t+1) = T*(a(:,t)+Kstar(:,:,t)*v(:,t));
|
||||
Kstar(:,:,t) = Pstar(:,:,t)*ZZ'*iFstar(:,:,t); %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0)
|
||||
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); % DK (2012), 5.16
|
||||
Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kstar(:,:,t)')*T'+QQ; % (5.17) DK (2012) with L_0 plugged in
|
||||
a(:,t+1) = T*(a(:,t)+Kstar(:,:,t)*v(:,t)); % (5.13) DK (2012)
|
||||
end
|
||||
end
|
||||
else
|
||||
%see notes in kalman_filter_d.m for details of computations
|
||||
iFinf(di,di,t) = inv(Finf);
|
||||
PZI = Pinf(:,:,t)*ZZ'*iFinf(di,di,t);
|
||||
atilde(:,t) = a(:,t) + PZI*v(di,t);
|
||||
Kinf(:,di,t) = T*PZI;
|
||||
Linf(:,:,t) = T - Kinf(:,di,t)*ZZ;
|
||||
Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);
|
||||
Kstar(:,di,t) = (T*Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t);
|
||||
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*ZZ'*Kinf(:,di,t)'-T*Pinf(:,:,t)*ZZ'*Kstar(:,di,t)' + QQ;
|
||||
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*ZZ'*Kinf(:,di,t)';
|
||||
Kinf(:,di,t) = Pinf(:,:,t)*ZZ'*iFinf(di,di,t); %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z'
|
||||
atilde(:,t) = a(:,t) + Kinf(:,di,t)*v(di,t);
|
||||
Linf(:,:,t) = T - T*Kinf(:,di,t)*ZZ; %L^(0) in DK (2012), eq. 5.12
|
||||
Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); %(5.7) DK(2012)
|
||||
Kstar(:,di,t) = (Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t); %(5.12) DK(2012) with Kstar=T^{-1}*K^(1)
|
||||
Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kinf(:,di,t)'-Pinf(:,:,t)*ZZ'*Kstar(:,di,t)')*T' + QQ; %(5.14) DK(2012)
|
||||
Pinf(:,:,t+1) = T*(Pinf(:,:,t)-Pinf(:,:,t)*ZZ'*Kinf(:,di,t)')*T'; %(5.14) DK(2012)
|
||||
end
|
||||
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,
|
||||
for jnk=2:nk
|
||||
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
|
||||
end
|
||||
end
|
||||
|
|
@ -163,12 +168,12 @@ Pinf = Pinf(:,:,1:d);
|
|||
notsteady = 1;
|
||||
while notsteady && t<smpl
|
||||
t = t+1;
|
||||
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
|
||||
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); % make sure P is symmetric
|
||||
di = data_index{t};
|
||||
if isempty(di)
|
||||
atilde(:,t) = a(:,t);
|
||||
atilde(:,t) = a(:,t);
|
||||
L(:,:,t) = T;
|
||||
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
|
||||
P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %p. 111, DK(2012)
|
||||
else
|
||||
ZZ = Z(di,:);
|
||||
v(di,t) = Y(di,t) - ZZ*a(:,t);
|
||||
|
|
@ -220,34 +225,37 @@ end
|
|||
% $$$ PK(jnk,:,:,t+jnk) = Pf;
|
||||
% $$$ end
|
||||
% $$$ end
|
||||
%% backward pass
|
||||
t = smpl+1;
|
||||
while t>d+1
|
||||
t = t-1;
|
||||
di = data_index{t};
|
||||
if isempty(di)
|
||||
% in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
|
||||
r(:,t) = L(:,:,t)'*r(:,t+1);
|
||||
else
|
||||
ZZ = Z(di,:);
|
||||
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1);
|
||||
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %DK (2012), eq. 4.38
|
||||
end
|
||||
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
|
||||
etahat(:,t) = QRt*r(:,t);
|
||||
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); %DK (2012), eq. 4.35
|
||||
etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63
|
||||
end
|
||||
if d
|
||||
if d %diffuse periods
|
||||
% initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23
|
||||
r0 = zeros(mm,d+1);
|
||||
r0(:,d+1) = r(:,d+1);
|
||||
r1 = zeros(mm,d+1);
|
||||
for t = d:-1:1
|
||||
r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
|
||||
r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf
|
||||
di = data_index{t};
|
||||
if isempty(di)
|
||||
r1(:,t) = Linf(:,:,t)'*r1(:,t+1);
|
||||
else
|
||||
r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ...
|
||||
+ Linf(:,:,t)'*r1(:,t+1);
|
||||
r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ...
|
||||
+ Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1)
|
||||
end
|
||||
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
|
||||
etahat(:,t) = QRt*r0(:,t);
|
||||
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23
|
||||
etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -260,7 +268,6 @@ if decomp_flag
|
|||
eta_tm1t = QRt*Z(di,:)'*iF(di,di,t)*v(di,t);
|
||||
AAA = P(:,:,t)*Z(di,:)'*ZRQinv(di,di)*bsxfun(@times,Z(di,:)*R,eta_tm1t');
|
||||
% calculate decomposition
|
||||
Ttok = eye(mm,mm);
|
||||
decomp(1,:,:,t+1) = AAA;
|
||||
for h = 2:nk
|
||||
AAA = T*AAA;
|
||||
|
|
|
|||
Loading…
Reference in New Issue