PosteriorIRF

 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

0001 function PosteriorIRF(type)
0002 % Builds posterior IRFs after the MH algorithm.
0003 %
0004 % INPUTS
0005 %   o type       [char]     string specifying the joint density of the
0006 %                           deep parameters ('prior','posterior').
0007 %
0008 % OUTPUTS
0009 %   None                    (oo_ and plots).
0010 %
0011 % SPECIAL REQUIREMENTS
0012 %   None
0013 
0014 % PARALLEL CONTEXT
0015 % This funtion has been parallelized in two different points. Then we have two core
0016 % functions associated with it(the _core1 and _core2).
0017 % See also the comments random_walk_metropolis_hastings.m funtion.
0018 
0019 % Copyright (C) 2006-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 
0037 global options_ estim_params_ oo_ M_ bayestopt_ dataset_
0038 % Set the number of periods
0039 if isempty(options_.irf) || ~options_.irf 
0040     options_.irf = 40;
0041 end
0042 % Set varlist if necessary
0043 varlist = options_.varlist;
0044 if isempty(varlist)
0045     varlist = options_.varobs;
0046 end
0047 options_.varlist = varlist;
0048 nvar = size(varlist,1);
0049 IndxVariables = [];
0050 for i=1:nvar
0051     idx = strmatch(deblank(varlist(i,:)),M_.endo_names,'exact');
0052     if isempty(idx)
0053         disp(['PosteriorIRF :: ' deblank(varlist(i,:)) 'is not a declared endogenous variable!'])
0054     else
0055         IndxVariables = [IndxVariables,idx];
0056     end
0057 end
0058 % Set various parameters & Check or create directories
0059 nvx  = estim_params_.nvx;
0060 nvn  = estim_params_.nvn;
0061 ncx  = estim_params_.ncx;
0062 ncn  = estim_params_.ncn;
0063 np   = estim_params_.np ;
0064 npar = nvx+nvn+ncx+ncn+np;
0065 offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
0066 
0067 nvobs = dataset_.info.nvobs;
0068 gend = dataset_.info.ntobs;
0069 MaxNumberOfPlotPerFigure = 9;
0070 nn = sqrt(MaxNumberOfPlotPerFigure);
0071 MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
0072 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0073 if options_.dsge_var
0074     MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8);
0075 else
0076     MAX_nirfs_dsgevar = 0;
0077 end
0078 
0079 DirectoryName = CheckPath('Output',M_.dname);
0080 if strcmpi(type,'posterior')
0081     MhDirectoryName = CheckPath('metropolis',M_.dname);
0082 elseif strcmpi(type,'gsa')
0083     if options_.opt_gsa.pprior
0084         MhDirectoryName = CheckPath(['GSA' filesep 'prior'],M_.dname);
0085     else
0086         MhDirectoryName = CheckPath(['GSA' filesep 'mc'],M_.dname);
0087     end
0088 else
0089     MhDirectoryName = CheckPath('prior',M_.dname);
0090 end
0091 delete([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
0092 delete([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
0093 if strcmpi(type,'posterior')
0094     load([ MhDirectoryName filesep  M_.fname '_mh_history.mat'])
0095     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0096     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0097 elseif strcmpi(type,'gsa')
0098     RootDirectoryName = CheckPath('gsa',M_.dname);
0099     if options_.opt_gsa.pprior
0100         load([ RootDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
0101     else
0102         load([ RootDirectoryName filesep  M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
0103     end
0104     x=[lpmat0(istable,:) lpmat(istable,:)];
0105     clear lpmat istable
0106     NumberOfDraws=size(x,1);
0107     B=NumberOfDraws; options_.B = B;
0108 else% type = 'prior'
0109     NumberOfDraws = 500;
0110 end
0111 if ~strcmpi(type,'gsa')
0112     B = min([round(.5*NumberOfDraws),500]); options_.B = B;
0113 end
0114 try 
0115     delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
0116 catch 
0117     disp('No _IRFs (dsge) files to be deleted!')
0118 end
0119 try 
0120     delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
0121 catch 
0122     disp('No _IRFs (bvar-dsge) files to be deleted!')
0123 end
0124 irun = 0;
0125 IRUN = 0;
0126 irun2 = 0;
0127 NumberOfIRFfiles_dsge = 1;
0128 NumberOfIRFfiles_dsgevar = 1;
0129 ifil2 = 1;
0130 % Create arrays
0131 if B <= MAX_nruns
0132     stock_param = zeros(B, npar);
0133 else
0134     stock_param = zeros(MAX_nruns, npar);
0135 end
0136 if B >= MAX_nirfs_dsge
0137     stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
0138 else
0139     stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
0140 end
0141 if MAX_nirfs_dsgevar
0142     if B >= MAX_nirfs_dsgevar
0143         stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
0144     else
0145         stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B);
0146     end
0147     [mYY,mXY,mYX,mXX,Ydata,Xdata] = ...
0148         var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,...
0149                            options_.dsge_varlag,-1,options_.datafile,options_.varobs,...
0150                            options_.xls_sheet,options_.xls_range);
0151     NumberOfLags = options_.dsge_varlag;
0152     NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
0153     if options_.noconstant
0154         NumberOfParametersPerEquation = NumberOfLagsTimesNvobs;
0155     else
0156         NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1;
0157     end
0158     Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
0159 end
0160 
0161 % First block of code executed in parallel. The function devoted to do it is  PosteriorIRF_core1.m
0162 % function.
0163 
0164 b = 0;
0165 nosaddle = 0;
0166 
0167 localVars=[];
0168 
0169 % Save the local variables.
0170 
0171 localVars.IRUN = IRUN;
0172 localVars.irun = irun;
0173 localVars.irun2=irun2;
0174 localVars.nosaddle=nosaddle;
0175 
0176 localVars.type=type;
0177 if strcmpi(type,'posterior')
0178     while b<B 
0179         b = b + 1;
0180         x(b,:) = GetOneDraw(type);
0181     end
0182 end
0183 
0184 if ~strcmpi(type,'prior'),
0185     localVars.x=x;
0186 end
0187 
0188 b=0;
0189 if options_.dsge_var
0190     localVars.gend = gend;
0191     localVars.nvobs = nvobs;
0192     localVars.NumberOfParametersPerEquation = NumberOfParametersPerEquation;
0193     localVars.NumberOfLags = options_.dsge_varlag;
0194     localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
0195     localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
0196 end 
0197 localVars.nvar=nvar;
0198 localVars.IndxVariables=IndxVariables;
0199 localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
0200 localVars.MAX_nirfs_dsge=MAX_nirfs_dsge;
0201 localVars.MAX_nruns=MAX_nruns;
0202 localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
0203 localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
0204 localVars.ifil2=ifil2;
0205 localVars.MhDirectoryName=MhDirectoryName;
0206 
0207 % Like sequential execution!
0208 if isnumeric(options_.parallel),
0209     [fout] = PosteriorIRF_core1(localVars,1,B,0);
0210 else
0211     % Parallel execution!
0212     [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
0213     for j=1:totCPU-1,
0214         nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge);
0215         NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles;
0216         if MAX_nirfs_dsgevar,
0217             nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar);
0218         else
0219             nfiles=0;
0220         end
0221         NumberOfIRFfiles_dsgevar(j+1) =NumberOfIRFfiles_dsgevar(j)+nfiles;
0222         nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
0223         ifil2(j+1) =ifil2(j)+nfiles;
0224     end
0225     localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
0226     localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
0227     localVars.ifil2=ifil2;
0228     
0229     globalVars = struct('M_',M_, ...
0230                         'options_', options_, ...
0231                         'bayestopt_', bayestopt_, ...
0232                         'estim_params_', estim_params_, ...
0233                         'oo_', oo_, ...
0234                         'dataset_',dataset_);
0235     
0236     % which files have to be copied to run remotely
0237     NamFileInput(1,:) = {'',[M_.fname '_static.m']};
0238     NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
0239     if options_.steadystate_flag,
0240         NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
0241     end
0242     [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info);
0243     
0244 end
0245 
0246 % END first parallel section!
0247 
0248 if nosaddle
0249     disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) 
0250 end
0251 
0252 ReshapeMatFiles('irf_dsge',type)
0253 if MAX_nirfs_dsgevar
0254     ReshapeMatFiles('irf_bvardsge')
0255 end
0256 
0257 if strcmpi(type,'gsa')
0258     return
0259 end
0260 
0261 IRF_DSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
0262 NumberOfIRFfiles_dsge = length(IRF_DSGEs);
0263 
0264 IRF_BVARDSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
0265 NumberOfIRFfiles_dsgevar = length(IRF_BVARDSGEs);
0266 
0267 MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
0268 MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
0269 VarIRF = zeros(options_.irf,nvar,M_.exo_nbr);
0270 DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
0271 HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
0272 
0273 if options_.TeX
0274     for i=1:nvar
0275         if i==1
0276             varlist_TeX = M_.endo_names_tex(IndxVariables(i),:);
0277         else
0278             varlist_TeX = char(varlist_TeX,M_.endo_names_tex(IndxVariables(i),:));
0279         end
0280     end
0281 end
0282 
0283 fprintf('MH: Posterior (dsge) IRFs...\n');
0284 tit(M_.exo_names_orig_ord,:) = M_.exo_names;
0285 kdx = 0;
0286 
0287 for file = 1:NumberOfIRFfiles_dsge
0288     load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
0289     for i = 1:M_.exo_nbr
0290         for j = 1:nvar
0291             for k = 1:size(STOCK_IRF_DSGE,1)
0292                 kk = k+kdx;
0293                 [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),...
0294                  DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig);
0295             end
0296         end
0297     end
0298     kdx = kdx + size(STOCK_IRF_DSGE,1);
0299 
0300 end
0301 
0302 clear STOCK_IRF_DSGE;
0303 
0304 for i = 1:M_.exo_nbr
0305     for j = 1:nvar
0306         name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
0307         eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
0308         eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
0309         eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
0310         eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']);
0311         eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
0312         eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
0313     end
0314 end
0315 
0316 
0317 if MAX_nirfs_dsgevar
0318     MeanIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
0319     MedianIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
0320     VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
0321     DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr);
0322     HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);    
0323     fprintf('MH: Posterior (bvar-dsge) IRFs...\n');
0324     tit(M_.exo_names_orig_ord,:) = M_.exo_names;
0325     kdx = 0;
0326     for file = 1:NumberOfIRFfiles_dsgevar
0327         load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
0328         for i = 1:M_.exo_nbr
0329             for j = 1:nvar
0330                 for k = 1:size(STOCK_IRF_BVARDSGE,1)
0331                     kk = k+kdx;
0332                     [MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
0333                      HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
0334                         posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig);
0335                 end
0336             end
0337         end
0338         kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
0339     end
0340     clear STOCK_IRF_BVARDSGE; 
0341     for i = 1:M_.exo_nbr
0342         for j = 1:nvar
0343             name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
0344             eval(['oo_.PosteriorIRF.bvardsge.Mean.' name ' = MeanIRFdsgevar(:,j,i);']);
0345             eval(['oo_.PosteriorIRF.bvardsge.Median.' name ' = MedianIRFdsgevar(:,j,i);']);
0346             eval(['oo_.PosteriorIRF.bvardsge.Var.' name ' = VarIRFdsgevar(:,j,i);']);
0347             eval(['oo_.PosteriorIRF.bvardsge.Distribution.' name ' = DistribIRFdsgevar(:,:,j,i);']);
0348             eval(['oo_.PosteriorIRF.bvardsge.HPDinf.' name ' = HPDIRFdsgevar(:,1,j,i);']);
0349             eval(['oo_.PosteriorIRF.bvardsge.HPDsup.' name ' = HPDIRFdsgevar(:,2,j,i);']);
0350         end
0351     end
0352 end
0353 %%
0354 %%      Finally I build the plots.
0355 %%
0356 
0357 
0358 % Second block of code executed in parallel, with the exception of file
0359 % .tex generation always run in sequentially. This portion of code is execute in parallel by
0360 % PosteriorIRF_core2.m function.
0361 
0362 
0363 % Save the local variables.
0364 localVars=[];
0365 
0366 Check=options_.TeX;
0367 if (Check)
0368     localVars.varlist_TeX=varlist_TeX;
0369 end
0370 
0371 
0372 localVars.nvar=nvar;
0373 localVars.MeanIRF=MeanIRF;
0374 localVars.tit=tit;
0375 localVars.nn=nn;
0376 localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
0377 localVars.HPDIRF=HPDIRF;
0378 localVars.varlist=varlist;
0379 localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
0380 if options_.dsge_var
0381     localVars.HPDIRFdsgevar=HPDIRFdsgevar;
0382     localVars.MeanIRFdsgevar = MeanIRFdsgevar;
0383 end    
0384 
0385 %%% The files .TeX are genereted in sequential way always!
0386 
0387 if options_.TeX
0388     fidTeX = fopen([DirectoryName filesep M_.fname '_BayesianIRF.TeX'],'w');
0389     fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
0390     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
0391     fprintf(fidTeX,' \n');
0392     titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
0393     
0394     for i=1:M_.exo_nbr
0395         NAMES = [];
0396         TEXNAMES = [];
0397         
0398         for j=1:nvar
0399             if max(abs(MeanIRF(:,j,i))) > 10^(-6)  
0400                 
0401                 name = deblank(varlist(j,:));
0402                 texname = deblank(varlist_TeX(j,:));
0403 
0404                 if j==1
0405                     NAMES = name;
0406                     TEXNAMES = ['$' texname '$'];
0407                 else
0408                     NAMES = char(NAMES,name);
0409                     TEXNAMES = char(TEXNAMES,['$' texname '$']);
0410                 end
0411             end
0412             
0413         end
0414         fprintf(fidTeX,'\\begin{figure}[H]\n');
0415         for jj = 1:size(TEXNAMES,1)
0416             fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
0417         end
0418         fprintf(fidTeX,'\\centering \n');
0419         fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
0420         if options_.relative_irf
0421             fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
0422         else
0423             fprintf(fidTeX,'\\caption{Bayesian IRF.}');
0424         end
0425         fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
0426         fprintf(fidTeX,'\\end{figure}\n');
0427         fprintf(fidTeX,' \n');
0428     end
0429     
0430     fprintf(fidTeX,'%% End of TeX file.\n');
0431     fclose(fidTeX);
0432     
0433 end
0434 
0435 % The others file format are generated in parallel by PosteriorIRF_core2!
0436 
0437 
0438 % Comment for testing!
0439 if ~exist('OCTAVE_VERSION')
0440     if isnumeric(options_.parallel)  || (M_.exo_nbr*ceil(size(varlist,1)/MaxNumberOfPlotPerFigure))<8,
0441         [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
0442     else
0443         isRemoteOctave = 0;
0444         for indPC=1:length(options_.parallel),
0445             isRemoteOctave = isRemoteOctave + (findstr(options_.parallel(indPC).MatlabOctavePath, 'octave'));
0446         end
0447         if isRemoteOctave
0448             [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
0449         else
0450             globalVars = struct('M_',M_, ...
0451                 'options_', options_);
0452             
0453             [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
0454         end
0455     end
0456 end
0457 % END parallel code!
0458 
0459 
0460 fprintf('MH: Posterior IRFs, done!\n');