Skip to content
GitLab
Commit 7668bc4a authored by Michel Juillard's avatar Michel Juillard
Browse files

-new smoother function kalman/smoother/kalman_smoother.m 

-fixing bugs in dynare_estimation_1.m
parent b261eb0b
Hide whitespace changes
Inline Side-by-side
matlab/DsgeSmoother.m
View file @ 7668bc4a
......@@ -30,7 +30,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2006-2009 Dynare Team
% Copyright (C) 2006-2010 Dynare Team
%
% This file is part of Dynare.
%
......@@ -70,8 +70,9 @@ set_all_parameters(xparam1);
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
[T,R,SteadyState] = dynare_resolve;
bayestopt_.mf = bayestopt_.mf2;
[T,R,SteadyState] = dynare_resolve(bayestopt_.smoother_var_list,...
bayestopt_.smoother_restrict_columns,[]);
bayestopt_.mf = bayestopt_.smoother_mf;
if options_.noconstant
constant = zeros(nobs,1);
else
......@@ -96,7 +97,7 @@ else
end
start = options_.presample+1;
np = size(T,1);
mf = bayestopt_.mf;
mf = bayestopt_.smoother_mf;
% ------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
% ------------------------------------------------------------------------------
......@@ -184,13 +185,16 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
Pinf(1:nk,1:nk) = diag(dd);
end
end
kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol;
data1 = Y-trend;
% -----------------------------------------------------------------------------
% 4. Kalman smoother
% -----------------------------------------------------------------------------
if any(any(H ~= 0)) % should be replaced by a flag
if kalman_algo == 1
[alphahat,epsilonhat,etahat,ahat,aK] = ...
DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
kalman_smoother(ST,Z,R1,Q,H,Pinf,Pstar,data1,nobs,np,smpl);
if all(alphahat(:)==0)
kalman_algo = 2;
if ~estim_params_.ncn
......@@ -255,13 +259,14 @@ if any(any(H ~= 0)) % should be replaced by a flag
end
end
else
H = 0;
if kalman_algo == 1
if missing_value
[alphahat,etahat,ahat,aK] = missing_DiffuseKalmanSmoother1(T,R,Q, ...
Pinf,Pstar,Y,trend,nobs,np,smpl,mf,data_index);
else
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q, ...
Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = ...
kalman_smoother(T,R,Q,H,Pstar,data1,start,mf,kalman_tol,riccati_tol);
end
if all(alphahat(:)==0)
kalman_algo = 2;
......
matlab/dynare_config.m
View file @ 7668bc4a
......@@ -40,6 +40,7 @@ dynareroot = strrep(which('dynare'),'dynare.m','');
addpath([dynareroot '/distributions/'])
addpath([dynareroot '/kalman/'])
addpath([dynareroot '/kalman/likelihood'])
addpath([dynareroot '/kalman/smoother'])
addpath([dynareroot '/AIM/'])
addpath([dynareroot '/partial_information/'])
......
matlab/dynare_estimation_1.m
View file @ 7668bc4a
......@@ -350,32 +350,40 @@ missing_value = ~(number_of_observations == gend*n_varobs);
initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
if options_.mode_compute == 0
if options_.mode_compute == 0 && length(options_.mode_file) == 0
if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
oo_.Smoother.variance = P;
if options_.filter_covariance
oo_.Smoother.variance = P;
end
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead,i_endo,:);
oo_.FilteredVariablesKStepAheadVariances = PK(options_.filter_step_ahead,i_endo,i_endo,:);
oo_.FilteredVariablesShockDecomposition = decomp(options_.filter_step_ahead,i_endo,:,:);
oo_.FilteredVariablesKStepAhead = ...
aK(options_.filter_step_ahead,i_endo,:);
if ~isempty(PK)
oo_.FilteredVariablesKStepAheadVariances = ...
PK(options_.filter_step_ahead,i_endo,i_endo,:);
end
if ~isempty(decomp)
oo_.FilteredVariablesShockDecomposition = ...
decomp(options_.filter_step_ahead,i_endo,:,:);
end
end
for i=bayestopt_.smoother_saved_var_list
i1 = bayestop_.smoother_var_list(i);
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(dr.order_var(i1),:)) ' = atT(i,:)'';']);
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i1),:)) ' = squeeze(aK(1,i,:));']);
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(dr.order_var(i1),:)) ' = updated_variables(i,:)'';']);
for i=bayestopt_.smoother_saved_var_list'
i1 = dr.order_var(bayestopt_.smoother_var_list(i));
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = atT(i,:)'';']);
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ' = squeeze(aK(1,i,:));']);
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ' = updated_variables(i,:)'';']);
end
for i=1:M_.exo_nbr
eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
end
end
if length(options_.mode_file) == 0
return;
end
return;
end
%% Estimation of the posterior mode or likelihood mode
......@@ -1105,22 +1113,26 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.integration_order = d;
oo_.Smoother.variance = P;
i_endo_nbr = 1:M_.endo_nbr;
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
i_endo_nbr,:);
i_endo,:);
if isfield(options_,'kalman_algo')
if options_.kalman_algo > 2
oo_.FilteredVariablesKStepAheadVariances = PK(options_.filter_step_ahead,i_endo_nbr,i_endo_nbr,:);
if ~isempty(PK)
oo_.FilteredVariablesKStepAheadVariances = ...
PK(options_.filter_step_ahead,i_endo,i_endo_nbr,:);
end
if ~isempty(decomp)
oo_.FilteredVariablesShockDecomposition = ...
decomp(options_.filter_step_ahead,i_endo_nbr,:,:);
decomp(options_.filter_step_ahead,i_endo,:,:);
end
end
end
for i=1:M_.endo_nbr
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = atT(i,:)'';']);
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ' = squeeze(aK(1,i,:));']);
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(dr.order_var(i),:)) ...
for i=bayestopt_.smoother_saved_var_list'
i1 = dr.order_var(bayestopt_.smoother_var_list(i));
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = atT(i,:)'';']);
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ' = squeeze(aK(1,i,:));']);
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ...
' = updated_variables(i,:)'';']);
end
[nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
......
matlab/kalman/smoother/kalman_smoother.m 0 → 100644
View file @ 7668bc4a
function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = kalman_smoother(T,R,Q,H,P0,Y,start,mf,kalman_tol,riccati_tol)
% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = kalman_smoother(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
% Computes the kalman smoother of a stationary state space model.
%
% 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.
% P0 [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
% alphahat: smoothed state variables (a_{t|T})
% etahat: smoothed shocks
% atilde: matrix of updated variables (a_{t|t})
% aK: 3D array of k step ahead filtered state variables (a_{t+k|t})
% SPECIAL REQUIREMENTS
% 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
%
% 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/>.
global options_
option_filter_covariance = options_.filter_covariance;
option_filter_decomposition = options_.filter_decomposition;
nk = options_.nk;
smpl = size(Y,2); % Sample size.
mm = size(T,2); % Number of state variables.
pp = size(Y,1); % Maximum number of
% observed variables.
rr = size(Q,1);
v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
atilde = zeros(mm,smpl);
K = zeros(mm,pp,smpl);
aK = zeros(nk,mm,smpl+nk);
iF = zeros(pp,pp,smpl);
P = zeros(mm,mm,smpl+1);
QQ = R*Q*R';
QRt = Q*R';
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
oldK = 0;
if option_filter_covariance
PK = zeros(nk,mm,mm,smpl+nk);
else
PK = [];
end
if option_filter_decomposition
decomp = zeros(nk,mm,rr,smpl+nk);
else
decomp = [];
end
P(:,:,1) = P0;
t = 0;
notsteady = 1;
F_singular = 1;
while notsteady & t<smpl
t = t+1;
v(:,t) = Y(:,t)-a(mf,t);
F = P(mf,mf,t) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
atilde(:,t) = a(:,t);
a(:,t+1) = T*a(:,t);
P(:,:,t+1) = T*P(:,:,t)*T'+QQ;
end
else
F_singular = 0;
iF(:,:,t) = inv(F);
K1 = P(:,mf,t)*iF(:,:,t);
atilde(:,t) = a(:,t) + K1*v(:,t);
K(:,:,t) = T*K1;
a(:,t+1) = T*atilde(:,t);
P(:,:,t+1) = (T*P(:,:,t)-K(:,:,t)*P(mf,:,t))*T'+QQ;
end
aK(1,:,t+1) = a(:,t+1);
if option_filter_covariance
Pf = P(:,:,t);
Pf = T*Pf*T' + QQ;
PK(1,:,:,t+1) = Pf;
end
for jnk=2:nk,
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
if option_filter_covariance
Pf = T*Pf*T' + QQ;
PK(jnk,:,:,t+jnk) = Pf;
end
end
notsteady = max(max(abs(K(:,:,t)-oldK))) > riccati_tol;
oldK = K(:,:,t);
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;
while t < smpl
t = t+1;
v(:,t) = Y(:,t)-a(mf,t);
atilde(:,t) = a(:,t) + K1*v(:,t);
a(:,t+1) = T*atilde(:,t);
aK(1,:,t+1) = a(:,t+1);
if option_filter_covariance
Pf = P(:,:,t);
Pf = T*Pf*T' + QQ;
PK(1,:,:,t+1) = Pf;
end
for jnk=2:nk,
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
if option_filter_covariance
Pf = T*Pf*T' + QQ;
PK(jnk,:,:,t+jnk) = Pf;
end
end
end
K= cat(3,K(:,:,1:t),repmat(K(:,:,t0),[1 1 smpl-t0+1]));
P = cat(3,P(:,:,1:t),repmat(P(:,:,t0),[1 1 smpl-t0+1]));
iF = cat(3,iF(:,:,1:t),repmat(iF(:,:,t0),[1 1 smpl-t0+1]));
end
t = smpl+1;
while t>1
t = t-1;
r(:,t) = T'*r(:,t+1);
r(mf,t) = r(mf,t)+iF(:,:,t)*v(:,t) - K(:,:,t)'*r(:,t+1);
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
end
epsilonhat = Y-alphahat(mf,:);
if option_filter_decomposition
ZRQinv = inv(QQ(mf,mf));
for t = 1:smpl
% calculate eta_tm1t
eta = QRt(:,mf)*iF(:,:,t)*v(:,t);
AAA = P(:,mf,t)*ZRQinv*bsxfun(@times,R(mf,:),eta');
% calculate decomposition
Ttok = eye(mm,mm);
decomp(1,:,:,t+1) = AAA;
for h = 2:nk
AAA = T*AAA;
decomp(h,:,:,t+h) = AAA;
end
end
end
if ~option_filter_covariance
P = [];
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment