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