stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

v4.1: Files related to kalman filter (likelihood evaluation of a state space model with or without missing observations). 

* Added a new tolerance parameter specific to the iteration on the riccati equation.
 * Added a kalman filter routine allowing for missing observations.
 * I do not distinguish anymore models with and without measurement errors (the same m file is used for both models to evaluate the likelihood). For a model without measurement errors H hat to be set to 0 scalar.

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2148 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2008-10-13 18:54:42 +00:00
parent 61e307d737
commit f58cdcf676
6 changed files with 283 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/dynare_config.m
View File

@ -38,6 +38,8 @@ end
dynareroot = strrep(which('dynare.m'),'dynare.m','');
addpath([dynareroot '/distributions/'])
addpath([dynareroot '/kalman/'])
% Add path to distribution-related function if under Matlab
% without the statistics toolbox
if ~exist('OCTAVE_VERSION') && isempty(ver('stats'))

1
matlab/global_initialization.m
View File

@ -145,6 +145,7 @@ function global_initialization()
options_.first_obs = 1;
options_.kalman_algo = 0;
options_.kalman_tol = 1e-12;
options_.riccati_tol = 1e-6;
options_.lik_algo = 1;
options_.lik_init = 1;
options_.load_mh_file = 0;

24
matlab/kalman/build_selection_matrix.m Normal file
View File

@ -0,0 +1,24 @@
function Z = build_selection_matrix(mf,m,p)
% Builds the selection matrix of the measurement equation from the vector
% of indices mf.
% Copyright (C) 2004-2008 Dynare Team
%
% 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/>.
Z = zeros(p,m);
for i=1:p
Z(i,mf(i)) = 1;
end

96
matlab/kalman/kalman_filter.m Normal file
View File

@ -0,0 +1,96 @@
function [LIK, lik] = kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
% Computes the likelihood of a stationnary state space model without measurement error, in the case of a non-singular var-cov matrix
%
% 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 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, likelihood
% lik [double] vector, density of observations in each period.
%
% 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.
% Copyright (C) 2004-2008 Dynare Team
%
% 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/>.
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,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
oldK = 0;
number_of_observations = smpl*pp;
lik(smpl+1) = number_of_observations*log(2*pi); % Log likelihood constant.
notsteady = 1; % Steady state flag.
F_singular = 1;
while notsteady & t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+QQ;
end
else
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;
end
notsteady = max(max(abs(K-oldK)))>riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
reste = smpl-t;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
lik(t) = lik(t) + reste*log(dF);
LIK = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);% Minus the log-likelihood.

110
matlab/kalman/missing_observations_kalman_filter.m Normal file
View File

@ -0,0 +1,110 @@
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)
% 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] saclar.
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% 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.
% Copyright (C) 2004-2008 Dynare Team
%
% 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/>.
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,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
oldK = 0;
lik(smpl+1) = number_of_observations*log(2*pi); % Log likelihood constant.
notsteady = 1; % Steady state flag.
F_singular = 1;
while notsteady & t<smpl
t = t+1;
if isempty(data_index{t})
a = T*a;
P = T*P*transpose(T)+QQ;
else
MF = mf(data_index{t}); % Set the selection for observed variables.
v = Y(data_index{t},t)-a(MF);
if ~isscalar(H) % => Errors in the measurement equation.
F = P(MF,MF) + H(data_index{t},data_index{t});
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;
end
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+QQ;
end
else
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;
end
notsteady = ~( (t>no_more_missing_observations) && (max(max(abs(K-oldK)))<riccati_tol) );
oldK = K;
end
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
reste = smpl-t;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
lik(t) = lik(t) + reste*log(dF);
LIK = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);% Minus the log-likelihood.

50
matlab/kalman/steady_state_kalman_gain.m Normal file
View File

@ -0,0 +1,50 @@
function [K,iF,P] = steady_state_kalman_gain(T,R,Q,H,mf)
% Given the invariant state space representation of a model, this
% function computes the gain matrix and the covariance matrix of the
% state vector at the steady state of the kalman filter.
%
% INPUTS
% T [double] m*m transition matrix of the state vector.
% R [double] m*q matrix (q is the number of structural innovations).
% Q [double] q*q covariance matrix of the structural innovations.
% H [double] p*p covariance matrix of the measurement error.
% mf [integer] p*1 vector, indices for the observed variables
%
% OUTPUTS
% K [double] kalman gain matrix.
% P [double] covariance matrix of the state vector.
%
% SPECIAL REQUIREMENTS
% Needs a solver for Riccati equations (dare.m)
% Copyright (C) 2004-2008 Dynare Team
%
% 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/>.
m = length(T);
p = length(mf);
Z = build_selection_matrix(mf,m,p);
if isempty(H)
H = zeros(p,p);
end
QQ = R*Q*transpose(R);
P = dare(T,transpose(Z),QQ,H);
iF = inv(Z*P*transpose(Z)+H);
K = T*P*transpose(Z)*iF;