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

correcting bugs in diffuse Kalman smoother

time-shift
Michel Juillard 2011-01-13 11:38:12 +01:00
parent 2e0a36ae9f
commit 6ea053e47d
3 changed files with 61 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
matlab/DiffuseKalmanSmoother1_Z.m
View File

@ -101,7 +101,7 @@ while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
% The univariate diffuse kalman filter should be used.
return
else
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
if rcond(Fstar(:,:,t)) < kalman_tol
if ~all(abs(Fstar(:,:,t))<kalman_tol)
% The univariate diffuse kalman filter should be used.
@ -133,7 +133,8 @@ while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
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)*Fstar(:,:,t)*Kstar(:,:,t)' + QQ;
% Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*Fstar(:,:,t)*Kstar(:,:,t)' + QQ;
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
end

56
matlab/DiffuseKalmanSmoother3.m
View File

@ -68,25 +68,18 @@ a = zeros(mm,smpl);
a1 = zeros(mm,smpl+1);
aK = zeros(nk,mm,smpl+nk);
PK = zeros(nk,mm,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);
Fi = zeros(pp,smpl_diff);
Fstar = zeros(pp,smpl);
Finf = zeros(pp,smpl);
Fi = zeros(pp,smpl);
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);
Linf = zeros(mm,mm,pp,smpl);
L0 = zeros(mm,mm,pp,smpl);
Kstar = zeros(mm,pp,smpl);
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;
Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
crit = options_.kalman_tol;
@ -132,29 +125,7 @@ while newRank & t < smpl
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
P0=Pinf(:,:,t);
% newRank = any(diag(P0(mf,mf))>crit);
% if newRank==0, id = i; end,
if ~isempty(options_.diffuse_d),
newRank = (icc<options_.diffuse_d);
%if newRank & any(diag(P0(mf,mf))>crit)==0;
if newRank & (any(diag(P0(mf,mf))>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(P0(mf,mf))>crit);
newRank = (any(diag(P0(mf,mf))>crit) | rank(P0,crit1));
if newRank==0,
options_.diffuse_d = icc;
end
end,
% if newRank==0,
% options_.diffuse_d=i; %this is buggy
% 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
@ -170,13 +141,20 @@ while newRank & t < smpl
for jnk=2:nk
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
if newRank
oldRank = rank(P0,crit);
else
oldRank = 0;
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);
if newRank
newRank = rank(P0,crit1);
end
if oldRank ~= newRank
disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!')
end
end

59
matlab/DiffuseKalmanSmootherH1_Z.m
View File

@ -97,25 +97,48 @@ 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
Finf = Z*Pinf(:,:,t)*Z';
if rcond(Finf) < kalman_tol
if ~all(abs(Finf(:)) < kalman_tol)
% The univariate diffuse kalman filter should be used.
return
else
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H;
if rcond(Fstar(:,:,t)) < kalman_tol
if ~all(abs(Fstar(:,:,t))<kalman_tol)
% The univariate diffuse kalman filter should be used.
return
else
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)*Z'*iFstar(:,:,t);
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*Z'*Kstar(:,:,t)')*T'+QQ;
a(:,:,t+1) = T*(a(:,:,t)+Kstar(:,:,t)*v(:,t));
end
end
else
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*dynare_squeeze(aK(jnk-1,:,t+jnk-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)'-T*Pinf(:,:,t)*Z'*Kstar(:,:,t)' + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
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*dynare_squeeze(aK(jnk-1,:,t+jnk-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);