Make default jscale depend on number of estimated parameters

Related to https://git.dynare.org/Dynare/dynare/-/issues/1849
unit-tests
Johannes Pfeifer 2022-05-31 21:23:57 +02:00
parent 6185767f63
commit 9de0688cfe
6 changed files with 64 additions and 20 deletions

View File

@ -6381,9 +6381,8 @@ block decomposition of the model (see :opt:`block`).
.. option:: mh_jscale = DOUBLE
The scale parameter of the jumping distributions covariance
matrix (Metropolis-Hastings or TaRB-algorithm). The default
value is rarely satisfactory. This option must be tuned to
The scale parameter of the jumping distribution's covariance
matrix (Metropolis-Hastings or TaRB-algorithm). This option must be tuned to
obtain, ideally, an acceptance ratio of 25%-33%. Basically, the
idea is to increase the variance of the jumping distribution if
the acceptance ratio is too high, and decrease the same
@ -6401,7 +6400,7 @@ block decomposition of the model (see :opt:`block`).
subsequent runs via the ``posterior_sampler_options`` option
:ref:`scale_file <scale-file>`. Both ``mode_compute=6`` and
``scale_file`` will overwrite any value specified in
``estimated_params`` with the tuned value. Default: ``0.2``.
``estimated_params`` with the tuned value. Default: ``2.38/sqrt(n)``.
Note also that for the Random Walk Metropolis Hastings
algorithm, it is possible to use option :opt:`mh_tune_jscale
@ -6409,7 +6408,7 @@ block decomposition of the model (see :opt:`block`).
of ``mh_jscale``. In this case, the ``mh_jscale`` option must
not be used.
.. option:: mh_init_scale = DOUBLE
.. option:: mh_init_scale = DOUBLE (deprecated)
The scale to be used for drawing the initial value of the
Metropolis-Hastings chain. Generally, the starting points
@ -6421,8 +6420,8 @@ block decomposition of the model (see :opt:`block`).
at the beginning of Dynare execution, i.e. the default will not
take into account potential changes in ``mh_jscale`` introduced
by either ``mode_compute=6`` or the
``posterior_sampler_options`` option :ref:`scale_file
<scale-file>`. If ``mh_init_scale`` is too wide during
``posterior_sampler_options`` option :ref:`scale_file<scale-file>`.
If ``mh_init_scale`` is too wide during
initalization of the posterior sampler so that 100 tested draws
are inadmissible (e.g. Blanchard-Kahn conditions are always
violated), Dynare will request user input of a new
@ -6434,6 +6433,26 @@ block decomposition of the model (see :opt:`block`).
most 10 times, at which point Dynare will abort with an error
message.
.. option:: mh_init_scale_factor = DOUBLE
The multiple of ``mh_jscale`` used for drawing the initial value of the
Metropolis-Hastings chain. Generally, the starting points
should be overdispersed for the *Brooks and Gelman (1998)*
convergence diagnostics to be meaningful. Default:
``2``
If ``mh_init_scale_factor`` is too wide during
initalization of the posterior sampler so that 100 tested draws
are inadmissible (e.g. Blanchard-Kahn conditions are always
violated), Dynare will request user input of a new
``mh_init_scale_factor`` value with which the next 100 draws will be
drawn and tested. If the :opt:`nointeractive` option has been
invoked, the program will instead automatically decrease
``mh_init_scale_factor`` by 10 percent after 100 futile draws and try
another 100 draws. This iterative procedure will take place at
most 10 times, at which point Dynare will abort with an error
message.
.. option:: mh_tune_jscale [= DOUBLE]
Automatically tunes the scale parameter of the jumping
@ -6450,7 +6469,7 @@ block decomposition of the model (see :opt:`block`).
.. option:: mh_tune_guess = DOUBLE
Specifies the initial value for the :opt:`mh_tune_jscale
<mh_tune_jscale [= DOUBLE]>` option. Default: ``0.2``. Must not
<mh_tune_jscale [= DOUBLE]>` option. Default: ``2.39/sqrt(n)``. Must not
be set if :opt:`mh_tune_jscale <mh_tune_jscale [= DOUBLE]>` is
not used.

View File

@ -55,6 +55,9 @@ n = length(Parameters);
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.

View File

@ -396,9 +396,9 @@ options_.logged_steady_state = 0;
options_.mh_conf_sig = 0.90;
options_.prior_interval = 0.90;
options_.mh_drop = 0.5;
options_.mh_jscale = 0.2;
options_.mh_jscale = [];
options_.mh_tune_jscale.status = false;
options_.mh_tune_jscale.guess = .2;
options_.mh_tune_jscale.guess = [];
options_.mh_tune_jscale.target = .33;
options_.mh_tune_jscale.maxiter = 200000;
options_.mh_tune_jscale.rho = .7;
@ -406,7 +406,7 @@ options_.mh_tune_jscale.stepsize = 1000;
options_.mh_tune_jscale.c1 = .02;
options_.mh_tune_jscale.c2 = .05;
options_.mh_tune_jscale.c3 = 4;
options_.mh_init_scale = 2*options_.mh_jscale;
options_.mh_init_scale_factor = 2;
options_.mh_initialize_from_previous_mcmc.status = false;
options_.mh_initialize_from_previous_mcmc.directory = '';
options_.mh_initialize_from_previous_mcmc.record = '';

View File

@ -534,8 +534,11 @@ if options_.analytic_derivation
end
end
% If jscale isn't specified for an estimated parameter, use global option options_.jscale, set to 0.2, by default.
% If jscale isn't specified for an estimated parameter, use global option options_.jscale, set to optimal value for a normal distribution by default.
% Note that check_posterior_sampler_options and mode_compute=6 may overwrite the setting
if isempty(options_.mh_jscale)
options_.mh_jscale=2.38/sqrt(length(xparam1));
end
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;

View File

@ -197,7 +197,14 @@ if ~options_.load_mh_file && ~options_.mh_recover
if isempty(d)
candidate = prior_draw();
else
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
if isfield(options_,'mh_init_scale')
if trial==1
fprintf('\nposterior_sampler_initialization: the mh_init_scale-option is deprecated. You should use the mh_init_scale_factor-option instead.\n')
end
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
else
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale_factor*options_.mh_jscale, npar);
end
end
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2(j,new_estimated_parameters) = candidate(new_estimated_parameters);
@ -219,13 +226,25 @@ if ~options_.load_mh_file && ~options_.mh_recover
if init_iter > 100 && validate == 0
disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.'])
if options_.nointeractive
disp(['Estimation::mcmc: I reduce mh_init_scale by 10 percent:'])
options_.mh_init_scale = .9*options_.mh_init_scale;
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale))
disp(['Estimation::mcmc: I reduce mh_init_scale_factor by 10 percent:'])
if isfield(options_,'mh_init_scale')
options_.mh_init_scale = .9*options_.mh_init_scale;
fprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.\n',options_.mh_init_scale)
else
options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor;
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is now equal to %f.\n',options_.mh_init_scale_factor)
end
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is now equal to %f.\n',options_.mh_init_scale_factor)
else
disp(['Estimation::mcmc: You should reduce mh_init_scale...'])
disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... ');
if isfield(options_,'mh_init_scale')
disp(['Estimation::mcmc: You should reduce mh_init_scale...'])
fprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.\n',options_.mh_init_scale)
options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... ');
else
disp(['Estimation::mcmc: You should reduce mh_init_scale_factor...'])
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is equal to %f.\n',options_.mh_init_scale_factor)
options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... ');
end
end
trial = trial+1;
end

@ -1 +1 @@
Subproject commit c48248fc0d65ed70fb13dd508ed0254f113cbcb6
Subproject commit c721718bf621a1c16803216d6a4156931b7cc373