Harmonize filters/likelihood with smoothers by using new option diffuse_kalman_tol

time-shift
Marco Ratto 2015-04-08 15:49:12 +02:00
parent fe48de4406
commit ca8f0ea006
4 changed files with 16 additions and 18 deletions

View File

@ -328,6 +328,7 @@ mm = length(T);
pp = DynareDataset.vobs; pp = DynareDataset.vobs;
rr = length(Q); rr = length(Q);
kalman_tol = DynareOptions.kalman_tol; kalman_tol = DynareOptions.kalman_tol;
diffuse_kalman_tol = DynareOptions.kalman_tol;
riccati_tol = DynareOptions.riccati_tol; riccati_tol = DynareOptions.riccati_tol;
Y = transpose(DynareDataset.data)-trend; Y = transpose(DynareDataset.data)-trend;
@ -405,13 +406,13 @@ switch DynareOptions.lik_init
if no_missing_data_flag if no_missing_data_flag
[dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ... [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ... zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ... kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr); T,R,Q,H,Z,mm,pp,rr);
else else
[dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ... [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ... Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ... zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ... kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr); T,R,Q,H,Z,mm,pp,rr);
end end
diffuse_periods = length(dlik); diffuse_periods = length(dlik);
@ -448,7 +449,7 @@ switch DynareOptions.lik_init
DatasetInfo.missing.no_more_missing_observations, ... DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ... Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ... zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ... kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H1,Z,mmm,pp,rr); T,R,Q,H1,Z,mmm,pp,rr);
diffuse_periods = size(dlik,1); diffuse_periods = size(dlik,1);
end end

View File

@ -1,4 +1,4 @@
function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr) function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, kalman_tol, diffuse_kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr)
% Computes the diffuse likelihood of a state space model. % Computes the diffuse likelihood of a state space model.
% %
% INPUTS % INPUTS
@ -59,14 +59,13 @@ dlik = zeros(smpl,1); % Initialization of the vector gathering the densitie
dLIK = Inf; % Default value of the log likelihood. dLIK = Inf; % Default value of the log likelihood.
oldK = Inf; oldK = Inf;
s = 0; s = 0;
crit1=1.e-6;
while rank(Pinf,crit1) && (t<=last) while rank(Pinf,diffuse_kalman_tol) && (t<=last)
s = t-start+1; s = t-start+1;
v = Y(:,t)-Z*a; v = Y(:,t)-Z*a;
Finf = Z*Pinf*Z'; Finf = Z*Pinf*Z';
if rcond(Finf) < kalman_tol if rcond(Finf) < diffuse_kalman_tol
if ~all(abs(Finf(:)) < kalman_tol) if ~all(abs(Finf(:)) < diffuse_kalman_tol)
% The univariate diffuse kalman filter should be used instead. % The univariate diffuse kalman filter should be used instead.
return return
else else

View File

@ -1,7 +1,7 @@
function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,number_of_observations,no_more_missing_observations, ... function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,number_of_observations,no_more_missing_observations, ...
Y, start, last, ... Y, start, last, ...
a, Pinf, Pstar, ... a, Pinf, Pstar, ...
kalman_tol, riccati_tol, presample, ... kalman_tol, diffuse_kalman_tol, riccati_tol, presample, ...
T, R, Q, H, Z, mm, pp, rr) T, R, Q, H, Z, mm, pp, rr)
% Computes the diffuse likelihood of a state space model when some observations are missing. % Computes the diffuse likelihood of a state space model when some observations are missing.
% %
@ -66,14 +66,13 @@ t = start; % Initialization of the time index.
dlik = zeros(smpl,1); % Initialization of the vector gathering the densities. dlik = zeros(smpl,1); % Initialization of the vector gathering the densities.
dLIK = Inf; % Default value of the log likelihood. dLIK = Inf; % Default value of the log likelihood.
oldK = Inf; oldK = Inf;
crit1=1.e-6;
if isequal(H,0) if isequal(H,0)
H = zeros(pp,pp); H = zeros(pp,pp);
end end
s = 0; s = 0;
while rank(Pinf,crit1) && (t<=last) while rank(Pinf,diffuse_kalman_tol) && (t<=last)
s = t-start+1; s = t-start+1;
d_index = data_index{t}; d_index = data_index{t};
if isempty(d_index) if isempty(d_index)
@ -84,8 +83,8 @@ while rank(Pinf,crit1) && (t<=last)
ZZ = Z(d_index,:); ZZ = Z(d_index,:);
v = Y(d_index,t)-ZZ*a; v = Y(d_index,t)-ZZ*a;
Finf = ZZ*Pinf*ZZ'; Finf = ZZ*Pinf*ZZ';
if rcond(Finf) < kalman_tol if rcond(Finf) < diffuse_kalman_tol
if ~all(abs(Finf(:)) < kalman_tol) if ~all(abs(Finf(:)) < diffuse_kalman_tol)
% The univariate diffuse kalman filter should be used. % The univariate diffuse kalman filter should be used.
return return
else else

View File

@ -1,4 +1,4 @@
function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index, number_of_observations, no_more_missing_observations, Y, start, last, a, Pinf, Pstar, kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr) function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index, number_of_observations, no_more_missing_observations, Y, start, last, a, Pinf, Pstar, kalman_tol, diffuse_kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr)
% Computes the likelihood of a state space model (initialization with diffuse steps, univariate approach). % Computes the likelihood of a state space model (initialization with diffuse steps, univariate approach).
%@info: %@info:
@ -110,8 +110,7 @@ dLIK = Inf; % Default value of the log likelihood.
oldK = Inf; oldK = Inf;
llik = zeros(smpl,pp); llik = zeros(smpl,pp);
crit1 = 1.e-6; newRank = rank(Pinf,diffuse_kalman_tol);
newRank = rank(Pinf,crit1);
l2pi = log(2*pi); l2pi = log(2*pi);
while newRank && (t<=last) while newRank && (t<=last)
@ -139,7 +138,7 @@ while newRank && (t<=last)
end end
end end
if newRank if newRank
oldRank = rank(Pinf,crit1); oldRank = rank(Pinf,diffuse_kalman_tol);
else else
oldRank = 0; oldRank = 0;
end end
@ -147,7 +146,7 @@ while newRank && (t<=last)
Pstar = T*Pstar*T'+QQ; Pstar = T*Pstar*T'+QQ;
Pinf = T*Pinf*T'; Pinf = T*Pinf*T';
if newRank if newRank
newRank = rank(Pinf,crit1); newRank = rank(Pinf,diffuse_kalman_tol);
end end
if oldRank ~= newRank if oldRank ~= newRank
disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!') disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!')