Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Show changes
Commits on Source (22)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1903 additions and 63 deletions
matlab/GetPosteriorParametersStatistics.m
View file @ f3843b49
......@@ -72,7 +72,7 @@ skipline()
try
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
catch
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
end
num_draws=NumberOfDraws*options_.mh_nblck;
......
matlab/McMCDiagnostics.m
View file @ f3843b49
......@@ -69,6 +69,46 @@ if issue_an_error_message
error('Estimation::mcmc::diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...')
end
% compute inefficiency factor
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
param_name=[];
for jj=1:npar
param_name = strvcat(param_name,get_the_name(jj,options_.TeX,M_,estim_params_,options_));
end
fprintf('\nMCMC Inefficiency factors per block.\n');
IFAC_header={'Parameter'};
print_string=['%',num2str(size(param_name,2)+3),'s'];
print_string_header=['%',num2str(size(param_name,2)+3),'s'];
for j=1:nblck,
IFAC_header=[IFAC_header, {['block ' int2str(j)]}];
print_string=[print_string,' \t %12.3f'];
print_string_header=[print_string_header,' \t %12s'];
end
print_string=[print_string,'\n'];
print_string_header=[print_string_header,'\n'];
fprintf(print_string_header,IFAC_header{1,:});
for jj=1:npar
Draws = GetAllPosteriorDraws(jj,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
Draws = reshape(Draws,[NumberOfDraws nblck]);
Nc = min(1000, NumberOfDraws/2);
for ll=1:nblck
Ifac(ll,jj) = mcmc_ifac(Draws(:,ll), Nc);
end
tmp = num2cell(Ifac(:,jj));
fprintf(print_string,param_name(jj,:),tmp{:})
end
skipline()
record.InefficiencyFactorsPerBlock = Ifac;
update_last_mh_history_file(MetropolisFolder, ModelName, record);
PastDraws = sum(record.MhDraws,1);
LastFileNumber = PastDraws(2);
LastLineNumber = record.MhDraws(end,3);
......@@ -436,4 +476,5 @@ if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fprintf(fidTeX,'\n');
fprintf(fidTeX,'% End Of TeX file.');
fclose(fidTeX);
end
\ No newline at end of file
end
matlab/check_posterior_sampler_options.m 0 → 100644
View file @ f3843b49
function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_)
% function posterior_sampler_options = check_posterior_sampler_options(posterior_sampler_options, options_)
% initialization of posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampling_method;
posterior_sampler_options.proposal_distribution = options_.proposal_distribution;
bounds = posterior_sampler_options.bounds;
invhess = posterior_sampler_options.invhess;
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
if ~options_.cova_compute
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.')
end
if options_.load_mh_file && options_.use_mh_covariance_matrix,
[junk, invhess] = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='MH';
posterior_sampler_options.serial_bar_title='Metropolis-Hastings';
posterior_sampler_options.save_tmp_file=1;
end
% check specific options for slice sampler
if strcmp(options_.posterior_sampling_method,'slice')
% by default, slice sampler should trigger
% mode_compute=0 and
% mh_replic=100 (much smaller than the default mh_replic=20000 of RWMH)
% moreover slice must be associated to:
% options_.mh_posterior_mode_estimation = 0;
% this is done below, but perhaps preprocessing should do this?
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;
posterior_sampler_options = set_default_option(posterior_sampler_options,'rotated',0);
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',[]);
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 = set_default_option(posterior_sampler_options,'W1',0.8*(bounds.ub-bounds.lb));
if ~isempty(options_.optim_opt)
options_list = read_key_value_string(options_.optim_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
posterior_sampler_options.W1 = options_list{i,2}*(bounds.ub-bounds.lb);
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};
otherwise
warning(['slice_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
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 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('I cannot start rotated slice sampler because')
disp('there is no previous MCMC to load ')
disp('or the Hessian at the mode is not computed.')
error('Rotated slice cannot start')
end
if isempty(invhess)
error('oops! This error should not occur, please contact developers.')
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
[V1 D]=eig(invhess);
posterior_sampler_options.V1=V1;
posterior_sampler_options.WR=sqrt(diag(D))*3;
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
options_.mh_posterior_mode_estimation = 0;
end
matlab/dynare_config.m
View file @ f3843b49
......@@ -83,6 +83,10 @@ if ~isoctave
% Replacements for functions of the stats toolbox
addpath([dynareroot '/missing/stats/'])
end
if isempty(ver('econ')),
% Replacements for functions of the econ toolbox
addpath([dynareroot '/missing/econ/'])
end
end
% ordeig() doesn't exist in Octave
......
matlab/dynare_estimation_1.m
View file @ f3843b49
......@@ -352,6 +352,7 @@ end
oo_.posterior.optimization.mode = xparam1;
oo_.posterior.optimization.Variance = [];
invhess=[];
if ~options_.mh_posterior_mode_estimation
if options_.cova_compute
invhess = inv(hh);
......@@ -366,6 +367,8 @@ else
xparam1 = bayestopt_.p5;
idNaN = isnan(xparam1);
xparam1(idNaN) = bayestopt_.p1(idNaN);
outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub);
xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars);
end
if ~options_.cova_compute
......@@ -415,30 +418,32 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end
% runs MCMC
if options_.mh_replic
if options_.load_mh_file && options_.use_mh_covariance_matrix
invhess = compute_mh_covariance_matrix;
end
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
if options_.cova_compute
posterior_sampler_options = options_.posterior_sampler_options;
posterior_sampler_options.bounds = bounds;
posterior_sampler_options.invhess = invhess;
[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 ...
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_);
else
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.')
posterior_sampler(objective_function,options_.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
end
options_.analytic_derivation = ana_deriv_old;
end
%% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_);
if options_.mh_posterior_mode_estimation
CutSample(M_, options_, estim_params_);
return
else
if ~options_.nodiagnostic && options_.mh_replic>0
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
end
%% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_);
%% Estimation of the marginal density from the Mh draws:
if options_.mh_replic
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
% Store posterior statistics by parameter name
oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames);
if ~options_.nograph
......
matlab/dynare_estimation_init.m
View file @ f3843b49
......@@ -590,3 +590,29 @@ else
estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=Sigma_e_non_zero_rows;
end
% slice posterior sampler does not require mode or hessian to run
if strcmp(options_.posterior_sampling_method,'slice')
options_.mh_posterior_mode_estimation=1;
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
matlab/global_initialization.m
View file @ f3843b49
......@@ -458,15 +458,17 @@ options_.prefilter = 0;
options_.presample = 0;
options_.prior_trunc = 1e-10;
options_.smoother = 0;
options_.student_degrees_of_freedom = 3;
options_.posterior_max_subsample_draws = 1200;
options_.sub_draws = [];
options_.use_mh_covariance_matrix = 0;
options_.gradient_method = 2; %used by csminwel and newrat
options_.gradient_epsilon = 1e-6; %used by csminwel and newrat
options_.posterior_sampling_method = 'random_walk_metropolis_hastings';
options_.proposal_distribution = 'rand_multivariate_normal';
options_.student_degrees_of_freedom = 3;
options_.posterior_sampler_options = []; %extended set of options for individual posterior samplers
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.student_degrees_of_freedom = 3;
options_.posterior_sampler_options.tarb.proposal_distribution = 'rand_multivariate_normal';
options_.posterior_sampler_options.tarb.student_degrees_of_freedom = 3;
options_.trace_plot_ma = 200;
options_.mh_autocorrelation_function_size = 30;
options_.plot_priors = 1;
......
matlab/marginal_density.m
View file @ f3843b49
function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_)
function [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_)
% function marginal = marginal_density()
% Computes the marginal density
%
......@@ -58,6 +58,9 @@ lpost_mode = posterior_kernel_at_the_mode;
xparam1 = posterior_mean;
hh = inv(SIGMA);
fprintf(' Done!\n');
if ~isfield(oo_,'posterior_mode')
oo_=fill_mh_mode(posterior_mode',NaN(npar,1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
end
% save the posterior mean and the inverse of the covariance matrix
% (usefull if the user wants to perform some computations using
......@@ -120,4 +123,87 @@ while check_coverage
end
end
oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));
\ No newline at end of file
oo_.MarginalDensity.ModifiedHarmonicMean = mean(marginal(:,2));
return
function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
%function oo_=fill_mh_mode(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_, field_name)
%
% INPUTS
% o xparam1 [double] (p*1) vector of estimate parameters.
% o stdh [double] (p*1) vector of estimate parameters.
% o M_ Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o estim_params_ Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o options_ Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o bayestopt_ Matlab's structure describing the priors (initialized by dynare, see @ref{bayesopt_}).
% o oo_ Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%
% OUTPUTS
% o oo_ Matlab's structure gathering the results
%
% SPECIAL REQUIREMENTS
% None.
nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters).
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters).
np = estim_params_.np ; % Number of deep parameters.
nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
if np
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
name = bayestopt_.name{ip};
eval(['oo_.' field_name '_mode.parameters.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.parameters.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if nvx
ip = 1;
for i=1:nvx
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
eval(['oo_.' field_name '_mode.shocks_std.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.shocks_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if nvn
ip = nvx+1;
for i=1:nvn
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
eval(['oo_.' field_name '_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.measurement_errors_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
if ncx
ip = nvx+nvn+1;
for i=1:ncx
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
eval(['oo_.' field_name '_mode.shocks_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.shocks_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
if ncn
ip = nvx+nvn+ncx+1;
for i=1:ncn
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
eval(['oo_.' field_name '_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.' field_name '_std_at_mode.measurement_errors_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
return
\ No newline at end of file
matlab/mcmc_ifac.m 0 → 100644
View file @ f3843b49
function Ifac = mcmc_ifac(X, Nc)
% function Ifac = mcmc_ifac(X, Nc)
% Compute inefficiency factor of a MCMC sample X
%
% INPUTS
% X: time series
% Nc: # of lags
%
% OUTPUTS
% Ifac: inefficiency factor of MCMC sample
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
Nc = floor(min(Nc, length(X)/2));
if mod(Nc,2),
Nc=Nc-1;
end
AcorrXSIM = autocorr(X(:), Nc);
%
%Calculate the Parzen Weight
Parzen=zeros(Nc+1,1);
for i=1: Nc/2+1
Parzen(i)=1 - 6*(i/Nc)^2+ 6*(i/Nc)^3;
end
for i=(Nc/2)+1: Nc+1
Parzen(i)=2 * (1-(i/Nc))^3;
end
Parzen=Parzen';
Ifac= 1+2*sum(Parzen(:).* AcorrXSIM);
matlab/missing/econ/autocorr.m 0 → 100644
View file @ f3843b49
function acf = autocorr(y, ar)
% function acf = autocorr(y, ar)
% autocorrelation function of y
%
% INPUTS
% y: time series
% ar: # of lags
%
% OUTPUTS
% acf: autocorrelation for lags 1 to ar
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
y=y(:);
acf = NaN(ar+1,1);
acf(1)=1;
m = mean(y);
sd = std(y,1);
for i=1:ar,
acf(i+1) = (y(i+1:end)-m)'*(y(1:end-i)-m)./((size(y,1))*sd^2);
end
matlab/posterior_sampler.m 0 → 100644
View file @ f3843b49
function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% Random Walk Metropolis-Hastings algorithm.
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
% function (posterior kernel).
% o ProposalFun [char] string specifying the name of the proposal
% density
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% o sampler_options structure
% .invhess [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 dataset_ data structure
% o dataset_info dataset info structure
% o options_ options structure
% o M_ model structure
% o estim_params_ estimated parameters structure
% o bayestopt_ estimation options structure
% o oo_ outputs structure
%
% ALGORITHM
% Random-Walk Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS
% None.
%
% PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed
% in parallel. The code suitable to be executed in
% parallel on multi core or cluster machine (in general a 'for' cycle)
% has been removed from this function and been placed in the random_walk_metropolis_hastings_core.m funtion.
%
% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
% to manage the parallel computations.
%
% 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 functions.
% Copyright (C) 2006-2015 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/>.
% In Metropolis, we set penalty to Inf so as to reject all parameter sets triggering an error during target density computation
global objective_function_penalty_base
objective_function_penalty_base = Inf;
vv = sampler_options.invhess;
% Initialization of the random walk metropolis-hastings chains.
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
% Load last mh history file
load_last_mh_history_file(MetropolisFolder, ModelName);
% Only for test parallel results!!!
% To check the equivalence between parallel and serial computation!
% First run in serial mode, and then comment the follow line.
% save('recordSerial.mat','-struct', 'record');
% For parallel runs after serial runs with the abobe line active.
% TempRecord=load('recordSerial.mat');
% record.Seeds=TempRecord.Seeds;
% Snapshot of the current state of computing. It necessary for the parallel
% execution (i.e. to execute in a corretct way a portion of code remotely or
% on many cores). The mandatory variables for local/remote parallel
% 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, ...
'ProposalFun', ProposalFun, ...
'xparam1', xparam1, ...
'vv', vv, ...
'sampler_options', sampler_options, ...
'mh_bounds', mh_bounds, ...
'ix2', ix2, ...
'ilogpo2', ilogpo2, ...
'ModelName', ModelName, ...
'fline', fline, ...
'npar', npar, ...
'nruns', nruns, ...
'NewFile', NewFile, ...
'MAX_nruns', MAX_nruns, ...
'd', d, ...
'InitSizeArray',InitSizeArray, ...
'record', record, ...
'dataset_', dataset_, ...
'dataset_info', dataset_info, ...
'options_', options_, ...
'M_',M_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_,...
'varargin',[]);
% User doesn't want to use parallel computing, or wants to compute a
% single chain compute Random walk Metropolis-Hastings algorithm sequentially.
if isnumeric(options_.parallel) || (nblck-fblck)==0,
fout = posterior_sampler_core(localVars, fblck, nblck, 0);
record = fout.record;
% Parallel in Local or remote machine.
else
% Global variables for parallel routines.
globalVars = struct();
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[ModelName '_static.m']};
NamFileInput(2,:) = {'',[ModelName '_dynamic.m']};
if options_.steadystate_flag,
NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']};
end
if (options_.load_mh_file~=0) && any(fline>1) ,
NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']};
end
if exist([ModelName '_optimal_mh_scale_parameter.mat'])
NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_optimal_mh_scale_parameter.mat']};
end
% from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
if options_.TaRB.use_TaRB
[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,
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.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)));
record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)));
record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)));
end
end
irun = fout(1).irun;
NewFile = fout(1).NewFile;
record.MCMCConcludedSuccessfully = 1; %set indicator for successful run
update_last_mh_history_file(MetropolisFolder, ModelName, record);
% Provide diagnostic output
skipline()
disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.'])
disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.'])
disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
disp(['Estimation::mcmc: Current acceptance ratio per chain: '])
for i=1:nblck
if i<10
disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%'])
else
disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%'])
end
end
if max(record.FunctionEvalPerIteration)>1
disp(['Estimation::mcmc: Current function evaluations per iteration: '])
for i=1:nblck
if i<10
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])
else
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])
end
end
end
\ No newline at end of file
matlab/posterior_sampler_core.m 0 → 100644
View file @ f3843b49
function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% Contains the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for'
% cycle are completely independent to be suitable for parallel execution.
%
% INPUTS
% o myimput [struc] The mandatory variables for local/remote
% parallel computing obtained from random_walk_metropolis_hastings.m
% function.
% o fblck and nblck [integer] The Metropolis-Hastings chains.
% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed.
% The integer whoaim is the integer that
% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
% cluster.
% o ThisMatlab [integer] Allows us to distinguish between the
% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab,
% ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS
% o myoutput [struc]
% If executed without parallel, this is the original output of 'for b =
% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or
% remote machine. In this case:
% record;
% irun;
% NewFile;
% OutputFileName
%
% ALGORITHM
% Portion of Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS.
% None.
%
% PARALLEL CONTEXT
% See the comments in the random_walk_metropolis_hastings.m funtion.
% Copyright (C) 2006-2015 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/>.
if nargin<4,
whoiam=0;
end
% reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
TargetFun=myinputs.TargetFun;
ProposalFun=myinputs.ProposalFun;
xparam1=myinputs.xparam1;
mh_bounds=myinputs.mh_bounds;
last_draw=myinputs.ix2;
last_posterior=myinputs.ilogpo2;
fline=myinputs.fline;
npar=myinputs.npar;
nruns=myinputs.nruns;
NewFile=myinputs.NewFile;
MAX_nruns=myinputs.MAX_nruns;
sampler_options=myinputs.sampler_options;
d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
dataset_ = myinputs.dataset_;
dataset_info = myinputs.dataset_info;
bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
M_ = myinputs.M_;
oo_ = myinputs.oo_;
% Necessary only for remote computing!
if whoiam
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
end
MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
save_tmp_file = sampler_options.save_tmp_file;
options_.lik_algo = 1;
OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
sampler_options.n = npar;
sampler_options.ProposalDensity = 'multivariate_normal_pdf';
elseif strcmpi(ProposalFun,'rand_multivariate_student')
sampler_options.n = options_.student_degrees_of_freedom;
sampler_options.ProposalDensity = 'multivariate_student_pdf';
end
%
% Now I run the (nblck-fblck+1) MCMC chains
%
sampler_options.xparam1 = xparam1;
if ~isempty(d),
sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
end
block_iter=0;
for curr_block = fblck:nblck,
LastSeeds=[];
block_iter=block_iter+1;
try
% This will not work if the master uses a random number generator not
% available in the slave (different Matlab version or
% Matlab/Octave cluster). Therefore the trap.
%
% Set the random number generator type (the seed is useless but needed by the function)
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
% Set the state of the RNG
set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
catch
% If the state set by master is incompatible with the slave, we only reseed
set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
end
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
OpenOldFile(curr_block) = 0;
else
x2 = zeros(InitSizeArray(curr_block),npar);
logpo2 = zeros(InitSizeArray(curr_block),1);
end
%Prepare waiting bars
if whoiam
refresh_rate = sampler_options.parallel_bar_refresh_rate;
bar_title = sampler_options.parallel_bar_title;
prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
else
refresh_rate = sampler_options.serial_bar_refresh_rate;
bar_title = sampler_options.serial_bar_title;
hh = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
set(hh,'Name',bar_title);
end
accepted_draws_this_chain = 0;
accepted_draws_this_file = 0;
feval_this_chain = 0;
feval_this_file = 0;
draw_index_current_file = fline(curr_block); %get location of first draw in current block
draw_iter = 1;
while draw_iter <= nruns(curr_block)
[par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
x2(draw_index_current_file,:) = par;
last_draw(curr_block,:) = par;
logpo2(draw_index_current_file) = logpost;
last_posterior(curr_block) = logpost;
feval_this_chain = feval_this_chain + sum(neval);
feval_this_file = feval_this_file + sum(neval);
accepted_draws_this_chain = accepted_draws_this_chain + accepted;
accepted_draws_this_file = accepted_draws_this_file + accepted;
prtfrc = draw_iter/nruns(curr_block);
if mod(draw_iter, refresh_rate)==0
if accepted_draws_this_chain/draw_iter==1 && sum(neval)>1
dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Function eval per draw %4.3f', feval_this_chain/draw_iter)]);
else
dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]);
end
if save_tmp_file
[LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
save([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds');
end
end
if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done
[LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
if save_tmp_file,
delete([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
end
save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);
fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
fprintf(fidlog,[' Acceptance ratio......: ' num2str(accepted_draws_this_file/length(logpo2)) '\n']);
fprintf(fidlog,[' Feval per iteration...: ' num2str(feval_this_file/length(logpo2)) '\n']);
fprintf(fidlog,[' Posterior mean........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']);
fprintf(fidlog,[' Minimum value.........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']);
fprintf(fidlog,[' Maximum value.........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
fprintf(fidlog,' \n');
fclose(fidlog);
accepted_draws_this_file = 0;
feval_this_file = 0;
if draw_iter == nruns(curr_block) % I record the last draw...
record.LastParameters(curr_block,:) = x2(end,:);
record.LastLogPost(curr_block) = logpo2(end);
end
% size of next file in chain curr_block
InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns);
% initialization of next file if necessary
if InitSizeArray(curr_block)
x2 = zeros(InitSizeArray(curr_block),npar);
logpo2 = zeros(InitSizeArray(curr_block),1);
NewFile(curr_block) = NewFile(curr_block) + 1;
draw_index_current_file = 0;
end
end
draw_iter=draw_iter+1;
draw_index_current_file = draw_index_current_file + 1;
end% End of the simulations for one mh-block.
record.AcceptanceRatio(curr_block) = accepted_draws_this_chain/(draw_iter-1);
record.FunctionEvalPerIteration(curr_block) = feval_this_chain/(draw_iter-1);
dyn_waitbar_close(hh);
[record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state();
OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']};
end% End of the loop over the mh-blocks.
myoutput.record = record;
myoutput.irun = draw_index_current_file;
myoutput.NewFile = NewFile;
myoutput.OutputFileName = OutputFileName;
\ No newline at end of file
matlab/posterior_sampler_initialization.m 0 → 100644
View file @ f3843b49
This diff is collapsed.
matlab/posterior_sampler_iteration.m 0 → 100644
View file @ f3843b49
function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,varargin)
% function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
% posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
posterior_sampling_method = sampler_options.posterior_sampling_method;
mh_bounds = sampler_options.bounds;
switch posterior_sampling_method
case 'slice'
[par, logpost, neval] = slice_sampler(TargetFun,last_draw, [mh_bounds.lb mh_bounds.ub], sampler_options,varargin{:});
accepted = 1;
case 'random_walk_metropolis_hastings'
neval = 1;
ProposalFun = sampler_options.proposal_distribution;
proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition;
n = sampler_options.n;
par = feval(ProposalFun, last_draw, proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
try
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
else
logpost = -inf;
end
r = logpost-last_posterior;
if (logpost > -inf) && (log(rand) < r)
accepted = 1;
else
accepted = 0;
par = last_draw;
logpost = last_posterior;
end
case 'tailored_random_block_metropolis_hastings'
options_=varargin{3};
bayestopt_=varargin{6};
npar=length(last_draw);
%% randomize indices for blocking in this iteration
global objective_function_penalty_base
indices=randperm(npar)';
blocks=[1; (1+cumsum((rand(length(indices)-1,1)>(1-options_.TaRB.new_block_probability))))];
nblocks=blocks(end,1); %get number of blocks this iteration
current_draw=last_draw'; %get starting point for current draw for updating
blocked_draws_counter=0;
accepted_draws_counter=0;
for block_iter=1:nblocks
blocked_draws_counter=blocked_draws_counter+1;
nxopt=length(indices(blocks==block_iter,1)); %get size of current block
par_start_current_block=current_draw(indices(blocks==block_iter,1));
[xopt_current_block, fval, exitflag, hess_mat_optimizer, options_, Scale] = dynare_minimize_objective(@TaRB_optimizer_wrapper,par_start_current_block,options_.TaRB.mode_compute,options_,[mh_bounds.lb(indices(blocks==block_iter,1),1) mh_bounds.ub(indices(blocks==block_iter,1),1)],bayestopt_.name,bayestopt_,[],...
current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper
varargin{:}); %inputs for objective
objective_function_penalty_base=Inf; %reset penalty that may have been changed by optimizer
%% covariance for proposal density
hessian_mat = reshape(hessian('TaRB_optimizer_wrapper',xopt_current_block, ...
options_.gstep,...
current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper
varargin{:}),nxopt,nxopt);
if any(any(isnan(hessian_mat))) || any(any(isinf(hessian_mat)))
inverse_hessian_mat=eye(nxopt)*1e-4; %use diagonal
else
inverse_hessian_mat=inv(hessian_mat); %get inverse Hessian
if any(any((isnan(inverse_hessian_mat)))) || any(any((isinf(inverse_hessian_mat))))
inverse_hessian_mat=eye(nxopt)*1e-4; %use diagonal
end
end
[proposal_covariance_Cholesky_decomposition_upper,negeigenvalues]=chol(inverse_hessian_mat);
%if not positive definite, use generalized Cholesky of Eskow/Schnabel
if negeigenvalues~=0
proposal_covariance_Cholesky_decomposition_upper=chol_SE(inverse_hessian_mat,0);
end
proposal_covariance_Cholesky_decomposition_upper=proposal_covariance_Cholesky_decomposition_upper*diag(bayestopt_.jscale(indices(blocks==block_iter,1),:));
%get proposal draw
if strcmpi(sampler_options.proposal_distribution,'rand_multivariate_normal')
n = nxopt;
elseif strcmpi(sampler_options.proposal_distribution,'rand_multivariate_student')
n = options_.student_degrees_of_freedom;
end
proposed_par = feval(sampler_options.proposal_distribution, xopt_current_block', proposal_covariance_Cholesky_decomposition_upper, n);
% check whether draw is valid and compute posterior
if all( proposed_par(:) > mh_bounds.lb(indices(blocks==block_iter,1),:) ) && all( proposed_par(:) < mh_bounds.ub(indices(blocks==block_iter,1),:) )
try
logpost = - feval('TaRB_optimizer_wrapper', proposed_par(:),...
current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper
varargin{:});
catch
logpost = -inf;
end
else
logpost = -inf;
end
%get ratio of proposal densities, required because proposal depends
%on current mode via Hessian and is thus not symmetric anymore
if strcmpi(sampler_options.proposal_distribution,'rand_multivariate_normal')
proposal_density_proposed_move_forward=multivariate_normal_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
proposal_density_proposed_move_backward=multivariate_normal_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
elseif strcmpi(sampler_options.proposal_distribution,'rand_multivariate_student')
proposal_density_proposed_move_forward=multivariate_student_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
proposal_density_proposed_move_backward=multivariate_student_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n);
end
accprob=logpost-last_posterior+ log(proposal_density_proposed_move_backward)-log(proposal_density_proposed_move_forward); %Formula (6), Chib/Ramamurthy
if (logpost > -inf) && (log(rand) < accprob)
current_draw(indices(blocks==block_iter,1))=proposed_par;
last_posterior=logpost;
accepted_draws_counter =accepted_draws_counter +1;
else %no updating
%do nothing, keep old value
end
end
accepted=accepted_draws_counter/blocked_draws_counter;
par = current_draw;
neval=1;
case 'independent_metropolis_hastings'
neval = 1;
ProposalFun = sampler_options.proposal_distribution;
ProposalDensity = sampler_options.ProposalDensity;
proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition;
n = sampler_options.n;
xparam1 = sampler_options.xparam1;
par = feval(ProposalFun, sampler_options.xparam1, proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
try
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
else
logpost = -inf;
end
r = logpost - last_posterior + ...
log(feval(ProposalDensity, last_draw, xparam1, proposal_covariance_Cholesky_decomposition, n)) - ...
log(feval(ProposalDensity, par, xparam1, proposal_covariance_Cholesky_decomposition, n));
if (logpost > -inf) && (log(rand) < r)
accepted = 1;
else
accepted = 0;
par = last_draw;
logpost = last_posterior;
end
end
\ No newline at end of file
matlab/rotated_slice_sampler.m 0 → 100644
View file @ f3843b49
function [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
% ----------------------------------------------------------
% ROTATED SLICE SAMPLER - with stepping out (Neal, 2003)
% extension of the orthogonal univarite sampler (slice_sampler.m)
% copyright M. Ratto (European Commission)
%
% objective_function(theta,varargin): -log of any unnormalized pdf
% with varargin (optional) a vector of auxiliaty parameters
% to be passed to f( ).
% ----------------------------------------------------------
%
% INPUTS
% objective_function: objective function (expressed as minus the log of a density)
% theta: last value of theta
% thetaprior: bounds of the theta space
% sampler_options: posterior sampler options
% varargin: optional input arguments to objective function
%
% OUTPUTS
% theta: new theta sample
% fxsim: value of the objective function for the new sample
% neval: number of function evaluations
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
theta=theta(:);
npar = length(theta);
neval = zeros(npar,1);
W1=[];
if isfield(sampler_options,'WR'),
W1 = sampler_options.WR;
end
if ~isempty(sampler_options.mode),
mm = sampler_options.mode;
n = length(mm);
for j=1:n,
distance(j)=sqrt(sum((theta-mm(j).m).^2));
end
[m, im] = min(distance);
r=im;
V1 = mm(r).m;
jj=0;
for j=1:n,
if j~=r,
jj=jj+1;
tmp=mm(j).m-mm(r).m;
%tmp=mm(j).m-theta;
V1(:,jj)=tmp/norm(tmp);
end
end
resul=randperm(n-1,n-1);
V1 = V1(:,resul);
%V1 = V1(:, randperm(n-1));
% %d = chol(mm(r).invhess);
% %V1 = transpose(feval(sampler_options.proposal_distribution, transpose(mm(r).m), d, npar));
%
% V1=eye(npar);
% V1=V1(:,randperm(npar));
% for j=1:2,
% V1(:,j)=mm(r(j)).m-theta;
% V1(:,j)=V1(:,j)/norm(V1(:,j));
% end
% % Gram-Schmidt
% for j=2:npar,
% for k=1:j-1,
% V1(:,j)=V1(:,j)-V1(:,k)'*V1(:,j)*V1(:,k);
% end
% V1(:,j)=V1(:,j)/norm(V1(:,j));
% end
% for j=1:n,
% distance(j)=sqrt(sum((theta-mm(j).m).^2));
% end
% [m, im] = min(distance);
% if im==r,
% fxsim=[];
% return,
% else
% theta1=theta;
% end
else
V1 = sampler_options.V1;
end
npar=size(V1,2);
for it=1:npar,
theta0 = theta;
neval(it) = 0;
xold = 0;
% XLB = thetaprior(3);
% XUB = thetaprior(4);
tb=sort([(thetaprior(:,1)-theta)./V1(:,it) (thetaprior(:,2)-theta)./V1(:,it)],2);
XLB=max(tb(:,1));
XUB=min(tb(:,2));
if isempty(W1),
W = (XUB-XLB); %*0.8;
else
W = W1(it);
end
% -------------------------------------------------------
% 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1))
% THIS DEFINES THE SLICE S={x: z < ln(f(x))}
% -------------------------------------------------------
fxold = -feval(objective_function,theta,varargin{:});
%I have to be sure that the rotation is for L,R or for Fxold, theta(it)
neval(it) = neval(it) + 1;
Z = fxold + log(rand(1,1));
% -------------------------------------------------------------
% 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE
% STEPPING-OUT PROCEDURE
% -------------------------------------------------------------
u = rand(1,1);
L = max(XLB,xold-W*u);
R = min(XUB,L+W);
%[L R]=slice_rotation(L, R, alpha);
while(L > XLB)
xsim = L;
theta = theta0+xsim*V1(:,it);
fxl = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (fxl <= Z)
break;
end
L = max(XLB,L-W);
end
while(R < XUB)
xsim = R;
theta = theta0+xsim*V1(:,it);
fxr = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (fxr <= Z)
break;
end
R = min(XUB,R+W);
end
% ------------------------------------------------------
% 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA)
% ------------------------------------------------------
fxsim = Z-1;
while (fxsim < Z)
u = rand(1,1);
xsim = L + u*(R - L);
theta = theta0+xsim*V1(:,it);
fxsim = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (xsim > xold)
R = xsim;
else
L = xsim;
end
end
end
% if ~isempty(sampler_options.mode),
% dist1=sqrt(sum((theta-mm(r).m).^2));
% if dist1>distance(r),
% theta=theta1;
% fxsim=[];
% end
% end
end
matlab/slice_sampler.m 0 → 100644
View file @ f3843b49
function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
% function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
% ----------------------------------------------------------
% UNIVARIATE SLICE SAMPLER - stepping out (Neal, 2003)
% W: optimal value in the range (3,10)*std(x)
% - see C.Planas and A.Rossi (2014)
% objective_function(theta,varargin): -log of any unnormalized pdf
% with varargin (optional) a vector of auxiliaty parameters
% to be passed to f( ).
% ----------------------------------------------------------
%
% INPUTS
% objective_function: objective function (expressed as minus the log of a density)
% theta: last value of theta
% thetaprior: bounds of the theta space
% sampler_options: posterior sampler options
% varargin: optional input arguments to objective function
%
% OUTPUTS
% theta: new theta sample
% fxsim: value of the objective function for the new sample
% neval: number of function evaluations
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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/>.
if sampler_options.rotated %&& ~isempty(sampler_options.V1),
[theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin{:});
if isempty(sampler_options.mode), % jumping
return,
else
nevalR=sum(neval);
end
end
theta=theta(:);
npar = length(theta);
W1 = sampler_options.W1;
neval = zeros(npar,1);
for it=1:npar,
neval(it) = 0;
W = W1(it);
xold = theta(it);
% XLB = thetaprior(3);
% XUB = thetaprior(4);
XLB = thetaprior(it,1);
XUB = thetaprior(it,2);
% -------------------------------------------------------
% 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1))
% THIS DEFINES THE SLICE S={x: z < ln(f(x))}
% -------------------------------------------------------
fxold = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
Z = fxold + log(rand(1,1));
% -------------------------------------------------------------
% 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE
% STEPPING-OUT PROCEDURE
% -------------------------------------------------------------
u = rand(1,1);
L = max(XLB,xold-W*u);
R = min(XUB,L+W);
while(L > XLB)
xsim = L;
theta(it) = xsim;
fxl = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (fxl <= Z)
break;
end
L = max(XLB,L-W);
end
while(R < XUB)
xsim = R;
theta(it) = xsim;
fxr = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (fxr <= Z)
break;
end
R = min(XUB,R+W);
end
% ------------------------------------------------------
% 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA)
% ------------------------------------------------------
fxsim = Z-1;
while (fxsim < Z)
u = rand(1,1);
xsim = L + u*(R - L);
theta(it) = xsim;
fxsim = -feval(objective_function,theta,varargin{:});
neval(it) = neval(it) + 1;
if (xsim > xold)
R = xsim;
else
L = xsim;
end
end
end
if sampler_options.rotated && ~isempty(sampler_options.mode), % jumping
neval=sum(neval)+nevalR;
end
preprocessor/DynareBison.yy
View file @ f3843b49
......@@ -83,14 +83,14 @@ class ParsingDriver;
#define yylex driver.lexer->lex
}
%token AIM_SOLVER ANALYTIC_DERIVATION AR AUTOCORR TARB_MODE_COMPUTE
%token BAYESIAN_IRF BETA_PDF BLOCK USE_CALIBRATION USE_TARB TARB_NEW_BLOCK_PROBABILITY SILENT_OPTIMIZER
%token AIM_SOLVER ANALYTIC_DERIVATION AR AUTOCORR POSTERIOR_SAMPLING_METHOD
%token BAYESIAN_IRF BETA_PDF BLOCK USE_CALIBRATION SILENT_OPTIMIZER
%token BVAR_DENSITY BVAR_FORECAST NODECOMPOSITION DR_DISPLAY_TOL HUGE_NUMBER
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA TARB_OPTIM
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION
%token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
%token CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED STUDENT_DEGREES_OF_FREEDOM
%token CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED
%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
%token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
......@@ -152,7 +152,7 @@ class ParsingDriver;
%token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD WEIBULL WEIBULL_PDF
%token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
%token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION
%token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS
%token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS POSTERIOR_SAMPLER_OPTIONS
%token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON
%token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION
%token MS_ESTIMATION MS_SIMULATION MS_COMPUTE_MDD MS_COMPUTE_PROBABILITIES MS_FORECAST
......@@ -180,7 +180,7 @@ class ParsingDriver;
%type <string_val> filename symbol vec_of_vec_value vec_value_list date_expr
%type <string_val> vec_value_1 vec_value signed_inf signed_number_w_inf
%type <string_val> range vec_value_w_inf vec_value_1_w_inf
%type <string_val> integer_range signed_integer_range
%type <string_val> integer_range signed_integer_range sub_sampling_options list_sub_sampling_option
%type <string_pair_val> named_var
%type <symbol_type_val> change_type_arg
%type <vector_string_val> change_type_var_list subsamples_eq_opt prior_eq_opt options_eq_opt calibration_range
......@@ -1757,14 +1757,11 @@ estimation_options : o_datafile
| o_distribution_approximation
| o_dirname
| o_huge_number
| o_use_tarb
| o_tarb_mode_compute
| o_tarb_new_block_probability
| o_tarb_optim
| o_silent_optimizer
| o_proposal_distribution
| o_student_degrees_of_freedom
| o_no_posterior_kernel_density
| o_posterior_sampling_method
| o_posterior_sampler_options
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
......@@ -1777,15 +1774,48 @@ optim_options : list_optim_option
| optim_options COMMA list_optim_option;
;
list_tarb_optim_option : QUOTED_STRING COMMA QUOTED_STRING
{ driver.tarb_optim_options_string($1, $3); }
| QUOTED_STRING COMMA signed_number
{ driver.tarb_optim_options_num($1, $3); }
;
list_sub_sampling_option : QUOTED_STRING COMMA QUOTED_STRING
{
$1->insert(0, "''");
$1->append("'', ''");
$1->append(*$3);
$1->append("''");
$$ = $1;
}
| QUOTED_STRING COMMA signed_number
{
$1->insert(0, "''");
$1->append("'',");
$1->append(*$3);
$$ = $1;
}
;
tarb_optim_options : list_tarb_optim_option
| tarb_optim_options COMMA list_tarb_optim_option;
;
sub_sampling_options : list_sub_sampling_option
{ $$ = $1; }
| sub_sampling_options COMMA list_sub_sampling_option
{
$1->append(",");
$1->append(*$3);
$$ = $1;
}
;
list_sampling_option : QUOTED_STRING COMMA QUOTED_STRING
{ driver.sampling_options_string($1, $3); }
| QUOTED_STRING COMMA signed_number
{ driver.sampling_options_num($1, $3); }
| QUOTED_STRING COMMA '(' sub_sampling_options ')'
{
$4->insert(0,"(");
$4->append(")");
driver.sampling_options_string($1, $4);
}
;
sampling_options : list_sampling_option
| sampling_options COMMA list_sampling_option;
;
varobs : VAROBS { driver.check_varobs(); } varobs_list ';';
......@@ -2679,6 +2709,8 @@ o_est_first_obs : FIRST_OBS EQUAL vec_int
| FIRST_OBS EQUAL vec_int_number
{ driver.option_vec_int("first_obs", $3); }
;
o_posterior_sampling_method : POSTERIOR_SAMPLING_METHOD EQUAL QUOTED_STRING
{ driver.option_str("posterior_sampler_options.posterior_sampling_method", $3); } ;
o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); };
o_data_first_obs : FIRST_OBS EQUAL date_expr { driver.option_date("firstobs", $3); } ;
o_data_last_obs : LAST_OBS EQUAL date_expr { driver.option_date("lastobs", $3); } ;
......@@ -2735,12 +2767,11 @@ o_posterior_max_subsample_draws : POSTERIOR_MAX_SUBSAMPLE_DRAWS EQUAL INT_NUMBER
o_mh_drop : MH_DROP EQUAL non_negative_number { driver.option_num("mh_drop", $3); };
o_mh_jscale : MH_JSCALE EQUAL non_negative_number { driver.option_num("mh_jscale", $3); };
o_optim : OPTIM EQUAL '(' optim_options ')';
o_tarb_optim : TARB_OPTIM EQUAL '(' tarb_optim_options ')';
o_proposal_distribution : PROPOSAL_DISTRIBUTION EQUAL symbol { driver.option_str("proposal_distribution", $3); };
o_posterior_sampler_options : POSTERIOR_SAMPLER_OPTIONS EQUAL '(' sampling_options ')' ;
o_proposal_distribution : PROPOSAL_DISTRIBUTION EQUAL symbol { driver.option_str("posterior_sampler_options.posterior_sampling_method.proposal_distribution", $3); };
o_no_posterior_kernel_density : NO_POSTERIOR_KERNEL_DENSITY
{ driver.option_num("moments_posterior_density.indicator", "0"); }
;
o_student_degrees_of_freedom : STUDENT_DEGREES_OF_FREEDOM EQUAL INT_NUMBER { driver.option_num("student_degrees_of_freedom", $3); };
o_mh_init_scale : MH_INIT_SCALE EQUAL non_negative_number { driver.option_num("mh_init_scale", $3); };
o_mode_file : MODE_FILE EQUAL filename { driver.option_str("mode_file", $3); };
o_mode_compute : MODE_COMPUTE EQUAL INT_NUMBER { driver.option_num("mode_compute", $3); };
......@@ -2999,9 +3030,6 @@ o_equations : EQUATIONS EQUAL vec_int
| EQUATIONS EQUAL vec_int_number
{ driver.option_vec_int("ms.equations",$3); }
;
o_use_tarb : USE_TARB { driver.option_num("TaRB.use_TaRB", "1"); };
o_tarb_mode_compute : TARB_MODE_COMPUTE EQUAL INT_NUMBER { driver.option_num("TaRB.mode_compute", $3); };
o_tarb_new_block_probability : TARB_NEW_BLOCK_PROBABILITY EQUAL non_negative_number {driver.option_num("TaRB.new_block_probability",$3); };
o_silent_optimizer : SILENT_OPTIMIZER { driver.option_num("silent_optimizer", "1"); };
o_instruments : INSTRUMENTS EQUAL '(' symbol_list ')' {driver.option_symbol_list("instruments"); };
......
preprocessor/DynareFlex.ll
View file @ f3843b49
......@@ -390,7 +390,6 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;}
<DYNARE_STATEMENT>proposal_distribution {return token::PROPOSAL_DISTRIBUTION;}
<DYNARE_STATEMENT>no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;}
<DYNARE_STATEMENT>student_degrees_of_freedom {return token::STUDENT_DEGREES_OF_FREEDOM;}
<DYNARE_STATEMENT>alpha {
yylval->string_val = new string(yytext);
......@@ -588,10 +587,8 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>outvars {return token::OUTVARS;}
<DYNARE_STATEMENT>huge_number {return token::HUGE_NUMBER;}
<DYNARE_STATEMENT>dr_display_tol {return token::DR_DISPLAY_TOL;}
<DYNARE_STATEMENT>use_tarb {return token::USE_TARB;}
<DYNARE_STATEMENT>tarb_mode_compute {return token::TARB_MODE_COMPUTE;}
<DYNARE_STATEMENT>tarb_new_block_probability {return token::TARB_NEW_BLOCK_PROBABILITY;}
<DYNARE_STATEMENT>tarb_optim {return token::TARB_OPTIM;}
<DYNARE_STATEMENT>posterior_sampling_method {return token::POSTERIOR_SAMPLING_METHOD;}
<DYNARE_STATEMENT>posterior_sampler_options {return token::POSTERIOR_SAMPLER_OPTIONS;}
<DYNARE_STATEMENT>silent_optimizer {return token::SILENT_OPTIMIZER;}
<DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
<DYNARE_STATEMENT>occbin {return token::OCCBIN;}
......
preprocessor/ParsingDriver.cc
View file @ f3843b49
......@@ -1655,29 +1655,30 @@ ParsingDriver::optim_options_num(string *name, string *value)
}
void
ParsingDriver::tarb_optim_options_helper(const string &name)
ParsingDriver::sampling_options_helper(const string &name)
{
if (options_list.string_options.find("TaRB.optim_opt") == options_list.string_options.end())
options_list.string_options["TaRB.optim_opt"] = "";
if (options_list.string_options.find("posterior_sampler_options.sampling_opt") ==
options_list.string_options.end())
options_list.string_options["posterior_sampler_options.sampling_opt"] = "";
else
options_list.string_options["TaRB.optim_opt"] += ",";
options_list.string_options["TaRB.optim_opt"] += "''" + name + "'',";
options_list.string_options["posterior_sampler_options.sampling_opt"] += ",";
options_list.string_options["posterior_sampler_options.sampling_opt"] += "''" + name + "'',";
}
void
ParsingDriver::tarb_optim_options_string(string *name, string *value)
ParsingDriver::sampling_options_string(string *name, string *value)
{
tarb_optim_options_helper(*name);
options_list.string_options["TaRB.optim_opt"] += "''" + *value + "''";
sampling_options_helper(*name);
options_list.string_options["posterior_sampler_options.sampling_opt"] += "''" + *value + "''";
delete name;
delete value;
}
void
ParsingDriver::tarb_optim_options_num(string *name, string *value)
ParsingDriver::sampling_options_num(string *name, string *value)
{
tarb_optim_options_helper(*name);
options_list.string_options["TaRB.optim_opt"] += *value;
sampling_options_helper(*name);
options_list.string_options["posterior_sampler_options.sampling_opt"] += *value;
delete name;
delete value;
}
......
preprocessor/ParsingDriver.hh
View file @ f3843b49
......@@ -98,7 +98,7 @@ private:
//! Creates option "optim_opt" in OptionsList if it doesn't exist, else add a comma, and adds the option name
void optim_options_helper(const string &name);
void tarb_optim_options_helper(const string &name);
void sampling_options_helper(const string &name);
//! Stores temporary symbol table
SymbolList symbol_list;
......@@ -450,10 +450,10 @@ public:
void optim_options_string(string *name, string *value);
//! Adds an optimization option (numeric value)
void optim_options_num(string *name, string *value);
//! Adds a TaRB optimization option (string value)
void tarb_optim_options_string(string *name, string *value);
//! Adds a TaRB optimization option (numeric value)
void tarb_optim_options_num(string *name, string *value);
//! Adds an sampling option (string value)
void sampling_options_string(string *name, string *value);
//! Adds an sampling option (numeric value)
void sampling_options_num(string *name, string *value);
//! Check that no observed variable has yet be defined
void check_varobs();
//! Add a new observed variable
......