Computes the diffuse likelihood of a state space model.
INPUTS
Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables.
start [integer] scalar, first observation.
last [integer] scalar, last observation.
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).
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.
Q [double] rr*rr covariance matrix of the structural innovations.
H [double] pp*pp covariance matrix of the measurement errors (if H is equal to zero (scalar) there is no measurement error).
Z [double] pp*mm matrix, selection matrix or pp linear independant combinations of the state vector.
mm [integer] scalar, number of state variables.
pp [integer] scalar, number of observed variables.
rr [integer] scalar, number of structural innovations.
OUTPUTS
LIK [double] scalar, minus loglikelihood
lik [double] smpl*1 vector, log density of each vector of observations.
a [double] mm*1 vector, current estimate of the state vector.
Pstar [double] mm*mm matrix, covariance matrix of the state vector.
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).0001 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) 0002 % Computes the diffuse likelihood of a state space model. 0003 % 0004 % INPUTS 0005 % Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables. 0006 % start [integer] scalar, first observation. 0007 % last [integer] scalar, last observation. 0008 % a [double] mm*1 vector, levels of the state variables. 0009 % Pinf [double] mm*mm matrix used to initialize the covariance matrix of the state vector. 0010 % Pstar [double] mm*mm matrix used to initialize the covariance matrix of the state vector. 0011 % kalman_tol [double] scalar, tolerance parameter (rcond). 0012 % riccati_tol [double] scalar, tolerance parameter (riccati iteration). 0013 % presample [integer] scalar, presampling if strictly positive. 0014 % T [double] mm*mm matrix, transition matrix in the state equations. 0015 % R [double] mm*rr matrix relating the structural innovations to the state vector. 0016 % Q [double] rr*rr covariance matrix of the structural innovations. 0017 % H [double] pp*pp covariance matrix of the measurement errors (if H is equal to zero (scalar) there is no measurement error). 0018 % Z [double] pp*mm matrix, selection matrix or pp linear independant combinations of the state vector. 0019 % mm [integer] scalar, number of state variables. 0020 % pp [integer] scalar, number of observed variables. 0021 % rr [integer] scalar, number of structural innovations. 0022 % 0023 % OUTPUTS 0024 % LIK [double] scalar, minus loglikelihood 0025 % lik [double] smpl*1 vector, log density of each vector of observations. 0026 % a [double] mm*1 vector, current estimate of the state vector. 0027 % Pstar [double] mm*mm matrix, covariance matrix of the state vector. 0028 % 0029 % REFERENCES 0030 % See "Filtering and Smoothing of State Vector for Diffuse State Space 0031 % Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series 0032 % Analysis, vol. 24(1), pp. 85-98). 0033 0034 % Copyright (C) 2004-2011 Dynare Team 0035 % 0036 % This file is part of Dynare. 0037 % 0038 % Dynare is free software: you can redistribute it and/or modify 0039 % it under the terms of the GNU General Public License as published by 0040 % the Free Software Foundation, either version 3 of the License, or 0041 % (at your option) any later version. 0042 % 0043 % Dynare is distributed in the hope that it will be useful, 0044 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0045 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0046 % GNU General Public License for more details. 0047 % 0048 % You should have received a copy of the GNU General Public License 0049 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0050 0051 % Get sample size. 0052 smpl = last-start+1; 0053 0054 % Initialize some variables. 0055 dF = 1; 0056 QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 0057 t = start; % Initialization of the time index. 0058 dlik = zeros(smpl,1); % Initialization of the vector gathering the densities. 0059 dLIK = Inf; % Default value of the log likelihood. 0060 oldK = Inf; 0061 s = 0; 0062 0063 while rank(Pinf,kalman_tol) && (t<=last) 0064 s = t-start+1; 0065 v = Y(:,t)-Z*a; 0066 Finf = Z*Pinf*Z'; 0067 if rcond(Finf) < kalman_tol 0068 if ~all(abs(Finf(:)) < kalman_tol) 0069 % The univariate diffuse kalman filter should be used instead. 0070 return 0071 else 0072 Fstar = Z*Pstar*Z' + H; 0073 if rcond(Fstar) < kalman_tol 0074 if ~all(abs(Fstar(:))<kalman_tol) 0075 % The univariate diffuse kalman filter should be used. 0076 return 0077 else 0078 a = T*a; 0079 Pstar = T*Pstar*transpose(T)+QQ; 0080 Pinf = T*Pinf*transpose(T); 0081 end 0082 else 0083 iFstar = inv(Fstar); 0084 dFstar = det(Fstar); 0085 Kstar = Pstar*Z'*iFstar; 0086 dlik(s)= log(dFstar) + v'*iFstar*v; 0087 Pinf = T*Pinf*transpose(T); 0088 Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ; 0089 a = T*(a+Kstar*v); 0090 end 0091 end 0092 else 0093 dlik(s)= log(det(Finf)); 0094 iFinf = inv(Finf); 0095 Kinf = Pinf*Z'*iFinf; 0096 Fstar = Z*Pstar*Z' + H; 0097 Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; 0098 Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ; 0099 Pinf = T*(Pinf-Pinf*Z'*Kinf')*T'; 0100 a = T*(a+Kinf*v); 0101 end 0102 t = t+1; 0103 end 0104 0105 if t>last 0106 warning(['There isn''t enough information to estimate the initial conditions of the nonstationary variables']); 0107 dLIK = NaN; 0108 return 0109 end 0110 0111 dlik = dlik(1:s); 0112 dlik = .5*(dlik + pp*log(2*pi)); 0113 0114 dLIK = sum(dlik(1+presample:end));