Computes the likelihood of a state space model (initialization with diffuse steps, univariate approach).
0001 function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index, number_of_observations, no_more_missing_observations, Y, start, last, a, Pinf, Pstar, kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr) 0002 % Computes the likelihood of a state space model (initialization with diffuse steps, univariate approach). 0003 0004 %@info: 0005 %! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{Pstar}, @var{llik} ] =} univariate_kalman_filter_d (@var{data_index}, @var{number_of_observations},@var{no_more_missing_observations}, @var{Y}, @var{start}, @var{last}, @var{a}, @var{Pinf}, @var{Pstar}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{R},@var{Q},@var{H},@var{Z},@var{mm},@var{pp},@var{rr}) 0006 %! @anchor{univariate_kalman_filter_d} 0007 %! @sp 1 0008 %! Computes the likelihood of a state space model (initialization with diffuse steps, univariate approach). 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item data_index 0014 %! Matlab's cell, 1*T cell of column vectors of indices (in the vector of observed variables). 0015 %! @item number_of_observations 0016 %! Integer scalar, effective number of observations. 0017 %! @item no_more_missing_observations 0018 %! Integer scalar, date after which there is no more missing observation (it is then possible to switch to the steady state kalman filter). 0019 %! @item Y 0020 %! Matrix (@var{pp}*T) of doubles, data. 0021 %! @item start 0022 %! Integer scalar, first period. 0023 %! @item last 0024 %! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T). 0025 %! @item a 0026 %! Vector (@var{mm}*1) of doubles, initial mean of the state vector. 0027 %! @item Pinf 0028 %! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector (non stationary part). 0029 %! @item Pstar 0030 %! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector (stationary part). 0031 %! @item kalman_tol 0032 %! Double scalar, tolerance parameter (rcond, inversibility of the covariance matrix of the prediction errors). 0033 %! @item riccati_tol 0034 %! Double scalar, tolerance parameter (iteration over the Riccati equation). 0035 %! @item presample 0036 %! Integer scalar, presampling if strictly positive (number of initial iterations to be discarded when evaluating the likelihood). 0037 %! @item T 0038 %! Matrix (@var{mm}*@var{mm}) of doubles, transition matrix of the state equation. 0039 %! @item R 0040 %! Matrix (@var{mm}*@var{rr}) of doubles, second matrix of the state equation relating the structural innovations to the state variables. 0041 %! @item Q 0042 %! Matrix (@var{rr}*@var{rr}) of doubles, covariance matrix of the structural innovations (noise in the state equation). 0043 %! @item H 0044 %! Vector (@var{pp}) of doubles, diagonal of covariance matrix of the measurement errors (corelation among measurement errors is handled by a model transformation). 0045 %! @item Z 0046 %! Matrix (@var{pp}*@var{mm}) of doubles, matrix relating the states to the observed variables. 0047 %! @item mm 0048 %! Integer scalar, number of state variables. 0049 %! @item pp 0050 %! Integer scalar, number of observed variables. 0051 %! @item rr 0052 %! Integer scalar, number of structural innovations. 0053 %! @item Zflag 0054 %! Integer scalar, equal to 0 if Z is a vector of indices targeting the obseved variables in the state vector, equal to 1 if Z is a @var{pp}*@var{mm} matrix. 0055 %! @item diffuse_periods 0056 %! Integer scalar, number of diffuse filter periods in the initialization step. 0057 %! @end table 0058 %! @sp 2 0059 %! @strong{Outputs} 0060 %! @sp 1 0061 %! @table @ @var 0062 %! @item dLIK 0063 %! Double scalar, value of (minus) the likelihood for the first d observations. 0064 %! @item likk 0065 %! Column vector (d*1) of doubles, values of the density of each (fist) observations. 0066 %! @item a 0067 %! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample. 0068 %! @item Pstar 0069 %! Matrix (@var{mm}*@var{mm}) of doubles, covariance of the state vector at the end of the subsample. 0070 %! @item llik 0071 %! Matrix of doubles (d*@var{pp}) used by @ref{DsgeLikelihood_hh} 0072 %! @end table 0073 %! @sp 2 0074 %! @strong{This function is called by:} 0075 %! @sp 1 0076 %! @ref{dsge_likelihood}, @ref{DsgeLikelihood_hh} 0077 %! @sp 2 0078 %! @strong{This function calls:} 0079 %! @sp 1 0080 %! @ref{univariate_kalman_filter_ss} 0081 %! @end deftypefn 0082 %@eod: 0083 0084 % Copyright (C) 2004-2011 Dynare Team 0085 % 0086 % This file is part of Dynare. 0087 % 0088 % Dynare is free software: you can redistribute it and/or modify 0089 % it under the terms of the GNU General Public License as published by 0090 % the Free Software Foundation, either version 3 of the License, or 0091 % (at your option) any later version. 0092 % 0093 % Dynare is distributed in the hope that it will be useful, 0094 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0095 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0096 % GNU General Public License for more details. 0097 % 0098 % You should have received a copy of the GNU General Public License 0099 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0100 0101 % Get sample size. 0102 smpl = last-start+1; 0103 0104 % Initialize some variables. 0105 dF = 1; 0106 QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 0107 t = start; % Initialization of the time index. 0108 dlikk= zeros(smpl,1); % Initialization of the vector gathering the densities. 0109 dLIK = Inf; % Default value of the log likelihood. 0110 oldK = Inf; 0111 llik = NaN(smpl,pp); 0112 0113 newRank = rank(Pinf,kalman_tol); 0114 l2pi = log(2*pi); 0115 0116 while newRank && (t<=last) 0117 s = t-start+1; 0118 d_index = data_index{t}; 0119 for i=1:length(d_index) 0120 Zi = Z(d_index(i),:); 0121 prediction_error = Y(d_index(i),t) - Zi*a; 0122 Fstar = Zi*Pstar*Zi' + H(d_index(i)); 0123 Finf = Zi*Pinf*Zi'; 0124 Kstar = Pstar*Zi'; 0125 if Finf>kalman_tol && newRank 0126 Kinf = Pinf*Zi'; 0127 Kinf_Finf = Kinf/Finf; 0128 a = a + Kinf_Finf*prediction_error; 0129 Pstar = Pstar + Kinf*(Kinf_Finf'*(Fstar/Finf)) - Kstar*Kinf_Finf' - Kinf_Finf*Kstar'; 0130 Pinf = Pinf - Kinf*Kinf_Finf'; 0131 llik(s,d_index(i)) = log(Finf) + l2pi; 0132 dlikk(s) = dlikk(s) + llik(s,d_index(i)); 0133 elseif Fstar>kalman_tol 0134 llik(s,d_index(i)) = log(Fstar) + prediction_error*prediction_error/Fstar + l2pi; 0135 dlikk(s) = dlikk(s) + llik(s,d_index(i)); 0136 a = a+Kstar*(prediction_error/Fstar); 0137 Pstar = Pstar-Kstar*(Kstar'/Fstar); 0138 end 0139 end 0140 if newRank 0141 oldRank = rank(Pinf,kalman_tol); 0142 else 0143 oldRank = 0; 0144 end 0145 a = T*a; 0146 Pstar = T*Pstar*T'+QQ; 0147 Pinf = T*Pinf*T'; 0148 if newRank 0149 newRank = rank(Pinf,kalman_tol); 0150 end 0151 if oldRank ~= newRank 0152 disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!') 0153 end 0154 t = t+1; 0155 end 0156 0157 if (t>last) 0158 error(['univariate_diffuse_kalman_filter:: There isn''t enough information to estimate the initial conditions of the nonstationary variables']); 0159 LIK = NaN; 0160 return 0161 end 0162 0163 % Divide by two. 0164 dlikk = .5*dlikk(1:s); 0165 llik = .5*llik(1:s,:); 0166 0167 dLIK = sum(dlikk(1+presample:end));