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

New functions. posterior draws used after the metropolis (forecasts, irfs, moments) are built and saved just after the call to metropolis.m 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@964 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2006-10-10 19:37:26 +00:00
parent cf650ea461
commit 897a084881
5 changed files with 252 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/PosteriorIRF.m
View File

@ -26,7 +26,7 @@ ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
%%
%
MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
nn = sqrt(MaxNumberOfPlotPerFigure);
DirectoryName = CheckPath('Output');

8
matlab/dynare_estimation.m
View File

@ -817,6 +817,9 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
end
end
pindx = estim_params_.param_vals(:,1);
save([M_.fname '_params'],'pindx');
if (any(bayestopt_.pshape >0 ) & options_.mh_replic) | ...
(any(bayestopt_.pshape >0 ) & options_.load_mh_file) %% not ML estimation
bounds = prior_bounds(bayestopt_);
@ -1323,5 +1326,8 @@ if ~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
end
if options_.forecast > 0 & options_.mh_replic == 0 & ~options_.load_mh_file
forecast(var_list_);
forecast(var_list);
end
pindx = estim_params_.param_vals(:,1);
save([M_.fname '_pindx','pindx']);

132
matlab/prior_draw.m Normal file
View File

@ -0,0 +1,132 @@
function pdraw = prior_draw(init,cc)
% Build one draw from the prior distribution.
%
% INPUTS
% o SampleSize [integer] Size of the sample to build
%
% OUTPUTS
% None.
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% Gnu Public License.
global M_ options_ estim_params_
persistent fname npar bounds pshape pmean pstd a b p3 p4
if init
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];
pshape = bayestopt_.pshape;
pmean = bayestopt_.pmean;
pstd = bayestopt_.pstdev;
p1 = bayestopt_.p1;
p2 = bayestopt_.p2;
p3 = bayestopt_.p3;
p4 = bayestopt_.p4;
a = zeros(npar,1);
b = zeros(npar,1);
if nargin == 2
bounds = cc;
else
bounds = [-Inf Inf];
end
for i = 1:npar
if pshape(i) == 3
b(i) = pstd(i)^2/(pmean(i)-p3(i));
a(i) = (pmean(i)-p3(i))/b(i);
elseif pshape(i) == 1
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
a(i) = (1-mu)*mu^2/stdd^2 - mu;
b(i) = a*(1/mu - 1);
elseif pshape(i) ==
end
pdraw = zeros(npar,1);
end
condition = 1;
pdraw = zeros(npar,1);
return
end
for i = 1:npar
switch pshape(i)
case 5% Uniform prior.
pdraw(i) = rand*(p4(i)-p3(i)) + p3(i);
case 3% Gaussian prior.
while condition
tmp = randn*pstd(i) + pmean(i);
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
pdraw(i) = tmp;
break
end
end
case 2% Gamma prior.
while condition
g = gamma_draw(a(i),b(i),p3(i));
if g >= bounds(i,1) && g <= bounds(i,2)
pdraw(i) = g;
break
end
end
case 1% Beta distribution (TODO: generalized beta distribution)
while condition
y1 = gamma_draw(a(i),1,0);
y2 = gamma_draw(b(i),1,0);
tmp = y1/(y1+y2);
if tmp >= bounds(i,1) && tmp <= bounds(i,2)
pdraw(i) = tmp;
break
end
end
case 4% INV-GAMMA1 distribution
case 6% INV-GAMMA2 distribution
otherwise
disp('prior_draw:: Error!')
disp('Unknown prior distribution.')
pdraw(i) = NaN;
end
end
function gamma_draw(a,b,c)
% Bauwens, Lubrano & Richard (page 316)
if a >30
z = randn;
g = b*(z+sqrt(4*a-1))^2/4 + c;
else
x = -1;
while x<0
u1 = rand;
y = tan(pi*u1);
x = y*sqrt(2*a-1)+a-1;
end
while condition
u2 = rand;
if log(u2) <= log(1+y^2)+(a-1)*log(x/(a-1))-y*sqrt(2*a-1);
break
end
end
g = x*b+c;
end

112
matlab/selec_posterior_draws.m Normal file
View File

@ -0,0 +1,112 @@
function SampleAddress = selec_posterior_draws(SampleSize,info,filepath,filename)
% Selects a sample of draws from the posterior distribution.
%
% INPUTS
% o SampleSize [integer] Size of the sample to build
% o info [integer] Ff 1 then posterior draws are saved on disk.
% OUTPUTS
% o SampleAddress [integer] A (SampleSize*4) array, each line specifies the
% location of a posterior draw:
% Column 2 --> Chain number
% Column 3 --> (mh) File number
% Column 4 --> (mh) line number
%
% ALGORITHM
% None.
%
% SPECIAL REQUIREMENTS
% None.
%
% SPECIAL THANKS
% o Antoine.
%
% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
% 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;
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')
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')
end
end

1
matlab/simk.m
View File

@ -143,7 +143,6 @@ for iter = 1:options_.maxit
ix = indnv(jwci,iz) ;
iy__ = indnv(icc1,iz) ;
temp = zeros(size(w,1),size(iz,1)) ;
temp(:,ix) = w ;
temp(:,iy__) = temp(:,iy__)-w0*c(j1i,1:ncc1) ;
w = temp ;
jwci = iz ;