Commit cd3b5bf3 authored by Johannes Pfeifer's avatar Johannes Pfeifer Committed by Stéphane Adjemian
Browse files

Fix bug in Bayesian estimation where the filter_step_ahead command produced no output

- Improves documentation of resulting matrices
- Adds unit test for filter_step_ahead option
parent 39807643
......@@ -4924,6 +4924,9 @@ Default value is @code{0}.
forecast covariance matrices.
@item filter_step_ahead = [@var{INTEGER1}:@var{INTEGER2}]
See below.
@item filter_step_ahead = [@var{INTEGER1} @var{INTEGER2} @dots{}]
@anchor{filter_step_ahead}
@vindex oo_.FilteredVariablesKStepAhead
@vindex oo_.FilteredVariablesKStepAheadVariances
......@@ -5186,7 +5189,7 @@ Variable set by the @code{estimation} command, if it is used with the
@defvr {MATLAB/Octave variable} oo_.FilteredVariablesKStepAhead
Variable set by the @code{estimation} command, if it is used with the
@code{filter_step_ahead} option.
@code{filter_step_ahead} option. The k-steps are stored along the rows while the columns indicate the respective variables. The third dimension of the array provides the observation for which the forecast has been made. For example, if @code{filter_step_ahead=[1 2 4]} and @code{nobs=200}, the element (3,5,204) stores the four period ahead filtered value of variable 5 computed at time t=200 for time t=204. The periods at the beginning and end of the sample for which no forecasts can be made, e.g. entries (1,5,1) and (1,5,204) in the example, are set to zero.
@end defvr
@defvr {MATLAB/Octave variable} oo_.FilteredVariablesKStepAheadVariances
......@@ -5194,6 +5197,16 @@ Variable set by the @code{estimation} command, if it is used with the
@code{filter_step_ahead} option.
@end defvr
@defvr {MATLAB/Octave variable} oo_.Filtered_Variables_X_step_ahead
Variable set by the @code{estimation} command, if it is used with the @code{filter_step_ahead} option in the context of Bayesian estimation. Fields are of the form:
@example
@code{oo_.Filtered_Variables_X_step_ahead.@var{VARIABLE_NAME}}
@end example
The nth entry stores the k-step ahead filtered variable computed at time n for time n+k.
@end defvr
@defvr {MATLAB/Octave variable} oo_.PosteriorIRF.dsge
Variable set by the @code{estimation} command, if it is used with the
@code{bayesian_irf} option. Fields are of the form:
......
......@@ -56,29 +56,67 @@ if options_.TeX
end
Mean = zeros(n2,nvar);
Median = zeros(n2,nvar);
Std = zeros(n2,nvar);
Var = zeros(n2,nvar);
Distrib = zeros(9,n2,nvar);
HPD = zeros(2,n2,nvar);
fprintf(['Estimation::mcmc: ' tit1 '\n']);
stock1 = zeros(n1,n2,B);
k = 0;
filter_step_ahead_indicator=0;
for file = 1:ifil
load([DirectoryName '/' M_.fname var_type int2str(file)]);
if size(size(stock),2) == 4
stock = squeeze(stock(1,:,1:n2,:));
filter_step_ahead_indicator=1;
stock_filter_step_ahead=zeros(n1,n2,size(stock,4),length(options_.filter_step_ahead));
for ii=1:length(options_.filter_step_ahead)
K_step_ahead=options_.filter_step_ahead(ii);
stock_filter_step_ahead(:,:,:,ii)=stock(ii,:,1+K_step_ahead:n2+K_step_ahead,:);
end
stock = squeeze(stock(1,:,1+1:1+n2,:)); %1 step ahead starts at entry 2
end
k = k(end)+(1:size(stock,3));
stock1(:,:,k) = stock;
if filter_step_ahead_indicator
stock1_filter_step_ahead(:,:,k,:) = stock_filter_step_ahead;
end
end
clear stock
if filter_step_ahead_indicator
clear stock_filter_step_ahead
filter_steps=length(options_.filter_step_ahead);
Mean_filter_step_ahead = zeros(filter_steps,nvar,n2);
Median_filter_step_ahead = zeros(filter_steps,nvar,n2);
Var_filter_step_ahead = zeros(filter_steps,nvar,n2);
Distrib_filter_step_ahead = zeros(9,filter_steps,nvar,n2);
HPD_filter_step_ahead = zeros(2,filter_steps,nvar,n2);
end
tmp =zeros(B,1);
for i = 1:nvar
for j = 1:n2
[Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
if filter_step_ahead_indicator
for K_step = 1:length(options_.filter_step_ahead)
[Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
end
end
end
end
clear stock1
if filter_step_ahead_indicator %write matrices corresponding to ML
clear stock1_filter_step_ahead
FilteredVariablesKStepAhead=zeros(length(options_.filter_step_ahead),nvar,n2+max(options_.filter_step_ahead));
FilteredVariablesKStepAheadVariances=zeros(length(options_.filter_step_ahead),nvar,n2+max(options_.filter_step_ahead));
for K_step = 1:length(options_.filter_step_ahead)
FilteredVariablesKStepAhead(K_step,:,1+options_.filter_step_ahead(K_step):n2+options_.filter_step_ahead(K_step))=Mean_filter_step_ahead(K_step,:,:);
FilteredVariablesKStepAheadVariances(K_step,:,1+options_.filter_step_ahead(K_step):n2+options_.filter_step_ahead(K_step))=Mean_filter_step_ahead(K_step,:,:);
end
oo_.FilteredVariablesKStepAhead=FilteredVariablesKStepAhead;
oo_.FilteredVariablesKStepAheadVariances=FilteredVariablesKStepAheadVariances;
end
for i = 1:nvar
name = deblank(names1(SelecVariables(i),:));
eval(['oo_.' name3 '.Mean.' name ' = Mean(:,i);']);
......@@ -87,6 +125,17 @@ for i = 1:nvar
eval(['oo_.' name3 '.deciles.' name ' = Distrib(:,:,i);']);
eval(['oo_.' name3 '.HPDinf.' name ' = HPD(1,:,i);']);
eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i);']);
if filter_step_ahead_indicator
for K_step = 1:length(options_.filter_step_ahead)
name4=['Filtered_Variables_',num2str(K_step),'_step_ahead'];
eval(['oo_.' name4 '.Mean.' name ' = squeeze(Mean_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.Median.' name ' = squeeze(Median_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.Var.' name ' = squeeze(Var_filter_step_ahead(K_step,i,:));']);
eval(['oo_.' name4 '.deciles.' name ' = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));']);
eval(['oo_.' name4 '.HPDinf.' name ' = squeeze(HPD_filter_step_ahead(1,K_step,i,:));']);
eval(['oo_.' name4 '.HPDsup.' name ' = squeeze(HPD_filter_step_ahead(2,K_step,i,:));']);
end
end
end
%%
%% Finally I build the plots.
......
......@@ -300,7 +300,7 @@ if options_.filtered_vars
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'UpdatedVariables',DirectoryName, ...
'_update');
pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast (filtered variables)',...
pm3(endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',...
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
end
......
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
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*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*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;
initval;
k = 6;
m = mst;
P = 2.25;
c = 0.45;
e = 1;
W = 4;
R = 1.02;
d = 0.85;
n = 0.19;
l = 0.86;
y = 0.6;
gy_obs = exp(gam);
gp_obs = exp(-gam);
dA = exp(gam);
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
estimated_params;
alp, 0.356;
bet, 0.993;
gam, 0.0085;
end;
varobs gp_obs gy_obs;
estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, filter_step_ahead = [1 4 8 12], forecast=20,smoother,filtered_vars) m P c;
/*
* The following lines were used to generate the data file. If you want to
* generate another random data file, comment the "estimation" line and uncomment
* the following lines.
*/
//stoch_simul(periods=200, order=1);
//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs'));
% computes the steady state of fs2000 analyticaly
% largely inspired by the program of F. Schorfheide
% 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/>.
function [ys,check] = fs2000_steadystate(ys,exe)
global M_
alp = M_.params(1);
bet = M_.params(2);
gam = M_.params(3);
mst = M_.params(4);
rho = M_.params(5);
psi = M_.params(6);
del = M_.params(7);
check = 0;
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;
ys =[
m
P
c
e
W
R
k
d
n
l
gy_obs
gp_obs
y
dA ];
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
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*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*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;
initval;
k = 6;
m = mst;
P = 2.25;
c = 0.45;
e = 1;
W = 4;
R = 1.02;
d = 0.85;
n = 0.19;
l = 0.86;
y = 0.6;
gy_obs = exp(gam);
gp_obs = exp(-gam);
dA = exp(gam);
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
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.223;
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;
estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=2000, mh_nblocks=1, mh_jscale=0.8,filter_step_ahead = [1 4 8 12], forecast=20,smoother,filtered_vars) m P c;
/*
* The following lines were used to generate the data file. If you want to
* generate another random data file, comment the "estimation" line and uncomment
* the following lines.
*/
//stoch_simul(periods=200, order=1);
//datatomfile('fsdat_simul', char('gy_obs', 'gp_obs'));
% computes the steady state of fs2000 analyticaly
% largely inspired by the program of F. Schorfheide
% 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/>.
function [ys,check] = fs2000_steadystate(ys,exe)
global M_
alp = M_.params(1);
bet = M_.params(2);
gam = M_.params(3);
mst = M_.params(4);
rho = M_.params(5);
psi = M_.params(6);
del = M_.params(7);
check = 0;
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;
ys =[
m
P
c
e
W
R
k
d
n
l
gy_obs
gp_obs
y
dA ];
% Generated data, used by fs2000.mod
gy_obs =[
1.0030045
1.0002599
0.99104664
1.0321162
1.0223545
1.0043614
0.98626929
1.0092127
1.0357197
1.0150827
1.0051548
0.98465775
0.99132132
0.99904153
1.0044641
1.0179198
1.0113462
0.99409421
0.99904293
1.0448336
0.99932433
1.0057004
0.99619787
1.0267504
1.0077645
1.0058026
1.0025891
0.9939097
0.99604693
0.99908569
1.0151094
0.99348134
1.0039124
1.0145805
0.99800868
0.98578138
1.0065771
0.99843919
0.97979062
0.98413351
0.96468174
1.0273857
1.0225211
0.99958667
1.0111157
1.0099585
0.99480311
1.0079265
0.98924573
1.0070613
1.0075706
0.9937151
1.0224711
1.0018891
0.99051863
1.0042944
1.0184055
0.99419508
0.99756624
1.0015983
0.9845772
1.0004407
1.0116237
0.9861885
1.0073094
0.99273355
1.0013224
0.99777979
1.0301686
0.96809556
0.99917088
0.99949253
0.96590004
1.0083938
0.96662298
1.0221454
1.0069792
1.0343996
1.0066531
1.0072525
0.99743563
0.99723703
1.000372
0.99013917
1.0095223
0.98864268
0.98092242
0.98886488
1.0030341
1.01894
0.99155059
0.99533235