From 9d91625c10374adaec0b75a137af917ef7ed0583 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Wed, 2 Nov 2011 11:10:58 +0100 Subject: [PATCH] fixing bug related to measurement errors --- matlab/DsgeLikelihood.m | 73 ++++++++++++------- .../likelihood/univariate_kalman_filter.m | 6 +- .../likelihood/univariate_kalman_filter_d.m | 4 +- .../likelihood/univariate_kalman_filter_ss.m | 5 +- 4 files changed, 53 insertions(+), 35 deletions(-) diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m index 406a265cf..4678313f4 100644 --- a/matlab/DsgeLikelihood.m +++ b/matlab/DsgeLikelihood.m @@ -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), ... diff --git a/matlab/kalman/likelihood/univariate_kalman_filter.m b/matlab/kalman/likelihood/univariate_kalman_filter.m index 3a088adce..e2f280a66 100644 --- a/matlab/kalman/likelihood/univariate_kalman_filter.m +++ b/matlab/kalman/likelihood/univariate_kalman_filter.m @@ -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 diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_d.m b/matlab/kalman/likelihood/univariate_kalman_filter_d.m index 7b9d28a13..a17c0cf1d 100644 --- a/matlab/kalman/likelihood/univariate_kalman_filter_d.m +++ b/matlab/kalman/likelihood/univariate_kalman_filter_d.m @@ -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 diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m index 51a893c88..785282863 100644 --- a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m +++ b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m @@ -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