diff --git a/matlab/CutSample.m b/matlab/CutSample.m
index a3efbb01ad43ead1a02b7d031206fe048b3b5e97..5ec0881a8677eda180b302911d4f6a6a35c701aa 100644
--- a/matlab/CutSample.m
+++ b/matlab/CutSample.m
@@ -1,12 +1,14 @@
function CutSample(M_, options_, estim_params_)
-
-% function CutSample()
-% Takes a subset from metropolis
+% function CutSample(M_, options_, estim_params_)
+% Takes a subset from metropolis draws by storing the required information
+% like the first MH-file to be loaded and the first line in that file to be
+% loaded into the record structure saved on harddisk into the
+% _mh_history-file
%
% INPUTS
-% options_ [structure]
-% estim_params_ [structure]
-% M_ [structure]
+% M_ [structure] Dynare model structure
+% options_ [structure] Dynare options structure
+% estim_params_ [structure] Parameter structure
%
% OUTPUTS
% none
@@ -14,7 +16,7 @@ function CutSample(M_, options_, estim_params_)
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2015 Dynare Team
%
% This file is part of Dynare.
%
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index a813fbc0abb6fefd1ea752b955fd7f5437bd9067..58282eba5ff81098613902b2cd74851af5282f5b 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -450,7 +450,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if ~options_.nodiagnostic && options_.mh_replic>0
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
end
- %% Here i discard first half of the draws:
+ %% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_);
%% Estimation of the marginal density from the Mh draws:
if options_.mh_replic
@@ -467,7 +467,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
else
load([M_.fname '_results'],'oo_');
end
- error_flag = metropolis_draw(1);
+ [error_flag,junk,options_]= metropolis_draw(1,options_,estim_params_,M_);
if options_.bayesian_irf
if error_flag
error('Estimation::mcmc: I cannot compute the posterior IRFs!')
diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m
index 771ba0a43310ddb2542ea04ccd6f005e44425efc..f0fcf49bd969524ef1cd883983d125c4ac67980a 100644
--- a/matlab/metropolis_draw.m
+++ b/matlab/metropolis_draw.m
@@ -1,18 +1,30 @@
-function [xparams, logpost]=metropolis_draw(init)
+function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
% function [xparams, logpost]=metropolis_draw(init)
% Builds draws from metropolis
%
% INPUTS:
-% init: scalar equal to 1 (first call) or 0
+% init: scalar equal to
+% 1: first call to store the required information
+% on files, lines, and chains to be read
+% in persistent variables to make them available
+% for future calls
+% 0: load a parameter draw
+% Additional Inputs required for initialization
+% options_ [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+% M_ [structure] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+% estim_params_ [structure] Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%
% OUTPUTS:
-% xparams: vector of estimated parameters
+% xparams: if init==0: vector of estimated parameters
+% if init==1: error flaog
% logpost: log of posterior density
+% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%
% SPECIAL REQUIREMENTS
-% none
+%
+% Requires CutSample to be run before in order to set up mh_history-file
-% Copyright (C) 2003-2011 Dynare Team
+% Copyright (C) 2003-2015 Dynare Team
%
% This file is part of Dynare.
%
@@ -29,22 +41,24 @@ function [xparams, logpost]=metropolis_draw(init)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-global options_ estim_params_ M_
persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
xparams = 0;
logpost = 0;
if init
+ %get number of parameters
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;
+ %get path of metropolis files
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName];
+ %load mh_history-file with info on what to load
load_last_mh_history_file(MetropolisFolder, FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
@@ -52,7 +66,7 @@ if init
LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
- MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ?
mh_nblck = options_.mh_nblck;
% set sub_draws option if empty
if isempty(options_.sub_draws)
@@ -67,20 +81,21 @@ if init
end
end
return
-end
-
-ChainNumber = ceil(rand*mh_nblck);
-DrawNumber = ceil(rand*NumberOfDraws);
+else %not initialization, return one draw
+ %get random draw from random chain
+ ChainNumber = ceil(rand*mh_nblck);
+ DrawNumber = ceil(rand*NumberOfDraws);
-if DrawNumber <= MAX_nruns-FirstLine+1
- MhFilNumber = FirstMhFile;
- MhLine = FirstLine+DrawNumber-1;
-else
- DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
- MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
- MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
-end
-
-load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
-xparams = x2(MhLine,:);
-logpost= logpo2(MhLine);
\ No newline at end of file
+ if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line
+ MhFilNumber = FirstMhFile;
+ MhLine = FirstLine+DrawNumber-1;
+ else %draw in other file
+ DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
+ MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
+ MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
+ end
+ %load parameters and posterior
+ load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
+ xparams = x2(MhLine,:);
+ logpost= logpo2(MhLine);
+end
\ No newline at end of file