2017-04-27 10:44:27 +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,rescale_prediction_error_covariance,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.
2017-05-16 15:10:20 +02:00
% Y [double] pp*smpl matrix of data.
2011-09-19 17:01:24 +02:00
% start [integer] scalar, index of the first observation.
% last [integer] scalar, index of the last observation.
2021-01-18 17:33:21 +01:00
% a [double] pp*1 vector, levels of the predicted initial state variables (E_{0}(alpha_1)).
2011-09-19 17:01:24 +02:00
% P [double] pp*pp matrix, covariance matrix of the initial state vector.
% kalman_tol [double] scalar, tolerance parameter (rcond).
2017-05-16 15:10:20 +02:00
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
2011-09-19 17:01:24 +02:00
% presample [integer] scalar, presampling if strictly positive.
% T [double] mm*mm transition matrix of the state equation.
2017-05-16 15:10:20 +02:00
% Q [double] rr*rr covariance matrix of the structural innovations.
2011-09-19 17:01:24 +02:00
% R [double] mm*rr matrix, mapping structural innovations to state variables.
2017-05-16 15:10:20 +02:00
% 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.
2011-09-19 17:01:24 +02:00
% mm [integer] scalar, dimension of the state vector.
% pp [integer] scalar, number of observed variables.
2017-05-16 15:10:20 +02:00
% rr [integer] scalar, number of structural innovations.
%
% 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.
2021-01-18 17:33:21 +01:00
% a [double] mm*1 vector, current estimate of the state vector tomorrow (E_{T}(alpha_{T+1})).
2017-05-16 15:10:20 +02:00
% P [double] mm*mm matrix, covariance matrix of the states.
%
2008-10-13 20:54:42 +02:00
%
% NOTES
2017-05-16 15:10:20 +02:00
% The vector "lik" is used to evaluate the jacobian of the likelihood.
2009-12-16 18:17:34 +01:00
2021-05-02 17:07:12 +02:00
% Copyright (C) 2004-2021 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
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://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
2012-06-19 16:35:44 +02:00
if isequal ( H , 0 )
H = zeros ( pp , pp ) ;
end
2011-09-19 17:01:24 +02:00
% Get sample size.
smpl = last - start + 1 ;
% Initialize some variables.
dF = 1 ;
2021-05-02 17:07:12 +02:00
isqvec = false ;
if ndim ( Q ) > 2
Qvec = Q ;
Q = Q ( : , : , 1 ) ;
isqvec = true ;
end
2011-09-19 17:01:24 +02:00
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 ;
2017-04-27 10:44:27 +02:00
F_singular = true ;
2012-07-01 15:19:36 +02:00
s = 0 ;
2020-01-03 13:19:57 +01:00
rescale_prediction_error_covariance0 = rescale_prediction_error_covariance ;
2012-07-01 15:19:36 +02:00
2015-11-08 18:59:30 +01:00
while notsteady && t < = last
2011-09-19 17:01:24 +02:00
s = t - start + 1 ;
2009-12-16 18:17:34 +01:00
d_index = data_index { t } ;
2021-05-02 17:07:12 +02:00
if isqvec
QQ = R * Qvec ( : , : , t + 1 ) * transpose ( R ) ;
end
2009-12-16 18:17:34 +01:00
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
2017-04-27 10:44:27 +02:00
badly_conditioned_F = false ;
if rescale_prediction_error_covariance
sig = sqrt ( diag ( F ) ) ;
if any ( diag ( F ) < kalman_tol ) || rcond ( F ./ ( sig * sig ' ) ) < kalman_tol
badly_conditioned_F = true ;
end
else
if rcond ( F ) < kalman_tol
2020-01-03 13:19:57 +01:00
sig = sqrt ( diag ( F ) ) ;
if any ( diag ( F ) < kalman_tol ) || rcond ( F ./ ( sig * sig ' ) ) < kalman_tol
badly_conditioned_F = true ;
else
rescale_prediction_error_covariance = 1 ;
end
% badly_conditioned_F = true;
2017-04-27 10:44:27 +02:00
end
end
if badly_conditioned_F
2009-12-16 18:17:34 +01:00
if ~ all ( abs ( F ( : ) ) < kalman_tol )
2017-04-27 10:44:27 +02:00
% Use univariate filter.
2009-12-16 18:17:34 +01:00
return
2008-10-13 20:54:42 +02:00
else
2017-04-27 10:44:27 +02:00
% Pathological case, discard draw
2016-06-19 17:13:04 +02:00
return
2008-10-16 15:40:40 +02:00
end
2009-12-16 18:17:34 +01:00
else
2017-04-27 10:44:27 +02:00
F_singular = false ;
2017-05-16 12:42:01 +02:00
if rescale_prediction_error_covariance
2017-04-27 10:44:27 +02:00
log_dF = log ( det ( F ./ ( sig * sig ' ) ) ) + 2 * sum ( log ( sig ) ) ;
iF = inv ( F ./ ( sig * sig ' ) ) ./ ( sig * sig ' ) ;
2020-01-03 13:19:57 +01:00
rescale_prediction_error_covariance = rescale_prediction_error_covariance0 ;
2017-04-27 10:44:27 +02:00
else
log_dF = log ( det ( F ) ) ;
iF = inv ( F ) ;
end
2015-10-13 17:15:01 +02:00
lik ( s ) = log_dF + transpose ( v ) * iF * v + length ( d_index ) * log ( 2 * pi ) ;
2011-09-19 17:01:24 +02:00
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 ) ;
2021-05-02 17:07:12 +02:00
if t > = no_more_missing_observations && ~ isqvec
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.
2012-05-01 14:07:03 +02:00
if t < = last
2017-04-27 10:44:27 +02:00
[ tmp , lik ( s + 1 : end ) ] = kalman_filter_ss ( Y , t , last , a , T , K , iF , log_dF , Z , pp , Zflag ) ;
2011-09-19 17:01:24 +02:00
end
2009-05-28 20:14:07 +02:00
2011-09-19 17:01:24 +02:00
% Compute minus the log-likelihood.
2012-01-22 22:40:46 +01:00
if presample > = diffuse_periods
LIK = sum ( lik ( 1 + presample - diffuse_periods : end ) ) ;
else
LIK = sum ( lik ) ;
end