From 0f72e6ed63ef226b52b5b87f8a5e2210979aab60 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 19 Jun 2016 12:41:51 +0200 Subject: [PATCH] Account for measurement error in non-MCMC forecasts --- matlab/dyn_forecast.m | 17 +++++++- matlab/forcst.m | 18 ++++++-- matlab/forecast_graphs.m | 91 ++++++++++++++++++++++++++++++++++------ matlab/simultxdet.m | 23 ++++++++-- 4 files changed, 129 insertions(+), 20 deletions(-) diff --git a/matlab/dyn_forecast.m b/matlab/dyn_forecast.m index bea652c3a..56f977318 100644 --- a/matlab/dyn_forecast.m +++ b/matlab/dyn_forecast.m @@ -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 diff --git a/matlab/forcst.m b/matlab/forcst.m index fec6201de..c2e46ddf3 100644 --- a/matlab/forcst.m +++ b/matlab/forcst.m @@ -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,:); diff --git a/matlab/forecast_graphs.m b/matlab/forecast_graphs.m index 9053b4843..0b6708eb6 100644 --- a/matlab/forecast_graphs.m +++ b/matlab/forecast_graphs.m @@ -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 \ No newline at end of file + 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); diff --git a/matlab/simultxdet.m b/matlab/simultxdet.m index 2de1bb4b7..913d98101 100644 --- a/matlab/simultxdet.m +++ b/matlab/simultxdet.m @@ -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