2011-09-19 17:01:24 +02:00
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)
2008-10-13 20:54:42 +02:00
% Computes the likelihood of a state space model in the case with missing observations.
%
% INPUTS
2008-10-16 20:03:33 +02:00
% data_index [cell] 1*smpl cell of column vectors of indices.
% number_of_observations [integer] scalar.
2011-09-19 17:01:24 +02:00
% 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.
2008-10-13 20:54:42 +02:00
%
2011-09-19 17:01:24 +02:00
% OUTPUTS
2010-12-10 22:31:30 +01:00
% LIK [double] scalar, MINUS loglikelihood
2008-10-13 20:54:42 +02:00
% lik [double] vector, density of observations in each period.
2011-09-19 17:01:24 +02:00
% a [double] mm*1 vector, estimated level of the states.
% P [double] mm*mm matrix, covariance matrix of the states.
2008-10-13 20:54:42 +02:00
%
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
2009-12-16 18:17:34 +01:00
2011-02-16 18:31:48 +01:00
% Copyright (C) 2004-2011 Dynare Team
2008-10-13 20:54:42 +02:00
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
2009-12-16 18:17:34 +01:00
2011-09-19 17:01:24 +02:00
% Set defaults
if nargin < 20
Zflag = 0 ;
diffuse_periods = 0 ;
end
if nargin < 21
diffuse_periods = 0 ;
end
if isempty ( Zflag )
Zflag = 0 ;
end
if isempty ( diffuse_periods )
diffuse_periods = 0 ;
end
% Get sample size.
smpl = last - start + 1 ;
% Initialize some variables.
dF = 1 ;
QQ = R * Q * transpose ( R ) ; % Variance of R times the vector of structural innovations.
t = start ; % Initialization of the time index.
lik = zeros ( smpl , 1 ) ; % Initialization of the vector gathering the densities.
LIK = Inf ; % Default value of the log likelihood.
2010-04-16 12:19:09 +02:00
oldK = Inf ;
2011-09-19 17:01:24 +02:00
notsteady = 1 ;
2009-12-16 18:17:34 +01:00
F_singular = 1 ;
2011-09-19 17:01:24 +02:00
while notsteady & t < = last
s = t - start + 1 ;
2009-12-16 18:17:34 +01:00
d_index = data_index { t } ;
if isempty ( d_index )
a = T * a ;
P = T * P * transpose ( T ) + QQ ;
else
2011-09-19 17:01:24 +02:00
% Compute the prediction error and its variance
if Zflag
z = Z ( d_index , : ) ;
v = Y ( d_index , t ) - z * a ;
F = z * P * z ' + H ( d_index , d_index ) ;
else
z = Z ( d_index ) ;
v = Y ( d_index , t ) - a ( z ) ;
F = P ( z , z ) + H ( d_index , d_index ) ;
2009-12-16 18:17:34 +01:00
end
if rcond ( F ) < kalman_tol
if ~ all ( abs ( F ( : ) ) < kalman_tol )
return
2008-10-13 20:54:42 +02:00
else
2009-12-16 18:17:34 +01:00
a = T * a ;
P = T * P * transpose ( T ) + QQ ;
2008-10-16 15:40:40 +02:00
end
2009-12-16 18:17:34 +01:00
else
F_singular = 0 ;
dF = det ( F ) ;
iF = inv ( F ) ;
2011-09-19 17:01:24 +02:00
lik ( s ) = log ( dF ) + transpose ( v ) * iF * v + length ( d_index ) * log ( 2 * pi ) ;
if Zflag
K = P * z ' * iF ;
P = T * ( P - K * z * P ) * transpose ( T ) + QQ ;
else
K = P ( : , z ) * iF ;
P = T * ( P - K * P ( z , : ) ) * transpose ( T ) + QQ ;
end
a = T * ( a + K * v ) ;
if t > = no_more_missing_observations
2010-06-24 18:18:57 +02:00
notsteady = max ( abs ( K ( : ) - oldK ) ) > riccati_tol ;
oldK = K ( : ) ;
end
2009-05-28 20:14:07 +02:00
end
2008-10-13 20:54:42 +02:00
end
2011-09-19 17:01:24 +02:00
t = t + 1 ;
2009-12-16 18:17:34 +01:00
end
if F_singular
error ( ' The variance of the forecast error remains singular until the end of the sample' )
end
2011-09-19 17:01:24 +02:00
% Divide by two.
lik ( 1 : s ) = . 5 * lik ( 1 : s ) ;
2009-12-16 18:17:34 +01:00
2011-09-19 17:01:24 +02:00
% Call steady state Kalman filter if needed.
if t < last
[ tmp , lik ( s + 1 : end ) ] = kalman_filter_ss ( Y , t , last , a , T , K , iF , dF , Z , pp , Zflag ) ;
end
2009-05-28 20:14:07 +02:00
2011-09-19 17:01:24 +02:00
% Compute minus the log-likelihood.
if presample
if presample > = diffuse_periods
lik = lik ( 1 + ( presample - diffuse_periods ) : end ) ;
end
end
LIK = sum ( lik ) ;