PosteriorIRF.m 15.6 KB
Newer Older
1
function PosteriorIRF(type)
2
3
4
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
5
%   o type       [char]     string specifying the joint density of the
6
7
8
%                           deep parameters ('prior','posterior').
%
% OUTPUTS
9
10
11
12
%   None                    (oo_ and plots).
%
% SPECIAL REQUIREMENTS
%   None
sebastien's avatar
sebastien committed
13

14
15
16
% PARALLEL CONTEXT
% This funtion has been parallelized in two different points. Then we have two core
% 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
sebastien's avatar
sebastien committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
%
% 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 <http://www.gnu.org/licenses/>.
35

36

37
38
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info

39
% Set the number of periods
40
if isempty(options_.irf) || ~options_.irf
41
42
43
44
45
    options_.irf = 40;
end
% Set varlist if necessary
varlist = options_.varlist;
if isempty(varlist)
46
    varlist = options_.varobs;
47
48
end
options_.varlist = varlist;
49
nvar = length(varlist);
50
51
IndxVariables = [];
for 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
58
    end
end
59
60
61

% Get index of shocks for requested IRFs
irf_shocks_indx = getIrfShocksIndx();
62

63
64
65
66
67
68
69
% 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;
70
71
offset = npar-np; clear('nvx','nvn','ncx','ncn','np');

72
73
nvobs = dataset_.vobs;
gend = dataset_.nobs;
74
75
76
77
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);
78
if options_.dsge_var
79
80
81
82
    MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8);
else
    MAX_nirfs_dsgevar = 0;
end
83

84
DirectoryName = CheckPath('Output',M_.dname);
85
if strcmpi(type,'posterior')
86
    MhDirectoryName = CheckPath('metropolis',M_.dname);
87
elseif 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
93
else
94
    MhDirectoryName = CheckPath('prior',M_.dname);
95
end
96
97
98
99
100
101
102
103
104

%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']);


105
if strcmpi(type,'posterior')
Stéphane Adjemian's avatar
Stéphane Adjemian committed
106
107
    B = options_.sub_draws;
    options_.B = B;
108
109
110
    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
111
elseif strcmpi(type,'gsa')
112
    RootDirectoryName = CheckPath('gsa',M_.dname);
113
114
115
116
117
    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
118
119
    x=[lpmat0(istable,:) lpmat(istable,:)];
    clear lpmat istable
Stéphane Adjemian's avatar
Stéphane Adjemian committed
120
    B=size(x,1); options_.B = B;
121
else% type = 'prior'
Stéphane Adjemian's avatar
Stéphane Adjemian committed
122
123
    B = options_.prior_draws;
    options_.B = B;
124
end
125

126
127
128
129
130
131
132
133
irun = 0;
IRUN = 0;
irun2 = 0;
NumberOfIRFfiles_dsge = 1;
NumberOfIRFfiles_dsgevar = 1;
ifil2 = 1;
% Create arrays
if B <= MAX_nruns
134
    stock_param = zeros(B, npar);
135
else
136
    stock_param = zeros(MAX_nruns, npar);
137
138
end
if B >= MAX_nirfs_dsge
139
    stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
140
else
141
    stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B);
142
143
144
145
146
147
148
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
149
    NumberOfLags = options_.dsge_varlag;
150
151
152
153
154
155
156
157
    NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
    if options_.noconstant
        NumberOfParametersPerEquation = NumberOfLagsTimesNvobs;
    else
        NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1;
    end
    Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
158

159
160
% First block of code executed in parallel. The function devoted to do it is  PosteriorIRF_core1.m
% function.
161

162
b = 0;
163
164
165
166
167

localVars=[];

% Save the local variables.

Marco Ratto's avatar
Marco Ratto committed
168
localVars.IRUN = IRUN;
169
170
localVars.irun = irun;
localVars.irun2=irun2;
171
localVars.npar = npar;
172

173
localVars.type=type;
174
if strcmpi(type,'posterior')
175
    while b<B
176
        b = b + 1;
177
        x(b,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
178
    end
179
180
end

181
if ~strcmpi(type,'prior')
182
    localVars.x=x;
183
184
185
end

b=0;
186
187
188
189
190
191
192
if options_.dsge_var
    localVars.gend = gend;
    localVars.nvobs = nvobs;
    localVars.NumberOfParametersPerEquation = NumberOfParametersPerEquation;
    localVars.NumberOfLags = options_.dsge_varlag;
    localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
    localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
193
end
194
195
196
197
198
199
200
201
localVars.nvar=nvar;
localVars.IndxVariables=IndxVariables;
localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
localVars.MAX_nirfs_dsge=MAX_nirfs_dsge;
localVars.MAX_nruns=MAX_nruns;
localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
202
localVars.MhDirectoryName=MhDirectoryName;
203

204
% Like sequential execution!
205
if isnumeric(options_.parallel)
206
    [fout] = PosteriorIRF_core1(localVars,1,B,0);
207
    nosaddle = fout.nosaddle;
208
else
209
    % Parallel execution!
210
    [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
211
    for j=1:totCPU-1
212
213
        nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge);
        NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles;
214
        if MAX_nirfs_dsgevar
215
216
217
218
            nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar);
        else
            nfiles=0;
        end
219
220
221
        NumberOfIRFfiles_dsgevar(j+1) =NumberOfIRFfiles_dsgevar(j)+nfiles;
        nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
        ifil2(j+1) =ifil2(j)+nfiles;
222
    end
223
224
225
    localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
    localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
    localVars.ifil2=ifil2;
226

227
    globalVars = struct('M_',M_, ...
228
229
230
                        'options_', options_, ...
                        'bayestopt_', bayestopt_, ...
                        'estim_params_', estim_params_, ...
231
                        'oo_', oo_, ...
232
233
234
                        'dataset_',dataset_, ...
                        'dataset_info',dataset_info);

235
    % which files have to be copied to run remotely
236
237
    NamFileInput(1,:) = {'',[M_.fname '.static.m']};
    NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']};
238
    if M_.set_auxiliary_variables
239
        NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
240
    end
241
242
    if options_.steadystate_flag
        if options_.steadystate_flag == 1
243
244
            NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
        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
252
        nosaddle = nosaddle + fout(j).nosaddle;
    end
253

254
end
255

256
% END first parallel section!
257

258
if nosaddle
259
    disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
260
261
end

262
ReshapeMatFiles('irf_dsge',type)
263
264
265
266
267
if MAX_nirfs_dsgevar
    ReshapeMatFiles('irf_bvardsge')
end

if strcmpi(type,'gsa')
268
    return
269
270
end

271
IRF_DSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']);
272
273
NumberOfIRFfiles_dsge = length(IRF_DSGEs);

274
IRF_BVARDSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']);
275
276
277
278
279
280
281
282
283
NumberOfIRFfiles_dsgevar = length(IRF_BVARDSGEs);

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
284
    varlist_TeX = cell(nvar, 1);
285
    for i=1:nvar
286
        varlist_TeX(i) = {M_.endo_names_tex{IndxVariables(i)}};
287
    end
288
289
end

290
fprintf('Estimation::mcmc: Posterior (dsge) IRFs...\n');
291
tit(M_.exo_names_orig_ord) = M_.exo_names;
292
293
294
kdx = 0;

for file = 1:NumberOfIRFfiles_dsge
295
    load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']);
296
    for i = irf_shocks_indx
297
298
299
300
301
302
        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,j,i,:)),0,options_.mh_conf_sig);
            end
303
304
        end
    end
305
    kdx = kdx + size(STOCK_IRF_DSGE,1);
306

307
308
309
310
end

clear STOCK_IRF_DSGE;

311
for i = irf_shocks_indx
312
    for j = 1:nvar
313
        name = sprintf('%s_%s', M_.endo_names{IndxVariables(j)}, tit{i});
314
315
316
317
318
319
        oo_.PosteriorIRF.dsge.Mean.(name) = MeanIRF(:,j,i);
        oo_.PosteriorIRF.dsge.Median.(name) = MedianIRF(:,j,i);
        oo_.PosteriorIRF.dsge.Var.(name) = VarIRF(:,j,i);
        oo_.PosteriorIRF.dsge.deciles.(name) = DistribIRF(:,:,j,i);
        oo_.PosteriorIRF.dsge.HPDinf.(name) = HPDIRF(:,1,j,i);
        oo_.PosteriorIRF.dsge.HPDsup.(name) = HPDIRF(:,2,j,i);
320
    end
321
322
323
324
325
326
327
328
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);
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
333
    kdx = 0;
    for file = 1:NumberOfIRFfiles_dsgevar
334
        load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']);
335
        for i = irf_shocks_indx
336
337
338
339
340
341
342
343
344
345
346
            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,options_.mh_conf_sig);
                end
            end
        end
        kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
    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
352
353
354
355
356
            oo_.PosteriorIRF.bvardsge.Mean.(name) = MeanIRFdsgevar(:,j,i);
            oo_.PosteriorIRF.bvardsge.Median.(name) = MedianIRFdsgevar(:,j,i);
            oo_.PosteriorIRF.bvardsge.Var.(name) = VarIRFdsgevar(:,j,i);
            oo_.PosteriorIRF.bvardsge.deciles.(name) = DistribIRFdsgevar(:,:,j,i);
            oo_.PosteriorIRF.bvardsge.HPDinf.(name) = HPDIRFdsgevar(:,1,j,i);
            oo_.PosteriorIRF.bvardsge.HPDsup.(name) = HPDIRFdsgevar(:,2,j,i);
357
358
359
        end
    end
end
360
361
362
%
%      Finally I build the plots.
%
363
364


365
366
367
% Second block of code executed in parallel, with the exception of file
% .tex generation always run in sequentially. This portion of code is execute in parallel by
% PosteriorIRF_core2.m function.
368

369
if ~options_.nograph && ~options_.no_graph.posterior
370
371
    % Save the local variables.
    localVars=[];
372

373
374
375
376
    Check=options_.TeX;
    if (Check)
        localVars.varlist_TeX=varlist_TeX;
    end
377

378

379
380
381
382
383
384
385
386
387
388
389
390
    localVars.nvar=nvar;
    localVars.MeanIRF=MeanIRF;
    localVars.tit=tit;
    localVars.nn=nn;
    localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
    localVars.HPDIRF=HPDIRF;
    localVars.varlist=varlist;
    localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
    if options_.dsge_var
        localVars.HPDIRFdsgevar=HPDIRFdsgevar;
        localVars.MeanIRFdsgevar = MeanIRFdsgevar;
    end
391

392
393
394
    % The files .TeX are genereted in sequential way always!

    subplotnum = 0;
395
    titTeX(M_.exo_names_orig_ord) = M_.exo_names_tex;
396
397
398
399
400
    if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
        fidTeX = fopen([DirectoryName filesep 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');
401
        titTeX(M_.exo_names_orig_ord) = M_.exo_names_tex;
402
403
404
405
406
407
408
409
410
411
412

        for ii=irf_shocks_indx
            figunumber = 0;

            for jj=1:nvar
                if max(abs(MeanIRF(:,jj,ii))) >= options_.impulse_responses.plot_threshold
                    subplotnum = subplotnum+1;

                    if subplotnum == 1
                        fprintf(fidTeX,'\\begin{figure}[H]\n');
                    end
413
                end
414
415
416
417
                if subplotnum == MaxNumberOfPlotPerFigure || (jj == nvar  && subplotnum> 0)
                    figunumber = figunumber+1;

                    fprintf(fidTeX,'\\centering \n');
418
                    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);
419
420
421
                    if options_.relative_irf
                        fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
                    else
422
                        fprintf(fidTeX,'\\caption{Bayesian IRF: Orthogonalized shock to $%s$.}\n',titTeX{ii});
423
                    end
424
                    fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%d}\n', tit{ii},figunumber);
425
426
427
                    fprintf(fidTeX,'\\end{figure}\n');
                    fprintf(fidTeX,' \n');
                    subplotnum = 0;
428
                end
429
            end
430
        end
431
432
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
433
    end
434

435
    % The others file format are generated in parallel by PosteriorIRF_core2!
436

437
    if ~isoctave
438
        if isnumeric(options_.parallel)  || (M_.exo_nbr*ceil(length(varlist)/MaxNumberOfPlotPerFigure))<8
439
440
            [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
        else
441
442
443
444
445
446
447
448
449
            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_);
450

451
452
                [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
            end
453
        end
454
455
    else
        [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0);
456
    end
457
    % END parallel code!
458

459
end
460

461
fprintf('Estimation::mcmc: Posterior IRFs, done!\n');