various bug fixes for filter and smoother with missing observations
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2282 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
parent
5bf9ed894c
commit
bf65cb4f10
|
@ -74,6 +74,7 @@ end
|
||||||
|
|
||||||
Fstar = zeros(pp,smpl_diff);
|
Fstar = zeros(pp,smpl_diff);
|
||||||
Finf = zeros(pp,smpl_diff);
|
Finf = zeros(pp,smpl_diff);
|
||||||
|
Fi = zeros(pp,smpl_diff);
|
||||||
Ki = zeros(mm,pp,smpl);
|
Ki = zeros(mm,pp,smpl);
|
||||||
Li = zeros(mm,mm,pp,smpl);
|
Li = zeros(mm,mm,pp,smpl);
|
||||||
Linf = zeros(mm,mm,pp,smpl_diff);
|
Linf = zeros(mm,mm,pp,smpl_diff);
|
||||||
|
|
|
@ -83,6 +83,7 @@ end
|
||||||
|
|
||||||
Fstar = zeros(pp,smpl_diff);
|
Fstar = zeros(pp,smpl_diff);
|
||||||
Finf = zeros(pp,smpl_diff);
|
Finf = zeros(pp,smpl_diff);
|
||||||
|
Fi = zeros(pp,smpl_diff);
|
||||||
Ki = zeros(mm,pp,smpl);
|
Ki = zeros(mm,pp,smpl);
|
||||||
Li = zeros(mm,mm,pp,smpl);
|
Li = zeros(mm,mm,pp,smpl);
|
||||||
Linf = zeros(mm,mm,pp,smpl_diff);
|
Linf = zeros(mm,mm,pp,smpl_diff);
|
||||||
|
|
|
@ -53,7 +53,7 @@ function [LIK, lik] = missing_observations_diffuse_kalman_filter(T,R,Q,H,Pinf,Ps
|
||||||
oldK = 0;
|
oldK = 0;
|
||||||
lik = zeros(smpl+1,1);
|
lik = zeros(smpl+1,1);
|
||||||
LIK = Inf;
|
LIK = Inf;
|
||||||
lik(smpl+1) = number_of_observations*pp*log(2*pi);
|
lik(smpl+1) = number_of_observations*log(2*pi);
|
||||||
notsteady = 1;
|
notsteady = 1;
|
||||||
reste = 0;
|
reste = 0;
|
||||||
|
|
||||||
|
|
|
@ -75,6 +75,7 @@ end
|
||||||
|
|
||||||
Fstar = zeros(pp,smpl_diff);
|
Fstar = zeros(pp,smpl_diff);
|
||||||
Finf = zeros(pp,smpl_diff);
|
Finf = zeros(pp,smpl_diff);
|
||||||
|
Fi = zeros(pp,smpl_diff);
|
||||||
Ki = zeros(mm,pp,smpl);
|
Ki = zeros(mm,pp,smpl);
|
||||||
Li = zeros(mm,mm,pp,smpl);
|
Li = zeros(mm,mm,pp,smpl);
|
||||||
Linf = zeros(mm,mm,pp,smpl_diff);
|
Linf = zeros(mm,mm,pp,smpl_diff);
|
||||||
|
@ -211,34 +212,34 @@ while notsteady & t<smpl
|
||||||
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
|
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
|
||||||
end
|
end
|
||||||
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
|
P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ;
|
||||||
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
|
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
|
||||||
end
|
|
||||||
P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
|
|
||||||
Fi_s = Fi(:,t);
|
|
||||||
Ki_s = Ki(:,:,t);
|
|
||||||
L_s =Li(:,:,:,t);
|
|
||||||
if t<smpl
|
|
||||||
t_steady = t+1;
|
|
||||||
P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
|
|
||||||
Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
|
|
||||||
Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
|
|
||||||
Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
|
|
||||||
end
|
|
||||||
while t<smpl
|
|
||||||
t=t+1;
|
|
||||||
a(:,t) = a1(:,t);
|
|
||||||
di = data_index{t}';
|
|
||||||
for i=di
|
|
||||||
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
|
|
||||||
a1(:,t+1) = T*a(:,t);
|
|
||||||
for jnk=1:nk,
|
|
||||||
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
% $$$ P_s=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
|
||||||
|
% $$$ Fi_s = Fi(:,t);
|
||||||
|
% $$$ Ki_s = Ki(:,:,t);
|
||||||
|
% $$$ L_s =Li(:,:,:,t);
|
||||||
|
% $$$ if t<smpl
|
||||||
|
% $$$ t_steady = t+1;
|
||||||
|
% $$$ P = cat(3,P(:,:,1:t),repmat(P(:,:,t),[1 1 smpl-t_steady+1]));
|
||||||
|
% $$$ Fi = cat(2,Fi(:,1:t),repmat(Fi_s,[1 1 smpl-t_steady+1]));
|
||||||
|
% $$$ Li = cat(4,Li(:,:,:,1:t),repmat(L_s,[1 1 smpl-t_steady+1]));
|
||||||
|
% $$$ Ki = cat(3,Ki(:,:,1:t),repmat(Ki_s,[1 1 smpl-t_steady+1]));
|
||||||
|
% $$$ end
|
||||||
|
% $$$ while t<smpl
|
||||||
|
% $$$ t=t+1;
|
||||||
|
% $$$ a(:,t) = a1(:,t);
|
||||||
|
% $$$ di = data_index{t}';
|
||||||
|
% $$$ for i=di
|
||||||
|
% $$$ 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
|
||||||
|
% $$$ a1(:,t+1) = T*a(:,t);
|
||||||
|
% $$$ for jnk=1:nk,
|
||||||
|
% $$$ aK(jnk,:,t+jnk) = T^jnk*a(:,t);
|
||||||
|
% $$$ end
|
||||||
|
% $$$ end
|
||||||
ri=zeros(mm,1);
|
ri=zeros(mm,1);
|
||||||
t = smpl+1;
|
t = smpl+1;
|
||||||
while t>d+1
|
while t>d+1
|
||||||
|
|
|
@ -84,6 +84,7 @@ end
|
||||||
|
|
||||||
Fstar = zeros(pp,smpl_diff);
|
Fstar = zeros(pp,smpl_diff);
|
||||||
Finf = zeros(pp,smpl_diff);
|
Finf = zeros(pp,smpl_diff);
|
||||||
|
Fi = zeros(pp,smpl);
|
||||||
Ki = zeros(mm,pp,smpl);
|
Ki = zeros(mm,pp,smpl);
|
||||||
Li = zeros(mm,mm,pp,smpl);
|
Li = zeros(mm,mm,pp,smpl);
|
||||||
Linf = zeros(mm,mm,pp,smpl_diff);
|
Linf = zeros(mm,mm,pp,smpl_diff);
|
||||||
|
@ -214,39 +215,39 @@ while notsteady & t<smpl
|
||||||
PK(jnk,:,:,t+jnk) = Pf;
|
PK(jnk,:,:,t+jnk) = Pf;
|
||||||
end
|
end
|
||||||
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
|
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
|
||||||
notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
|
% 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);
|
|
||||||
di = data_index{t}';
|
|
||||||
for i=di
|
|
||||||
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
|
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);
|
||||||
|
% $$$ di = data_index{t}';
|
||||||
|
% $$$ for i=di
|
||||||
|
% $$$ 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=zeros(mm,1);
|
ri=zeros(mm,1);
|
||||||
t = smpl+1;
|
t = smpl+1;
|
||||||
while t > d+1
|
while t > d+1
|
||||||
|
|
Loading…
Reference in New Issue