function PosteriorIRF(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).
%
% SPECIAL REQUIREMENTS
% None
% PARALLEL CONTEXT
% This funtion has been parallelized in two different points. Then we have two core
% functions associated with it(the _core1 and _core2).
% See also the comments posterior_sampler.m funtion.
% Copyright (C) 2006-2018 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 .
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
% Set the number of periods
if isempty(options_.irf) || ~options_.irf
options_.irf = 40;
end
% Set varlist if necessary
varlist = options_.varlist;
if isempty(varlist)
varlist = options_.varobs;
end
options_.varlist = varlist;
nvar = length(varlist);
IndxVariables = [];
for i=1:nvar
idx = strmatch(varlist{i}, M_.endo_names, 'exact');
if isempty(idx)
disp(['PosteriorIRF :: ' varlist{i} 'is not a declared endogenous variable!'])
else
IndxVariables = [IndxVariables, idx];
end
end
% Get index of shocks for requested IRFs
irf_shocks_indx = getIrfShocksIndx();
% Set various parameters & Check or create directories
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; clear('nvx','nvn','ncx','ncn','np');
nvobs = dataset_.vobs;
gend = dataset_.nobs;
MaxNumberOfPlotPerFigure = 9;
nn = sqrt(MaxNumberOfPlotPerFigure);
MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
if options_.dsge_var
MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8);
else
MAX_nirfs_dsgevar = 0;
end
DirectoryName = CheckPath('Output',M_.dname);
if strcmpi(type,'posterior')
MhDirectoryName = CheckPath('metropolis',M_.dname);
elseif strcmpi(type,'gsa')
if options_.opt_gsa.pprior
MhDirectoryName = CheckPath(['GSA' filesep 'prior'],M_.dname);
else
MhDirectoryName = CheckPath(['GSA' filesep 'mc'],M_.dname);
end
else
MhDirectoryName = CheckPath('prior',M_.dname);
end
%delete old stale files before creating new ones
delete_stale_file([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
delete_stale_file([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
delete_stale_file([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']);
delete_stale_file([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat']);
delete_stale_file([MhDirectoryName filesep M_.fname '_param_irf*.mat']);
if strcmpi(type,'posterior')
B = options_.sub_draws;
options_.B = B;
if round((1-options_.mh_conf_sig)*B)<2
fprintf('\nPosteriorIRF:: options_.mh_conf_sig times options_.sub_draws is too small to generate HPDIs. I am omitting them.\n')
end
elseif strcmpi(type,'gsa')
RootDirectoryName = CheckPath('gsa',M_.dname);
if options_.opt_gsa.pprior
load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
else
load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
end
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat istable
B=size(x,1); options_.B = B;
else% type = 'prior'
B = options_.prior_draws;
options_.B = B;
end
irun = 0;
IRUN = 0;
irun2 = 0;
NumberOfIRFfiles_dsge = 1;
NumberOfIRFfiles_dsgevar = 1;
ifil2 = 1;
% Create arrays
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,nvar,M_.exo_nbr,MAX_nirfs_dsge);
else
stock_irf_dsge = zeros(options_.irf,nvar,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
NumberOfLags = options_.dsge_varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
if options_.noconstant
NumberOfParametersPerEquation = NumberOfLagsTimesNvobs;
else
NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1;
end
Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m
% function.
b = 0;
localVars=[];
% Save the local variables.
localVars.IRUN = IRUN;
localVars.irun = irun;
localVars.irun2=irun2;
localVars.npar = npar;
localVars.type=type;
if strcmpi(type,'posterior')
while b= options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
if subplotnum == 1
fprintf(fidTeX,'\\begin{figure}[H]\n');
end
end
if subplotnum == MaxNumberOfPlotPerFigure || (jj == nvar && subplotnum> 0)
figunumber = figunumber+1;
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s/%s_Bayesian_IRF_%s_%d}\n',options_.figures.textwidth*min(subplotnum/nn,1),DirectoryName,M_.fname,tit{ii},figunumber);
if options_.relative_irf
fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
else
fprintf(fidTeX,'\\caption{Bayesian IRF: Orthogonalized shock to $%s$.}\n',titTeX{ii});
end
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%d}\n', tit{ii},figunumber);
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
subplotnum = 0;
end
end
end
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
% The others file format are generated in parallel by PosteriorIRF_core2!
if ~isoctave
if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(length(varlist)/MaxNumberOfPlotPerFigure))<8
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
else
isRemoteOctave = 0;
for indPC=1:length(options_.parallel)
isRemoteOctave = isRemoteOctave + (findstr(options_.parallel(indPC).MatlabOctavePath, 'octave'));
end
if isRemoteOctave
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
else
globalVars = struct('M_',M_, ...
'options_', options_);
[fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
end
end
else
[fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
end
% END parallel code!
end
fprintf('Estimation::mcmc: Posterior IRFs, done!\n');