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

Merge pull request #1324 from JohannesPfeifer/kalman 

Three bugfixes for missing_DiffuseKalmanSmootherH1_Z.m that led to wrong results
time-shift
MichelJuillard 2016-11-13 19:13:10 +01:00 committed by GitHub
parent 920a9911ee 3ebc982b02
commit 31d2c7e332
28 changed files with 404 additions and 140 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
doc/dynare.texi
View File

@ -5919,7 +5919,7 @@ Use the Univariate Diffuse Kalman Filter
Default value is @code{0}. In case of missing observations of single or all series, Dynare treats those missing values as unobserved states and uses the Kalman filter to infer their value (see e.g. @cite{Durbin and Koopman (2012), Ch. 4.10})
This procedure has the advantage of being capable of dealing with observations where the forecast error variance matrix becomes singular for some variable(s).
If this happens, the respective observation enters with a weight of zero in the log-likelihood, i.e. this observation for the respective variable(s) is dropped
from the likelihood computations (for details see @cite{Durbin and Koopman (2012), Ch. 6.4 and 7.2.5}). If the use of a multivariate Kalman filter is specified and a
from the likelihood computations (for details see @cite{Durbin and Koopman (2012), Ch. 6.4 and 7.2.5} and @cite{Koopman and Durbin (2000)}). If the use of a multivariate Kalman filter is specified and a
singularity is encountered, Dynare by default automatically switches to the univariate Kalman filter for this parameter draw. This behavior can be changed via the
@ref{use_univariate_filters_if_singularity_is_detected} option.
@ -5954,12 +5954,16 @@ Triggers the computation k-step ahead filtered values, i.e. @math{E_{t}{y_{t+k}}
@anchor{filter_decomposition} Triggers the computation of the shock
decomposition of the above k-step ahead filtered values. Stores results in @code{oo_.FilteredVariablesShockDecomposition}.
@item smoothed_state_uncertainty
@anchor{smoothed_state_uncertainty} Triggers the computation of the variance of smoothed estimates, i.e.
@code{Var_T(y_t)}. Stores results in @code{oo_.Smoother.State_uncertainty}.
@item diffuse_filter
@anchor{diffuse_filter}
Uses the diffuse Kalman filter (as described in
@cite{Durbin and Koopman (2012)} and @cite{Koopman and Durbin
(2003)}) to estimate models with non-stationary observed variables.
(2003)} for the multivariate and @cite{Koopman and Durbin
(2000)} for the univariate filter) to estimate models with non-stationary observed variables.
When @code{diffuse_filter} is used the @code{lik_init} option of
@code{estimation} has no effect.
@ -6512,6 +6516,19 @@ two dimensions. The third dimension of the array provides the
observation for which the forecast has been made.
@end defvr
@defvr {MATLAB/Octave variable} oo_.Smoother.State_uncertainty
@anchor{oo_.Smoother.State_uncertainty}
Three-dimensional array set by the @code{estimation} command (if used with the
@code{smoother}) without Metropolis,
or by the @code{calib_smoother} command, if the @code{o_smoothed_state_uncertainty} option
has been requested.
Contains the series of covariance matrices for the state estimate given the full data
from the Kalman smoother. The @code{M_.endo_nbr} times @code{M_.endo_nbr} times
@code{T} array contains the variables in declaration order along the first
two dimensions. The third dimension of the array provides the
observation for which the smoothed estimate has been made.
@end defvr
@defvr {MATLAB/Octave variable} oo_.Smoother.SteadyState
@anchor{oo_.Smoother.SteadyState}
Variable set by the @code{estimation} command (if used with the
@ -14128,6 +14145,11 @@ Dynamics and Control}, 32(11), 3397--3414
@item
Koop, Gary (2003), @i{Bayesian Econometrics}, John Wiley & Sons
@item
Koopman, S. J. and J. Durbin (2000): ``Fast Filtering and Smoothing for
Multivariate State Space Models,'' @i{Journal of Time
Series Analysis}, 21(3), 281--296
@item
Koopman, S. J. and J. Durbin (2003): ``Filtering and Smoothing of
State Vector for Diffuse State Space Models,'' @i{Journal of Time

15
matlab/DsgeSmoother.m
View File

@ -1,4 +1,4 @@
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition,state_uncertainty] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
% Estimation of the smoothed variables and innovations.
%
% INPUTS
@ -25,6 +25,8 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% o decomp (K*m*r*(T+K)) 4D array of shock decomposition of k-step ahead
% filtered variables (decision-rule order)
% o trend_addition [double] (n*T) pure trend component; stored in options_.varobs order
% o state_uncertainty [double] (K,K,T) array, storing the uncertainty
% about the smoothed state (decision-rule order)
%
% Notes:
% m: number of endogenous variables (M_.endo_nbr)
@ -233,10 +235,10 @@ ST = T;
R1 = R;
if kalman_algo == 1 || kalman_algo == 3
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH1_Z(ST, ...
Z,R1,Q,H,Pinf,Pstar, ...
data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition);
options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition,options_.smoothed_state_uncertainty);
if isinf(alphahat)
if kalman_algo == 1
kalman_algo = 2;
@ -273,11 +275,11 @@ if kalman_algo == 2 || kalman_algo == 4
end
end
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(ST, ...
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty] = missing_DiffuseKalmanSmootherH3_Z(ST, ...
Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ...
options_.filter_decomposition);
options_.filter_decomposition,options_.smoothed_state_uncertainty);
end
@ -299,4 +301,7 @@ if expanded_state_vector_for_univariate_filter && (kalman_algo == 2 || kalman_al
if ~isempty(P)
P = P(k,k,:);
end
if ~isempty(state_uncertainty)
state_uncertainty = state_uncertainty(k,k,:);
end
end

8
matlab/dynare_estimation_1.m
View File

@ -163,8 +163,8 @@ end
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value);
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value);
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
end
%reset qz_criterium
options_.qz_criterium=qz_criterium_old;
@ -483,8 +483,8 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && options_.partial_information == 0 % to be fixed
%% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
[oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
[oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
if ~options_.nograph,
[nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);

4
matlab/evaluate_smoother.m
View File

@ -97,9 +97,9 @@ if ischar(parameters)
end
end
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = ...
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty] = ...
DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
if nargout==2
Smoothed_variables_declaration_order_deviation_form=atT(oo_.dr.inv_order_var(bayestopt_.smoother_var_list),:);

1
matlab/global_initialization.m
View File

@ -413,6 +413,7 @@ options_.bayesian_th_moments = 0;
options_.diffuse_filter = 0;
options_.filter_step_ahead = [];
options_.filtered_vars = 0;
options_.smoothed_state_uncertainty = 0;
options_.first_obs = NaN;
options_.nobs = NaN;
options_.kalman_algo = 0;

2
matlab/imcforecast.m
View File

@ -123,7 +123,9 @@ if estimated_model
%store qz_criterium
qz_criterium_old=options_.qz_criterium;
options_=select_qz_criterium_value(options_);
options_state_uncertainty_old=options_.state_uncertainty;
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam,gend,data,data_index,missing_value);
options_.state_uncertainty=options_state_uncertainty_old;
%get constant part
if options_.noconstant
constant = zeros(size(ys,1),options_cond_fcst.periods+1);

10
matlab/kalman/likelihood/kalman_filter_d.m
View File

@ -87,7 +87,7 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
else
iFstar = inv(Fstar);
dFstar = det(Fstar);
Kstar = Pstar*Z'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1}
Kstar = Pstar*Z'*iFstar; %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0)
dlik(s)= log(dFstar) + v'*iFstar*v; %set w_t to bottom case in bottom equation page 172, DK (2012)
Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
Pstar = T*(Pstar-Pstar*Z'*Kstar')*T'+QQ; % (5.17) DK (2012)
@ -96,8 +96,8 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
end
else %F_{\infty,t} positive definite
%To compare to DK (2012), this block makes use of the following transformation
%Kstar=K^{(1)}*T^{-1}=M_{*}*F^{(1)}+M_{\infty}*F^{(2)}
% =P_{*}*Z'*F^{(1)}+P_{\infty}*Z'*((-1)*(F_{\infty}^{-1})*F_{*}*(F_{\infty}^{-1}))
%Kstar=T^{-1}*K^{(1)}=M_{*}*F^{(1)}+M_{\infty}*F^{(2)}
% =P_{*}*Z'*F^{(1)}+P_{\infty}*Z'*((-1)*(-F_{\infty}^{-1})*F_{*}*(F_{\infty}^{-1}))
% =[P_{*}*Z'-Kinf*F_{*})]*F^{(1)}
%Make use of L^{0}'=(T-K^{(0)}*Z)'=(T-T*M_{\infty}*F^{(1)}*Z)'
% =(T-T*P_{\infty*Z'*F^{(1)}*Z)'=(T-T*Kinf*Z)'
@ -106,9 +106,9 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
% =T*[(P_{\infty}*(-K^{(1)*Z}))+P_{*}*(I-Z'*Kinf')*T'+RQR]
dlik(s)= log(det(Finf)); %set w_t to top case in bottom equation page 172, DK (2012)
iFinf = inv(Finf);
Kinf = Pinf*Z'*iFinf; %define Kinf=K_0*T^{-1} with M_{\infty}=Pinf*Z'
Kinf = Pinf*Z'*iFinf; %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z'
Fstar = Z*Pstar*Z' + H; %(5.7) DK(2012)
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; %(5.12) DK(2012)
Kstar = (Pstar*Z'-Kinf*Fstar)*iFinf; %(5.12) DK(2012); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix
Pstar = T*(Pstar-Pstar*Z'*Kinf'-Pinf*Z'*Kstar')*T'+QQ; %(5.14) DK(2012)
Pinf = T*(Pinf-Pinf*Z'*Kinf')*T'; %(5.14) DK(2012)
a = T*(a+Kinf*v); %(5.13) DK(2012)

6
matlab/kalman/likelihood/missing_observations_kalman_filter_d.m
View File

@ -106,10 +106,10 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
else
iFstar = inv(Fstar);
dFstar = det(Fstar);
Kstar = Pstar*ZZ'*iFstar; %(5.15) of DK (2012) with Kstar=K^(0)*T^{-1}
Kstar = Pstar*ZZ'*iFstar; %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0)
dlik(s) = log(dFstar) + v'*iFstar*v + length(d_index)*log(2*pi); %set w_t to bottom case in bottom equation page 172, DK (2012)
Pinf = T*Pinf*transpose(T); % (5.16) DK (2012)
Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ; % (5.17) DK (2012)
Pstar = T*(Pstar-Pstar*ZZ'*Kstar')*T'+QQ; % (5.17) DK (2012) with L_0 plugged in
a = T*(a+Kstar*v); % (5.13) DK (2012)
end
end
@ -119,7 +119,7 @@ while rank(Pinf,diffuse_kalman_tol) && (t<=last)
Kinf = Pinf*ZZ'*iFinf;
%see notes in kalman_filter_d.m for details of computations
Fstar = ZZ*Pstar*ZZ' + H(d_index,d_index); %(5.7) DK(2012)
Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf; %(5.12) DK(2012)
Kstar = (Pstar*ZZ'-Kinf*Fstar)*iFinf; %(5.12) DK(2012); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix
Pstar = T*(Pstar-Pstar*ZZ'*Kinf'-Pinf*ZZ'*Kstar')*T'+QQ; %(5.14) DK(2012)
Pinf = T*(Pinf-Pinf*ZZ'*Kinf')*T'; %(5.14) DK(2012)
a = T*(a+Kinf*v); %(5.13) DK(2012)

162
matlab/missing_DiffuseKalmanSmootherH1_Z.m
View File

@ -1,6 +1,6 @@
function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag)
function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag)
% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
% Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix.
%
% INPUTS
@ -20,6 +20,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal
% kalman_tol tolerance for reciprocal condition number
% diffuse_kalman_tol tolerance for reciprocal condition number (for Finf) and the rank of Pinf
% decomp_flag if true, compute filter decomposition
% state_uncertainty_flag if true, compute uncertainty about smoothed
% state estimate
%
% OUTPUTS
% alphahat: smoothed variables (a_{t|T})
@ -33,6 +35,7 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal
% PK: 4D array of k-step ahead forecast error variance
% matrices (meaningless for periods 1:d)
% decomp: decomposition of the effect of shocks on filtered values
% V: 3D array of state uncertainty matrices
%
% Notes:
% Outputs are stored in decision-rule order, i.e. to get variables in order of declaration
@ -43,6 +46,8 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal
% 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).
% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
% Second Edition, Ch. 5
% Copyright (C) 2004-2016 Dynare Team
%
@ -78,76 +83,92 @@ aK = zeros(nk,mm,smpl+nk);
PK = zeros(nk,mm,mm,smpl+nk);
iF = zeros(pp,pp,smpl);
Fstar = zeros(pp,pp,smpl);
iFstar = zeros(pp,pp,smpl);
iFinf = zeros(pp,pp,smpl);
K = zeros(mm,pp,smpl);
L = zeros(mm,mm,smpl);
Linf = zeros(mm,mm,smpl);
Lstar = zeros(mm,mm,smpl);
Kstar = zeros(mm,pp,smpl);
Kinf = zeros(mm,pp,smpl);
P = zeros(mm,mm,smpl+1);
Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
steady = smpl;
Pstar = zeros(spstar(1),spstar(2),smpl+1);
Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl+1);
Pinf(:,:,1) = Pinf1;
rr = size(Q,1);
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
epsilonhat = zeros(rr,smpl);
epsilonhat = zeros(rr,smpl);
r = zeros(mm,smpl+1);
Finf_singular = zeros(1,smpl);
if state_uncertainty_flag
V = zeros(mm,mm,smpl);
N = zeros(mm,mm,smpl+1);
else
V=[];
end
t = 0;
while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
t = t+1;
di = data_index{t};
if isempty(di)
atilde(:,t) = a(:,t);
%no observations, propagate forward without updating based on
%observations
atilde(:,t) = a(:,t);
a(:,t+1) = T*atilde(:,t);
Linf(:,:,t) = T;
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T' + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
else
ZZ = Z(di,:);
v(di,t)= Y(di,t) - ZZ*a(:,t);
Finf = ZZ*Pinf(:,:,t)*ZZ';
if rcond(Finf) < diffuse_kalman_tol
if ~all(abs(Finf(:)) < diffuse_kalman_tol)
ZZ = Z(di,:); %span selector matrix
v(di,t)= Y(di,t) - ZZ*a(:,t); %get prediction error v^(0) in (5.13) DK (2012)
Finf = ZZ*Pinf(:,:,t)*ZZ'; % (5.7) in DK (2012)
if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0
% The univariate diffuse kalman filter should be used.
alphahat = Inf;
return
else
Fstar(:,:,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);
if rcond(Fstar(:,:,t)) < kalman_tol
else %rank of F_{\infty,t} is 0
Finf_singular(1,t) = 1;
Fstar(:,:,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); % (5.7) in DK (2012)
if rcond(Fstar(:,:,t)) < kalman_tol %F_{*} is singular
if ~all(abs(Fstar(:,:,t))<kalman_tol)
% The univariate diffuse kalman filter should be used.
alphahat = Inf;
return
else
else %rank 0
a(:,t+1) = T*a(:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
end
else
iFstar = inv(Fstar(:,:,t));
Kstar(:,:,t) = Pstar(:,:,t)*ZZ'*iFstar(:,:,t);
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T);
Pstar(:,:,t+1) = T*(Pstar(:,:,t)-Pstar(:,:,t)*ZZ'*Kstar(:,:,t)')*T'+QQ;
a(:,t+1) = T*(a(:,t)+Kstar(:,:,t)*v(:,t));
iFstar(:,:,t) = inv(Fstar(:,:,t));
Kstar(:,:,t) = Pstar(:,:,t)*ZZ'*iFstar(:,:,t); %(5.15) of DK (2012) with Kstar=T^{-1}*K^(0)
Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); % DK (2012), 5.16
Lstar(:,:,t) = T - T*Kstar(:,di,t)*ZZ; %L^(0) in DK (2012), eq. 5.12
Pstar(:,:,t+1) = T*Pstar(:,:,t)*Lstar(:,:,t)'+QQ; % (5.17) DK (2012)
a(:,t+1) = T*(a(:,t)+Kstar(:,:,t)*v(:,t)); % (5.13) DK (2012)
end
end
else
%see notes in kalman_filter_d.m for details of computations
iFinf(di,di,t) = inv(Finf);
PZI = Pinf(:,:,t)*ZZ'*iFinf(di,di,t);
atilde(:,t) = a(:,t) + PZI*v(di,t);
Kinf(:,di,t) = T*PZI;
Linf(:,:,t) = T - Kinf(:,di,t)*ZZ;
Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di);
Kstar(:,di,t) = (T*Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*ZZ'*Kinf(:,di,t)'-T*Pinf(:,:,t)*ZZ'*Kstar(:,di,t)' + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*ZZ'*Kinf(:,di,t)';
Kinf(:,di,t) = Pinf(:,:,t)*ZZ'*iFinf(di,di,t); %define Kinf=T^{-1}*K_0 with M_{\infty}=Pinf*Z'
atilde(:,t) = a(:,t) + Kinf(:,di,t)*v(di,t);
Linf(:,:,t) = T - T*Kinf(:,di,t)*ZZ; %L^(0) in DK (2012), eq. 5.12
Fstar(di,di,t) = ZZ*Pstar(:,:,t)*ZZ' + H(di,di); %(5.7) DK(2012)
Kstar(:,di,t) = (Pstar(:,:,t)*ZZ'-Kinf(:,di,t)*Fstar(di,di,t))*iFinf(di,di,t); %(5.12) DK(2012) with Kstar=T^{-1}*K^(1); note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf", but it is correct in their appendix
Pstar(:,:,t+1) = T*Pstar(:,:,t)*Linf(:,:,t)'-T*Kinf(:,di,t)*Finf*Kstar(:,di,t)'*T' + QQ; %(5.14) DK(2012)
Pinf(:,:,t+1) = T*Pinf(:,:,t)*Linf(:,:,t)'; %(5.14) DK(2012)
end
a(:,t+1) = T*atilde(:,t);
a(:,t+1) = T*atilde(:,t);
aK(1,:,t+1) = a(:,t+1);
% isn't a meaningless as long as we are in the diffuse part? MJ
for jnk=2:nk,
for jnk=2:nk
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
end
@ -155,20 +176,21 @@ end
d = t;
P(:,:,d+1) = Pstar(:,:,d+1);
iFinf = iFinf(:,:,1:d);
iFstar= iFstar(:,:,1:d);
Linf = Linf(:,:,1:d);
Fstar = Fstar(:,:,1:d);
Lstar = Lstar(:,:,1:d);
Kstar = Kstar(:,:,1:d);
Pstar = Pstar(:,:,1:d);
Pinf = Pinf(:,:,1:d);
notsteady = 1;
while notsteady && t<smpl
t = t+1;
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); % make sure P is symmetric
di = data_index{t};
if isempty(di)
atilde(:,t) = a(:,t);
atilde(:,t) = a(:,t);
L(:,:,t) = T;
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %p. 111, DK(2012)
else
ZZ = Z(di,:);
v(di,t) = Y(di,t) - ZZ*a(:,t);
@ -220,34 +242,79 @@ end
% $$$ PK(jnk,:,:,t+jnk) = Pf;
% $$$ end
% $$$ end
%% backward pass; r_T and N_T, stored in entry (smpl+1) were initialized at 0
t = smpl+1;
while t>d+1
t = t-1;
di = data_index{t};
if isempty(di)
r(:,t) = L(:,:,t)'*r(:,t+1);
% in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
r(:,t) = L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0
if state_uncertainty_flag
N(:,:,t)=L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0
end
else
ZZ = Z(di,:);
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1);
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38
if state_uncertainty_flag
N(:,:,t)=ZZ'*iF(di,di,t)*ZZ+L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42
end
end
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); %DK (2012), eq. 4.35
etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63
if state_uncertainty_flag
V(:,:,t) = P(:,:,t)-P(:,:,t)*N(:,:,t)*P(:,:,t); %DK (2012), eq. 4.43
end
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
end
if d
if d %diffuse periods
% initialize r_d^(0) and r_d^(1) as below DK (2012), eq. 5.23
r0 = zeros(mm,d+1);
r0(:,d+1) = r(:,d+1);
r1 = zeros(mm,d+1);
r0(:,d+1) = r(:,d+1); %set r0_{d}, i.e. shifted by one period
r1 = zeros(mm,d+1); %set r1_{d}, i.e. shifted by one period
if state_uncertainty_flag
%N_0 at (d+1) is N(d+1), so we can use N for continuing and storing N_0-recursion
N_1=zeros(mm,mm,d+1); %set N_1_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26
N_2=zeros(mm,mm,d+1); %set N_2_{d}=0, i.e. shifted by one period, below DK (2012), eq. 5.26
end
for t = d:-1:1
r0(:,t) = Linf(:,:,t)'*r0(:,t+1);
di = data_index{t};
if isempty(di)
r1(:,t) = Linf(:,:,t)'*r1(:,t+1);
else
r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*r0(:,t+1)) ...
+ Linf(:,:,t)'*r1(:,t+1);
if ~Finf_singular(1,t)
r0(:,t) = Linf(:,:,t)'*r0(:,t+1); % DK (2012), eq. 5.21 where L^(0) is named Linf
r1(:,t) = Z(di,:)'*(iFinf(di,di,t)*v(di,t)-Kstar(:,di,t)'*T'*r0(:,t+1)) ...
+ Linf(:,:,t)'*r1(:,t+1); % DK (2012), eq. 5.21, noting that i) F^(1)=(F^Inf)^(-1)(see 5.10), ii) where L^(0) is named Linf, and iii) Kstar=T^{-1}*K^(1)
if state_uncertainty_flag
L_1=(-T*Kstar(:,di,t)*Z(di,:)); % noting that Kstar=T^{-1}*K^(1)
N(:,:,t)=Linf(:,:,t)'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.19, noting that L^(0) is named Linf
N_1(:,:,t)=Z(di,:)'*iFinf(di,di,t)*Z(di,:)+Linf(:,:,t)'*N_1(:,:,t+1)*Linf(:,:,t)...
+L_1'*N(:,:,t+1)*Linf(:,:,t); % DK (2012), eq. 5.29; note that, compared to DK (2003) this drops the term (L_1'*N(:,:,t+1)*Linf(:,:,t))' in the recursion due to it entering premultiplied by Pinf when computing V, and Pinf*Linf'*N=0
N_2(:,:,t)=Z(di,:)'*(-iFinf(di,di,t)*Fstar(di,di,t)*iFinf(di,di,t))*Z(di,:) ...
+ Linf(:,:,t)'*N_2(:,:,t+1)*Linf(:,:,t)...
+ Linf(:,:,t)'*N_1(:,:,t+1)*L_1...
+ L_1'*N_1(:,:,t+1)'*Linf(:,:,t)...
+ L_1'*N(:,:,t+1)*L_1; % DK (2012), eq. 5.29
end
else
r0(:,t) = Z(di,:)'*iFstar(di,di,t)*v(di,t)-Lstar(:,di,t)'*r0(:,t+1); % DK (2003), eq. (14)
r1(:,t) = T'*r1(:,t+1); % DK (2003), eq. (14)
if state_uncertainty_flag
N(:,:,t)=Z(di,:)'*iFstar(di,di,t)*Z(di,:)...
+Lstar(:,:,t)'*N(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14)
N_1(:,:,t)=T'*N_1(:,:,t+1)*Lstar(:,:,t); % DK (2003), eq. (14)
N_2(:,:,t)=T'*N_2(:,:,t+1)*T'; % DK (2003), eq. (14)
end
end
end
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23
etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135
if state_uncertainty_flag
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
- Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)...
- Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30
end
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t);
etahat(:,t) = QRt*r0(:,t);
end
end
@ -260,7 +327,6 @@ if decomp_flag
eta_tm1t = QRt*Z(di,:)'*iF(di,di,t)*v(di,t);
AAA = P(:,:,t)*Z(di,:)'*ZRQinv(di,di)*bsxfun(@times,Z(di,:)*R,eta_tm1t');
% calculate decomposition
Ttok = eye(mm,mm);
decomp(1,:,:,t+1) = AAA;
for h = 2:nk
AAA = T*AAA;

151
matlab/missing_DiffuseKalmanSmootherH3_Z.m
View File

@ -1,5 +1,5 @@
function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag)
% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag)
function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
% Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
% Univariate treatment of multivariate time series.
%
@ -20,6 +20,8 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm
% kalman_tol tolerance for zero divider
% diffuse_kalman_tol tolerance for zero divider
% decomp_flag if true, compute filter decomposition
% state_uncertainty_flag if true, compute uncertainty about smoothed
% state estimate
%
% OUTPUTS
% alphahat: smoothed state variables (a_{t|T})
@ -33,12 +35,22 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp] = missing_DiffuseKalmanSm
% PK: 4D array of k-step ahead forecast error variance
% matrices (meaningless for periods 1:d)
% decomp: decomposition of the effect of shocks on filtered values
% V: 3D array of state uncertainty matrices
%
% Notes:
% Outputs are stored in decision-rule order, i.e. to get variables in order of declaration
% as in M_.endo_names, ones needs code along the lines of:
% variables_declaration_order(dr.order_var,:) = alphahat
%
% Algorithm:
%
% Uses the univariate filter as described in Durbin/Koopman (2012): "Time
% Series Analysis by State Space Methods", Oxford University Press,
% Second Edition, Ch. 6.4 + 7.2.5
% and
% Koopman/Durbin (2000): "Fast Filtering and Smoothing for Multivariatze State Space
% Models", in Journal of Time Series Analysis, vol. 21(3), pp. 281-296.
%
% 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
@ -84,21 +96,29 @@ Finf = zeros(pp,smpl);
Fi = zeros(pp,smpl);
Ki = zeros(mm,pp,smpl);
Kstar = zeros(mm,pp,smpl);
Kinf = zeros(spstar(1),pp,smpl);
P = zeros(mm,mm,smpl+1);
P1 = P;
PK = zeros(nk,mm,mm,smpl+nk);
Pstar = zeros(spstar(1),spstar(2),smpl); Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl); Pinf(:,:,1) = Pinf1;
Pstar = zeros(spstar(1),spstar(2),smpl);
Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl);
Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
steady = smpl;
rr = size(Q,1); % number of structural shocks
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
etahat = zeros(rr,smpl);
epsilonhat = zeros(rr,smpl);
r = zeros(mm,smpl);
if state_uncertainty_flag
V = zeros(mm,mm,smpl);
N = zeros(mm,mm,smpl);
else
V=[];
end
t = 0;
icc=0;
@ -111,23 +131,28 @@ while newRank && t < smpl
di = data_index{t}';
for i=di
Zi = Z(i,:);
v(i,t) = Y(i,t)-Zi*a(:,t);
Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i);
Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
Kstar(:,i,t) = Pstar(:,:,t)*Zi';
if Finf(i,t) > diffuse_kalman_tol && newRank
v(i,t) = Y(i,t)-Zi*a(:,t); % nu_{t,i} in 6.13 in DK (2012)
Fstar(i,t) = Zi*Pstar(:,:,t)*Zi' +H(i); % F_{*,t} in 5.7 in DK (2012), relies on H being diagonal
Finf(i,t) = Zi*Pinf(:,:,t)*Zi'; % F_{\infty,t} in 5.7 in DK (2012)
Kstar(:,i,t) = Pstar(:,:,t)*Zi'; % KD (2000), eq. (15)
if Finf(i,t) > diffuse_kalman_tol && newRank % F_{\infty,t,i} = 0, use upper part of bracket on p. 175 DK (2012) for w_{t,i}
icc=icc+1;
Kinf(:,i,t) = Pinf(:,:,t)*Zi';
Kinf(:,i,t) = Pinf(:,:,t)*Zi'; % KD (2000), eq. (15)
Kinf_Finf = Kinf(:,i,t)/Finf(i,t);
a(:,t) = a(:,t) + Kinf_Finf*v(i,t);
a(:,t) = a(:,t) + Kinf_Finf*v(i,t); % KD (2000), eq. (16)
Pstar(:,:,t) = Pstar(:,:,t) + ...
Kinf(:,i,t)*Kinf_Finf'*(Fstar(i,t)/Finf(i,t)) - ...
Kstar(:,i,t)*Kinf_Finf' - ...
Kinf_Finf*Kstar(:,i,t)';
Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t);
Kinf_Finf*Kstar(:,i,t)'; % KD (2000), eq. (16)
Pinf(:,:,t) = Pinf(:,:,t) - Kinf(:,i,t)*Kinf(:,i,t)'/Finf(i,t); % KD (2000), eq. (16)
elseif Fstar(i,t) > kalman_tol
a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t);
Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t);
a(:,t) = a(:,t) + Kstar(:,i,t)*v(i,t)/Fstar(i,t); % KD (2000), eq. (17)
Pstar(:,:,t) = Pstar(:,:,t) - Kstar(:,i,t)*Kstar(:,i,t)'/Fstar(i,t); % KD (2000), eq. (17)
% Pinf is passed through unaltered, see eq. (17) of
% Koopman/Durbin (2000)
else
% do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
% p. 157, DK (2012)
end
end
if newRank
@ -142,8 +167,7 @@ while newRank && t < smpl
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
P0=Pinf(:,:,t+1);
if newRank,
if newRank
newRank = rank(Pinf(:,:,t+1),diffuse_kalman_tol);
end
if oldRank ~= newRank
@ -169,15 +193,18 @@ while notsteady && t<smpl
di = data_index{t}';
for i=di
Zi = Z(i,:);
v(i,t) = Y(i,t) - Zi*a(:,t);
Fi(i,t) = Zi*P(:,:,t)*Zi' + H(i);
Ki(:,i,t) = P(:,:,t)*Zi';
v(i,t) = Y(i,t) - Zi*a(:,t); % nu_{t,i} in 6.13 in DK (2012)
Fi(i,t) = Zi*P(:,:,t)*Zi' + H(i); % F_{t,i} in 6.13 in DK (2012), relies on H being diagonal
Ki(:,i,t) = P(:,:,t)*Zi'; % K_{t,i}*F_(i,t) in 6.13 in DK (2012)
if Fi(i,t) > kalman_tol
a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t);
a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t); %filtering according to (6.13) in DK (2012)
P(:,:,t) = P(:,:,t) - Ki(:,i,t)*Ki(:,i,t)'/Fi(i,t); %filtering according to (6.13) in DK (2012)
else
% do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
% p. 157, DK (2012)
end
end
a1(:,t+1) = T*a(:,t);
a1(:,t+1) = T*a(:,t); %transition according to (6.14) in DK (2012)
Pf = P(:,:,t);
aK(1,:,t+1) = a1(:,t+1);
for jnk=1:nk
@ -187,7 +214,7 @@ while notsteady && t<smpl
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
end
P(:,:,t+1) = T*P(:,:,t)*T' + QQ;
P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %transition according to (6.14) in DK (2012)
% notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
end
% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
@ -221,43 +248,93 @@ end
% $$$ PK(jnk,:,:,t+jnk) = Pf;
% $$$ end
% $$$ end
%% do backward pass
ri=zeros(mm,1);
if state_uncertainty_flag
Ni=zeros(mm,mm);
end
t = smpl+1;
while t > d+1
t = t-1;
di = flipud(data_index{t})';
for i = di
if Fi(i,t) > kalman_tol
ri = Z(i,:)'/Fi(i,t)*v(i,t)+ri-Ki(:,i,t)'*ri/Fi(i,t)*Z(i,:)';
Li = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
ri = Z(i,:)'/Fi(i,t)*v(i,t)+Li'*ri; % DK (2012), 6.15, equation for r_{t,i-1}
if state_uncertainty_flag
Ni = Z(i,:)'/Fi(i,t)*Z(i,:)+Li'*Ni*Li; % KD (2000), eq. (23)
end
end
end
r(:,t) = ri;
r(:,t) = ri; % DK (2012), below 6.15, r_{t-1}=r_{t,0}
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
etahat(:,t) = QRt*r(:,t);
ri = T'*ri;
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
if state_uncertainty_flag
N(:,:,t) = Ni; % DK (2012), below 6.15, N_{t-1}=N_{t,0}
V(:,:,t) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); % KD (2000), eq. (7) with N_{t-1} stored in N(:,:,t)
Ni = T'*Ni*T; % KD (2000), eq. (23), equation for N_{t-1,p_{t-1}}
end
end
if d
r0 = zeros(mm,d);
r0(:,d) = ri;
r1 = zeros(mm,d);
if state_uncertainty_flag
%N_0 at (d+1) is N(d+1), so we can use N for continuing and storing N_0-recursion
N_0=zeros(mm,mm,d); %set N_1_{d}=0, below KD (2000), eq. (24)
N_0(:,:,d) = Ni;
N_1=zeros(mm,mm,d); %set N_1_{d}=0, below KD (2000), eq. (24)
N_2=zeros(mm,mm,d); %set N_2_{d}=0, below KD (2000), eq. (24)
end
for t = d:-1:1
di = flipud(data_index{t})';
for i = di
if Finf(i,t) > diffuse_kalman_tol
if Finf(i,t) > diffuse_kalman_tol
% recursions need to be from highest to lowest term in order to not
% overwrite lower terms still needed in thisstep
r1(:,t) = Z(i,:)'*v(i,t)/Finf(i,t) + ...
(Kinf(:,i,t)'*Fstar(i,t)/Finf(i,t)-Kstar(:,i,t)')*r0(:,t)/Finf(i,t)*Z(i,:)' + ...
r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)';
r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)';
r1(:,t)-Kinf(:,i,t)'*r1(:,t)/Finf(i,t)*Z(i,:)'; % KD (2000), eq. (25) for r_1
r0(:,t) = r0(:,t)-Kinf(:,i,t)'*r0(:,t)/Finf(i,t)*Z(i,:)'; % KD (2000), eq. (25) for r_0
if state_uncertainty_flag
Linf = eye(mm) - Kinf(:,i,t)'/Finf(i,t);
L0 = (Kinf(:,i,t)*(Fstar(i,t)/Finf(i,t))-Kstar(:,i,t))*Z(i,:)/Finf(i,t);
N_2(:,:,t)=Z(i,:)'/Finf(i,t)^2*Z(i,:)*Fstar(i,t) ...
+ Linf'*N_2(:,:,t)*Linf...
+ Linf'*N_1(:,:,t)*L0...
+ L0'*N_1(:,:,t)'*Linf...
+ L0'*N_0(:,:,t)*L0; % DK (2012), eq. 5.29
N_1(:,:,t)=Z(i,:)'/Finf(i,t)*Z(i,:)+Linf'*N_1(:,:,t)*Linf...
+L0'*N_0(:,:,t)*Linf; % DK (2012), eq. 5.29; note that, compared to DK (2003) this drops the term (L_1'*N(:,:,t+1)*Linf(:,:,t))' in the recursion due to it entering premultiplied by Pinf when computing V, and Pinf*Linf'*N=0
N_0(:,:,t)=Linf'*N_0(:,:,t)*Linf; % DK (2012), eq. 5.19, noting that L^(0) is named Linf
end
elseif Fstar(i,t) > kalman_tol % step needed whe Finf == 0
r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+r0(:,t)-(Kstar(:,i,t)'*r0(:,t))/Fstar(i,t)*Z(i,:)';
L_i=eye(mm) - Kstar(:,i,t)*Z(i,:)*Fstar(i,t);
r0(:,t) = Z(i,:)'/Fstar(i,t)*v(i,t)+L_i'*r0(:,t); % propagate r0 and keep r1 fixed
if state_uncertainty_flag
N_0(:,:,t)=Z(i,:)'/Fstar(i,t)*Z(i,:)+L_i'*N_0(:,:,t)*L_i; % propagate N_0 and keep N_1 and N_2 fixed
end
end
end
alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t);
alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); % KD (2000), eq. (26)
r(:,t) = r0(:,t);
etahat(:,t) = QRt*r(:,t);
etahat(:,t) = QRt*r(:,t); % KD (2000), eq. (27)
if state_uncertainty_flag
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)...
-(Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t))'...
- Pinf(:,:,t)*N_1(:,:,t)*Pstar(:,:,t)...
- Pinf(:,:,t)*N_2(:,:,t)*Pinf(:,:,t); % DK (2012), eq. 5.30
end
if t > 1
r0(:,t-1) = T'*r0(:,t);
r1(:,t-1) = T'*r1(:,t);
r0(:,t-1) = T'*r0(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
r1(:,t-1) = T'*r1(:,t); % KD (2000), below eq. (25) r_{t-1,p_{t-1}}=T'*r_{t,0}
if state_uncertainty_flag
N_0(:,:,t-1)= T'*N_0(:,t)*T; % KD (2000), below eq. (25) N_{t-1,p_{t-1}}=T'*N_{t,0}*T
N_1(:,:,t-1)= T'*N_1(:,t)*T; % KD (2000), below eq. (25) N^1_{t-1,p_{t-1}}=T'*N^1_{t,0}*T
N_2(:,:,t-1)= T'*N_2(:,t)*T; % KD (2000), below eq. (25) N^2_{t-1,p_{t-1}}=T'*N^2_{t,0}*T
end
end
end
end

33
matlab/pm3.m
View File

@ -85,6 +85,7 @@ fprintf(['Estimation::mcmc: ' tit1 '\n']);
k = 0;
filter_step_ahead_indicator=0;
filter_covar_indicator=0;
state_uncert_indicator=0;
for file = 1:ifil
loaded_file=load([DirectoryName '/' M_.fname var_type int2str(file)]);
@ -118,6 +119,13 @@ for file = 1:ifil
end
k = k(end)+(1:size(stock,2));
stock1(:,k) = stock;
elseif strcmp(var_type,'_state_uncert')
if file==1 %on first run, initialize variable for storing filter_step_ahead
stock1_state_uncert=NaN(n1,n2,size(stock,3),B);
end
state_uncert_indicator=1;
k = k(end)+(1:size(stock,4));
stock1_state_uncert(:,:,:,k) = stock;
else
if file==1 %on first run, initialize variable for storing filter_step_ahead
stock1 = zeros(n1,n2,B);
@ -138,8 +146,7 @@ if filter_step_ahead_indicator
if options_.estimation.moments_posterior_density.indicator
Density_filter_step_ahead = zeros(options_.estimation.moments_posterior_density.gridpoints,2,filter_steps,nvar,n2);
end
end
if filter_covar_indicator
elseif filter_covar_indicator
draw_dimension=4;
oo_.FilterCovariance.Mean = squeeze(mean(stock1_filter_covar,draw_dimension));
oo_.FilterCovariance.Median = squeeze(median(stock1_filter_covar,draw_dimension));
@ -161,6 +168,28 @@ if filter_covar_indicator
oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2));
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);
return
elseif state_uncert_indicator
draw_dimension=4;
oo_.Smoother.State_uncertainty.Mean = squeeze(mean(stock1_state_uncert,draw_dimension));
oo_.Smoother.State_uncertainty.Median = squeeze(median(stock1_state_uncert,draw_dimension));
oo_.Smoother.State_uncertainty.var = squeeze(var(stock1_state_uncert,0,draw_dimension));
if size(stock1_state_uncert,draw_dimension)>2
hpd_interval = quantile(stock1_state_uncert,[(1-options_.mh_conf_sig)/2 (1-options_.mh_conf_sig)/2+options_.mh_conf_sig],draw_dimension);
else
size_matrix=size(stock1_state_uncert);
hpd_interval=NaN([size_matrix(1:3),2]);
end
if size(stock1_state_uncert,draw_dimension)>9
post_deciles =quantile(stock1_state_uncert,[0.1:0.1:0.9],draw_dimension);
else
size_matrix=size(stock1_state_uncert);
post_deciles=NaN([size_matrix(1:3),9]);
end
oo_.Smoother.State_uncertainty.post_deciles=post_deciles;
oo_.Smoother.State_uncertainty.HPDinf=squeeze(hpd_interval(:,:,:,1));
oo_.Smoother.State_uncertainty.HPDsup=squeeze(hpd_interval(:,:,:,2));
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);
return
end
if strcmp(var_type,'_trend_coeff') %two dimensional arrays

28
matlab/prior_posterior_statistics.m
View File

@ -117,6 +117,10 @@ if options_.filter_covariance
MAX_filter_covariance = min(B,ceil(MaxNumberOfBytes/(endo_nbr^2*(gend+1))/8));
end
if options_.smoothed_state_uncertainty
MAX_n_smoothed_state_uncertainty = min(B,ceil(MaxNumberOfBytes/((endo_nbr*endo_nbr)*gend)/8));
end
varlist = options_.varlist;
if isempty(varlist)
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
@ -129,13 +133,13 @@ for i=1:nvar
end
end
n_variables_to_fill=12;
n_variables_to_fill=13;
irun = ones(n_variables_to_fill,1);
ifil = zeros(n_variables_to_fill,1);
run_smoother = 0;
if options_.smoother || options_.forecast || options_.filter_step_ahead
if options_.smoother || options_.forecast || options_.filter_step_ahead || options_.smoothed_state_uncertainty
run_smoother = 1;
end
@ -144,11 +148,17 @@ if options_.filter_covariance
filter_covariance=1;
end
smoothed_state_uncertainty=0;
if options_.smoothed_state_uncertainty
smoothed_state_uncertainty=1;
end
% Store the variable mandatory for local/remote parallel computing.
localVars.type=type;
localVars.run_smoother=run_smoother;
localVars.filter_covariance=filter_covariance;
localVars.smoothed_state_uncertainty=smoothed_state_uncertainty;
localVars.gend=gend;
localVars.Y=Y;
localVars.data_index=data_index;
@ -181,6 +191,9 @@ end
if options_.filter_covariance
localVars.MAX_filter_covariance = MAX_filter_covariance;
end
if options_.smoothed_state_uncertainty
localVars.MAX_n_smoothed_state_uncertainty = MAX_n_smoothed_state_uncertainty ;
end
localVars.MAX_n_smoothed_constant=MAX_n_smoothed_constant;
localVars.MAX_n_smoothed_trend=MAX_n_smoothed_trend;
localVars.MAX_n_trend_coeff=MAX_n_trend_coeff;
@ -246,6 +259,10 @@ else
ifil(10,j+1) =ifil(10,j)+nfiles;
nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_trend);
ifil(11,j+1) =ifil(11,j)+nfiles;
if smoothed_state_uncertainty
nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_state_uncertainty);
ifil(13,j+1) =ifil(13,j)+nfiles;
end
end
end
localVars.ifil = ifil;
@ -297,7 +314,12 @@ if options_.smoother
pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',...
'',M_.endo_names(1:M_.orig_endo_nbr, :),M_.endo_names_tex,M_.endo_names,...
varlist,'Trend',DirectoryName,'_smoothed_trend');
if smoothed_state_uncertainty
pm3(endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',...
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'StateUncertainty',DirectoryName,'_state_uncert');
end
if nvn
for obs_iter=1:length(options_.varobs)
meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names(strmatch(options_.varobs{obs_iter},M_.endo_names,'exact'),:)];

32
matlab/prior_posterior_statistics_core.m
View File

@ -21,6 +21,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
% _trend_coeff;
% _smoothed_trend;
% _smoothed_constant;
% _state_uncert;
%
% ALGORITHM
% Portion of prior_posterior.m function.
@ -58,6 +59,7 @@ end
type=myinputs.type;
run_smoother=myinputs.run_smoother;
filter_covariance=myinputs.filter_covariance;
smoothed_state_uncertainty=myinputs.smoothed_state_uncertainty;
gend=myinputs.gend;
Y=myinputs.Y;
data_index=myinputs.data_index;
@ -83,8 +85,12 @@ if naK
MAX_naK=myinputs.MAX_naK;
end
if filter_covariance
MAX_filter_covariance=myinputs.MAX_filter_covariance;
MAX_filter_covariance=myinputs.MAX_filter_covariance;
end
if smoothed_state_uncertainty
MAX_n_smoothed_state_uncertainty=myinputs.MAX_n_smoothed_state_uncertainty;
end
exo_nbr=myinputs.exo_nbr;
maxlag=myinputs.maxlag;
MAX_nsmoo=myinputs.MAX_nsmoo;
@ -177,6 +183,9 @@ stock_ys = NaN(MAX_nruns,endo_nbr);
if filter_covariance
stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance);
end
if smoothed_state_uncertainty
stock_smoothed_uncert = zeros(endo_nbr,endo_nbr,gend,MAX_n_smoothed_state_uncertainty);
end
for b=fpar:B
if strcmpi(type,'prior')
@ -195,7 +204,7 @@ for b=fpar:B
if run_smoother
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition] = ...
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition,state_uncertainty] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value);
stock_trend_coeff(options_.varobs_id,irun(9))=trend_coeff;
@ -322,6 +331,9 @@ for b=fpar:B
if filter_covariance
stock_filter_covariance(dr.order_var,dr.order_var,:,irun(8)) = P;
end
if smoothed_state_uncertainty
stock_smoothed_uncert(dr.order_var,dr.order_var,:,irun(13)) = state_uncertainty;
end
else
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
end
@ -330,7 +342,7 @@ for b=fpar:B
stock_ys(irun(5),:) = SteadyState';
irun = irun + ones(12,1);
irun = irun + ones(13,1);
if run_smoother && (irun(1) > MAX_nsmoo || b == B),
@ -463,6 +475,17 @@ for b=fpar:B
irun(irun_index) = 1;
end
irun_index=13;
if run_smoother && smoothed_state_uncertainty && (irun(irun_index) > MAX_n_smoothed_state_uncertainty || b == B)
stock = stock_smoothed_uncert(:,:,:,1:irun(irun_index)-1);
ifil(irun_index) = ifil(irun_index) + 1;
save([DirectoryName '/' M_.fname '_state_uncert' int2str(ifil(irun_index)) '.mat'],'stock');
if RemoteFlag==1,
OutputFileName_state_uncert = [OutputFileName_state_uncert; {[DirectoryName filesep], [M_.fname '_state_uncert' int2str(ifil(irun_index)) '.mat']}];
end
irun(irun_index) = 1;
end
dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
end
@ -480,7 +503,8 @@ if RemoteFlag==1,
OutputFileName_filter_covar;
OutputFileName_trend_coeff;
OutputFileName_smoothed_trend;
OutputFileName_smoothed_constant];
OutputFileName_smoothed_constant;
OutputFileName_state_uncert];
end
dyn_waitbar_close(h);

13
matlab/store_smoother_results.m
View File

@ -1,4 +1,4 @@
function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend)
function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty)
% oo_=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend)
% Writes the smoother results into respective fields in oo_
%
@ -24,6 +24,8 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da
% decomp [4D array] (K*m*r*(T+K)) 4D array of shock decomposition of k-step ahead
% filtered variables (decision-rule order)
% Trend [double] [nvarobs*T] matrix of trends in observables
% state_uncertainty [double] (K,K,T) array, storing the uncertainty
% about the smoothed state (decision-rule order)
%
% Outputs:
% oo_ [structure] storing the results:
@ -38,7 +40,8 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da
% oo_.UpdatedVariables: structure storing the updated variables
% oo_.SmoothedShocks: structure storing the smoothed shocks
% oo_.SmoothedMeasurementErrors: structure storing the smoothed measurement errors
%
% oo_.Smoother.State_uncertainty: smoothed state uncertainty (declaration order)
% yf [double] (nvarobs*T) matrix storing the smoothed observed variables (order of options_.varobs)
%
% Notes:
@ -117,6 +120,9 @@ if options_.filter_covariance
oo_.Smoother.Variance = P;
end
if options_.smoothed_state_uncertainty
oo_.Smoother.State_uncertainty=state_uncertainty;
end
%get indices of smoothed variables
i_endo_in_bayestopt_smoother_varlist = bayestopt_.smoother_saved_var_list;
i_endo_in_dr_matrices=bayestopt_.smoother_var_list(i_endo_in_bayestopt_smoother_varlist);
@ -207,6 +213,9 @@ end
if options_.filter_covariance
oo_.Smoother.Variance(oo_.dr.order_var,oo_.dr.order_var,:)=oo_.Smoother.Variance;
end
if options_.smoothed_state_uncertainty
oo_.Smoother.State_uncertainty(oo_.dr.order_var,oo_.dr.order_var,:)=state_uncertainty;
end
%% get smoothed shocks
for exo_iter=1:M_.exo_nbr

7
preprocessor/DynareBison.yy
View File

@ -153,7 +153,7 @@ class ParsingDriver;
%token <string_val> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX DIRICHLET
%token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD WEIBULL WEIBULL_PDF
%token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION
%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION SMOOTHED_STATE_UNCERTAINTY
%token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS POSTERIOR_SAMPLER_OPTIONS
%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON
%token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION
@ -1777,6 +1777,7 @@ estimation_options : o_datafile
| o_partial_information
| o_filter_covariance
| o_filter_decomposition
| o_smoothed_state_uncertainty
| o_selected_variables_only
| o_conditional_variance_decomposition
| o_cova_compute
@ -2581,6 +2582,7 @@ calib_smoother_option : o_filtered_vars
| o_filter_decomposition
| o_diffuse_kalman_tol
| o_diffuse_filter
| o_smoothed_state_uncertainty
;
extended_path : EXTENDED_PATH ';'
@ -3129,6 +3131,9 @@ o_filter_covariance : FILTER_COVARIANCE
o_filter_decomposition : FILTER_DECOMPOSITION
{ driver.option_num("filter_decomposition","1");}
;
o_smoothed_state_uncertainty : SMOOTHED_STATE_UNCERTAINTY
{ driver.option_num("smoothed_state_uncertainty","1");}
;
o_selected_variables_only : SELECTED_VARIABLES_ONLY
{ driver.option_num("selected_variables_only","1");}
;

1
preprocessor/DynareFlex.ll
View File

@ -573,6 +573,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>k_order_solver {return token::K_ORDER_SOLVER; }
<DYNARE_STATEMENT>filter_covariance {return token::FILTER_COVARIANCE; }
<DYNARE_STATEMENT>filter_decomposition {return token::FILTER_DECOMPOSITION; }
<DYNARE_STATEMENT>smoothed_state_uncertainty {return token::SMOOTHED_STATE_UNCERTAINTY; }
<DYNARE_STATEMENT>selected_variables_only {return token::SELECTED_VARIABLES_ONLY; }
<DYNARE_STATEMENT>pruning {return token::PRUNING; }
<DYNARE_STATEMENT>save_draws {return token::SAVE_DRAWS; }

3
tests/TeX/fs2000_corr_ME.mod
View File

@ -168,7 +168,8 @@ end;
write_latex_prior_table;
estimation(mode_compute=8,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=4000, mh_nblocks=1, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20,moments_varendo,contemporaneous_correlation,conditional_variance_decomposition=[1 2 4]) m P c e W R k d y gy_obs;
estimation(mode_compute=8,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=4000, mh_nblocks=1, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20,
moments_varendo,contemporaneous_correlation,conditional_variance_decomposition=[1 2 4],smoothed_state_uncertainty) m P c e W R k d y gy_obs;
trace_plot(options_,M_,estim_params_,'PosteriorDensity',1);
trace_plot(options_,M_,estim_params_,'StructuralShock',1,'e_a')

4
tests/kalman_filter_smoother/algo1.mod
View File

@ -31,8 +31,8 @@ stderr e_z, uniform_pdf,,, 0.01, 0.1;
end;
varobs dw dx dy z;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

4
tests/kalman_filter_smoother/algo2.mod
View File

@ -31,8 +31,8 @@ stderr e_z, uniform_pdf,,, 0.01, 0.1;
end;
varobs dw dx dy z;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo1_mode,kalman_algo=2,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo1_mode,kalman_algo=2,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

4
tests/kalman_filter_smoother/algo3.mod
View File

@ -34,8 +34,8 @@ stderr e_z, inv_gamma_pdf,0.01, inf;
end;
varobs w x y;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

4
tests/kalman_filter_smoother/algo4.mod
View File

@ -34,8 +34,8 @@ stderr e_z, inv_gamma_pdf,0.01, inf;
end;
varobs w x y;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter,kalman_algo=4,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter,kalman_algo=4,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

2
tests/kalman_filter_smoother/algo4a.mod
View File

@ -33,7 +33,7 @@ end;
varobs dw dx y z;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty);
//estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter);
//checking smoother consistency

2
tests/kalman_filter_smoother/algo4b.mod
View File

@ -33,7 +33,7 @@ end;
varobs dw dx y z;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty);
//estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algo3_mode,diffuse_filter);
//checking smoother consistency

2
tests/kalman_filter_smoother/algoH1.mod
View File

@ -34,7 +34,7 @@ end;
varobs dw dx dy z;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

2
tests/kalman_filter_smoother/algoH2.mod
View File

@ -35,7 +35,7 @@ end;
varobs dw dx dy z;
//estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,kalman_algo=2);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algoH1_mode,kalman_algo=2,filtered_vars);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,mode_compute=0,mode_file=algoH1_mode,kalman_algo=2,filtered_vars,smoothed_state_uncertainty);
//checking smoother consistency
X = oo_.SmoothedVariables;

2
tests/kalman_filter_smoother/algoH3.mod
View File

@ -37,7 +37,7 @@ end;
varobs w x y;
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter);
estimation(datafile=data,first_obs=1000,nobs=200,mh_replic=0,diffuse_filter,smoothed_state_uncertainty);
stoch_simul(irf=0);

10
tests/kalman_filter_smoother/fs2000_smoother_only.mod
View File

@ -100,11 +100,11 @@ check;
varobs gp_obs gy_obs;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=1) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=2) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=3) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=4) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=1, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=2, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=3, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear, smoother,kalman_algo=4, smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
/*
* The following lines were used to generate the data file. If you want to

6
tests/kalman_filter_smoother/fs2000_smoother_only_ns.mod
View File

@ -107,9 +107,9 @@ P_obs (log(mst)-gam);
Y_obs (gam);
end;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=3) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=4) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=3,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
estimation(order=1, datafile=fsdat_simul, mode_compute=0,nobs=192, loglinear,diffuse_filter, smoother,kalman_algo=4,smoothed_state_uncertainty) m P c e W R k d n l gy_obs gp_obs y dA;
/*
* The following lines were used to generate the data file. If you want to