diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 13f79ed6d03ea02dae3620a8c271fcbadc6a4be5..eeecc56608821f0f29310b27b7c3fa830d13ff0f 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -326,11 +326,11 @@ end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
-Z = BayesInfo.mf;
+Z = BayesInfo.mf; %selector for observed variables
no_missing_data_flag = ~DatasetInfo.missing.state;
-mm = length(T);
-pp = DynareDataset.vobs;
-rr = length(Q);
+mm = length(T); %number of states
+pp = DynareDataset.vobs; %number of observables
+rr = length(Q); %number of shocks
kalman_tol = DynareOptions.kalman_tol;
diffuse_kalman_tol = DynareOptions.diffuse_kalman_tol;
riccati_tol = DynareOptions.riccati_tol;
@@ -351,6 +351,7 @@ if (kalman_algo == 2) || (kalman_algo == 4)
H = diag(H);
mmm = mm;
else
+ %Augment state vector (follows Section 6.4.3 of DK (2012))
Z = [Z, eye(pp)];
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
diff --git a/matlab/kalman/likelihood/kalman_filter_d.m b/matlab/kalman/likelihood/kalman_filter_d.m
index 3b1161b7597cf30ae997391694f5020349995a47..4a2c94cf148e83e8a3ae028af7d4635ef82e61f8 100644
--- a/matlab/kalman/likelihood/kalman_filter_d.m
+++ b/matlab/kalman/likelihood/kalman_filter_d.m
@@ -8,8 +8,10 @@ function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, k
% a [double] mm*1 vector, levels of the state variables.
% Pinf [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% Pstar [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
-% kalman_tol [double] scalar, tolerance parameter (rcond).
-% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
+% kalman_tol [double] scalar, tolerance parameter (rcond) of F_star.
+% diffuse_kalman_tol [double] scalar, tolerance parameter (rcond) of Pinf to signify end of diffuse filtering and Finf.
+% riccati_tol [double] scalar, tolerance parameter (riccati iteration);
+% not used in this filter as usually diffuse phase will be left before convergence of filter to steady state.
% presample [integer] scalar, presampling if strictly positive.
% T [double] mm*mm matrix, transition matrix in the state equations.
% R [double] mm*rr matrix relating the structural innovations to the state vector.
@@ -28,10 +30,13 @@ function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, k
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
-% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
-% Analysis, vol. 24(1), pp. 85-98).
+% Models", S.J. Koopman and J. Durbin (2003), in Journal of Time Series
+% Analysis, vol. 24(1), pp. 85-98.
+% and
+% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
+% Second Edition, Ch. 5 and 7.2
-% Copyright (C) 2004-2013 Dynare Team
+% Copyright (C) 2004-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -62,48 +67,58 @@ s = 0;
while rank(Pinf,diffuse_kalman_tol) && (t<=last)
s = t-start+1;
- v = Y(:,t)-Z*a;
- Finf = Z*Pinf*Z';
- if rcond(Finf) < diffuse_kalman_tol
- if ~all(abs(Finf(:)) < diffuse_kalman_tol)
+ v = Y(:,t)-Z*a; %get prediction error v^(0) in (5.13) DK (2012)
+ Finf = Z*Pinf*Z'; % (5.7) in DK (2012)
+ %do case distinction based on whether F_{\infty,t} has full rank or 0 rank
+ if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
+ if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0
% The univariate diffuse kalman filter should be used instead.
return
- else
- Fstar = Z*Pstar*Z' + H;
- if rcond(Fstar) < kalman_tol
+ else %rank of F_{\infty,t} is 0
+ Fstar = Z*Pstar*Z' + H; % (5.7) in DK (2012)
+ if rcond(Fstar) < kalman_tol %F_{*} is singular
if ~all(abs(Fstar(:))<kalman_tol)
% The univariate diffuse kalman filter should be used.
return
- else
+ else %rank 0
a = T*a;
Pstar = T*Pstar*transpose(T)+QQ;
- Pinf = T*Pinf*transpose(T);
+ Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
end
else
iFstar = inv(Fstar);
dFstar = det(Fstar);
- Kstar = Pstar*Z'*iFstar;
- dlik(s)= log(dFstar) + v'*iFstar*v;
- Pinf = T*Pinf*transpose(T);
- Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ;
- a = T*(a+Kstar*v);
+ Kstar = Pstar*Z'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1}
+ dlik(s)= log(dFstar) + v'*iFstar*v; %set w_t to bottom case in bottom equation page 172, DK (2012)
+ Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
+ Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ; % (5.17) DK (2012)
+ a = T*(a+Kstar*v); % (5.13) DK (2012)
end
end
- else
- dlik(s)= log(det(Finf));
+ else %F_{\infty,t} positive definite
+ %To compare to DK (2012), this block makes use of the following transformation
+ %Kstar=K^{(1)}*T^{-1}=M_{*}*F^{(1)}+M_{\infty}*F^{(2)}
+ % =P_{*}*Z'*F^{(1)}+P_{\infty}*Z'*((-1)*(F_{\infty}^{-1})*F_{*}*(F_{\infty}^{-1}))
+ % =[P_{*}*Z'-Kinf*F_{*})]*F^{(1)}
+ %Make use of L^{0}'=(T-K^{(0)}*Z)'=(T-T*M_{\infty}*F^{(1)}*Z)'
+ % =(T-T*P_{\infty*Z'*F^{(1)}*Z)'=(T-T*Kinf*Z)'
+ % = (T*(I-*Kinf*Z))'=(I-Z'*Kinf')*T'
+ %P_{*}=T*P_{\infty}*L^{(1)}+T*P_{*}*L^{(0)}+RQR
+ % =T*[(P_{\infty}*(-K^{(1)*Z}))+P_{*}*(I-Z'*Kinf')*T'+RQR]
+ dlik(s)= log(det(Finf)); %set w_t to top case in bottom equation page 172, DK (2012)
iFinf = inv(Finf);
- Kinf = Pinf*Z'*iFinf;
- Fstar = Z*Pstar*Z' + H;
- Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf;
- Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ;
- Pinf = T*(Pinf-Pinf*Z'*Kinf')*T';
- a = T*(a+Kinf*v);
+ Kinf = Pinf*Z'*iFinf; %define Kinf=K_0*T^{-1} with M_{\infty}=Pinf*Z'
+ Fstar = Z*Pstar*Z' + H; %(5.7) DK(2012)
+ Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; %(5.12) DK(2012)
+ Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ; %(5.14) DK(2012)
+ Pinf = T*(Pinf-Pinf*Z'*Kinf')*T'; %(5.14) DK(2012)
+ a = T*(a+Kinf*v); %(5.13) DK(2012)
end
t = t+1;
end
if t>last
- warning(['There isn''t enough information to estimate the initial conditions of the nonstationary variables']);
+ warning(['kalman_filter_d: There isn''t enough information to estimate the initial conditions of the nonstationary variables. The diffuse Kalman filter never left the diffuse stage.']);
dLIK = NaN;
return
end
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
index 008e1aa8960c6a33a771b2597761b766d3071ac5..0b582a2d44b19664e57070e9ba723f89deac5ea6 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
@@ -35,10 +35,14 @@ function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,n
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
-% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
-% Analysis, vol. 24(1), pp. 85-98).
+% Models", S.J. Koopman and J. Durbin (2003), in Journal of Time Series
+% Analysis, vol. 24(1), pp. 85-98.
+% and
+% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
+% Second Edition, Ch. 5 and 7.2
-% Copyright (C) 2004-2012 Dynare Team
+%
+% Copyright (C) 2004-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -76,54 +80,57 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
s = t-start+1;
d_index = data_index{t};
if isempty(d_index)
+ %no observations, propagate forward without updating based on
+ %observations
a = T*a;
Pstar = T*Pstar*transpose(T)+QQ;
Pinf = T*Pinf*transpose(T);
else
- ZZ = Z(d_index,:);
- v = Y(d_index,t)-ZZ*a;
- Finf = ZZ*Pinf*ZZ';
- if rcond(Finf) < diffuse_kalman_tol
- if ~all(abs(Finf(:)) < diffuse_kalman_tol)
+ ZZ = Z(d_index,:); %span selector matrix
+ v = Y(d_index,t)-ZZ*a; %get prediction error v^(0) in (5.13) DK (2012)
+ Finf = ZZ*Pinf*ZZ'; % (5.7) in DK (2012)
+ if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
+ if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0
% The univariate diffuse kalman filter should be used.
return
- else
- Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index);
- if rcond(Fstar) < kalman_tol
+ else %rank of F_{\infty,t} is 0
+ Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index); % (5.7) in DK (2012)
+ if rcond(Fstar) < kalman_tol %F_{*} is singular
if ~all(abs(Fstar(:))<kalman_tol)
% The univariate diffuse kalman filter should be used.
return
- else
+ else %rank 0
a = T*a;
Pstar = T*Pstar*transpose(T)+QQ;
- Pinf = T*Pinf*transpose(T);
+ Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
end
else
iFstar = inv(Fstar);
dFstar = det(Fstar);
- Kstar = Pstar*ZZ'*iFstar;
- dlik(s) = log(dFstar) + v'*iFstar*v + length(d_index)*log(2*pi);
- Pinf = T*Pinf*transpose(T);
- Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ;
- a = T*(a+Kstar*v);
+ Kstar = Pstar*ZZ'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1}
+ dlik(s) = log(dFstar) + v'*iFstar*v + length(d_index)*log(2*pi); %set w_t to bottom case in bottom equation page 172, DK (2012)
+ Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
+ Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ; % (5.17) DK (2012)
+ a = T*(a+Kstar*v); % (5.13) DK (2012)
end
end
else
- dlik(s) = log(det(Finf));
+ dlik(s) = log(det(Finf)); %set w_t to top case in bottom equation page 172, DK (2012)
iFinf = inv(Finf);
Kinf = Pinf*ZZ'*iFinf;
- Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index);
- Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf;
- Pstar = T*(Pstar-Pstar*ZZ'*Kinf'-Pinf*ZZ'*Kstar')*T'+QQ;
- Pinf = T*(Pinf-Pinf*ZZ'*Kinf')*T';
- a = T*(a+Kinf*v);
+ %see notes in kalman_filter_d.m for details of computations
+ Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index); %(5.7) DK(2012)
+ Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf; %(5.12) DK(2012)
+ Pstar = T*(Pstar-Pstar*ZZ'*Kinf'-Pinf*ZZ'*Kstar')*T'+QQ; %(5.14) DK(2012)
+ Pinf = T*(Pinf-Pinf*ZZ'*Kinf')*T'; %(5.14) DK(2012)
+ a = T*(a+Kinf*v); %(5.13) DK(2012)
end
end
t = t+1;
end
if t==(last+1)
- warning(['There isn''t enough information to estimate the initial conditions of the nonstationary variables']);
+ warning(['kalman_filter_d: There isn''t enough information to estimate the initial conditions of the nonstationary variables. The diffuse Kalman filter never left the diffuse stage.']);
dLIK = NaN;
return
end
diff --git a/matlab/kalman/likelihood/univariate_kalman_filter.m b/matlab/kalman/likelihood/univariate_kalman_filter.m
index 2f79f42a5caf4392c38c9033771c7b5dd00c42c9..b0c8fec15cec23a14deeaf60bc6db9c388c16fb6 100644
--- a/matlab/kalman/likelihood/univariate_kalman_filter.m
+++ b/matlab/kalman/likelihood/univariate_kalman_filter.m
@@ -49,7 +49,7 @@ function [LIK, lik,a,P] = univariate_kalman_filter(data_index,number_of_observat
%! @item rr
%! Integer scalar, number of structural innovations.
%! @item Zflag
-%! Integer scalar, equal to 0 if Z is a vector of indices targeting the obseved variables in the state vector, equal to 1 if Z is a @var{pp}*@var{mm} matrix.
+%! Integer scalar, equal to 0 if Z is a vector of indices targeting the observed variables in the state vector, equal to 1 if Z is a @var{pp}*@var{mm} matrix.
%! @item diffuse_periods
%! Integer scalar, number of diffuse filter periods in the initialization step.
%! @end table
@@ -76,8 +76,15 @@ function [LIK, lik,a,P] = univariate_kalman_filter(data_index,number_of_observat
%! @ref{univariate_kalman_filter_ss}
%! @end deftypefn
%@eod:
+%
+% Algorithm:
+%
+% Uses the univariate filter as described in Durbin/Koopman (2012): "Time
+% Series Analysis by State Space Methods", Oxford University Press,
+% Second Edition, Ch. 6.4 + 7.2.5
-% Copyright (C) 2004-2012 Dynare Team
+
+% Copyright (C) 2004-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -154,7 +161,7 @@ else
LIK={inf,DLIK,Hess};
end
-while notsteady && t<=last
+while notsteady && t<=last %loop over t
s = t-start+1;
d_index = data_index{t};
if Zflag
@@ -163,22 +170,22 @@ while notsteady && t<=last
z = Z(d_index);
end
oldP = P(:);
- for i=1:rows(z)
+ for i=1:rows(z) %loop over i
if Zflag
- prediction_error = Y(d_index(i),t) - z(i,:)*a;
- PZ = P*z(i,:)';
- Fi = z(i,:)*PZ + H(d_index(i));
+ prediction_error = Y(d_index(i),t) - z(i,:)*a; % nu_{t,i} in 6.13 in DK (2012)
+ PZ = P*z(i,:)'; % Z_{t,i}*P_{t,i}*Z_{t,i}'
+ Fi = z(i,:)*PZ + H(d_index(i)); % F_{t,i} in 6.13 in DK (2012), relies on H being diagonal
else
- prediction_error = Y(d_index(i),t) - a(z(i));
- PZ = P(:,z(i));
- Fi = PZ(z(i)) + H(d_index(i));
+ prediction_error = Y(d_index(i),t) - a(z(i)); % nu_{t,i} in 6.13 in DK (2012)
+ PZ = P(:,z(i)); % Z_{t,i}*P_{t,i}*Z_{t,i}'
+ Fi = PZ(z(i)) + H(d_index(i)); % F_{t,i} in 6.13 in DK (2012), relies on H being diagonal
end
if Fi>kalman_tol
- Ki = PZ/Fi;
+ Ki = PZ/Fi; %K_{t,i} in 6.13 in DK (2012)
if t>=no_more_missing_observations
K(:,i) = Ki;
end
- lik(s,i) = log(Fi) + prediction_error*prediction_error/Fi + l2pi;
+ lik(s,i) = log(Fi) + (prediction_error*prediction_error)/Fi + l2pi; %Top equation p. 175 in DK (2012)
if analytic_derivation,
if analytic_derivation==2,
[Da,DP,DLIKt,D2a,D2P, Hesst] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady,D2a,D2Yss,D2P);
@@ -193,8 +200,11 @@ while notsteady && t<=last
end
dlik(s,:)=dlik(s,:)+DLIKt';
end
- a = a + Ki*prediction_error;
- P = P - PZ*Ki';
+ a = a + Ki*prediction_error; %filtering according to (6.13) in DK (2012)
+ P = P - PZ*Ki'; %filtering according to (6.13) in DK (2012)
+ else
+ % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
+ % p. 157, DK (2012)
end
end
if analytic_derivation,
@@ -204,8 +214,8 @@ while notsteady && t<=last
[Da,DP] = univariate_computeDstate(k,a,P,T,Da,DP,DT,DOm,notsteady);
end
end
- a = T*a;
- P = T*P*T' + QQ;
+ a = T*a; %transition according to (6.14) in DK (2012)
+ P = T*P*T' + QQ; %transition according to (6.14) in DK (2012)
if t>=no_more_missing_observations
notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_d.m b/matlab/kalman/likelihood/univariate_kalman_filter_d.m
index 4222f763f7d2f8ebba4780adc272629eb9587738..292a59c6921523faff222659dbafb8c29fa64771 100644
--- a/matlab/kalman/likelihood/univariate_kalman_filter_d.m
+++ b/matlab/kalman/likelihood/univariate_kalman_filter_d.m
@@ -80,6 +80,12 @@ function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index,
%! @ref{univariate_kalman_filter_ss}
%! @end deftypefn
%@eod:
+%
+% Algorithm:
+%
+% Uses the diffuse univariate filter as described in Durbin/Koopman (2012): "Time
+% Series Analysis by State Space Methods", Oxford University Press,
+% Second Edition, Ch. 5, 6.4 + 7.2.5
% Copyright (C) 2004-2016 Dynare Team
%
@@ -119,23 +125,26 @@ while newRank && (t<=last)
d_index = data_index{t};
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));
- Finf = Zi*Pinf*Zi';
+ prediction_error = Y(d_index(i),t) - Zi*a; % nu_{t,i} in 6.13 in DK (2012)
+ Fstar = Zi*Pstar*Zi' + H(d_index(i)); % F_{*,t} in 5.7 in DK (2012), relies on H being diagonal
+ Finf = Zi*Pinf*Zi'; % F_{\infty,t} in 5.7 in DK (2012), relies on H being diagonal
Kstar = Pstar*Zi';
- if Finf>diffuse_kalman_tol && newRank
+ if Finf>diffuse_kalman_tol && newRank % F_{\infty,t,i} = 0, use upper part of bracket on p. 175 DK (2012) for w_{t,i}
Kinf = Pinf*Zi';
Kinf_Finf = Kinf/Finf;
a = a + Kinf_Finf*prediction_error;
Pstar = Pstar + Kinf*(Kinf_Finf'*(Fstar/Finf)) - Kstar*Kinf_Finf' - Kinf_Finf*Kstar';
Pinf = Pinf - Kinf*Kinf_Finf';
- llik(s,d_index(i)) = log(Finf) + l2pi;
+ llik(s,d_index(i)) = log(Finf) + l2pi;
dlikk(s) = dlikk(s) + llik(s,d_index(i));
elseif Fstar>kalman_tol
- llik(s,d_index(i)) = log(Fstar) + prediction_error*prediction_error/Fstar + l2pi;
+ llik(s,d_index(i)) = log(Fstar) + (prediction_error*prediction_error/Fstar) + l2pi;
dlikk(s) = dlikk(s) + llik(s,d_index(i));
a = a+Kstar*(prediction_error/Fstar);
Pstar = Pstar-Kstar*(Kstar'/Fstar);
+ else
+ % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
+ % p. 157, DK (2012)
end
end
if newRank
diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
index 0f6171055c408c1cd9cfb8fd92a05ef103432d3e..cbc1d5775351fc9abf61bdda94c0d17905ef1c45 100644
--- a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
+++ b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
@@ -54,8 +54,10 @@ function [LIK,likk,a] = univariate_kalman_filter_ss(Y,start,last,a,P,kalman_tol,
%! @sp 1
%! @end deftypefn
%@eod:
+%
+% Algorithm: See univariate_kalman_filter.m
-% Copyright (C) 2011-2012 Dynare Team
+% Copyright (C) 2011-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -144,6 +146,9 @@ while t<=last
end
dlikk(s,:)=dlikk(s,:)+DLIKt';
end
+ else
+ % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
+ % p. 157, DK (2012)
end
end
if analytic_derivation,