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