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

Replace global variables in CutSample.m and metropolis_draw.m by function arguments 

Makes them usable outside of their current environment
time-shift
Johannes Pfeifer 2015-04-04 19:36:57 +02:00
parent 5131f3907f
commit ca5c714e29
3 changed files with 49 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
matlab/CutSample.m
View File

@ -1,12 +1,14 @@
function CutSample(M_, options_, estim_params_) function CutSample(M_, options_, estim_params_)
% function CutSample(M_, options_, estim_params_)
% function CutSample() % Takes a subset from metropolis draws by storing the required information
% Takes a subset from metropolis % 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 % INPUTS
% options_ [structure] % M_ [structure] Dynare model structure
% estim_params_ [structure] % options_ [structure] Dynare options structure
% M_ [structure] % estim_params_ [structure] Parameter structure
% %
% OUTPUTS % OUTPUTS
% none % none
@ -14,7 +16,7 @@ function CutSample(M_, options_, estim_params_)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2005-2012 Dynare Team % Copyright (C) 2005-2015 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %

4
matlab/dynare_estimation_1.m
View File

@ -450,7 +450,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if ~options_.nodiagnostic && options_.mh_replic>0 if ~options_.nodiagnostic && options_.mh_replic>0
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_); oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
end end
%% Here i discard first half of the draws: %% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_); CutSample(M_, options_, estim_params_);
%% Estimation of the marginal density from the Mh draws: %% Estimation of the marginal density from the Mh draws:
if options_.mh_replic if options_.mh_replic
@ -467,7 +467,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
else else
load([M_.fname '_results'],'oo_'); load([M_.fname '_results'],'oo_');
end end
error_flag = metropolis_draw(1); [error_flag,junk,options_]= metropolis_draw(1,options_,estim_params_,M_);
if options_.bayesian_irf if options_.bayesian_irf
if error_flag if error_flag
error('Estimation::mcmc: I cannot compute the posterior IRFs!') error('Estimation::mcmc: I cannot compute the posterior IRFs!')

61
matlab/metropolis_draw.m
View File

@ -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) % function [xparams, logpost]=metropolis_draw(init)
% Builds draws from metropolis % Builds draws from metropolis
% %
% INPUTS: % 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: % OUTPUTS:
% xparams: vector of estimated parameters % xparams: if init==0: vector of estimated parameters
% if init==1: error flaog
% logpost: log of posterior density % logpost: log of posterior density
% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% %
% SPECIAL REQUIREMENTS % 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. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % 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 persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
xparams = 0; xparams = 0;
logpost = 0; logpost = 0;
if init if init
%get number of parameters
nvx = estim_params_.nvx; nvx = estim_params_.nvx;
nvn = estim_params_.nvn; nvn = estim_params_.nvn;
ncx = estim_params_.ncx; ncx = estim_params_.ncx;
ncn = estim_params_.ncn; ncn = estim_params_.ncn;
np = estim_params_.np ; np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np; npar = nvx+nvn+ncx+ncn+np;
%get path of metropolis files
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname; FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName]; BaseName = [MetropolisFolder filesep FileName];
%load mh_history-file with info on what to load
load_last_mh_history_file(MetropolisFolder, FileName); load_last_mh_history_file(MetropolisFolder, FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
@ -52,7 +66,7 @@ if init
LastMhFile = TotalNumberOfMhFiles; LastMhFile = TotalNumberOfMhFiles;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); 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; mh_nblck = options_.mh_nblck;
% set sub_draws option if empty % set sub_draws option if empty
if isempty(options_.sub_draws) if isempty(options_.sub_draws)
@ -67,20 +81,21 @@ if init
end end
end end
return return
end else %not initialization, return one draw
%get random draw from random chain
ChainNumber = ceil(rand*mh_nblck);
DrawNumber = ceil(rand*NumberOfDraws);
ChainNumber = ceil(rand*mh_nblck); if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line
DrawNumber = ceil(rand*NumberOfDraws); MhFilNumber = FirstMhFile;
MhLine = FirstLine+DrawNumber-1;
if DrawNumber <= MAX_nruns-FirstLine+1 else %draw in other file
MhFilNumber = FirstMhFile; DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
MhLine = FirstLine+DrawNumber-1; MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
else MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); end
MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); %load parameters and posterior
MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns; load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
end xparams = x2(MhLine,:);
logpost= logpo2(MhLine);
load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2'); end
xparams = x2(MhLine,:);
logpost= logpo2(MhLine);