1function [LIK, LIKK, a, P] = kalman_filter(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,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P) 2% Computes the likelihood of a stationary state space model. 3 4%@info: 5%! @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}) 6%! @anchor{kalman_filter} 7%! @sp 1 8%! Computes the likelihood of a stationary state space model, given initial condition for the states (mean and variance). 9%! @sp 2 10%! @strong{Inputs} 11%! @sp 1 12%! @table @ @var 13%! @item Y 14%! Matrix (@var{pp}*T) of doubles, data. 15%! @item start 16%! Integer scalar, first period. 17%! @item last 18%! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T). 19%! @item a 20%! Vector (@var{mm}*1) of doubles, initial mean of the state vector. 21%! @item P 22%! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector. 23%! @item kalman_tol 24%! Double scalar, tolerance parameter (rcond, inversibility of the covariance matrix of the prediction errors). 25%! @item riccati_tol 26%! Double scalar, tolerance parameter (iteration over the Riccati equation). 27%! @item presample 28%! Integer scalar, presampling if strictly positive (number of initial iterations to be discarded when evaluating the likelihood). 29%! @item T 30%! Matrix (@var{mm}*@var{mm}) of doubles, transition matrix of the state equation. 31%! @item Q 32%! Matrix (@var{rr}*@var{rr}) of doubles, covariance matrix of the structural innovations (noise in the state equation). 33%! @item R 34%! Matrix (@var{mm}*@var{rr}) of doubles, second matrix of the state equation relating the structural innovations to the state variables. 35%! @item H 36%! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar). 37%! @item Z 38%! 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}). 39%! @item mm 40%! Integer scalar, number of state variables. 41%! @item pp 42%! Integer scalar, number of observed variables. 43%! @item rr 44%! Integer scalar, number of structural innovations. 45%! @item Zflag 46%! 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. 47%! @item diffuse_periods 48%! Integer scalar, number of diffuse filter periods in the initialization step. 49%! @end table 50%! @sp 2 51%! @strong{Outputs} 52%! @sp 1 53%! @table @ @var 54%! @item LIK 55%! Double scalar, value of (minus) the likelihood. 56%! @item LIKK 57%! Column vector of doubles, values of the density of each observation. 58%! @item a 59%! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample. 60%! @item P 61%! Matrix (@var{mm}*@var{mm}) of doubles, covariance of the state vector at the end of the (sub)sample. 62%! @end table 63%! @sp 2 64%! @strong{This function is called by:} 65%! @sp 1 66%! @ref{DsgeLikelihood} 67%! @sp 2 68%! @strong{This function calls:} 69%! @sp 1 70%! @ref{kalman_filter_ss} 71%! @end deftypefn 72%@eod: 73 74% Copyright (C) 2004-2017 Dynare Team 75% 76% This file is part of Dynare. 77% 78% Dynare is free software: you can redistribute it and/or modify 79% it under the terms of the GNU General Public License as published by 80% the Free Software Foundation, either version 3 of the License, or 81% (at your option) any later version. 82% 83% Dynare is distributed in the hope that it will be useful, 84% but WITHOUT ANY WARRANTY; without even the implied warranty of 85% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 86% GNU General Public License for more details. 87% 88% You should have received a copy of the GNU General Public License 89% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 90 91% Set defaults. 92if nargin<17 93 Zflag = 0; 94end 95 96if nargin<18 97 diffuse_periods = 0; 98end 99 100if nargin<19 101 analytic_derivation = 0; 102end 103 104if isempty(Zflag) 105 Zflag = 0; 106end 107 108if isempty(diffuse_periods) 109 diffuse_periods = 0; 110end 111 112% Get sample size. 113smpl = last-start+1; 114 115% Initialize some variables. 116dF = 1; 117QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations. 118t = start; % Initialization of the time index. 119likk = zeros(smpl,1); % Initialization of the vector gathering the densities. 120LIK = Inf; % Default value of the log likelihood. 121oldK = Inf; 122notsteady = 1; 123F_singular = true; 124asy_hess=0; 125 126if analytic_derivation == 0 127 DLIK=[]; 128 Hess=[]; 129 LIKK=[]; 130else 131 k = size(DT,3); % number of structural parameters 132 DLIK = zeros(k,1); % Initialization of the score. 133 Da = zeros(mm,k); % Derivative State vector. 134 dlikk = zeros(smpl,k); 135 136 if Zflag==0 137 C = zeros(pp,mm); 138 for ii=1:pp, C(ii,Z(ii))=1; end % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT) 139 else 140 C=Z; 141 end 142 dC = zeros(pp,mm,k); % either selection matrix or schur have zero derivatives 143 if analytic_derivation==2 144 Hess = zeros(k,k); % Initialization of the Hessian 145 D2a = zeros(mm,k,k); % State vector. 146 d2C = zeros(pp,mm,k,k); 147 else 148 asy_hess=D2T; 149 Hess=[]; 150 D2a=[]; 151 D2T=[]; 152 D2Yss=[]; 153 end 154 if asy_hess 155 Hess = zeros(k,k); % Initialization of the Hessian 156 end 157 LIK={inf,DLIK,Hess}; 158 LIKK={likk,dlikk}; 159end 160 161rescale_prediction_error_covariance0=rescale_prediction_error_covariance; 162while notsteady && t<=last 163 s = t-start+1; 164 if Zflag 165 v = Y(:,t)-Z*a; 166 F = Z*P*Z' + H; 167 else 168 v = Y(:,t)-a(Z); 169 F = P(Z,Z) + H; 170 end 171 badly_conditioned_F = false; 172 if rescale_prediction_error_covariance 173 sig=sqrt(diag(F)); 174 if any(diag(F)<kalman_tol) || rcond(F./(sig*sig'))<kalman_tol 175 badly_conditioned_F = true; 176 end 177 else 178 if rcond(F)<kalman_tol 179 sig=sqrt(diag(F)); 180 if any(diag(F)<kalman_tol) || rcond(F./(sig*sig'))<kalman_tol 181 badly_conditioned_F = true; 182 else 183 rescale_prediction_error_covariance=1; 184 end 185 end 186 end 187 if badly_conditioned_F 188 if ~all(abs(F(:))<kalman_tol) 189 % Use univariate filter (will remove observations with zero variance prediction error) 190 return 191 else 192 % Pathological case, discard draw. 193 return 194 end 195 else 196 F_singular = false; 197 if rescale_prediction_error_covariance 198 log_dF = log(det(F./(sig*sig')))+2*sum(log(sig)); 199 iF = inv(F./(sig*sig'))./(sig*sig'); 200 rescale_prediction_error_covariance=rescale_prediction_error_covariance0; 201 else 202 log_dF = log(det(F)); 203 iF = inv(F); 204 end 205 likk(s) = log_dF+transpose(v)*iF*v; 206 if Zflag 207 K = P*Z'*iF; 208 Ptmp = T*(P-K*Z*P)*transpose(T)+QQ; 209 else 210 K = P(:,Z)*iF; 211 Ptmp = T*(P-K*P(Z,:))*transpose(T)+QQ; 212 end 213 tmp = (a+K*v); 214 if analytic_derivation 215 if analytic_derivation==2 216 [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); 217 else 218 [Da,DP,DLIKt,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady); 219 end 220 if t>presample 221 DLIK = DLIK + DLIKt; 222 if analytic_derivation==2 || asy_hess 223 Hess = Hess + Hesst; 224 end 225 end 226 dlikk(s,:)=DLIKt; 227 end 228 a = T*tmp; 229 P = Ptmp; 230 notsteady = max(abs(K(:)-oldK))>riccati_tol; 231 oldK = K(:); 232 end 233 t = t+1; 234end 235 236if F_singular 237 error('The variance of the forecast error remains singular until the end of the sample') 238end 239 240% Add observation's densities constants and divide by two. 241likk(1:s) = .5*(likk(1:s) + pp*log(2*pi)); 242if analytic_derivation 243 DLIK = DLIK/2; 244 dlikk = dlikk/2; 245 if analytic_derivation==2 || asy_hess 246 if asy_hess==0 247 Hess = Hess + tril(Hess,-1)'; 248 end 249 Hess = -Hess/2; 250 end 251end 252 253% Call steady state Kalman filter if needed. 254if t <= last 255 if analytic_derivation 256 if analytic_derivation==2 257 [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, log_dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, D2a, D2T, D2Yss); 258 else 259 [tmp, tmp2] = kalman_filter_ss(Y, t, last, a, T, K, iF, log_dF, Z, pp, Zflag, analytic_derivation, Da, DT, DYss, asy_hess); 260 end 261 likk(s+1:end) = tmp2{1}; 262 dlikk(s+1:end,:) = tmp2{2}; 263 DLIK = DLIK + tmp{2}; 264 if analytic_derivation==2 || asy_hess 265 Hess = Hess + tmp{3}; 266 end 267 else 268 [tmp, likk(s+1:end)] = kalman_filter_ss(Y, t, last, a, T, K, iF, log_dF, Z, pp, Zflag); 269 end 270end 271 272% Compute minus the log-likelihood. 273if presample>diffuse_periods 274 LIK = sum(likk(1+(presample-diffuse_periods):end)); 275else 276 LIK = sum(likk); 277end 278 279if analytic_derivation 280 if analytic_derivation==2 || asy_hess 281 LIK={LIK, DLIK, Hess}; 282 else 283 LIK={LIK, DLIK}; 284 end 285 LIKK={likk, dlikk}; 286else 287 LIKK=likk; 288end 289