Diffuse smoother: bug correction when singularity occurs
parent
4ad491b0a8
commit
32bf5e9a0e
|
@ -63,6 +63,7 @@ global options_
|
||||||
d = 0;
|
d = 0;
|
||||||
decomp = [];
|
decomp = [];
|
||||||
nk = options_.nk;
|
nk = options_.nk;
|
||||||
|
kalman_tol = options_.kalman_tol;
|
||||||
spinf = size(Pinf1);
|
spinf = size(Pinf1);
|
||||||
spstar = size(Pstar1);
|
spstar = size(Pstar1);
|
||||||
v = zeros(pp,smpl);
|
v = zeros(pp,smpl);
|
||||||
|
@ -91,14 +92,35 @@ etahat = zeros(rr,smpl);
|
||||||
r = zeros(mm,smpl+1);
|
r = zeros(mm,smpl+1);
|
||||||
|
|
||||||
t = 0;
|
t = 0;
|
||||||
while rank(Pinf(:,:,t+1),crit1) & t<smpl
|
while rank(Pinf(:,:,t+1),kalman_tol) && (t<smpl)
|
||||||
t = t+1;
|
t = t+1;
|
||||||
v(:,t)= Y(:,t) - Z*a(:,t);
|
v(:,t)= Y(:,t) - Z*a(:,t);
|
||||||
F = Z*Pinf(:,:,t)*Z';
|
Finf = Z*Pinf(:,:,t)*Z' ;
|
||||||
if rcond(F) < crit
|
if rcond(Finf) < kalman_tol
|
||||||
|
if ~all(abs(Finf(:)) < kalman_tol)
|
||||||
|
% The univariate diffuse kalman filter should be used.
|
||||||
return
|
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
|
end
|
||||||
iFinf(:,:,t) = inv(F);
|
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(Finf);
|
||||||
PZI = Pinf(:,:,t)*Z'*iFinf(:,:,t);
|
PZI = Pinf(:,:,t)*Z'*iFinf(:,:,t);
|
||||||
atilde(:,t) = a(:,t) + PZI*v(:,t);
|
atilde(:,t) = a(:,t) + PZI*v(:,t);
|
||||||
Kinf(:,:,t) = T*PZI;
|
Kinf(:,:,t) = T*PZI;
|
||||||
|
@ -113,6 +135,8 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
|
||||||
Kstar(:,:,t) = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
|
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;
|
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)';
|
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
|
||||||
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
d = t;
|
d = t;
|
||||||
P(:,:,d+1) = Pstar(:,:,d+1);
|
P(:,:,d+1) = Pstar(:,:,d+1);
|
||||||
|
|
|
@ -203,7 +203,7 @@ if any(any(H ~= 0)) % should be replaced by a flag
|
||||||
nobs,np,smpl,mf);
|
nobs,np,smpl,mf);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
elseif options_.kalman_algo == 2
|
elseif kalman_algo == 2
|
||||||
if ~estim_params_.ncn
|
if ~estim_params_.ncn
|
||||||
[alphahat,epsilonhat,etahat,ahat,aK] = ...
|
[alphahat,epsilonhat,etahat,ahat,aK] = ...
|
||||||
DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
|
DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
|
||||||
|
@ -289,7 +289,7 @@ else
|
||||||
data1,nobs,np,smpl);
|
data1,nobs,np,smpl);
|
||||||
end
|
end
|
||||||
if all(alphahat(:)==0)
|
if all(alphahat(:)==0)
|
||||||
options_.kalman_algo = 4;
|
kalman_algo = 4;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
if kalman_algo == 4
|
if kalman_algo == 4
|
||||||
|
|
Loading…
Reference in New Issue