Commit f2ca6d0a authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Changed kalman filter routines to allow for arbitrary initial conditions...

Changed kalman filter routines to allow for arbitrary initial conditions (needed for the introduction of breaks on the estimated
parameters and also for the estimation of the initial states).

Added specialized routines for steady state  kalman filter.

Completed header of DsgeLikelihood (missing refs to the routines called by DsgeLikelihood).
parent 9c22fc1b
......@@ -103,7 +103,7 @@ function [fval,cost_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,Dynar
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resol}
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{schur_statespace_transformation}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @end deftypefn
%@eod:
......
function [LIK, lik] = kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
function [LIK, lik, 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)
% Computes the likelihood of a stationnary state space model.
%
% INPUTS
% INPUTS
% Y [double] pp*smpl matrix of data.
% start [integer] scalar, index of the first observation (column of Y).
% last [integer] scalar, index of the last observation (column of Y).
% a [double] mm*1 vector, initial level of the state vector.
% P [double] mm*mm 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.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix with stationary variables
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% mf [integer] pp*1 vector of indices.
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
%
% OUTPUTS
% LIK [double] scalar, MINUS loglikelihood
% lik [double] vector, density of observations in each period.
% Z [integer] pp*1 vector of indices for the observed variables, if Zflag=0, pp*mm matrix if Zflag>0.
% mm [integer] scalar, dimension of the state vector.
% pp [integer] scalar, number of observed variables.
% rr [integer] scalar, number of structural innovations.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
% OUTPUTS
% LIK [double]
% lik [double] (last-start+1)*1 vector, density of observations in each periods.
% a [double] mm*1 vector, estimated level of the states.
% P [double] mm*mm matrix, covariance matrix of the states.
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
%
% Copyright (C) 2004-2011 Dynare Team
%
% stephane DOT adjemian AT ens DOT fr
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
......@@ -42,23 +47,46 @@ function [LIK, lik] = kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
smpl = size(Y,2); % Sample size.
mm = size(T,2); % Number of state variables.
pp = size(Y,1); % Maximum number of observed variables.
a = zeros(mm,1); % State vector.
dF = 1; % det(F).
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
% Set defaults.
if nargin<17
Zflag = 0;
diffuse_periods = 0;
end
if nargin<18
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.
oldK = Inf;
notsteady = 1; % Steady state flag.
notsteady = 1;
F_singular = 1;
while notsteady && t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
while notsteady && t<=last
s = t-start+1;
if Zflag
v = Y(:,t)-Z*a;
F = Z*P*Z' + H;
else
v = Y(:,t)-a(Z);
F = P(Z,Z) + H;
end
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
......@@ -70,31 +98,37 @@ while notsteady && t<smpl
F_singular = 0;
dF = det(F);
iF = inv(F);
lik(t) = log(dF)+transpose(v)*iF*v;
K = P(:,mf)*iF;
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+QQ;
notsteady = max(abs(K(:)-oldK)) > riccati_tol;
lik(s) = log(dF)+transpose(v)*iF*v;
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);
notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
end
t = t+1;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
lik(t0:smpl) = lik(t0:smpl) + log(dF);
end
% Add observation's densities constants and divide by two.
lik(1:s) = .5*(lik(1:s) + pp*log(2*pi));
% adding log-likelihhod constants
lik = (lik + pp*log(2*pi))/2;
% 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
LIK = sum(lik(start:end)); % Minus the log-likelihood.
\ No newline at end of file
% Compute minus the log-likelihood.
if presample
if presample>=diffuse_periods
lik = lik(1+(presample-diffuse_periods):end);
end
end
LIK = sum(lik);
\ No newline at end of file
function [LIK, lik] = diffuse_kalman_filter(T,R,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol,riccati_tol)
function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, kalman_tol, riccati_tol, presample, T, R, Q, H, Z, mm, pp, rr)
% Computes the diffuse likelihood of a state space model.
%
% INPUTS
% T [double] mm*mm transition matrix
% R [double] mm*rr matrix
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp covariance matrix of the measurement errors (if H is equal to zero (scalar) there is no measurement error).
% INPUTS
% Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables.
% start [integer] scalar, first observation.
% last [integer] scalar, last observation.
% a [double] mm*1 vector, levels of the state variables.
% Pinf [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% Pstar [double] mm*mm matrix used to initialize the covariance matrix of the state vector.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z [double] pp*mm matrix, selection matrix or pp linear independant combinations of the 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 matrix, transition matrix in the state equations.
% R [double] mm*rr matrix relating the structural innovations to the state vector.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp covariance matrix of the measurement errors (if H is equal to zero (scalar) there is no measurement error).
% Z [double] pp*mm matrix, selection matrix or pp linear independant combinations of the state vector.
% mm [integer] scalar, number of state variables.
% pp [integer] scalar, number of observed variables.
% rr [integer] scalar, number of structural innovations.
%
% OUTPUTS
% LIK: MINUS loglikelihood
% lik: density vector in each period
% OUTPUTS
% LIK [double] scalar, minus loglikelihood
% lik [double] smpl*1 vector, log density of each vector of observations.
% a [double] mm*1 vector, current estimate of the state vector.
% Pstar [double] mm*mm matrix, covariance matrix of the state vector.
%
% REFERENCES
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
% Copyright (C) 2004-2010 Dynare Team
% Copyright (C) 2004-2011 Dynare Team
%
% This file is part of Dynare.
%
......@@ -40,26 +48,24 @@ function [LIK, lik] = diffuse_kalman_filter(T,R,Q,H,Pinf,Pstar,Y,start,Z,kalman_
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
[pp,smpl] = size(Y);
mm = size(T,2);
a = zeros(mm,1);
% Get sample size.
smpl = last-start+1;
% Initialize some variables.
dF = 1;
QQ = R*Q*transpose(R);
t = 0;
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = start; % Initialization of the time index.
dlik = zeros(smpl,1); % Initialization of the vector gathering the densities.
dLIK = Inf; % Default value of the log likelihood.
oldK = Inf;
lik = zeros(smpl,1);
LIK = Inf;
% lik(smpl+1) = smpl*pp*log(2*pi);
notsteady = 1;
reste = 0;
while rank(Pinf,kalman_tol) && (t<smpl)
t = t+1;
while rank(Pinf,kalman_tol) && (t<=last)
s = t-start+1;
v = Y(:,t)-Z*a;
Finf = Z*Pinf*Z' ;
Finf = Z*Pinf*Z';
if rcond(Finf) < kalman_tol
if ~all(abs(Finf(:)) < kalman_tol)
% The univariate diffuse kalman filter should be used.
% The univariate diffuse kalman filter should be used instead.
return
else
Fstar = Z*Pstar*Z' + H;
......@@ -76,72 +82,42 @@ while rank(Pinf,kalman_tol) && (t<smpl)
iFstar = inv(Fstar);
dFstar = det(Fstar);
Kstar = Pstar*Z'*iFstar;
lik(t) = log(dFstar) + v'*iFstar*v;
dlik(s)= log(dFstar) + v'*iFstar*v;
Pinf = T*Pinf*transpose(T);
Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ;
a = T*(a+Kstar*v);
a = T*(a+Kstar*v);
end
end
else
lik(t) = log(det(Finf));
iFinf = inv(Finf);
Kinf = Pinf*Z'*iFinf;
Fstar = Z*Pstar*Z' + H;
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf;
Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ;
Pinf = T*(Pinf-Pinf*Z'*Kinf')*T';
a = T*(a+Kinf*v);
dlik(s)= log(det(Finf));
iFinf = inv(Finf);
Kinf = Pinf*Z'*iFinf;
Fstar = Z*Pstar*Z' + H;
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf;
Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ;
Pinf = T*(Pinf-Pinf*Z'*Kinf')*T';
a = T*(a+Kinf*v);
end
t = t+1;
end
if t == smpl
if t>last
warning(['There isn''t enough information to estimate the initial conditions of the nonstationary variables']);
LIK = NaN;
dLIK = NaN;
return
end
F_singular = 1;
while notsteady && (t<smpl)
t = t+1;
v = Y(:,t)-Z*a;
F = Z*Pstar*Z' + H;
dF = det(F);
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
Pstar = T*Pstar*T'+QQ;
end
else
F_singular = 0;
iF = inv(F);
lik(t) = log(dF)+v'*iF*v;
K = Pstar*Z'*iF;
a = T*(a+K*v);
Pstar = T*(Pstar-K*Z*Pstar)*T'+QQ;
notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
end
end
if F_singular == 1
warning(['The variance of the forecast error remains singular until the end of the sample'])
LIK = NaN;
return
end
dlik = dlik(1:s);
dlik = .5*(dlik + pp*log(2*pi));
if t < smpl
t0 = t+1;
while t<smpl
t = t+1;
v = Y(:,t)-Z*a;
a = T*(a+K*v);
lik(t) = v'*iF*v;
if presample
if presample>=length(dlik)
dLIK = 0;
dlik = [];
else
dlik = dlik(1+presample:end);
dLIK = sum(dlik);% Minus the log-likelihood (for the initial periods).
end
lik(t0:smpl) = lik(t0:smpl) + log(dF);
end
lik = (lik + pp*log(2*pi))/2;
LIK = sum(lik(start:end)); % Minus the log-likelihood.
\ No newline at end of file
else
dLIK = sum(dlik);
end
\ No newline at end of file
function [LIK, lik, a] = kalman_filter_ss(Y,start,last,a,T,K,iF,dF,Z,pp,Zflag)
% Computes the likelihood of a stationnary state space model (steady state kalman filter).
%
% INPUTS
% Y [double] pp*smpl matrix of data.
% start [integer] scalar, index of the first observation (column of Y).
% last [integer] scalar, index of the last observation (column of Y).
% a [double] mm*1 vector, initial level of the state vector.
% P [double] mm*mm matrix, covariance matrix of the initial state vector.
% T [double] mm*mm transition matrix of the state equation.
% K [double] mm*pp matrix, steady state kalman gain.
% iF [double] pp*pp matrix, inverse of the steady state covariance matrix of the predicted errors.
% dF [double] scalar, determinant of the steady state covariance matrix of the predicted errors.
% Z [integer] pp*1 vector of indices for the observed variables, if Zflag=0.
% pp [integer] scalar, number of observed variables.
%
%
% OUTPUTS
% LIK [double] scalar, minus log likelihood.
% lik [double] (last-start+1)*1 vector, density of each observation.
% a [double] mm*1 vector, estimate of the state vector.
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2011 Dynare Team
% stephane DOT adjemian AT ens DOT fr
%
% 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/>.
% Get sample size.
smpl = last-start+1;
% Initialize some variables.
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.
while t <= last
if Zflag
v = Y(:,t)-Z*a;
else
v = Y(:,t)-a(Z);
end
a = T*(a+K*v);
lik(t-start+1) = transpose(v)*iF*v;
t = t+1;
end
% Adding constant determinant of F (prediction error covariance matrix)
lik = lik + log(dF);
% Add log-likelihhod constants and divide by two
lik = .5*(lik + pp*log(2*pi));
% Sum the observation's densities (minus the likelihood)
LIK = sum(lik);
\ No newline at end of file
function [LIK, lik] = missing_observations_kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations)
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)
% Computes the likelihood of a state space model in the case with missing observations.
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix of the initial state vector.
% Y [double] pp*smpl matrix of detrended data, where pp is the maximum number of observed variables
% trend [double] pp*smpl matrix
% start [integer] scalar, likelihood evaluation starts at 'start'.
% mf [integer] pp*1 vector of indices.
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
% data_index [cell] 1*smpl cell of column vectors of indices.
% number_of_observations [integer] scalar.
% no_more_missing_observations [integer] scalar.
% 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.
%
% OUTPUTS
% OUTPUTS
% LIK [double] scalar, MINUS loglikelihood
% lik [double] vector, density of observations in each period.
% a [double] mm*1 vector, estimated level of the states.
% P [double] mm*mm matrix, covariance matrix of the states.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
......@@ -46,34 +49,53 @@ function [LIK, lik] = missing_observations_kalman_filter(T,R,Q,H,P,Y,start,mf,k
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
smpl = size(Y,2); % Sample size.
mm = size(T,2); % Number of state variables.
pp = size(Y,1); % Maximum number of observed variables.
a = zeros(mm,1); % State vector.
dF = 1; % det(F).
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
% 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.
oldK = Inf;
notsteady = 1; % Steady state flag.
notsteady = 1;
F_singular = 1;
while notsteady && t<smpl
t = t+1;
while notsteady & t<=last
s = t-start+1;
d_index = data_index{t};
if isempty(d_index)
a = T*a;
P = T*P*transpose(T)+QQ;
else
MF = mf(d_index); % Set the selection for observed variables.
v = Y(d_index,t)-a(MF);
if ~isscalar(H) % => Errors in the measurement equation.
F = P(MF,MF) + H(d_index,d_index);
else% =>
% case 1. No errors in the measurement (H=0) and more than one variable is observed in this state space model.
% case 2. Errors in the measurement equation, but only one variable is observed in this state-space model.
F = P(MF,MF)+H;
% 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);
end
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
......@@ -86,33 +108,40 @@ while notsteady && t<smpl
F_singular = 0;
dF = det(F);
iF = inv(F);
lik(t) = log(dF) + transpose(v)*iF*v + length(d_index)*log(2*pi);
K = P(:,MF)*iF;
a = T*(a+K*v);
P = T*(P-K*P(MF,:))*transpose(T)+QQ;
if t>no_more_missing_observations
lik(s) = log(dF) + transpose(v)*iF*v + length(d_index)*log(2*pi);
if Zflag