Commit 1be52aaa authored by sebastien's avatar sebastien
Browse files

Beautification: removed tabulation characters which were left in previous beautification pass

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3300 ac1d8469-bf42-47a9-8791-bf33cf982152
parent a8614965
......@@ -54,31 +54,31 @@ function [alphahat,etahat,atilde, aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf1,Pstar
global options_
nk = options_.nk;
spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
aK = zeros(nk,mm,smpl+1);
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;
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);
r = zeros(mm,smpl+1);
steady = smpl;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
Z = zeros(pp,mm);
for i=1:pp;
......@@ -88,24 +88,24 @@ end
t = 0;
while rank(Pinf(:,:,t+1),crit1) & t<smpl
t = t+1;
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
if rcond(Pinf(mf,mf,t)) < crit
return
return
end
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
PZI = Pinf(:,mf,t)*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);
Kinf(:,:,t) = T*PZI;
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);
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t);
Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t);
Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
end
d = t;
P(:,:,d+1) = Pstar(:,:,d+1);
......@@ -121,7 +121,7 @@ while notsteady & t<smpl
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)) < crit
return
return
end
iF(:,:,t) = inv(P(mf,mf,t));
PZI = P(:,mf,t)*iF(:,:,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^(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);
......@@ -161,17 +161,17 @@ t = smpl+1;
while t>d+1
t = t-1;
r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
end
if d
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);
r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
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);
etahat(:,t) = QRt*r0(:,t);
end
end
......@@ -63,32 +63,32 @@ 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);
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;
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);
r = zeros(mm,smpl+1);
steady = smpl;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
t = 0;
while rank(Pinf(:,:,t+1),crit1) & t<smpl
......@@ -96,23 +96,23 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
v(:,t)= Y(:,t) - Z*a(:,t);
F = Z*Pinf(:,:,t)*Z';
if rcond(F) < crit
return
return
end
iFinf(:,:,t) = inv(F);
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);
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);
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
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)';
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
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);
......@@ -129,7 +129,7 @@ while notsteady & t<smpl
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
F = Z*P(:,:,t)*Z';
if rcond(F) < crit
return
return
end
iF(:,:,t) = inv(F);
PZI = P(:,:,t)*Z'*iF(:,:,t);
......@@ -139,9 +139,9 @@ while notsteady & t<smpl
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
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);
......@@ -163,17 +163,17 @@ while t<smpl
a(:,t+1) = T*atilde(:,t);
Pf = P(:,:,t);
for jnk=1:nk,
Pf = T*Pf*T' + QQ;
Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t);
PK(jnk,:,:,t+jnk) = Pf;
PK(jnk,:,:,t+jnk) = Pf;
end
end
t = smpl+1;
while t>d+1
t = t-1;
r(:,t) = Z'*iF(:,:,t)*v(:,t) + L(:,:,t)'*r(:,t+1);
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
end
if d
r0 = zeros(mm,d+1);
......@@ -182,8 +182,8 @@ if d
for t = d:-1:1
r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
r1(:,t) = Z'*(iFinf(:,:,t)*v(:,t)-Kstar(:,:,t)'*r0(:,t+1)) + Linf(:,:,t)'*r1(:,t+1);
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);
etahat(:,t) = QRt*r0(:,t);
end
end
......
......@@ -61,11 +61,11 @@ 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);
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);
PK = zeros(nk,mm,mm,smpl+nk);
......@@ -75,29 +75,29 @@ else
smpl_diff=rank(Pinf1);
end
Fstar = zeros(pp,smpl_diff);
Finf = zeros(pp,smpl_diff);
Fi = 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;
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;
Fstar = zeros(pp,smpl_diff);
Finf = zeros(pp,smpl_diff);
Fi = 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;
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);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
steady = smpl;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
Z = zeros(pp,mm);
for i=1:pp;
......@@ -106,7 +106,7 @@ end
t = 0;
icc=0;
newRank = rank(Pinf(:,:,1),crit1);
newRank = rank(Pinf(:,:,1),crit1);
while newRank & t < smpl
t = t+1;
a(:,t) = a1(:,t);
......@@ -115,21 +115,21 @@ while newRank & t < smpl
Pstar1(:,:,t) = Pstar(:,:,t);
Pinf1(:,:,t) = Pinf(:,:,t);
for i=1:pp
v(i,t) = Y(i,t)-a(mf(i),t)-trend(i,t);
Fstar(i,t) = Pstar(mf(i),mf(i),t);
Finf(i,t) = Pinf(mf(i),mf(i),t);
Kstar(:,i,t) = Pstar(:,mf(i),t);
v(i,t) = Y(i,t)-a(mf(i),t)-trend(i,t);
Fstar(i,t) = Pstar(mf(i),mf(i),t);
Finf(i,t) = Pinf(mf(i),mf(i),t);
Kstar(:,i,t) = Pstar(:,mf(i),t);
if Finf(i,t) > crit & newRank, % original MJ: if Finf(i,t) > crit
icc=icc+1;
Kinf(:,i,t) = Pinf(:,mf(i),t);
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))*Z(i,:)/Finf(i,t);
a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
Pstar(:,:,t) = Pstar(:,:,t) + ...
Kinf(:,i,t) = Pinf(:,mf(i),t);
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))*Z(i,:)/Finf(i,t);
a(:,t) = a(:,t) + Kinf(:,i,t)*v(i,t)/Finf(i,t);
Pstar(:,:,t) = Pstar(:,:,t) + ...
Kinf(:,i,t)*transpose(Kinf(:,i,t))*Fstar(i,t)/(Finf(i,t)*Finf(i,t)) - ...
(Kstar(:,i,t)*transpose(Kinf(:,i,t)) +...
Kinf(:,i,t)*transpose(Kstar(:,i,t)))/Finf(i,t);
Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t);
Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*transpose(Kinf(:,i,t))/Finf(i,t);
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
Pinf(:,:,t)=tril(Pinf(:,:,t))+transpose(tril(Pinf(:,:,t),-1));
% new terminiation criteria by M. Ratto
......@@ -158,23 +158,23 @@ while newRank & t < smpl
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. [stphane,11-03-2004].
%% rank(Pinf)=0. [stphane,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)*transpose(Kstar(:,i,t))/Fstar(i,t);
a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*transpose(Kstar(:,i,t))/Fstar(i,t);
Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1));
end
end
a1(:,t+1) = T*a(:,t);
a1(:,t+1) = T*a(:,t);
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
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);
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
......@@ -255,8 +255,8 @@ while t>d+1
end
end
r(:,t) = ri;
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
ri = T'*ri;
end
if d
......@@ -274,9 +274,9 @@ if d
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);
alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
r(:,t) = r0(:,t);
etahat(:,t) = QRt*r(:,t);
etahat(:,t) = QRt*r(:,t);
if t > 1
r0(:,t-1) = T'*r0(:,t);
r1(:,t-1) = T'*r1(:,t);
......
......@@ -54,32 +54,32 @@ function [alphahat,epsilonhat,etahat,atilde, aK] = DiffuseKalmanSmootherH1(T,R,Q
global options_
nk = options_.nk;
spinf = size(Pinf1);
spstar = size(Pstar1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
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);
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;
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);
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+1);
r = zeros(mm,smpl+1);
Z = zeros(pp,mm);
for i=1:pp;
......@@ -89,23 +89,23 @@ end
t = 0;
while rank(Pinf(:,:,t+1),crit1) & t<smpl
t = t+1;
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
if rcond(Pinf(mf,mf,t)) < crit
return
return
end
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
PZI = Pinf(:,mf,t)*iFinf(:,:,t);
atilde(:,t) = a(:,t) + PZI*v(:,t);
Kinf(:,:,t) = T*PZI;
a(:,t+1) = T*atilde(:,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(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t) + H;
Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Pstar(mf,mf,t) + H;
Kstar(:,:,t) = (T*Pstar(:,mf,t)-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)-T*Pstar(:,mf,t)*transpose(Kinf(:,:,t))-Kinf(:,:,t)*Pinf(mf,mf,t)*transpose(Kstar(:,:,t)) + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf,t)*transpose(Kinf(:,:,t));
end
d = t;
P(:,:,d+1) = Pstar(:,:,d+1);
......@@ -121,7 +121,7 @@ while notsteady & t<smpl
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
return
end
iF(:,:,t) = inv(P(mf,mf,t) + H);
PZI = P(:,mf,t)*iF(:,:,t);
......@@ -167,7 +167,7 @@ if d
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);