1function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
2%   Generates and stores Posterior IRFs
3%   PARALLEL CONTEXT
4%   This function perfoms in parallel execution a portion of the PosteriorIRF.m code.
5%   This is a special kind of parallel function. Unlike of other parallel functions,
6%   that running in parallel a 'for' cycle, this function run in parallel a
7%   'while' loop! The parallelization of 'while' loop (when possible) is a more
8%   sophisticated procedure.
9%
10%   See also the comment in posterior_sampler_core.m funtion.
11%
12% INPUTS
13%   See the comment in posterior_sampler_core.m funtion.
14%
15% OUTPUTS
16% o myoutput  [struc]
17%  Contained:
18%  OutputFileName_dsge, OutputFileName_param and OutputFileName_bvardsge.
19%
20% ALGORITHM
21%   Portion of PosteriorIRF.m function. Specifically the 'while' cycle.
22%
23% SPECIAL REQUIREMENTS.
24%   None.
25%
26% Copyright (C) 2006-2019 Dynare Team
27%
28% This file is part of Dynare.
29%
30% Dynare is free software: you can redistribute it and/or modify
31% it under the terms of the GNU General Public License as published by
32% the Free Software Foundation, either version 3 of the License, or
33% (at your option) any later version.
34%
35% Dynare is distributed in the hope that it will be useful,
36% but WITHOUT ANY WARRANTY; without even the implied warranty of
37% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
38% GNU General Public License for more details.
39%
40% You should have received a copy of the GNU General Public License
41% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
42
43
44global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
45
46if nargin<4
47    whoiam=0;
48end
49
50% Reshape 'myinputs' for local computation.
51% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
52
53IRUN = myinputs.IRUN;
54irun =myinputs.irun;
55irun2=myinputs.irun2;
56npar=myinputs.npar;
57type=myinputs.type;
58if ~strcmpi(type,'prior')
59    x=myinputs.x;
60end
61
62nvar=myinputs.nvar;
63IndxVariables=myinputs.IndxVariables;
64MAX_nirfs_dsgevar=myinputs.MAX_nirfs_dsgevar;
65MAX_nirfs_dsge=myinputs.MAX_nirfs_dsge;
66MAX_nruns=myinputs.MAX_nruns;
67
68NumberOfIRFfiles_dsge=myinputs.NumberOfIRFfiles_dsge;
69NumberOfIRFfiles_dsgevar=myinputs.NumberOfIRFfiles_dsgevar;
70ifil2=myinputs.ifil2;
71
72if options_.dsge_var
73    nvobs=myinputs.nvobs;
74    NumberOfParametersPerEquation = myinputs.NumberOfParametersPerEquation;
75    NumberOfLagsTimesNvobs = myinputs.NumberOfLagsTimesNvobs;
76    Companion_matrix = myinputs.Companion_matrix;
77    stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
78    bounds = prior_bounds(bayestopt_,options_.prior_trunc);
79end
80
81
82if whoiam
83    Parallel=myinputs.Parallel;
84end
85
86
87% MhDirectoryName = myinputs.MhDirectoryName;
88if strcmpi(type,'posterior')
89    MhDirectoryName = CheckPath('metropolis',M_.dname);
90elseif strcmpi(type,'gsa')
91    if options_.opt_gsa.pprior
92        MhDirectoryName = CheckPath(['gsa' filesep 'prior'],M_.dname);
93    else
94        MhDirectoryName = CheckPath(['gsa' filesep 'mc'],M_.dname);
95    end
96else
97    MhDirectoryName = CheckPath('prior',M_.dname);
98end
99
100RemoteFlag = 0;
101
102if whoiam
103    if Parallel(ThisMatlab).Local==0
104        RemoteFlag =1;
105    end
106    prct0={0,whoiam,Parallel(ThisMatlab)};
107else
108    prct0=0;
109end
110if strcmpi(type,'posterior')
111    h = dyn_waitbar(prct0,'Bayesian (posterior) IRFs...');
112elseif strcmpi(type,'gsa')
113    h = dyn_waitbar(prct0,'GSA (prior) IRFs...');
114else
115    h = dyn_waitbar(prct0,'Bayesian (prior) IRFs...');
116end
117
118
119OutputFileName_bvardsge = {};
120OutputFileName_dsge = {};
121OutputFileName_param = {};
122
123
124fpar = fpar-1;
125fpar0=fpar;
126nosaddle=0;
127
128if whoiam
129    ifil2=ifil2(whoiam);
130    NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam);
131    NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam);
132end
133
134% Parallel 'while' very good!!!
135stock_param=zeros(MAX_nruns,npar);
136stock_irf_dsge=zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge);
137while fpar<B
138    fpar = fpar + 1;
139    irun = irun+1;
140    irun2 = irun2+1;
141    if strcmpi(type,'prior')
142        deep = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
143    else
144        deep = x(fpar,:);
145    end
146    stock_param(irun2,:) = deep;
147    set_parameters(deep);
148    [dr,info,M_,options_,oo_] =compute_decision_rules(M_,options_,oo_);
149    oo_.dr = dr;
150    if info(1)
151        nosaddle = nosaddle + 1;
152        fpar = fpar - 1;
153        irun = irun-1;
154        irun2 = irun2-1;
155        if info(1) == 1
156            errordef = 'Static variables are not uniquely defined';
157        elseif info(1) == 2
158            errordef = 'Dll problem';
159        elseif info(1) == 3
160            errordef = 'No stable trajectory';
161        elseif info(1) == 4
162            errordef = 'Indeterminacy';
163        elseif info(1) == 5
164            errordef = 'Rank condition  is not satisfied';
165        else
166            errordef = get_error_message(info, options_);
167        end
168        if strcmpi(type,'prior')
169            disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')'])
170            continue
171        else
172            error(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ') with sample ' type])
173        end
174    end
175    SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
176    SS = transpose(chol(SS));
177    irf_shocks_indx = getIrfShocksIndx();
178    for i=irf_shocks_indx
179        if SS(i,i) > 1e-13
180            if options_.order>1 && options_.relative_irf % normalize shock to 0.01 before IRF generation for GIRFs; multiply with 100 later
181                y=irf(M_,options_,dr,SS(M_.exo_names_orig_ord,i)./SS(i,i)/100, options_.irf, options_.drop,options_.replic,options_.order);
182            else
183                y=irf(M_,options_,dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order);
184            end
185            if options_.relative_irf && options_.order==1 %multiply with 100 for backward compatibility
186                y = 100*y/SS(i,i);
187            end
188            for j = 1:nvar
189                if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12
190                    stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:));
191                end
192            end
193        end
194    end
195    if MAX_nirfs_dsgevar
196        IRUN = IRUN+1;
197        [~,~,~,~,~,~,~,PHI,SIGMAu,iXX] =  dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
198        dsge_prior_weight = M_.params(strmatch('dsge_prior_weight', M_.param_names));
199        DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
200        SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));
201        explosive_var  = 1;
202        while explosive_var
203            % draw from the marginal posterior of SIGMA
204            SIGMAu_draw = rand_inverse_wishart(dataset_.vobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
205                                               SIGMA_inv_upper_chol);
206            % draw from the conditional posterior of PHI
207            PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.vobs, PHI, ...
208                                          chol(SIGMAu_draw)', chol(iXX)');
209            Companion_matrix(1:dataset_.vobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
210            % Check for stationarity
211            explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
212        end
213        % Get the mean
214        mu = zeros(1,dataset_.vobs);
215        % Get rotation
216        if dsge_prior_weight > 0
217            Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
218            A0 = Atheta(bayestopt_.mfys,:);
219            OMEGAstar = qr2(A0');
220        end
221        SIGMAu_chol = chol(SIGMAu_draw)';
222        SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
223        PHIpower = eye(NumberOfLagsTimesNvobs);
224        irfs = zeros (options_.irf,dataset_.vobs*M_.exo_nbr);
225        tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
226        irfs(1,:) = tmp3(:)';
227        for t = 2:options_.irf
228            PHIpower = Companion_matrix*PHIpower;
229            tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
230            irfs(t,:)  = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
231        end
232        tmp_dsgevar = kron(ones(options_.irf,1),mu);
233        for j = 1:(dataset_.vobs*M_.exo_nbr)
234            if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
235                tmp_dsgevar(:,j) = (irfs(:,j));
236            end
237        end
238        if IRUN < MAX_nirfs_dsgevar
239            stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
240        else
241            stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
242            save([MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat'], 'stock_irf_bvardsge');
243            if RemoteFlag==1
244                OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
245            end
246            NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
247            IRUN =0;
248        end
249    end
250    if irun == MAX_nirfs_dsge || irun == B || fpar == B
251        if fpar == B
252            stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
253            if MAX_nirfs_dsgevar && (fpar == B || IRUN == B)
254                stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
255                save([MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat'], 'stock_irf_bvardsge');
256                NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
257                if RemoteFlag==1
258                    OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}];
259                end
260                irun = 0;
261            end
262        end
263        save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge');
264        if RemoteFlag==1
265            OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}];
266        end
267        NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
268        irun = 0;
269    end
270    if irun2 == MAX_nruns || fpar == B
271        if fpar == B
272            stock_param = stock_param(1:irun2,:);
273        end
274        stock = stock_param;
275        save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock');
276        if RemoteFlag==1
277            OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}];
278        end
279        ifil2 = ifil2 + 1;
280        irun2 = 0;
281    end
282    dyn_waitbar((fpar-fpar0)/(B-fpar0),h);
283end
284
285dyn_waitbar_close(h);
286
287if whoiam==0
288    if nosaddle
289        disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
290    end
291end
292
293% Copy the rusults of computation on the call machine (specifically in the
294% directory on call machine that contain the model).
295
296myoutput.OutputFileName = [OutputFileName_dsge;
297                    OutputFileName_param;
298                    OutputFileName_bvardsge];
299
300myoutput.nosaddle = nosaddle;
301