Commit 9d91625c authored by MichelJuillard's avatar MichelJuillard

fixing bug related to measurement errors

parent 9ac3e7a5
......@@ -145,6 +145,7 @@ ys = [];
trend_coeff = [];
exit_flag = 1;
info = 0;
singularity_flag = 0;
% Set flag related to analytical derivatives.
if nargout > 9
......@@ -315,28 +316,25 @@ Y = DynareDataset.data-trend;
kalman_algo = DynareOptions.kalman_algo;
% resetting measurement errors covariance matrix for univariate filters
no_correlation_flag = 1;
if (kalman_algo == 2) || (kalman_algo == 4)
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
else
no_correlation_flag = 0;
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H = zeros(nobs,1);
mmm = mm+pp;
end
end
if no_correlation_flag
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
mmm = mm+pp;
end
end
......@@ -386,10 +384,33 @@ switch DynareOptions.lik_init
if isinf(dLIK)
% Go to univariate diffuse filter if singularity problem.
kalman_algo = 4;
singularity_flag = 1;
end
end
if (kalman_algo==4)
% Univariate Diffuse Kalman Filter
if singularity_flag
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H = zeros(nobs,1);
mmm = mm+pp;
end
end
% no need to test again for correlation elements
singularity_flag = 0;
end
[dLIK,tmp,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
......@@ -398,7 +419,7 @@ switch DynareOptions.lik_init
diffuse_periods = length(tmp);
end
case 4% Start from the solution of the Riccati equation.
if kalman_algo ~= 2
if kalman_algo ~= 2
kalman_algo = 1;
end
if isequal(H,0)
......@@ -476,8 +497,6 @@ end
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
singularity_flag = 0;
if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
if no_missing_data_flag
if DynareOptions.block == 1
......@@ -520,24 +539,22 @@ if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) )
if singularity_flag
if isequal(H,0)
H = zeros(nobs,1);
mmm = mm;
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
mmm = mm;
else
no_correlation_flag = 0;
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
H = zeros(nobs,1);
mmm = mm+pp;
end
end
if no_correlation_flag
mmm = mm;
else
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blckdiag(Pinf,zeros(pp));
mmm = mm+pp;
end
end
LIK = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
......
......@@ -39,7 +39,7 @@ function [LIK, likk,a,P] = univariate_kalman_filter(data_index,number_of_observa
%! @item R
%! Matrix (@var{mm}*@var{rr}) of doubles, second matrix of the state equation relating the structural innovations to the state variables.
%! @item H
%! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar).
%! Vector (@var{pp}) of doubles, diagonal of covariance matrix of the measurement errors (corelation among measurement errors is handled by a model transformation).
%! @item Z
%! Matrix (@var{pp}*@var{mm}) of doubles or vector of integers, matrix relating the states to the observed variables or vector of indices (depending on the value of @var{Zflag}).
%! @item mm
......@@ -132,10 +132,10 @@ while notsteady && t<=last
for i=1:rows(z)
if Zflag
prediction_error = Y(d_index(i),t) - z(i,:)*a;
Fi = z(i,:)*P*z(i,:)' + H(d_index(i),d_index(i));
Fi = z(i,:)*P*z(i,:)' + H(d_index(i));
else
prediction_error = Y(d_index(i),t) - a(z(i));
Fi = P(z(i),z(i)) + H(d_index(i),d_index(i));
Fi = P(z(i),z(i)) + H(d_index(i));
end
if Fi>kalman_tol
if Zflag
......
......@@ -41,7 +41,7 @@ function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index,
%! @item Q
%! Matrix (@var{rr}*@var{rr}) of doubles, covariance matrix of the structural innovations (noise in the state equation).
%! @item H
%! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar).
%! Vector (@var{pp}) of doubles, diagonal of covariance matrix of the measurement errors (corelation among measurement errors is handled by a model transformation).
%! @item Z
%! Matrix (@var{pp}*@var{mm}) of doubles, matrix relating the states to the observed variables.
%! @item mm
......@@ -119,7 +119,7 @@ while newRank && (t<=last)
for i=1:length(d_index)
Zi = Z(d_index(i),:);
prediction_error = Y(d_index(i),t) - Zi*a;
Fstar = Zi*Pstar*Zi' + H(d_index(i),d_index(i));
Fstar = Zi*Pstar*Zi' + H(d_index(i));
Finf = Zi*Pinf*Zi';
Kstar = Pstar*Zi';
if Finf>kalman_tol && newRank
......
......@@ -25,6 +25,7 @@ function [LIK,likk,a] = univariate_kalman_filter_ss(Y,start,last,a,P,kalman_tol,
%! @item T
%! Matrix (@var{mm}*@var{mm}) of doubles, transition matrix of the state equation.
%! @item H
%! Vector (@var{pp}) of doubles, diagonal of covariance matrix of the measurement errors (corelation among measurement errors is handled by a model transformation).
%! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar).
%! @item Z
%! Matrix (@var{pp}*@var{mm}) of doubles or vector of integers, matrix relating the states to the observed variables or vector of indices (depending on the value of @var{Zflag}).
......@@ -89,10 +90,10 @@ while t<=last
for i=1:pp
if Zflag
prediction_error = Y(i,t) - Z(i,:)*a;
Fi = Z(i,:)*PP*Z(i,:)' + H(i,i);
Fi = Z(i,:)*PP*Z(i,:)' + H(i);
else
prediction_error = Y(i,t) - a(Z(i));
Fi = PP(Z(i),Z(i)) + H(i,i);
Fi = PP(Z(i),Z(i)) + H(i);
end
if Fi>kalman_tol
if Zflag
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment