diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index 61f7361d9..c07fc9231 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -1,10 +1,18 @@ -function PosteriorIRF(type,dispString) +function oo_=PosteriorIRF(type,options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString) +% PosteriorIRF(type,options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString) % Builds posterior IRFs after the MH algorithm. % % INPUTS -% o type [char] string specifying the joint density of the -% deep parameters ('prior','posterior'). -% o dispString [char] string to display in the console. +% o type [char] string specifying the joint density of the +% deep parameters ('prior','posterior'). +% o options_ [structure] storing the options +% o estim_params_ [structure] storing information about estimated parameters +% o oo_ [structure] storing the results +% o M_ [structure] storing the model information +% o bayestopt_ [structure] storing information about priors +% o dataset_ [structure] storing the dataset +% o dataset_info [structure] Various information about the dataset +% o dispString [char] string to display in the console. % % OUTPUTS % None (oo_ and plots). @@ -34,9 +42,6 @@ function PosteriorIRF(type,dispString) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info - % Set the number of periods if isempty(options_.irf) || ~options_.irf options_.irf = 40; @@ -62,13 +67,7 @@ end irf_shocks_indx = getIrfShocksIndx(M_, options_); % 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; -offset = npar-np; clear('nvx','nvn','ncx','ncn','np'); +npar = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np ; nvobs = dataset_.vobs; gend = dataset_.nobs; @@ -118,7 +117,8 @@ elseif strcmpi(type,'gsa') end x=[lpmat0(istable,:) lpmat(istable,:)]; clear lpmat istable - B=size(x,1); options_.B = B; + B=size(x,1); + options_.B = B; else% type = 'prior' B = options_.prior_draws; options_.B = B; @@ -130,23 +130,7 @@ irun2 = 0; NumberOfIRFfiles_dsge = 1; NumberOfIRFfiles_dsgevar = 1; ifil2 = 1; -% Create arrays -if B <= MAX_nruns - stock_param = zeros(B, npar); -else - stock_param = zeros(MAX_nruns, npar); -end -if B >= MAX_nirfs_dsge - stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge); -else - stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B); -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 NumberOfLags = options_.dsge_varlag; NumberOfLagsTimesNvobs = NumberOfLags*nvobs; if options_.noconstant @@ -154,10 +138,9 @@ if MAX_nirfs_dsgevar else NumberOfParametersPerEquation = NumberOfLagsTimesNvobs+1; end - Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs); end -% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m +% First block of code executed in parallel. The function devoted to do it is PosteriorIRF_core1.m % function. b = 0; @@ -183,7 +166,6 @@ if ~strcmpi(type,'prior') localVars.x=x; end -b=0; if options_.dsge_var localVars.gend = gend; localVars.nvobs = nvobs; @@ -202,6 +184,16 @@ localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; localVars.ifil2=ifil2; localVars.MhDirectoryName=MhDirectoryName; +%store main structures +localVars.options_=options_; +localVars.estim_params_= estim_params_; +localVars.M_= M_; +localVars.oo_= oo_; +localVars.bayestopt_= bayestopt_; +localVars.dataset_= dataset_; +localVars.dataset_info= dataset_info; + + % Like sequential execution! if isnumeric(options_.parallel) [fout] = PosteriorIRF_core1(localVars,1,B,0); @@ -225,14 +217,6 @@ else 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_info',dataset_info); - % which files have to be copied to run remotely NamFileInput(1,:) = {'',[M_.fname '.static.m']}; NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']}; @@ -246,7 +230,7 @@ else NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '.steadystate.m']}; end end - [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info); + [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, [], options_.parallel_info); nosaddle=0; for j=1:length(fout) nosaddle = nosaddle + fout(j).nosaddle; @@ -260,9 +244,9 @@ if nosaddle disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) end -ReshapeMatFiles('irf_dsge',type) +ReshapeMatFiles(M_.fname,M_.dname,M_.exo_nbr,M_.endo_nbr,options_,'irf_dsge',type) if MAX_nirfs_dsgevar - ReshapeMatFiles('irf_bvardsge') + ReshapeMatFiles(M_.fname,M_.dname,M_.exo_nbr,M_.endo_nbr,options_,'irf_bvardsge') end if strcmpi(type,'gsa') @@ -293,21 +277,21 @@ tit = M_.exo_names; kdx = 0; for file = 1:NumberOfIRFfiles_dsge - load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']); + temp=load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']); for i = irf_shocks_indx for j = 1:nvar - for k = 1:size(STOCK_IRF_DSGE,1) + for k = 1:size(temp.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); + DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(temp.STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig); end end end - kdx = kdx + size(STOCK_IRF_DSGE,1); + kdx = kdx + size(temp.STOCK_IRF_DSGE,1); end -clear STOCK_IRF_DSGE; +clear temp; for i = irf_shocks_indx for j = 1:nvar @@ -332,20 +316,20 @@ if MAX_nirfs_dsgevar tit = M_.exo_names; kdx = 0; for file = 1:NumberOfIRFfiles_dsgevar - load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']); + temp=load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']); for i = irf_shocks_indx for j = 1:nvar - for k = 1:size(STOCK_IRF_BVARDSGE,1) + for k = 1:size(temp.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); + posterior_moments(squeeze(temp.STOCK_IRF_BVARDSGE(k,j,i,:)),0,options_.mh_conf_sig); end end end - kdx = kdx + size(STOCK_IRF_BVARDSGE,1); + kdx = kdx + size(temp.STOCK_IRF_BVARDSGE,1); end - clear STOCK_IRF_BVARDSGE; + clear temp; for i = irf_shocks_indx for j = 1:nvar name = sprintf('%s_%s', M_.endo_names{IndxVariables(j)}, tit{i}); @@ -358,25 +342,24 @@ if MAX_nirfs_dsgevar end end end -% -% Finally I build the plots. -% +%% Finally I build the plots. % 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 +% .tex generation always run in sequentially. This portion of code is executed in parallel by % PosteriorIRF_core2.m function. if ~options_.nograph && ~options_.no_graph.posterior % Save the local variables. localVars=[]; - + localVars.dname=M_.dname; + localVars.fname=M_.fname; + localVars.options_=options_; Check=options_.TeX; if (Check) localVars.varlist_TeX=varlist_TeX; end - localVars.nvar=nvar; localVars.MeanIRF=MeanIRF; localVars.tit=tit; @@ -445,9 +428,7 @@ if ~options_.nograph && ~options_.no_graph.posterior if isRemoteOctave [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); else - globalVars = struct('M_',M_, ... - 'options_', options_); - + globalVars = []; [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info); end end @@ -455,7 +436,6 @@ if ~options_.nograph && ~options_.no_graph.posterior [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); end % END parallel code! - end -fprintf('%s: Posterior IRFs, done!\n',dispString); +fprintf('%s: Posterior IRFs, done!\n',dispString); \ No newline at end of file diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index a0a47dd5a..2eac57416 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -1,4 +1,5 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) +%myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) % Generates and stores Posterior IRFs % PARALLEL CONTEXT % This function perfoms in parallel execution a portion of the PosteriorIRF.m code. @@ -40,9 +41,6 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - -global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info - if nargin<4 whoiam=0; end @@ -50,6 +48,14 @@ end % Reshape 'myinputs' for local computation. % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: +options_= myinputs.options_; +estim_params_= myinputs.estim_params_; +M_= myinputs.M_; +oo_= myinputs.oo_; +bayestopt_= myinputs.bayestopt_; +dataset_= myinputs.dataset_; +dataset_info= myinputs.dataset_info; + IRUN = myinputs.IRUN; irun =myinputs.irun; irun2=myinputs.irun2; @@ -78,13 +84,10 @@ if options_.dsge_var bounds = prior_bounds(bayestopt_,options_.prior_trunc); end - if whoiam Parallel=myinputs.Parallel; end - -% MhDirectoryName = myinputs.MhDirectoryName; if strcmpi(type,'posterior') MhDirectoryName = CheckPath('metropolis',M_.dname); elseif strcmpi(type,'gsa') @@ -115,12 +118,10 @@ else h = dyn_waitbar(prct0,'Bayesian (prior) IRFs...'); end - OutputFileName_bvardsge = {}; OutputFileName_dsge = {}; OutputFileName_param = {}; - fpar = fpar-1; fpar0=fpar; nosaddle=0; diff --git a/matlab/PosteriorIRF_core2.m b/matlab/PosteriorIRF_core2.m index 73ea4c493..09b8650d9 100644 --- a/matlab/PosteriorIRF_core2.m +++ b/matlab/PosteriorIRF_core2.m @@ -1,5 +1,5 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam,ThisMatlab) -% function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) +% myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) % Generates the Posterior IRFs plot from the IRFs generated in % PosteriorIRF_core1 % @@ -47,8 +47,6 @@ 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 . -global options_ M_ - if nargin<4 whoiam=0; end @@ -56,6 +54,7 @@ end % Reshape 'myinputs' for local computation. % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: +options_=myinputs.options_; Check=options_.TeX; if (Check) varlist_TeX=myinputs.varlist_TeX; @@ -67,6 +66,9 @@ tit=myinputs.tit; nn=myinputs.nn; MAX_nirfs_dsgevar=myinputs.MAX_nirfs_dsgevar; HPDIRF=myinputs.HPDIRF; +dname=myinputs.dname; +fname=myinputs.fname; + if options_.dsge_var HPDIRFdsgevar=myinputs.HPDIRFdsgevar; MeanIRFdsgevar=myinputs.MeanIRFdsgevar; @@ -82,7 +84,7 @@ end % To save the figures where the function is computed! -DirectoryName = CheckPath('Output',M_.dname); +DirectoryName = CheckPath('Output',dname); RemoteFlag = 0; if whoiam @@ -117,7 +119,6 @@ for i=fpar:npar h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i))); %white below HPDIinf and minimum of HPDIinf end plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); box on axis tight xlim([1 options_.irf]); @@ -135,7 +136,6 @@ for i=fpar:npar plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1) plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1) end - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); box on axis tight xlim([1 options_.irf]); @@ -155,9 +155,9 @@ for i=fpar:npar if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0) figunumber = figunumber+1; - dyn_saveas(hh_fig,[DirectoryName '/' M_.fname '_Bayesian_IRF_' tit{i} '_' int2str(figunumber)],options_.nodisplay,options_.graph_format); + dyn_saveas(hh_fig,[DirectoryName '/' fname '_Bayesian_IRF_' tit{i} '_' int2str(figunumber)],options_.nodisplay,options_.graph_format); if RemoteFlag==1 - OutputFileName = [OutputFileName; {[DirectoryName,filesep], [M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}]; + OutputFileName = [OutputFileName; {[DirectoryName,filesep], [fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}]; end subplotnum = 0; end @@ -165,7 +165,6 @@ for i=fpar:npar if whoiam fprintf('Done! \n'); waitbarString = [ 'Exog. shocks ' int2str(i) '/' int2str(npar) ' done.']; - % 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 diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m index 38e529e59..58107602e 100644 --- a/matlab/ReshapeMatFiles.m +++ b/matlab/ReshapeMatFiles.m @@ -1,10 +1,15 @@ -function ReshapeMatFiles(type, type2) -% function ReshapeMatFiles(type, type2) +function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2) +% function ReshapeMatFiles(fname, dname, exo_nbr, endo_nbr, options_, type, type2) % Reshapes and sorts (along the mcmc simulations) the mat files generated by DYNARE. % 4D-arrays are splitted along the first dimension. % 3D-arrays are splitted along the second dimension. % % INPUTS: +% fname: [string] filename +% dname: [string] directory name +% exo_nbr: [integer] number of exogenous variables +% endo_nbr: [integer] number of endogenous variables +% options_: [struct] options structure % type: statistics type in the repertory: % dgse % irf_bvardsge @@ -25,7 +30,7 @@ function ReshapeMatFiles(type, type2) % SPECIAL REQUIREMENTS % none -% Copyright © 2003-2017 Dynare Team +% Copyright © 2003-2023 Dynare Team % % This file is part of Dynare. % @@ -42,43 +47,41 @@ function ReshapeMatFiles(type, type2) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ options_ - -if nargin==1 - MhDirectoryName = [ CheckPath('metropolis',M_.dname) filesep ]; +if nargin==6 + MhDirectoryName = [ CheckPath('metropolis',dname) filesep ]; else if strcmpi(type2,'posterior') - MhDirectoryName = [CheckPath('metropolis',M_.dname) filesep ]; + MhDirectoryName = [CheckPath('metropolis',dname) filesep ]; elseif strcmpi(type2,'gsa') if options_.opt_gsa.morris==1 - MhDirectoryName = [CheckPath('gsa/screen',M_.dname) filesep ]; + MhDirectoryName = [CheckPath('gsa/screen',dname) filesep ]; elseif options_.opt_gsa.morris==2 - MhDirectoryName = [CheckPath('gsa/identif',M_.dname) filesep ]; + MhDirectoryName = [CheckPath('gsa/identif',dname) filesep ]; elseif options_.opt_gsa.pprior - MhDirectoryName = [CheckPath(['gsa' filesep 'prior'],M_.dname) filesep ]; + MhDirectoryName = [CheckPath(['gsa' filesep 'prior'],dname) filesep ]; else - MhDirectoryName = [CheckPath(['gsa' filesep 'mc'],M_.dname) filesep ]; + MhDirectoryName = [CheckPath(['gsa' filesep 'mc'],dname) filesep ]; end else - MhDirectoryName = [CheckPath('prior',M_.dname) filesep ]; + MhDirectoryName = [CheckPath('prior',dname) filesep ]; end end switch type case 'irf_dsge' CAPtype = 'IRF_DSGE'; - TYPEsize = [ options_.irf , size(options_.varlist,1) , M_.exo_nbr ]; + TYPEsize = [ options_.irf , size(options_.varlist,1) , exo_nbr ]; TYPEarray = 4; case 'irf_bvardsge' CAPtype = 'IRF_BVARDSGE'; - TYPEsize = [ options_.irf , length(options_.varobs) , M_.exo_nbr ]; + TYPEsize = [ options_.irf , length(options_.varobs) , exo_nbr ]; TYPEarray = 4; case 'smooth' CAPtype = 'SMOOTH'; - TYPEsize = [ M_.endo_nbr , options_.nobs ]; + TYPEsize = [ endo_nbr , options_.nobs ]; TYPEarray = 3; case 'filter' CAPtype = 'FILTER'; - TYPEsize = [ M_.endo_nbr , options_.nobs + 1 ];% TO BE CHECKED! + TYPEsize = [ endo_nbr , options_.nobs + 1 ];% TO BE CHECKED! TYPEarray = 3; case 'error' CAPtype = 'ERROR'; @@ -86,22 +89,22 @@ switch type TYPEarray = 3; case 'innov' CAPtype = 'INNOV'; - TYPEsize = [ M_.exo_nbr , options_.nobs ]; + TYPEsize = [ exo_nbr , options_.nobs ]; TYPEarray = 3; case 'forcst' CAPtype = 'FORCST'; - TYPEsize = [ M_.endo_nbr , options_.forecast ]; + TYPEsize = [ endo_nbr , options_.forecast ]; TYPEarray = 3; case 'forcst1' CAPtype = 'FORCST1'; - TYPEsize = [ M_.endo_nbr , options_.forecast ]; + TYPEsize = [ endo_nbr , options_.forecast ]; TYPEarray = 3; otherwise disp('ReshapeMatFiles :: Unknown argument!') return end -TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*.mat']); +TYPEfiles = dir([MhDirectoryName fname '_' type '*.mat']); NumberOfTYPEfiles = length(TYPEfiles); B = options_.B; @@ -116,36 +119,29 @@ switch TYPEarray for f1=1:NumberOfTYPEfiles-foffset eval(['STOCK_' CAPtype ' = zeros(NumberOfPeriodsPerTYPEfiles,TYPEsize(2),TYPEsize(3),B);']) for f2 = 1:NumberOfTYPEfiles - load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']); + load([MhDirectoryName fname '_' type int2str(f2) '.mat']); eval(['STOCK_' CAPtype '(:,:,1:+size(stock_' type ',3),idx+1:idx+size(stock_' type ',4))=stock_' ... type '(jdx+1:jdx+NumberOfPeriodsPerTYPEfiles,:,:,:);']) eval(['idx = idx + size(stock_' type ',4);']) end - %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',4);']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]); jdx = jdx + NumberOfPeriodsPerTYPEfiles; idx = 0; end if reste eval(['STOCK_' CAPtype ' = zeros(reste,TYPEsize(2),TYPEsize(3),B);']) for f2 = 1:NumberOfTYPEfiles - load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']); + load([MhDirectoryName fname '_' type int2str(f2) '.mat']); eval(['STOCK_' CAPtype '(:,:,:,idx+1:idx+size(stock_' type ',4))=stock_' type '(jdx+1:jdx+reste,:,:,:);']) eval(['idx = idx + size(stock_' type ',4);']) end - %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',4);']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(NumberOfTYPEfiles-foffset+1) '.mat'],['STOCK_' CAPtype]); end else - load([MhDirectoryName M_.fname '_' type '1.mat']); - %eval(['STOCK_' CAPtype ' = sort(stock_' type ',4);']) + load([MhDirectoryName fname '_' type '1.mat']); eval(['STOCK_' CAPtype ' = stock_' type ';']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]); end - % Original file format may be useful in some cases... - % for file = 1:NumberOfTYPEfiles - % delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat']) - % end case 3 if NumberOfTYPEfiles>1 NumberOfPeriodsPerTYPEfiles = ceil( TYPEsize(2)/NumberOfTYPEfiles ); @@ -155,31 +151,24 @@ switch TYPEarray for f1=1:NumberOfTYPEfiles-1 eval(['STOCK_' CAPtype ' = zeros(TYPEsize(1),NumberOfPeriodsPerTYPEfiles,B);']) for f2 = 1:NumberOfTYPEfiles - load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']); + load([MhDirectoryName fname '_' type int2str(f2) '.mat']); eval(['STOCK_' CAPtype '(:,:,idx+1:idx+size(stock_ ' type ',3))=stock_' type '(:,jdx+1:jdx+NumberOfPeriodsPerTYPEfiles,:);']) eval(['idx = idx + size(stock_' type ',3);']) end - %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',3);']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(f1) '.mat'],['STOCK_' CAPtype]); jdx = jdx + NumberOfPeriodsPerTYPEfiles; idx = 0; end eval(['STOCK_' CAPtype ' = zeros(TYPEsize(1),reste,B);']) for f2 = 1:NumberOfTYPEfiles - load([MhDirectoryName M_.fname '_' type int2str(f2) '.mat']); + load([MhDirectoryName fname '_' type int2str(f2) '.mat']); eval(['STOCK_' CAPtype '(:,:,idx+1:idx+size(stock_' type ',3))=stock_' type '(:,jdx+1:jdx+reste,:);']) eval(['idx = idx + size(stock_' type ',3);']) end - %eval(['STOCK_' CAPtype ' = sort(STOCK_' CAPtype ',3);']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(NumberOfTYPEfiles) '.mat'],['STOCK_' CAPtype]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(NumberOfTYPEfiles) '.mat'],['STOCK_' CAPtype]); else - load([MhDirectoryName M_.fname '_' type '1.mat']); - %eval(['STOCK_' CAPtype ' = sort(stock_' type ',3);']) + load([MhDirectoryName fname '_' type '1.mat']); eval(['STOCK_' CAPtype ' = stock_' type ';']) - save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]); + save([MhDirectoryName fname '_' CAPtype 's' int2str(1) '.mat'],['STOCK_' CAPtype ]); end - % Original file format may be useful in some cases... - % for file = 1:NumberOfTYPEfiles - % delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat']) - % end end \ No newline at end of file diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 4c63a2ccf..94797144e 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -225,70 +225,70 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation current_optimizer = optimizer_vec{optim_iter}; [xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); - fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval); + fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval); if isnumeric(current_optimizer) if current_optimizer==5 - newratflag = new_rat_hess_info.newratflag; - new_rat_hess_info = new_rat_hess_info.new_rat_hess_info; + newratflag = new_rat_hess_info.newratflag; + new_rat_hess_info = new_rat_hess_info.new_rat_hess_info; elseif current_optimizer==6 %save scaling factor - save([M_.dname filesep 'Output' filesep M_.fname '_optimal_mh_scale_parameter.mat'],'Scale'); - options_.mh_jscale = Scale; - bayestopt_.jscale(:) = options_.mh_jscale; - end + save([M_.dname filesep 'Output' filesep M_.fname '_optimal_mh_scale_parameter.mat'],'Scale'); + options_.mh_jscale = Scale; + bayestopt_.jscale(:) = options_.mh_jscale; + end end if ~isnumeric(current_optimizer) || ~isequal(current_optimizer,6) %always already computes covariance matrix - if options_.cova_compute == 1 %user did not request covariance not to be computed - if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood') - ana_deriv_old = options_.analytic_derivation; - options_.analytic_derivation = 2; - [~,~,~,~,hh] = feval(objective_function,xparam1, ... - dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); - options_.analytic_derivation = ana_deriv_old; + if options_.cova_compute == 1 %user did not request covariance not to be computed + if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood') + ana_deriv_old = options_.analytic_derivation; + options_.analytic_derivation = 2; + [~,~,~,~,hh] = feval(objective_function,xparam1, ... + dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); + options_.analytic_derivation = ana_deriv_old; elseif ~isnumeric(current_optimizer) || ~(isequal(current_optimizer,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood')) - % enter here if i) not mode_compute_5, ii) if mode_compute_5 and newratflag==1; - % with flag==0 or 2 and dsge_likelihood, we force to use - % the hessian from outer product gradient of optimizer 5 below - if options_.hessian.use_penalized_objective - penalized_objective_function = str2func('penalty_objective_function'); - hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); - else - hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); - end - hh = reshape(hh, nx, nx); - elseif isnumeric(current_optimizer) && isequal(current_optimizer,5) - % other numerical hessian options available with optimizer - % 5 and dsge_likelihood - % - % if newratflag == 0 - % compute outer product gradient of optimizer 5 - % - % if newratflag == 2 - % compute 'mixed' outer product gradient of optimizer 5 - % with diagonal elements computed with numerical second order derivatives - % - % uses univariate filters, so to get max # of available - % densities for outer product gradient - kalman_algo0 = options_.kalman_algo; - compute_hessian = 1; - if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4)) - options_.kalman_algo=2; - if options_.lik_init == 3 - options_.kalman_algo=4; - end - elseif newratflag==0 % hh already contains outer product gradient with univariate filter - compute_hessian = 0; - end - if compute_hessian - crit = options_.newrat.tolerance.f; - newratflag = newratflag>0; - hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx); - end - options_.kalman_algo = kalman_algo0; + % enter here if i) not mode_compute_5, ii) if mode_compute_5 and newratflag==1; + % with flag==0 or 2 and dsge_likelihood, we force to use + % the hessian from outer product gradient of optimizer 5 below + if options_.hessian.use_penalized_objective + penalized_objective_function = str2func('penalty_objective_function'); + hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); + else + hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); end + hh = reshape(hh, nx, nx); + elseif isnumeric(current_optimizer) && isequal(current_optimizer,5) + % other numerical hessian options available with optimizer + % 5 and dsge_likelihood + % + % if newratflag == 0 + % compute outer product gradient of optimizer 5 + % + % if newratflag == 2 + % compute 'mixed' outer product gradient of optimizer 5 + % with diagonal elements computed with numerical second order derivatives + % + % uses univariate filters, so to get max # of available + % densities for outer product gradient + kalman_algo0 = options_.kalman_algo; + compute_hessian = 1; + if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4)) + options_.kalman_algo=2; + if options_.lik_init == 3 + options_.kalman_algo=4; + end + elseif newratflag==0 % hh already contains outer product gradient with univariate filter + compute_hessian = 0; + end + if compute_hessian + crit = options_.newrat.tolerance.f; + newratflag = newratflag>0; + hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,0,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx); + end + options_.kalman_algo = kalman_algo0; end end - parameter_names = bayestopt_.name; + end + parameter_names = bayestopt_.name; end if options_.cova_compute || current_optimizer==5 || current_optimizer==6 save([M_.dname filesep 'Output' filesep M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval'); @@ -472,7 +472,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... if error_flag error('%s: I cannot compute the posterior IRFs!',dispString) end - PosteriorIRF('posterior',dispString); + oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString); end if options_.moments_varendo if error_flag @@ -504,7 +504,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... error('%s: I cannot compute the posterior statistics!',dispString) end if options_.order==1 && ~options_.particle.status - prior_posterior_statistics('posterior',dataset_,dataset_info,dispString); %get smoothed and filtered objects and forecasts + oo_=prior_posterior_statistics('posterior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString); %get smoothed and filtered objects and forecasts else error('%s: Particle Smoothers are not yet implemented.',dispString) end diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m index 4a3389d22..9adfc643a 100644 --- a/matlab/dynare_sensitivity.m +++ b/matlab/dynare_sensitivity.m @@ -410,9 +410,9 @@ if options_gsa.rmse end end - prior_posterior_statistics('gsa',dataset_, dataset_info,'gsa::mcmc'); + oo_=prior_posterior_statistics('gsa',dataset_, dataset_info,M_,oo_,options_,estim_params_,bayestopt_,'gsa::mcmc'); if options_.bayesian_irf - PosteriorIRF('gsa','gsa::mcmc'); + oo_=PosteriorIRF('gsa',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,'gsa::mcmc'); end options_gsa.load_rmse=0; % else diff --git a/matlab/pm3.m b/matlab/pm3.m index 3800cd515..d08b06830 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -1,4 +1,4 @@ -function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryName,var_type,dispString) +function oo_=pm3(M_,options_,oo_,n1,n2,ifil,B,tit1,tit2,tit_tex,names1,names2,name3,DirectoryName,var_type,dispString) % Computes, stores and plots the posterior moment statistics. % @@ -8,8 +8,7 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa % ifil [scalar] number of moment files to load % B [scalar] number of subdraws % tit1 [string] Figure title -% tit2 [string] not used -% tit3 [string] Save name for figure +% tit2 [string] Save name for figure % tit_tex [cell array] TeX-Names for Variables % names1 [cell array] Names of all variables in the moment matrix from % which names2 is selected @@ -20,6 +19,10 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa % var_type [string] suffix of the filename from which to load moment % matrix % dispString [string] string to be displayes in the command window +% +% OUTPUTS +% oo_ [structure] storing the results + % PARALLEL CONTEXT % See also the comment in posterior_sampler.m funtion. @@ -42,8 +45,6 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ M_ oo_ - nn = 3; MaxNumberOfPlotsPerFigure = nn^2; % must be square varlist = names2; @@ -284,9 +285,7 @@ if strcmp(var_type,'_trend_coeff') || max(max(abs(Mean(:,:))))<=10^(-6) || all(a fprintf(['%s: ' tit1 ', done!\n'],dispString); return %not do plots end -%% %% Finally I build the plots. -%% if ~options_.nograph && ~options_.no_graph.posterior % Block of code executed in parallel, with the exception of file @@ -297,8 +296,6 @@ if ~options_.nograph && ~options_.no_graph.posterior % % %%% The file .TeX! are not saved in parallel. - - % Store the variable mandatory for local/remote parallel computing. localVars=[]; @@ -313,8 +310,14 @@ if ~options_.nograph && ~options_.no_graph.posterior end localVars.MaxNumberOfPlotsPerFigure=MaxNumberOfPlotsPerFigure; localVars.name3=name3; - localVars.tit3=tit3; + localVars.tit2=tit2; localVars.Mean=Mean; + localVars.TeX=options_.TeX; + localVars.nodisplay=options_.nodisplay; + localVars.graph_format=options_.graph_format; + localVars.dname=M_.dname; + localVars.fname=M_.fname; + % Like sequential execution! nvar0=nvar; @@ -332,16 +335,13 @@ if ~options_.nograph && ~options_.no_graph.posterior if isRemoteOctave fout = pm3_core(localVars,1,nvar,0); else - globalVars = struct('M_',M_, ... - 'options_', options_, ... - 'oo_', oo_); + globalVars = []; [fout, nvar0, totCPU] = masterParallel(options_.parallel, 1, nvar, [],'pm3_core', localVars,globalVars, options_.parallel_info); end end else % For the time being in Octave enviroment the pm3.m is executed only in % serial modality, to avoid problem with the plots. - fout = pm3_core(localVars,1,nvar,0); end @@ -365,8 +365,8 @@ if ~options_.nograph && ~options_.no_graph.posterior if subplotnum == MaxNumberOfPlotsPerFigure || i == nvar fprintf(fidTeX,'\\begin{figure}[H]\n'); fprintf(fidTeX,'\\centering \n'); - fprintf(fidTeX,['\\includegraphics[width=%2.2f\\textwidth]{%s/Output/%s_' name3 '_%s}\n'],options_.figures.textwidth*min(subplotnum/nn,1),M_.dname,M_.fname, tit3{i}); - fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,tit3{i}); + fprintf(fidTeX,['\\includegraphics[width=%2.2f\\textwidth]{%s/Output/%s_' name3 '_%s}\n'],options_.figures.textwidth*min(subplotnum/nn,1),M_.dname,M_.fname, tit2{i}); + fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,tit2{i}); fprintf(fidTeX,'\\caption{%s}\n',tit1); fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,' \n'); diff --git a/matlab/pm3_core.m b/matlab/pm3_core.m index 343c70897..f820ed233 100644 --- a/matlab/pm3_core.m +++ b/matlab/pm3_core.m @@ -1,13 +1,23 @@ function myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab) - +% myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab) % PARALLEL CONTEXT % Core functionality for pm3.m function, which can be parallelized. % INPUTS -% See the comment in posterior_sampler_core.m funtion. - +% o myimput [struc] The mandatory variables for local/remote +% parallel computing obtained from prior_posterior_statistics.m +% function. +% o fpar and nvar [integer] first variable and number of variables +% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed. +% The integer whoaim is the integer that +% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the +% cluster. +% o ThisMatlab [integer] Allows us to distinguish between the +% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab, +% ... Then it is the index number of this slave machine in the cluster. +% % OUTPUTS -% o myoutput [struc] +% o myoutput [struct] Contains file names % % % SPECIAL REQUIREMENTS. @@ -45,16 +55,18 @@ varlist=myinputs.varlist; MaxNumberOfPlotsPerFigure=myinputs.MaxNumberOfPlotsPerFigure; name3=myinputs.name3; -tit3=myinputs.tit3; +tit2=myinputs.tit2; Mean=myinputs.Mean; +options_.TeX=myinputs.TeX; +options_.nodisplay=myinputs.nodisplay; +options_.graph_format=myinputs.graph_format; +fname=myinputs.fname; +dname=myinputs.dname; if whoiam Parallel=myinputs.Parallel; end - -global options_ M_ oo_ - if options_.TeX varlist_TeX=myinputs.varlist_TeX; end @@ -64,8 +76,6 @@ if whoiam h = dyn_waitbar(prct0,'Parallel plots pm3 ...'); end - - figunumber = 0; subplotnum = 0; hh_fig = dyn_figure(options_.nodisplay,'Name',[tit1 ' ' int2str(figunumber+1)]); @@ -110,14 +120,14 @@ for i=fpar:nvar if whoiam if Parallel(ThisMatlab).Local==0 - DirectoryName = CheckPath('Output',M_.dname); + DirectoryName = CheckPath('Output',dname); end end if subplotnum == MaxNumberOfPlotsPerFigure || i == nvar - dyn_saveas(hh_fig,[M_.dname '/Output/' M_.fname '_' name3 '_' tit3{i}],options_.nodisplay,options_.graph_format); + dyn_saveas(hh_fig,[dname '/Output/' fname '_' name3 '_' tit2{i}],options_.nodisplay,options_.graph_format); if RemoteFlag==1 - OutputFileName = [OutputFileName; {[M_.dname, filesep, 'Output',filesep], [M_.fname '_' name3 '_' deblank(tit3(i,:)) '.*']}]; + OutputFileName = [OutputFileName; {[dname, filesep, 'Output',filesep], [fname '_' name3 '_' deblank(tit2(i,:)) '.*']}]; end subplotnum = 0; figunumber = figunumber+1; @@ -127,12 +137,8 @@ for i=fpar:nvar end if whoiam - % waitbarString = [ 'Variable ' int2str(i) '/' int2str(nvar) ' done.']; - % fMessageStatus((i-fpar+1)/(nvar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab)); dyn_waitbar((i-fpar+1)/(nvar-fpar+1),h); end - - end if whoiam diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 4bdddfd2d..58f428c5a 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -1,17 +1,19 @@ -function prior_posterior_statistics(type,dataset,dataset_info,dispString) -% function prior_posterior_statistics(type,dataset,dataset_info,dispString) +function oo_=prior_posterior_statistics(type,dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString) +% oo_=prior_posterior_statistics(type,dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_,dispString)) % Computes Monte Carlo filter smoother and forecasts % % INPUTS -% type: posterior -% prior -% gsa -% dataset: data structure -% dataset_info: dataset structure -% dispString: string to display in the command window -% +% type [string] posterior, prior, or gsa +% o dataset_ [structure] storing the dataset +% o dataset_info [structure] Various information about the dataset +% o M_ [structure] storing the model information +% o oo_ [structure] storing the results +% o options_ [structure] storing the options +% o estim_params_ [structure] storing information about estimated parameters +% o bayestopt_ [structure] storing information about priors +% dispString: [string] display info in the command window % OUTPUTS -% none +% oo_: [structure] storing the results % % SPECIAL REQUIREMENTS % none @@ -37,36 +39,23 @@ function prior_posterior_statistics(type,dataset,dataset_info,dispString) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -if nargin < 4 +if nargin < 9 dispString = 'prior_posterior_statistics'; end - -global options_ estim_params_ oo_ M_ bayestopt_ - localVars=[]; -Y = transpose(dataset.data); -gend = dataset.nobs; -data_index = dataset_info.missing.aindex; -missing_value = dataset_info.missing.state; -mean_varobs = dataset_info.descriptive.mean; +Y = transpose(dataset_.data); +gend = dataset_.nobs; - -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; +npar = estim_params_.nvx+nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np; naK = length(options_.filter_step_ahead); MaxNumberOfBytes=options_.MaxNumberOfBytes; endo_nbr=M_.endo_nbr; exo_nbr=M_.exo_nbr; meas_err_nbr=length(M_.Correlation_matrix_ME); -iendo = 1:endo_nbr; horizon = options_.forecast; -IdObs = bayestopt_.mfys; if horizon i_last_obs = gend+(1-M_.maximum_endo_lag:0); end @@ -130,13 +119,6 @@ varlist = options_.varlist; if isempty(varlist) varlist = sort(M_.endo_names(1:M_.orig_endo_nbr)); end -nvar = length(varlist); -SelecVariables = []; -for i=1:nvar - if ~isempty(strmatch(varlist{i}, M_.endo_names, 'exact')) - SelecVariables = [SelecVariables; strmatch(varlist{i}, M_.endo_names, 'exact')]; - end -end n_variables_to_fill=13; @@ -171,17 +153,17 @@ localVars.filter_covariance=filter_covariance; localVars.smoothed_state_uncertainty=smoothed_state_uncertainty; localVars.gend=gend; localVars.Y=Y; -localVars.data_index=data_index; -localVars.missing_value=missing_value; +localVars.data_index=dataset_info.missing.aindex; +localVars.missing_value=dataset_info.missing.state; localVars.varobs=options_.varobs; -localVars.mean_varobs=mean_varobs; +localVars.mean_varobs=dataset_info.descriptive.mean; localVars.irun=irun; localVars.endo_nbr=endo_nbr; localVars.nvn=nvn; localVars.naK=naK; localVars.horizon=horizon; -localVars.iendo=iendo; -localVars.IdObs=IdObs; +localVars.iendo=1:endo_nbr; +localVars.IdObs=bayestopt_.mfys; if horizon localVars.i_last_obs=i_last_obs; localVars.MAX_nforc1=MAX_nforc1; @@ -212,6 +194,13 @@ localVars.MAX_momentsno = MAX_momentsno; localVars.ifil=ifil; localVars.DirectoryName = DirectoryName; +localVars.M_=M_; +localVars.oo_=oo_; +localVars.options_=options_; +localVars.estim_params_=estim_params_; +localVars.bayestopt_=bayestopt_; + + if strcmpi(type,'posterior') record=load_last_mh_history_file(DirectoryName, M_.fname); [nblck, npar] = size(record.LastParameters); @@ -290,11 +279,7 @@ else end end localVars.ifil = ifil; - globalVars = struct('M_',M_, ... - 'options_', options_, ... - 'bayestopt_', bayestopt_, ... - 'estim_params_', estim_params_, ... - 'oo_', oo_); + globalVars = []; % which files have to be copied to run remotely NamFileInput(1,:) = {'',[M_.fname '.static.m']}; NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']}; @@ -330,28 +315,28 @@ if ~isnumeric(options_.parallel) end if options_.smoother - pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',... - '',varlist, M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(1),B,'Smoothed variables',... + varlist, M_.endo_names_tex,M_.endo_names,... varlist,'SmoothedVariables',DirectoryName,'_smooth',dispString); - pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',... - '',M_.exo_names,M_.exo_names_tex,M_.exo_names,... + oo_=pm3(M_,options_,oo_,exo_nbr,gend,ifil(2),B,'Smoothed shocks',... + M_.exo_names,M_.exo_names_tex,M_.exo_names,... M_.exo_names,'SmoothedShocks',DirectoryName,'_inno',dispString); - pm3(endo_nbr,1,ifil(9),B,'Trend_coefficients',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,1,ifil(9),B,'Trend_coefficients',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'TrendCoeff',DirectoryName,'_trend_coeff',dispString); - pm3(endo_nbr,gend,ifil(10),B,'Smoothed constant',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(10),B,'Smoothed constant',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'Constant',DirectoryName,'_smoothed_constant',dispString); - pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(11),B,'Smoothed trend',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'Trend',DirectoryName,'_smoothed_trend',dispString); - pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(1),B,'Updated Variables',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'UpdatedVariables',DirectoryName, ... '_update',dispString); if smoothed_state_uncertainty - pm3(endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'StateUncertainty',DirectoryName,'_state_uncert',dispString); end @@ -360,24 +345,24 @@ if options_.smoother meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}]; texnames{obs_iter,1}=['\sigma^{ME}_' M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}]; end - pm3(meas_err_nbr,gend,ifil(3),B,'Smoothed measurement errors',... - '',meas_error_names,texnames,meas_error_names,... - meas_error_names,'SmoothedMeasurementErrors',DirectoryName,'_error',dispString) + oo_=pm3(M_,options_,oo_,meas_err_nbr,gend,ifil(3),B,'Smoothed measurement errors',... + meas_error_names,texnames,meas_error_names,... + meas_error_names,'SmoothedMeasurementErrors',DirectoryName,'_error',dispString); end end if options_.filtered_vars - pm3(endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead',dispString); end if options_.forecast - pm3(endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'MeanForecast',DirectoryName,'_forc_mean',dispString); - pm3(endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'PointForecast',DirectoryName,'_forc_point',dispString); if ~isequal(M_.H,0) && ~isempty(intersect(options_.varobs,varlist)) texnames = cell(length(options_.varobs), 1); @@ -387,15 +372,15 @@ if options_.forecast texnames{obs_iter}=M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}; end varlist_forecast_ME=intersect(options_.varobs,varlist); - pm3(meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',... - '',varlist_forecast_ME,texnames,obs_names,... - varlist_forecast_ME,'PointForecastME',DirectoryName,'_forc_point_ME',dispString) + oo_=pm3(M_,options_,oo_,meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',... + varlist_forecast_ME,texnames,obs_names,... + varlist_forecast_ME,'PointForecastME',DirectoryName,'_forc_point_ME',dispString); end end if options_.filter_covariance - pm3(endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',... - '',varlist,M_.endo_names_tex,M_.endo_names,... + oo_=pm3(M_,options_,oo_,endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',... + varlist,M_.endo_names_tex,M_.endo_names,... varlist,'FilterCovariance',DirectoryName,'_filter_covar',dispString); end diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m index babb70c3b..74ee479e9 100644 --- a/matlab/prior_posterior_statistics_core.m +++ b/matlab/prior_posterior_statistics_core.m @@ -47,14 +47,17 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ oo_ M_ bayestopt_ estim_params_ - if nargin<4 whoiam=0; end % Reshape 'myinputs' for local computation. % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: +M_=myinputs.M_; +oo_=myinputs.oo_; +options_=myinputs.options_; +estim_params_=myinputs.estim_params_; +bayestopt_=myinputs.bayestopt_; type=myinputs.type; run_smoother=myinputs.run_smoother; diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod index eaaddf7de..13f686415 100644 --- a/tests/estimation/fs2000.mod +++ b/tests/estimation/fs2000.mod @@ -132,4 +132,4 @@ oo_.MarginalDensity.LaplaceApproximation = Laplace; %reset correct Laplace %test prior sampling options_.prior_draws=100; options_.no_graph.posterior=0; -prior_posterior_statistics('prior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts +oo_=prior_posterior_statistics('prior',dataset_,dataset_info,M_,oo_,options_,estim_params_,bayestopt_); %get smoothed and filtered objects and forecasts