1function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_) 2% function [xparams, logpost]=metropolis_draw(init) 3% Builds draws from metropolis 4% 5% INPUTS: 6% init: scalar equal to 7% 1: first call to store the required information 8% on files, lines, and chains to be read 9% in persistent variables to make them available 10% for future calls 11% 0: load a parameter draw 12% Additional Inputs required for initialization 13% options_ [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}). 14% M_ [structure] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}). 15% estim_params_ [structure] Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}). 16% 17% OUTPUTS: 18% xparams: if init==0: vector of estimated parameters 19% if init==1: error flaog 20% logpost: log of posterior density 21% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}). 22% 23% SPECIAL REQUIREMENTS 24% 25% Requires CutSample to be run before in order to set up mh_history-file 26 27% Copyright (C) 2003-2017 Dynare Team 28% 29% This file is part of Dynare. 30% 31% Dynare is free software: you can redistribute it and/or modify 32% it under the terms of the GNU General Public License as published by 33% the Free Software Foundation, either version 3 of the License, or 34% (at your option) any later version. 35% 36% Dynare is distributed in the hope that it will be useful, 37% but WITHOUT ANY WARRANTY; without even the implied warranty of 38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39% GNU General Public License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 43 44persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns 45 46xparams = 0; 47logpost = 0; 48 49if init 50 %get number of parameters 51 nvx = estim_params_.nvx; 52 nvn = estim_params_.nvn; 53 ncx = estim_params_.ncx; 54 ncn = estim_params_.ncn; 55 np = estim_params_.np ; 56 npar = nvx+nvn+ncx+ncn+np; 57 %get path of metropolis files 58 MetropolisFolder = CheckPath('metropolis',M_.dname); 59 FileName = M_.fname; 60 BaseName = [MetropolisFolder filesep FileName]; 61 %load mh_history-file with info on what to load 62 load_last_mh_history_file(MetropolisFolder, FileName); 63 FirstMhFile = record.KeepedDraws.FirstMhFile; 64 FirstLine = record.KeepedDraws.FirstLine; 65 TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); 66 LastMhFile = TotalNumberOfMhFiles; 67 TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); 68 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); 69 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ? 70 mh_nblck = options_.mh_nblck; 71 % set sub_draws option if empty 72 if isempty(options_.sub_draws) 73 options_.sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws)); 74 else 75 if options_.sub_draws>NumberOfDraws*mh_nblck 76 skipline() 77 disp(['Estimation::mcmc: The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!']) 78 disp('Estimation::mcmc: You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).') 79 skipline() 80 xparams = 1; % xparams is interpreted as an error flag 81 end 82 end 83 return 84else %not initialization, return one draw 85 %get random draw from random chain 86 ChainNumber = ceil(rand*mh_nblck); 87 DrawNumber = ceil(rand*NumberOfDraws); 88 89 if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line 90 MhFilNumber = FirstMhFile; 91 MhLine = FirstLine+DrawNumber-1; 92 else %draw in other file 93 DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); 94 MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); 95 MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns; 96 end 97 %load parameters and posterior 98 load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2'); 99 xparams = x2(MhLine,:); 100 logpost= logpo2(MhLine); 101end