From 2a632d2efc4f268b22d700689b6c16175c9fe01b Mon Sep 17 00:00:00 2001 From: adjemian Date: Tue, 17 Jun 2008 21:23:02 +0000 Subject: [PATCH] Computation of 2nd order posterior moments (depending on the size of the model the posterior distribution of dr is saved on disk). git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1879 ac1d8469-bf42-47a9-8791-bf33cf982152 --- .../dsge_posterior_theoretical_covariance.m | 40 ++-- ...erior_theoretical_variance_decomposition.m | 40 +++- matlab/selec_posterior_draws.m | 218 ++++++++++-------- matlab/size_of_the_reduced_form_model.m | 7 + 4 files changed, 189 insertions(+), 116 deletions(-) create mode 100644 matlab/size_of_the_reduced_form_model.m diff --git a/matlab/dsge_posterior_theoretical_covariance.m b/matlab/dsge_posterior_theoretical_covariance.m index e6e70bdd9..dc45a3a41 100644 --- a/matlab/dsge_posterior_theoretical_covariance.m +++ b/matlab/dsge_posterior_theoretical_covariance.m @@ -24,7 +24,7 @@ function dsge_posterior_theoretical_covariance() global M_ options_ oo_ type = 'posterior';% To be defined as a input argument later... -NumberOfSimulations = 800;% To be defined in a global structure... +NumberOfSimulations = 1000;% To be defined in a global structure... % Set varlist (vartan) [ivar,vartan] = set_stationary_variables_list; @@ -40,8 +40,13 @@ fname = [ MhDirectoryName '/' M_.fname]; DrawsFiles = dir([fname '_' type '_draws*' ]); if ~rows(DrawsFiles) if strcmpi(type,'posterior') - SampleAddress = selec_posterior_draws(NumberOfSimulations,1); - else% (samples from the prior) To be done later... + drsize = size_of_the_reduced_form_model(oo_.dr); + if drsize*NumberOfSimulations>101 %Too big model! + drsize=0; + end + SampleAddress = selec_posterior_draws(NumberOfSimulations,drsize); + else + % (samples from the prior) To be done later... end DrawsFiles = dir([fname '_' type '_draws*']); end @@ -49,10 +54,8 @@ end % Get the number of stationary endogenous variables. nvar = length(ivar); - - nar = options_.ar;% Saves size of the auto-correlation function. -options_.ar = 0;% Set the size of the auto-correlation function. +options_.ar = 0;% Set the size of the auto-correlation function to zero. NumberOfDrawsFiles = rows(DrawsFiles); MaXNumberOfCovarLines = ceil(options_.MaxNumberOfBytes/(nvar*(nvar+1)/2)/8); @@ -74,11 +77,15 @@ linea = 0; for file = 1:NumberOfDrawsFiles load([MhDirectoryName '/' DrawsFiles(file).name]); NumberOfDraws = rows(pdraws); + isdrsaved = cols(pdraws)-1; for linee = 1:NumberOfDraws linea = linea+1; - draw = pdraws(linee,:); - set_parameters(draw); - [dr,info] = resol(oo_.steady_state,0); + if isdrsaved + dr = pdraws{linee,2}; + else + set_parameters(pdraws{linee,1}); + [dr,info] = resol(oo_.steady_state,0); + end tmp = th_autocovariances(dr,ivar); for i=1:nvar for j=i:nvar @@ -90,13 +97,14 @@ for file = 1:NumberOfDrawsFiles CovarFileNumber = CovarFileNumber + 1; linea = 0; test = CovarFileNumber-NumberOfCovarFiles; - if ~(CovarFileNumber-NumberOfCovarFiles)% Prepare the last round... + if ~test% Prepare the last round... Covariance_matrix = zeros(NumberOfLinesInTheLastCovarFile,nvar*(nvar+1)/2); NumberOfCovarLines = NumberOfLinesInTheLastCovarFile; - elseif CovarFileNumber-NumberOfCovarFiles<0; + CovarFileNumber = CovarFileNumber - 1; + elseif test<0; Covariance_matrix = zeros(MaXNumberOfCovarLines,nvar*(nvar+1)/2); else - clear('Covariance_matrix'); + clear('Covariance_matrix'); end end end @@ -108,7 +116,7 @@ for i=1:nvar for j=i:nvar i1 = 1; tmp = zeros(NumberOfSimulations,1); - for file = 1:NumberOfDrawsFiles + for file = 1:CovarFileNumber load([fname '_Posterior2ndOrderMoments' int2str(file)]); i2 = i1 + rows(Covariance_matrix) - 1; tmp(i1:i2) = Covariance_matrix(:,idx(i,j,nvar)); @@ -124,7 +132,7 @@ for i=1:nvar eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = hpd_interval(2);']); eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = post_deciles;']); eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = density;']); - end + end end @@ -133,10 +141,12 @@ end function k = idx(i,j,n) k = (i-1)*n+j-i*(i-1)/2; - function r = rows(M) r = size(M,1); +function c = cols(M) + c = size(M,2); + function name = fieldname(i,j,vlist) n1 = deblank(vlist(i,:)); n2 = deblank(vlist(j,:)); diff --git a/matlab/dsge_posterior_theoretical_variance_decomposition.m b/matlab/dsge_posterior_theoretical_variance_decomposition.m index 16d8babed..e8914c061 100644 --- a/matlab/dsge_posterior_theoretical_variance_decomposition.m +++ b/matlab/dsge_posterior_theoretical_variance_decomposition.m @@ -41,7 +41,11 @@ fname = [ MhDirectoryName '/' M_.fname]; DrawsFiles = dir([fname '_' type '_draws*' ]); if ~rows(DrawsFiles) if strcmpi(type,'posterior') - SampleAddress = selec_posterior_draws(NumberOfSimulations,1); + drsize = size_of_the_reduced_form_model(oo_.dr); + if drsize*NumberOfSimulations>101%Big model! + drsize=0; + end + SampleAddress = selec_posterior_draws(NumberOfSimulations,drsize); else% (samples from the prior) To be done later... end DrawsFiles = dir([fname '_' type '_draws*']); @@ -75,12 +79,16 @@ DecompFileNumber = 1; linea = 0; for file = 1:NumberOfDrawsFiles load([MhDirectoryName '/' DrawsFiles(file).name]); + isdrsaved = cols(pdraws)-1; NumberOfDraws = rows(pdraws); for linee = 1:NumberOfDraws linea = linea+1; - draw = pdraws(linee,:); - set_parameters(draw); - [dr,info] = resol(oo_.steady_state,0); + if isdrsaved + dr = pdraws{linee,2}; + else + set_parameters(pdraws{linee,1}); + [dr,info] = resol(oo_.steady_state,0); + end tmp = th_autocovariances(dr,ivar); for i=1:nvar Decomposition_array(linea,i) = tmp{1}(i,i); @@ -95,10 +103,11 @@ for file = 1:NumberOfDrawsFiles DecompFileNumber = DecompFileNumber + 1; linea = 0; test = DecompFileNumber-NumberOfDecompFiles; - if ~(DecompFileNumber-NumberOfDecompFiles)% Prepare the last round... + if ~test% Prepare the last round... Decomposition_array = zeros(NumberOfLinesInTheLastDecompFile,nvar*(nexo+1)); NumberOfDecompLines = NumberOfLinesInTheLastDecompFile; - elseif DecompFileNumber-NumberOfDecompFiles<0; + DecompFileNumber = DecompFileNumber - 1; + elseif test<0; Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*(nexo+1)); else clear('Decomposition_array'); @@ -114,7 +123,7 @@ for i=1:nvar for j=1:nexo i1 = 1; tmp = zeros(NumberOfSimulations,1); - for file = 1:NumberOfDrawsFiles + for file = 1:DecompFileNumber load([fname '_PosteriorVarianceDecomposition' int2str(file)]); i2 = i1 + rows(Decomposition_array) - 1; tmp(i1:i2) = Decomposition_array(:,nvar+(i-1)*nexo+j); @@ -149,4 +158,19 @@ for i=1:nvar eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.deciles.' name ' = post_deciles;']); eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']); end -end \ No newline at end of file +end + + +function k = idx(i,j,n) + k = (i-1)*n+j-i*(i-1)/2; + +function r = rows(M) + r = size(M,1); + +function c = cols(M) + c = size(M,2); + +function name = fieldname(i,j,vlist) + n1 = deblank(vlist(i,:)); + n2 = deblank(vlist(j,:)); + name = [n1 '.' n2]; \ No newline at end of file diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m index de735881e..9c34346dc 100644 --- a/matlab/selec_posterior_draws.m +++ b/matlab/selec_posterior_draws.m @@ -1,13 +1,13 @@ -function SampleAddress = selec_posterior_draws(SampleSize,info) - -% function SampleAddress = selec_posterior_draws(SampleSize,info) -% Selects a sample of draws from the posterior distribution and if info=1 -% saves the draws in _pdraws mat files (metropolis folder). This routine is more -% efficient than metropolis_draw.m because here an _mh file cannot be opened twice. +function SampleAddress = selec_posterior_draws(SampleSize,drsize) +% Selects a sample of draws from the posterior distribution and if nargin>1 +% saves the draws in _pdraws mat files (metropolis folder). If drsize>0 +% the dr structure, associated to the parameters, is also saved in _pdraws. +% This routine is more efficient than metropolis_draw.m because here an +% _mh file cannot be opened twice. % % INPUTS -% o SampleSize [integer] Size of the sample to build -% o info [integer] If 1 then posterior draws are saved on disk. +% o SampleSize [integer] Size of the sample to build. +% o drsize [double] structure dr is drsize megaoctets. % % OUTPUTS % o SampleAddress [integer] A (SampleSize*4) array, each line specifies the @@ -21,92 +21,124 @@ function SampleAddress = selec_posterior_draws(SampleSize,info) % % part of DYNARE, copyright Dynare Team (2006-2008) % Gnu Public License. - - global M_ options_ estim_params_ - 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; - - MhDirectoryName = CheckPath('metropolis'); - fname = [ MhDirectoryName '/' M_.fname]; - load([ fname '_mh_history']); - FirstMhFile = record.KeepedDraws.FirstMhFile; - FirstLine = record.KeepedDraws.FirstLine; - TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); - LastMhFile = TotalNumberOfMhFiles; - TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); - NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); - MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); - mh_nblck = options_.mh_nblck; - - SampleAddress = zeros(SampleSize,4); - for i = 1:SampleSize - ChainNumber = ceil(rand*mh_nblck); - DrawNumber = ceil(rand*NumberOfDraws); - SampleAddress(i,1) = DrawNumber; - SampleAddress(i,2) = ChainNumber; - if DrawNumber <= MAX_nruns-FirstLine+1 - MhFileNumber = FirstMhFile; - MhLineNumber = FirstLine+DrawNumber-1; - else - DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); - MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); - MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns; + global M_ options_ estim_params_ oo_ + + % Number of parameters: + npar = estim_params_.nvx; + npar = npar + estim_params_.nvn; + npar = npar + estim_params_.ncx; + npar = npar + estim_params_.ncn; + npar = npar + estim_params_.np; + + % Select one task: + switch nargin + case 1 + info = 0; + case 2 + if drsize>0 + info=2; + MAX_mega_bytes = 10;% Should be an option... + else + info=1; + end + drawsize = drsize+npar*8/1048576; + otherwise + error(['selec_posterior_draws:: Unexpected number of input arguments!']) end - SampleAddress(i,3) = MhFileNumber; - SampleAddress(i,4) = MhLineNumber; - end - SampleAddress = sortrows(SampleAddress,[3 2]); - if nargin == 2 & info - if SampleSize <= MAX_nruns% The posterior draws are saved in one file. - pdraws = zeros(SampleSize,npar); - old_mhfile = 0; - old_mhblck = 0; - for i = 1:SampleSize - mhfile = SampleAddress(i,3); - mhblck = SampleAddress(i,2); - if (mhfile ~= old_mhfile) | (mhblck ~= old_mhblck) - load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') + + % Get informations about the mcmc: + MhDirectoryName = CheckPath('metropolis'); + fname = [ MhDirectoryName '/' M_.fname]; + load([ fname '_mh_history']); + FirstMhFile = record.KeepedDraws.FirstMhFile; + FirstLine = record.KeepedDraws.FirstLine; + TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); + LastMhFile = TotalNumberOfMhFiles; + TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); + NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); + MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); + mh_nblck = options_.mh_nblck; + + % Randomly select draws in the posterior distribution: + SampleAddress = zeros(SampleSize,4); + for i = 1:SampleSize + ChainNumber = ceil(rand*mh_nblck); + DrawNumber = ceil(rand*NumberOfDraws); + SampleAddress(i,1) = DrawNumber; + SampleAddress(i,2) = ChainNumber; + if DrawNumber <= MAX_nruns-FirstLine+1 + MhFileNumber = FirstMhFile; + MhLineNumber = FirstLine+DrawNumber-1; + else + DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); + MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); + MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns; end - pdraws(i,:) = x2(SampleAddress(i,4),:); - old_mhfile = mhfile; - old_mhblck = mhblck; - end - clear('x2') - save([fname '_posterior_draws'],'pdraws') - else% The posterior draws are saved in xx files. - NumberOfFiles = ceil(SampleSize/MAX_nruns); - NumberOfLines = SampleSize - (NumberOfFiles-1)*MAX_nruns; - linee = 0; - fnum = 1; - pdraws = zeros(MAX_nruns,npar); - old_mhfile = 0; - old_mhblck = 0; - for i=1:SampleSize - linee = linee+1; - mhfile = SampleAddress(i,3); - mhblck = SampleAddress(i,2); - if (mhfile ~= old_mhfile) | (mhblck ~= old_mhblck) - load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') - end - pdraws(linee,:) = x2(SampleAddress(i,4),:); - old_mhfile = mhfile; - old_mhblck = mhblck; - if fnum < NumberOfFiles && linee == MAX_nruns - linee = 0; - save([fname '_posterior_draws' num2str(fnum)],'pdraws') - fnum = fnum+1; - if fnum < NumberOfFiles - pdraws = zeros(MAX_nruns,npar); - else - pdraws = zeros(NumberOfLines,npar); - end - end - end - save([fname '_posterior_draws' num2str(fnum)],'pdraws') + SampleAddress(i,3) = MhFileNumber; + SampleAddress(i,4) = MhLineNumber; end - end \ No newline at end of file + SampleAddress = sortrows(SampleAddress,[3 2]); + + % Selected draws in the posterior distribution, and if drsize>0 + % reduced form solutions, are saved on disk. + if info + if SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file. + pdraws = cell(SampleSize,info); + old_mhfile = 0; + old_mhblck = 0; + for i = 1:SampleSize + mhfile = SampleAddress(i,3); + mhblck = SampleAddress(i,2); + if (mhfile ~= old_mhfile) | (mhblck ~= old_mhblck) + load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') + end + pdraws(i,1) = {x2(SampleAddress(i,4),:)}; + if info-1 + set_parameters(pdraws{i,1}); + [dr,info] = resol(oo_.steady_state,0); + pdraws(i,2) = { dr }; + end + old_mhfile = mhfile; + old_mhblck = mhblck; + end + clear('x2') + save([fname '_posterior_draws'],'pdraws') + else% The posterior draws are saved in xx files. + NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize); + NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes); + NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile; + linee = 0; + fnum = 1; + pdraws = cell(NumberOfDrawsPerFile,info); + old_mhfile = 0; + old_mhblck = 0; + for i=1:SampleSize + linee = linee+1; + mhfile = SampleAddress(i,3); + mhblck = SampleAddress(i,2); + if (mhfile ~= old_mhfile) | (mhblck ~= old_mhblck) + load([fname '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2') + end + pdraws(linee,1) = {x2(SampleAddress(i,4),:)}; + if info-1 + set_parameters(pdraws{linee,1}); + [dr,info] = resol(oo_.steady_state,0); + pdraws(linee,2) = { dr }; + end + old_mhfile = mhfile; + old_mhblck = mhblck; + if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile + linee = 0; + save([fname '_posterior_draws' num2str(fnum)],'pdraws') + fnum = fnum+1; + if fnum < NumberOfFiles + pdraws = cell(NumberOfDrawsPerFile,info); + else + pdraws = cell(NumberOfLines,info); + end + end + end + save([fname '_posterior_draws' num2str(fnum)],'pdraws') + end + end \ No newline at end of file diff --git a/matlab/size_of_the_reduced_form_model.m b/matlab/size_of_the_reduced_form_model.m new file mode 100644 index 000000000..b876cb67d --- /dev/null +++ b/matlab/size_of_the_reduced_form_model.m @@ -0,0 +1,7 @@ +function mega = size_of_the_reduced_form_model(dr) + names = fieldnames(dr); + number_of_scalars = 0; + for field=1:size(names,1) + number_of_scalars = number_of_scalars + prod(size(getfield(dr,names{field}))); + end + mega = 8 * number_of_scalars / 1048576 ; \ No newline at end of file