Computes the likelihood of a stationnary state space model.
0001 function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P) 0002 % Computes the likelihood of a stationnary state space model. 0003 0004 %@info: 0005 %! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} DsgeLikelihood (@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods}) 0006 %! @anchor{kalman_filter} 0007 %! @sp 1 0008 %! Computes the likelihood of a stationary state space model, given initial condition for the states (mean and variance). 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item Y 0014 %! Matrix (@var{pp}*T) of doubles, data. 0015 %! @item start 0016 %! Integer scalar, first period. 0017 %! @item last 0018 %! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T). 0019 %! @item a 0020 %! Vector (@var{mm}*1) of doubles, initial mean of the state vector. 0021 %! @item P 0022 %! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector. 0023 %! @item kalman_tol 0024 %! Double scalar, tolerance parameter (rcond, inversibility of the covariance matrix of the prediction errors). 0025 %! @item riccati_tol 0026 %! Double scalar, tolerance parameter (iteration over the Riccati equation). 0027 %! @item presample 0028 %! Integer scalar, presampling if strictly positive (number of initial iterations to be discarded when evaluating the likelihood). 0029 %! @item T 0030 %! Matrix (@var{mm}*@var{mm}) of doubles, transition matrix of the state equation. 0031 %! @item Q 0032 %! Matrix (@var{rr}*@var{rr}) of doubles, covariance matrix of the structural innovations (noise in the state equation). 0033 %! @item R 0034 %! Matrix (@var{mm}*@var{rr}) of doubles, second matrix of the state equation relating the structural innovations to the state variables. 0035 %! @item H 0036 %! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar). 0037 %! @item Z 0038 %! Matrix (@var{pp}*@var{mm}) of doubles or vector of integers, matrix relating the states to the observed variables or vector of indices (depending on the value of @var{Zflag}). 0039 %! @item mm 0040 %! Integer scalar, number of state variables. 0041 %! @item pp 0042 %! Integer scalar, number of observed variables. 0043 %! @item rr 0044 %! Integer scalar, number of structural innovations. 0045 %! @item Zflag 0046 %! 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. 0047 %! @item diffuse_periods 0048 %! Integer scalar, number of diffuse filter periods in the initialization step. 0049 %! @end table 0050 %! @sp 2 0051 %! @strong{Outputs} 0052 %! @sp 1 0053 %! @table @ @var 0054 %! @item LIK 0055 %! Double scalar, value of (minus) the likelihood. 0056 %! @item likk 0057 %! Column vector of doubles, values of the density of each observation. 0058 %! @item a 0059 %! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample. 0060 %! @item P 0061 %! Matrix (@var{mm}*@var{mm}) of doubles, covariance of the state vector at the end of the (sub)sample. 0062 %! @end table 0063 %! @sp 2 0064 %! @strong{This function is called by:} 0065 %! @sp 1 0066 %! @ref{DsgeLikelihood} 0067 %! @sp 2 0068 %! @strong{This function calls:} 0069 %! @sp 1 0070 %! @ref{kalman_filter_ss} 0071 %! @end deftypefn 0072 %@eod: 0073 0074 % Copyright (C) 2004-2011 Dynare Team 0075 % 0076 % This file is part of Dynare. 0077 % 0078 % Dynare is free software: you can redistribute it and/or modify 0079 % it under the terms of the GNU General Public License as published by 0080 % the Free Software Foundation, either version 3 of the License, or 0081 % (at your option) any later version. 0082 % 0083 % Dynare is distributed in the hope that it will be useful, 0084 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0085 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0086 % GNU General Public License for more details. 0087 % 0088 % You should have received a copy of the GNU General Public License 0089 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0090 0091 % AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr 0092 0093 % Set defaults. 0094 if nargin<17 0095 Zflag = 0; 0096 end 0097 0098 if nargin<18 0099 diffuse_periods = 0; 0100 end 0101 0102 if nargin<19 0103 analytic_derivation = 0; 0104 end 0105 0106 if isempty(Zflag) 0107 Zflag = 0; 0108 end 0109 0110 if isempty(diffuse_periods) 0111 diffuse_periods = 0; 0112 end 0113 0114 % Get sample size. 0115 smpl = last-start+1; 0116 0117 % Initialize some variables. 0118 dF = 1; 0119 QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 0120 t = start; % Initialization of the time index. 0121 likk = zeros(smpl,1); % Initialization of the vector gathering the densities. 0122 LIK = Inf; % Default value of the log likelihood. 0123 oldK = Inf; 0124 notsteady = 1; 0125 F_singular = 1; 0126 0127 if analytic_derivation == 0, 0128 DLIK=[]; 0129 Hess=[]; 0130 else 0131 k = size(DT,3); % number of structural parameters 0132 DLIK = zeros(k,1); % Initialization of the score. 0133 Da = zeros(mm,k); % Derivative State vector. 0134 0135 if Zflag==0, 0136 C = zeros(pp,mm); 0137 for ii=1:pp; C(ii,Z(ii))=1;end % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT) 0138 else 0139 C=Z; 0140 end 0141 dC = zeros(pp,mm,k); % either selection matrix or schur have zero derivatives 0142 if analytic_derivation==2, 0143 Hess = zeros(k,k); % Initialization of the Hessian 0144 D2a = zeros(mm,k,k); % State vector. 0145 d2C = zeros(pp,mm,k,k); 0146 else 0147 Hess=[]; 0148 D2a=[]; 0149 D2T=[]; 0150 D2Yss=[]; 0151 end 0152 LIK={inf,DLIK,Hess}; 0153 end 0154 0155 while notsteady && t<=last 0156 s = t-start+1; 0157 if Zflag 0158 v = Y(:,t)-Z*a; 0159 F = Z*P*Z' + H; 0160 else 0161 v = Y(:,t)-a(Z); 0162 F = P(Z,Z) + H; 0163 end 0164 if rcond(F) < kalman_tol 0165 if ~all(abs(F(:))<kalman_tol) 0166 return 0167 else 0168 a = T*a; 0169 P = T*P*transpose(T)+QQ; 0170 end 0171 else 0172 F_singular = 0; 0173 dF = det(F); 0174 iF = inv(F); 0175 likk(s) = log(dF)+transpose(v)*iF*v; 0176 if Zflag 0177 K = P*Z'*iF; 0178 Ptmp = T*(P-K*Z*P)*transpose(T)+QQ; 0179 else 0180 K = P(:,Z)*iF; 0181 Ptmp = T*(P-K*P(Z,:))*transpose(T)+QQ; 0182 end 0183 tmp = (a+K*v); 0184 if analytic_derivation, 0185 if analytic_derivation==2, 0186 [Da,DP,DLIKt,D2a,D2P, Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P); 0187 else 0188 [Da,DP,DLIKt] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady); 0189 end 0190 if t>presample 0191 DLIK = DLIK + DLIKt; 0192 if analytic_derivation==2, 0193 Hess = Hess + Hesst; 0194 end 0195 end 0196 end 0197 a = T*tmp; 0198 P=Ptmp; 0199 notsteady = max(abs(K(:)-oldK))>riccati_tol; 0200 oldK = K(:); 0201 end 0202 t = t+1; 0203 end 0204 0205 if F_singular 0206 error('The variance of the forecast error remains singular until the end of the sample') 0207 end 0208 0209 % Add observation's densities constants and divide by two. 0210 likk(1:s) = .5*(likk(1:s) + pp*log(2*pi)); 0211 0212 % Call steady state Kalman filter if needed. 0213 if t <= last 0214 if analytic_derivation, 0215 if analytic_derivation==2, 0216 [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ... 0217 analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss); 0218 else 0219 [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ... 0220 analytic_derivation,Da,DT,DYss); 0221 end 0222 DLIK = DLIK + tmp{2}; 0223 if analytic_derivation==2, 0224 Hess = Hess + tmp{3}; 0225 end 0226 else 0227 [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag); 0228 end 0229 end 0230 0231 % Compute minus the log-likelihood. 0232 if presample 0233 if presample>=diffuse_periods 0234 likk = likk(1+(presample-diffuse_periods):end); 0235 end 0236 end 0237 LIK = sum(likk); 0238 0239 if analytic_derivation, 0240 DLIK = DLIK/2; 0241 if analytic_derivation==2, 0242 Hess = Hess + tril(Hess,-1)'; 0243 Hess = -Hess/2; 0244 LIK={LIK, DLIK, Hess}; 0245 else 0246 LIK={LIK, DLIK}; 0247 end 0248 end