diff --git a/README.md b/README.md
index 9d5ec4aeaaba24f97227a791f1a0e653488288d0..5edd2eadc2db551270e01afe5c71ac877a63303c 100644
--- a/README.md
+++ b/README.md
@@ -54,6 +54,7 @@ a 32-bit Octave.
1. [**General Instructions**](#general-instructions)
1. [**Debian or Ubuntu**](#debian-or-ubuntu)
1. [**Fedora, CentOS or RHEL**](#fedora-centos-or-rhel)
+1. [**Arch Linux**](#arch-linux)
1. [**Windows**](#windows)
1. [**macOS**](#macos)
1. [**Docker**](#docker)
@@ -264,7 +265,8 @@ make -f makefile.gf FFLAGS="-O2 -std=legacy" PROGRAM=x13as
sudo cp x13as /usr/bin/
```
-If you use MATLAB, we strongly advise to also rename or exclude the GCC libraries shipped with MATLAB to avoid possible conflicts with GCC libraries shipped by Fedora, see e.g. [Matlab on Fedora 33](https://mutschler.eu/linux/install-guides/fedora-post-install/#matlab) or [MATLAB-ArchWiki](https://wiki.archlinux.org/index.php/MATLAB) for instructions.
+If you use MATLAB, we strongly advise to also rename or exclude the GCC libraries shipped with MATLAB to avoid possible conflicts with GCC libraries shipped by Fedora, see e.g. the [OpenGL Section of the ArchWiki article on MATLAB](https://wiki.archlinux.org/title/MATLAB#OpenGL_acceleration) for instructions.
+
Now use the following commands if using MATLAB (adapt them for Octave, see above):
```sh
@@ -291,6 +293,43 @@ bison --version # bison (GNU Bison) 3.6.4
```
Now configure dynare as above.
+## Arch Linux
+
+The following steps show how to install Dynare on Arch Linux from source.
+- Install all needed dependencies:
+```sh
+pacman -S git make meson boost blas gsl libmatio gcc-fortran gcc-libs
+```
+- Compile and install SLICOT:
+```sh
+wget https://github.com/SLICOT/SLICOT-Reference/archive/refs/tags/v5.9.tar.gz
+tar xf v5.9.tar.gz
+cd SLICOT-Reference-5.9
+make -f makefile_Unix FORTRAN=gfortran OPTS="-O2 -fno-underscoring -fdefault-integer-8" LOADER=gfortran lib
+mkdir -p /usr/local/lib
+cp slicot.a /usr/local/lib/libslicot64_pic.a
+cd ..
+```
+- Install `x13as-bin` from the AUR ([link to PKGBUILD](https://aur.archlinux.org/packages/x13as-bin)), either by using your favorite AUR-helper or by following the [Arch wiki](https://wiki.archlinux.org/title/Arch_User_Repository)
+- Prepare the Dynare sources:
+```sh
+git clone --recurse-submodules https://git.dynare.org/Dynare/dynare.git
+cd dynare
+```
+- Configure Dynare from the source directory (adjust `matlab_path` if you use a different version than R2024a):
+```sh
+meson setup -Dmatlab_path=/usr/local/MATLAB/R2024a -Dbuildtype=debugoptimized build-matlab
+```
+- Compile:
+```sh
+meson compile -C build-matlab
+```
+- Optionally, run the testsuite:
+```sh
+meson test -C build-matlab
+```
+If you run into errors with GCC libraries/`libstdc++`, you need to prevent MATLAB from using the libraries it is shipping. For example, you could follow the [OpenGL Section of the ArchWiki article on MATLAB](https://wiki.archlinux.org/title/MATLAB#OpenGL_acceleration) by setting/exporting `LD_PRELOAD` and `LD_LIBRARY` (the exports could also be placed in your `.zprofile`). If necessary, add `java.opts` as described in the section.
+
## Windows
- Install [MSYS2](http://www.msys2.org)
diff --git a/doc/manual/source/bibliography.rst b/doc/manual/source/bibliography.rst
index 90ec693d0b9e1fd54fd56a67c18b82c76e1ad44f..9e17ff619a4dba24360f6f1d9716e3e807fa29ce 100644
--- a/doc/manual/source/bibliography.rst
+++ b/doc/manual/source/bibliography.rst
@@ -14,6 +14,7 @@ Bibliography
* Backus, David K., Patrick J. Kehoe, and Finn E. Kydland (1992): “International Real Business Cycles,” *Journal of Political Economy*, 100(4), 745–775.
* Baxter, Marianne and Robert G. King (1999): “Measuring Business Cycles: Approximate Band-pass Filters for Economic Time Series,” *Review of Economics and Statistics*, 81(4), 575–593.
* Bini, Dario A., Guy Latouche, and Beatrice Meini (2002): “Solving matrix polynomial equations arising in queueing problems,” *Linear Algebra and its Applications*, 340, 225–244.
+* Boehl, Gregor (2022): “DIME MCMC: A Swiss Army Knife for Bayesian Inference”, *SSRN No. 4250395*
* Born, Benjamin and Johannes Pfeifer (2014): “Policy risk and the business cycle”, *Journal of Monetary Economics*, 68, 68-85.
* Boucekkine, Raouf (1995): “An alternative methodology for solving nonlinear forward-looking models,” *Journal of Economic Dynamics and Control*, 19, 711–734.
* Brayton, Flint and Peter Tinsley (1996): “A Guide to FRB/US: A Macroeconomic Model of the United States,” *Finance and Economics Discussion Series*, 1996-42.
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 8c0b4b72bf7c2a3c316d736930456573fa7640f9..68c9871d4909fbe6556b2e78d2ba563b2c2537ee 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -7462,6 +7462,12 @@ observed variables.
density is recentered to the previous draw in every
step.
+ ``'dime_mcmc'``
+
+ Instructs Dynare to use the Differential-Independence Mixture Ensemble ("DIME") MCMC sampler of *Boehl (2022)* instead of the standard Random-Walk Metropolis-Hastings. DIME is robust for odd shaped, multimodal, black-box distributions and shown to require significantly less likelihood evaluations than alternative samplers. Many chains run simultaneously, thereby further increasing sampling speed. The algorithm is based on a gradient-free global multi-start optimizer and does not require any posterior mode density maximization prior to MCMC sampling. DIME proposals are generated from an endogenous and adaptive proposal distribution, thereby providing close-to-optimal proposal distributions for black box target distributions without manual fine-tuning. Does not yet support ``moments_varendo``, ``bayesian_irf``, and ``smoother``.
+
+ Note that, since DIME is using parameter transformations, setting parameter bounds is often counterproductive. The ``prior_trunc`` option is disabled and set to zero.
+
``'tailored_random_block_metropolis_hastings'``
Instructs Dynare to use the Tailored randomized block
@@ -7578,6 +7584,42 @@ observed variables.
the draws when the current ``_mh*_blck`` file is
full. Default: ``0``
+ ``'dime_mcmc'``
+
+ Available options are:
+
+ ``'niter'``
+
+ Number of iterations. Default value is 1500.
+
+ ``'nchain'``
+
+ Number of chains/particles. A range between :math:`4 ndim` and :math:`6 ndim` is recommended. For very difficult problems, double the number of chains. Default value is: :math:`5 ndim`.
+
+ ``'tune'``
+
+ Number of ensemble iterations to keep after burnin. Default value is chosen to obtain at least 50,000 samples.
+
+ ``'aimh_prob'``
+
+ Probability to draw a global transition kernel. By default this is set to :math:`0.1`. It is usually not necessary to change this value.
+
+ ``'df_proposal_dist'``
+
+ Degrees of freedom of the multivariate t distribution used for global kernel proposals. Defaults to :math:`10`.
+
+ ``'rho'``
+
+ Decay parameter for the mean and covariances of the global transistion kernel. Defaults to :math:`0.999`.
+
+ ``'gamma'``
+
+ Mean stretch factor for the proposal vector. By default, it is :math:`2.38 / \sqrt{2\,\mathrm{ndim}}` as recommended by *ter Braak (2006)*
+
+ ``'sigma'``
+
+ Standard deviation of the Gaussian used to stretch the proposal vector. This is normally negligible. Defaults to :math:`1e-5`.
+
``'independent_metropolis_hastings'``
Takes the same options as in the case of ``random_walk_metropolis_hastings``.
diff --git a/license.txt b/license.txt
index 45842121d569505d897dd14c9757fb9a3424cf2b..51b37da97b0b9067955540c6c0ec037bc28f0454 100644
--- a/license.txt
+++ b/license.txt
@@ -118,6 +118,10 @@ Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin
2013-2017 Dynare Team
License: GPL-3+
+Files: matlab/estimation/smc/logsumexp.m
+Copyright: 2020 Nicholas J. Higham
+License: BSD-2-clause
+
Files: matlab/one_sided_hp_filter.m
Copyright: 2010-2015 Alexander Meyer-Gohde
2015-2017 Dynare Team
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 5af50faa2733af6b162b2bdf53e44d596f64a64a..9e6b685ac83f528558a918dc5f387b7d959fa18b 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -498,6 +498,13 @@ options_.posterior_sampler_options.hssmc.particles = 20000;
options_.posterior_sampler_options.hssmc.scale = 0.5;
options_.posterior_sampler_options.hssmc.acpt = 1.00;
options_.posterior_sampler_options.hssmc.target = 0.25;
+% DIME MCMC Sampler
+options_.posterior_sampler_options.dime.niter = 1500;
+options_.posterior_sampler_options.dime.parallel = false;
+options_.posterior_sampler_options.dime.aimh_prob = .1;
+options_.posterior_sampler_options.dime.sigma = 1e-5;
+options_.posterior_sampler_options.dime.df_proposal_dist = 10;
+options_.posterior_sampler_options.dime.rho = 0.999;
% DSMH: Dynamic Striated Metropolis-Hastings algorithm
options_.posterior_sampler_options.dsmh.H = 25 ;
options_.posterior_sampler_options.dsmh.N = 20 ;
diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index dcb6e80dc7eb25919b2388bc5c3581dc2988f915..fd6cdce1cdbf6b93c02d86fed1db3587e46692e2 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -42,6 +42,16 @@ if ishssmc(options_)
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
draws = transpose(posterior.particles(column,:));
+elseif isdime(options_)
+ posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
+ tune = posterior.tune;
+ chains = posterior.chains(end-tune:end,:,:);
+ if column>0
+ chains = reshape(chains, [], size(chains, 3));
+ draws = chains(:,column);
+ else
+ draws = posterior.lprobs;
+ end
else
iline = FirstLine;
linee = 1;
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index cd791de522c7007e206bf6dd883783c9f2a369c9..299ccdaca07cc4fafb3cc257e63b3c69336d3bc3 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -69,6 +69,8 @@ if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
% Set function handle for GetAllPosteriorDraws
getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
+elseif isdime(options_)
+ getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
else
if ~isfield(oo_,'MarginalDensity') || (issmc(options_) && ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean'))
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
@@ -81,6 +83,11 @@ end
if ishssmc(options_)
num_draws = options_.posterior_sampler_options.hssmc.particles;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
+elseif isdime(options_)
+ nchain = options_.posterior_sampler_options.current_options.nchain;
+ tune = options_.posterior_sampler_options.current_options.tune;
+ num_draws = nchain*tune;
+ hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
else
num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
diff --git a/matlab/estimation/check_posterior_sampler_options.m b/matlab/estimation/check_posterior_sampler_options.m
index 29de7cf6891de044ffe52e957f037e37262597fa..e93510075baf1f6956c47beb1a5c148f9eefc3a7 100644
--- a/matlab/estimation/check_posterior_sampler_options.m
+++ b/matlab/estimation/check_posterior_sampler_options.m
@@ -405,6 +405,49 @@ if init
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
+ case 'dime_mcmc'
+
+ % default options
+ posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dime);
+
+ % 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 'nchain'
+ posterior_sampler_options.nchain = options_list{i,2};
+ case 'niter'
+ posterior_sampler_options.niter = options_list{i,2};
+ case 'parallel'
+ posterior_sampler_options.parallel = options_list{i,2};
+ case 'tune'
+ posterior_sampler_options.tune = options_list{i,2};
+ case 'aimh_prob'
+ posterior_sampler_options.aimh_prob = options_list{i,2};
+ case 'sigma'
+ posterior_sampler_options.sigma = options_list{i,2};
+ case 'gamma'
+ posterior_sampler_options.gamma = options_list{i,2};
+ case 'df_proposal_dist'
+ posterior_sampler_options.df_proposal_dist = options_list{i,2};
+ case 'rho'
+ posterior_sampler_options.rho = options_list{i,2};
+ otherwise
+ warning(['dime: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+
+ if ~isfield(posterior_sampler_options,'tune')
+ posterior_sampler_options.tune = ceil(50000 / posterior_sampler_options.nchain);
+ dprintf('dime: setting number of ensemble iterations to keep (`tune`) to %d.', posterior_sampler_options.tune)
+ end
+ options_.mode_compute = 0;
+ options_.cova_compute = 0;
+ options_.mh_replic = 0;
+ options_.mh_posterior_mode_estimation = false;
+
case 'dsmh'
% default options
diff --git a/matlab/estimation/dynare_estimation_1.m b/matlab/estimation/dynare_estimation_1.m
index 444280a79a80e026feb8fd3357e55a92c0800342..52a51a4cf5862903152be09b763afef77b5c07d7 100644
--- a/matlab/estimation/dynare_estimation_1.m
+++ b/matlab/estimation/dynare_estimation_1.m
@@ -31,6 +31,8 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
+dispString = 'Estimation::mcmc';
+
if issmc(options_)
options_.mode_compute = 0;
options_.mh_replic = 0;
@@ -38,8 +40,10 @@ if issmc(options_)
options_.load_mh_file = false;
options_.load_results_after_load_mh = false;
end
-
-dispString = 'Estimation::mcmc';
+if isdime(options_) && options_.prior_trunc
+ options_.prior_trunc = 0;
+ fprintf('%s: DIME requires no prior truncation. Resetting options_.prior_trunc=0.\n', dispString);
+end
if ~exist([M_.dname filesep 'Output'],'dir')
mkdir(M_.dname,'Output');
@@ -396,6 +400,10 @@ if ishssmc(options_)
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = posterior_sampler_options;
oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+elseif isdime(options_)
+ [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
+ options_.posterior_sampler_options.current_options = posterior_sampler_options;
+ dime(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
elseif isdsmh(options_)
dsmh(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
end
@@ -456,9 +464,12 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
oo_.convergence=oo_load_mh.oo_.convergence;
end
end
+ elseif isdime(options_) && ~options_.nodiagnostic
+ % provide plot of log densities over iterations
+ oo_.lprob = trace_plot_dime(options_, M_);
end
% Estimation of the marginal density from the Mh draws:
- if ishssmc(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
+ if ishssmc(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
if ~issmc(options_)
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
@@ -469,8 +480,9 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
end
% Store posterior mean in a vector and posterior variance in
% a matrix
- [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
- = GetPosteriorMeanVariance(options_, M_);
+ if ~isdime(options_)
+ [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] = GetPosteriorMeanVariance(options_, M_);
+ end
elseif options_.load_mh_file && options_.load_results_after_load_mh
% load fields from previous MCMC run stored in results-file
field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
@@ -537,7 +549,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
end
end
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
- if ~ishssmc(options_)
+ if ~ishssmc(options_) && ~isdime(options_)
if error_flag
error('%s: I cannot compute the posterior statistics!',dispString)
end
diff --git a/matlab/estimation/dynare_estimation_init.m b/matlab/estimation/dynare_estimation_init.m
index 8378285dcef48104ed10faa21c451afb73407a1a..13c1c9e8328ec63c1e0fd124b7b2fbd2f421279f 100644
--- a/matlab/estimation/dynare_estimation_init.m
+++ b/matlab/estimation/dynare_estimation_init.m
@@ -372,6 +372,9 @@ end
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
+% set default number of chains for DIME
+options_.posterior_sampler_options.dime.nchain = 5 * length(xparam1);
+
% Build the dataset
if ~isempty(options_.datafile)
[~,name] = fileparts(options_.datafile);
diff --git a/matlab/estimation/generate_trace_plots.m b/matlab/estimation/generate_trace_plots.m
index 409ab22bdfa78eb09635cd34672a809640509446..b6e6c8b62f70c3495822dab56d2e9341f9e1a58a 100644
--- a/matlab/estimation/generate_trace_plots.m
+++ b/matlab/estimation/generate_trace_plots.m
@@ -30,6 +30,9 @@ function generate_trace_plots(chain_number)
global M_ options_ estim_params_
+if issmc(options_)
+ error('generate_trace_plots:: SMC methods do not support trace plots')
+end
% Get informations about the posterior draws:
MetropolisFolder = CheckPath('metropolis', M_.dname);
@@ -65,4 +68,4 @@ for ii=1:size(estim_params_.corrx, 1)
parameter_name_1 = M_.exo_names{estim_params_.corrx(ii,1)};
parameter_name_2 = M_.exo_names{estim_params_.corrx(ii,2)};
trace_plot(options_, M_, estim_params_, 'StructuralShock', chain_number, parameter_name_1, parameter_name_2)
-end
\ No newline at end of file
+end
diff --git a/matlab/estimation/smc/dime.m b/matlab/estimation/smc/dime.m
new file mode 100644
index 0000000000000000000000000000000000000000..bb2e5886e5f4bd639e979dba661b23fc14bf8a23
--- /dev/null
+++ b/matlab/estimation/smc/dime.m
@@ -0,0 +1,245 @@
+function dime(TargetFun, init_x, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, dr , steady_state, exo_steady_state, exo_det_steady_state)
+
+% Differential-Independence Mixture Ensemble ("DIME") MCMC sampling
+% as proposed in "Ensemble MCMC Sampling for Robust Bayesian Inference"
+% (Gregor Boehl, 2022, SSRN No. 4250395):
+%
+% https://gregorboehl.com/live/dime_mcmc_boehl.pdf
+%
+% INPUTS
+% - TargetFun [char] string specifying the name of the objective function (posterior kernel).
+% - init_x [double] p×1 vector of parameters to be estimated (initial values, not used).
+% - mh_bounds [double] p×2 matrix defining lower and upper bounds for the parameters.
+% - dataset_ [dseries] sample
+% - dataset_info [struct] informations about the dataset
+% - options_ [struct] Dynare's options
+% - M_ [struct] model description
+% - estim_params_ [struct] estimated parameters
+% - bayestopt_ [struct] describing the priors
+% - steady_state [double] steady state of endogenous variables
+% - exo_steady_state [double] steady state of exogenous variables
+% - exo_det_steady_state [double] steady state of exogenous deterministic variables
+%
+% SPECIAL REQUIREMENTS
+% None.
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+ % initialize
+ ndim = length(init_x);
+ prop_cov = eye(ndim);
+ prop_mean = zeros(1,ndim);
+ naccepted = 1;
+ cumlweight = -inf;
+ fixPSD = 1e-16*prop_cov;
+
+ % get options & assign problem specific default value
+ opts = options_.posterior_sampler_options.current_options;
+ if isfield(opts,'gamma')
+ g0 = opts.gamma;
+ else
+ g0 = 2.38 / sqrt(2 * ndim);
+ end
+ nchain = opts.nchain;
+ dft = opts.df_proposal_dist;
+ tune = opts.tune;
+ if tune > opts.niter
+ error('dime: number of ensemble iterations to keep (tune=%d) exceeds number of iterations (niter=%d)!', tune, opts.niter)
+ end
+
+ % temporary workaround to deal with parallelization
+ if opts.parallel
+ % check if parallel is possible
+ installed_toolboxes = ver;
+ toolbox_installed = any(strcmp("Parallel Computing Toolbox", {installed_toolboxes.Name}));
+ if ~toolbox_installed || feature('numcores') == 1
+ error('dime: parallel processing option is chosen but Parallel Computing Toolbox is not installed or machine has only one core')
+ end
+ end
+
+ % Set location for the simulated particles.
+ SimulationFolder = CheckPath('dime', M_.dname);
+ % Define prior distribution
+ Prior = dprior(bayestopt_, options_.prior_trunc);
+ bounds = check_bounds(Prior, mh_bounds);
+
+ % Set function handle for the objective
+ eval(sprintf('%s = @(x) %s(x, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, mh_bounds, dr , steady_state, exo_steady_state, exo_det_steady_state, []);', 'funobj', func2str(TargetFun)));
+
+ % Initialization of the sampler (draws from the prior distribution with finite logged likelihood)
+ t0 = tic;
+ [x, ~, lprob] = ...
+ smc_samplers_initialization(funobj, 'dime', nchain, Prior, SimulationFolder, opts.niter);
+ tt = toc(t0);
+ x = ptransform(x', bounds, true);
+
+ disp_verbose(sprintf('Estimation:dime: log-posterior standard deviation (post.std) of a multivariate normal would be %.2f.\n', sqrt(0.5*ndim)), options_.verbosity);
+
+ % preallocate
+ chains = zeros(opts.niter, nchain, ndim);
+ lprobs = zeros(opts.niter, nchain);
+ isplit = fix(nchain/2) + 1;
+
+ for iter = 1:opts.niter
+
+ for complementary_ensemble = [false true]
+
+ % define current ensemble
+ if complementary_ensemble
+ idcur = 1:isplit;
+ idref = (isplit+1):nchain;
+ else
+ idcur = (isplit+1):nchain;
+ idref = 1:isplit;
+ end
+ cursize = length(idcur);
+ refsize = nchain - cursize;
+
+ % get differential evolution proposal
+ % draw the indices of the complementary chains
+ i1 = (0:cursize-1) + randi([1 cursize-1], 1, cursize);
+ i2 = (0:cursize-1) + randi([1 cursize-2], 1, cursize);
+ i2(i2 >= i1) = i2(i2 >= i1) + 1;
+
+ % add small noise and calculate proposal
+ f = opts.sigma*randn(cursize, 1);
+ q = x(idcur,:) + g0 .* (x(idref(mod(i1, refsize) + 1),:) - x(idref(mod(i2, refsize) + 1),:)) + f;
+ factors = zeros(cursize,1);
+
+ % get AIMH proposal
+ xchnge = rand(cursize,1) <= opts.aimh_prob;
+
+ % draw alternative candidates and calculate their proposal density
+ prop_cov_chol = chol(prop_cov*(dft - 2)/dft + fixPSD);
+ xcand = rand_multivariate_student(prop_mean, prop_cov_chol, dft, sum(xchnge));
+ lprop_old = log(multivariate_student_pdf(x(idcur(xchnge),:), prop_mean, prop_cov_chol, dft));
+ lprop_new = log(multivariate_student_pdf(xcand, prop_mean, prop_cov_chol, dft));
+
+ % update proposals and factors
+ q(xchnge,:) = xcand;
+ factors(xchnge) = lprop_old - lprop_new;
+
+ % Metropolis-Hasings
+ newlprob = log_prob_fun(funobj, Prior, bounds, opts.parallel, q);
+ lnpdiff = factors + newlprob - lprob(idcur);
+ accepted = lnpdiff > log(rand(cursize,1));
+ naccepted = naccepted + sum(accepted);
+
+ % update chains
+ x(idcur(accepted),:) = q(accepted,:);
+ lprob(idcur(accepted)) = newlprob(accepted);
+ end
+
+ % store
+ chains(iter,:,:) = ptransform(x, bounds, false);
+ lprobs(iter,:) = lprob;
+
+ % log weight of current ensemble
+ lweight = logsumexp(lprob) + log(naccepted) - log(nchain);
+
+ % calculate stats for current ensemble
+ ncov = cov(x);
+ nmean = mean(x);
+
+ % update AIMH proposal distribution
+ newcumlweight = logsumexp([cumlweight lweight]);
+ prop_cov = exp(cumlweight - newcumlweight) .* prop_cov + exp(lweight - newcumlweight) .* ncov;
+ prop_mean = exp(cumlweight - newcumlweight) .* prop_mean + exp(lweight - newcumlweight) .* nmean;
+ cumlweight = newcumlweight + log(opts.rho);
+
+ % be informative
+ tt = toc(t0);
+ if iter == 1 | ~mod(iter,15)
+ disp_verbose(' #iter. post.mode post.std %accept lapsed', options_.verbosity)
+ disp_verbose(' ------ --------- -------- ------- ------', options_.verbosity)
+ end
+ counter = sprintf('%u/%u', iter, opts.niter);
+ disp_verbose(sprintf('%11s %10.3f %1.2e %5.2f%% %5.2fs', counter, max(lprob), std(lprob), 100*naccepted/nchain, tt), options_.verbosity)
+ naccepted = 0;
+ end
+
+ save(sprintf('%s%schains.mat', SimulationFolder, filesep()), 'chains', 'lprobs', 'tune')
+end
+
+function lprobs = log_prob_fun(loglikefun, Prior, bounds, runs_in_parallel, x)
+
+ dim = length(x);
+ x = ptransform(x, bounds, false);
+ lprobs = zeros(dim,1);
+ if runs_in_parallel
+ parfor i=1:dim
+ par = x(i,:)';
+ loglikelihood = -loglikefun(par);
+ % temporary workaround for bug #1930
+ logprior = Prior.density(par);
+ lprobs(i) = loglikelihood + logprior;
+ end
+ else
+ for i=1:dim
+ lprobs(i) = -loglikefun(x(i,:)');
+ end
+ end
+end
+
+function nbounds = check_bounds(Prior, bounds)
+
+ % ensure that prior is continous at bounds
+ nbounds = bounds;
+ prior_at_mean = Prior.density(Prior.mean);
+ dim = length(Prior.mean);
+ for j = 1:dim
+ if ~isinf(bounds.ub(j))
+ x = Prior.mean;
+ x(j) = bounds.ub(j);
+ if Prior.density(x - 1e-16) - prior_at_mean > log(1e-16)
+ nbounds.ub(j) = nbounds.ub(j) + 1e-2;
+ end
+ x = Prior.mean;
+ x(j) = bounds.lb(j);
+ if Prior.density(x + 1e-16) - prior_at_mean > log(1e-16)
+ nbounds.lb(j) = nbounds.lb(j) - 1e-2;
+ end
+ end
+ end
+end
+
+function x = ptransform(x, bounds, con2unc)
+
+ % parameter transformation
+ lb = bounds.lb;
+ ub = bounds.ub;
+ dim = size(x,2);
+ for j = 1:dim
+ % one-sided
+ if ~isinf(lb(j)) && isinf(ub(j))
+ if con2unc
+ x(:,j) = log(x(:,j) - lb(j));
+ else
+ x(:,j) = lb(j) + exp(x(:,j));
+ end
+ % two-sided
+ elseif ~isinf(lb(j)) && ~isinf(ub(j))
+ if con2unc
+ y = (x(:,j) - lb(j))/(ub(j) - lb(j));
+ x(:,j) = log(y ./ (1 - y));
+ else
+ x(:,j) = lb(j) + (ub(j) - lb(j)) ./ (1 + exp(-x(:,j)));
+ end
+ end
+ end
+end
diff --git a/matlab/estimation/smc/isdime.m b/matlab/estimation/smc/isdime.m
new file mode 100644
index 0000000000000000000000000000000000000000..11f6a4ba9b92dbdf41dec34b9562ff9e1e40acc1
--- /dev/null
+++ b/matlab/estimation/smc/isdime.m
@@ -0,0 +1,25 @@
+function bool = isdime(options_)
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+bool = false;
+if isfield(options_, 'posterior_sampler_options')
+ if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'dime_mcmc')
+ bool = true;
+ end
+end
diff --git a/matlab/estimation/smc/logsumexp.m b/matlab/estimation/smc/logsumexp.m
new file mode 100644
index 0000000000000000000000000000000000000000..cad3165568cdaee5daf6eae798e432b60f9cae85
--- /dev/null
+++ b/matlab/estimation/smc/logsumexp.m
@@ -0,0 +1,56 @@
+function [lse,sm] = logsumexp(x)
+%LOGSUMEXP Log-sum-exp function.
+% lse = LOGSUMEXP(x) returns the log-sum-exp function evaluated at
+% the vector x, defined by lse = log(sum(exp(x)).
+% [lse,sm] = LOGSUMEXP(x) also returns the softmax function evaluated
+% at x, defined by sm = exp(x)/sum(exp(x)).
+% The functions are computed in a way that avoids overflow and
+% optimizes numerical stability.
+%
+% Reference:
+% P. Blanchard, D. J. Higham, and N. J. Higham.
+% Accurately computing the log-sum-exp and softmax functions.
+% IMA J. Numer. Anal., Advance access, 2020.
+%
+% Taken from https://de.mathworks.com/matlabcentral/fileexchange/84892-logsumexp-softmax
+%
+% Copyright (c) 2020, Nicholas J. Higham
+% All rights reserved.
+%
+% Redistribution and use in source and binary forms, with or without
+% modification, are permitted provided that the following conditions are met:
+%
+% * Redistributions of source code must retain the above copyright notice, this
+% list of conditions and the following disclaimer.
+%
+% * Redistributions in binary form must reproduce the above copyright notice,
+% this list of conditions and the following disclaimer in the documentation
+% and/or other materials provided with the distribution.
+%
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+if ~isvector(x), error('Input x must be a vector.'), end
+
+n = length(x);
+s = 0; e = zeros(n,1);
+[xmax,k] = max(x); a = xmax;
+s = 0;
+for i = 1:n
+ e(i) = exp(x(i)-xmax);
+ if i ~= k
+ s = s + e(i);
+ end
+end
+lse = a + log1p(s);
+if nargout > 1
+ sm = e/(1+s);
+end
diff --git a/matlab/estimation/smc/trace_plot_dime.m b/matlab/estimation/smc/trace_plot_dime.m
new file mode 100644
index 0000000000000000000000000000000000000000..02719db9903a88c8de026f82b799f881500e6698
--- /dev/null
+++ b/matlab/estimation/smc/trace_plot_dime.m
@@ -0,0 +1,70 @@
+function lprobs_sample = trace_plot_dime(options_, M_)
+
+% Plot the history of the densities of an ensemble to visually inspect convergence.
+%
+% INPUTS
+% - options_ [struct] dynare's options
+% - M_ [struct] model description
+% - oo_ [struct] outputs
+%
+% SPECIAL REQUIREMENTS
+% None.
+
+% Copyright © 2022-2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+
+lprobs = GetAllPosteriorDraws(options_, M_.dname, [], 0);
+tune = options_.posterior_sampler_options.current_options.tune;
+[niter, nchain] = size(lprobs);
+
+if max(lprobs(end,:)) > 0
+ ulim = max(lprobs(end,:))*1.2;
+else
+ ulim = max(lprobs(end,:))/1.2;
+end
+if min(lprobs(end,:)) > 0
+ llim = min(lprobs(end,:))/5;
+else
+ llim = min(lprobs(end,:))*5;
+end
+
+graphFolder = CheckPath('graphs',M_.dname);
+latexFolder = CheckPath('latex',M_.dname);
+hh_fig = dyn_figure(options_.nodisplay,'Name','DIME Convergence Diagnostics');
+
+hold on
+lines_tune = plot(1:niter-tune, lprobs(1:end-tune,:), 'Color', '#D95319');
+lines_sample = plot(niter-tune:niter, lprobs(end-tune:end,:), 'Color', '#0072BD');
+for i = 1:nchain
+ lines_tune(i).Color(4) = min(1,10/nchain);
+ lines_sample(i).Color(4) = min(1,10/nchain);
+end
+ylim([llim, ulim])
+hold off
+
+dyn_saveas(hh_fig,[graphFolder '/' M_.fname '_trace_lprob'],options_.nodisplay,options_.graph_format);
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+ fprintf(fidTeX,'\\begin{figure}[H]\n');
+ fprintf(fidTeX,'\\centering \n');
+ fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s_trace_lprob}\n',[graphFolder '/' M_.fname]);
+ fprintf(fidTeX,'\\caption{Ensemble traces of posterior densities for DIME.\n');
+ fprintf(fidTeX,'\\label{Fig:DIME_trace}\n');
+ fprintf(fidTeX,'\\end{figure}\n');
+ fprintf(fidTeX,'\n');
+end
+
+lprobs_sample = lprobs(end-tune:end,:);
diff --git a/matlab/issmc.m b/matlab/issmc.m
index a79d444db9239152acc6d962a0b7345b8309138d..0f80730a55635be528e9fc6cf648921ff57b5957 100644
--- a/matlab/issmc.m
+++ b/matlab/issmc.m
@@ -17,4 +17,4 @@ function bool = issmc(options_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-bool = ishssmc(options_) || isdsmh(options_);
+bool = ishssmc(options_) || isdsmh(options_) || isdime(options_);
diff --git a/meson.build b/meson.build
index 8bd67e600f125b63617bacfe8c81ed5415ad4ced..e2466bca8b58519a22f87ea9d72553a05c046283 100644
--- a/meson.build
+++ b/meson.build
@@ -845,6 +845,8 @@ mod_and_m_tests = [
'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'estimation/fs2000.mod' ],
'extra' : [ 'estimation/fsdat_simul.m' ] },
+ { 'test' : [ 'estimation/dime_mcmc/fs2000.mod' ],
+ 'extra' : [ 'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'estimation/hssmc/fs2000.mod' ],
'extra' : [ 'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'gsa/ls2003a.mod',
diff --git a/tests/estimation/dime_mcmc/fs2000.mod b/tests/estimation/dime_mcmc/fs2000.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3491ce5c5a832eda88c090faa3ec5449f265a9d5
--- /dev/null
+++ b/tests/estimation/dime_mcmc/fs2000.mod
@@ -0,0 +1,89 @@
+// 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;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+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;
+W = l/n;
+-(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;
+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);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+options_.solve_tolf = 1e-12;
+
+estimation(order=1, datafile='../fsdat_simul.m', nobs=192, loglinear,
+ posterior_sampling_method='dime_mcmc',
+ posterior_sampler_options=('tune',200, 'niter',400),
+bayesian_irf, smoother, moments_varendo,consider_all_endogenous
+);