Computes the likelihood of a state space model in the case with missing observations.
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 data.
start [integer] scalar, index of the first observation.
last [integer] scalar, index of the last observation.
a [double] pp*1 vector, initial level of the state vector.
P [double] pp*pp matrix, covariance matrix of the initial 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 transition matrix of the state equation.
Q [double] rr*rr covariance matrix of the structural innovations.
R [double] mm*rr matrix, mapping structural innovations to state variables.
H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
Z [integer] pp*1 vector of indices for the observed variables.
mm [integer] scalar, dimension of the state vector.
pp [integer] scalar, number of observed variables.
rr [integer] scalar, number of structural innovations.
OUTPUTS
LIK [double] scalar, MINUS loglikelihood
lik [double] vector, density of observations in each period.
a [double] mm*1 vector, estimated level of the states.
P [double] mm*mm matrix, covariance matrix of the states.
NOTES
The vector "lik" is used to evaluate the jacobian of the likelihood.0001 function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,number_of_observations,no_more_missing_observations,Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods) 0002 % Computes the likelihood of a state space model in the case with missing observations. 0003 % 0004 % INPUTS 0005 % data_index [cell] 1*smpl cell of column vectors of indices. 0006 % number_of_observations [integer] scalar. 0007 % no_more_missing_observations [integer] scalar. 0008 % Y [double] pp*smpl matrix of data. 0009 % start [integer] scalar, index of the first observation. 0010 % last [integer] scalar, index of the last observation. 0011 % a [double] pp*1 vector, initial level of the state vector. 0012 % P [double] pp*pp matrix, covariance matrix of the initial state vector. 0013 % kalman_tol [double] scalar, tolerance parameter (rcond). 0014 % riccati_tol [double] scalar, tolerance parameter (riccati iteration). 0015 % presample [integer] scalar, presampling if strictly positive. 0016 % T [double] mm*mm transition matrix of the state equation. 0017 % Q [double] rr*rr covariance matrix of the structural innovations. 0018 % R [double] mm*rr matrix, mapping structural innovations to state variables. 0019 % H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors. 0020 % Z [integer] pp*1 vector of indices for the observed variables. 0021 % mm [integer] scalar, dimension of the state vector. 0022 % pp [integer] scalar, number of observed variables. 0023 % rr [integer] scalar, number of structural innovations. 0024 % 0025 % OUTPUTS 0026 % LIK [double] scalar, MINUS loglikelihood 0027 % lik [double] vector, density of observations in each period. 0028 % a [double] mm*1 vector, estimated level of the states. 0029 % P [double] mm*mm matrix, covariance matrix of the states. 0030 % 0031 % 0032 % NOTES 0033 % The vector "lik" is used to evaluate the jacobian of the likelihood. 0034 0035 % Copyright (C) 2004-2011 Dynare Team 0036 % 0037 % This file is part of Dynare. 0038 % 0039 % Dynare is free software: you can redistribute it and/or modify 0040 % it under the terms of the GNU General Public License as published by 0041 % the Free Software Foundation, either version 3 of the License, or 0042 % (at your option) any later version. 0043 % 0044 % Dynare is distributed in the hope that it will be useful, 0045 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0046 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0047 % GNU General Public License for more details. 0048 % 0049 % You should have received a copy of the GNU General Public License 0050 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0051 0052 % Set defaults 0053 if nargin<20 0054 Zflag = 0; 0055 diffuse_periods = 0; 0056 end 0057 0058 if nargin<21 0059 diffuse_periods = 0; 0060 end 0061 0062 if isempty(Zflag) 0063 Zflag = 0; 0064 end 0065 0066 if isempty(diffuse_periods) 0067 diffuse_periods = 0; 0068 end 0069 0070 % Get sample size. 0071 smpl = last-start+1; 0072 0073 % Initialize some variables. 0074 dF = 1; 0075 QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 0076 t = start; % Initialization of the time index. 0077 lik = zeros(smpl,1); % Initialization of the vector gathering the densities. 0078 LIK = Inf; % Default value of the log likelihood. 0079 oldK = Inf; 0080 notsteady = 1; 0081 F_singular = 1; 0082 0083 while notsteady & t<=last 0084 s = t-start+1; 0085 d_index = data_index{t}; 0086 if isempty(d_index) 0087 a = T*a; 0088 P = T*P*transpose(T)+QQ; 0089 else 0090 % Compute the prediction error and its variance 0091 if Zflag 0092 z = Z(d_index,:); 0093 v = Y(d_index,t)-z*a; 0094 F = z*P*z' + H(d_index,d_index); 0095 else 0096 z = Z(d_index); 0097 v = Y(d_index,t) - a(z); 0098 F = P(z,z) + H(d_index,d_index); 0099 end 0100 if rcond(F) < kalman_tol 0101 if ~all(abs(F(:))<kalman_tol) 0102 return 0103 else 0104 a = T*a; 0105 P = T*P*transpose(T)+QQ; 0106 end 0107 else 0108 F_singular = 0; 0109 dF = det(F); 0110 iF = inv(F); 0111 lik(s) = log(dF) + transpose(v)*iF*v + length(d_index)*log(2*pi); 0112 if Zflag 0113 K = P*z'*iF; 0114 P = T*(P-K*z*P)*transpose(T)+QQ; 0115 else 0116 K = P(:,z)*iF; 0117 P = T*(P-K*P(z,:))*transpose(T)+QQ; 0118 end 0119 a = T*(a+K*v); 0120 if t>=no_more_missing_observations 0121 notsteady = max(abs(K(:)-oldK))>riccati_tol; 0122 oldK = K(:); 0123 end 0124 end 0125 end 0126 t = t+1; 0127 end 0128 0129 if F_singular 0130 error('The variance of the forecast error remains singular until the end of the sample') 0131 end 0132 0133 % Divide by two. 0134 lik(1:s) = .5*lik(1:s); 0135 0136 % Call steady state Kalman filter if needed. 0137 if t<=last 0138 [tmp, lik(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag); 0139 end 0140 0141 % Compute minus the log-likelihood. 0142 if presample>=diffuse_periods 0143 LIK = sum(lik(1+presample-diffuse_periods:end)); 0144 else 0145 LIK = sum(lik); 0146 end