diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index ac479c526cea3d607ff104e1a0917fdc0313e516..fc1c3bcf58ce1b158458f244062f46433da3188f 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -1,11 +1,11 @@ function PosteriorIRF(type) -% Builds posterior IRFs after the MH algorithm. -% -% INPUTS +% Builds posterior IRFs after the MH algorithm. +% +% INPUTS % o type [char] string specifying the joint density of the -% deep parameters ('prior','posterior'). -% -% OUTPUTS +% deep parameters ('prior','posterior'). +% +% OUTPUTS % None (oo_ and plots). % % SPECIAL REQUIREMENTS @@ -34,9 +34,10 @@ function PosteriorIRF(type) % along with Dynare. If not, see <http://www.gnu.org/licenses/>. -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ +global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info + % Set the number of periods -if isempty(options_.irf) || ~options_.irf +if isempty(options_.irf) || ~options_.irf options_.irf = 40; end % Set varlist if necessary @@ -58,7 +59,7 @@ end % Get index of shocks for requested IRFs irf_shocks_indx = getIrfShocksIndx(); - + % Set various parameters & Check or create directories nvx = estim_params_.nvx; nvn = estim_params_.nvn; @@ -68,8 +69,8 @@ np = estim_params_.np ; npar = nvx+nvn+ncx+ncn+np; offset = npar-np; clear('nvx','nvn','ncx','ncn','np'); -nvobs = dataset_.info.nvobs; -gend = dataset_.info.ntobs; +nvobs = dataset_.vobs; +gend = dataset_.nobs; MaxNumberOfPlotPerFigure = 9; nn = sqrt(MaxNumberOfPlotPerFigure); MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8); @@ -111,14 +112,14 @@ else% type = 'prior' B = options_.prior_draws; options_.B = B; end -try +try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']) -catch +catch disp('No _IRFs (dsge) files to be deleted!') end -try +try delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat']) -catch +catch disp('No _IRFs (bvar-dsge) files to be deleted!') end irun = 0; @@ -144,10 +145,6 @@ if MAX_nirfs_dsgevar else stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B); end - [mYY,mXY,mYX,mXX,Ydata,Xdata] = ... - var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,... - options_.dsge_varlag,-1,options_.datafile,options_.varobs,... - options_.xls_sheet,options_.xls_range); NumberOfLags = options_.dsge_varlag; NumberOfLagsTimesNvobs = NumberOfLags*nvobs; if options_.noconstant @@ -174,7 +171,7 @@ localVars.npar = npar; localVars.type=type; if strcmpi(type,'posterior') - while b<B + while b<B b = b + 1; x(b,:) = GetOneDraw(type); end @@ -192,7 +189,7 @@ if options_.dsge_var localVars.NumberOfLags = options_.dsge_varlag; localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs; localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs); -end +end localVars.nvar=nvar; localVars.IndxVariables=IndxVariables; localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar; @@ -225,14 +222,15 @@ else localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge; localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; localVars.ifil2=ifil2; - + globalVars = struct('M_',M_, ... 'options_', options_, ... 'bayestopt_', bayestopt_, ... 'estim_params_', estim_params_, ... 'oo_', oo_, ... - 'dataset_',dataset_); - + 'dataset_',dataset_, ... + 'dataset_info',dataset_info); + % which files have to be copied to run remotely NamFileInput(1,:) = {'',[M_.fname '_static.m']}; NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']}; @@ -244,13 +242,13 @@ else for j=1:length(fout), nosaddle = nosaddle + fout(j).nosaddle; end - + end % END first parallel section! if nosaddle - disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) + disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) end ReshapeMatFiles('irf_dsge',type) @@ -323,7 +321,7 @@ if MAX_nirfs_dsgevar 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); - HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr); + HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr); fprintf('Estimation::mcmc: Posterior (bvar-dsge) IRFs...\n'); tit(M_.exo_names_orig_ord,:) = M_.exo_names; kdx = 0; @@ -341,7 +339,7 @@ if MAX_nirfs_dsgevar end kdx = kdx + size(STOCK_IRF_BVARDSGE,1); end - clear STOCK_IRF_BVARDSGE; + clear STOCK_IRF_BVARDSGE; for i = irf_shocks_indx for j = 1:nvar name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))]; @@ -384,7 +382,7 @@ localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure; if options_.dsge_var localVars.HPDIRFdsgevar=HPDIRFdsgevar; localVars.MeanIRFdsgevar = MeanIRFdsgevar; -end +end % The files .TeX are genereted in sequential way always! @@ -394,14 +392,14 @@ if options_.TeX fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,' \n'); titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex; - + for i=irf_shocks_indx NAMES = []; TEXNAMES = []; - + for j=1:nvar - if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold - + if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold + name = deblank(varlist(j,:)); texname = deblank(varlist_TeX(j,:)); @@ -413,7 +411,7 @@ if options_.TeX TEXNAMES = char(TEXNAMES,['$' texname '$']); end end - + end fprintf(fidTeX,'\\begin{figure}[H]\n'); for jj = 1:size(TEXNAMES,1) @@ -430,10 +428,10 @@ if options_.TeX fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,' \n'); end - + fprintf(fidTeX,'%% End of TeX file.\n'); fclose(fidTeX); - + end % The others file format are generated in parallel by PosteriorIRF_core2! @@ -453,7 +451,7 @@ if ~isoctave else globalVars = struct('M_',M_, ... 'options_', options_); - + [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info); end end diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 1eff387ad64115fcce0e12710c2136ebc612bbc5..0d0a1a752e4ba9921bc39a74ffb61344b9e02167 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -41,7 +41,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) % along with Dynare. If not, see <http://www.gnu.org/licenses/>. -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ +global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info if nargin<4, whoiam=0; @@ -189,24 +189,24 @@ while fpar<B end if MAX_nirfs_dsgevar IRUN = IRUN+1; - [fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); - DSGE_PRIOR_WEIGHT = floor(dataset_.info.ntobs*(1+dsge_prior_weight)); - SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.info.ntobs*(dsge_prior_weight+1))); + DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight)); + SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1))); explosive_var = 1; while explosive_var % draw from the marginal posterior of SIGMA - SIGMAu_draw = rand_inverse_wishart(dataset_.info.nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ... + SIGMAu_draw = rand_inverse_wishart(dataset_.vobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ... SIGMA_inv_upper_chol); % draw from the conditional posterior of PHI - PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.info.nvobs, PHI, ... + PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.vobs, PHI, ... chol(SIGMAu_draw)', chol(iXX)'); - Companion_matrix(1:dataset_.info.nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:)); + Companion_matrix(1:dataset_.vobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:)); % Check for stationarity explosive_var = any(abs(eig(Companion_matrix))>1.000000001); end % Get the mean - mu = zeros(1,dataset_.info.nvobs); + mu = zeros(1,dataset_.vobs); % Get rotation if dsge_prior_weight > 0 Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); @@ -216,24 +216,24 @@ while fpar<B SIGMAu_chol = chol(SIGMAu_draw)'; SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar'; PHIpower = eye(NumberOfLagsTimesNvobs); - irfs = zeros (options_.irf,dataset_.info.nvobs*M_.exo_nbr); - tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA; + irfs = zeros (options_.irf,dataset_.vobs*M_.exo_nbr); + tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA; irfs(1,:) = tmp3(:)'; for t = 2:options_.irf PHIpower = Companion_matrix*PHIpower; - tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA; + tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA; irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu); end tmp_dsgevar = kron(ones(options_.irf,1),mu); - for j = 1:(dataset_.info.nvobs*M_.exo_nbr) + for j = 1:(dataset_.vobs*M_.exo_nbr) if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 tmp_dsgevar(:,j) = (irfs(:,j)); end end if IRUN < MAX_nirfs_dsgevar - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr); + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr); else - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr); + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr); instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];, eval(['save ' instr]); @@ -242,7 +242,7 @@ while fpar<B end NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; IRUN =0; - stock_irf_dsgevar = zeros(options_.irf,dataset_.info.nvobs,M_.exo_nbr,MAX_nirfs_dsgevar); + stock_irf_dsgevar = zeros(options_.irf,dataset_.vobs,M_.exo_nbr,MAX_nirfs_dsgevar); end end if irun == MAX_nirfs_dsge || irun == B || fpar == B @@ -279,20 +279,6 @@ while fpar<B ifil2 = ifil2 + 1; irun2 = 0; end -% if isoctave -% if (whoiam==0), -% printf(['Posterior IRF %3.f%% done\r'],(fpar/B*100)); -% end -% elseif ~whoiam, -% waitbar(fpar/B,h); -% end -% if whoiam, -% if ~isoctave -% fprintf('Done! \n'); -% end -% waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(B) ' done.']; -% fMessageStatus((fpar-fpar0)/(B-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); -% end dyn_waitbar((fpar-fpar0)/(B-fpar0),h); end @@ -310,9 +296,5 @@ end myoutput.OutputFileName = [OutputFileName_dsge; OutputFileName_param; OutputFileName_bvardsge]; - -myoutput.nosaddle = nosaddle; - - - +myoutput.nosaddle = nosaddle; diff --git a/matlab/PosteriorIRF_core2.m b/matlab/PosteriorIRF_core2.m index f33fbc2832ac0529a3146fe4bed8713adc552b90..c20702eb5f4658996bb7aa9a0690a6e49ba7724d 100644 --- a/matlab/PosteriorIRF_core2.m +++ b/matlab/PosteriorIRF_core2.m @@ -5,7 +5,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) % Perform in parallel execution a portion of the PosteriorIRF.m code. % See also the comment in random_walk_metropolis_hastings_core.m funtion. % -% INPUTS +% INPUTS % See the comment in random_walk_metropolis_hastings_core.m funtion. % % OUTPUTS @@ -13,8 +13,8 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) % Contained: % OutputFileName (i.e. the figures without the file .txt). % -% ALGORITHM -% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle. +% ALGORITHM +% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle. % % SPECIAL REQUIREMENTS. % None. @@ -36,7 +36,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see <http://www.gnu.org/licenses/>. -global options_ M_ +global options_ M_ if nargin<4, whoiam=0; @@ -87,7 +87,7 @@ OutputFileName={}; subplotnum = 0; for i=fpar:npar, figunumber = 0; - + for j=1:nvar if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold subplotnum = subplotnum+1; @@ -96,7 +96,7 @@ for i=fpar:npar, elseif subplotnum == 1 && ~options_.relative_irf hh = dyn_figure(options_,'Name',['Orthogonalized shock to ' tit(i,:)]); end - + set(0,'CurrentFigure',hh) subplot(nn,nn,subplotnum); if ~MAX_nirfs_dsgevar @@ -108,12 +108,12 @@ for i=fpar:npar, set(h2,'FaceColor',[1 1 1]); set(h2,'BaseValue',min(HPDIRF(:,1,j,i))); plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); box on axis tight xlim([1 options_.irf]); hold off - else + else h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); set(h1,'FaceColor',[.9 .9 .9]); set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); @@ -136,9 +136,9 @@ for i=fpar:npar, else if options_.debug fprintf('POSTERIOR_IRF: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(varlist(j,:)),tit(i,:),options_.impulse_responses.plot_threshold) - end + end end - + if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0) figunumber = figunumber+1; dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_); @@ -154,9 +154,8 @@ for i=fpar:npar, % fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); dyn_waitbar((i-fpar+1)/(npar-fpar+1),[],waitbarString); end -end% loop over exo_var +end% loop over exo_var myoutput.OutputFileName = OutputFileName; -