diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index d9617ce09..dc6e85ae7 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -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'); diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index 6cb33a5a8..0665c32f7 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -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']); \ No newline at end of file diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m new file mode 100644 index 000000000..e8e201557 --- /dev/null +++ b/matlab/prior_draw.m @@ -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 \ No newline at end of file diff --git a/matlab/selec_posterior_draws.m b/matlab/selec_posterior_draws.m new file mode 100644 index 000000000..8028eaa5d --- /dev/null +++ b/matlab/selec_posterior_draws.m @@ -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 \ No newline at end of file diff --git a/matlab/simk.m b/matlab/simk.m index 71d55e888..40ce95a74 100644 --- a/matlab/simk.m +++ b/matlab/simk.m @@ -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 ;