stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bug fix and efficiency changes in smoothing routines.

time-shift
Stéphane Adjemian (Charybdis) 2010-02-12 12:31:32 +01:00
parent a828c3d2c1
commit a83c0e68c7
13 changed files with 131 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
matlab/DiffuseKalmanSmoother1.m
View File

@ -59,7 +59,7 @@ spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+1);
aK = zeros(nk,mm,smpl+nk);
iF = zeros(pp,pp,smpl);
Fstar = zeros(pp,pp,smpl);
iFinf = zeros(pp,pp,smpl);
@ -71,11 +71,11 @@ 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;
crit1 = 1.e-8;
steady = smpl;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
@ -99,7 +99,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
a(:,t+1) = T*atilde(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t);
@ -129,9 +129,9 @@ while notsteady & t<smpl
K(:,:,t) = T*PZI;
L(:,:,t) = T-K(:,:,t)*Z;
a(:,t+1) = T*atilde(:,t);
aK(1,:,t+1) = a(:,t+1);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-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);
@ -154,7 +154,7 @@ while t<smpl
a(:,t+1) = T*a(:,t) + K_s*v(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
t = smpl+1;
@ -174,4 +174,4 @@ if d
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
etahat(:,t) = QRt*r0(:,t);
end
end
end

16
matlab/DiffuseKalmanSmoother1_Z.m
View File

@ -105,8 +105,8 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
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);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
@ -138,10 +138,13 @@ while notsteady & t<smpl
L(:,:,t) = T-K(:,:,t)*Z;
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
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);
@ -162,10 +165,13 @@ while t<smpl
atilde(:,t) = a(:,t) + PZI*v(:,t);
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
end
t = smpl+1;

23
matlab/DiffuseKalmanSmoother3.m
View File

@ -165,13 +165,14 @@ while newRank & t < smpl
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
P0=Pinf(:,:,t+1);
P0 = Pinf(:,:,t+1);
if newRank,
%newRank = any(diag(P0(mf,mf))>crit);
newRank = rank(P0,crit1);
@ -207,12 +208,15 @@ while notsteady & t<smpl
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
end
end
a1(:,t+1) = T*a(:,t);
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
@ -239,10 +243,13 @@ while t<smpl
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
end
ri=zeros(mm,1);

15
matlab/DiffuseKalmanSmoother3_Z.m
View File

@ -91,7 +91,6 @@ Ki = zeros(mm,pp,smpl);
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;
@ -162,9 +161,10 @@ while newRank & t < smpl
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);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squueze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
@ -205,11 +205,14 @@ while notsteady & t<smpl
end
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1)
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);

23
matlab/DiffuseKalmanSmootherH1.m
View File

@ -98,8 +98,9 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
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);
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t) + H;
@ -118,25 +119,26 @@ Pinf = Pinf(:,:,1:d);
notsteady = 1;
while notsteady & t<smpl
t = t+1;
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
if rcond(P(mf,mf,t) + H) < crit
return
end
return
end
iF(:,:,t) = inv(P(mf,mf,t) + H);
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*atilde(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-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
if t<smpl
PZI_s = PZI;
PZI_s = PZI;
K_s = K(:,:,t);
iF_s = iF(:,:,t);
P_s = P(:,:,t+1);
@ -151,8 +153,9 @@ while t<smpl
v(:,t) = Y(:,t) - a(mf,t) - trend(:,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);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
t = smpl+1;

18
matlab/DiffuseKalmanSmootherH1_Z.m
View File

@ -108,8 +108,8 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
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);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
@ -141,10 +141,13 @@ while notsteady & t<smpl
L(:,:,t) = T-K(:,:,t)*Z;
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
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);
@ -165,10 +168,13 @@ while t<smpl
atilde(:,t) = a(:,t) + PZI*v(:,t);
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
end
t = smpl+1;

17
matlab/DiffuseKalmanSmootherH3.m
View File

@ -163,9 +163,10 @@ while newRank & t < smpl
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1)
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
@ -206,8 +207,9 @@ while notsteady & t<smpl
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
@ -233,8 +235,9 @@ while t<smpl
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
ri=zeros(mm,1);

20
matlab/DiffuseKalmanSmootherH3_Z.m
View File

@ -91,7 +91,6 @@ 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;
@ -163,9 +162,10 @@ while newRank & t < smpl
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);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
@ -207,10 +207,13 @@ while notsteady & t<smpl
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)) ;
end
end
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
@ -238,11 +241,14 @@ while t<smpl
end
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
end
ri=zeros(mm,1);

15
matlab/DiffuseKalmanSmootherH3corr.m
View File

@ -122,8 +122,9 @@ while newRank & t < smpl
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
aK(1,:,t+1) = a1(:,t+1)
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
@ -159,8 +160,9 @@ while notsteady & t<smpl
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
@ -186,8 +188,9 @@ while t<smpl
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end

10
matlab/missing_DiffuseKalmanSmoother1.m
View File

@ -60,7 +60,7 @@ spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+1);
aK = zeros(nk,mm,smpl+nk);
iF = zeros(pp,pp,smpl);
Fstar = zeros(pp,pp,smpl);
iFinf = zeros(pp,pp,smpl);
@ -113,9 +113,9 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf1,t)* ...
transpose(Kinf(:,di,t));
end
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
d = t;
@ -150,7 +150,7 @@ while notsteady & t<smpl
end
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-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);

9
matlab/missing_DiffuseKalmanSmoother1_Z.m
View File

@ -121,7 +121,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
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);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
d = t;
@ -157,10 +157,13 @@ while notsteady & t<smpl
end
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
end

18
matlab/missing_DiffuseKalmanSmoother3.m
View File

@ -91,11 +91,11 @@ Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
crit = options_.kalman_tol;
crit1 = 1.e-6;
crit1 = 1.e-6;
steady = smpl;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
@ -167,16 +167,17 @@ while newRank & t < smpl
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
end
end
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
P0=Pinf(:,:,t+1);
if newRank,
%newRank = any(diag(P0(mf,mf))>crit);
newRank = rank(P0,crit1);
newRank = rank(P0,crit1);
end
end
@ -212,10 +213,13 @@ while notsteady & t<smpl
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a1(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);

15
matlab/missing_DiffuseKalmanSmoother3_Z.m
View File

@ -92,7 +92,6 @@ Ki = zeros(mm,pp,smpl);
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;
@ -164,9 +163,10 @@ while newRank & t < smpl
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);
a1(:,t+1) = T*a(:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
@ -209,10 +209,13 @@ while notsteady & t<smpl
end
a1(:,t+1) = T*a(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
PK(jnk,:,:,t+jnk) = Pf;
if jnk>1
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
end
end
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);