adaptive_metropolis_hastings.m 6.21 KB
Newer Older
1
function record=adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
2
3
4
5
6
7
8
9
10
11
12
13
%function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
% Random walk Metropolis-Hastings algorithm. 
% 
% INPUTS 
%   o TargetFun  [char]     string specifying the name of the objective
%                           function (posterior kernel).
%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
%   o vv         [double]   (p*p) matrix, posterior covariance matrix (at the mode).
%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters. 
%   o varargin              list of argument following mh_bounds
%  
% OUTPUTS 
14
%   o record     [struct]   structure describing the iterations
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
%
% ALGORITHM 
%   Metropolis-Hastings.       
%
% SPECIAL REQUIREMENTS
%   None.
%
% PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed
% in parallel. The code sutable to be executed in
% parallel on multi core or cluster machine (in general a 'for' cycle),
% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
% Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in
% parallel and called name_function.m and name_function_core.m. 
% In addition in parallel package we have second set of functions used
% to manage the parallel computation.
%
% This function was the first function to be parallelized, later other
% functions have been parallelized using the same methodology.
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management funtions.

% Copyright (C) 2006-2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

global M_ options_ bayestopt_ estim_params_ oo_


old_options = options_;

accept_target = options_.amh.accept_target;
60
m_directory = [M_.fname '/metropolis/']; 
61

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
if options_.load_mh_file == 0
    delete([m_directory 'adaptive_metropolis_proposal_*.mat']);
    nP = 0;
else
    D = dir([m_directory 'adaptive_metropolis_proposal_*.mat']);
    nP = size(D,1);
end;

if nP == 0
    jscale = options_.mh_jscale;
    bayestopt_.jscale = jscale;
    save([m_directory 'adaptive_metropolis_proposal_0'],'vv','jscale');
    nP = 1;
else
    tmp = load([m_directory 'adaptive_metropolis_proposal_' ...
                int2str(nP-1)],'vv','jscale');
    vv = tmp.vv;
    bayestopt_.jscale = tmp.jscale;
end

if options_.amh.cova_steps
    bayestopt_.jscale = tune_scale_parameter(TargetFun, ...
                                              ProposalFun,xparam1,vv,mh_bounds,varargin{:});
end
86
87
88
89
90

for i=1:options_.amh.cova_steps
    options_.mh_replic = options_.amh.cova_replic;
    random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
                                    xparam1,vv,mh_bounds,varargin{:});
91
92
    tot_draws = total_draws(M_);
    options_.mh_drop = (tot_draws-options_.amh.cova_replic)/tot_draws;
93
94
    CutSample(M_,options_,estim_params_);
    [junk,vv] = compute_mh_covariance_matrix();
95
96
97
98
99
    jscale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin{:});
    bayestopt_.jscale = jscale;
    save([m_directory 'adaptive_metropolis_proposal_' ...
          int2str(nP)],'vv','jscale');
    nP = nP + 1;
100
101
102
end

options_.mh_replic = old_options.mh_replic;
103
options_.mh_drop = old_options.mh_drop;
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
                                         xparam1,vv,mh_bounds,varargin{:});
                


function selected_scale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
global options_ bayestopt_

selected_scale = [];

maxit = options_.amh.scale_tuning_maxit;
accept_target = options_.amh.accept_target;
test_runs = options_.amh.scale_tuning_test_runs;
tolerance = options_.amh.scale_tuning_tolerance;
Scales = zeros(maxit,1);
AvRates = zeros(maxit,1);
120
Scales(1) = bayestopt_.jscale;
121
122
123
124
125
126
127
128
129

for i=1:maxit
    options_.mh_replic = options_.amh.scale_tuning_blocksize;
    bayestopt_.jscale = Scales(i);
    record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
                                             xparam1,vv, ...
                                             mh_bounds,varargin{:});
    AvRates(i) = mean(record.AcceptationRates);    

130
131
132
133
134
135
136
137
138
139
140
    if i < test_runs
        i_kept_runs = 1:i;
    else
        i_kept_runs = i_kept_runs+1;
    end
    
    kept_runs_s = Scales(i_kept_runs);
    kept_runs_a = AvRates(i_kept_runs);
    
    if i > test_runs
        a_mean = mean(kept_runs_a);
141
        if (a_mean > (1-tolerance)*accept_target) && ...
142
143
144
                (a_mean < (1+tolerance)*accept_target) && ...
                all(kept_runs_a > (1-test_runs*tolerance)*accept_target) && ...
                all(kept_runs_a < (1+test_runs*tolerance)*accept_target)
145
146
147
148
149
150
            selected_scale = mean(Scales((i-test_runs+1):i));
            disp(['Selected scale: ' num2str(selected_scale)])
            return
        end
    end
    
151
152
153
154
    mean_kept_runs_a = mean(kept_runs_a);
    if (mean_kept_runs_a/accept_target < 1-test_runs*tolerance) ...
            || (mean_kept_runs_a/accept_target > 1+test_runs*tolerance)
        damping_factor = 1
155
    else
156
        damping_factor = 1/3
157
    end
158
159
160
161
    Scales(i+1) = mean(kept_runs_s)*(mean(kept_runs_a)/ ...
                                     accept_target)^damping_factor;


162
163
    options_.load_mh_file = 1;
    
164
165
166
167
168
169
170
171
172
173
    disp(100*kept_runs_s')
    disp(100*kept_runs_a')
    disp(['Selected scale ' num2str(Scales(i+1))])    
end

error('AMH scale tuning: tuning didn''t converge')

function y = total_draws(M_)
load([M_.fname '/metropolis/' M_.fname '_mh_history'])
y = sum(record.MhDraws(:,1));