stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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
time-shift
adjemian 2008-06-17 21:23:02 +00:00
parent 120d412152
commit 2a632d2efc
4 changed files with 189 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
matlab/dsge_posterior_theoretical_covariance.m
View File

@ -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,:));

40
matlab/dsge_posterior_theoretical_variance_decomposition.m
View File

@ -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
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];

218
matlab/selec_posterior_draws.m
View File

@ -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
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

7
matlab/size_of_the_reduced_form_model.m Normal file
View File

@ -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 ;