1function PosteriorIRF(type) 2% Builds posterior IRFs after the MH algorithm. 3% 4% INPUTS 5% o type [char] string specifying the joint density of the 6% deep parameters ('prior','posterior'). 7% 8% OUTPUTS 9% None (oo_ and plots). 10% 11% SPECIAL REQUIREMENTS 12% None 13 14% PARALLEL CONTEXT 15% This funtion has been parallelized in two different points. Then we have two core 16% functions associated with it(the _core1 and _core2). 17% See also the comments posterior_sampler.m funtion. 18 19% Copyright (C) 2006-2018 Dynare Team 20% 21% This file is part of Dynare. 22% 23% Dynare is free software: you can redistribute it and/or modify 24% it under the terms of the GNU General Public License as published by 25% the Free Software Foundation, either version 3 of the License, or 26% (at your option) any later version. 27% 28% Dynare is distributed in the hope that it will be useful, 29% but WITHOUT ANY WARRANTY; without even the implied warranty of 30% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 31% GNU General Public License for more details. 32% 33% You should have received a copy of the GNU General Public License 34% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 35 36 37global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info 38 39% Set the number of periods 40if isempty(options_.irf) || ~options_.irf 41 options_.irf = 40; 42end 43% Set varlist if necessary 44varlist = options_.varlist; 45if isempty(varlist) 46 varlist = options_.varobs; 47end 48options_.varlist = varlist; 49nvar = length(varlist); 50IndxVariables = []; 51for i=1:nvar 52 idx = strmatch(varlist{i}, M_.endo_names, 'exact'); 53 if isempty(idx) 54 disp(['PosteriorIRF :: ' varlist{i} 'is not a declared endogenous variable!']) 55 else 56 IndxVariables = [IndxVariables, idx]; 57 end 58end 59 60% Get index of shocks for requested IRFs 61irf_shocks_indx = getIrfShocksIndx(); 62 63% Set various parameters & Check or create directories 64nvx = estim_params_.nvx; 65nvn = estim_params_.nvn; 66ncx = estim_params_.ncx; 67ncn = estim_params_.ncn; 68np = estim_params_.np ; 69npar = nvx+nvn+ncx+ncn+np; 70offset = npar-np; clear('nvx','nvn','ncx','ncn','np'); 71 72nvobs = dataset_.vobs; 73gend = dataset_.nobs; 74MaxNumberOfPlotPerFigure = 9; 75nn = sqrt(MaxNumberOfPlotPerFigure); 76MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8); 77MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); 78if options_.dsge_var 79 MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8); 80else 81 MAX_nirfs_dsgevar = 0; 82end 83 84DirectoryName = CheckPath('Output',M_.dname); 85if strcmpi(type,'posterior') 86 MhDirectoryName = CheckPath('metropolis',M_.dname); 87elseif strcmpi(type,'gsa') 88 if options_.opt_gsa.pprior 89 MhDirectoryName = CheckPath(['GSA' filesep 'prior'],M_.dname); 90 else 91 MhDirectoryName = CheckPath(['GSA' filesep 'mc'],M_.dname); 92 end 93else 94 MhDirectoryName = CheckPath('prior',M_.dname); 95end 96 97%delete old stale files before creating new ones 98delete_stale_file([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']); 99delete_stale_file([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']); 100delete_stale_file([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']); 101delete_stale_file([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat']); 102delete_stale_file([MhDirectoryName filesep M_.fname '_param_irf*.mat']); 103 104 105if strcmpi(type,'posterior') 106 B = options_.sub_draws; 107 options_.B = B; 108 if round((1-options_.mh_conf_sig)*B)<2 109 fprintf('\nPosteriorIRF:: options_.mh_conf_sig times options_.sub_draws is too small to generate HPDIs. I am omitting them.\n') 110 end 111elseif strcmpi(type,'gsa') 112 RootDirectoryName = CheckPath('gsa',M_.dname); 113 if options_.opt_gsa.pprior 114 load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable') 115 else 116 load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable') 117 end 118 x=[lpmat0(istable,:) lpmat(istable,:)]; 119 clear lpmat istable 120 B=size(x,1); options_.B = B; 121else% type = 'prior' 122 B = options_.prior_draws; 123 options_.B = B; 124end 125 126irun = 0; 127IRUN = 0; 128irun2 = 0; 129NumberOfIRFfiles_dsge = 1; 130NumberOfIRFfiles_dsgevar = 1; 131ifil2 = 1; 132% Create arrays 133if B <= MAX_nruns 134 stock_param = zeros(B, npar); 135else 136 stock_param = zeros(MAX_nruns, npar); 137end 138if B >= MAX_nirfs_dsge 139 stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge); 140else 141 stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B); 142end 143if MAX_nirfs_dsgevar 144 if B >= MAX_nirfs_dsgevar 145 stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar); 146 else 147 stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B); 148 end 149 NumberOfLags = options_.dsge_varlag; 150 NumberOfLagsTimesNvobs = NumberOfLags*nvobs; 151 if options_.noconstant 152 NumberOfParametersPerEquation = NumberOfLagsTimesNvobs; 153 else 154 NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1; 155 end 156 Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs); 157end 158 159% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m 160% function. 161 162b = 0; 163 164localVars=[]; 165 166% Save the local variables. 167 168localVars.IRUN = IRUN; 169localVars.irun = irun; 170localVars.irun2=irun2; 171localVars.npar = npar; 172 173localVars.type=type; 174if strcmpi(type,'posterior') 175 while b<B 176 b = b + 1; 177 x(b,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_); 178 end 179end 180 181if ~strcmpi(type,'prior') 182 localVars.x=x; 183end 184 185b=0; 186if options_.dsge_var 187 localVars.gend = gend; 188 localVars.nvobs = nvobs; 189 localVars.NumberOfParametersPerEquation = NumberOfParametersPerEquation; 190 localVars.NumberOfLags = options_.dsge_varlag; 191 localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs; 192 localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs); 193end 194localVars.nvar=nvar; 195localVars.IndxVariables=IndxVariables; 196localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar; 197localVars.MAX_nirfs_dsge=MAX_nirfs_dsge; 198localVars.MAX_nruns=MAX_nruns; 199localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge; 200localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; 201localVars.ifil2=ifil2; 202localVars.MhDirectoryName=MhDirectoryName; 203 204% Like sequential execution! 205if isnumeric(options_.parallel) 206 [fout] = PosteriorIRF_core1(localVars,1,B,0); 207 nosaddle = fout.nosaddle; 208else 209 % Parallel execution! 210 [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B); 211 for j=1:totCPU-1 212 nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge); 213 NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles; 214 if MAX_nirfs_dsgevar 215 nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar); 216 else 217 nfiles=0; 218 end 219 NumberOfIRFfiles_dsgevar(j+1) =NumberOfIRFfiles_dsgevar(j)+nfiles; 220 nfiles = ceil(nBlockPerCPU(j)/MAX_nruns); 221 ifil2(j+1) =ifil2(j)+nfiles; 222 end 223 localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge; 224 localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; 225 localVars.ifil2=ifil2; 226 227 globalVars = struct('M_',M_, ... 228 'options_', options_, ... 229 'bayestopt_', bayestopt_, ... 230 'estim_params_', estim_params_, ... 231 'oo_', oo_, ... 232 'dataset_',dataset_, ... 233 'dataset_info',dataset_info); 234 235 % which files have to be copied to run remotely 236 NamFileInput(1,:) = {'',[M_.fname '.static.m']}; 237 NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']}; 238 if M_.set_auxiliary_variables 239 NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']}; 240 end 241 if options_.steadystate_flag 242 if options_.steadystate_flag == 1 243 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']}; 244 else 245 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '.steadystate.m']}; 246 end 247 end 248 [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info); 249 nosaddle=0; 250 for j=1:length(fout) 251 nosaddle = nosaddle + fout(j).nosaddle; 252 end 253 254end 255 256% END first parallel section! 257 258if nosaddle 259 disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) 260end 261 262ReshapeMatFiles('irf_dsge',type) 263if MAX_nirfs_dsgevar 264 ReshapeMatFiles('irf_bvardsge') 265end 266 267if strcmpi(type,'gsa') 268 return 269end 270 271IRF_DSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']); 272NumberOfIRFfiles_dsge = length(IRF_DSGEs); 273 274IRF_BVARDSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']); 275NumberOfIRFfiles_dsgevar = length(IRF_BVARDSGEs); 276 277MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr); 278MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr); 279VarIRF = zeros(options_.irf,nvar,M_.exo_nbr); 280DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr); 281HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr); 282 283if options_.TeX 284 varlist_TeX = cell(nvar, 1); 285 for i=1:nvar 286 varlist_TeX(i) = {M_.endo_names_tex{IndxVariables(i)}}; 287 end 288end 289 290fprintf('Estimation::mcmc: Posterior (dsge) IRFs...\n'); 291tit(M_.exo_names_orig_ord) = M_.exo_names; 292kdx = 0; 293 294for file = 1:NumberOfIRFfiles_dsge 295 load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']); 296 for i = irf_shocks_indx 297 for j = 1:nvar 298 for k = 1:size(STOCK_IRF_DSGE,1) 299 kk = k+kdx; 300 [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),... 301 DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig); 302 end 303 end 304 end 305 kdx = kdx + size(STOCK_IRF_DSGE,1); 306 307end 308 309clear STOCK_IRF_DSGE; 310 311for i = irf_shocks_indx 312 for j = 1:nvar 313 name = sprintf('%s_%s', M_.endo_names{IndxVariables(j)}, tit{i}); 314 oo_.PosteriorIRF.dsge.Mean.(name) = MeanIRF(:,j,i); 315 oo_.PosteriorIRF.dsge.Median.(name) = MedianIRF(:,j,i); 316 oo_.PosteriorIRF.dsge.Var.(name) = VarIRF(:,j,i); 317 oo_.PosteriorIRF.dsge.deciles.(name) = DistribIRF(:,:,j,i); 318 oo_.PosteriorIRF.dsge.HPDinf.(name) = HPDIRF(:,1,j,i); 319 oo_.PosteriorIRF.dsge.HPDsup.(name) = HPDIRF(:,2,j,i); 320 end 321end 322 323 324if MAX_nirfs_dsgevar 325 MeanIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr); 326 MedianIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr); 327 VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr); 328 DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr); 329 HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr); 330 fprintf('Estimation::mcmc: Posterior (bvar-dsge) IRFs...\n'); 331 tit(M_.exo_names_orig_ord) = M_.exo_names; 332 kdx = 0; 333 for file = 1:NumberOfIRFfiles_dsgevar 334 load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']); 335 for i = irf_shocks_indx 336 for j = 1:nvar 337 for k = 1:size(STOCK_IRF_BVARDSGE,1) 338 kk = k+kdx; 339 [MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),... 340 HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ... 341 posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig); 342 end 343 end 344 end 345 kdx = kdx + size(STOCK_IRF_BVARDSGE,1); 346 end 347 clear STOCK_IRF_BVARDSGE; 348 for i = irf_shocks_indx 349 for j = 1:nvar 350 name = sprintf('%s_%s', M_.endo_names{IndxVariables(j)}, tit{i}); 351 oo_.PosteriorIRF.bvardsge.Mean.(name) = MeanIRFdsgevar(:,j,i); 352 oo_.PosteriorIRF.bvardsge.Median.(name) = MedianIRFdsgevar(:,j,i); 353 oo_.PosteriorIRF.bvardsge.Var.(name) = VarIRFdsgevar(:,j,i); 354 oo_.PosteriorIRF.bvardsge.deciles.(name) = DistribIRFdsgevar(:,:,j,i); 355 oo_.PosteriorIRF.bvardsge.HPDinf.(name) = HPDIRFdsgevar(:,1,j,i); 356 oo_.PosteriorIRF.bvardsge.HPDsup.(name) = HPDIRFdsgevar(:,2,j,i); 357 end 358 end 359end 360% 361% Finally I build the plots. 362% 363 364 365% Second block of code executed in parallel, with the exception of file 366% .tex generation always run in sequentially. This portion of code is execute in parallel by 367% PosteriorIRF_core2.m function. 368 369if ~options_.nograph && ~options_.no_graph.posterior 370 % Save the local variables. 371 localVars=[]; 372 373 Check=options_.TeX; 374 if (Check) 375 localVars.varlist_TeX=varlist_TeX; 376 end 377 378 379 localVars.nvar=nvar; 380 localVars.MeanIRF=MeanIRF; 381 localVars.tit=tit; 382 localVars.nn=nn; 383 localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar; 384 localVars.HPDIRF=HPDIRF; 385 localVars.varlist=varlist; 386 localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure; 387 if options_.dsge_var 388 localVars.HPDIRFdsgevar=HPDIRFdsgevar; 389 localVars.MeanIRFdsgevar = MeanIRFdsgevar; 390 end 391 392 % The files .TeX are genereted in sequential way always! 393 394 % The files .TeX are generated in sequential way always! 395 subplotnum = 0; 396 titTeX(M_.exo_names_orig_ord) = M_.exo_names_tex; 397 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 398 fidTeX = fopen([DirectoryName filesep M_.fname '_BayesianIRF.tex'],'w'); 399 fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n'); 400 fprintf(fidTeX,['%% ' datestr(now,0) '\n']); 401 fprintf(fidTeX,' \n'); 402 titTeX(M_.exo_names_orig_ord) = M_.exo_names_tex; 403 404 for ii=irf_shocks_indx 405 figunumber = 0; 406 407 for jj=1:nvar 408 if max(abs(MeanIRF(:,jj,ii))) >= options_.impulse_responses.plot_threshold 409 subplotnum = subplotnum+1; 410 411 if subplotnum == 1 412 fprintf(fidTeX,'\\begin{figure}[H]\n'); 413 end 414 name = varlist{jj}; 415 texname = varlist_TeX{jj}; 416 fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],name,['$' texname '$']); 417 end 418 if subplotnum == MaxNumberOfPlotPerFigure || (jj == nvar && subplotnum> 0) 419 figunumber = figunumber+1; 420 421 fprintf(fidTeX,'\\centering \n'); 422 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); 423 if options_.relative_irf 424 fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']); 425 else 426 fprintf(fidTeX,'\\caption{Bayesian IRF: Orthogonalized shock to $%s$.}\n',titTeX{ii}); 427 end 428 fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%d}\n', tit{ii},figunumber); 429 fprintf(fidTeX,'\\end{figure}\n'); 430 fprintf(fidTeX,' \n'); 431 subplotnum = 0; 432 end 433 end 434 end 435 fprintf(fidTeX,'%% End of TeX file.\n'); 436 fclose(fidTeX); 437 end 438 439 % The others file format are generated in parallel by PosteriorIRF_core2! 440 441 442 % Comment for testing! 443 if ~isoctave 444 if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(length(varlist)/MaxNumberOfPlotPerFigure))<8 445 [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); 446 else 447 isRemoteOctave = 0; 448 for indPC=1:length(options_.parallel) 449 isRemoteOctave = isRemoteOctave + (findstr(options_.parallel(indPC).MatlabOctavePath, 'octave')); 450 end 451 if isRemoteOctave 452 [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); 453 else 454 globalVars = struct('M_',M_, ... 455 'options_', options_); 456 457 [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info); 458 end 459 end 460 else 461 [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); 462 end 463 % END parallel code! 464 465end 466 467fprintf('Estimation::mcmc: Posterior IRFs, done!\n'); 468