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

Account for measurement error in non-MCMC forecasts

time-shift
Johannes Pfeifer 2016-06-19 12:41:51 +02:00 committed by Stéphane Adjemian (Hermes)
parent 46450c6f1f
commit 0f72e6ed63
4 changed files with 129 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
matlab/dyn_forecast.m
View File

@ -127,7 +127,11 @@ switch task
end
if M.exo_det_nbr == 0
[yf,int_width] = forcst(oo.dr,y0,horizon,var_list,M,oo,options);
if isequal(M.H,0)
[yf,int_width] = forcst(oo.dr,y0,horizon,var_list,M,oo,options);
else
[yf,int_width,int_width_ME] = forcst(oo.dr,y0,horizon,var_list,M,oo,options);
end
else
exo_det_length = size(oo.exo_det_simul,1)-M.maximum_lag;
if horizon > exo_det_length
@ -139,8 +143,13 @@ else
elseif horizon < exo_det_length
ex = zeros(exo_det_length,M.exo_nbr);
end
[yf,int_width] = simultxdet(y0,ex,oo.exo_det_simul,...
if isequal(M.H,0)
[yf,int_width] = simultxdet(y0,ex,oo.exo_det_simul,...
options.order,var_list,M,oo,options);
else
[yf,int_width,int_width_ME] = simultxdet(y0,ex,oo.exo_det_simul,...
options.order,var_list,M,oo,options);
end
end
if ~isscalar(trend) %add trend back to forecast
@ -162,6 +171,10 @@ for i=1:n_var
forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))';
forecast.HPDinf.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width(1:horizon,i);
forecast.HPDsup.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width(1:horizon,i);
if ~isequal(M.H,0) && ismember(var_list(i,:),options.varobs)
forecast.HPDinf_ME.(vname)= yf(i,maximum_lag+(1:horizon))' - int_width_ME(1:horizon,i);
forecast.HPDsup_ME.(vname) = yf(i,maximum_lag+(1:horizon))' + int_width_ME(1:horizon,i);
end
end
for i=1:M.exo_det_nbr

18
matlab/forcst.m
View File

@ -1,5 +1,5 @@
function [yf,int_width]=forcst(dr,y0,horizon,var_list,M_,oo_,options_)
% function [yf,int_width]=forecst(dr,y0,horizon,var_list,M_,oo_,options_)
function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,options_)
% function [yf,int_width,int_width_ME]=forecst(dr,y0,horizon,var_list,M_,oo_,options_)
% computes mean forecast for a given value of the parameters
% computes also confidence band for the forecast
%
@ -16,6 +16,8 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list,M_,oo_,options_)
% yf: mean forecast
% int_width: distance between upper bound and
% mean forecast
% int_width_ME:distance between upper bound and
% mean forecast when considering measurement error
%
% SPECIAL REQUIREMENTS
% none
@ -81,12 +83,22 @@ for i=1:horizon
sigma_u = A*sigma_u*A';
sigma_y = sigma_y+sigma_u;
end
if nargout==3
var_yf_ME=var_yf;
[loc_H,loc_varlist]=ismember(options_.varobs',options_.varlist);
loc_varlist(loc_varlist==0)=[];
var_yf_ME(:,loc_varlist)=var_yf(:,loc_varlist)+repmat(diag(M_.H(loc_H,loc_H))',horizon,1);
int_width_ME = zeros(horizon,nvar);
end
fact = norminv((1-options_.forecasts.conf_sig)/2,0,1);
int_width = zeros(horizon,M_.endo_nbr);
int_width = zeros(horizon,nvar);
for i=1:nvar
int_width(:,i) = -fact*sqrt(var_yf(:,i));
if nargout==3
int_width_ME(:,i) = -fact*sqrt(var_yf_ME(:,i));
end
end
yf = yf(ivar,:);

91
matlab/forecast_graphs.m
View File

@ -1,5 +1,5 @@
function forecast_graphs(var_list,M_, oo_,options_)
% function forecast_graphs(var_list)
% function forecast_graphs(var_list,M_, oo_,options_)
% Plots the classical forecasts created by dyn_forecast.m
%
% Inputs:
@ -55,12 +55,8 @@ if ~exist([dname '/graphs'],'dir')
mkdir(dname,'graphs');
end
%% write LaTeX-Header
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([M_.dname '/graphs/',fname,'_forcst.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by Dynare''s graph_decomp.m.\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
fidTeX=write_LaTeX_header([M_.dname '/graphs/',fname,'_forcst.tex']);
end
m = 1;
@ -74,7 +70,7 @@ for j= 1:nvar
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s/graphs/forcst%d}\n',dname,n_fig);
fprintf(fidTeX,'\\label{Fig:forcst:%d}\n',n_fig);
fprintf(fidTeX,'\\caption{Mean forecasts and %2.0f%% confidence intervals}\n',options_.forecasts.conf_sig*100);
fprintf(fidTeX,'\\caption{Mean forecasts and %2.0f\\%% confidence intervals}\n',options_.forecasts.conf_sig*100);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
@ -116,9 +112,80 @@ if m > 1
end
end
%% write LaTeX-Footer
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,' \n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
write_LaTeX_footer(fidTeX);
end
if isfield(oo_.forecast,'HPDinf_ME')
var_names=fieldnames(oo_.forecast.HPDinf_ME);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX=write_LaTeX_header([M_.dname '/graphs/',fname,'_forcst_ME.tex']);
end
m = 1;
n_fig = 1;
hh=dyn_figure(options_,'Name','Forecasts including ME (I)');
for j= 1:length(var_names)
if m > nc*nr;
dyn_saveas(hh,[ dname '/graphs/forcst_ME' int2str(n_fig)],options_);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s/graphs/forcst_ME%d}\n',options_.figures.textwidth*min((m-1)/nc,1),dname,n_fig);
fprintf(fidTeX,'\\label{Fig:forcst_ME:%d}\n',n_fig);
fprintf(fidTeX,'\\caption{Mean forecasts and %2.0f\\%% confidence intervals accounting for measurement error}\n',options_.conf_sig*100);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
n_fig =n_fig+1;
eval(['hh=dyn_figure(options_,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
m = 1;
end
subplot(nr,nc,m);
vn = var_names{j};
obs = 0;
plot([NaN(obs,1); yf.(vn)],'b-');
hold on
plot([NaN(obs,1); oo_.forecast.HPDinf_ME.(vn)],'b-');
hold on
plot([NaN(obs,1); oo_.forecast.HPDsup_ME.(vn)],'b-');
title(vn,'Interpreter','none');
xlim([1 obs+length(oo_.forecast.HPDsup_ME.(vn))])
hold off
m = m + 1;
end
if m > 1
dyn_saveas(hh,[dname '/graphs/forcst_ME' int2str(n_fig)],options_);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s/graphs/forcst_ME%d}\n',options_.figures.textwidth*min((m-1)/nc,1),dname,n_fig);
fprintf(fidTeX,'\\label{Fig:forcst_ME:%d}\n',n_fig);
fprintf(fidTeX,'\\caption{Mean forecasts and %2.0f\\%% confidence intervals accounting for measurement error}\n',options_.conf_sig*100);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
write_LaTeX_footer(fidTeX);
end
end
function fidTeX=write_LaTeX_header(fname)
%% write LaTeX-Header
fidTeX = fopen(fname,'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by Dynare''s forecast_graphs.m.\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
function write_LaTeX_footer(fidTeX)
%% write LaTeX-Footer
fprintf(fidTeX,' \n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);

23
matlab/simultxdet.m
View File

@ -1,4 +1,4 @@
function [y_,int_width]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
function [y_,int_width,int_width_ME]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
%function [y_,int_width]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
%
% Simulates a stochastic model in the presence of deterministic exogenous shocks
@ -9,12 +9,20 @@ function [y_,int_width]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_
% ex_det: matrix of deterministic exogenous shocks, starting at period 1-M_.maximum_lag
% iorder: order of approximation
% var_list: list of endogenous variables to simulate
% int_width_ME:distance between upper bound and
% mean forecast when considering measurement error
% OUTPUTS:
% yf: mean forecast
% int_width: distance between upper bound and
% mean forecast
% int_width_ME:distance between upper bound and
% mean forecast when considering measurement error
%
% The forecast horizon is equal to size(ex, 1).
% The condition size(ex,1)+M_.maximum_lag=size(ex_det,1) must be verified
% for consistency.
% Copyright (C) 2008-2012 Dynare Team
% Copyright (C) 2008-2016 Dynare Team
%
% This file is part of Dynare.
%
@ -124,6 +132,7 @@ sigma_u = B*M_.Sigma_e*B';
sigma_u1 = ghu1*M_.Sigma_e*ghu1';
sigma_y = 0;
var_yf=NaN(horizon,nvar); %initialize
for i=1:iter
sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
var_yf(i,:) = diag(sigma_y1)';
@ -135,8 +144,16 @@ for i=1:iter
end
fact = norminv((1-options_.forecasts.conf_sig)/2,0,1);
if nargout==3
var_yf_ME=var_yf;
var_yf_ME(:,options_.varobs_id)=var_yf(:,options_.varobs_id)+repmat(diag(M_.H)',horizon,1);
int_width_ME = zeros(horizon,M_.endo_nbr);
end
int_width = zeros(iter,endo_nbr);
int_width = zeros(iter,nvar);
for i=1:nvar
int_width(:,i) = fact*sqrt(var_yf(:,i));
if nargout==3
int_width_ME(:,i) = -fact*sqrt(var_yf_ME(:,i));
end
end