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

Bug corrections : TeX output, oo_ structure (posterior density), posterior IRF. 

Added : Posterior IRFs for BVAR-DSGE, conditional forecasts

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1119 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2006-12-15 11:37:24 +00:00
parent be18e12a43
commit 47ae595059
9 changed files with 865 additions and 430 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/DsgeLikelihood.m
View File

@ -192,4 +192,4 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
fval = (likelihood-lnprior);

332
matlab/GetPosteriorParametersStatistics.m
View File

@ -1,5 +1,22 @@
function GetPosteriorParametersStatistics()
% stephane.adjemian@ens.fr [09-09-2005]
function get_posterior_parameters_statistics()
% This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output).
%
% INPUTS
% None.
%
% OUTPUTS
% None.
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global estim_params_ M_ options_ bayestopt_ oo_
TeX = options_.TeX;
@ -12,6 +29,8 @@ np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
DirectoryName = CheckPath('metropolis');
OutputDirectoryName = CheckPath('Output');
load([ DirectoryName '/' M_.fname '_mh_history'])
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
@ -21,138 +40,199 @@ FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
clear record;
disp(' ')
disp(' ')
disp('ESTIMATION RESULTS')
disp(' ')
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
pnames=[' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
tit2 = sprintf('%10s %7s %10s %14s %4s %6s\n',' ','prior mean','post. mean','conf. interval','prior','pstdev');
pformat = '%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f';
disp(' ');disp(' ');disp('ESTIMATION RESULTS');disp(' ');
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
if np
disp(' ')
disp('parameters')
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
name = bayestopt_.name{ip};
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name, ...
bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
eval(['oo_.posterior_mean.parameters.' name ' = post_mean;']);
eval(['oo_.posterior_hpdinf.parameters.' name ' = hpd_interval(1);']);
eval(['oo_.posterior_hpdsup.parameters.' name ' = hpd_interval(2);']);
eval(['oo_.posterior_median.' name ' = post_median;']);
eval(['oo_.posterior_variance.' name ' = post_var;']);
eval(['oo_.posterior_deciles.' name ' = post_deciles;']);
eval(['oo_.posterior_density.' name ' = density;']);
ip = ip+1;
end
if TeX
fid = TeXBegin(OutputDirectoryName,M_.fname,1,'parameters');
end
disp(' ')
disp('parameters')
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
name = bayestopt_.name{ip};
disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
oo_ = Filloo(oo_,name,'parameters',post_mean,hpd_interval,post_median,post_var,post_deciles,density);
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,1,'parameters');
end
end
if nvx
ip = 1;
disp(' ')
disp('standard deviation of shocks')
disp(tit2)
for i=1:nvx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name,bayestopt_.pmean(ip),post_mean, ...
hpd_interval,pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
M_.Sigma_e(k,k) = post_mean*post_mean;
eval(['oo_.posterior_mean.shocks_std.' name ' = post_mean;']);
eval(['oo_.posterior_hpdinf.shocks_std.' name ' = hpd_interval(1);']);
eval(['oo_.posterior_hpdsup.shocks_std.' name ' = hpd_interval(2);']);
eval(['oo_.posterior_median.shocks_std.' name ' = post_median;']);
eval(['oo_.posterior_variance.shocks_std.' name ' = post_var;']);
eval(['oo_.posterior_deciles.shocks_std.' name ' = post_deciles;']);
eval(['oo_.posterior_density.shocks_std.' name ' = density;']);
ip = ip+1;
end
if TeX
fid = TeXBegin(OutputDirectoryName,M_.fname,2,'standard deviation of structural shocks');
end
ip = 1;
disp(' ')
disp('standard deviation of shocks')
disp(tit2)
for i=1:nvx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval,...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
M_.Sigma_e(k,k) = post_mean*post_mean;
oo_ = Filloo(oo_,name,'shocks_std',post_mean,hpd_interval,post_median,post_var,post_deciles,density);
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,2,'standard deviation of structural shocks');
end
end
if nvn
disp(' ')
disp('standard deviation of measurement errors')
disp(tit2)
ip = nvx+1;
for i=1:nvn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name,...
bayestopt_.pmean(ip), ...
post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
eval(['oo_.posterior_mean.measurement_errors_std.' name ' = post_mean;']);
eval(['oo_.posterior_hpdinf.measurement_errors_std.' name ' = hpd_interval(1);']);
eval(['oo_.posterior_hpdsup.measurement_errors_std.' name ' = hpd_interval(2);']);
eval(['oo_.posterior_median.measurement_errors_std.' name ' = post_median;']);
eval(['oo_.posterior_variance.measurement_errors_std.' name ' = post_var;']);
eval(['oo_.posterior_deciles.measurement_errors_std.' name ' = post_deciles;']);
eval(['oo_.posterior_density.measurement_errors_std.' name ' = density;']);
ip = ip+1;
end
if TeX
fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors')
end
disp(' ')
disp('standard deviation of measurement errors')
disp(tit2)
ip = nvx+1;
for i=1:nvn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
oo_ = Filloo(oo_,name,'measurement_errors_std',post_mean,hpd_interval,post_median,...
post_var,post_deciles,density);
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,3,'standard deviation of measurement errors');
end
end
if ncx
disp(' ')
disp('correlation of shocks')
disp(tit2)
ip = nvx+nvn+1;
for i=1:ncx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
eval(['oo_.posterior_mean.shocks_corr.' NAME ' = post_mean;']);
eval(['oo_.posterior_hpdinf.shocks_corr.' NAME ' = hpd_interval(1);']);
eval(['oo_.posterior_hpdsup.shocks_corr.' NAME ' = hpd_interval(2);']);
eval(['oo_.posterior_median.shocks_corr.' NAME ' = post_median;']);
eval(['oo_.posterior_variance.shocks_corr.' NAME ' = post_var;']);
eval(['oo_.posterior_deciles.shocks_corr.' NAME ' = post_deciles;']);
eval(['oo_.posterior_density.shocks_corr.' NAME ' = density;']);
M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
ip = ip+1;
end
if TeX
fid = TeXBegin(OutputDirectoryName,M_.fname,4,'correlation of structural shocks');
end
disp(' ')
disp('correlation of shocks')
disp(tit2)
ip = nvx+nvn+1;
for i=1:ncx
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
oo_ = Filloo(oo_,NAME,'shocks_corr',post_mean,hpd_interval,post_median,post_var,post_deciles,density);
M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,4,'correlation of structural shocks');
end
end
if ncn
disp(' ')
disp('correlation of measurement errors')
disp(tit2)
ip = nvx+nvn+ncx+1;
for i=1:ncn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
eval(['oo_.posterior_mean.measurement_errors_corr.' NAME ' = post_mean;']);
eval(['oo_.posterior_hpdinf.measurement_errors_corr.' NAME ' = hpd_interval(1);']);
eval(['oo_.posterior_hpdsup.measurement_errors_corr.' NAME ' = hpd_interval(2);']);
eval(['oo_.posterior_median.measurement_errors_corr.' NAME ' = post_median;']);
eval(['oo_.posterior_variance.measurement_errors_corr.' NAME ' = post_var;']);
eval(['oo_.posterior_deciles.measurement_errors_corr.' NAME ' = post_decile;']);
eval(['oo_.posterior_density.measurement_errors_corr.' NAME ' = density;']);
ip = ip+1;
end
end
if TeX
fid = TeXBegin(OutputDirectoryName,M_.fname,5,'correlation of measurement errors');
end
disp(' ')
disp('correlation of measurement errors')
disp(tit2)
ip = nvx+nvn+ncx+1;
for i=1:ncn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
oo_ = Filloo(oo_,NAME,'measurement_errors_corr',post_mean,hpd_interval,...
post_median,post_var,post_deciles,density);
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
if TeX
TeXEnd(fid,5,'correlation of measurement errors');
end
end
%
%% subfunctions:
%
function fid = TeXBegin(directory,fname,fnum,title)
TeXfile = [directory '/' fname '_Posterior_Mean_' int2str(fnum) '.TeX'];
fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by Dynare.\n');
fprintf(fidTeX,['%% RESULTS FROM METROPOLIS HASTINGS (' title ')\n']);
fprintf(fidTeX,['%% ' datestr(now,0)]);
fprintf(fidTeX,' \n');
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{table}\n');
fprintf(fidTeX,'\\centering\n');
fprintf(fidTeX,'\\begin{tabular}{l|lcccccc} \n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,[' & Prior distribution & Prior mean & Prior ' ...
's.d. & Posterior mean & Posterior s.d. & HPD inf & HPD sup\\\\ \n']);
fprintf(fidTeX,'\\hline \\\\ \n');
fid = fidTeX;
function TeXCore(fid,name,shape,priormean,priorstd,postmean,poststd,hpd)
fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n'],...
name,...
shape,...
priormean,...
priorstd,...
postmean,...
poststd,...
hpd(1),...
hpd(2));
function TeXEnd(fid,fnum,title)
fprintf(fid,'\\hline\\hline \n');
fprintf(fid,'\\end{tabular}\n ');
fprintf(fid,['\\caption{Results from Metropolis-Hastings (' title ')}\n ']);
fprintf(fid,['\\label{Table:MHPosterior:' int2str(fnum) '}\n']);
fprintf(fid,'\\end{table}\n');
fprintf(fid,'%% End of TeX file.\n');
fclose(fid);
function oo = Filloo(oo,name,type,postmean,hpdinterval,postmedian,postvar,postdecile,density)
eval(['oo.posterior_mean.' type '.' name ' = postmean;']);
eval(['oo.posterior_hpdinf.' type '.' name ' = hpdinterval(1);']);
eval(['oo.posterior_hpdsup.' type '.' name ' = hpdinterval(2);']);
eval(['oo.posterior_median.' type '.' name ' = postmedian;']);
eval(['oo.posterior_variance.' type '.' name ' = postvar;']);
eval(['oo.posterior_deciles.' type '.' name ' = postdecile;']);
eval(['oo.posterior_density.' type '.' name ' = density;']);

4
matlab/PlotPosteriorDistributions.m
View File

@ -88,8 +88,8 @@ for i=1:npar
else
j = i - (nvx+nvn+ncx+ncn);
name = deblank(M_.param_names(estim_params_.param_vals(j,1),:));
eval(['x1 = oo_.posterior_density.' name '(:,1);'])
eval(['f1 = oo_.posterior_density.' name '(:,2);'])
eval(['x1 = oo_.posterior_density.parameters.' name '(:,1);'])
eval(['f1 = oo_.posterior_density.parameters.' name '(:,2);'])
if options_.posterior_mode_estimation
eval(['pmode = oo_.posterior_mode.parameters.' name ';'])
end

747
matlab/PosteriorIRF.m
View File

@ -1,281 +1,466 @@
function posterior_irf(type)
% Metropolis-Hastings algorithm.
%
% INPUTS
% o type [char] string specifying the joint density of the
% deep parameters ('prior','posterior').
%
% OUTPUTS
% None (oo_ and plots).
%
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global options_ estim_params_ oo_ M_ dsge_prior_weight
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
%
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
DirectoryName = CheckPath('Output');
if strcmpi(type,'posterior')
MhDirectoryName = CheckPath('metropolis');
elseif strcmpi(type,'gsa')
MhDirectoryName = CheckPath('GSA');
else
MhDirectoryName = CheckPath('prior');
end
MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
if strcmpi(type,'posterior')
load([ MhDirectoryName '/' M_.fname '_mh_history'])
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
elseif strcmpi(type,'gsa')
load([ MhDirectoryName '/' M_.fname '_prior'],'lpmat','istable')
x=lpmat(istable,:);
clear lpmat istable
NumberOfDraws=size(x,1);
B=NumberOfDraws; options_.B = B;
else% type = 'prior'
NumberOfDraws = 500;
end
if ~strcmpi(type,'gsa')
B = min([round(.5*NumberOfDraws),500]); options_.B = B;
end
try delete([MhDirectoryName '\' M_.fname '_IRFs*']);
catch disp('No _IRFs files to be deleted!')
end
irun = 0;
irun2 = 0;
NumberOfIRFfiles = 1;
ifil2 = 1;
if strcmpi(type,'posterior')
h = waitbar(0,'Bayesian (posterior) IRFs...');
elseif strcmpi(type,'gsa')
h = waitbar(0,'GSA (prior) IRFs...');
else
h = waitbar(0,'Bayesian (prior) IRFs...');
end
if B <= MAX_nruns
stock_param = zeros(B, npar);
else
stock_param = zeros(MAX_nruns, npar);
end
if B >= MAX_nirfs
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
else
stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
end
for b=1:B
irun = irun+1;
irun2 = irun2+1;
if ~strcmpi(type,'gsa')
deep = GetOneDraw(type);
else
deep = x(b,:);
end
stock_param(irun2,:) = deep;
set_parameters(deep);
dr = resol(oo_.steady_state,0);
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
cs = transpose(chol(SS));
for i = 1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
end
for j = 1:M_.endo_nbr
if max(y(j,:)) - min(y(j,:)) > 1e-10
stock_irf(:,j,i,irun) = transpose(y(j,:));
end
end
end
end
if irun == MAX_nirfs | irun == B | b == B
if b == B
stock_irf = stock_irf(:,:,:,1:irun);
end
save([MhDirectoryName '/' M_.fname '_irf' int2str(NumberOfIRFfiles)],'stock_irf');
NumberOfIRFfiles = NumberOfIRFfiles+1;
irun = 0;
end
if irun2 == MAX_nruns | b == B
if b == B
stock_param = stock_param(1:irun2,:);
end
stock = stock_param;
save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2)],'stock');
ifil2 = ifil2 + 1;
irun2 = 0;
end
waitbar(b/B,h);
end
NumberOfIRFfiles = NumberOfIRFfiles-1;
ifil2 = ifil2-1;
close(h);
ReshapeMatFiles('irf',type)
if strcmpi(type,'gsa')
return
end
varlist = options_.varlist;
if isempty(varlist)
varlist = M_.endo_names;
SelecVariables = transpose(1:M_.endo_nbr);
nvar = M_.endo_nbr;
else
nvar = size(varlist,1);
SelecVariables = [];
for i=1:nvar
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
end
end
end
MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
StdIRF = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
if options_.TeX
varlist_TeX = [];
for i=1:nvar
varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(SelecVariables(i),:));
end
end
fprintf('MH: Posterior IRFs...\n');
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles
load([MhDirectoryName '/' M_.fname '_IRFs' int2str(file)]);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),DistribIRF(kk,:,j,i)] = ...
posterior_moments(squeeze(STOCK_IRF(k,SelecVariables(j),i,:)),0);
end
end
end
kdx = kdx + size(STOCK_IRF,1);
end
clear STOCK_IRF;
for i = 1:M_.exo_nbr
for j = 1:nvar
name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.Distribution.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
end
end
%%
%% Finally i build the plots.
%%
if options_.TeX
fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
end
%%
subplotnum = 0;
for i=1:M_.exo_nbr
NAMES = [];
if options_.TeX
TEXNAMES = [];
end
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
subplotnum = subplotnum+1;
if options_.nograph
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
end
else
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
end
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold on
for k = 1:9
plot(1:options_.irf,DistribIRF(:,k,j,i),'-g','linewidth',0.5)
end
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',1)
xlim([1 options_.irf]);
hold off
name = deblank(varlist(j,:));
NAMES = strvcat(NAMES,name);
if options_.TeX
texname = deblank(varlist_TeX(j,:));
TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
end
title(name,'Interpreter','none')
end
if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum>0)
figunumber = figunumber+1;
set(hh,'visible','on')
eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']);
set(hh,'visible','off')
if options_.nograph, close(hh), end
if options_.TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
if options_.relative_irf
fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
else
fprintf(fidTeX,'\\caption{Bayesian IRF.}');
end
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
subplotnum = 0;
end
end% loop over selected endo_var
end% loop over exo_var
%%
if options_.TeX
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
fprintf('MH: Posterior IRFs, done!\n');
function posterior_irf(type)
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
% o type [char] string specifying the joint density of the
% deep parameters ('prior','posterior').
%
% OUTPUTS
% None (oo_ and plots).
%
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global options_ estim_params_ oo_ M_ bayestopt_
if isempty(options_.irf) | ~options_.irf
options_.irf = 40;
end
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
nvobs = size(options_.varobs,1);
gend = options_.nobs;
%
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
DirectoryName = CheckPath('Output');
if strcmpi(type,'posterior')
MhDirectoryName = CheckPath('metropolis');
elseif strcmpi(type,'gsa')
MhDirectoryName = CheckPath('GSA');
else
MhDirectoryName = CheckPath('prior');
end
MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8)+50;
else
MAX_nirfs_dsgevar = 0;
end
if strcmpi(type,'posterior')
load([ MhDirectoryName '/' M_.fname '_mh_history'])
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
elseif strcmpi(type,'gsa')
load([ MhDirectoryName '/' M_.fname '_prior'],'lpmat','istable')
x=lpmat(istable,:);
clear lpmat istable
NumberOfDraws=size(x,1);
B=NumberOfDraws; options_.B = B;
else% type = 'prior'
NumberOfDraws = 500;
end
if ~strcmpi(type,'gsa')
B = min([round(.5*NumberOfDraws),500]); options_.B = B;
end
try delete([MhDirectoryName '/' M_.fname '_irf_dsge*.mat'])
catch disp('No _IRFs (dsge) files to be deleted!')
end
try delete([MhDirectoryName '/' M_.fname '_irf_bvardsge*.mat'])
catch disp('No _IRFs (bvar-dsge) files to be deleted!')
end
irun = 0;
IRUN = 0;
irun2 = 0;
NumberOfIRFfiles_dsge = 1;
NumberOfIRFfiles_dsgevar = 1;
ifil2 = 1;
if strcmpi(type,'posterior')
h = waitbar(0,'Bayesian (posterior) IRFs...');
elseif strcmpi(type,'gsa')
h = waitbar(0,'GSA (prior) IRFs...');
else
h = waitbar(0,'Bayesian (prior) IRFs...');
end
if B <= MAX_nruns
stock_param = zeros(B, npar);
else
stock_param = zeros(MAX_nruns, npar);
end
if B >= MAX_nirfs_dsge
stock_irf_dsge = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs_dsge);
else
stock_irf_dsge = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
end
if MAX_nirfs_dsgevar
if B >= MAX_nirfs_dsgevar
stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
else
stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B);
end
[mYY,mXY,mYX,mXX,Ydata,Xdata] = ...
var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);
NumberOfLags = options_.varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
COMP_draw = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
for b=1:B
irun = irun+1;
irun2 = irun2+1;
if ~strcmpi(type,'gsa')
deep = GetOneDraw(type);
else
deep = x(b,:);
end
stock_param(irun2,:) = deep;
set_parameters(deep);
dr = resol(oo_.steady_state,0);
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
SS = transpose(chol(SS));
for i = 1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
end
for j = 1:M_.endo_nbr
if max(y(j,:)) - min(y(j,:)) > 1e-10
stock_irf_dsge(:,j,i,irun) = transpose(y(j,:));
end
end
end
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
[fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
tmp1 = SIGMAu*gend*(dsge_prior_weight+1);
val = 1;
tmp1 = chol(inv(tmp1))';
while val;
% draw from the marginal posterior of sig
tmp2 = tmp1*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
SIGMAu_draw = inv(tmp2*tmp2');
% draw from the conditional posterior of PHI
VARvecPHI = kron(SIGMAu_draw,iXX);
PHI_draw = PHI(:) + chol(VARvecPHI)'*randn(nvobs*NumberOfLagsTimesNvobs,1);
COMP_draw(1:nvobs,:) = reshape(PHI_draw,NumberOfLagsTimesNvobs,nvobs)';
% Check for stationarity
tests = find(abs(eig(COMP_draw))>0.9999999999);
if isempty(tests)
val=0;
end
end
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
A0 = Atheta(bayestopt_.mfys,:);
[OMEGAstar,SIGMAtr] = qr2(A0');
end
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = COMP_draw*PHIpower;
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)';
end
for j = 1:(nvobs*M_.exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
end
if irun == MAX_nirfs_dsge | irun == B | b == B
if b == B
stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
if MAX_nirfs_dsgevar & (b == B | IRUN == B)
stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
irun = 0;
end
end
save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge)],'stock_irf_dsge');
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
irun = 0;
end
if irun2 == MAX_nruns | b == B
if b == B
stock_param = stock_param(1:irun2,:);
end
stock = stock_param;
save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2)],'stock');
ifil2 = ifil2 + 1;
irun2 = 0;
end
waitbar(b/B,h);
end
NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge-1;
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar-1;
ifil2 = ifil2-1;
close(h);
ReshapeMatFiles('irf_dsge')
if MAX_nirfs_dsgevar
ReshapeMatFiles('irf_bvardsge')
end
if strcmpi(type,'gsa')
return
end
varlist = options_.varlist;
if isempty(varlist)
if MAX_nirfs_dsgevar
varlist = options_.varobs;
nvar = size(varlist,1);
SelecVariables = [];
for i=1:nvar
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
end
end
else
varlist = M_.endo_names;
SelecVariables = transpose(1:M_.endo_nbr);
nvar = M_.endo_nbr;
end
else
nvar = size(varlist,1);
SelecVariables = [];
for i=1:nvar
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
end
end
if MAX_nirfs_dsgevar% Here I test if declared variables are observed variables.
SelecVariables = [];
varlistbis = [];
for i=1:nvar
if ~isempty(strmatch(varlist(i,:),options_.varobs,'exact'))
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
varlistbis = strvcat(varlistbis,varlist(i,:));
else
disp(' ')
disp(['Warning :: ' varlist(i,:) 'is not an observed variable!'])
disp(['This variable will not be considered for the IRFs.'])
disp(' ')
end
end
varlist = varlistbis;
end
end
MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
VarIRF = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
if options_.TeX
varlist_TeX = [];
for i=1:nvar
varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(SelecVariables(i),:));
end
end
fprintf('MH: Posterior (dsge) IRFs...\n');
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsge
load([MhDirectoryName '/' M_.fname '_IRF_DSGEs' int2str(file)]);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF_DSGE,1)
kk = k+kdx;
[MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),DistribIRF(kk,:,j,i)] = ...
posterior_moments(squeeze(STOCK_IRF_DSGE(k,SelecVariables(j),i,:)),0);
end
end
end
kdx = kdx + size(STOCK_IRF_DSGE,1);
end
clear STOCK_IRF_DSGE;
for i = 1:M_.exo_nbr
for j = 1:nvar
name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
end
end
if MAX_nirfs_dsgevar
MeanIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
MedianIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
fprintf('MH: Posterior (bvar-dsge) IRFs...\n');
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
for file = 1:NumberOfIRFfiles_dsgevar
load([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs' int2str(file)]);
for i = 1:M_.exo_nbr
for j = 1:nvar
for k = 1:size(STOCK_IRF_BVARDSGE,1)
kk = k+kdx;
[MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0);%SelecVariables(j)
end
end
end
kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
end
clear STOCK_IRF_BVARDSGE;
for i = 1:M_.exo_nbr
for j = 1:nvar
name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
eval(['oo_.PosteriorIRF.bvardsge.Mean.' name ' = MeanIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.Median.' name ' = MedianIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.Var.' name ' = VarIRFdsgevar(:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.Distribution.' name ' = DistribIRFdsgevar(:,:,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.HPDinf.' name ' = HPDIRFdsgevar(:,1,j,i);']);
eval(['oo_.PosteriorIRF.bvardsge.HPDsup.' name ' = HPDIRFdsgevar(:,2,j,i);']);
end
end
end
%%
%% Finally I build the plots.
%%
if options_.TeX
fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
end
%%
subplotnum = 0;
for i=1:M_.exo_nbr
NAMES = [];
if options_.TeX
TEXNAMES = [];
end
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
subplotnum = subplotnum+1;
if options_.nograph
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
end
else
if subplotnum == 1 & options_.relative_irf
hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
elseif subplotnum == 1 & ~options_.relative_irf
hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
end
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
if ~MAX_nirfs_dsgevar
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min(HPDIRF(:,1,j,i)));
hold on
h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i)));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min(HPDIRF(:,1,j,i)));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
hold off
else
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
hold on;
h2 = area(1:options_.irf,HPDIRF(:,1,j,i));
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
plot(1:options_.irf,MeanIRFdsgevar(:,j,i),'--k','linewidth',2)
plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1)
plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1)
box on
axis tight
xlim([1 options_.irf]);
hold off
end
name = deblank(varlist(j,:));
NAMES = strvcat(NAMES,name);
if options_.TeX
texname = deblank(varlist_TeX(j,:));
TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
end
title(name,'Interpreter','none')
end
if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum> 0)
figunumber = figunumber+1;
set(hh,'visible','on')
eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']);
set(hh,'visible','off')
if options_.nograph, close(hh), end
if options_.TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
if options_.relative_irf
fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
else
fprintf(fidTeX,'\\caption{Bayesian IRF.}');
end
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
subplotnum = 0;
end
end% loop over selected endo_var
end% loop over exo_var
%%
if options_.TeX
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
fprintf('MH: Posterior IRFs, done!\n');

26
matlab/ReshapeMatFiles.m
View File

@ -26,10 +26,14 @@ else
end
end
switch type
case 'irf'
CAPtype = 'IRF';
case 'irf_dsge'
CAPtype = 'IRF_DSGE';
TYPEsize = [ options_.irf , M_.endo_nbr , M_.exo_nbr ];
TYPEarray = 4;
TYPEarray = 4;
case 'irf_bvardsge'
CAPtype = 'IRF_BVARDSGE';
TYPEsize = [ options_.irf , size(options_.varobs,1) , M_.exo_nbr ];
TYPEarray = 4;
case 'smooth'
CAPtype = 'SMOOTH';
TYPEsize = [ M_.endo_nbr , options_.nobs ];
@ -59,7 +63,7 @@ end
return
end
TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*']);
TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*.mat']);
NumberOfTYPEfiles = length(TYPEfiles);
B = options_.B;
@ -98,9 +102,10 @@ end
%eval(['STOCK_' CAPtype ' = sort(stock_' type ',4);'])
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1)],['STOCK_' CAPtype ]);
end
for file = 1:NumberOfTYPEfiles
delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
end
% Original file format may be useful in some cases...
% for file = 1:NumberOfTYPEfiles
% delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
% end
case 3
if NumberOfTYPEfiles>1
NumberOfPeriodsPerTYPEfiles = ceil( TYPEsize(2)/NumberOfTYPEfiles );
@ -132,7 +137,8 @@ end
%eval(['STOCK_' CAPtype ' = sort(stock_' type ',3);'])
save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1)],['STOCK_' CAPtype ]);
end
for file = 1:NumberOfTYPEfiles
delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
end
% Original file format may be useful in some cases...
% for file = 1:NumberOfTYPEfiles
% delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
% end
end

20
matlab/dynare_estimation.m
View File

@ -650,7 +650,7 @@ OutputDirectoryName = CheckPath('Output');
if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior mode) Latex output
if np
filename = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_1.TeX'];
filename = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_1.TeX'];
fidTeX = fopen(filename,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (parameters)\n');
@ -670,7 +670,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
M_.param_names_tex(estim_params_.param_vals(i,1),:),...
deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
bayestopt_.pmean(ip),...
estim_params_.param_vals(i,6),...
bayestopt_.pstdev(ip),...
xparam1(ip),...
stdh(ip));
ip = ip + 1;
@ -685,7 +685,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX);
end
if nvx
TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_2.TeX'];
TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_2.TeX'];
fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of structural shocks)\n');
@ -706,7 +706,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
deblank(M_.exo_names_tex(k,:)),...
deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
bayestopt_.pmean(ip),...
estim_params_.var_exo(i,7),...
bayestopt_.pstdev(ip),...
xparam1(ip), ...
stdh(ip));
ip = ip+1;
@ -721,7 +721,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX);
end
if nvn
TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_3.TeX'];
TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_3.TeX'];
fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of measurement errors)\n');
@ -741,7 +741,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
deblank(M_.endo_names_tex(idx,:)), ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
bayestopt_.pmean(ip), ...
estim_params_.var_endo(i,7),...
bayestopt_.pstdev(ip),...
xparam1(ip),...
stdh(ip));
ip = ip+1;
@ -755,7 +755,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX);
end
if ncx
TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_4.TeX'];
TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_4.TeX'];
fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of structural shocks)\n');
@ -776,7 +776,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
[deblank(M_.exo_names_tex(k1,:)) ',' deblank(M_.exo_names_tex(k2,:))], ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
bayestopt_.pmean(ip), ...
estim_params_.corrx(i,8), ...
bayestopt_.pstdev(ip), ...
xparam1(ip), ...
stdh(ip));
ip = ip+1;
@ -790,7 +790,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fclose(fidTeX);
end
if ncn
TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_5.TeX'];
TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_5.TeX'];
fidTeX = fopen(TeXfile,'w');
fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of measurement errors)\n');
@ -811,7 +811,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
[deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))], ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pmean(ip), ...
estim_params_.corrn(i,8), ...
bayestopt_.pstdev(ip), ...
xparam1(ip), ...
stdh(ip));
ip = ip+1;

119
matlab/imcforecast.m Executable file
View File

@ -0,0 +1,119 @@
function icforecast(ptype,cV,cS,cL,H,mcValue,B,ci)
% stephane.adjemian@ens.fr
global options_ oo_ M_
xparam = get_posterior_parameters(ptype);
gend = options_.nobs;
% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1 & ~options_.logdata
rawdata = log(rawdata);
end
if options_.prefilter == 1
bayestopt_.mean_varobs = mean(rawdata,1);
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else
data = transpose(rawdata);
end
set_parameters(xparam);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data);
InitState(:,1) = atT(:,end);
[T,R,ys,info] = dynare_resolve;
sQ = sqrt(M_.Sigma_e);
NumberOfStates = length(InitState);
FORCS1 = zeros(NumberOfStates,H+1,B);
for b=1:B
FORCS1(:,1,b) = InitState;
end
EndoSize = size(M_.endo_names,1);
ExoSize = size(M_.exo_names,1);
n1 = size(cV,1);
n2 = size(cS,1);
if n1 ~= n2
disp('imcforecast :: Error!')
disp(['imcforecast :: The number of variables doesn''t match the number of shocks'])
return
end
idx = [];
jdx = [];
for i = 1:n1
idx = [idx ; oo_.dr.inv_order_var(strmatch(deblank(cV(i,:)),M_.endo_names,'exact'))];
jdx = [jdx ; strmatch(deblank(cS(i,:)),M_.exo_names,'exact')];
end
mv = zeros(n1,NumberOfStates);
mu = zeros(ExoSize,n2);
for i=1:n1
mv(i,idx(i)) = 1;
mu(jdx(i),i) = 1;
end
if (size(mcValue,2) == 1);
mcValue = mcValue*ones(1,cL);
else
cL = size(mcValue,2);
end
randn('state',0);
for b=1:B
shocks = sQ*randn(ExoSize,H);
shocks(jdx,:) = zeros(length(jdx),H);
FORCS1(:,:,b) = mcforecast3(cL,H,mcValue,shocks,FORCS1(:,:,b),T,R,mv, mu);
end
mFORCS1 = mean(FORCS1,3);
tt = (1-ci)/2;
t1 = round(B*tt);
t2 = round(B*(1-tt));
forecasts.controled_variables = cV;
forecasts.instruments = cS;
for i = 1:EndoSize
eval(['forecasts.cond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS1(i,:)'';']);
tmp = sort(squeeze(FORCS1(i,:,:))');
eval(['forecasts.cond.ci.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ...
' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
end
clear FORCS1;
FORCS2 = zeros(NumberOfStates,H+1,B);
for b=1:B
FORCS2(:,1,b) = InitState;
end
randn('state',0);
for b=1:B
shocks = sQ*randn(ExoSize,H);
shocks(jdx,:) = zeros(length(jdx),H);
FORCS2(:,:,b) = mcforecast3(0,H,mcValue,shocks,FORCS2(:,:,b),T,R,mv, mu);
end
mFORCS2 = mean(FORCS2,3);
for i = 1:EndoSize
eval(['forecasts.uncond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS2(i,:)'';']);
tmp = sort(squeeze(FORCS2(i,:,:))');
eval(['forecasts.uncond.ci.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ...
' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
end
save('conditional_forecasts.mat','forecasts');

13
matlab/mcforecast3.m Executable file
View File

@ -0,0 +1,13 @@
function forcs = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
% stephane.adjemian@ens.fr [06-11-2006]
if cL
e = zeros(size(mcValue,1),cL);
for t = 1:cL
e(:,t) = inv(mv*R*mu)*(mcValue(:,t)-mv*T*forcs(:,t)-mv*R*shocks(:,t));
forcs(:,t+1) = T*forcs(:,t)+R*(mu*e(:,t)+shocks(:,t));
end
end
for t = cL+1:H
forcs(:,t+1) = T*forcs(:,t)+R*shocks(:,t);
end

32
matlab/plot_icforecast.m Executable file
View File

@ -0,0 +1,32 @@
function plot_icforecast(Variable)
% stephane.adjemian@ens.fr
load conditional_forecasts;
eval(['ci1 = forecasts.cond.ci.' Variable ';'])
eval(['m1 = forecasts.cond.mean.' Variable ';'])
eval(['ci2 = forecasts.uncond.ci.' Variable ';'])
eval(['m2 = forecasts.uncond.mean.' Variable ';'])
H = length(m1);
% area(1:H,ci1(2,:),'FaceColor',[.9 .9 .9],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
h1 = area(1:H,ci1(2,1:H))
set(h1,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
set(h1,'FaceColor',[.9 .9 .9])
hold on
% area(1:H,ci1(1,:),'FaceColor',[1 1 1],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
h2 = area(1:H,ci1(1,1:H));
set(h2,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
set(h2,'FaceColor',[1 1 1])
plot(1:H,m1,'-k','linewidth',3)
plot(1:H,m2,'--k','linewidth',3)
plot(1:H,ci2(1,:),'--k','linewidth',1)
plot(1:H,ci2(2,:),'--k','linewidth',1)
axis tight
hold off