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

Implement heteroskedastic filter and smoother

time-shift
Marcoo Ratto 2021-05-02 17:07:12 +02:00 committed by Sébastien Villemot
parent 0d549bd87a
commit 1645f38269
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
24 changed files with 570 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
doc/manual/source/the-model-file.rst
View File

@ -2555,6 +2555,59 @@ blocks.
See above for the meaning of the ``overwrite`` option.
.. block:: heteroskedastic_shocks ;
heteroskedastic_shocks(overwrite);
|br| In *estimation context*, it implements heteroskedastic filters, where the standard error of shocks may unexpectedly change in every period.
The standard deviation of shocks may be either provided directly or set/modified in each observed period by a scale factor.
If ``std0`` is the usual standard error for ``shock1``, then:
* using a scale factor in period ``t`` implies: ``std(shock1|t)=std0(shock1)*scale(t)``
* using a provided value in period ``t`` implies: ``std(shock1|t)=value(t)``.
The block has a similar syntax as the ``shocks`` block in a perfect foresight context.
It should contain one or more occurrences of the following
group of three lines (for setting values)::
var VARIABLE_NAME;
periods INTEGER[:INTEGER] [[,] INTEGER[:INTEGER]]...;
values DOUBLE | (EXPRESSION) [[,] DOUBLE | (EXPRESSION) ]...;
OR (for setting scale factors)::
var VARIABLE_NAME;
periods INTEGER[:INTEGER] [[,] INTEGER[:INTEGER]]...;
scales DOUBLE | (EXPRESSION) [[,] DOUBLE | (EXPRESSION) ]...;
NOTE: ``scales`` and ``values`` cannot be simultaneously set for the same shock in the same period, but it is
possible to set ``values`` for some periods and ``scales`` for other periods for the same shock.
*Example*
::
heteroskedastic_shocks;
var e1;
periods 86:87 88:97;
scales 0.5 0;
var e2;
periods 86:87 88:97;
values 0.04 0.01;
end;
var e3;
periods 86:87;
values 0.04;
end;
var e3;
periods 88:97;
scales 0;
end;
.. specvar:: Sigma_e
|br| This special variable specifies directly the covariance
@ -6429,6 +6482,11 @@ block decomposition of the model (see :opt:`block`).
integrated processes (their first difference or k-difference
must be stationary).
.. option:: heteroskedastic_filter
Runs filter, likelihood, and smoother using heteroskedastic definitions provided in
a ``heteroskedastic_shocks`` block.
.. option:: selected_variables_only
Only run the classical smoother on the variables listed just

16
matlab/DsgeSmoother.m
View File

@ -171,7 +171,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo ~= 4
kalman_algo = 3;
else
if ~all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
if ~all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is not diagonal...
%Augment state vector (follows Section 6.4.3 of DK (2012))
expanded_state_vector_for_univariate_filter=1;
T = blkdiag(T,zeros(vobs));
@ -225,6 +225,10 @@ end
ST = T;
R1 = R;
if options_.heteroskedastic_filter
Q=get_Qvec_heteroskedastic_filter(Q,smpl,M_);
end
if kalman_algo == 1 || kalman_algo == 3
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
@ -255,7 +259,15 @@ if kalman_algo == 2 || kalman_algo == 4
Z = [Z, eye(vobs)];
ST = blkdiag(ST,zeros(vobs));
np = size(ST,1);
Q = blkdiag(Q,H);
if options_.heteroskedastic_filter
Qvec=Q;
Q=NaN(size(Qvec,1)+size(H,1),size(Qvec,1)+size(H,1),smpl+1);
for kv=1:size(Qvec,3)
Q(:,:,kv) = blkdiag(Qvec(:,:,kv),H);
end
else
Q = blkdiag(Q,H);
end
R1 = blkdiag(R,eye(vobs));
if kalman_algo == 4
%recompute Schur state space transformation with

1
matlab/default_option_values.m
View File

@ -355,6 +355,7 @@ options_.dataset.nobs = NaN;
options_.dataset.xls_sheet = [];
options_.dataset.xls_range = [];
options_.Harvey_scale_factor = 10;
options_.heteroskedastic_filter = false;
options_.MaxNumberOfBytes = 1e8;
options_.MaximumNumberOfMegaBytes = 111;
options_.analytic_derivation = 0; % Not a boolean, can also take values -1 or 2

39
matlab/dsge_likelihood.m
View File

@ -257,6 +257,7 @@ kalman_tol = DynareOptions.kalman_tol;
diffuse_kalman_tol = DynareOptions.diffuse_kalman_tol;
riccati_tol = DynareOptions.riccati_tol;
Y = transpose(DynareDataset.data)-trend;
smpl = size(Y,2);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
@ -267,6 +268,10 @@ kalman_algo = DynareOptions.kalman_algo;
diffuse_periods = 0;
expanded_state_vector_for_univariate_filter=0;
singular_diffuse_filter = 0;
if DynareOptions.heteroskedastic_filter
Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model);
end
switch DynareOptions.lik_init
case 1% Standard initialization with the steady state of the state equation.
if kalman_algo~=2
@ -304,6 +309,11 @@ switch DynareOptions.lik_init
Z(i,BayesInfo.mf(i))=1;
end
Zflag = 1;
if DynareOptions.heteroskedastic_filter
QQ=Qvec;
else
QQ=Q;
end
% Run diffuse kalman filter on first periods.
if (kalman_algo==3)
% Multivariate Diffuse Kalman Filter
@ -315,13 +325,13 @@ switch DynareOptions.lik_init
[dLIK,dlik,a_0_given_tm1,Pstar] = kalman_filter_d(Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
T,R,QQ,H,Z,mm,pp,rr);
else
[dLIK,dlik,a_0_given_tm1,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
T,R,QQ,H,Z,mm,pp,rr);
end
diffuse_periods = length(dlik);
if isinf(dLIK)
@ -359,7 +369,15 @@ switch DynareOptions.lik_init
Pinf = blkdiag(Pinf,zeros(pp));
H1 = zeros(pp,1);
mmm = mm+pp;
if DynareOptions.heteroskedastic_filter
clear QQ
for kv=1:size(Qvec,3)
QQ(:,:,kv) = blkdiag(Qvec(:,:,kv),H);
end
Qvec=QQ;
else
QQ = Q;
end
end
end
@ -372,7 +390,7 @@ switch DynareOptions.lik_init
Y, 1, size(Y,2), ...
a_0_given_tm1, Pinf, Pstar, ...
kalman_tol, diffuse_kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H1,Z,mmm,pp,rr);
T,R,QQ,H1,Z,mmm,pp,rr);
diffuse_periods = size(dlik,1);
end
if isnan(dLIK)
@ -579,6 +597,9 @@ end
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if DynareOptions.heteroskedastic_filter
Q=Qvec;
end
singularity_has_been_detected = false;
% First test multivariate filter if specified; potentially abort and use univariate filter instead
@ -689,7 +710,14 @@ if (kalman_algo==2) || (kalman_algo==4)
Z = [Z1, eye(pp)];
Zflag=1;
T = blkdiag(T,zeros(pp));
Q = blkdiag(Q,H);
if DynareOptions.heteroskedastic_filter
clear Q
for kv=1:size(Qvec,3)
Q(:,:,kv) = blkdiag(Qvec(:,:,kv),H);
end
else
Q = blkdiag(Q,H);
end
R = blkdiag(R,eye(pp));
Pstar = blkdiag(Pstar,H);
Pinf = blkdiag(Pinf,zeros(pp));
@ -866,3 +894,4 @@ if isfield(M_,'filter_initial_state') && ~isempty(M_.filter_initial_state)
end
end
end

34
matlab/dynare_estimation_init.m
View File

@ -652,3 +652,37 @@ if options_.mh_replic || options_.load_mh_file
[current_options, options_] = check_posterior_sampler_options([], options_, bounds);
options_.posterior_sampler_options.current_options = current_options;
end
%% set heteroskedastic shocks
if options_.heteroskedastic_filter
if options_.fast_kalman_filter
error('estimation option conflict: "heteroskedastic_filter" incompatible with "fast_kalman_filter"')
end
if options_.block
error('estimation option conflict: "heteroskedastic_filter" incompatible with block kalman filter')
end
if options_.analytic_derivation
error(['estimation option conflict: analytic_derivation isn''t available ' ...
'for heteroskedastic_filter'])
end
M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs);
M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs);
xname = fieldnames(M_.heteroskedastic_shocks.Qhet);
for k=1:length(xname)
inx = strcmp(xname{k},M_.exo_names);
isqscale=false;
if isfield(M_.heteroskedastic_shocks.Qhet.(xname{k}),'scale')
M_.heteroskedastic_shocks.Qscale(inx,M_.heteroskedastic_shocks.Qhet.(xname{k}).time_scale)=M_.heteroskedastic_shocks.Qhet.(xname{k}).scale.^2;
isqscale=true;
end
if isfield(M_.heteroskedastic_shocks.Qhet.(xname{k}),'value')
if isqscale && ~isempty(intersect(M_.heteroskedastic_shocks.Qhet.(xname{k}).time_value,M_.heteroskedastic_shocks.Qhet.(xname{k}).time_scale))
fprintf('\ndynare_estimation_init: With the option "heteroskedastic_shocks" you cannot define\n')
fprintf('dynare_estimation_init: the scale and the value for the same shock \n')
fprintf('dynare_estimation_init: in the same period!\n')
error('Scale and value defined for the same shock in the same period with "heteroskedastic_shocks".')
end
M_.heteroskedastic_shocks.Qvalue(inx,M_.heteroskedastic_shocks.Qhet.(xname{k}).time_value)=M_.heteroskedastic_shocks.Qhet.(xname{k}).value.^2;
end
end
end

26
matlab/initial_estimation_checks.m
View File

@ -42,6 +42,22 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy
maximum_number_non_missing_observations=max(sum(~isnan(DynareDataset.data(2:end,:)),2));
init_number_non_missing_observations=sum(~isnan(DynareDataset.data(1,:)),2);
if DynareOptions.heteroskedastic_filter
if DynareOptions.order>1
error('initial_estimation_checks:: heteroskedastic shocks are only supported with the Kalman filter/smoother')
end
observations_by_period=sum(~isnan(DynareDataset.data),2);
base_shocks = find(diag(Model.Sigma_e));
zero_shocks = ~ismember(1:Model.exo_nbr,base_shocks);
non_zero_shocks_by_period=repmat(length(base_shocks),size(observations_by_period));
% check periods for which base shocks are scaled to zero
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(Model.heteroskedastic_shocks.Qscale(base_shocks,1:DynareOptions.nobs)==0,1)';
% check periods for which base shocks are set to zero
non_zero_shocks_by_period = non_zero_shocks_by_period-sum(Model.heteroskedastic_shocks.Qvalue(base_shocks,1:DynareOptions.nobs)==0,1)';
% check periods for which other shocks are set to a positive number
non_zero_shocks_by_period = non_zero_shocks_by_period+sum(Model.heteroskedastic_shocks.Qvalue(zero_shocks,1:DynareOptions.nobs)>0,1)';
end
if DynareOptions.order>1
if any(any(isnan(DynareDataset.data)))
error('initial_estimation_checks:: particle filtering does not support missing observations')
@ -95,8 +111,14 @@ if init_number_non_missing_observations>Model.exo_nbr+non_zero_ME
end
end
if maximum_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
if DynareOptions.heteroskedastic_filter
if any(observations_by_period>(non_zero_shocks_by_period+non_zero_ME))
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance: Check heteroskedastic block and shocks calibration!'])
end
else
if maximum_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
end
end
if init_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
if DynareOptions.no_init_estimation_check_first_obs

35
matlab/kalman/get_Qvec_heteroskedastic_filter.m Normal file
View File

@ -0,0 +1,35 @@
function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,Model)
isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal...
Qvec=repmat(Q,[1 1 smpl+1]);
for k=1:smpl
inx = ~isnan(Model.heteroskedastic_shocks.Qvalue(:,k));
if any(inx)
if isqdiag
Qvec(inx,inx,k)=diag(Model.heteroskedastic_shocks.Qvalue(inx,k));
else
inx = find(inx);
for s=1:length(inx)
if Q(inx(s),inx(s))>1.e-14
tmpscale = sqrt(Model.heteroskedastic_shocks.Qvalue(inx(s),k)./Q(inx(s),inx(s)));
Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale;
Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale;
else
Qvec(inx(s),inx(s),k)=Model.heteroskedastic_shocks.Qvalue(inx(s),k);
end
end
end
end
inx = ~isnan(Model.heteroskedastic_shocks.Qscale(:,k));
if any(inx)
if isqdiag
Qvec(inx,inx,k)=Qvec(inx,inx,k).*diag(Model.heteroskedastic_shocks.Qscale(inx,k));
else
inx = find(inx);
for s=1:length(inx)
tmpscale = sqrt(Model.heteroskedastic_shocks.Qscale(inx(s),k));
Qvec(inx(s),:,k) = Qvec(inx(s),:,k).*tmpscale;
Qvec(:,inx(s),k) = Qvec(:,inx(s),k).*tmpscale;
end
end
end
end

15
matlab/kalman/likelihood/kalman_filter.m
View File

@ -71,7 +71,7 @@ function [LIK, LIKK, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_t
%! @end deftypefn
%@eod:
% Copyright (C) 2004-2017 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -114,6 +114,12 @@ smpl = last-start+1;
% Initialize some variables.
dF = 1;
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = start; % Initialization of the time index.
likk = zeros(smpl,1); % Initialization of the vector gathering the densities.
@ -161,6 +167,9 @@ end
rescale_prediction_error_covariance0=rescale_prediction_error_covariance;
while notsteady && t<=last
s = t-start+1;
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if Zflag
v = Y(:,t)-Z*a;
F = Z*P*Z' + H;
@ -227,7 +236,9 @@ while notsteady && t<=last
end
a = T*tmp;
P = Ptmp;
notsteady = max(abs(K(:)-oldK))>riccati_tol;
if ~isqvec
notsteady = max(abs(K(:)-oldK))>riccati_tol;
end
oldK = K(:);
end
t = t+1;

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

@ -37,7 +37,7 @@ function [dLIK,dlik,a,Pstar] = kalman_filter_d(Y, start, last, a, Pinf, Pstar, k
% Durbin/Koopman (2012): "Time Series Analysis by State Space Methods", Oxford University Press,
% Second Edition, Ch. 5 and 7.2
% Copyright (C) 2004-2018 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -59,6 +59,12 @@ smpl = last-start+1;
% Initialize some variables.
dF = 1;
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
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.
@ -70,6 +76,9 @@ while rank(Z*Pinf*Z',diffuse_kalman_tol) && (t<=last)
s = t-start+1;
v = Y(:,t)-Z*a; %get prediction error v^(0) in (5.13) DK (2012)
Finf = Z*Pinf*Z'; % (5.7) in DK (2012)
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
%do case distinction based on whether F_{\infty,t} has full rank or 0 rank
if rcond(Finf) < diffuse_kalman_tol %F_{\infty,t} = 0
if ~all(abs(Finf(:)) < diffuse_kalman_tol) %rank-deficient but not rank 0

13
matlab/kalman/likelihood/missing_observations_kalman_filter.m
View File

@ -32,7 +32,7 @@ function [LIK, lik, a, P] = missing_observations_kalman_filter(data_index,numbe
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2004-2017 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -76,6 +76,12 @@ smpl = last-start+1;
% Initialize some variables.
dF = 1;
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
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.
@ -89,6 +95,9 @@ rescale_prediction_error_covariance0=rescale_prediction_error_covariance;
while notsteady && t<=last
s = t-start+1;
d_index = data_index{t};
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if isempty(d_index)
a = T*a;
P = T*P*transpose(T)+QQ;
@ -147,7 +156,7 @@ while notsteady && t<=last
P = T*(P-K*P(z,:))*transpose(T)+QQ;
end
a = T*(a+K*v);
if t>=no_more_missing_observations
if t>=no_more_missing_observations && ~isqvec
notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
end

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

@ -42,7 +42,7 @@ function [dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(data_index,n
% Second Edition, Ch. 5 and 7.2
%
% Copyright (C) 2004-2017 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -65,6 +65,12 @@ smpl = last-start+1;
% Initialize some variables.
dF = 1;
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
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.
@ -79,6 +85,9 @@ s = 0;
while rank(Pinf,diffuse_kalman_tol) && (t<=last)
s = t-start+1;
d_index = data_index{t};
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if isempty(d_index)
%no observations, propagate forward without updating based on
%observations

13
matlab/kalman/likelihood/univariate_kalman_filter.m
View File

@ -84,7 +84,7 @@ function [LIK, lik,a,P] = univariate_kalman_filter(data_index,number_of_observat
% Second Edition, Ch. 6.4 + 7.2.5
% Copyright (C) 2004-2017 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -116,6 +116,12 @@ end
smpl = last-start+1;
% Initialize some variables.
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = start; % Initialization of the time index.
lik = zeros(smpl,pp); % Initialization of the matrix gathering the densities at each time and each observable
@ -164,6 +170,9 @@ end
while notsteady && t<=last %loop over t
s = t-start+1;
d_index = data_index{t};
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if Zflag
z = Z(d_index,:);
else
@ -216,7 +225,7 @@ while notsteady && t<=last %loop over t
end
a = T*a; %transition according to (6.14) in DK (2012)
P = T*P*T' + QQ; %transition according to (6.14) in DK (2012)
if t>=no_more_missing_observations
if t>=no_more_missing_observations && ~isqvec
notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
end

11
matlab/kalman/likelihood/univariate_kalman_filter_d.m
View File

@ -87,7 +87,7 @@ function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index,
% Series Analysis by State Space Methods", Oxford University Press,
% Second Edition, Ch. 5, 6.4 + 7.2.5
% Copyright (C) 2004-2018 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -109,6 +109,12 @@ smpl = last-start+1;
% Initialize some variables.
dF = 1;
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = start; % Initialization of the time index.
dlikk= zeros(smpl,1); % Initialization of the vector gathering the densities.
@ -159,6 +165,9 @@ while newRank && (t<=last)
oldRank = 0;
end
a = T*a;
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
Pstar = T*Pstar*T'+QQ;
Pinf = T*Pinf*T';
if newRank

21
matlab/missing_DiffuseKalmanSmootherH1_Z.m
View File

@ -99,6 +99,12 @@ Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl+1);
Pinf(:,:,1) = Pinf1;
rr = size(Q,1);
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
@ -117,6 +123,9 @@ t = 0;
while rank(Pinf(:,:,t+1),diffuse_kalman_tol) && t<smpl
t = t+1;
di = data_index{t};
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if isempty(di)
%no observations, propagate forward without updating based on
%observations
@ -189,6 +198,9 @@ while notsteady && t<smpl
t = t+1;
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); % make sure P is symmetric
di = data_index{t};
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
if isempty(di)
atilde(:,t) = a(:,t);
L(:,:,t) = T;
@ -265,6 +277,9 @@ while t>d+1
end
end
alphahat(:,t) = a(:,t) + P(:,:,t)*r(:,t); %DK (2012), eq. 4.35
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
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
@ -312,6 +327,9 @@ if d %diffuse periods
end
end
alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); % DK (2012), eq. 5.23
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
etahat(:,t) = QRt*r0(:,t); % DK (2012), p. 135
if state_uncertainty_flag
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N(:,:,t)*Pstar(:,:,t)...
@ -328,6 +346,9 @@ if decomp_flag
for t = max(d,1):smpl
di = data_index{t};
% calculate eta_tm1t
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
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

23
matlab/missing_DiffuseKalmanSmootherH3_Z.m
View File

@ -57,7 +57,7 @@ function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalma
% Models", S.J. Koopman and J. Durbin (2003), in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98.
% Copyright (C) 2004-2018 Dynare Team
% Copyright (C) 2004-2021 Dynare Team
%
% This file is part of Dynare.
%
@ -110,6 +110,12 @@ Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
rr = size(Q,1); % number of structural shocks
isqvec = false;
if ndim(Q)>2
Qvec = Q;
Q=Q(:,:,1);
isqvec = true;
end
QQ = R*Q*transpose(R);
QRt = Q*transpose(R);
alphahat = zeros(mm,smpl);
@ -176,6 +182,9 @@ while newRank && t < smpl
for jnk=2:nk
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
end
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T';
if newRank
@ -220,6 +229,9 @@ while notsteady && t<smpl
% p. 157, DK (2012)
end
end
if isqvec
QQ = R*Qvec(:,:,t+1)*transpose(R);
end
a1(:,t+1) = T*a(:,t); %transition according to (6.14) in DK (2012)
P(:,:,t+1) = T*P(:,:,t)*T' + QQ; %transition according to (6.14) in DK (2012)
Pf = P(:,:,t+1);
@ -288,6 +300,9 @@ while t > d+1
end
r(:,t) = ri; % DK (2012), below 6.15, r_{t-1}=r_{t,0}
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t);
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
etahat(:,t) = QRt*r(:,t);
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
if state_uncertainty_flag
@ -339,6 +354,9 @@ if d
end
alphahat(:,t) = a1(:,t) + Pstar1(:,:,t)*r0(:,t) + Pinf1(:,:,t)*r1(:,t); % KD (2000), eq. (26)
r(:,t) = r0(:,t);
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
etahat(:,t) = QRt*r(:,t); % KD (2000), eq. (27)
if state_uncertainty_flag
V(:,:,t)=Pstar(:,:,t)-Pstar(:,:,t)*N_0(:,:,t)*Pstar(:,:,t)...
@ -371,6 +389,9 @@ if decomp_flag
end
% calculate eta_tm1t
if isqvec
QRt = Qvec(:,:,t)*transpose(R);
end
eta_tm1t = QRt*ri_d;
% calculate decomposition
Ttok = eye(mm,mm);

2
preprocessor

@ -1 +1 @@
Subproject commit 97e9609c28fb4bda15b699da033f02fe02ca6a0e
Subproject commit b12b254f2c90c92ba8218d07a9860762d41be0fd

4
scripts/dynare.el
View File

@ -82,7 +82,7 @@
;; Keywords that may appear in blocks, and that begin a statement which will be
;; closed by a semicolon
(defvar dynare-statements-like
'("stderr" "values" "restriction" "exclusion" "upper_cholesky" "lower_cholesky")
'("stderr" "values" "scales" "restriction" "exclusion" "upper_cholesky" "lower_cholesky")
"Dynare statements-like keywords.")
;; Those keywords that makes the lexer enter the DYNARE_BLOCK start condition
@ -91,7 +91,7 @@
;; referenced in another eval-when-compile statement in dynare-calculate-indentation
(eval-when-compile
(defvar dynare-blocks
'("model" "steady_state_model" "initval" "endval" "histval" "shocks"
'("model" "steady_state_model" "initval" "endval" "histval" "shocks" "heteroskedastic_shocks"
"shock_groups" "init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
"estimated_params_init" "estimated_params_bounds" "osr_params_bounds"
"observation_trends" "deterministic_trends" "optim_weights" "homotopy_setup"

6
tests/Makefile.am
View File

@ -50,6 +50,8 @@ MODFILES = \
estimation/MH_recover/fs2000_recover.mod \
estimation/MH_recover/fs2000_recover_2.mod \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/heteroskedastic_shocks/fs2000_het.mod \
estimation/heteroskedastic_shocks/fs2000_het_corr.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/tune_mh_jscale/fs2000.mod \
estimation/method_of_moments/AnScho/AnScho_MoM.mod \
@ -566,6 +568,8 @@ XFAIL_MODFILES = ramst_xfail.mod \
optimal_policy/Ramsey/ramsey_ex_wrong_ss_file_xfail.mod \
estimation/fs2000_mixed_ML_xfail.mod \
estimation/fs2000_stochastic_singularity_xfail.mod \
estimation/heteroskedastic_shocks/fs2000_het_XFAIL.mod \
estimation/heteroskedastic_shocks/fs2000_het_XFAIL2.mod \
identification/ident_unit_root/ident_unit_root_xfail.mod \
identification/LindeTrabandt/LindeTrabandt2019_xfail.mod \
steady_state/Linear_steady_state_xfail.mod \
@ -1179,6 +1183,8 @@ EXTRA_DIST = \
kalman/likelihood_from_dynare/fs2000ns_estimation_check.inc \
identification/as2007/G_QT.mat \
estimation/fsdat_simul.m \
estimation/heteroskedastic_shocks/fs2000_het_model.inc \
estimation/heteroskedastic_shocks/fs2000_het_check.inc \
ep/mean_preserving_spread.m \
ep/rbcii_steady_state.m \
decision_rules/example1_results_dyn_432.mat \

33
tests/estimation/heteroskedastic_shocks/fs2000_het.mod Normal file
View File

@ -0,0 +1,33 @@
// See fs2000.mod in the examples/ directory for details on the model
// tests heteroskedastic filter/smoother
// includes lagged exogenous variable introduced by preprocessor
@#include "fs2000_het_model.inc"
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
stoch_simul(order=1,loglinear);
forecast;
options_.solve_tolf = 1e-12;
heteroskedastic_shocks;
var e_b;
periods 100:120;
values 0.01;
var e_a;
periods 100:120;
scales 0;
end;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_compute=5,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous,heteroskedastic_filter);
@#define mode_file_name="fs2000_het_mode"
@#include "fs2000_het_check.inc"

27
tests/estimation/heteroskedastic_shocks/fs2000_het_XFAIL.mod Normal file
View File

@ -0,0 +1,27 @@
// See fs2000.mod in the examples/ directory for details on the model
// Tests that stochastic singularity is correctly filtered out
@#include "fs2000_het_model.inc"
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
options_.solve_tolf = 1e-12;
heteroskedastic_shocks;
var e_b;
periods 100:120;
values 0.01;
var e_a;
periods 99:120;
scales 0;
end;
estimation(order=1,datafile='../fsdat_simul',nobs=192,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous,heteroskedastic_filter);

32
tests/estimation/heteroskedastic_shocks/fs2000_het_XFAIL2.mod Normal file
View File

@ -0,0 +1,32 @@
// See fs2000.mod in the examples/ directory for details on the model
// Tests that setting scale and value in same period for same shock is correctly filtered out
@#include "fs2000_het_model.inc"
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
options_.solve_tolf = 1e-12;
heteroskedastic_shocks;
var e_b;
periods 100:120;
values 0.01;
var e_a;
periods 100:120;
scales 0;
// Wrongly set scale on top of value for the same shock and period
var e_a;
periods 110;
values 0.01;
end;
estimation(order=1,datafile='../fsdat_simul',nobs=192,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous,heteroskedastic_filter);

55
tests/estimation/heteroskedastic_shocks/fs2000_het_check.inc Normal file
View File

@ -0,0 +1,55 @@
options_.debug=true;
options_.lik_init=1;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_file=@{mode_file_name},mode_compute=0,kalman_algo=1,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous);
fval_algo_1=oo_.likelihood_at_initial_parameters;
if isfield(oo_,'SmoothedMeasurementErrors')
SmoothedMeasurementErrors(:,:,1)=cell2mat(struct2cell(oo_.SmoothedMeasurementErrors));
end
SmoothedShocks(:,:,1)=cell2mat(struct2cell(oo_.SmoothedShocks));
SmoothedVariables(:,:,1)=cell2mat(struct2cell(oo_.SmoothedVariables));
options_.lik_init=1;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_file=@{mode_file_name},mode_compute=0,kalman_algo=2,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous);
fval_algo_2=oo_.likelihood_at_initial_parameters;
if isfield(oo_,'SmoothedMeasurementErrors')
SmoothedMeasurementErrors(:,:,2)=cell2mat(struct2cell(oo_.SmoothedMeasurementErrors));
end
SmoothedShocks(:,:,2)=cell2mat(struct2cell(oo_.SmoothedShocks));
SmoothedVariables(:,:,2)=cell2mat(struct2cell(oo_.SmoothedVariables));
options_.lik_init=1;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_file=@{mode_file_name},mode_compute=0,kalman_algo=3,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous);
fval_algo_3=oo_.likelihood_at_initial_parameters;
if isfield(oo_,'SmoothedMeasurementErrors')
SmoothedMeasurementErrors(:,:,3)=cell2mat(struct2cell(oo_.SmoothedMeasurementErrors));
end
SmoothedShocks(:,:,3)=cell2mat(struct2cell(oo_.SmoothedShocks));
SmoothedVariables(:,:,3)=cell2mat(struct2cell(oo_.SmoothedVariables));
options_.lik_init=1;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_file=@{mode_file_name},mode_compute=0,kalman_algo=4,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous);
fval_algo_4=oo_.likelihood_at_initial_parameters;
if isfield(oo_,'SmoothedMeasurementErrors')
SmoothedMeasurementErrors(:,:,4)=cell2mat(struct2cell(oo_.SmoothedMeasurementErrors));
end
SmoothedShocks(:,:,4)=cell2mat(struct2cell(oo_.SmoothedShocks));
SmoothedVariables(:,:,4)=cell2mat(struct2cell(oo_.SmoothedVariables));
if isfield(oo_,'SmoothedMeasurementErrors')
if max(max(abs(SmoothedMeasurementErrors-repmat(SmoothedMeasurementErrors(:,:,1),[1,1,4]))))>1e-6
error('SmoothedMeasurementErrors do not match')
end
end
if max(max(abs(SmoothedShocks-repmat(SmoothedShocks(:,:,1),[1,1,4]))))>1e-6
error('SmoothedShocks do not match')
end
if max(max(abs(SmoothedVariables-repmat(SmoothedVariables(:,:,1),[1,1,4]))))>1e-6
error('SmoothedVariables do not match')
end
if max(abs([fval_algo_2,fval_algo_3,fval_algo_4]-fval_algo_1))>1e-5
error('Likelihoods do not match')
end

37
tests/estimation/heteroskedastic_shocks/fs2000_het_corr.mod Normal file
View File

@ -0,0 +1,37 @@
// See fs2000.mod in the examples/ directory for details on the model
// tests heteroskedastic filter/smoother
// includes lagged exogenous variable introduced by preprocessor
@#include "fs2000_het_model.inc"
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
corr e_a, e_m = 0.01;
var gp_obs; stderr 0.001;
var gy_obs; stderr 0.001;
corr gp_obs, gy_obs = 0.01;
end;
steady;
check;
stoch_simul(order=1,loglinear);
forecast;
options_.solve_tolf = 1e-12;
heteroskedastic_shocks;
var e_b;
periods 100:120;
values 0.01;
var e_a;
periods 100:120;
scales 0;
end;
estimation(order=1,datafile='../fsdat_simul',nobs=192,mode_compute=5,loglinear,mh_replic=0,smoother,filtered_vars,forecast=8,filter_step_ahead=[1:8],consider_all_endogenous,heteroskedastic_filter);
@#define mode_file_name="fs2000_het_corr_mode"
@#include "fs2000_het_check.inc"

69
tests/estimation/heteroskedastic_shocks/fs2000_het_model.inc Normal file
View File

@ -0,0 +1,69 @@
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m e_b;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*exp(e_b)*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*exp(e_b)*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.1;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;