1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
60
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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
function Scale = calibrate_mh_scale_parameter(ObjectiveFunction, CovarianceMatrix, Parameters, MhBounds, options, varargin)
% function Scale = calibrate_mh_scale_parameter(ObjectiveFunction, CovarianceMatrix, Parameters, MhBounds, options, varargin)
% Tune the MH scale parameter so that the overall acceptance ratio is close to AcceptanceTarget.
%
% INPUTS
% - ObjectiveFunction [fhandle] Function (posterior kernel).
% - CovarianceMatrix [double] n*n matrix, covariance matrix of the jumping distribution.
% - Parameters [double] n*1 vector, parameter values.
% - MhBounds [double] n*2 matrix, bounds on the possible values for the parameters.
% - options [structure] content of options_.tune_mh_jscale.
% - varargin [cell] Additional arguments to be passed to ObjectiveFunction.
%
% OUTPUTS
% - Scale [double] scalar, optimal scale parameter for the jumping distribution.
%
% Note: program terminates if c3 consecutive runs of stepsize draws occured where
% i) the overall acceptance rate was less than c1 from target and
% ii) less than c2 over the last stepsize=2000 draws.
% Adjustment between steps takes place using a weighted average with the exponent being rho
% Copyright © 2020-2023 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 <https://www.gnu.org/licenses/>.
% Fire up the wait bar
hh_fig = dyn_waitbar(0,'Tuning of the scale parameter...');
set(hh_fig,'Name','Tuning of the scale parameter.');
% Intilialize various counters.
j = 1; jj = 1; isux = 0; jsux = 0; i = 0;
% Evaluate the objective function.
logpo0 = - feval(ObjectiveFunction, Parameters, varargin{:});
logpo1 = logpo0;
% Get the dimension of the problem.
n = length(Parameters);
% Initialize the correction on the scale factor.
correction = 1.0;
% Set the initial value of the scale parameter
if isempty(options.guess)
options.guess=2.38/sqrt(length(Parameters));
end
Scale = options.guess;
% Transposition of some arrays.
MhBounds = MhBounds';
Parameters = Parameters';
% Compute the Cholesky of the covariance matrix, return an error if the
% matrix is not positive definite.
try
dd = chol(CovarianceMatrix);
catch
error('The covariance matrix has to be a symmetric positive definite matrix!')
end
% Set parameters related to the proposal distribution
if options.rwmh.proposal_distribution=='rand_multivariate_normal'
nu = n;
elseif options.rwmh.proposal_distribution=='rand_multivariate_student'
nu = options.rwmh.student_degrees_of_freedom;
end
% Random Walk Metropolis Hastings iterations...
while j<=options.maxiter
% Obtain a proposal (jump)
proposal = feval(options.rwmh.proposal_distribution, Parameters, Scale*dd, nu);
% If out of boundaries set the posterior kernel equal to minus infinity
% so that the proposal will be rejected with probability one.
if all(proposal > MhBounds(1,:)) && all(proposal < MhBounds(2,:))
logpo0 = -feval(ObjectiveFunction, proposal(:), varargin{:});
else
logpo0 = -inf;
end
% Move if the proposal is enough likely...
if logpo0>-inf && log(rand)<logpo0-logpo1
Parameters = proposal;
logpo1 = logpo0;
isux = isux + 1;
jsux = jsux + 1;
end% ... otherwise I don't move.
prtfrc = j/options.maxiter;
% Update the waitbar
if ~mod(j, 10)
dyn_waitbar(prtfrc, hh_fig, sprintf('Acceptance ratio [during last %u]: %f [%f]', options.stepsize, isux/j, jsux/jj));
end
% Adjust the value of the scale parameter.
if ~mod(j, options.stepsize)
r1 = jsux/jj; % Local acceptance ratio
r2 = isux/j; % Overall acceptance ratio
% Set correction for the scale factor
c1 = r1/options.target;
if abs(c1-1)>.05
correction = correction^options.rho*c1^(1-options.rho);
else
correction = c1;
end
% Apply correction
if c1>0
Scale = Scale*correction;
else
Scale = Scale/10;
end
% Update some counters.
jsux = 0; jj = 0;
if abs(r2-options.target)<options.c2 && abs(r1-options.target)<options.c1
i = i+1;
else
i = 0;
end
% Test convergence.
if i>options.c3
break
end
end
j = j+1;
jj = jj + 1;
end
dyn_waitbar_close(hh_fig);
|
