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