stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

provisions for reworked posterior sampling options: 

- handle sub lists of individual samplers
- split checks in dynare_estimation_init.m and before running posterior_sampler.m [invhess checks]
- posterior sampler options checks moved from initial_estimation_checks.m to check_posterior_sampler_options.m
- added use_mh_covariance_matrix to imh and rwmh
- slice re-sets mode_compute=0 cova_compute=0
- updated test function
time-shift
Marco Ratto 2016-05-13 21:35:59 +02:00 committed by Johannes Pfeifer
parent 8b85ca19bf
commit 7b3c42c6e1
8 changed files with 520 additions and 331 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

598
matlab/check_posterior_sampler_options.m
View File

@ -1,193 +1,405 @@
function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_) function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds)
% function posterior_sampler_options = check_posterior_sampler_options(posterior_sampler_options, options_) % function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds)
% initialization of posterior samplers % initialization of posterior samplers
% %
% INPUTS % INPUTS
% posterior_sampler_options: posterior sampler options % posterior_sampler_options: posterior sampler options
% options_: structure storing the options % options_: structure storing the options
% OUTPUTS % OUTPUTS
% posterior_sampler_options: checked posterior sampler options % posterior_sampler_options: checked posterior sampler options
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2015 Dynare Team % Copyright (C) 2015 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
% Dynare is free software: you can redistribute it and/or modify % Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by % it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or % the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version. % (at your option) any later version.
% %
% Dynare is distributed in the hope that it will be useful, % Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of % but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details. % GNU General Public License for more details.
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampling_method;
posterior_sampler_options.proposal_distribution = options_.proposal_distribution; init=0;
bounds = posterior_sampler_options.bounds; if isempty(posterior_sampler_options),
invhess = posterior_sampler_options.invhess; init=1;
end
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice') if init,
if ~options_.cova_compute % set default options and user defined options
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.') posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method;
end posterior_sampler_options.bounds = bounds;
if options_.load_mh_file && options_.use_mh_covariance_matrix,
[junk, invhess] = compute_mh_covariance_matrix; switch posterior_sampler_options.posterior_sampling_method
posterior_sampler_options.invhess = invhess;
end case 'random_walk_metropolis_hastings'
posterior_sampler_options.parallel_bar_refresh_rate=50; posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3; posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='MH'; posterior_sampler_options.parallel_bar_title='RWMH';
posterior_sampler_options.serial_bar_title='Metropolis-Hastings'; posterior_sampler_options.serial_bar_title='RW Metropolis-Hastings';
posterior_sampler_options.save_tmp_file=1; posterior_sampler_options.save_tmp_file=1;
end
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.rwmh);
% check specific options for slice sampler
if strcmp(options_.posterior_sampling_method,'slice') % user defined options
% by default, slice sampler should trigger if ~isempty(options_.posterior_sampler_options.sampling_opt)
% mode_compute=0 and options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
% mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH) for i=1:rows(options_list)
% moreover slice must be associated to: switch options_list{i,1}
% options_.mh_posterior_mode_estimation = 0;
% this is done below, but perhaps preprocessing should do this? case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
posterior_sampler_options.parallel_bar_refresh_rate=1; strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
posterior_sampler_options.serial_bar_refresh_rate=1; error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
posterior_sampler_options.parallel_bar_title='SLICE'; 'rand_multivariate_student or rand_multivariate_normal as options']);
posterior_sampler_options.serial_bar_title='SLICE'; else
posterior_sampler_options.save_tmp_file=1; posterior_sampler_options.proposal_distribution=options_list{i,2};
posterior_sampler_options = set_default_option(posterior_sampler_options,'rotated',0); end
posterior_sampler_options = set_default_option(posterior_sampler_options,'slice_initialize_with_mode',0);
posterior_sampler_options = set_default_option(posterior_sampler_options,'use_mh_covariance_matrix',0);
posterior_sampler_options = set_default_option(posterior_sampler_options,'WR',[]); case 'student_degrees_of_freedom'
if ~isfield(posterior_sampler_options,'mode'), if options_list{i,2} <= 0
posterior_sampler_options.mode = []; error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else % multimodal case else
posterior_sampler_options.rotated = 1; posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end end
posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options = set_default_option(posterior_sampler_options,'W1',0.8*(bounds.ub-bounds.lb)); case 'use_mh_covariance_matrix'
% indicates to use the covariance matrix from previous iterations to
if ~isempty(options_.optim_opt) % define the covariance of the proposal distribution
options_list = read_key_value_string(options_.optim_opt); % default = 0
for i=1:rows(options_list) posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
switch options_list{i,1} options_.use_mh_covariance_matrix = options_list{i,2};
case 'rotated'
% triggers rotated slice iterations using a covariance otherwise
% matrix from initial burn-in iterations warning(['imh_sampler: Unknown option (' options_list{i,1} ')!'])
% must be associated with: end
% <use_mh_covariance_matrix> or <slice_initialize_with_mode> end
% default = 0 end
posterior_sampler_options.rotated = options_list{i,2};
case 'tailored_random_block_metropolis_hastings'
case 'mode' posterior_sampler_options.parallel_bar_refresh_rate=5;
% for multimodal posteriors, provide the list of modes as a posterior_sampler_options.serial_bar_refresh_rate=1;
% matrix, ordered by column, i.e. [x1 x2 x3] for three posterior_sampler_options.parallel_bar_title='TaRB-MH';
% modes x1 x2 x3 posterior_sampler_options.serial_bar_title='TaRB Metropolis-Hastings';
% MR note: not sure this is possible with the posterior_sampler_options.save_tmp_file=1;
% read_key_value_string ???
% if this is not possible <mode_files> does to job in any case % default options
% This will automatically trigger <rotated> posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.tarb);
% default = []
tmp_mode = options_list{i,2}; % user defined options
for j=1:size(tmp_mode,2), if ~isempty(options_.posterior_sampler_options.sampling_opt)
posterior_sampler_options.mode(j).m = tmp_mode(:,j); options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
end for i=1:rows(options_list)
case 'mode_files' switch options_list{i,1}
% for multimodal posteriors provide a list of mode files,
% one per mode. With this info, the code will automatically case 'proposal_distribution'
% set the <mode> option. The mode files need only to if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
% contain the xparam1 variable. strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
% This will automatically trigger <rotated> error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
% default = [] 'rand_multivariate_student or rand_multivariate_normal as options']);
posterior_sampler_options.mode_files = options_list{i,2}; else
posterior_sampler_options.proposal_distribution=options_list{i,2};
case 'slice_initialize_with_mode' end
% the default for slice is to set mode_compute = 0 in the
% preprocessor and start the chain(s) from a random location in the prior.
% This option first runs the optimizer and then starts the case 'student_degrees_of_freedom'
% chain from the mode. Associated with optios <rotated>, it will if options_list{i,2} <= 0
% use invhess from the mode to perform rotated slice error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
% iterations. else
% default = 0 posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
posterior_sampler_options.slice_initialize_with_mode = options_list{i,2}; end
case 'initial_step_size' case 'mode_compute'
% sets the initial size of the interval in the STEPPING-OUT PROCEDURE posterior_sampler_options.mode_compute=options_list{i,2};
% the initial_step_size must be a real number in [0, 1],
% and it sets the size as a proportion of the prior bounds, case 'new_block_probability'
% i.e. the size will be initial_step_size*(UB-LB) if options_list{i,2}<0 || options_list{i,2}>1
% slice sampler requires prior_truncation > 0! error('check_posterior_sampler_options:: The tarb new_block_probability must be between 0 and 1!')
% default = 0.8 else
posterior_sampler_options.W1 = options_list{i,2}*(bounds.ub-bounds.lb); posterior_sampler_options.new_block_probability=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% in association with <rotated> indicates to use the otherwise
% covariance matrix from previous iterations to define the
% rotated slice end
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2}; end
otherwise end
warning(['slice_sampler: Unknown option (' options_list{i,1} ')!'])
end case 'independent_metropolis_hastings'
end posterior_sampler_options.parallel_bar_refresh_rate=50;
end posterior_sampler_options.serial_bar_refresh_rate=3;
if options_.load_mh_file, posterior_sampler_options.parallel_bar_title='IMH';
posterior_sampler_options.slice_initialize_with_mode = 0; posterior_sampler_options.serial_bar_title='Ind. Metropolis-Hastings';
else posterior_sampler_options.save_tmp_file=1;
if ~posterior_sampler_options.slice_initialize_with_mode,
posterior_sampler_options.invhess=[]; % default options
end posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.imh);
end
% user defined options
if posterior_sampler_options.rotated, if ~isempty(options_.posterior_sampler_options.sampling_opt)
if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode), % rotated unimodal options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix) switch options_list{i,1}
skipline()
disp('I cannot start rotated slice sampler because') case 'proposal_distribution'
disp('there is no previous MCMC to load ') if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
disp('or the Hessian at the mode is not computed.') strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error('Rotated slice cannot start') error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
end 'rand_multivariate_student or rand_multivariate_normal as options']);
if isempty(invhess) else
error('oops! This error should not occur, please contact developers.') posterior_sampler_options.proposal_distribution=options_list{i,2};
end end
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix,
[junk, invhess] = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess; case 'student_degrees_of_freedom'
end if options_list{i,2} <= 0
[V1 D]=eig(invhess); error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
posterior_sampler_options.V1=V1; else
posterior_sampler_options.WR=sqrt(diag(D))*3; posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end end
end
case 'use_mh_covariance_matrix'
if ~isempty(posterior_sampler_options.mode_files), % multimodal case % indicates to use the covariance matrix from previous iterations to
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains % define the covariance of the proposal distribution
for j=1:length( modes ), % default = 0
load(modes{j}, 'xparam1') posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
mode(j).m=xparam1; options_.use_mh_covariance_matrix = options_list{i,2};
end
posterior_sampler_options.mode = mode; otherwise
posterior_sampler_options.rotated = 1; warning(['imh_sampler: Unknown option (' options_list{i,1} ')!'])
posterior_sampler_options.WR=[]; end
end end
options_.mh_posterior_mode_estimation = 0; end
end
case 'slice'
posterior_sampler_options.parallel_bar_refresh_rate=1;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='SLICE';
posterior_sampler_options.serial_bar_title='SLICE';
posterior_sampler_options.save_tmp_file=1;
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options,options_.posterior_sampler_options.slice);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'rotated'
% triggers rotated slice iterations using a covariance
% matrix from initial burn-in iterations
% must be associated with:
% <use_mh_covariance_matrix> or <slice_initialize_with_mode>
% default = 0
posterior_sampler_options.rotated = options_list{i,2};
case 'mode'
% for multimodal posteriors, provide the list of modes as a
% matrix, ordered by column, i.e. [x1 x2 x3] for three
% modes x1 x2 x3
% MR note: not sure this is possible with the
% read_key_value_string ???
% if this is not possible <mode_files> does to job in any case
% This will automatically trigger <rotated>
% default = []
tmp_mode = options_list{i,2};
for j=1:size(tmp_mode,2),
posterior_sampler_options.mode(j).m = tmp_mode(:,j);
end
case 'mode_files'
% for multimodal posteriors provide a list of mode files,
% one per mode. With this info, the code will automatically
% set the <mode> option. The mode files need only to
% contain the xparam1 variable.
% This will automatically trigger <rotated>
% default = []
posterior_sampler_options.mode_files = options_list{i,2};
case 'slice_initialize_with_mode'
% the default for slice is to set mode_compute = 0 in the
% preprocessor and start the chain(s) from a random location in the prior.
% This option first runs the optimizer and then starts the
% chain from the mode. Associated with optios <rotated>, it will
% use invhess from the mode to perform rotated slice
% iterations.
% default = 0
posterior_sampler_options.slice_initialize_with_mode = options_list{i,2};
case 'initial_step_size'
% sets the initial size of the interval in the STEPPING-OUT PROCEDURE
% the initial_step_size must be a real number in [0, 1],
% and it sets the size as a proportion of the prior bounds,
% i.e. the size will be initial_step_size*(UB-LB)
% slice sampler requires prior_truncation > 0!
% default = 0.8
if options_list{i,2}<=0 || options_list{i,2}>=1
error('check_posterior_sampler_options:: slice initial_step_size must be between 0 and 1')
else
posterior_sampler_options.initial_step_size=options_list{i,2};
end
case 'use_mh_covariance_matrix'
% in association with <rotated> indicates to use the
% covariance matrix from previous iterations to define the
% rotated slice
% default = 0
posterior_sampler_options.use_mh_covariance_matrix = options_list{i,2};
options_.use_mh_covariance_matrix = options_list{i,2};
otherwise
warning(['slice_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
% slice posterior sampler does not require mode or hessian to run
% needs to be set to 1 to skip parts in dynare_estimation_1.m
% requiring posterior maximization/calibrated smoother before MCMC
options_.mh_posterior_mode_estimation=1;
if ~ posterior_sampler_options.slice_initialize_with_mode
% by default, slice sampler should trigger
% mode_compute=0 and
% mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH)
options_.mode_compute = 0;
options_.cova_compute = 0;
else
if (isequal(options_.mode_compute,0) && isempty(options_.mode_file) )
skipline()
disp('check_posterior_sampler_options:: You have specified the option "slice_initialize_with_mode"')
disp('check_posterior_sampler_options:: to initialize the slice sampler using mode information')
disp('check_posterior_sampler_options:: but no mode file nor posterior maximization is selected,')
error('check_posterior_sampler_options:: The option "slice_initialize_with_mode" is inconsistent with mode_compute=0 or empty mode_file.')
else
options_.mh_posterior_mode_estimation=0;
end
end
if any(isinf(bounds.lb)) || any(isinf(bounds.ub)),
skipline()
disp('some priors are unbounded and prior_trunc is set to zero')
error('The option "slice" is inconsistent with prior_trunc=0.')
end
% moreover slice must be associated to:
% options_.mh_posterior_mode_estimation = 0;
% this is done below, but perhaps preprocessing should do this?
if ~isfield(posterior_sampler_options,'mode'),
posterior_sampler_options.mode = [];
else % multimodal case
posterior_sampler_options.rotated = 1;
end
% posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
if options_.load_mh_file,
posterior_sampler_options.slice_initialize_with_mode = 0;
else
if ~posterior_sampler_options.slice_initialize_with_mode,
posterior_sampler_options.invhess=[];
end
end
if ~isempty(posterior_sampler_options.mode_files), % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
for j=1:length( modes ),
load(modes{j}, 'xparam1')
mode(j).m=xparam1;
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
otherwise
end
return
end
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
skipline()
disp('check_posterior_sampler_options:: I cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed')
disp('check_posterior_sampler_options:: or there is no previous MCMC to load ')
error('check_posterior_sampler_options:: MCMC cannot start')
end
end
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix,
[junk, invhess] = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
% check specific options for slice sampler
if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
invhess = posterior_sampler_options.invhess;
if posterior_sampler_options.rotated,
if isempty(posterior_sampler_options.mode_files) && isempty(posterior_sampler_options.mode), % rotated unimodal
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
skipline()
disp('check_posterior_sampler_options:: I cannot start rotated slice sampler because')
disp('check_posterior_sampler_options:: there is no previous MCMC to load ')
disp('check_posterior_sampler_options:: or the Hessian at the mode is not computed.')
error('check_posterior_sampler_options:: Rotated slice cannot start')
end
if isempty(invhess)
error('check_posterior_sampler_options:: This error should not occur, please contact developers.')
end
% % % if options_.load_mh_file && options_.use_mh_covariance_matrix,
% % % [junk, invhess] = compute_mh_covariance_matrix;
% % % posterior_sampler_options.invhess = invhess;
% % % end
[V1 D]=eig(invhess);
posterior_sampler_options.V1=V1;
posterior_sampler_options.WR=sqrt(diag(D))*3;
end
end
% needs to be re-set to zero otherwise posterior analysis is filtered
% out in dynare_estimation_1.m
options_.mh_posterior_mode_estimation = 0;
else
end
return
function posterior_sampler_options = add_fields_(posterior_sampler_options, sampler_options)
fnam = fieldnames(sampler_options);
for j=1:length(fnam)
posterior_sampler_options.(fnam{j}) = sampler_options.(fnam{j});
end

9
matlab/dynare_estimation_1.m
View File

@ -420,15 +420,16 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if options_.mh_replic if options_.mh_replic
ana_deriv_old = options_.analytic_derivation; ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0; options_.analytic_derivation = 0;
posterior_sampler_options = options_.posterior_sampler_options; posterior_sampler_options = options_.posterior_sampler_options.current_options;
posterior_sampler_options.bounds = bounds;
posterior_sampler_options.invhess = invhess; posterior_sampler_options.invhess = invhess;
[posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_); [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_);
if strcmpi(options_.posterior_sampling_method,'adaptive_metropolis_hastings'), % keep old form only for this ... % store current options in global
options_.posterior_sampler_options.current_options = posterior_sampler_options;
if strcmpi(options_.posterior_sampler_options.posterior_sampling_method,'adaptive_metropolis_hastings'), % keep old form only for this ...
invhess = posterior_sampler_options.invhess; invhess = posterior_sampler_options.invhess;
feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
else else
posterior_sampler(objective_function,options_.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
end end
options_.analytic_derivation = ana_deriv_old; options_.analytic_derivation = ana_deriv_old;
end end

27
matlab/dynare_estimation_init.m
View File

@ -591,28 +591,7 @@ else
end end
% slice posterior sampler does not require mode or hessian to run if options_.mh_replic
if strcmp(options_.posterior_sampling_method,'slice') [current_options, options_] = check_posterior_sampler_options([], options_, bounds);
options_.mh_posterior_mode_estimation=1; options_.posterior_sampler_options.current_options = current_options;
options_.posterior_sampler_options = set_default_option(options_.posterior_sampler_options,'slice_initialize_with_mode',0);
if options_.posterior_sampler_options.slice_initialize_with_mode
if (isequal(options_.mode_compute,0) && isempty(options_.mode_file) )
skipline()
disp('You have specified the option "slice_initialize_with_mode"')
disp('to initialize the slice sampler using mode information')
disp('but no mode file nor posterior maximization is selected,')
error('The option "slice_initialize_with_mode" is inconsistent with mode_compute=0 or empty mode_file.')
else
options_.mh_posterior_mode_estimation=0;
end
end
if any(isinf(bounds.lb)) || any(isinf(bounds.ub)),
skipline()
disp('some priors are unbounded and prior_trunc is set to zero')
error('The option "slice" is inconsistent with prior_trunc=0.')
end
end end

26
matlab/global_initialization.m
View File

@ -435,11 +435,6 @@ options_.MCMC_jumping_covariance='hessian';
options_.use_calibration_initialization = 0; options_.use_calibration_initialization = 0;
options_.endo_vars_for_moment_computations_in_estimation=[]; options_.endo_vars_for_moment_computations_in_estimation=[];
% Tailored Random Block Metropolis-Hastings
options_.TaRB.use_TaRB = 0;
options_.TaRB.mode_compute=4;
options_.TaRB.new_block_probability=0.25; %probability that next parameter belongs to new block
% Run optimizer silently % Run optimizer silently
options_.silent_optimizer=0; options_.silent_optimizer=0;
@ -460,15 +455,32 @@ options_.prior_trunc = 1e-10;
options_.smoother = 0; options_.smoother = 0;
options_.posterior_max_subsample_draws = 1200; options_.posterior_max_subsample_draws = 1200;
options_.sub_draws = []; options_.sub_draws = [];
options_.use_mh_covariance_matrix = 0; % options_.use_mh_covariance_matrix = 0;
options_.gradient_method = 2; %used by csminwel and newrat options_.gradient_method = 2; %used by csminwel and newrat
options_.gradient_epsilon = 1e-6; %used by csminwel and newrat options_.gradient_epsilon = 1e-6; %used by csminwel and newrat
options_.posterior_sampler_options = []; %extended set of options for individual posterior samplers options_.posterior_sampler_options.sampling_opt = []; %extended set of options for individual posterior samplers
% Random Walk Metropolis-Hastings
options_.posterior_sampler_options.posterior_sampling_method = 'random_walk_metropolis_hastings'; options_.posterior_sampler_options.posterior_sampling_method = 'random_walk_metropolis_hastings';
options_.posterior_sampler_options.rwmh.proposal_distribution = 'rand_multivariate_normal'; options_.posterior_sampler_options.rwmh.proposal_distribution = 'rand_multivariate_normal';
options_.posterior_sampler_options.rwmh.student_degrees_of_freedom = 3; options_.posterior_sampler_options.rwmh.student_degrees_of_freedom = 3;
options_.posterior_sampler_options.rwmh.use_mh_covariance_matrix=0;
% Tailored Random Block Metropolis-Hastings
options_.posterior_sampler_options.tarb.proposal_distribution = 'rand_multivariate_normal'; options_.posterior_sampler_options.tarb.proposal_distribution = 'rand_multivariate_normal';
options_.posterior_sampler_options.tarb.student_degrees_of_freedom = 3; options_.posterior_sampler_options.tarb.student_degrees_of_freedom = 3;
options_.posterior_sampler_options.tarb.mode_compute=4;
options_.posterior_sampler_options.tarb.new_block_probability=0.25; %probability that next parameter belongs to new block
% Slice
options_.posterior_sampler_options.slice.proposal_distribution = '';
options_.posterior_sampler_options.slice.rotated=0;
options_.posterior_sampler_options.slice.slice_initialize_with_mode=0;
options_.posterior_sampler_options.slice.use_mh_covariance_matrix=0;
options_.posterior_sampler_options.slice.WR=[];
options_.posterior_sampler_options.slice.mode_files=[];
options_.posterior_sampler_options.slice.initial_step_size=0.8;
% Independent Metropolis-Hastings
options_.posterior_sampler_options.imh.proposal_distribution = 'rand_multivariate_normal';
options_.posterior_sampler_options.imh.use_mh_covariance_matrix=0;
options_.trace_plot_ma = 200; options_.trace_plot_ma = 200;
options_.mh_autocorrelation_function_size = 30; options_.mh_autocorrelation_function_size = 30;
options_.plot_priors = 1; options_.plot_priors = 1;

14
matlab/initial_estimation_checks.m
View File

@ -53,20 +53,6 @@ if maximum_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+Est
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!']) error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
end end
if ~(strcmpi(DynareOptions.proposal_distribution, 'rand_multivariate_student') || ...
strcmpi(DynareOptions.proposal_distribution, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
end
if DynareOptions.student_degrees_of_freedom <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
end
if DynareOptions.TaRB.use_TaRB && (DynareOptions.TaRB.new_block_probability<0 || DynareOptions.TaRB.new_block_probability>1)
error(['initial_estimation_checks:: The tarb_new_block_probability must be between 0 and 1!'])
end
if (any(BayesInfo.pshape >0 ) && DynareOptions.mh_replic) && DynareOptions.mh_nblck<1 if (any(BayesInfo.pshape >0 ) && DynareOptions.mh_replic) && DynareOptions.mh_nblck<1
error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.']) error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.'])
end end

16
matlab/posterior_sampler.m
View File

@ -89,11 +89,6 @@ load_last_mh_history_file(MetropolisFolder, ModelName);
% on many cores). The mandatory variables for local/remote parallel % on many cores). The mandatory variables for local/remote parallel
% computing are stored in the localVars struct. % computing are stored in the localVars struct.
if options_.TaRB.use_TaRB
options_.silent_optimizer=1; %locally set optimizer to silent mode
sampler_options.posterior_sampling_method='tailored_random_block_metropolis_hastings';
end
localVars = struct('TargetFun', TargetFun, ... localVars = struct('TargetFun', TargetFun, ...
'ProposalFun', ProposalFun, ... 'ProposalFun', ProposalFun, ...
'xparam1', xparam1, ... 'xparam1', xparam1, ...
@ -120,6 +115,9 @@ localVars = struct('TargetFun', TargetFun, ...
'oo_', oo_,... 'oo_', oo_,...
'varargin',[]); 'varargin',[]);
if strcmp(sampler_options.posterior_sampling_method,'tailored_random_block_metropolis_hastings');
localVars.options_.silent_optimizer=1; %locally set optimizer to silent mode
end
% User doesn't want to use parallel computing, or wants to compute a % User doesn't want to use parallel computing, or wants to compute a
% single chain compute Random walk Metropolis-Hastings algorithm sequentially. % single chain compute Random walk Metropolis-Hastings algorithm sequentially.
@ -145,11 +143,7 @@ else
end end
% from where to get back results % from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
if options_.TaRB.use_TaRB [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info);
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'TaRB_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info);
end
for j=1:totCPU, for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1; offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j))); record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
@ -190,4 +184,4 @@ if max(record.FunctionEvalPerIteration)>1
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])
end end
end end
end end

4
matlab/posterior_sampler_initialization.m
View File

@ -202,7 +202,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(['Estimation::mcmc: Write details about the MCMC... ']); fprintf(['Estimation::mcmc: Write details about the MCMC... ']);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns); AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns; AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
record.Sampler = options_.posterior_sampling_method; record.Sampler = options_.posterior_sampler_options.posterior_sampling_method;
record.Nblck = NumberOfBlocks; record.Nblck = NumberOfBlocks;
record.MhDraws = zeros(1,3); record.MhDraws = zeros(1,3);
record.MhDraws(1,1) = nruns(1); record.MhDraws(1,1) = nruns(1);
@ -446,4 +446,4 @@ elseif options_.mh_recover
end end
write_mh_history_file(MetropolisFolder, ModelName, record); write_mh_history_file(MetropolisFolder, ModelName, record);
end end
end end

157
tests/estimation/fs2000_slice.mod
View File

@ -1,76 +1,81 @@
// See fs2000.mod in the examples/ directory for details on the model // See fs2000.mod in the examples/ directory for details on the model
var m P c e W R k d n l gy_obs gp_obs y dA; var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m; varexo e_a e_m;
parameters alp bet gam mst rho psi del; parameters alp bet gam mst rho psi del;
alp = 0.33; alp = 0.33;
bet = 0.99; bet = 0.99;
gam = 0.003; gam = 0.003;
mst = 1.011; mst = 1.011;
rho = 0.7; rho = 0.7;
psi = 0.787; psi = 0.787;
del = 0.02; del = 0.02;
model; model;
dA = exp(gam+e_a); dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; -P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n; W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; -(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; 1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m; P*c = m;
m-1+d = l; m-1+d = l;
e = exp(e_a); e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1); gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA; gp_obs = (P/P(-1))*m(-1)/dA;
end; end;
initval; initval;
k = 6; k = 6;
m = mst; m = mst;
P = 2.25; P = 2.25;
c = 0.45; c = 0.45;
e = 1; e = 1;
W = 4; W = 4;
R = 1.02; R = 1.02;
d = 0.85; d = 0.85;
n = 0.19; n = 0.19;
l = 0.86; l = 0.86;
y = 0.6; y = 0.6;
gy_obs = exp(gam); gy_obs = exp(gam);
gp_obs = exp(-gam); gp_obs = exp(-gam);
dA = exp(gam); dA = exp(gam);
end; end;
shocks; shocks;
var e_a; stderr 0.014; var e_a; stderr 0.014;
var e_m; stderr 0.005; var e_m; stderr 0.005;
end; end;
steady; steady;
check; check;
estimated_params; estimated_params;
alp, beta_pdf, 0.356, 0.02; alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002; bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003; gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007; mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223; rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05; psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005; del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf; stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf; stderr e_m, inv_gamma_pdf, 0.008862, inf;
end; end;
varobs gp_obs gy_obs; varobs gp_obs gy_obs;
options_.posterior_sampling_method = 'slice'; //options_.posterior_sampling_method = 'slice';
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=50,mh_nblocks=2,mh_drop=0.2,mode_compute=0,cova_compute=0); estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=50,mh_nblocks=2,mh_drop=0.2, //mode_compute=0,cova_compute=0,
// continue with rotated slice posterior_sampling_method='slice'
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=100,mh_nblocks=2,mh_drop=0.5,load_mh_file,mode_compute=0,optim=('rotated',1,'use_mh_covariance_matrix',1)); );
// continue with rotated slice
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=100,mh_nblocks=2,mh_drop=0.5,load_mh_file,//mode_compute=0,
posterior_sampling_method='slice',
posterior_sampler_options=('rotated',1,'use_mh_covariance_matrix',1)
);