Compare commits

...

33 Commits

Author SHA1 Message Date
Stéphane Adjemian (Ryûk) a32db79a75
[WIP] Add Sequential Monte Carlo samplers. 2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 2fbbe66c0a
Add member to dprior class.
Name of the parameter.
2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 61498e644a
One file per method. 2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 3606b10f05
Add methods for computing moments.
- prior mean
 - prior mode
 - prior median
 - prior variance
2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 5077969aad
Add members to @dprior class. 2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 3d50844ae4
Make last input argument optional. 2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 3c3353b7ed
Add methods to dprior (density and densities).
Will be used as a replacement for priordens.
2023-12-12 18:18:38 +01:00
Stéphane Adjemian (Ryûk) 03a68ddb89
Cosmetic changes. 2023-12-12 18:18:38 +01:00
Sébastien Villemot b1aa88e8da
Add clang-tidy configuration file
[skip ci]
2023-12-12 17:32:43 +01:00
Sébastien Villemot d3aac5e2d7
Fix typo
[skip ci]
2023-12-12 17:28:52 +01:00
Sébastien Villemot 62b31aa279
Merge branch 'cosmetics' of git.dynare.org:JohannesPfeifer/dynare
Ref. !2218
2023-12-12 17:15:38 +01:00
Sébastien Villemot 9d6a25e368
dseries: update to CI configuration 2023-12-12 17:10:42 +01:00
Sébastien Villemot 43b24facb9
Configuration file: new default location; new default value for GlobalInitFile
Under Linux and macOS, the default location for the configuration file is now
dynare/dynare.ini under the configuration directories as defined by the XDG
specification. Under Windows, the default configuration file is now
%APPDATA%\dynare\dynare.ini.

There is now a default value for the global initialization file (GlobalInitFile
option of the configuration file): the global_init.m in the Dynare
configuration directory.
2023-12-12 17:09:55 +01:00
Sébastien Villemot cc15281b1f
Emacs mode: add missing keywords
– “bvar_irf” added to statements
– “target”, “auxname_target_nonstationary”, “component”, “growth”, “auxname”,
  “kind” added to statements-like
– “filter_initial_state” added to blocks
2023-12-12 10:54:03 +01:00
Sébastien Villemot c99230825f
Windows package: bump dependencies 2023-12-12 10:54:03 +01:00
Sébastien Villemot b7805cc667 Merge branch 'remove_files' into 'master'
Remove various unused files

See merge request Dynare/dynare!2217
2023-12-11 20:38:06 +00:00
Johannes Pfeifer ec76bda254 Remove obsolete Sylvester options
dr_block has been removed
2023-12-11 18:04:43 +01:00
Johannes Pfeifer 021b9dbb25 identication.checks.m: remove wrong condition 2023-12-11 18:04:43 +01:00
Johannes Pfeifer daecd1f720 DsgeSmoother.m: remove unnecessary space 2023-12-11 18:04:42 +01:00
Johannes Pfeifer 5a3d545db2 var_sample_moments.m: cosmetic changes 2023-12-11 18:04:42 +01:00
Johannes Pfeifer ed80c4ff3f load_last_mh_history_file.m: cosmetic changes 2023-12-11 18:04:42 +01:00
Johannes Pfeifer 678bd7aca9 dyn_forecast.m: cosmetic header fix 2023-12-11 18:04:42 +01:00
Johannes Pfeifer 97f6a4219b smirnov_test.m: update call to histc under Matlab 2023-12-11 18:04:41 +01:00
Johannes Pfeifer 31c91080e1 Remove shiftS.m, which is a duplicate of the one in dseries 2023-12-11 18:01:34 +01:00
Johannes Pfeifer 62e8b275a0 Remove further unused function from matlab folder 2023-12-11 18:01:33 +01:00
Johannes Pfeifer 435b103cf5 Remove unused functions, mostly related to old analytical derivatives 2023-12-11 18:01:33 +01:00
Sébastien Villemot d844043877
Merge branch 'plot_shock_decomposition' of git.dynare.org:JohannesPfeifer/dynare
Ref. !2215
2023-12-08 15:49:19 +01:00
Sébastien Villemot a31c76403d
Windows and macOS packages: move meson native/cross files to OS-specific directory 2023-12-08 14:34:14 +01:00
Sébastien Villemot 4ef9245a95
MEX files: remove calls to virtual method during construction
Such calls may bypass virtual dispatch.

Automatically detected by clang-tidy with
clang-analyzer-optin.cplusplus.VirtualCall check.
2023-12-07 18:34:38 +01:00
Sébastien Villemot 91c677ca7f
MEX files: drop C++ preprocessor directives now obsolete
Dynare++ is no longer distributed as a standalone binary.
2023-12-07 17:57:01 +01:00
Sébastien Villemot 56289c72d0
Drop obsolete Dynare++ example
[skip ci]
2023-12-07 16:04:33 +01:00
Sébastien Villemot 1f5f668313
Move docker folder below the scripts folder
[skip ci]
2023-12-07 16:02:33 +01:00
Johannes Pfeifer 54c4e9df09 plot_shock_decomposition.m: filter out case where data to plot is empty and writing the Excel file will crash 2023-12-07 13:50:18 +01:00
99 changed files with 2942 additions and 2248 deletions

9
.clang-tidy Normal file
View File

@ -0,0 +1,9 @@
# NB: to use clang-tidy on the MEX source code, make sure that you have
# libomp-dev installed (the LLVM implementation of OpenMP)
# TODO: add the following check families:
# - performance-*
# - bugprone-*
# - cppcoreguidelines-
Checks: 'modernize-*,-modernize-use-trailing-return-type,-clang-diagnostic-unqualified-std-cast-call'

View File

@ -479,7 +479,7 @@ If you want a certain version (e.g. 5.x) , then add `--single-branch --branch 5.
```sh
export BUILDDIR=build-matlab
export MATLABPATH=/Applications/MATLAB_R2023b.app
arch -$ARCH meson setup --native-file scripts/homebrew-native-$ARCH.ini -Dmatlab_path=$MATLABPATH -Dbuildtype=debugoptimized -Dfortran_args="['-B','$DYNAREDIR/slicot/lib']" $BUILDDIR
arch -$ARCH meson setup --native-file macOS/homebrew-native-$ARCH.ini -Dmatlab_path=$MATLABPATH -Dbuildtype=debugoptimized -Dfortran_args="['-B','$DYNAREDIR/slicot/lib']" $BUILDDIR
```
where you need to adapt the path to MATLAB.
Similarly, if you want to compile for Octave, replace the `-Dmatlab_path` option by `-Dbuild_for=octave`, and change the build directory to `build-octave`.
@ -514,4 +514,4 @@ e.g. by adding this to your mod file. Alternatively, you can create a `startup.m
## Docker
We offer a variety of pre-configured Docker containers for Dynare, pre-configured with Octave and MATLAB including all recommended toolboxes.
These are readily available for your convenience on [Docker Hub](https://hub.docker.com/r/dynare/dynare).
The docker folder contains [information and instructions](docker/README.md) to interact, built and customize the containers.
The `scripts/docker` folder contains [information and instructions](scripts/docker/README.md) to interact, built and customize the containers.

View File

@ -15,11 +15,16 @@ related to the model (and hence not placed in the model file). At the
moment, it is only used when using Dynare to run parallel
computations.
On Linux and macOS, the default location of the configuration file is
``$HOME/.dynare``, while on Windows it is ``%APPDATA%\dynare.ini``
(typically ``c:\Users\USERNAME\AppData\dynare.ini``). You
can specify a non standard location using the ``conffile`` option of
the ``dynare`` command (see :ref:`dyn-invoc`).
On Linux and macOS, the configuration file is searched by default under
``dynare/dynare.ini`` in the configuration directories defined by the XDG
specification (typically ``$HOME/.config/dynare/dynare.ini`` for the
user-specific configuration and ``/etc/xdg/dynare/dynare.ini`` for the
system-wide configuration, the former having precedence over the latter). Under
Windows, the configuration file is searched by default in
``%APPDATA%\dynare\dynare.ini`` (typically
``c:\Users\USERNAME\AppData\Roaming\dynare\dynare.ini``). You can specify a non
standard location using the ``conffile`` option of the ``dynare`` command (see
:ref:`dyn-invoc`).
The parsing of the configuration file is case-sensitive and it should
take the following form, with each option/choice pair placed on a
@ -76,8 +81,15 @@ processing. Currently, there is only one option available.
.. option:: GlobalInitFile = PATH_AND_FILE
The location of the global initialization file to be run at
the end of ``global_initialization.m``.
The location of a global initialization file that can be used to
customize some Dynare internals (typically default option values). This
is a MATLAB/Octave script.
If this option is not specified, Dynare will look for a
``global_init.m`` file in its configuration directory (typically
``$HOME/.config/dynare/global_init.m`` under Linux and macOS, and
``c:\Users\USERNAME\AppData\Roaming\dynare\global_init.m`` under
Windows).
*Example*

View File

@ -4787,31 +4787,6 @@ Computing the stochastic solution
and this option has no effect. More references can be found
`here <https://archives.dynare.org/DynareWiki/PartialInformation>`__ .
.. option:: sylvester = OPTION
Determines the algorithm used to solve the Sylvester equation
for block decomposed model. Possible values for OPTION are:
``default``
Uses the default solver for Sylvester equations
(``gensylv``) based on Ondra Kameniks algorithm (see
`here
<https://www.dynare.org/assets/team-presentations/sylvester.pdf>`__
for more information).
``fixed_point``
Uses a fixed point algorithm to solve the Sylvester
equation (``gensylv_fp``). This method is faster than
the default one for large scale models.
|br| Default value is ``default``.
.. option:: sylvester_fixed_point_tol = DOUBLE
The convergence criterion used in the fixed point
Sylvester solver. Its default value is ``1e-12``.
.. option:: dr = OPTION
@ -7479,6 +7454,18 @@ observed variables.
Chain draws than the MH-algorithm. Its relative (in)efficiency can be investigated via
the reported inefficiency factors.
``'hssmc'``
Instructs Dynare to use the *Herbst and Schorfheide (2014)*
version of the Sequential Monte-Carlo sampler instead of the
standard Random-Walk Metropolis-Hastings.
``'dsmh'``
Instructs Dynare to use the Dynamic Striated Metropolis Hastings
sampler proposed by *Waggoner, Wu and Zha (2016)* instead of the
standard Random-Walk Metropolis-Hastings.
.. option:: posterior_sampler_options = (NAME, VALUE, ...)
A list of NAME and VALUE pairs. Can be used to set options for
@ -7923,15 +7910,6 @@ observed variables.
See :opt:`aim_solver`.
.. option:: sylvester = OPTION
See :opt:`sylvester <sylvester = OPTION>`.
.. option:: sylvester_fixed_point_tol = DOUBLE
See :opt:`sylvester_fixed_point_tol <sylvester_fixed_point_tol
= DOUBLE>` .
.. option:: lyapunov = OPTION
Determines the algorithm used to solve the Lyapunov equation to
@ -9475,16 +9453,6 @@ adding prior information comes at the cost of a loss in efficiency of the estima
See :opt:`lyapunov_doubling_tol <lyapunov_doubling_tol = DOUBLE>`.
Default: ``1e-16``.
.. option:: sylvester = OPTION
See :opt:`sylvester <sylvester = OPTION>`.
Default: ``default``, i.e. uses ``gensylv``.
.. option:: sylvester_fixed_point_tol = DOUBLE
See :opt:`sylvester_fixed_point_tol <sylvester_fixed_point_tol = DOUBLE>`.
Default: ``1e-12``.
.. option:: qz_criterium = DOUBLE
See :opt:`qz_criterium <qz_criterium = DOUBLE>`.

View File

@ -1,6 +0,0 @@
Assuming that the dynare++ binary is in your PATH, you can run the example by using the following command
in a Command Prompt Window:
... > dynare++ example1.mod
Please, read the manual (doc\dynare++\dynare++-tutorial.pdf) for a description of the generated output.

View File

@ -1,45 +0,0 @@
/*
* This Dynare++ mod-file implements the RBC model with time-to-build
* described in Kamenik (2011): "DSGE Models with Dynare++. A Tutorial."
* Note that Dynare++ uses the same stock-at-the-end-of-period timing convention
* as the regular Dynare
*/
var Y, C, K, A, H, B;
varexo EPS, NU;
parameters beta, rho, alpha, delta, theta, psi, tau;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 1/(1.03^0.25);
delta = 0.025;
psi = 0;
theta = 2.95;
model;
C*theta*H^(1+psi) = (1-alpha)*Y;
beta*exp(B)*C/exp(B(1))/C(1)*
(exp(B(1))*alpha*Y(1)/K(1)+1-delta) = 1;
Y = exp(A)*K^alpha*H^(1-alpha);
K = exp(B(-1))*(Y(-1)-C(-1)) + (1-delta)*K(-1);
A = rho*A(-1) + tau*B(-1) + EPS;
B = tau*A(-1) + rho*B(-1) + NU;
end;
initval;
A = 0;
B = 0;
H = ((1-alpha)/(theta*(1-(delta*alpha)/(1/beta-1+delta))))^(1/(1+psi));
Y = (alpha/(1/beta-1+delta))^(alpha/(1-alpha))*H;
K = alpha/(1/beta-1+delta)*Y;
C = Y - delta*K;
end;
vcov = [0.0002 0.00005;
0.00005 0.0001
];
order = 7;

View File

@ -54,7 +54,7 @@ LIB64="$ROOTDIR"/macOS/deps/"$PKG_ARCH"/lib64
## - the macOS linker is different from GNU ld and does not have the equivalent of -Bstatic/-Bdynamic
## - libgfortran.spec does not include --as-needed on macOS, hence it will link the library anyways
## Also, it does not seem possible to override libgfortran.spec with the --specs option.
GCC_VERSION=$(sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" "$ROOTDIR"/scripts/homebrew-native-"$PKG_ARCH".ini)
GCC_VERSION=$(sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" "$ROOTDIR"/macOS/homebrew-native-"$PKG_ARCH".ini)
QUADMATH_DIR=$(mktemp -d)
ln -s "$BREWDIR"/opt/gcc/lib/gcc/"$GCC_VERSION"/libquadmath.a "$QUADMATH_DIR"
@ -67,7 +67,7 @@ cd "$ROOTDIR"
# NB: the addition of -Wl,-ld_classic is a workaround for https://github.com/mesonbuild/meson/issues/12282 (see also the native file)
common_meson_opts=(-Dbuild_for=matlab -Dbuildtype=release -Dprefer_static=true -Dfortran_args="[ '-B', '$LIB64/Slicot/' ]" \
-Dc_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dcpp_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dfortran_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" \
--native-file scripts/homebrew-native-$PKG_ARCH.ini)
--native-file macOS/homebrew-native-$PKG_ARCH.ini)
# Build for MATLAB ⩾ R2018b (x86_64) and MATLAB ⩾ R2023b (arm64)
arch -"$PKG_ARCH" meson setup "${common_meson_opts[@]}" -Dmatlab_path="$MATLAB_PATH" build-matlab --wipe

View File

@ -22,7 +22,7 @@ DEPS_ARCH ?= x86_64 # use x86_64 by default
BREWDIR := $(if $(filter arm64,$(DEPS_ARCH)),/opt/homebrew,/usr/local)
GCC_VERSION = $(shell sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" ../../scripts/homebrew-native-$(DEPS_ARCH).ini)
GCC_VERSION = $(shell sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" ../homebrew-native-$(DEPS_ARCH).ini)
ROOT_PATH = $(realpath .)

View File

@ -160,8 +160,6 @@ options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false);
options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15); % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10); % convergence criterion used in the fixed point Lyapunov solver
options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm
options_mom_ = set_default_option(options_mom_,'sylvester_fp',false); % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12); % convergence criterion used in the fixed point Sylvester solver
% mode check plot
options_mom_.mode_check.nolik = false; % we don't do likelihood (also this initializes mode_check substructure)

146
matlab/@dprior/admissible.m Normal file
View File

@ -0,0 +1,146 @@
function b = admissible(o, d)
% Return true iff d is an admissible draw in a distribution characterized by o.
%
% INPUTS
% - o [dprior] Distribution specification for a n×1 vector of independent continuous random variables
% - d [double] n×1 vector.
%
% OUTPUTS
% - b [logical] scalar.
%
% REMARKS
% None.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> d = Prior.draw()
% >> Prior.admissible(d)
% ans =
%
% logical
%
% 1
% Copyright © 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/>.
b = false;
if ~isequal(length(d), length(o.lb))
return
end
if all(d>=o.lb & d<=o.ub)
b = true;
end
return % --*-- Unit tests --*--
%@test:1
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
ndraws = 10;
% Call the tested routine
try
% Instantiate dprior object
o = dprior(BayesInfo, prior_trunc, false);
% Do simulations in a loop and estimate recursively the mean and the variance.
for i = 1:ndraws
draw = o.draw();
if ~o.admissible(draw)
error()
end
end
t(1) = true;
catch
t(1) = false;
end
T = all(t);
%@eof:1

View File

@ -1,7 +1,15 @@
function [tlogpostkern,loglik] = tempered_likelihood(TargetFun,xparam1,lambda,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_)
% function [tlogpostkern,loglik] = tempered_likelihood(TargetFun,xparam1,lambda,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_)
function lpd = densities(o, X)
% Copyright © 2022 Dynare Team
% Evaluate the logged prior densities at X (for each column).
%
% INPUTS
% - o [dprior]
% - X [double] m×n matrix, n points where the prior density is evaluated.
%
% OUTPUTS
% - lpd [double] 1×n, values of the logged prior density at X.
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
@ -18,7 +26,10 @@ function [tlogpostkern,loglik] = tempered_likelihood(TargetFun,xparam1,lambda,da
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
logpostkern = -feval(TargetFun,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
logprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
loglik = logpostkern-logprior ;
tlogpostkern = lambda*loglik + logprior;
n = columns(X);
lpd = NaN(1, n);
parfor i=1:n
lpd(i) = density(o, X(:,i));
end

384
matlab/@dprior/density.m Normal file
View File

@ -0,0 +1,384 @@
function [lpd, dlpd, d2lpd, info] = density(o, x)
% Evaluate the logged prior density at x.
%
% INPUTS
% - o [dprior]
% - x [double] m×1 vector, point where the prior density is evaluated.
%
% OUTPUTS
% - lpd [double] scalar, value of the logged prior density at x.
% - dlpd [double] m×1 vector, first order derivatives.
% - d2lpd [double] m×1 vector, second order derivatives.
%
% REMARKS
% Second order derivatives holder, d2lpd, has the same rank and shape than dlpd because the priors are
% independent (we would have to use a matrix if non orthogonal priors were allowed in Dynare).
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> lpd = Prior.dsensity(x)
% Copyright © 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/>.
lpd = 0.0;
if nargout>1
dlpd = zeros(1, length(x));
if nargout>2
d2lpd = dlpd;
if nargout>3
info = [];
end
end
end
if o.isuniform
if any(x(o.iduniform)-o.p3(o.iduniform)<0) || any(x(o.iduniform)-o.p4(o.iduniform)>0)
lpd = -Inf ;
if nargout==4
info = o.iduniform((x(o.iduniform)-o.p3(o.iduniform)<0) || (x(o.iduniform)-o.p4(o.iduniform)>0));
end
return
end
lpd = lpd - sum(log(o.p4(o.iduniform)-o.p3(o.iduniform))) ;
if nargout>1
dlpd(o.iduniform) = zeros(length(o.iduniform), 1);
if nargout>2
d2lpd(o.iduniform) = zeros(length(o.iduniform), 1);
end
end
end
if o.isgaussian
switch nargout
case 1
lpd = lpd + sum(lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian)));
case 2
[tmp, dlpd(o.idgaussian)] = lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian));
lpd = lpd + sum(tmp);
case {3,4}
[tmp, dlpd(o.idgaussian), d2lpd(o.idgaussian)] = lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian));
lpd = lpd + sum(tmp);
end
end
if o.isgamma
switch nargout
case 1
lpd = lpd + sum(lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma)));
if isinf(lpd), return, end
case 2
[tmp, dlpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 3
[tmp, dlpd(o.idgamma), d2lpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 4
[tmp, dlpd(o.idgamma), d2lpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
lpd = lpd + sum(tmp);
if isinf(lpd)
info = o.idgamma(isinf(tmp));
return
end
end
end
if o.isbeta
switch nargout
case 1
lpd = lpd + sum(lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta)));
if isinf(lpd), return, end
case 2
[tmp, dlpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 3
[tmp, dlpd(o.idbeta), d2lpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 4
[tmp, dlpd(o.idbeta), d2lpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
lpd = lpd + sum(tmp);
if isinf(lpd)
info = o.idbeta(isinf(tmp));
return
end
end
end
if o.isinvgamma1
switch nargout
case 1
lpd = lpd + sum(lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)));
if isinf(lpd), return, end
case 2
[tmp, dlpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 3
[tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 4
[tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
lpd = lpd + sum(tmp);
if isinf(lpd)
info = o.idinvgamma1(isinf(tmp));
return
end
end
end
if o.isinvgamma2
switch nargout
case 1
lpd = lpd + sum(lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2)));
if isinf(lpd), return, end
case 2
[tmp, dlpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 3
[tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 4
[tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
lpd = lpd + sum(tmp);
if isinf(lpd)
info = o.idinvgamma2(isinf(tmp));
return
end
end
end
if o.isweibull
switch nargout
case 1
lpd = lpd + sum(lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull)));
if isinf(lpd), return, end
case 2
[tmp, dlpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 3
[tmp, dlpd(o.idweibull), d2lpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
lpd = lpd + sum(tmp);
if isinf(lpd), return, end
case 4
[tmp, dlpd(o.idweibull), d2lpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
lpd = lpd + sum(tmp);
if isinf(lpd)
info = o.idweibull(isinf(tmp));
return
end
end
end
return % --*-- Unit tests --*--
%@test:1
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
% Call the tested routine
try
Prior = dprior(BayesInfo, prior_trunc, false);
% Compute density at the prior mode
lpdstar = Prior.density(p5);
% Draw random deviates in a loop and evaluate the density.
LPD = NaN(10000,1);
parfor i = 1:10000
x = Prior.draw();
LPD(i) = Prior.density(x);
end
t(1) = true;
catch
t(1) = false;
end
if t(1)
t(2) = all(LPD<=lpdstar);
end
T = all(t);
%@eof:1
%@test:2
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
% Call the tested routine
try
Prior = dprior(BayesInfo, prior_trunc, false);
mu = NaN(14,1);
std = NaN(14,1);
for i=1:14
% Define conditional density (it's also a marginal since priors are orthogonal)
f = @(x) exp(Prior.densities(substitute(p5, i, x)));
% TODO: Check the version of Octave we use (integral is available as a compatibility wrapper in latest Octave version)
m = integral(f, p3(i), p4(i));
density = @(x) f(x)/m; % rescaling is required since the probability mass depends on the conditioning.
% Compute the conditional expectation
mu(i) = integral(@(x) x.*density(x), p3(i), p4(i));
std(i) = sqrt(integral(@(x) ((x-mu(i)).^2).*density(x), p3(i), p4(i)));
end
t(1) = true;
catch
t(1) = false;
end
if t(1)
t(2) = all(abs(mu-.4)<1e-6);
t(3) = all(abs(std-.2)<1e-6);
end
T = all(t);
%@eof:2

View File

@ -1,26 +1,48 @@
classdef dprior
classdef dprior < handle
% Copyright © 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/>.
properties
p6 = []; % Prior first hyperparameter.
p7 = []; % Prior second hyperparameter.
p1 = []; % Prior mean.
p2 = []; % Prior stddev.
p3 = []; % Lower bound of the prior support.
p4 = []; % Upper bound of the prior support.
p5 = []; % Prior mode.
p6 = []; % Prior first hyperparameter.
p7 = []; % Prior second hyperparameter.
p11 = []; % Prior median
lb = []; % Truncated prior lower bound.
ub = []; % Truncated prior upper bound.
uniform_index = []; % Index for the uniform priors.
gaussian_index = []; % Index for the gaussian priors.
gamma_index = []; % Index for the gamma priors.
beta_index = []; % Index for the beta priors.
inverse_gamma_1_index = []; % Index for the inverse gamma type 1 priors.
inverse_gamma_2_index = []; % Index for the inverse gamma type 2 priors.
weibull_index = []; % Index for the weibull priors.
uniform_draws = false;
gaussian_draws = false;
gamma_draws = false;
beta_draws = false;
inverse_gamma_1_draws = false;
inverse_gamma_2_draws = false;
weibull_draws = false;
name = {}; % Name of the parameter
iduniform = []; % Index for the uniform priors.
idgaussian = []; % Index for the gaussian priors.
idgamma = []; % Index for the gamma priors.
idbeta = []; % Index for the beta priors.
idinvgamma1 = []; % Index for the inverse gamma type 1 priors.
idinvgamma2 = []; % Index for the inverse gamma type 2 priors.
idweibull = []; % Index for the weibull priors.
isuniform = false;
isgaussian = false;
isgamma = false;
isbeta = false;
isinvgamma1 = false;
isinvgamma2 = false;
isweibull = false;
end
methods
@ -38,10 +60,17 @@ classdef dprior
%
% REQUIREMENTS
% None.
o.p6 = bayestopt_.p6;
o.p7 = bayestopt_.p7;
o.p3 = bayestopt_.p3;
o.p4 = bayestopt_.p4;
if ~nargin % Empty object
return
end
if isfield(bayestopt_, 'p1'), o.p1 = bayestopt_.p1; end
if isfield(bayestopt_, 'p2'), o.p2 = bayestopt_.p2; end
if isfield(bayestopt_, 'p3'), o.p3 = bayestopt_.p3; end
if isfield(bayestopt_, 'p4'), o.p4 = bayestopt_.p4; end
if isfield(bayestopt_, 'p5'), o.p5 = bayestopt_.p5; end
if isfield(bayestopt_, 'p6'), o.p6 = bayestopt_.p6; end
if isfield(bayestopt_, 'p7'), o.p7 = bayestopt_.p7; end
if isfield(bayestopt_, 'p11'), o.p11 = bayestopt_.p11; end
bounds = prior_bounds(bayestopt_, PriorTrunc);
o.lb = bounds.lb;
o.ub = bounds.ub;
@ -50,138 +79,38 @@ classdef dprior
else
prior_shape = bayestopt_.pshape;
end
o.beta_index = find(prior_shape==1);
if ~isempty(o.beta_index)
o.beta_draws = true;
o.idbeta = find(prior_shape==1);
if ~isempty(o.idbeta)
o.isbeta = true;
end
o.gamma_index = find(prior_shape==2);
if ~isempty(o.gamma_index)
o.gamma_draws = true;
o.idgamma = find(prior_shape==2);
if ~isempty(o.idgamma)
o.isgamma = true;
end
o.gaussian_index = find(prior_shape==3);
if ~isempty(o.gaussian_index)
o.gaussian_draws = true;
o.idgaussian = find(prior_shape==3);
if ~isempty(o.idgaussian)
o.isgaussian = true;
end
o.inverse_gamma_1_index = find(prior_shape==4);
if ~isempty(o.inverse_gamma_1_index)
o.inverse_gamma_1_draws = true;
o.idinvgamma1 = find(prior_shape==4);
if ~isempty(o.idinvgamma1)
o.isinvgamma1 = true;
end
o.uniform_index = find(prior_shape==5);
if ~isempty(o.uniform_index)
o.uniform_draws = true;
o.iduniform = find(prior_shape==5);
if ~isempty(o.iduniform)
o.isuniform = true;
end
o.inverse_gamma_2_index = find(prior_shape==6);
if ~isempty(o.inverse_gamma_2_index)
o.inverse_gamma_2_draws = true;
o.idinvgamma2 = find(prior_shape==6);
if ~isempty(o.idinvgamma2)
o.isinvgamma2 = true;
end
o.weibull_index = find(prior_shape==8);
if ~isempty(o.weibull_index)
o.weibull_draws = true;
o.idweibull = find(prior_shape==8);
if ~isempty(o.idweibull)
o.isweibull = true;
end
end
function p = draw(o)
% Return a random draw from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
%
% REMARKS
% None.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> d = Prior.draw()
p = NaN(rows(o.lb), 1);
if o.uniform_draws
p(o.uniform_index) = rand(length(o.uniform_index),1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
while ~isempty(out_of_bound)
p(o.uniform_index) = rand(length(o.uniform_index), 1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
end
end
if o.gaussian_draws
p(o.gaussian_index) = randn(length(o.gaussian_index), 1).*o.p7(o.gaussian_index) + o.p6(o.gaussian_index);
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
while ~isempty(out_of_bound)
p(o.gaussian_index(out_of_bound)) = randn(length(o.gaussian_index(out_of_bound)), 1).*o.p7(o.gaussian_index(out_of_bound)) + o.p6(o.gaussian_index(out_of_bound));
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
end
end
if o.gamma_draws
p(o.gamma_index) = gamrnd(o.p6(o.gamma_index), o.p7(o.gamma_index))+o.p3(o.gamma_index);
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
while ~isempty(out_of_bound)
p(o.gamma_index(out_of_bound)) = gamrnd(o.p6(o.gamma_index(out_of_bound)), o.p7(o.gamma_index(out_of_bound)))+o.p3(o.gamma_index(out_of_bound));
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
end
end
if o.beta_draws
p(o.beta_index) = (o.p4(o.beta_index)-o.p3(o.beta_index)).*betarnd(o.p6(o.beta_index), o.p7(o.beta_index))+o.p3(o.beta_index);
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
while ~isempty(out_of_bound)
p(o.beta_index(out_of_bound)) = (o.p4(o.beta_index(out_of_bound))-o.p3(o.beta_index(out_of_bound))).*betarnd(o.p6(o.beta_index(out_of_bound)), o.p7(o.beta_index(out_of_bound)))+o.p3(o.beta_index(out_of_bound));
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
end
end
if o.inverse_gamma_1_draws
p(o.inverse_gamma_1_index) = ...
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index)/2, 2./o.p6(o.inverse_gamma_1_index)))+o.p3(o.inverse_gamma_1_index);
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
while ~isempty(out_of_bound)
p(o.inverse_gamma_1_index(out_of_bound)) = ...
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_1_index(out_of_bound))))+o.p3(o.inverse_gamma_1_index(out_of_bound));
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
end
end
if o.inverse_gamma_2_draws
p(o.inverse_gamma_2_index) = ...
1./gamrnd(o.p7(o.inverse_gamma_2_index)/2, 2./o.p6(o.inverse_gamma_2_index))+o.p3(o.inverse_gamma_2_index);
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
while ~isempty(out_of_bound)
p(o.inverse_gamma_2_index(out_of_bound)) = ...
1./gamrnd(o.p7(o.inverse_gamma_2_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_2_index(out_of_bound)))+o.p3(o.inverse_gamma_2_index(out_of_bound));
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
end
end
if o.weibull_draws
p(o.weibull_index) = wblrnd(o.p7(o.weibull_index), o.p6(o.weibull_index)) + o.p3(o.weibull_index);
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
while ~isempty(out_of_bound)
p(o.weibull_index(out_of_bound)) = wblrnd(o.p7(o.weibull_index(out_of_bound)), o.p6(o.weibull_index(out_of_bound)))+o.p3(o.weibull_index(out_of_bound));
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
end
end
end
function P = draws(o, n)
% Return n independent random draws from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - P [double] m×n matrix, random draw from the prior distribution.
%
% REMARKS
% If the Parallel Computing Toolbox is available, the main loop is run in parallel.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> Prior.draws(1e6)
P = NaN(rows(o.lb), 1);
parfor i=1:n
P(:,i) = draw(o);
end
end
end % dprior (constructor)
end % methods
end % classdef --*-- Unit tests --*--
%@test:1
@ -263,114 +192,22 @@ end % classdef --*-- Unit tests --*--
%$ try
%$ % Instantiate dprior object
%$ o = dprior(bayestopt_, prior_trunc, false);
%$ % Do simulations in a loop and estimate recursively the mean and the variance.
%$ for i = 1:ndraws
%$ draw = o.draw();
%$ m1 = m0 + (draw-m0)/i;
%$ m2 = m1*m1';
%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
%$ m0 = m1;
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = all(abs(m0-bayestopt_.p1)<3e-3);
%$ t(3) = all(all(abs(v0-diag(bayestopt_.p2.^2))<5e-3));
%$ end
%$ T = all(t);
%@eof:1
%@test:2
%$ % Fill global structures with required fields...
%$ prior_trunc = 1e-10;
%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
%$ p1 = .4*ones(14,1); % Prior mean
%$ p2 = .2*ones(14,1); % Prior std.
%$ p3 = NaN(14,1);
%$ p4 = NaN(14,1);
%$ p5 = NaN(14,1);
%$ p6 = NaN(14,1);
%$ p7 = NaN(14,1);
%$
%$ for i=1:14
%$ switch p0(i)
%$ case 1
%$ % Beta distribution
%$ p3(i) = 0;
%$ p4(i) = 1;
%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
%$ case 2
%$ % Gamma distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
%$ case 3
%$ % Normal distribution
%$ p3(i) = -Inf;
%$ p4(i) = Inf;
%$ p6(i) = p1(i);
%$ p7(i) = p2(i);
%$ p5(i) = p1(i);
%$ case 4
%$ % Inverse Gamma (type I) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
%$ case 5
%$ % Uniform distribution
%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
%$ p3(i) = p6(i);
%$ p4(i) = p7(i);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
%$ case 6
%$ % Inverse Gamma (type II) distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
%$ case 8
%$ % Weibull distribution
%$ p3(i) = 0;
%$ p4(i) = Inf;
%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
%$ otherwise
%$ error('This density is not implemented!')
%$ end
%$ end
%$
%$ bayestopt_.pshape = p0;
%$ bayestopt_.p1 = p1;
%$ bayestopt_.p2 = p2;
%$ bayestopt_.p3 = p3;
%$ bayestopt_.p4 = p4;
%$ bayestopt_.p5 = p5;
%$ bayestopt_.p6 = p6;
%$ bayestopt_.p7 = p7;
%$
%$ ndraws = 1e5;
%$
%$ % Call the tested routine
%$ try
%$ % Instantiate dprior object.
%$ o = dprior(bayestopt_, prior_trunc, false);
%$ X = o.draws(ndraws);
%$ m = mean(X, 2);
%$ v = var(X, 0, 2);
%$ % Instantiate dprior object
%$ o = dprior();
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = all(abs(m-bayestopt_.p1)<3e-3);
%$ t(3) = all(all(abs(v-bayestopt_.p2.^2)<5e-3));
%$ end
%$ T = all(t);
%@eof:2

197
matlab/@dprior/draw.m Normal file
View File

@ -0,0 +1,197 @@
function p = draw(o)
% Return a random draw from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
%
% REMARKS
% None.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> d = Prior.draw()
% Copyright © 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/>.
p = NaN(rows(o.lb), 1);
if o.isuniform
p(o.iduniform) = rand(length(o.iduniform),1).*(o.p4(o.iduniform)-o.p3(o.iduniform)) + o.p3(o.iduniform);
oob = find( (p(o.iduniform)>o.ub(o.iduniform)) | (p(o.iduniform)<o.lb(o.iduniform)));
while ~isempty(oob)
p(o.iduniform) = rand(length(o.iduniform), 1).*(o.p4(o.iduniform)-o.p3(o.iduniform)) + o.p3(o.iduniform);
oob = find( (p(o.iduniform)>o.ub(o.iduniform)) | (p(o.iduniform)<o.lb(o.iduniform)));
end
end
if o.isgaussian
p(o.idgaussian) = randn(length(o.idgaussian), 1).*o.p7(o.idgaussian) + o.p6(o.idgaussian);
oob = find( (p(o.idgaussian)>o.ub(o.idgaussian)) | (p(o.idgaussian)<o.lb(o.idgaussian)));
while ~isempty(oob)
p(o.idgaussian(oob)) = randn(length(o.idgaussian(oob)), 1).*o.p7(o.idgaussian(oob)) + o.p6(o.idgaussian(oob));
oob = find( (p(o.idgaussian)>o.ub(o.idgaussian)) | (p(o.idgaussian)<o.lb(o.idgaussian)));
end
end
if o.isgamma
p(o.idgamma) = gamrnd(o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
oob = find( (p(o.idgamma)>o.ub(o.idgamma)) | (p(o.idgamma)<o.lb(o.idgamma)));
while ~isempty(oob)
p(o.idgamma(oob)) = gamrnd(o.p6(o.idgamma(oob)), o.p7(o.idgamma(oob)))+o.p3(o.idgamma(oob));
oob = find( (p(o.idgamma)>o.ub(o.idgamma)) | (p(o.idgamma)<o.lb(o.idgamma)));
end
end
if o.isbeta
p(o.idbeta) = (o.p4(o.idbeta)-o.p3(o.idbeta)).*betarnd(o.p6(o.idbeta), o.p7(o.idbeta))+o.p3(o.idbeta);
oob = find( (p(o.idbeta)>o.ub(o.idbeta)) | (p(o.idbeta)<o.lb(o.idbeta)));
while ~isempty(oob)
p(o.idbeta(oob)) = (o.p4(o.idbeta(oob))-o.p3(o.idbeta(oob))).*betarnd(o.p6(o.idbeta(oob)), o.p7(o.idbeta(oob)))+o.p3(o.idbeta(oob));
oob = find( (p(o.idbeta)>o.ub(o.idbeta)) | (p(o.idbeta)<o.lb(o.idbeta)));
end
end
if o.isinvgamma1
p(o.idinvgamma1) = ...
sqrt(1./gamrnd(o.p7(o.idinvgamma1)/2, 2./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
oob = find( (p(o.idinvgamma1)>o.ub(o.idinvgamma1)) | (p(o.idinvgamma1)<o.lb(o.idinvgamma1)));
while ~isempty(oob)
p(o.idinvgamma1(oob)) = ...
sqrt(1./gamrnd(o.p7(o.idinvgamma1(oob))/2, 2./o.p6(o.idinvgamma1(oob))))+o.p3(o.idinvgamma1(oob));
oob = find( (p(o.idinvgamma1)>o.ub(o.idinvgamma1)) | (p(o.idinvgamma1)<o.lb(o.idinvgamma1)));
end
end
if o.isinvgamma2
p(o.idinvgamma2) = ...
1./gamrnd(o.p7(o.idinvgamma2)/2, 2./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
oob = find( (p(o.idinvgamma2)>o.ub(o.idinvgamma2)) | (p(o.idinvgamma2)<o.lb(o.idinvgamma2)));
while ~isempty(oob)
p(o.idinvgamma2(oob)) = ...
1./gamrnd(o.p7(o.idinvgamma2(oob))/2, 2./o.p6(o.idinvgamma2(oob)))+o.p3(o.idinvgamma2(oob));
oob = find( (p(o.idinvgamma2)>o.ub(o.idinvgamma2)) | (p(o.idinvgamma2)<o.lb(o.idinvgamma2)));
end
end
if o.isweibull
p(o.idweibull) = wblrnd(o.p7(o.idweibull), o.p6(o.idweibull)) + o.p3(o.idweibull);
oob = find( (p(o.idweibull)>o.ub(o.idweibull)) | (p(o.idweibull)<o.lb(o.idweibull)));
while ~isempty(oob)
p(o.idweibull(oob)) = wblrnd(o.p7(o.idweibull(oob)), o.p6(o.idweibull(oob)))+o.p3(o.idweibull(oob));
oob = find( (p(o.idweibull)>o.ub(o.idweibull)) | (p(o.idweibull)<o.lb(o.idweibull)));
end
end
return % --*-- Unit tests --*--
%@test:1
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
ndraws = 1e5;
m0 = BayesInfo.p1; %zeros(14,1);
v0 = diag(BayesInfo.p2.^2); %zeros(14);
% Call the tested routine
try
% Instantiate dprior object
o = dprior(BayesInfo, prior_trunc, false);
% Do simulations in a loop and estimate recursively the mean and the variance.
for i = 1:ndraws
draw = o.draw();
m1 = m0 + (draw-m0)/i;
m2 = m1*m1';
v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
m0 = m1;
end
t(1) = true;
catch
t(1) = false;
end
if t(1)
t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
end
T = all(t);
%@eof:1

133
matlab/@dprior/draws.m Normal file
View File

@ -0,0 +1,133 @@
function P = draws(o, n)
% Return n independent random draws from the prior distribution.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - P [double] m×n matrix, random draw from the prior distribution.
%
% REMARKS
% If the Parallel Computing Toolbox is available, the main loop is run in parallel.
%
% EXAMPLE
%
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
% >> Prior.draws(1e6)
% Copyright © 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/>.
P = NaN(rows(o.lb), 1);
parfor i=1:n
P(:,i) = draw(o);
end
return % --*-- Unit tests --*--
%@test:1
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
ndraws = 1e5;
% Call the tested routine
try
% Instantiate dprior object.
o = dprior(BayesInfo, prior_trunc, false);
X = o.draws(ndraws);
m = mean(X, 2);
v = var(X, 0, 2);
t(1) = true;
catch
t(1) = false;
end
if t(1)
t(2) = all(abs(m-BayesInfo.p1)<3e-3);
t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3));
end
T = all(t);
%@eof:1

31
matlab/@dprior/length.m Normal file
View File

@ -0,0 +1,31 @@
function n = length(o)
% Return the dimension of the random vector.
%
% INPUTS
% - o [dprior] Distribution specification for a n×1 vector of independent continuous random variables
%
% OUTPUTS
% - n [integer] scalar.
%
% REMARKS
% None.
% Copyright © 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/>.
n = length(o.lb);

42
matlab/@dprior/mean.m Normal file
View File

@ -0,0 +1,42 @@
function m = mean(o, resetmoments)
% Return the prior mean.
%
% INPUTS
% - o [dprior]
% - resetmoments [logical] Force the computation of the prior mean
%
% OUTPUTS
% - m [double] n×1 vector, prior mean
% Copyright © 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/>.
if nargin<2
resetmoments = false;
end
if any(isnan(o.p1))
resetmoments = true;
end
if resetmoments
o.p1 = NaN(size(o.p1));
o.moments('mean');
end
m = o.p1;

42
matlab/@dprior/median.m Normal file
View File

@ -0,0 +1,42 @@
function m = median(o, resetmoments)
% Return the prior median.
%
% INPUTS
% - o [dprior]
% - resetmoments [logical] Force the computation of the prior median
%
% OUTPUTS
% - m [double] n×1 vector, prior median
% Copyright © 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/>.
if nargin<2
resetmoments = false;
end
if any(isnan(o.p11))
resetmoments = true;
end
if resetmoments
o.p11 = NaN(size(o.p11));
o.moments('median');
end
m = o.p11;

42
matlab/@dprior/mode.m Normal file
View File

@ -0,0 +1,42 @@
function m = mode(o, resetmoments)
% Return the prior mode.
%
% INPUTS
% - o [dprior]
% - resetmoments [logical] Force the computation of the prior mode
%
% OUTPUTS
% - m [double] n×1 vector, prior mode
% Copyright © 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/>.
if nargin<2
resetmoments = false;
end
if any(isnan(o.p5))
resetmoments = true;
end
if resetmoments
o.p5 = NaN(size(o.p5));
o.moments('mode');
end
m = o.p5;

291
matlab/@dprior/moments.m Normal file
View File

@ -0,0 +1,291 @@
function o = moments(o, name)
% Compute the prior moments.
%
% INPUTS
% - o [dprior]
%
% OUTPUTS
% - o [dprior]
% Copyright © 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/>.
switch name
case 'mean'
m = o.p1;
case 'median'
m = o.p11;
case 'std'
m = o.p2;
case 'mode'
m = o.p5;
otherwise
error('%s is not an implemented moemnt.', name)
end
id = isnan(m);
if any(id)
% For some parameters the prior mean is not defined.
% We compute the first order moment from the
% hyperparameters, if the hyperparameters are not
% available an error is thrown.
if o.isuniform
jd = intersect(o.iduniform, find(id));
if ~isempty(jd)
if any(isnan(o.p3(jd))) || any(isnan(o.p4(jd)))
error('dprior::mean: Some hyperparameters are missing (uniform distribution).')
end
switch name
case 'mean'
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
case 'median'
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
case 'std'
m(jd) = (o.p4(jd)-o.p3(jd))/sqrt(12);
case 'mode' % Actually we have a continuum of modes with the uniform distribution.
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
end
end
end
if o.isgaussian
jd = intersect(o.idgaussian, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd)))
error('dprior::mean: Some hyperparameters are missing (gaussian distribution).')
end
switch name
case 'mean'
m(jd) = o.p6(jd);
case 'median'
m(jd) = o.p6(jd);
case 'std'
m(jd) = o.p7(jd);
case 'mode' % Actually we have a continuum of modes with the uniform distribution.
m(jd) = o.p6(jd);
end
end
end
if o.isgamma
jd = intersect(o.idgamma, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
error('dprior::mean: Some hyperparameters are missing (gamma distribution).')
end
% α → o.p6, β → o.p7
switch name
case 'mean'
m(jd) = o.p3(jd) + o.p6(jd).*o.p7(jd);
case 'median'
m(jd) = o.p3(jd) + gaminv(.5, o.p6(jd), o.p7(jda));
case 'std'
m(jd) = sqrt(o.p6(jd)).*o.p7(jd);
case 'mode'
m(jd) = 0;
hd = o.p6(jd)>1;
m(jd(hd)) = (o.p6(jd(hd))-1).*o.p7(jd(hd));
end
end
end
if o.isbeta
jd = intersect(o.idbeta, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd))) || any(isnan(o.p4(jd)))
error('dprior::mean: Some hyperparameters are missing (beta distribution).')
end
% α → o.p6, β → o.p7
switch name
case 'mean'
m(jd) = o.p3(jd) + (o.p6(jd)./(o.p6(jd)+o.p7(jd))).*(o.p4(jd)-o.p3(jd));
case 'median'
m(jd) = o.p3(jd) + betainv(.5, o.p6(jd), o.p7(jd)).*(o.p4(jd)-o.p3(jd));
case 'std'
m(jd) = (o.p4(jd)-o.p3(jd)).*sqrt(o.p6(jd).*o.p7(jd)./((o.p6(jd)+o.p7(jd)).^2.*(o.p6(jd)+o.p7(jd)+1)));
case 'mode'
h0 = true(jd, 1);
h1 = o.p6(jd)<=1 & o.p7(jd)>1; h0 = h0 & ~h1;
h2 = o.p7(jd)<=1 & o.p6(jd)>1; h0 = h0 & ~h2;
h3 = o.p6(jd)<1 & o.p7(jd)<1; h0 = h0 & ~h3;
h4 = ismembertol(o.p6(jd), 1) & ismembertol(o.p7(jd),1); h0 = h0 & ~h4;
m(jd(h1)) = o.p3(jd(h1)); % Standard β has a mode at 0
m(jd(h2)) = o.p4(jd(h2)); % Standard β has a mode at 1
m(jd(h3)) = o.p3(jd(h3)); % Standard β is bimodal, we pick the lowest mode (0)
m(jd(h4)) = o.p3(jd(h4)) + .5*(o.p4(jd(h4))-o.p3(jd(h4))); % Standard β is the uniform distribution (continuum of modes), we pick the mean as the mode
m(jd(h0)) = o.p3(jd(h0))+(o.p4(jd(h0))-o.p3(jd(h0))).*((o.p6(jd(h0))-1)./(o.p6(jd(h0))+o.p7(jd(h0))-2)); % β distribution is concave and has a unique interior mode.
end
end
end
if o.isinvgamma1
jd = intersect(o.idinvgamma1, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
error('dprior::mean: Some hyperparameters are missing (inverse gamma type 1 distribution).')
end
% s → o.p6, ν → o.p7
switch name
case 'mean'
m(jd) = o.p3(jd) + sqrt(.5*o.p6(jd)) .*(gamma(.5*(o.p7(jd)-1))./gamma(.5*o.p7(jd)));
case 'median'
m(jd) = o.p3(jd) + 1.0/sqrt(gaminv(.5, o.p7(jd)/2.0, 2.0/o.p6(jd)));
case 'std'
m(jd) = sqrt( o.p6(jd)./(o.p7(jd)-2)-(.5*o.p6(jd)).*(gamma(.5*(o.p7(jd)-1))./gamma(.5*o.p7(jd))).^2);
case 'mode'
m(jd) = sqrt((o.p7(jd)-1)./o.p6(jd));
end
end
end
if o.isinvgamma2
jd = intersect(o.idinvgamma2, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
error('dprior::mean: Some hyperparameters are missing (inverse gamma type 2 distribution).')
end
% s → o.p6, ν → o.p7
switch name
case 'mean'
m(jd) = o.p3(jd) + o.p6(jd)./(o.p7(jd)-2);
case 'median'
m(jd) = o.p3(jd) + 1.0/gaminv(.5, o.p7(jd)/2.0, 2.0/o.p6(jd));
case 'std'
m(jd) = sqrt(2./(o.p7(jd)-4)).*o.p6(jd)./(o.p7(jd)-2);
case 'mode'
m(jd) = o.p6(jd)./(o.p7(jd)+2);
end
end
end
if o.isweibull
jd = intersect(o.idweibull, find(id));
if ~isempty(jd)
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
error('dprior::mean: Some hyperparameters are missing (weibull distribution).')
end
% k → o.p6 (shape parameter), λ → o.p7 (scale parameter)
% See https://en.wikipedia.org/wiki/Weibull_distribution
switch name
case 'mean'
m(jd) = o.p3(jd) + o.p7(jd).*gamma(1+1./o.p6(jd));
case 'median'
m(jd) = o.p3(jd) + o.p7(jd).*log(2).^(1./o.p6(jd));
case 'std'
m(jd) = o.p7(jd).*sqrt(gamma(1+2./o.p6(jd))-gamma(1+1./o.p6(jd)).^2);
case 'mode'
m(jd) = 0;
hd = o.p6(jd)>1;
m(jd(hd)) = o.p3(jd(hd)) + o.p7(jd(hd)).*((o.p6(jd(hd))-1)./o.p6(jd(hd))).^(1./o.p6(jd(hd)));
end
end
end
switch name
case 'mean'
o.p1 = m;
case 'median'
o.p11 = m;
case 'std'
o.p2 = m;
case 'mode'
o.p5 = m;
end
end
return % --*-- Unit tests --*--
%@test:5
% Fill global structures with required fields...
prior_trunc = 1e-10;
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
p1 = .4*ones(14,1); % Prior mean
p2 = .2*ones(14,1); % Prior std.
p3 = NaN(14,1);
p4 = NaN(14,1);
p5 = NaN(14,1);
p6 = NaN(14,1);
p7 = NaN(14,1);
for i=1:14
switch p0(i)
case 1
% Beta distribution
p3(i) = 0;
p4(i) = 1;
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
case 2
% Gamma distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
case 3
% Normal distribution
p3(i) = -Inf;
p4(i) = Inf;
p6(i) = p1(i);
p7(i) = p2(i);
p5(i) = p1(i);
case 4
% Inverse Gamma (type I) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
case 5
% Uniform distribution
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
p3(i) = p6(i);
p4(i) = p7(i);
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
case 6
% Inverse Gamma (type II) distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
case 8
% Weibull distribution
p3(i) = 0;
p4(i) = Inf;
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
otherwise
error('This density is not implemented!')
end
end
BayesInfo.pshape = p0;
BayesInfo.p1 = p1;
BayesInfo.p2 = p2;
BayesInfo.p3 = p3;
BayesInfo.p4 = p4;
BayesInfo.p5 = p5;
BayesInfo.p6 = p6;
BayesInfo.p7 = p7;
% Call the tested routine
try
Prior = dprior(BayesInfo, prior_trunc, false);
t(1) = true;
catch
t(1) = false;
end
if t(1)
t(2) = all(Prior.mean()==.4);
t(3) = all(ismembertol(Prior.mean(true),.4));
t(4) = all(ismembertol(Prior.variance(),.04));
t(5) = all(ismembertol(Prior.variance(true),.04));
end
T = all(t);
%@eof:5

41
matlab/@dprior/subsref.m Normal file
View File

@ -0,0 +1,41 @@
function p = subsref(o, S)
% Overload subsref method.
% Copyright © 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/>.
switch S(1).type
case '.'
if ismember(S(1).subs, {'p1','p2','p3','p4','p5','p6','p7','lb','ub'})
p = builtin('subsref', o, S(1));
elseif ismember(S(1).subs, {'draw','length'})
p = feval(S(1).subs, o);
elseif ismember(S(1).subs, {'draws', 'density', 'densities', 'moments', 'admissible'})
p = feval(S(1).subs, o , S(2).subs{:});
elseif ismember(S(1).subs, {'mean', 'median', 'variance', 'mode'})
if (length(S)==2 && isempty(S(2).subs)) || length(S)==1
p = feval(S(1).subs, o);
else
p = feval(S(1).subs, o , S(2).subs{:});
end
else
error('dprior::subsref: unknown method (%s).', S(1).subs)
end
otherwise
error('dprior::subsref: %s indexing not implemented.', S(1).type)
end

42
matlab/@dprior/variance.m Normal file
View File

@ -0,0 +1,42 @@
function m = variance(o, resetmoments)
% Return the prior variance.
%
% INPUTS
% - o [dprior]
% - resetmoments [logical] Force the computation of the prior variance
%
% OUTPUTS
% - m [double] n×1 vector, prior variance
% Copyright © 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/>.
if nargin<2
resetmoments = false;
end
if any(isnan(o.p2))
resetmoments = true;
end
if resetmoments
o.p2 = NaN(size(o.p2));
o.moments('std');
end
m = o.p2.^2;

View File

@ -1,152 +0,0 @@
function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
%
% computes the asymptotic hessian matrix of the log-likelihood function of
% a state space model (notation as in kalman_filter.m in DYNARE
% Thanks to Nikolai Iskrev
%
% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
% dimension equal to the number of structural parameters
% Copyright © 2011-2017 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/>.
k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size.
pp = size(Y,1); % Maximum number of observed variables.
mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
oldK = 0;
notsteady = 1; % Steady state flag.
F_singular = 1;
lik = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
if nargout > 1
DLIK = zeros(k,1); % Initialization of the score.
end
AHess = zeros(k,k); % Initialization of the Hessian
Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k);
% for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady && t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+Om;
end
else
F_singular = 0;
iF = inv(F);
K = P(:,mf)*iF;
lik(t) = log(det(F))+transpose(v)*iF*v;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout > 1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
vecDPmf = reshape(DP(mf,mf,:),[],k);
% iPmf = inv(P(mf,mf));
if t>=start
AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
end
notsteady = max(max(abs(K-oldK))) > riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start && nargout >1
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
if t>=start
AHess = AHess + Dv'*iF*Dv;
end
a = T*(a+K*v);
lik(t) = transpose(v)*iF*v;
end
AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
if nargout > 1
for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
end
end
lik(t0:smpl) = lik(t0:smpl) + log(det(F));
% for ii = 1:k;
% for jj = 1:ii
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
% end
% end
end
AHess = -AHess;
if nargout > 1
DLIK = DLIK/2;
end
% adding log-likelihhod constants
lik = (lik + pp*log(2*pi))/2;
LIK = sum(lik(start:end)); % Minus the log-likelihood.
% end of main function
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
k = size(DT,3);
tmp = P-K*P(mf,:);
for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
% end of computeDKalman

View File

@ -275,7 +275,7 @@ if kalman_algo == 1 || kalman_algo == 3
kalman_algo = 2;
elseif kalman_algo == 3
fprintf('\nDsgeSmoother: Switching to univariate filter. This is usually due to co-integration in diffuse filter,\n')
fprintf(' otherwise it may be a sign of stochastic singularity.\n')
fprintf('otherwise it may be a sign of stochastic singularity.\n')
kalman_algo = 4;
else
error('This case shouldn''t happen')

View File

@ -1,23 +1,24 @@
function Draws = GetAllPosteriorDraws(dname, fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
% function Draws = GetAllPosteriorDraws(dname, fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
% Gets all posterior draws
function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
% Gets all posterior draws.
%
% INPUTS
% dname: name of directory with results
% fname: name of mod file
% column: column of desired parameter in draw matrix
% FirstMhFile: first mh file
% FirstLine: first line
% TotalNumberOfMhFile: total number of mh file
% NumberOfDraws: number of draws
% nblcks: total number of blocks
% blck: desired block to read
% - options_ [struct] Dynare's options.
% - dname [char] name of directory with results.
% - fname [char] name of mod file.
% - column [integer] scalar, parameter index.
% - FirstMhFile [integer] scalar, first MH file.
% - FirstLine [integer] scalar, first line in first MH file.
% - TotalNumberOfMhFile [integer] scalar, total number of MH file.
% - NumberOfDraws [integer] scalar, number of posterior draws.
% - nblcks [integer] scalar, total number of blocks.
% - blck: [integer] scalar, desired block to read.
%
% OUTPUTS
% Draws: draws from posterior distribution
% - draws: [double] NumberOfDraws×1 vector, draws from posterior distribution.
%
% SPECIAL REQUIREMENTS
% none
% REMARKS
% Only the first and third input arguments are required for SMC samplers.
% Copyright © 2005-2023 Dynare Team
%
@ -36,55 +37,61 @@ function Draws = GetAllPosteriorDraws(dname, fname, column, FirstMhFile, FirstLi
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
iline = FirstLine;
linee = 1;
DirectoryName = CheckPath('metropolis',dname);
if nblcks>1 && nargin<9
Draws = zeros(NumberOfDraws*nblcks,1);
iline0=iline;
if column>0
for blck = 1:nblcks
iline=iline0;
if ishssmc(options_)
% Load draws from the posterior distribution
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,:));
else
iline = FirstLine;
linee = 1;
DirectoryName = CheckPath('metropolis',dname);
if nblcks>1 && nargin<10
draws = zeros(NumberOfDraws*nblcks,1);
iline0=iline;
if column>0
for blck = 1:nblcks
iline=iline0;
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
NumberOfLines = size(x2(iline:end,:),1);
draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
linee = linee+NumberOfLines;
iline = 1;
end
end
else
for blck = 1:nblcks
iline=iline0;
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
NumberOfLines = size(logpo2(iline:end),1);
draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
linee = linee+NumberOfLines;
iline = 1;
end
end
end
else
if nblcks==1
blck=1;
end
if column>0
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
NumberOfLines = size(x2(iline:end,:),1);
Draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
linee = linee+NumberOfLines;
iline = 1;
end
end
else
for blck = 1:nblcks
iline=iline0;
else
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
NumberOfLines = size(logpo2(iline:end),1);
Draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
NumberOfLines = size(logpo2(iline:end,:),1);
draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
linee = linee+NumberOfLines;
iline = 1;
end
end
end
else
if nblcks==1
blck=1;
end
if column>0
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
NumberOfLines = size(x2(iline:end,:),1);
Draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
linee = linee+NumberOfLines;
iline = 1;
end
else
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'logpo2')
NumberOfLines = size(logpo2(iline:end,:),1);
Draws(linee:linee+NumberOfLines-1) = logpo2(iline:end);
linee = linee+NumberOfLines;
iline = 1;
end
end
end
end

View File

@ -1,18 +1,15 @@
function [mean,variance] = GetPosteriorMeanVariance(M_,drop)
function [mean, variance] = GetPosteriorMeanVariance(options_, M_)
%function [mean,variance] = GetPosteriorMeanVariance(M,drop)
% Computes the posterior mean and variance
% (+updates of oo_ & TeX output).
%
% INPUTS
% M_ [structure] Dynare model structure
% drop [double] share of draws to drop
% - options_ [struct] Dynare's options.
% - M_ [struct] Description of the model.
%
% OUTPUTS
% mean [double] mean parameter vector
% variance [double] variance
%
% SPECIAL REQUIREMENTS
% None.
% - mean [double] n×1 vector, posterior expectation.
% - variance [double] n×n matrix, posterior variance.
% Copyright © 2012-2023 Dynare Team
%
@ -31,37 +28,46 @@ function [mean,variance] = GetPosteriorMeanVariance(M_,drop)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName];
record=load_last_mh_history_file(MetropolisFolder, FileName);
NbrDraws = sum(record.MhDraws(:,1));
NbrFiles = sum(record.MhDraws(:,2));
NbrBlocks = record.Nblck;
mean = 0;
variance = 0;
NbrKeptDraws = 0;
for i=1:NbrBlocks
NbrDrawsCurrentBlock = 0;
for j=1:NbrFiles
o = load([BaseName '_mh' int2str(j) '_blck' int2str(i),'.mat']);
NbrDrawsCurrentFile = size(o.x2,1);
if NbrDrawsCurrentBlock + NbrDrawsCurrentFile <= drop*NbrDraws
if ishssmc(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
% Compute the posterior mean
mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
else
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName];
record=load_last_mh_history_file(MetropolisFolder, FileName);
NbrDraws = sum(record.MhDraws(:,1));
NbrFiles = sum(record.MhDraws(:,2));
NbrBlocks = record.Nblck;
mean = 0;
variance = 0;
NbrKeptDraws = 0;
for i=1:NbrBlocks
NbrDrawsCurrentBlock = 0;
for j=1:NbrFiles
o = load([BaseName '_mh' int2str(j) '_blck' int2str(i),'.mat']);
NbrDrawsCurrentFile = size(o.x2,1);
if NbrDrawsCurrentBlock + NbrDrawsCurrentFile <= options_.mh_drop*NbrDraws
NbrDrawsCurrentBlock = NbrDrawsCurrentBlock + NbrDrawsCurrentFile;
continue
elseif NbrDrawsCurrentBlock < options_.mh_drop*NbrDraws
FirstDraw = ceil(options_.mh_drop*NbrDraws - NbrDrawsCurrentBlock + 1);
x2 = o.x2(FirstDraw:end,:);
else
x2 = o.x2;
end
NbrKeptDrawsCurrentFile = size(x2,1);
%recursively compute mean and variance
mean = (NbrKeptDraws*mean + sum(x2)')/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
x2Demeaned = bsxfun(@minus,x2,mean');
variance = (NbrKeptDraws*variance + x2Demeaned'*x2Demeaned)/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
NbrDrawsCurrentBlock = NbrDrawsCurrentBlock + NbrDrawsCurrentFile;
continue
elseif NbrDrawsCurrentBlock < drop*NbrDraws
FirstDraw = ceil(drop*NbrDraws - NbrDrawsCurrentBlock + 1);
x2 = o.x2(FirstDraw:end,:);
else
x2 = o.x2;
NbrKeptDraws = NbrKeptDraws + NbrKeptDrawsCurrentFile;
end
NbrKeptDrawsCurrentFile = size(x2,1);
%recursively compute mean and variance
mean = (NbrKeptDraws*mean + sum(x2)')/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
x2Demeaned = bsxfun(@minus,x2,mean');
variance = (NbrKeptDraws*variance + x2Demeaned'*x2Demeaned)/(NbrKeptDraws+NbrKeptDrawsCurrentFile);
NbrDrawsCurrentBlock = NbrDrawsCurrentBlock + NbrDrawsCurrentFile;
NbrKeptDraws = NbrKeptDraws + NbrKeptDrawsCurrentFile;
end
end

View File

@ -1,5 +1,5 @@
function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output).
%
@ -34,10 +34,6 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
%if ~options_.mh_replic && options_.load_mh_file
% load([M_.fname '_results.mat'],'oo_');
%end
TeX = options_.TeX;
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@ -45,19 +41,20 @@ ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
MetropolisFolder = CheckPath('metropolis',M_.dname);
latexFolder = CheckPath('latex',M_.dname);
FileName = M_.fname;
record=load_last_mh_history_file(MetropolisFolder,FileName);
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);
mh_nblck = size(record.LastParameters,1);
clear record;
if ~issmc(options_)
MetropolisFolder = CheckPath('metropolis',M_.dname);
record=load_last_mh_history_file(MetropolisFolder,FileName);
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);
mh_nblck = size(record.LastParameters,1);
clear record;
end
header_width = row_header_width(M_, estim_params_, bayestopt_);
hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
@ -68,13 +65,26 @@ skipline(2)
disp('ESTIMATION RESULTS')
skipline()
if ~isfield(oo_,'MarginalDensity') || ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean')
[~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
% Set function handle for GetAllPosteriorDraws
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_);
end
fprintf('Log data density is %f.', oo_.MarginalDensity.ModifiedHarmonicMean);
% Set function handle for GetAllPosteriordraws
getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, M_.fname, i, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
end
fprintf('\nLog data density is %f.\n', oo_.MarginalDensity.ModifiedHarmonicMean);
num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
if ishssmc(options_)
num_draws = options_.posterior_sampler_options.hssmc.particles;
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);
end
if hpd_draws<2
fprintf('posterior_moments: There are not enough draws computes to compute HPD Intervals. Skipping their computation.\n')
@ -93,9 +103,9 @@ if np
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
else
@ -103,8 +113,8 @@ if np
name = bayestopt_.name{ip};
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
name = bayestopt_.name{ip};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
end
@ -137,8 +147,8 @@ if nvx
disp(tit2)
for i=1:nvx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
name = M_.exo_names{k};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@ -149,9 +159,8 @@ if nvx
name = M_.exo_names{k};
[post_mean, hpd_interval, post_var] = Extractoo(oo_, name, type);
catch
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
k = estim_params_.var_exo(i,1);
name = M_.exo_names{k};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
@ -181,8 +190,8 @@ if nvn
ip = nvx+1;
for i=1:nvn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_, name, type, post_mean, hpd_interval, post_median, post_var, post_deciles, density);
else
@ -190,8 +199,8 @@ if nvn
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
[post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
catch
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
@ -220,8 +229,8 @@ if ncx
ip = nvx+nvn+1;
for i=1:ncx
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1,options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws,1,options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@ -237,8 +246,8 @@ if ncx
NAME = sprintf('%s_%s', M_.exo_names{k1}, M_.exo_names{k2});
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
k1 = estim_params_.corrx(i,1);
k2 = estim_params_.corrx(i,2);
name = sprintf('%s,%s', M_.exo_names{k1}, M_.exo_names{k2});
@ -259,6 +268,7 @@ if ncx
TeXEnd(fid, 4, 'correlation of structural shocks');
end
end
if ncn
type = 'measurement_errors_corr';
if TeX
@ -270,8 +280,8 @@ if ncn
ip = nvx+nvn+ncx+1;
for i=1:ncn
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
@ -285,8 +295,8 @@ if ncn
NAME = sprintf('%s_%s', M_.endo_names{k1}, M_.endo_names{k2});
[post_mean,hpd_interval,post_var] = Extractoo(oo_, NAME, type);
catch
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, ip, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws, 1, options_.mh_conf_sig);
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, 1, options_.mh_conf_sig);
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
name = sprintf('%s,%s', M_.endo_names{k1}, M_.endo_names{k2});
@ -333,11 +343,8 @@ fprintf(fidTeX, ' & \\multicolumn{3}{c}{Prior} & \\multicolumn{4}{c}{Posterio
fprintf(fidTeX, ' \\cmidrule(r{.75em}){2-4} \\cmidrule(r{.75em}){5-8}\n');
fprintf(fidTeX, ' & Dist. & Mean & Stdev. & Mean & Stdev. & HPD inf & HPD sup\\\\\n');
fprintf(fidTeX, '\\midrule \\endhead \n');
fprintf(fidTeX, '\\bottomrule \\multicolumn{8}{r}{(Continued on next page)} \\endfoot \n');
fprintf(fidTeX, '\\bottomrule \\endlastfoot \n');
fid = fidTeX;
@ -375,4 +382,4 @@ hpd_interval = zeros(2,1);
post_mean = oo.posterior_mean.(type).(name);
hpd_interval(1) = oo.posterior_hpdinf.(type).(name);
hpd_interval(2) = oo.posterior_hpdsup.(type).(name);
post_var = oo.posterior_variance.(type).(name);
post_var = oo.posterior_variance.(type).(name);

View File

@ -1,246 +0,0 @@
function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% function Herbst_Schorfheide_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% SMC sampler from JAE 2014 .
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
% function (posterior kernel).
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% 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
%
% 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 posterior_sampler_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 © 2006-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/>.
% Create the tempering schedule
phi = bsxfun(@power,(bsxfun(@minus,1:1:options_.posterior_sampler_options.HSsmc.nphi,1)/(options_.posterior_sampler_options.HSsmc.nphi-1)),options_.posterior_sampler_options.HSsmc.lambda) ;
% tuning for MH algorithms matrices
zhat = 0 ; % normalization constant
csim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % scale parameter
ESSsim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % ESS
acptsim = zeros(options_.posterior_sampler_options.HSsmc.nphi,1) ; % average acceptance rate
% Step 0: Initialization of the sampler
[ param, tlogpost_i, loglik, bayestopt_] = ...
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.posterior_sampler_options.HSsmc.nparticles);
weights = ones(options_.posterior_sampler_options.HSsmc.nparticles,1)/options_.posterior_sampler_options.HSsmc.nparticles ;
% The Herbst and Schorfheide sampler starts here
for i=2:options_.posterior_sampler_options.HSsmc.nphi
% (a) Correction
% incremental weights
incwt = exp((phi(i)-phi(i-1))*loglik) ;
% update weights
weights = bsxfun(@times,weights,incwt) ;
sum_weights = sum(weights) ;
zhat = zhat + log(sum_weights) ;
% normalize weights
weights = weights/sum_weights ;
% (b) Selection
ESSsim(i) = 1/sum(weights.^2) ;
if (ESSsim(i) < options_.posterior_sampler_options.HSsmc.nparticles/2)
indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.HSsmc.nparticles) ;
param = param(:,indx_resmpl) ;
loglik = loglik(indx_resmpl) ;
tlogpost_i = tlogpost_i(indx_resmpl) ;
weights = ones(options_.posterior_sampler_options.HSsmc.nparticles,1)/options_.posterior_sampler_options.HSsmc.nparticles ;
end
% (c) Mutation
options_.posterior_sampler_options.HSsmc.c = options_.posterior_sampler_options.HSsmc.c*modified_logit(0.95,0.1,16.0,options_.posterior_sampler_options.HSsmc.acpt-options_.posterior_sampler_options.HSsmc.trgt) ;
% Calculate estimates of mean and variance
mu = param*weights ;
z = bsxfun(@minus,param,mu) ;
R = z*(bsxfun(@times,z',weights)) ;
Rchol = chol(R)' ;
% Mutation
if options_.posterior_sampler_options.HSsmc.option_mutation==1
[param,tlogpost_i,loglik,options_.posterior_sampler_options.HSsmc.acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,options_.posterior_sampler_options.HSsmc.c*Rchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
elseif options_.posterior_sampler_options.HSsmc.option_mutation==2
inv_R = inv(options_.posterior_sampler_options.HSsmc.c^2*R) ;
Rdiagchol = sqrt(diag(R)) ;
[param,tlogpost_i,loglik,options_.posterior_sampler_options.HSsmc.acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,options_.posterior_sampler_options.HSsmc.c*Rchol,options_.posterior_sampler_options.HSsmc.c*Rdiagchol,inv_R,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_) ;
end
acptsim(i) = options_.posterior_sampler_options.HSsmc.acpt ;
csim(i) = options_.posterior_sampler_options.HSsmc.c ;
% print information
fprintf(' Iteration = %5.0f / %5.0f \n', i, options_.posterior_sampler_options.HSsmc.nphi);
fprintf(' phi = %5.4f \n', phi(i));
fprintf(' Neff = %5.4f \n', ESSsim(i));
fprintf(' %accept. = %5.4f \n', acptsim(i));
end
indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.HSsmc.nparticles);
distrib_param = param(:,indx_resmpl);
fprintf(' Log_lik = %5.4f \n', zhat);
mean_xparam = mean(distrib_param,2);
npar = length(xparam1) ;
%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
%mat_var_cov = (mat_var_cov*mat_var_cov')/(options_.posterior_sampler_options.HSsmc.nparticles-1) ;
%std_xparam = sqrt(diag(mat_var_cov)) ;
lb95_xparam = zeros(npar,1) ;
ub95_xparam = zeros(npar,1) ;
for i=1:npar
temp = sortrows(distrib_param(i,:)') ;
lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.HSsmc.nparticles) ;
ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.HSsmc.nparticles) ;
end
TeX = options_.TeX;
str = sprintf(' Param. \t Lower Bound (95%%) \t Mean \t Upper Bound (95%%)');
for l=1:npar
[name,~] = get_the_name(l,TeX,M_,estim_params_,options_.varobs);
str = sprintf('%s\n %s \t\t %5.4f \t\t %7.5f \t\t %5.4f', str, name, lb95_xparam(l), mean_xparam(l), ub95_xparam(l));
end
disp([str])
disp('')
%% Plot parameters densities
[nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);
if TeX
fidTeX = fopen([M_.fname '_param_density.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by DSMH.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourier Transform approximation.
plt = 1 ;
%for plt = 1:nbplt,
if TeX
NAMES = [];
TeXNAMES = [];
end
hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
for k=1:npar %min(nstar,npar-(plt-1)*nstar)
subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
%kk = (plt-1)*nstar+k;
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.HSsmc.nparticles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
options_.posterior_sampler_options.HSsmc.nparticles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2));
hold on
if TeX
title(texname,'interpreter','latex')
else
title(name,'interpreter','none')
end
hold off
axis tight
drawnow
end
dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/floor(sqrt(npar)),1),M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameter densities based on the Herbst/Schorfheide sampler.}');
fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
%end
function [param,tlogpost_i,loglik,acpt] = mutation_RW(TargetFun,param,tlogpost_i,loglik,phi,i,cRchol,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
acpt = 0.0 ;
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
for j=1:options_.posterior_sampler_options.HSsmc.nparticles
validate= 0;
while validate==0
candidate = param(:,j) + cRchol*randn(size(param,1),1) ;
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
validate = 1;
if rand(1,1)<exp(tlogpostx-tlogpost_i(j)) % accept
acpt = acpt + 1 ;
param(:,j) = candidate;
loglik(j) = loglikx;
tlogpost_i(j) = tlogpostx;
end
end
end
end
end
acpt = acpt/options_.posterior_sampler_options.HSsmc.nparticles;
function [param,tlogpost_i,loglik,acpt] = mutation_Mixture(TargetFun,param,tlogpost_i,loglik,phi,i,cRchol,cRdiagchol,invR,mu,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
acpt = 0.0 ;
tlogpost_i = tlogpost_i + (phi(i)-phi(i-1))*loglik ;
for j=1:options_.posterior_sampler_options.HSsmc.nparticles
qx = 0 ;
q0 = 0 ;
alpind = rand(1,1) ;
validate= 0;
while validate==0
if alpind<options_.posterior_sampler_options.HSsmc.alp % RW, no need to modify qx and q0
candidate = param(:,j) + cRchol*randn(size(param,1),1);
elseif alpind<options_.posterior_sampler_options.HSsmc.alp + (1-options_.posterior_sampler_options.HSsmc.alp/2) % random walk with diagonal cov no need to modify qx and q0
candidate = param(:,j) + cRdiagchol*randn(size(param,1),1);
else % Proposal densities
candidate = mu + cRchol*randn(size(param,1),1);
qx = -.5*(candidate-mu)'*invR*(candidate-mu) ; % no need of the constants in the acceptation rule
q0 = -.5*(param(:,j)-mu)'*invR*(param(:,j)-mu) ;
end
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
[tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,phi(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
validate = 1;
if rand(1,1)<exp((tlogpostx-qx)-(tlogpost_i(j)-q0)) % accept
acpt = acpt + 1 ;
param(:,j) = candidate;
loglik(j) = loglikx;
tlogpost_i(j) = tlogpostx;
end
end
end
end
end
acpt = acpt/options_.posterior_sampler_options.HSsmc.nparticles;
function x = modified_logit(a,b,c,d)
tmp = exp(c*d) ;
x = a + b*tmp/(1 + tmp) ;

View File

@ -31,6 +31,7 @@ function oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
latexDirectoryName = CheckPath('latex',M_.dname);
graphDirectoryName = CheckPath('graphs',M_.dname);
@ -72,7 +73,7 @@ for i=1:npar
f1 = oo_.posterior_density.shocks_std.(name)(:,2);
oo_.prior_density.shocks_std.(name)(:,1) = x2;
oo_.prior_density.shocks_std.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.shocks_std.(name);
end
elseif i <= nvx+nvn
@ -81,7 +82,7 @@ for i=1:npar
f1 = oo_.posterior_density.measurement_errors_std.(name)(:,2);
oo_.prior_density.measurement_errors_std.(name)(:,1) = x2;
oo_.prior_density.measurement_errors_std.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.measurement_errors_std.(name);
end
elseif i <= nvx+nvn+ncx
@ -93,7 +94,7 @@ for i=1:npar
f1 = oo_.posterior_density.shocks_corr.(name)(:,2);
oo_.prior_density.shocks_corr.(name)(:,1) = x2;
oo_.prior_density.shocks_corr.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.shocks_corr.(name);
end
elseif i <= nvx+nvn+ncx+ncn
@ -105,7 +106,7 @@ for i=1:npar
f1 = oo_.posterior_density.measurement_errors_corr.(name)(:,2);
oo_.prior_density.measurement_errors_corr.(name)(:,1) = x2;
oo_.prior_density.measurement_errors_corr.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.measurement_errors_corr.(name);
end
else
@ -115,7 +116,7 @@ for i=1:npar
f1 = oo_.posterior_density.parameters.(name)(:,2);
oo_.prior_density.parameters.(name)(:,1) = x2;
oo_.prior_density.parameters.(name)(:,2) = f2;
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
pmod = oo_.posterior_mode.parameters.(name);
end
end
@ -130,7 +131,7 @@ for i=1:npar
set(hh_plt, 'color', [0.7 0.7 0.7]);
hold on;
plot(x1, f1, '-k', 'linewidth', 2);
if ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && ~options_.mh_posterior_mode_estimation
plot([pmod pmod], [0.0 1.1*top0], '--g', 'linewidth', 2);
end
box on
@ -160,4 +161,4 @@ for i=1:npar
end
subplotnum = 0;
end
end
end

View File

@ -1,113 +0,0 @@
function [ ix2, temperedlogpost, loglik, bayestopt_] = ...
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_, NumberOfParticles)
% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfParticles, bayestopt_] = ...
% SMC_samplers_initialization(TargetFun, xparam1, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,NumberOfParticles)
% Draw in prior distribution to initialize samplers.
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
% function (tempered posterior kernel and likelihood).
% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values).
% 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
%
% OUTPUTS
% o ix2 [double] (NumberOfParticles*npar) vector of starting points for different chains
% o ilogpo2 [double] (NumberOfParticles*1) vector of initial posterior values for different chains
% o iloglik2 [double] (NumberOfParticles*1) vector of initial likelihood values for different chains
% o ModelName [string] name of the mod-file
% o MetropolisFolder [string] path to the Metropolis subfolder
% o bayestopt_ [structure] estimation options structure
%
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2006-2022 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 outputs
ix2 = [];
ilogpo2 = [];
iloglik2 = [];
ModelName = [];
MetropolisFolder = [];
ModelName = M_.fname;
if ~isempty(M_.bvar)
ModelName = [ModelName '_bvar'];
end
MetropolisFolder = CheckPath('dsmh',M_.dname);
BaseName = [MetropolisFolder filesep ModelName];
npar = length(xparam1);
% Here we start a new DS Metropolis-Hastings, previous draws are discarded.
disp('Estimation:: Initialization...')
% Delete old dsmh files if any...
files = dir([BaseName '_dsmh*_blck*.mat']);
%if length(files)
% delete([BaseName '_dsmh*_blck*.mat']);
% disp('Estimation::smc: Old dsmh-files successfully erased!')
%end
% Delete old log file.
file = dir([ MetropolisFolder '/dsmh.log']);
%if length(file)
% delete([ MetropolisFolder '/dsmh.log']);
% disp('Estimation::dsmh: Old dsmh.log file successfully erased!')
% disp('Estimation::dsmh: Creation of a new dsmh.log file.')
%end
fidlog = fopen([MetropolisFolder '/dsmh.log'],'w');
fprintf(fidlog,'%% DSMH log file (Dynare).\n');
fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
fprintf(fidlog,' \n\n');
fprintf(fidlog,'%% Session 1.\n');
fprintf(fidlog,' \n');
prior_draw(bayestopt_,options_.prior_trunc);
% Find initial values for the NumberOfParticles chains...
options_=set_dynare_seed_local_options(options_,'default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::dsmh: Searching for initial values...');
ix2 = zeros(npar,NumberOfParticles);
temperedlogpost = zeros(NumberOfParticles,1);
loglik = zeros(NumberOfParticles,1);
%stderr = sqrt(bsxfun(@power,mh_bounds.ub-mh_bounds.lb,2)/12)/10;
for j=1:NumberOfParticles
validate = 0;
while validate == 0
candidate = prior_draw()';
% candidate = xparam1(:) + 0.001*randn(npar,1);%bsxfun(@times,stderr,randn(npar,1)) ;
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2(:,j) = candidate ;
[temperedlogpost(j),loglik(j)] = tempered_likelihood(TargetFun,candidate,0.0,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
if isfinite(loglik(j)) % if returned log-density is Inf or Nan (penalized value)
validate = 1;
end
end
end
end
fprintf(fidlog,' \n');
disp('Estimation:: Initial values found!')
skipline()

View File

@ -1,20 +1,17 @@
function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_,outputFolderName)
% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_,outputFolderName)
% initialization of posterior samplers
function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_, outputFolderName)
% Initialization of posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% fname: name of the mod-file
% dname: name of directory with metropolis folder
% options_: structure storing the options
% bounds: structure containing prior bounds
% bayestopt_: structure storing information about priors
% - posterior_sampler_options [struct] posterior sampler options
% - options_ [struct] options
% - bounds [struct] prior bounds
% - bayestopt_ [struct] information about priors
%
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
% options_: structure storing the options
% bayestopt_: structure storing information about priors
% outputFolderName: string of folder to store mat files
% - posterior_sampler_options [struct] checked posterior sampler options (updated)
% - options_ [struct] options (updated)
% - bayestopt_ [struct] information about priors (updated)
%
% SPECIAL REQUIREMENTS
% none
@ -40,15 +37,17 @@ if nargin < 7
outputFolderName = 'Output';
end
init=0;
init = false;
if isempty(posterior_sampler_options)
init=1;
init = true;
end
if init
% set default options and user defined options
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method;
posterior_sampler_options.bounds = bounds;
if ~issmc(options_)
posterior_sampler_options.bounds = bounds;
end
switch posterior_sampler_options.posterior_sampling_method
@ -227,7 +226,6 @@ if init
end
end
case 'slice'
posterior_sampler_options.parallel_bar_refresh_rate=1;
posterior_sampler_options.serial_bar_refresh_rate=1;
@ -342,50 +340,120 @@ if init
end
% moreover slice must be associated to:
% options_.mh_posterior_mode_estimation = false;
% this is done below, but perhaps preprocessing should do this?
% options_.mh_posterior_mode_estimation = false;
% this is done below, but perhaps preprocessing should do this?
if ~isempty(posterior_sampler_options.mode)
% multimodal case
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
% posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
if ~isempty(posterior_sampler_options.mode)
% multimodal case
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
% posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
if options_.load_mh_file
posterior_sampler_options.slice_initialize_with_mode = 0;
else
if ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
if options_.load_mh_file
posterior_sampler_options.slice_initialize_with_mode = 0;
else
if ~posterior_sampler_options.slice_initialize_with_mode
posterior_sampler_options.invhess=[];
end
end
if ~isempty(posterior_sampler_options.mode_files) % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
load(modes, 'xparams')
if size(xparams,2)<2
error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
end
for j=1:size(xparams,2)
mode(j).m=xparams(:,j);
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
case 'hssmc'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
% 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 'target'
posterior_sampler_options.target = options_list{i,2};
case 'steps'
posterior_sampler_options.steps = options_list{i,2};
case 'scale'
posterior_sampler_options.scale = options_list{i,2};
case 'particles'
posterior_sampler_options.particles = options_list{i,2};
case 'lambda'
posterior_sampler_options.lambda = options_list{i,2};
otherwise
warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
case 'dsmh'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
% 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 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
case 'number_of_particles'
posterior_sampler_options.particles = options_list{i,2};
otherwise
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
end
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
load(modes, 'xparams')
if size(xparams,2)<2
error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
end
for j=1:size(xparams,2)
mode(j).m=xparams(:,j);
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = true;
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
return
return
end
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
if ~options_.cova_compute && ~(options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix)
if strcmp('hessian',options_.MCMC_jumping_covariance)
skipline()
disp('check_posterior_sampler_options:: I cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed')

View File

@ -1,23 +1,19 @@
function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname,outputFolderName)
% function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname,outputFolderName)
function [mean, covariance, mode, kernel_at_the_mode] = compute_mh_covariance_matrix(names, fname, dname, outputFolderName)
% Estimation of the posterior covariance matrix, posterior mean, posterior mode and evaluation of the posterior kernel at the
% estimated mode, using the draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in
% a file <fname>_mh_mode.mat.
% estimated mode, using posterior draws from a metropolis-hastings.
%
% INPUTS
% o bayestopt_ [struct] characterizing priors
% o fname [string] name of model
% o dname [string] name of directory with metropolis folder
% o outputFolderName [string] name of directory to store results
% - names [cell] n×1 cell array of row char arrays, names of the estimated parameters.
% - fname [char] name of the model
% - dname [char] name of subfolder with output files
% - outputFolderName [char] name of directory to store results
%
% OUTPUTS
% o posterior_mean [double] (n*1) vector, posterior expectation of the parameters.
% o posterior_covariance [double] (n*n) matrix, posterior covariance of the parameters (computed from previous metropolis hastings).
% o posterior_mode [double] (n*1) vector, posterior mode of the parameters.
% o posterior_kernel_at_the_mode [double] scalar.
%
% SPECIAL REQUIREMENTS
% None.
% - mean [double] n×1 vector, posterior expectation of the parameters.
% - covariance [double] n×n matrix, posterior covariance of the parameters.
% - mode [double] n×1 vector, posterior mode of the parameters.
% - kernel_at_the_mode [double] scalar, value of the posterior kernel at the mode.
% Copyright © 2006-2023 Dynare Team
%
@ -35,9 +31,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin < 4
outputFolderName = 'Output';
end
MetropolisFolder = CheckPath('metropolis',dname);
BaseName = [MetropolisFolder filesep fname];
@ -49,29 +43,24 @@ TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
[nblck, n] = size(record.LastParameters);
posterior_kernel_at_the_mode = -Inf;
posterior_mean = zeros(n,1);
posterior_mode = NaN(n,1);
posterior_covariance = zeros(n,n);
kernel_at_the_mode = -Inf;
mean = zeros(n,1);
mode = NaN(n,1);
covariance = zeros(n,n);
offset = 0;
for b=1:nblck
first_line = FirstLine;
for n = FirstMhFile:TotalNumberOfMhFiles
load([ BaseName '_mh' int2str(n) '_blck' int2str(b) '.mat'],'x2','logpo2');
[tmp,idx] = max(logpo2);
if tmp>posterior_kernel_at_the_mode
posterior_kernel_at_the_mode = tmp;
posterior_mode = x2(idx,:);
[tmp, idx] = max(logpo2);
if tmp>kernel_at_the_mode
kernel_at_the_mode = tmp;
mode = x2(idx,:);
end
[posterior_mean,posterior_covariance,offset] = recursive_moments(posterior_mean,posterior_covariance,x2(first_line:end,:),offset);
[mean, covariance, offset] = recursive_moments(mean, covariance, x2(first_line:end,:), offset);
first_line = 1;
end
end
xparam1 = posterior_mode';
hh = inv(posterior_covariance);
fval = posterior_kernel_at_the_mode;
parameter_names = bayestopt_.name;
save([dname filesep outputFolderName filesep fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names');
mode = transpose(mode);

View File

@ -0,0 +1,65 @@
function [mu, covariance, mode, kernel_at_the_mode] = compute_posterior_covariance_matrix(names, fname, dname, options_, outputFolderName)
% Estimation of the posterior covariance matrix, posterior mean, posterior mode and evaluation of the posterior kernel at the
% estimated mode, using posterior draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in
% a file <fname>_mh_mode.mat, hssmc_mode.mat or dsmh__mode.mat under <dname>/<outputFolderName>/.
%
% INPUTS
% - names [cell] n×1 cell array of row char arrays, names of the estimated parameters.
% - fname [char] name of the model
% - dname [char] name of subfolder with output files
% - outputFolderName [char] name of directory to store results
%
% OUTPUTS
% - mean [double] n×1 vector, posterior expectation of the parameters.
% - covariance [double] n×n matrix, posterior covariance of the parameters.
% - mode [double] n×1 vector, posterior mode of the parameters.
% - kernel_at_the_mode [double] scalar, value of the posterior kernel at the mode.
% Copyright © 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/>.
if nargin<5
outputFolderName = 'Output';
end
if ishssmc(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
% Get the posterior mode
[kernel_at_the_mode, id] = max(posterior.tlogpostkernel);
mode = posterior.particles(:,id);
% Compute the posterior mean
mu = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
covariance = (posterior.particles-mu)*(posterior.particles-mu)'/length(posterior.tlogpostkernel);
else
[mu, covariance, mode, kernel_at_the_mode] = compute_mh_covariance_matrix(names, fname, dname, outputFolderName);
end
xparam1 = mode;
hh = inv(covariance);
fval = kernel_at_the_mode;
parameter_names = names;
if ishssmc(options_)
save(sprintf('%s/%s/hssmc_mode.mat', dname, outputFolderName), 'xparam1', 'hh', 'fval', 'parameter_names');
else
save(sprintf('%s/%s/%s_mh_mode.mat', dname, outputFolderName, fname), 'xparam1', 'hh', 'fval', 'parameter_names');
end

View File

@ -79,7 +79,7 @@ for jj = 1:npar
par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_.varobs);
param_name = vertcat(param_name, par_name_temp);
end
Draws = GetAllPosteriorDraws(M_.dname, M_.fname, jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
Draws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
Draws = reshape(Draws, [NumberOfDraws nblck]);
Nc = min(1000, NumberOfDraws/2);
for ll = 1:nblck

View File

@ -451,6 +451,7 @@ options_.nk = 1;
options_.noconstant = false;
options_.nodiagnostic = false;
options_.mh_posterior_mode_estimation = false;
options_.smc_posterior_mode_estimation = false;
options_.prefilter = 0;
options_.presample = 0;
options_.prior_trunc = 1e-10;
@ -463,7 +464,7 @@ options_.use_mh_covariance_matrix = false;
options_.gradient_method = 2; %used by csminwel and newrat
options_.gradient_epsilon = 1e-6; %used by csminwel and newrat
options_.posterior_sampler_options.sampling_opt = []; %extended set of options for individual posterior samplers
% Random Walk Metropolis-Hastings
% Random Walk Metropolis-Hastings
options_.posterior_sampler_options.posterior_sampling_method = 'random_walk_metropolis_hastings';
options_.posterior_sampler_options.rwmh.proposal_distribution = 'rand_multivariate_normal';
options_.posterior_sampler_options.rwmh.student_degrees_of_freedom = 3;
@ -491,23 +492,19 @@ options_.posterior_sampler_options.imh.proposal_distribution = 'rand_multivariat
options_.posterior_sampler_options.imh.use_mh_covariance_matrix=0;
options_.posterior_sampler_options.imh.save_tmp_file=0;
% Herbst and Schorfeide SMC Sampler
%options_.posterior_sampler = 'Herbst_Schorfheide' ;
options_.posterior_sampler_options.HSsmc.nphi= 25 ;
options_.posterior_sampler_options.HSsmc.lambda = 2 ;
options_.posterior_sampler_options.HSsmc.nparticles = 20000 ;
options_.posterior_sampler_options.HSsmc.c = 0.5 ;
options_.posterior_sampler_options.HSsmc.acpt = 0.25 ;
options_.posterior_sampler_options.HSsmc.trgt = 0.25 ;
options_.posterior_sampler_options.HSsmc.option_mutation = 1 ;
options_.posterior_sampler_options.HSsmc.alp = .9 ;
options_.posterior_sampler_options.hssmc.steps= 25;
options_.posterior_sampler_options.hssmc.lambda = 2;
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;
% DSMH: Dynamic Striated Metropolis-Hastings algorithm
%options_.posterior_sampler = 'DSMH' ;
options_.posterior_sampler_options.dsmh.H = 25 ;
options_.posterior_sampler_options.dsmh.N = 20 ;
options_.posterior_sampler_options.dsmh.G = 10 ;
options_.posterior_sampler_options.dsmh.K = 50 ;
options_.posterior_sampler_options.dsmh.lambda1 = 0.1 ;
options_.posterior_sampler_options.dsmh.nparticles = 20000 ;
options_.posterior_sampler_options.dsmh.particles = 20000 ;
options_.posterior_sampler_options.dsmh.alpha0 = 0.2 ;
options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
options_.posterior_sampler_options.dsmh.tau = 10 ;
@ -664,12 +661,6 @@ options_.use_dll = false;
% model evaluated using bytecode.dll
options_.bytecode = false;
% if true, use a fixed point method to solve Sylvester equation (for large scale models)
options_.sylvester_fp = false;
% convergence criteria to solve iteratively a sylvester equations
options_.sylvester_fixed_point_tol = 1e-12;
% if true, use a fixed point method to solve Lyapunov equation (for large scale models)
options_.lyapunov_fp = false;
% if true, use a doubling algorithm to solve Lyapunov equation (for large scale models)

View File

@ -17,4 +17,4 @@ function dprintf(str, varargin)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
disp(sprintf(str, varargin{:}));
disp(sprintf(str, varargin{:}));

View File

@ -11,7 +11,7 @@ function forecast = dyn_forecast(var_list,M_,options_,oo_,task,dataset_info)
% task: indicates how to initialize the forecast
% either 'simul' or 'smoother'
% dataset_info: Various informations about the dataset (descriptive statistics and missing observations).
%
% OUTPUTS
% forecast: structure containing fields
% Mean: point estimate

View File

@ -51,6 +51,9 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/discretionary_policy/' ; ...
'/distributions/' ; ...
'/ep/' ; ...
'/estimation/'; ...
'/estimation/smc/'; ...
'/estimation/resampler/'; ...
'/gsa/' ; ...
'/kalman/' ; ...
'/kalman/likelihood' ; ...

View File

@ -31,6 +31,14 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
if issmc(options_)
options_.mode_compute = 0;
options_.mh_replic = 0;
options_.mh_recover = false;
options_.load_mh_file = false;
options_.load_results_after_load_mh = false;
end
dispString = 'Estimation::mcmc';
if ~exist([M_.dname filesep 'Output'],'dir')
@ -88,7 +96,9 @@ if options_.order > 1
end
end
%% set objective function
%
% set objective function
%
if ~options_.dsge_var
if options_.particle.status
objective_function = str2func('non_linear_dsge_likelihood');
@ -147,7 +157,10 @@ if ~isempty(estim_params_)
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
%% perform initial estimation checks;
%
% perform initial estimation checks;
%
try
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);
catch % if check fails, provide info on using calibration if present
@ -164,8 +177,11 @@ catch % if check fails, provide info on using calibration if present
rethrow(e);
end
%% Run smoother if no estimation or mode-finding are requested
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
%
% Run smoother if no estimation or mode-finding are requested
%
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation && ~issmc(options_)
if options_.order==1 && ~options_.particle.status
if options_.smoother
if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
@ -210,11 +226,11 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
end
end
%% Estimation of the posterior mode or likelihood mode
%
% Estimation of the posterior mode or likelihood mode
%
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation && ~issmc(options_)
optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
@ -238,7 +254,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 2;
[~,~,~,~,hh] = feval(objective_function,xparam1, ...
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.analytic_derivation = ana_deriv_old;
elseif ~isnumeric(current_optimizer) || ~(isequal(current_optimizer,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood'))
% enter here if i) not mode_compute_5, ii) if mode_compute_5 and newratflag==1;
@ -292,16 +308,19 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
end
if ~options_.mh_posterior_mode_estimation && options_.cova_compute
if ~options_.mh_posterior_mode_estimation && options_.cova_compute && ~issmc(options_)
check_hessian_at_the_mode(hh, xparam1, M_, estim_params_, options_, bounds);
end
%% create mode_check_plots
if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
%
% create mode_check_plots
%
if options_.mode_check.status && ~options_.mh_posterior_mode_estimation && ~issmc(options_)
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
mode_check(objective_function,xparam1,hh,options_,M_,estim_params_,bayestopt_,bounds,false,...
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
options_.analytic_derivation = ana_deriv_old;
end
@ -309,43 +328,38 @@ oo_.posterior.optimization.mode = [];
oo_.posterior.optimization.Variance = [];
oo_.posterior.optimization.log_density=[];
invhess=[];
if ~options_.mh_posterior_mode_estimation
oo_.posterior.optimization.mode = xparam1;
if exist('fval','var')
oo_.posterior.optimization.log_density=-fval;
invhess = [];
if ~issmc(options_)
if ~options_.mh_posterior_mode_estimation
oo_.posterior.optimization.mode = xparam1;
if exist('fval','var')
oo_.posterior.optimization.log_density=-fval;
end
if options_.cova_compute
hsd = sqrt(diag(hh));
invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
stdh = sqrt(diag(invhess));
oo_.posterior.optimization.Variance = invhess;
end
else
variances = bayestopt_.p2.*bayestopt_.p2;
idInf = isinf(variances);
variances(idInf) = 1;
invhess = options_.mh_posterior_mode_estimation*diag(variances);
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
hsd = sqrt(diag(hh));
invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
stdh = sqrt(diag(invhess));
oo_.posterior.optimization.Variance = invhess;
end
else
variances = bayestopt_.p2.*bayestopt_.p2;
idInf = isinf(variances);
variances(idInf) = 1;
invhess = options_.mh_posterior_mode_estimation*diag(variances);
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
stdh = NaN(length(xparam1),1);
end
if options_.particle.status && isfield(options_.particle,'posterior_sampler')
if strcmpi(options_.particle.posterior_sampler,'Herbst_Schorfheide')
Herbst_Schorfheide_sampler(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 strcmpi(options_.particle.posterior_sampler,'DSMH')
DSMH_sampler(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)
end
end
if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
if ~issmc(options_) && any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
% display results table and store parameter estimates and standard errors in results
oo_ = display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Posterior', 'posterior');
% Laplace approximation to the marginal log density:
@ -366,56 +380,75 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
[~,~,~,~,~,~,~,oo_.dsge_var.posterior_mode.PHI_tilde,oo_.dsge_var.posterior_mode.SIGMA_u_tilde,oo_.dsge_var.posterior_mode.iXX,oo_.dsge_var.posterior_mode.prior] =...
feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
elseif ~issmc(options_) && ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
oo_=display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Maximum Likelihood', 'mle');
end
invhess = set_mcmc_jumping_covariance(invhess, nx, options_.MCMC_jumping_covariance, bayestopt_, 'dynare_estimation_1');
if ~issmc(options_)
invhess = set_mcmc_jumping_covariance(invhess, nx, options_.MCMC_jumping_covariance, bayestopt_, 'dynare_estimation_1');
end
if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
(any(bayestopt_.pshape >0 ) && options_.load_mh_file) %% not ML estimation
%reset bounds as lb and ub must only be operational during mode-finding
bounds = set_mcmc_prior_bounds(xparam1, bayestopt_, options_, 'dynare_estimation_1');
% Tunes the jumping distribution's scale parameter
if options_.mh_tune_jscale.status
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
options_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds,...
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
bayestopt_.jscale(:) = options_.mh_jscale;
fprintf('mh_jscale has been set equal to %s\n', num2str(options_.mh_jscale));
else
warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!')
%
% Run SMC sampler.
%
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 isdsmh(options_)
dsmh(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
end
%
% Run MCMC and compute posterior statistics.
%
if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(bayestopt_.pshape>0) && options_.load_mh_file) % not ML estimation
if ~issmc(options_)
% Reset bounds as lb and ub must only be operational during mode-finding
bounds = set_mcmc_prior_bounds(xparam1, bayestopt_, options_, 'dynare_estimation_1');
% Tune the jumping distribution's scale parameter
if options_.mh_tune_jscale.status
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
options_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds,...
dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
bayestopt_.jscale(:) = options_.mh_jscale;
fprintf('mh_jscale has been set equal to %s\n', num2str(options_.mh_jscale));
else
warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!')
end
end
end
% runs MCMC
if options_.mh_replic || options_.load_mh_file
posterior_sampler_options = options_.posterior_sampler_options.current_options;
posterior_sampler_options.invhess = invhess;
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, M_.fname, M_.dname, options_, bounds, bayestopt_);
% store current options in global
options_.posterior_sampler_options.current_options = posterior_sampler_options;
if options_.mh_replic
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString);
options_.analytic_derivation = ana_deriv_old;
% Run MCMC
if options_.mh_replic || options_.load_mh_file
posterior_sampler_options = options_.posterior_sampler_options.current_options;
posterior_sampler_options.invhess = invhess;
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, M_.fname, M_.dname, options_, bounds, bayestopt_);
% store current options in global
options_.posterior_sampler_options.current_options = posterior_sampler_options;
if options_.mh_replic
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 0;
posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString);
options_.analytic_derivation = ana_deriv_old;
end
end
% Discard first mh_drop percent of the draws:
CutSample(M_, options_, dispString);
end
%% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, dispString);
if options_.mh_posterior_mode_estimation
[~,~,posterior_mode,~] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname);
oo_=fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
if options_.mh_posterior_mode_estimation || (issmc(options_) && options_.smc_posterior_mode_estimation)
[~, covariance, posterior_mode, ~] = compute_posterior_covariance_matrix(bayestopt_.name, M_.fname, M_.dname, options_);
oo_ = fill_mh_mode(posterior_mode, sqrt(diag(covariance)), M_, options_, estim_params_, oo_, 'posterior');
%reset qz_criterium
options_.qz_criterium=qz_criterium_old;
options_.qz_criterium = qz_criterium_old;
return
else
%get stored results if required
if options_.load_mh_file && options_.load_results_after_load_mh
oo_load_mh=load([M_.dname filesep 'Output' filesep M_.fname '_results'],'oo_');
% Get stored results if required
if ~issmc(options_) && options_.load_mh_file && options_.load_results_after_load_mh
oo_load_mh = load(sprintf('%s/%s/%s_results', M_.dname, 'Output', M_.fname), 'oo_');
end
if ~options_.nodiagnostic
% Compute MCMC convergence diagnostics
if ~issmc(options_) && ~options_.nodiagnostic
if (options_.mh_replic>0 || (options_.load_mh_file && ~options_.load_results_after_load_mh))
oo_= mcmc_diagnostics(options_, estim_params_, M_,oo_);
elseif options_.load_mh_file && options_.load_results_after_load_mh
@ -424,9 +457,11 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end
end
end
%% Estimation of the marginal density from the Mh draws:
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
[~,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
% 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 ~issmc(options_)
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
% Store posterior statistics by parameter name
oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, prior_dist_names);
if ~options_.nograph
@ -435,13 +470,13 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
% Store posterior mean in a vector and posterior variance in
% a matrix
[oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
= GetPosteriorMeanVariance(M_,options_.mh_drop);
= GetPosteriorMeanVariance(options_, M_);
elseif options_.load_mh_file && options_.load_results_after_load_mh
%% load fields from previous MCMC run stored in results-file
% load fields from previous MCMC run stored in results-file
field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
'prior_density',...%fields set by PlotPosteriorDistributions
};
'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
'prior_density',...%fields set by PlotPosteriorDistributions
};
for field_iter=1:size(field_names,2)
if isfield(oo_load_mh.oo_,field_names{1,field_iter})
oo_.(field_names{1,field_iter})=oo_load_mh.oo_.(field_names{1,field_iter});
@ -456,7 +491,9 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
oo_.posterior.metropolis=oo_load_mh.oo_.posterior.metropolis;
end
end
[error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_);
if ~issmc(options_)
[error_flag, ~, options_]= metropolis_draw(1, options_, estim_params_, M_);
end
if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
if options_.bayesian_irf
if error_flag
@ -499,9 +536,9 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
error('%s: Particle Smoothers are not yet implemented.',dispString)
end
end
else
fprintf('%s: sub_draws was set to 0. Skipping posterior computations.',dispString);
end
else
fprintf('%s: sub_draws was set to 0. Skipping posterior computations.',dispString);
end
xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
@ -517,7 +554,7 @@ end
%Run and store classical smoother if needed
if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && ~options_.partial_information % to be fixed
|| ~options_.smoother ) && ~options_.partial_information % to be fixed
%% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothed variables
oo_=save_display_classical_smoother_results(xparam1,M_,oo_,options_,bayestopt_,dataset_,dataset_info,estim_params_);
end

View File

@ -1,8 +1,6 @@
function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
% function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, bounds] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
% performs initialization tasks before estimation or
% global sensitivity analysis
% Performs initialization tasks before estimation or global sensitivity analysis
%
% INPUTS
% var_list_: selected endogenous variables vector
@ -390,7 +388,7 @@ end
%set options for old interface from the ones for new interface
if ~isempty(dataset_)
options_.nobs = dataset_.nobs;
if options_.endogenous_prior
if options_.endogenous_prior
if ~isnan(dataset_info.missing.number_of_observations) && ~(dataset_info.missing.number_of_observations==0) %missing observations present
if dataset_info.missing.no_more_missing_observations<dataset_.nobs-10
fprintf('\ndynare_estimation_init: There are missing observations in the data.\n')
@ -537,15 +535,15 @@ end
if options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter
if ~isempty(options_.nk)
fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead. Disabling the option.\n')
fprintf('dynare_estimation_init: the inversion filter does not support filter_step_ahead. Disabling the option.\n')
options_.nk=[];
end
if options_.filter_covariance
fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
fprintf('dynare_estimation_init: the inversion filter does not support filter_covariance. Disabling the option.\n')
options_.filter_covariance=false;
end
if options_.smoothed_state_uncertainty
fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
fprintf('dynare_estimation_init: the inversion filter does not support smoothed_state_uncertainty. Disabling the option.\n')
options_.smoothed_state_uncertainty=false;
end
end

View File

@ -1,6 +1,15 @@
function x = bseastr(s1,s2)
function indices = kitagawa(weights, noise)
% Copyright © 2001-2017 Dynare Team
% Return indices for resampling.
%
% INPUTS
% - weights [double] n×1 vector of partcles' weights.
% - noise [double] scalar, uniform random deviates in [0,1]
%
% OUTPUTS
% - indices [integer] n×1 vector of indices in [1:n]
% Copyright © 2022-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -17,31 +26,20 @@ function x = bseastr(s1,s2)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
m = size(s1,1) ;
x = zeros(m,1) ;
s1=upper(deblank(s1));
s2=upper(deblank(s2));
n= length(weights);
for im = 1:m
key = s1(im,:) ;
h = size(s2,1) ;
l = 1 ;
while l <= h
mid = round((h+l)/2) ;
temp = s2(mid,:) ;
if ~ strcmp(key,temp)
for i = 1:min(length(key),length(temp))
if temp(i) > key(i)
h = mid - 1 ;
break
else
l = mid + 1 ;
break
end
end
else
x(im) = mid ;
break
end
if nargin<2, noise = rand; end
indices = NaN(n, 1);
cweights = cumsum(weights);
wweights = (transpose(0:n-1)+noise)*(1.0/n);
j = 1;
for i=1:n
while wweights(i)>cweights(j)
j = j+1;
end
indices(i) = j;
end

View File

@ -1,5 +1,5 @@
function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% Dynamic Striated Metropolis-Hastings algorithm.
%
% INPUTS
@ -33,7 +33,7 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management functions.
% Copyright © 2006-2023 Dynare Team
% Copyright © 2022-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -50,14 +50,15 @@ function DSMH_sampler(TargetFun,xparam1,mh_bounds,dataset_,dataset_info,options_
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
opts = options_.posterior_sampler_options.dsmh;
lambda = exp(bsxfun(@minus,options_.posterior_sampler_options.dsmh.H,1:1:options_.posterior_sampler_options.dsmh.H)/(options_.posterior_sampler_options.dsmh.H-1)*log(options_.posterior_sampler_options.dsmh.lambda1));
c = 0.055 ;
MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/10) ;
% Step 0: Initialization of the sampler
[ param, tlogpost_iminus1, loglik, bayestopt_] = ...
SMC_samplers_initialization(TargetFun, xparam1, mh_bounds, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,options_.posterior_sampler_options.dsmh.nparticles);
[param, tlogpost_iminus1, loglik, bayestopt_] = ...
smc_samplers_initialization(TargetFun, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
zhat = 1 ;
@ -78,20 +79,17 @@ end
weights = exp(loglik*(lambda(end)-lambda(end-1)));
weights = weights/sum(weights);
indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.nparticles);
indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.particles);
distrib_param = param(:,indx_resmpl);
mean_xparam = mean(distrib_param,2);
npar = length(xparam1);
%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
%mat_var_cov = (mat_var_cov*mat_var_cov')/(options_.HSsmc.nparticles-1) ;
%std_xparam = sqrt(diag(mat_var_cov)) ;
lb95_xparam = zeros(npar,1) ;
ub95_xparam = zeros(npar,1) ;
for i=1:npar
temp = sortrows(distrib_param(i,:)') ;
lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.nparticles) ;
ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.nparticles) ;
lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.particles) ;
ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.particles) ;
end
TeX = options_.TeX;
@ -130,9 +128,9 @@ for k=1:npar %min(nstar,npar-(plt-1)*nstar)
subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
%kk = (plt-1)*nstar+k;
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.nparticles,bandwidth,kernel_function);
optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
options_.posterior_sampler_options.dsmh.nparticles,optimal_bandwidth,kernel_function);
options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2));
hold on
if TeX

View File

@ -0,0 +1,136 @@
function mdd = hssmc(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
% Sequential Monte-Carlo sampler, Herbst and Schorfheide (JAE, 2014).
%
% INPUTS
% - TargetFun [char] string specifying the name of the objective function (posterior kernel).
% - xparam1 [double] p×1 vector of parameters to be estimated (initial values).
% - 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] estimated parameters
% - 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/>.
smcopt = options_.posterior_sampler_options.current_options;
% Set location for the simulated particles.
SimulationFolder = CheckPath('hssmc', M_.dname);
% Define prior distribution
Prior = dprior(bayestopt_, options_.prior_trunc);
% Set function handle for the objective
eval(sprintf('%s = @(x) %s(x, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, mh_bounds, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, []);', 'funobj', func2str(TargetFun)));
mlogit = @(x) .95 + .1/(1 + exp(-16*x)); % Update of the scale parameter
% Create the tempering schedule
phi = ((0:smcopt.steps-1)/(smcopt.steps-1)).^smcopt.lambda;
% Initialise the estimate of the marginal density of the data
mdd = .0;
% tuning for MH algorithms matrices
scl = zeros(smcopt.steps, 1); % scale parameter
ESS = zeros(smcopt.steps, 1); % ESS
acpt = zeros(smcopt.steps, 1); % average acceptance rate
% Initialization of the sampler (draws from the prior distribution with finite logged likelihood)
t0 = tic;
[particles, tlogpostkernel, loglikelihood] = ...
smc_samplers_initialization(funobj, 'hssmc', smcopt.particles, Prior, SimulationFolder, smcopt.steps);
tt = toc(t0);
dprintf('#Iter. lambda ESS Acceptance rate scale resample seconds')
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', 1, 0, 0, 0, 0, 'no', tt)
weights = ones(smcopt.particles, 1)/smcopt.particles;
resampled_particle_swarm = false;
for i=2:smcopt.steps % Loop over the weight on the liklihood (phi)
weights = weights.*exp((phi(i)-phi(i-1))*loglikelihood);
sweight = sum(weights);
weights = weights/sweight;
mdd = mdd + log(sweight);
ESS(i) = 1.0/sum(weights.^2);
if (2*ESS(i) < smcopt.particles) % Resampling
resampled_particle_swarm = true;
iresample = kitagawa(weights);
particles = particles(:,iresample);
loglikelihood = loglikelihood(iresample);
tlogpostkernel = tlogpostkernel(iresample);
weights = ones(smcopt.particles, 1)/smcopt.particles;
end
smcopt.scale = smcopt.scale*mlogit(smcopt.acpt-smcopt.target); % Adjust the scale parameter
scl(i) = smcopt.scale; % Scale parameter (for the jumping distribution in MH mutation step).
mu = particles*weights; % Weighted average of the particles.
z = particles-mu;
R = z*(z'.*weights); % Weighted covariance matrix of the particles.
t0 = tic;
acpt_ = false(smcopt.particles, 1);
tlogpostkernel = tlogpostkernel + (phi(i)-phi(i-1))*loglikelihood;
[acpt_, particles, loglikelihood, tlogpostkernel] = ...
randomwalk(funobj, chol(R, 'lower'), mu, scl(i), phi(i), acpt_, Prior, particles, loglikelihood, tlogpostkernel);
smcopt.acpt = sum(acpt_)/smcopt.particles; % Acceptance rate.
tt = toc(t0);
acpt(i) = smcopt.acpt;
if resampled_particle_swarm
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'yes', tt)
else
dprintf('%3u %5.4f %9.5E %5.4f %4.3f %3s %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'no', tt)
end
if i==smcopt.steps
iresample = kitagawa(weights);
particles = particles(:,iresample);
end
save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, smcopt.steps), 'particles', 'tlogpostkernel', 'loglikelihood')
resampled_particle_swarm = false;
end
end
function [accept, particles, loglikelihood, tlogpostkernel] = randomwalk(funobj, RR, mu, scale, phi, accept, Prior, particles, loglikelihood, tlogpostkernel)
parfor j=1:size(particles, 2)
notvalid= true;
while notvalid
candidate = particles(:,j) + scale*(RR*randn(size(mu)));
if Prior.admissible(candidate)
[tlogpost, loglik] = tempered_likelihood(funobj, candidate, phi, Prior);
if isfinite(loglik)
notvalid = false;
if rand<exp(tlogpost-tlogpostkernel(j))
accept(j) = true ;
particles(:,j) = candidate;
loglikelihood(j) = loglik;
tlogpostkernel(j) = tlogpost;
end
end
end
end
end
end

View File

@ -1,7 +1,6 @@
function indx = smc_resampling(weights,noise,number)
% function indx = smc_resampling(weights,noise,number)
function bool = isdsmh(options_)
% Copyright © 2022 Dynare Team
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
@ -18,13 +17,9 @@ function indx = smc_resampling(weights,noise,number)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
indx = zeros(number,1);
cumweights = cumsum(weights);
randvec = (transpose(1:number)-1+noise(:))/number;
j = 1;
for i=1:number
while (randvec(i)>cumweights(j))
j = j+1;
end
indx(i) = j;
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'dsmh')
bool = true;
end
end

View File

@ -0,0 +1,25 @@
function bool = ishssmc(options_)
% Copyright © 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/>.
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'hssmc')
bool = true;
end
end

View File

@ -0,0 +1,81 @@
function [particles, tlogpostkernel, loglikelihood, SimulationFolder] = smc_samplers_initialization(funobj, sampler, n, Prior, SimulationFolder, nsteps)
% Initialize SMC samplers by drawing initial particles in the prior distribution.
%
% INPUTS
% - TargetFun [char] string specifying the name of the objective function (posterior kernel).
% - sampler [char] name of the sampler.
% - n [integer] scalar, number of particles.
% - mh_bounds [double] p×2 matrix defining lower and upper bounds for the estimated 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] estimated parameters
% - oo_ [struct] outputs
%
% OUTPUTS
% - ix2 [double] p×n matrix of particles
% - ilogpo2 [double] n×1 vector of posterior kernel values for the particles
% - iloglik2 [double] n×1 vector of likelihood values for the particles
% - ModelName [string] name of the mod-file
% - MetropolisFolder [string] path to the Metropolis subfolder
% - bayestopt_ [structure] estimation options structure
%
% 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/>.
dprintf('Estimation:%s: Initialization...', sampler)
% Delete old mat files storign particles if any...
matfiles = sprintf('%s%sparticles*.mat', SimulationFolder, filesep());
files = dir(matfiles);
if ~isempty(files)
delete(matfiles);
dprintf('Estimation:%s: Old %s-files successfully erased.', sampler, sampler)
end
% Simulate a pool of particles characterizing the prior distribution (with the additional constraint that the likelihood is finite)
set_dynare_seed('default');
dprintf('Estimation:%s: Searching for initial values...', sampler);
particles = zeros(Prior.length(), n);
tlogpostkernel = zeros(n, 1);
loglikelihood = zeros(n, 1);
t0 = tic;
parfor j=1:n
notvalid = true;
while notvalid
candidate = Prior.draw();
if Prior.admissible(candidate)
particles(:,j) = candidate;
[tlogpostkernel(j), loglikelihood(j)] = tempered_likelihood(funobj, candidate, 0.0, Prior);
if isfinite(loglikelihood(j)) % if returned log-density is Inf or Nan (penalized value)
notvalid = false;
end
end
end
end
tt = toc(t0);
save(sprintf('%s%sparticles-1-%u.mat', SimulationFolder, filesep(), nsteps), 'particles', 'tlogpostkernel', 'loglikelihood')
dprintf('Estimation:%s: Initial values found (%.2fs)', sampler, tt)
skipline()

View File

@ -0,0 +1,35 @@
function [tlogpostkernel,loglikelihood] = tempered_likelihood(postkernelfun, xparam, lambda, Prior)
% Evaluate tempered likelihood (posterior kernel)
%
% INPUTS
% - postkernelfun [handle] Function handle for the opposite of the posterior kernel.
% - xparam [double] n×1 vector of parameters.
% - lambda [double] scalar between 0 and 1, weight on the posterior kernel.
% - Prior [dprior] Prior specification.
%
% OUTPUTS
% - tlogpostkernel [double] scalar, value of the tempered posterior kernel.
% - loglikelihood [double] scalar, value of the log likelihood.
% 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/>.
logpostkernel = -postkernelfun(xparam);
logprior = Prior.density(xparam);
loglikelihood = logpostkernel-logprior;
tlogpostkernel = lambda*loglikelihood + logprior;

View File

@ -1,22 +1,22 @@
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)
function oo_ = fill_mh_mode(xparam1, stdh, M_, options_, estim_params_, oo_, field_name)
% Fill oo_.<field_name>.mode and oo_.<field_name>.std_at_mode
%
% 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_}).
% - xparam1 [double] p×1 vector, estimated posterior mode.
% - stdh [double] p×1 vector, estimated posterior standard deviation.
% - M_ [struct] Description of the model.
% - estim_params_ [struct] Description of the estimated parameters.
% - options_ [struct] Dynare's options.
% - oo_ [struct] Estimation and simulation results.
%
% OUTPUTS
% o oo_ Matlab's structure gathering the results
% - oo_ Matlab's structure gathering the results
%
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2005-2021 Dynare Team
% Copyright © 2005-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -42,7 +42,8 @@ np = estim_params_.np ; % Number of deep parameters.
if np
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
name = bayestopt_.name{ip};
k = estim_params_.param_vals(i,1);
name = M_.param_names{k};
oo_.([field_name '_mode']).parameters.(name) = xparam1(ip);
oo_.([field_name '_std_at_mode']).parameters.(name) = stdh(ip);
ip = ip+1;
@ -90,4 +91,4 @@ if ncn
oo_.([field_name '_std_at_mode']).measurement_errors_corr.(name) = stdh(ip);
ip = ip+1;
end
end
end

View File

@ -1,84 +0,0 @@
function ftest (s1,s2)
% Copyright © 2001-2017 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/>.
global nvx nvy x y lag1
if size(s1,1) ~= 2
error ('Spécifiez deux fichiers pour la comparaison.') ;
end
for i = 1:2
if ~ isempty(find(abs(s1(i,:)) == 46))
error ('Entrez les noms de fichiers sans extensions.') ;
end
end
s1 = [s1 [' ';' ']] ;
file1 = [s1(1,1:min(find(abs(s1(1,:)) == 32))-1) '.BIN'] ;
file2 = [s1(2,1:min(find(abs(s1(2,:)) == 32))-1) '.BIN'] ;
fid=fopen(file1,'r') ;
n1 = fread(fid,1,'int') ;
n2 = fread(fid,1,'int') ;
n3 = fread(fid,1,'int') ;
lag1 = fread(fid,4,'int') ;
nvx = fread(fid,[n1,n3],'int') ;
x = fread(fid,[n1,n2],'float64') ;
fclose(fid) ;
nvx = char(nvx) ;
fid=fopen(file2,'r') ;
n1 = fread(fid,1,'int') ;
n2 = fread(fid,1,'int') ;
n3 = fread(fid,1,'int') ;
lag2 = fread(fid,4,'int') ;
nvy = fread(fid,[n1,n3],'int') ;
y = fread(fid,[n1,n2],'float64') ;
fclose(fid) ;
nvy = char(nvy) ;
if size(x,1) ~= size(y,1)
error ('FTEST: The two files don''t have the same number of variables.');
end
for i = 1:size(x,1)
if ~ strcmp(nvx(i,:),nvy(i,:))
error ('FTEST: The two files don''t have the same variables.') ;
end
end
if nnz(lag1 - lag2) > 0
error ('FTEST: Leads and lags aren''t the same in both files.') ;
end
j = zeros(size(s2,1),1);
for i=1:size(s2,1)
k = strmatch(s2(i,:),nvx,'exact') ;
if isempty(k)
t = ['FTEST: Variable ' s2(i) 'doesn''t exist'] ;
error (t) ;
else
j(i) =k;
end
end
y = y(j,:) ;
x = x(j,:) ;
%06/18/01 MJ replaced beastr by strmatch

View File

@ -1,52 +0,0 @@
function gcompare(s1,s2)
% GCOMPARE : GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
% This optional command plots the trajectories of a list of
% variables in two different simulations. One plot is drawn
% for each variable. The trajectories must have been previously
% saved by the instruction DYNASAVE. The simulation in file1
% is refered to as the base simulation.
% Copyright © 2001-2017 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/>.
global options_ M_
global nvx nvy x y lag1
ftest(s1,s2) ;
ix = [1-lag1(1):size(x,2)-lag1(1)]' ;
i = [lag1(1):size(ix,1)-lag1(2)+1]' ;
if options_.smpl == 0
i = [M_.maximum_lag:size(y,2)]' ;
else
i = [options_.smpl(1):options_.smpl(2)]' ;
end
for k = 1:size(x,1)
figure ;
plot (ix(i),x(k,i),ix(i),y(k,i)) ;
xlabel (['Periods']) ;
title (['Variable ' s2(k,:)]) ;
l = min(i) + 1;
ll = max(i) - 1 ;
text (l,x(k,l),s1(1,:)) ;
text (ll,y(k,ll),s1(2,:)) ;
end
% 06/18/01 MJ corrected treatment of options_.smpl
% 06/24/01 MJ removed color specification

View File

@ -1,255 +0,0 @@
function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
%
% computes the hessian matrix of the log-likelihood function of
% a state space model (notation as in kalman_filter.m in DYNARE
% Thanks to Nikolai Iskrev
%
% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
% dimension equal to the number of structural parameters
% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
% two dimensions equal to the number of structural parameters
% Copyright © 2011-2017 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/>.
k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size.
pp = size(Y,1); % Maximum number of observed variables.
mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
oldK = 0;
notsteady = 1; % Steady state flag.
F_singular = 1;
Hess = zeros(k,k); % Initialization of the Hessian
Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k);
D2a = zeros(mm,k,k); % State vector.
D2v = zeros(length(mf),k,k);
C = zeros(length(mf),mm);
for ii=1:length(mf); C(ii,mf(ii))=1;end % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
dC = zeros(length(mf),mm,k);
d2C = zeros(length(mf),mm,k,k);
s = zeros(pp,1); % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
ds = zeros(pp,1,k);
d2s = zeros(pp,1,k,k);
% for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady & t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+Om;
end
else
F_singular = 0;
iF = inv(F);
K = P(:,mf)*iF;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
[D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
tmp = (a+K*v);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii) - DYss(mf,ii);
% dai = da(:,:,ii);
dKi = DK(:,:,ii);
diFi = -iF*DF(:,:,ii)*iF;
dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
for jj = 1:ii
dFj = DF(:,:,jj);
diFj = -iF*DF(:,:,jj)*iF;
dKj = DK(:,:,jj);
d2Kij = D2K(:,:,jj,ii);
d2Fij = D2F(:,:,jj,ii);
d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
d2vij = -D2Yss(mf,jj,ii) - D2a(mf,jj,ii);
d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
end
Da(:,ii) = DT(:,:,ii)*tmp + T*dtmpi;
end
% vecDPmf = reshape(DP(mf,mf,:),[],k);
% iPmf = inv(P(mf,mf));
if t>=start
Hess = Hess + Hesst;
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
D2P = D2P1;
end
notsteady = max(max(abs(K-oldK))) > riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
tmp = (a+K*v);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
dKi = DK(:,:,ii);
diFi = -iF*DF(:,:,ii)*iF;
dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
for jj = 1:ii
dFj = DF(:,:,jj);
diFj = -iF*DF(:,:,jj)*iF;
dKj = DK(:,:,jj);
d2Kij = D2K(:,:,jj,ii);
d2Fij = D2F(:,:,jj,ii);
d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
d2vij = -D2Yss(mf,jj,ii) - D2a(mf,jj,ii);
d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
end
Da(:,ii) = DT(:,:,ii)*tmp + T*dtmpi;
end
if t>=start
Hess = Hess + Hesst;
end
a = T*(a+K*v);
end
% Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
% for ii = 1:k;
% for jj = 1:ii
% H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
% end
% end
end
Hess = Hess + tril(Hess,-1)';
Hess = -Hess/2;
% end of main function
function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
% computes (i,j) term in the Hessian
Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
% end of getHesst_ij
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
k = size(DT,3);
tmp = P-K*P(mf,:);
for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
% end of computeDKalman
function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0)
% computes the second derivatives of the Kalman matrices
% note: A=T in main func.
k = size(dA,3);
tmp = P0-K0*P0(mf,:);
[ns,no] = size(K0);
% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
% APC = A*P0*C';
% APA = A*P0*A';
d2K = zeros(ns,no,k,k);
d2S = zeros(no,no,k,k);
d2P1 = zeros(ns,ns,k,k);
for ii = 1:k
dAi = dA(:,:,ii);
dFi = dP0(mf,mf,ii);
d2Omi = d2Om(:,:,ii);
diFi = -iF*dFi*iF;
dKi = dK0(:,:,ii);
for jj = 1:k
dAj = dA(:,:,jj);
dFj = dP0(mf,mf,jj);
d2Omj = d2Om(:,:,jj);
dFj = dP0(mf,mf,jj);
diFj = -iF*dFj*iF;
dKj = dK0(:,:,jj);
d2Aij = d2A(:,:,jj,ii);
d2Pij = d2P0(:,:,jj,ii);
d2Omij = d2Om(:,:,jj,ii);
% second order
d2Fij = d2Pij(mf,mf) ;
% d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi';
d2APC = d2Pij(:,mf);
d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj;
d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ;
dtmpi = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii);
dtmpj = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj);
d2tmp = d2Pij - d2KCP;
d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi';
d2K(:,:,ii,jj) = d2Kij; %#ok<NASGU>
d2P1(:,:,ii,jj) = d2AtmpA + d2Omij; %#ok<*NASGU>
d2S(:,:,ii,jj) = d2Fij;
% d2iS(:,:,ii,jj) = d2iF;
end
end
% end of computeD2Kalman

View File

@ -34,11 +34,16 @@ end
% empirical cdfs.
xmix= [x1;x2];
bin = [-inf ; sort(xmix) ; inf];
ncount1 = histc (x1 , bin);
ncount1 = ncount1(:);
ncount2 = histc (x2 , bin);
ncount2 = ncount2(:);
if isoctave
ncount1 = histc(x1 , bin);
else
ncount1 = histcounts(x1 , bin);
end
if isoctave
ncount2 = histc(x2 , bin);
else
ncount2 = histcounts(x2 , bin);
end
cum1 = cumsum(ncount1)./sum(ncount1);
cum1 = cum1(1:end-1);

View File

@ -170,7 +170,7 @@ end
jweak = zeros(1,param_nbr);
jweak_pair = zeros(param_nbr,param_nbr);
if test_flag ~= 0 || test_flag ~= 0
if test_flag ~= 0
% these tests only apply to Jacobians, not to Gram matrices, i.e. Hessian-type or 'covariance' matrices
Pco = NaN(param_nbr,param_nbr);
for ii = 1:size(Xparnonzero,2)

View File

@ -1,9 +1,6 @@
function S = shiftS(S,n)
%function S = shiftS(S,n)
%
% Removes the first n elements of a one dimensional cell array.
function bool = issmc(options_)
% Copyright © 2013-2014 Dynare Team
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
@ -20,9 +17,4 @@ function S = shiftS(S,n)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if length(S) >= n+1
S = S(n+1:end);
else
S = {};
end
end
bool = ishssmc(options_) || isdsmh(options_);

View File

@ -7,10 +7,8 @@ function record = load_last_mh_history_file(MetropolisFolder, ModelName)
% Outputs:
% record [struct] structure storing the MH history
%
% Notes: The record structure is written to the caller workspace via an
% assignin statement.
% Copyright © 2013-2017 Dynare Team
% Copyright © 2013-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -36,7 +34,7 @@ mh_history_files = dir([BaseName '_mh_history_*.mat']);
% Consistency with older versions of Dynare.
if isequal(length(mh_history_files),0)
if exist([BaseName '_mh_history.mat'])
if exist([BaseName '_mh_history.mat'],'file')
format_mh_history_file = 1; % old Dynare format
else
error(['Estimation::load_mh_file: I cannot find any mh-history file in ' MetropolisFolder '!'])
@ -46,7 +44,7 @@ else
end
if format_mh_history_file %needed to preserve backward compatibility
load([BaseName '_mh_history.mat']);
load([BaseName '_mh_history.mat'],'record');
record.LastLogPost = record.LastLogLiK;
record.InitialLogPost = record.InitialLogLiK;
record.LastSeeds = record.Seeds;
@ -60,7 +58,7 @@ if format_mh_history_file %needed to preserve backward compatibility
record = rmfield(record,'AcceptationRates');
save([BaseName '_mh_history_0.mat'],'record');
else
load([BaseName '_mh_history_' num2str(length(mh_history_files)-1) '.mat']);
load([BaseName '_mh_history_' num2str(length(mh_history_files)-1) '.mat'],'record');
% add fields that have later been introduced
if ~isfield(record,'MCMCConcludedSuccessfully')
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...

View File

@ -60,7 +60,7 @@ xparam1 = posterior_mean;
hh = inv(SIGMA);
fprintf(' Done!\n');
if ~isfield(oo_,'posterior_mode') || (options_.mh_replic && isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice'))
oo_=fill_mh_mode(posterior_mode',NaN(npar,1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
oo_=fill_mh_mode(posterior_mode',NaN(npar,1),M_,options_,estim_params_,oo_,'posterior');
end
% save the posterior mean and the inverse of the covariance matrix

View File

@ -1,9 +1,18 @@
function [xparams,lpd,hessian_mat] = ...
maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,options_,M_,bayestopt_,estim_params_,oo_)
function [xparams, lpd, hessian_mat] = ...
maximize_prior_density(iparams, names, options_, M_, Prior, estim_params_, oo_)
% Maximizes the logged prior density using Chris Sims' optimization routine.
%
% INPUTS
% iparams [double] vector of initial parameters.
% - iparams [double] vector of initial parameters.
% - Prior [dprior] vector specifying prior densities shapes.
% - DynareOptions [struct] Options, AKA options_
% - DynareModel [struct] Model description, AKA M_
% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_
% - DynareResults [struct] Results, AKA oo_
%
%
% prior_shape [integer] vector specifying prior densities shapes.
% prior_hyperparameter_1 [double] vector, first hyperparameter.
% prior_hyperparameter_2 [double] vector, second hyperparameter.
@ -37,10 +46,18 @@ function [xparams,lpd,hessian_mat] = ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[xparams, lpd, exitflag, hessian_mat]=dynare_minimize_objective('minus_logged_prior_density', ...
iparams, options_.mode_compute, options_, [prior_inf_bound, prior_sup_bound], ...
bayestopt_.name, bayestopt_, [], ...
prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound, ...
options_,M_,estim_params_,oo_);
[xparams, lpd, exitflag, hessian_mat] = dynare_minimize_objective('minus_logged_prior_density', ...
iparams, ...
options_.mode_compute, ...
options_, ...
[Prior.p3, Prior.p4], ...
names, ...
[], ...
[], ...
Prior,
options_, ...
M_, ...
estim_params_, ...
oo_);
lpd = -lpd;

View File

@ -1,8 +0,0 @@
function [res,jac,domerr] = mcpath_function(func,z,jacflag,varargin)
domerr = 0;
if jacflag
[res,jac] = func(z,varargin{:});
else
res = func(z,varargin{:});
end

View File

@ -1,55 +0,0 @@
function metropolis_run_analysis(M_,basetopt,j)
%function metropolis_run_analysis(M_,basetopt,j
% analizes Metropolis runs
%
% INPUTS
% M_: (struct) Model structure
% basetopt: (struct) Estimated parameters structure
% j: (int) Index of estimated paramter
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Copyright © 2003-2017 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/>.
load([M_.fname '/metropolis/' M_.fname '_mh_history'])
nblck = record.Nblck;
ndraws = sum(record.MhDraws(:,1));
logPost = [];
params = [];
blck = 1;
for i=1:record.LastFileNumber
fname = [M_.fname '/metropolis/' M_.fname '_mh' int2str(i) '_blck' ...
int2str(blck) '.mat'];
if exist(fname,'file')
o=load(fname);
logPost = [logPost; o.logpo2];
params = [params; o.x2];
end
end
figure;
subplot(2,1,1)
plot(logPost)
subplot(2,1,2)
plot(params(:,j))
title(['parameter ' int2str(j)])

View File

@ -59,7 +59,7 @@ nblck = size(record.LastParameters,1);
clear record;
% Get all the posterior draws:
PosteriorDraws = GetAllPosteriorDraws(M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck);
PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck);
% Compute the autocorrelation function:
[~,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size);

View File

@ -1,24 +1,19 @@
function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
% [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
function [fval, info, exitflag, ~, ~] = minus_logged_prior_density(xparams, Prior, options_, M_, estim_params_, oo_)
% Evaluates minus the logged prior density.
%
%
% INPUTS
% xparams [double] vector of parameters.
% pshape [integer] vector specifying prior densities shapes.
% p6 [double] vector, first hyperparameter.
% p7 [double] vector, second hyperparameter.
% p3 [double] vector, prior's lower bound.
% p4 [double] vector, prior's upper bound.
% prior_sup_bound [double] vector, prior's upper bound.
% options_ [structure] describing the options
% M_ [structure] describing the model
% estim_params_ [structure] characterizing parameters to be estimated
% oo_ [structure] storing the results
% - xparams [double] vector of parameters.
% - Prior [dprior] vector specifying prior densities shapes.
% - DynareOptions [struct] Options, AKA options_
% - DynareModel [struct] Model description, AKA M_
% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_
% - DynareResults [struct] Results, AKA oo_
%
% OUTPUTS
% f [double] value of minus the logged prior density.
% info [double] vector: second entry stores penalty, first entry the error code.
%
% - fval [double] value of minus the logged prior density.
% - info [double] 4×1 vector, second entry stores penalty, first entry the error code, last entry a penalty (used for optimization).
% Copyright © 2009-2023 Dynare Team
%
% This file is part of Dynare.
@ -36,10 +31,7 @@ function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparam
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fake_1 = [];
fake_2 = [];
exit_flag = 1;
exitflag = true;
info = zeros(4,1);
%------------------------------------------------------------------------------
@ -47,74 +39,75 @@ info = zeros(4,1);
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if ~isequal(options_.mode_compute,1) && any(xparams<p3)
k = find(xparams<p3);
if ~isequal(options_.mode_compute, 1) && any(xparams<Prior.p3)
k = find(xparams<Prior.p3);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 41;
info(4) = sum((p3(k)-xparams(k)).^2);
info(4) = sum((Prior.p3(k)-xparams(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if ~isequal(options_.mode_compute,1) && any(xparams>p4)
k = find(xparams>p4);
if ~isequal(options_.mode_compute, 1) && any(xparams>Prior.p4)
k = find(xparams>Prior.p4);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 42;
info(4) = sum((xparams(k)-p4(k)).^2);
info(4) = sum((xparams(k)-Prior.p4(k)).^2);
return
end
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
M_ = set_all_parameters(xparams,estim_params_,M_);
M_ = set_all_parameters(xparams, estim_params_, M_);
Q = M_.Sigma_e;
H = M_.H;
% Test if Q is positive definite.
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covariances')
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_, 'calibrated_covariances')
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
[Q_is_positive_definite, penalty] = ispd(Q);
if ~Q_is_positive_definite
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the
% eigenvalues of this matrix in order to build the endogenous penalty.
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(estim_params_,'calibrated_covariances')
correct_flag=check_consistency_covariances(Q);
if isfield(estim_params_, 'calibrated_covariances')
correct_flag = check_consistency_covariances(Q);
if ~correct_flag
penalty = sum(Q(estim_params_.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 71;
info(4) = penalty;
return4
end
end
end
% Test if H is positive definite.
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covariances_ME')
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_, 'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H);
if ~H_is_positive_definite
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues
% of this matrix in order to build the endogenous penalty.
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 44;
info(4) = penalty;
return
end
if isfield(estim_params_,'calibrated_covariances_ME')
correct_flag=check_consistency_covariances(H);
if isfield(estim_params_, 'calibrated_covariances_ME')
correct_flag = check_consistency_covariances(H);
if ~correct_flag
penalty = sum(H(estim_params_.calibrated_covariances_ME.position).^2);
fval = Inf;
exit_flag = 0;
exitflag = false;
info(1) = 72;
info(4) = penalty;
return
@ -127,7 +120,7 @@ end
% 2. Check BK and steady state
%-----------------------------
[~,info] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
[~, info] = resol(0, M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
@ -137,14 +130,14 @@ if info(1)
%meaningful second entry of output that can be used
fval = Inf;
info(4) = info(2);
exit_flag = 0;
exitflag = false;
return
else
fval = Inf;
info(4) = 0.1;
exit_flag = 0;
exitflag = false;
return
end
end
fval = - priordens(xparams,pshape,p6,p7,p3,p4);
fval = - Prior.density(xparams);

@ -1 +1 @@
Subproject commit 80446e7cdcf392167e91b428b61d77e028e47674
Subproject commit a3874403fe5b403f74d727006cb4ccf198cf63c8

View File

@ -1,5 +1,5 @@
function optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
% optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
function optimize_prior(options_, M_, oo_, Prior, estim_params_, pnames)
% This routine computes the mode of the prior density using an optimization algorithm.
%
% INPUTS
@ -26,24 +26,25 @@ function optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
oo_.dr = set_state_space(oo_.dr, M_, options_);
% Initialize to the prior mean
oo_.dr = set_state_space(oo_.dr,M_);
xparam1 = bayestopt_.p1;
xparam1 = Prior.p1;
% Pertubation of the initial condition.
look_for_admissible_initial_condition = 1; scale = 1.0; iter = 0;
look_for_admissible_initial_condition = true; scale = 1.0; iter = 0;
while look_for_admissible_initial_condition
xinit = xparam1+scale*randn(size(xparam1));
if all(xinit(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
M_ = set_all_parameters(xinit,estim_params_,M_);
[oo_.dr,INFO,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if all(xinit>Prior.p3) && all(xinit<Prior.p4)
M_ = set_all_parameters(xinit, estim_params_, M_);
[dr, INFO, M_, oo_] = resol(0, M_, options_, oo_);
if ~INFO(1)
look_for_admissible_initial_condition = 0;
look_for_admissible_initial_condition = false;
end
else
if iter == 2000
if iter==2000
scale = scale/1.1;
iter = 0;
iter = 0;
else
iter = iter+1;
end
@ -52,23 +53,19 @@ end
% Maximization of the prior density
[xparams, lpd, hessian_mat] = ...
maximize_prior_density(xinit, bayestopt_.pshape, ...
bayestopt_.p6, ...
bayestopt_.p7, ...
bayestopt_.p3, ...
bayestopt_.p4,options_,M_,bayestopt_,estim_params_,oo_);
maximize_prior_density(xinit, pnames, options_, M_, Prior, estim_params_, oo_);
% Display the results.
% Display results.
skipline(2)
disp('------------------')
disp('PRIOR OPTIMIZATION')
disp('------------------')
skipline()
for i = 1:length(xparams)
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0,M_,estim_params_,options_.varobs) '.'])
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
dprintf('deep parameter %u: %s.', i, get_the_name(i, 0, M_, estim_params_, options_.varobs))
dprintf(' Initial condition ........ %s.', num2str(xinit(i)))
dprintf(' Prior mode ............... %s.', num2str(Prior.p5(i)))
dprintf(' Optimized prior mode ..... %s.', num2str(xparams(i)))
skipline()
end
skipline()
skipline()

View File

@ -536,6 +536,11 @@ b = size(z,3);
if ~isempty(options_.plot_shock_decomp.plot_init_date)
my_initial_date = max(initial_date,options_.plot_shock_decomp.plot_init_date);
a = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_init_date);
if isempty(a)
warning(['You set plot_init_date larger than the last observation. The last observation is %s,\n' ...
'while you requested as the first observation %s.Exiting because there is nothing to do.'],initial_date+b,options_.plot_shock_decomp.plot_init_date);
return;
end
end
if ~isempty(options_.plot_shock_decomp.plot_end_date)
if options_.plot_shock_decomp.plot_end_date<=(max(initial_date:initial_date+b-1))

View File

@ -1,8 +1,7 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_, dispString)
% function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_, dispString)
% Metropolis-Hastings initialization.
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_, dispString)
% Posterior sampler initialization.
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
@ -257,7 +256,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('%s: Initialization at the posterior mode.\n\n',dispString);
fprintf('%s: Initialization at the posterior mode.\n\n',dispString);
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
for i=1:length(ix2(1,:))
fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']);

View File

@ -1,51 +1,13 @@
function bounds = prior_bounds(bayestopt, prior_trunc)
function bounds = prior_bounds(bayestopt_, priortrunc)
%@info:
%! @deftypefn {Function File} {@var{bounds} =} prior_bounds (@var{bayesopt},@var{option})
%! @anchor{prior_bounds}
%! @sp 1
%! Returns bounds for the prior densities. For each estimated parameter the lower and upper bounds
%! are such that the defined intervals contains a probability mass equal to 1-2*@var{option}.prior_trunc. The
%! default value for @var{option}.prior_trunc is 1e-10 (set in @ref{global_initialization}).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item bayestopt
%! Matlab's structure describing the prior distribution (initialized by @code{dynare}).
%! @item option
%! Matlab's structure describing the options (initialized by @code{dynare}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item bounds
%! A structure with two fields lb and up (p*1 vectors of doubles, where p is the number of estimated parameters) for the lower and upper bounds.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{get_prior_info}, @ref{dynare_estimation_1}, @ref{dynare_estimation_init}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! None.
%! @end deftypefn
%@eod:
% function bounds = prior_bounds(bayestopt)
% computes bounds for prior density.
%
% INPUTS
% bayestopt [structure] characterizing priors (shape, mean, p1..p4)
% - bayestopt [struct] characterizing priors (shape, mean, p1..p4)
% - priortrunc [double] scalar, probability mass in the tails to be removed
%
% OUTPUTS
% bounds [double] structure specifying prior bounds (lb and ub fields)
%
% SPECIAL REQUIREMENTS
% none
% - bounds [struct] prior bounds (lb, lower bounds, and ub, upper bounds, fields are n×1 vectors)
% Copyright © 2003-2023 Dynare Team
%
@ -64,74 +26,78 @@ function bounds = prior_bounds(bayestopt, prior_trunc)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
pshape = bayestopt.pshape;
p3 = bayestopt.p3;
p4 = bayestopt.p4;
p6 = bayestopt.p6;
p7 = bayestopt.p7;
if nargin<2, priortrunc = 0.0; end
bounds.lb = zeros(length(p6),1);
bounds.ub = zeros(length(p6),1);
assert(priortrunc>=0 && priortrunc<=1, 'Second input argument must be non negative and not larger than one.')
pshape = bayestopt_.pshape;
p3 = bayestopt_.p3;
p4 = bayestopt_.p4;
p6 = bayestopt_.p6;
p7 = bayestopt_.p7;
bounds.lb = zeros(size(p6));
bounds.ub = zeros(size(p6));
for i=1:length(p6)
switch pshape(i)
case 1
if prior_trunc == 0
if priortrunc==0
bounds.lb(i) = p3(i);
bounds.ub(i) = p4(i);
else
bounds.lb(i) = betainv(prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
bounds.ub(i) = betainv(1-prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
bounds.lb(i) = betainv(priortrunc, p6(i), p7(i))*(p4(i)-p3(i))+p3(i);
bounds.ub(i) = betainv(1.0-priortrunc, p6(i), p7(i))*(p4(i)-p3(i))+p3(i);
end
case 2
if prior_trunc == 0
if priortrunc==0
bounds.lb(i) = p3(i);
bounds.ub(i) = Inf;
else
bounds.lb(i) = gaminv(prior_trunc,p6(i),p7(i))+p3(i);
bounds.ub(i) = gaminv(1-prior_trunc,p6(i),p7(i))+p3(i);
bounds.lb(i) = gaminv(priortrunc, p6(i), p7(i))+p3(i);
bounds.ub(i) = gaminv(1.0-priortrunc, p6(i), p7(i))+p3(i);
end
case 3
if prior_trunc == 0
bounds.lb(i) = max(-Inf,p3(i));
bounds.ub(i) = min(Inf,p4(i));
if priortrunc == 0
bounds.lb(i) = max(-Inf, p3(i));
bounds.ub(i) = min(Inf, p4(i));
else
bounds.lb(i) = max(norminv(prior_trunc,p6(i),p7(i)),p3(i));
bounds.ub(i) = min(norminv(1-prior_trunc,p6(i),p7(i)),p4(i));
bounds.lb(i) = max(norminv(priortrunc, p6(i), p7(i)), p3(i));
bounds.ub(i) = min(norminv(1-priortrunc, p6(i), p7(i)), p4(i));
end
case 4
if prior_trunc == 0
if priortrunc==0
bounds.lb(i) = p3(i);
bounds.ub(i) = Inf;
else
bounds.lb(i) = 1/sqrt(gaminv(1-prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
bounds.ub(i) = 1/sqrt(gaminv(prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
bounds.lb(i) = 1.0/sqrt(gaminv(1.0-priortrunc, p7(i)/2.0, 2.0/p6(i)))+p3(i);
bounds.ub(i) = 1.0/sqrt(gaminv(priortrunc, p7(i)/2.0, 2.0/p6(i)))+p3(i);
end
case 5
if prior_trunc == 0
if priortrunc == 0
bounds.lb(i) = p6(i);
bounds.ub(i) = p7(i);
else
bounds.lb(i) = p6(i)+(p7(i)-p6(i))*prior_trunc;
bounds.ub(i) = p7(i)-(p7(i)-p6(i))*prior_trunc;
bounds.lb(i) = p6(i)+(p7(i)-p6(i))*priortrunc;
bounds.ub(i) = p7(i)-(p7(i)-p6(i))*priortrunc;
end
case 6
if prior_trunc == 0
if priortrunc == 0
bounds.lb(i) = p3(i);
bounds.ub(i) = Inf;
else
bounds.lb(i) = 1/gaminv(1-prior_trunc, p7(i)/2, 2/p6(i))+p3(i);
bounds.ub(i) = 1/gaminv(prior_trunc, p7(i)/2, 2/p6(i))+ p3(i);
bounds.lb(i) = 1.0/gaminv(1.0-priortrunc, p7(i)/2.0, 2.0/p6(i))+p3(i);
bounds.ub(i) = 1.0/gaminv(priortrunc, p7(i)/2.0, 2.0/p6(i))+ p3(i);
end
case 8
if prior_trunc == 0
if priortrunc == 0
bounds.lb(i) = p3(i);
bounds.ub(i) = Inf;
else
bounds.lb(i) = p3(i)+wblinv(prior_trunc,p6(i),p7(i));
bounds.ub(i) = p3(i)+wblinv(1-prior_trunc,p6(i),p7(i));
bounds.lb(i) = p3(i)+wblinv(priortrunc, p6(i), p7(i));
bounds.ub(i) = p3(i)+wblinv(1.0-priortrunc, p6(i), p7(i));
end
otherwise
error(sprintf('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i)));
error('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i));
end
end
end

View File

@ -212,9 +212,9 @@ if strcmpi(type,'posterior')
mh_nblck = options_.mh_nblck;
if B==NumberOfDraws*mh_nblck
% we load all retained MH runs !
logpost=GetAllPosteriorDraws(M_.dname, M_.fname, 0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
logpost=GetAllPosteriorDraws(options_, M_.dname, M_.fname, 0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
for column=1:npar
x(:,column) = GetAllPosteriorDraws(M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
x(:,column) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
end
else
logpost=NaN(B,1);
@ -390,4 +390,4 @@ if ~isnumeric(options_.parallel)
if leaveSlaveOpen == 0
closeSlave(options_.parallel,options_.parallel_info.RemoteTmpFolder),
end
end
end

View File

@ -1,123 +0,0 @@
function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
% function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
%
% computes the derivative of the log-likelihood function of
% a state space model (notation as in kalman_filter.m in DYNARE
% thanks to Nikolai Iskrev
%
% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
% dimension equal to the number of structural parameters
% Copyright © 2009-2017 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/licen
k = size(DT,3); % number of structural parameters
smpl = size(Y,2); % Sample size.
mm = size(T,2); % Number of state variables.
a = zeros(mm,1); % State vector.
Om = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = 0; % Initialization of the time index.
oldK = 0;
notsteady = 1; % Steady state flag.
F_singular = 1;
DLIK = zeros(k,1); % Initialization of the score.
Da = zeros(mm,k); % State vector.
Dv = zeros(length(mf),k); % observation vector.
% for ii = 1:k
% DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
% end
while notsteady & t<smpl
t = t+1;
v = Y(:,t)-a(mf);
F = P(mf,mf) + H;
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+Om;
end
else
F_singular = 0;
iF = inv(F);
K = P(:,mf)*iF;
[DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
a = T*(a+K*v);
P = T*(P-K*P(mf,:))*transpose(T)+Om;
DP = DP1;
end
notsteady = max(max(abs(K-oldK))) > riccati_tol;
oldK = K;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
for ii = 1:k
tmp0(:,:,ii) = iF*DF(:,:,ii)*iF;
end
if t < smpl
t0 = t+1;
while t < smpl
t = t+1;
v = Y(:,t)-a(mf);
for ii = 1:k
Dv(:,ii) = -Da(mf,ii)-DYss(mf,ii);
Da(:,ii) = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
if t>=start
DLIK(ii,1) = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
end
end
a = T*(a+K*v);
end
for ii = 1:k
% DLIK(ii,1) = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
end
end
DLIK = DLIK/2;
% end of main function
function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
k = size(DT,3);
tmp = P-K*P(mf,:);
for ii = 1:k
DF(:,:,ii) = DP(mf,mf,ii) + DH(:,:,ii);
DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
DK(:,:,ii) = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
Dtmp = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
end
% end of computeDKalman

59
matlab/substitute.m Normal file
View File

@ -0,0 +1,59 @@
function v = substitute(v, i, x)
% Substitute a scalar in a vector.
%
% INPUTS
% - v [double] m×1 vector
% - i [integer] scalar, index for the scalar to be replaced
% - x [double] scalar or 1×n vector.
%
% OUTPUTS
% - v [double] m×1 vector or m×n matrix (with substituted value(s))
%
% REMARKS
% If x is a vector with n elements, then n substitutions are performed (returning n updated vectors in a matrix with n columns)
%
% EXAMPLES
% >> v = ones(2,1);
% >> substitude(v, 1, 0)
%
% ans = %
%
% 0
%
% 1
%
% >> substitute(v, 1, [3 4])
%
% ans =
%
% 3 4
% 1 1
% Copyright © 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/>.
assert(isvector(v), 'First input argument must be a vector.')
assert(isvector(x), 'Last input argument must be a scalar or a vector.')
assert(isscalar(i) && isint(i) && i>0 && i<=length(v), 'Second input argument must be a scalar integer')
if length(x)==1
v(i) = x;
else
v = repmat(v, 1, length(x));
v(i,:) = x;
end

View File

@ -67,7 +67,7 @@ n_nblocks_to_plot=length(blck);
if n_nblocks_to_plot==1
% Get all the posterior draws:
PosteriorDraws = GetAllPosteriorDraws(M_.dname,M_.fname,column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
else
PosteriorDraws=NaN(TotalNumberOfMhDraws,n_nblocks_to_plot);
save_string='';
@ -75,7 +75,7 @@ else
title_string_tex='';
end
for block_iter=1:n_nblocks_to_plot
PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
save_string=[save_string,'_',num2str(blck(block_iter))];
if options_.TeX
title_string_tex=[title_string_tex, ', ' num2str(blck(block_iter))];

View File

@ -1,4 +1,5 @@
function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)%datafile,varobs,xls_sheet,xls_range)
function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)
% dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)
% Computes the sample moments of a VAR model.
%
% The VAR(p) model is defined by:
@ -51,7 +52,7 @@ function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, datase
% SPECIAL REQUIREMENTS
% None.
% Copyright © 2007-2021 Dynare Team
% Copyright © 2007-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -80,14 +81,11 @@ if isequal(var_trend_order,-1)
elseif isequal(var_trend_order,0)
% Constant and no linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*nlag+1);
indx = NumberOfVariables*nlag+1;
elseif isequal(var_trend_order,1)
% Constant and linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*nlag+2);
indx = NumberOfVariables*nlag+1:NumberOfVariables*nlag+2;
else
error('Estimation::var_sample_moments: trend must be equal to -1,0 or 1!')
return
end
% I build matrices Y and X
@ -108,10 +106,4 @@ dataset_info.mYX=Y'*X;
dataset_info.mXY=X'*Y;
dataset_info.mXX=X'*X;
dataset_info.Ydata=Y;
dataset_info.Xdata=X;
% assignin('base', 'mYY', Y'*Y);
% assignin('base', 'mYX', Y'*X);
% assignin('base', 'mXY', X'*Y);
% assignin('base', 'mXX', X'*X);
% assignin('base', 'Ydata', Y);
% assignin('base', 'Xdata', X);
dataset_info.Xdata=X;

View File

@ -457,7 +457,8 @@ shared_module('k_order_welfare', k_order_welfare_src, kwargs : korder_mex_kwargs
# Unit tests
k_order_test_kwargs = { 'include_directories' : [ mex_incdir, korder_incdir ],
k_order_test_kwargs = { 'include_directories' : mex_kwargs['include_directories'] + korder_incdir,
'cpp_args' : mex_kwargs['cpp_args'],
'link_args' : exe_link_args,
'build_rpath' : exe_rpath,
'link_with' : korder_lib }
@ -804,6 +805,8 @@ mod_and_m_tests = [
'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'estimation/fs2000.mod' ],
'extra' : [ 'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'estimation/hssmc/fs2000.mod' ],
'extra' : [ 'estimation/fsdat_simul.m' ] },
{ 'test' : [ 'gsa/ls2003a.mod',
'gsa/ls2003.mod',
'gsa/ls2003scr.mod',

View File

@ -1,5 +1,5 @@
/*
* Copyright © 2008-2022 Dynare Team
* Copyright © 2008-2023 Dynare Team
*
* This file is part of Dynare.
*
@ -248,7 +248,7 @@ DynareStateNameList::DynareStateNameList(const KordpDynare& dynare, const Dynare
const DynareNameList& denl)
{
for (int i = 0; i < dynare.nYs; i++)
names.emplace_back(dnl.getName(i + dynare.nstat()));
for (int i = 0; i < dynare.nexog(); i++)
names.emplace_back(dnl.getName(i + dynare.nStat));
for (int i = 0; i < dynare.nExog; i++)
names.emplace_back(denl.getName(i));
}

View File

@ -90,12 +90,10 @@ public:
ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
// Create submatrix (with data sharing)
ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols);
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit ConstGeneralMatrix(const mxArray* p) :
data(p), rows {static_cast<int>(mxGetM(p))}, cols {static_cast<int>(mxGetN(p))}, ld {rows}
{
}
#endif
virtual ~ConstGeneralMatrix() = default;
ConstGeneralMatrix& operator=(const ConstGeneralMatrix& v) = delete;
@ -234,12 +232,10 @@ public:
// Create submatrix (with data sharing)
GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols);
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit GeneralMatrix(mxArray* p) :
data(p), rows {static_cast<int>(mxGetM(p))}, cols {static_cast<int>(mxGetN(p))}, ld {rows}
{
}
#endif
virtual ~GeneralMatrix() = default;
GeneralMatrix& operator=(const GeneralMatrix& m) = default;

View File

@ -125,7 +125,6 @@ SylvParams::getArrayNames() const
return names;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray*
SylvParams::DoubleParamItem::createMatlabArray() const
{
@ -136,11 +135,11 @@ mxArray*
SylvParams::IntParamItem::createMatlabArray() const
{
mxArray* res = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
# if MX_HAS_INTERLEAVED_COMPLEX
#if MX_HAS_INTERLEAVED_COMPLEX
*mxGetInt32s(res) = value;
# else
#else
*static_cast<int*>(mxGetData(res)) = value;
# endif
#endif
return res;
}
@ -225,4 +224,3 @@ SylvParams::createStructArray() const
return res;
}
#endif

View File

@ -25,9 +25,7 @@
#include <string>
#include <vector>
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
# include <dynmex.h>
#endif
#include <dynmex.h>
enum class status
{
@ -111,9 +109,7 @@ protected:
ParamItem<double>::operator=(val);
return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
[[nodiscard]] mxArray* createMatlabArray() const;
#endif
};
class IntParamItem : public ParamItem<int>
@ -132,9 +128,7 @@ protected:
ParamItem<int>::operator=(val);
return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
[[nodiscard]] mxArray* createMatlabArray() const;
#endif
};
class BoolParamItem : public ParamItem<bool>
@ -153,9 +147,7 @@ protected:
ParamItem<bool>::operator=(val);
return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
[[nodiscard]] mxArray* createMatlabArray() const;
#endif
};
class MethodParamItem : public ParamItem<solve_method>
@ -174,9 +166,7 @@ protected:
ParamItem<solve_method>::operator=(val);
return *this;
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
[[nodiscard]] mxArray* createMatlabArray() const;
#endif
};
public:
@ -224,9 +214,8 @@ public:
void print(const std::string& prefix) const;
void print(std::ostream& fdesc, const std::string& prefix) const;
[[nodiscard]] std::vector<const char*> getArrayNames() const;
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
[[nodiscard]] mxArray* createStructArray() const;
#endif
private:
void copy(const SylvParams& p);
};

View File

@ -107,14 +107,12 @@ Vector::Vector(const Vector& v, int off_arg, int skip, int l) : len(l), data {ne
copy(v.data + off_arg * v.s, v.s * skip);
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
Vector::Vector(mxArray* p) :
len {static_cast<int>(mxGetNumberOfElements(p))}, data {mxGetPr(p)}, destroy {false}
{
if (!mxIsDouble(p) || mxIsComplex(p) || mxIsSparse(p))
throw SYLV_MES_EXCEPTION("This is not a dense array of real doubles.");
}
#endif
bool
Vector::operator==(const Vector& y) const
@ -305,14 +303,12 @@ ConstVector::ConstVector(const double* d, int skip, int l) : len {l}, s {skip},
{
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
ConstVector::ConstVector(const mxArray* p) :
len {static_cast<int>(mxGetNumberOfElements(p))}, data {mxGetPr(p)}
{
if (!mxIsDouble(p))
throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles.");
}
#endif
bool
ConstVector::operator==(const ConstVector& y) const

View File

@ -28,9 +28,7 @@
#include <complex>
#include <utility>
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
# include <dynmex.h>
#endif
#include <dynmex.h>
class GeneralMatrix;
class ConstVector;
@ -72,9 +70,7 @@ public:
Vector(const Vector& v, int off_arg, int l);
Vector(Vector& v, int off_arg, int skip, int l);
Vector(const Vector& v, int off_arg, int skip, int l);
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit Vector(mxArray* p);
#endif
Vector& operator=(const Vector& v);
/* The move-assignment operator is not implemented, because moving pointers
across class instances would break the reference semantics that the
@ -198,9 +194,7 @@ public:
ConstVector(const ConstVector& v, int off_arg, int l);
ConstVector(const ConstVector& v, int off_arg, int skip, int l);
ConstVector(const double* d, int skip, int l);
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit ConstVector(const mxArray* p);
#endif
virtual ~ConstVector() = default;
ConstVector& operator=(const ConstVector& v) = delete;
ConstVector& operator=(ConstVector&& v) = delete;

View File

@ -216,7 +216,7 @@ public:
auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
ten_smart->zeros();
ten = ten_smart.get();
insert(std::move(ten_smart));
_Tparent::insert(std::move(ten_smart));
}
Symmetry sym {i, j};
@ -245,7 +245,7 @@ public:
auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
ten_smart->zeros();
ten = ten_smart.get();
insert(std::move(ten_smart));
_Tparent::insert(std::move(ten_smart));
}
Symmetry sym {0, j};

View File

@ -1,6 +1,6 @@
/*
* Copyright © 2004-2011 Ondra Kamenik
* Copyright © 2019 Dynare Team
* Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@ -63,11 +63,9 @@ public:
ConstGeneralMatrix(m, first_row, first_col, rows, cols)
{
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit ConstTwoDMatrix(const mxArray* p) : ConstGeneralMatrix(p)
{
}
#endif
~ConstTwoDMatrix() override = default;
ConstTwoDMatrix& operator=(const ConstTwoDMatrix& v) = delete;
@ -125,11 +123,9 @@ public:
{
}
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
explicit TwoDMatrix(mxArray* p) : GeneralMatrix(p)
{
}
#endif
~TwoDMatrix() override = default;
TwoDMatrix& operator=(const TwoDMatrix& m) = default;

@ -1 +1 @@
Subproject commit 378d00fc3a50a11c3dce615ea337ad55989c938e
Subproject commit cd86b1895d97d1b195849072a8d26fd372ca5d8c

View File

@ -55,39 +55,40 @@
;; Also include "end" in this list
(defvar dynare-statements
'("var" "varexo" "varexo_det" "trend_var" "log_trend_var"
"predetermined_variables" "parameters" "model_local_variable"
"model_info" "estimation" "set_time" "data" "varobs"
"varexobs" "unit_root_vars" "rplot" "osr_params" "osr" "dynatype"
"dynasave" "model_comparison" "change_type" "load_params_and_steady_state"
"save_params_and_steady_state" "write_latex_dynamic_model"
"write_latex_static_model" "write_latex_original_model"
"write_latex_steady_state_model" "steady" "check" "simul" "stoch_simul"
"var_model" "trend_component_model" "var_expectation_model" "pac_model"
"dsample" "planner_objective" "ramsey_model" "ramsey_policy"
"evaluate_planner_objective" "occbin_setup" "occbin_solver"
"occbin_write_regimes" "occbin_graph"
"predetermined_variables" "parameters" "model_local_variable" "model_info"
"estimation" "set_time" "data" "varobs" "varexobs" "unit_root_vars" "rplot"
"osr_params" "osr" "dynatype" "dynasave" "model_comparison" "change_type"
"load_params_and_steady_state" "save_params_and_steady_state"
"write_latex_dynamic_model" "write_latex_static_model"
"write_latex_original_model" "write_latex_steady_state_model" "steady"
"check" "simul" "stoch_simul" "var_model" "trend_component_model"
"var_expectation_model" "pac_model" "dsample" "planner_objective"
"ramsey_model" "ramsey_policy" "evaluate_planner_objective" "occbin_setup"
"occbin_solver" "occbin_write_regimes" "occbin_graph"
"discretionary_policy" "identification" "bvar_density" "bvar_forecast"
"dynare_sensitivity" "initval_file" "histval_file" "forecast"
"bvar_irf" "dynare_sensitivity" "initval_file" "histval_file" "forecast"
"shock_decomposition" "realtime_shock_decomposition"
"plot_shock_decomposition" "initial_condition_decomposition"
"squeeze_shock_decomposition" "sbvar"
"ms_estimation" "ms_simulation" "ms_compute_mdd" "ms_compute_probabilities"
"ms_forecast" "ms_irf" "ms_variance_decomposition" "conditional_forecast"
"plot_conditional_forecast" "method_of_moments"
"markov_switching" "svar" "svar_global_identification_check"
"external_function" "calib_smoother" "model_diagnostics" "extended_path"
"smoother2histval" "perfect_foresight_setup" "perfect_foresight_solver"
"squeeze_shock_decomposition" "sbvar" "ms_estimation" "ms_simulation"
"ms_compute_mdd" "ms_compute_probabilities" "ms_forecast" "ms_irf"
"ms_variance_decomposition" "conditional_forecast"
"plot_conditional_forecast" "method_of_moments" "markov_switching" "svar"
"svar_global_identification_check" "external_function" "calib_smoother"
"model_diagnostics" "extended_path" "smoother2histval"
"perfect_foresight_setup" "perfect_foresight_solver"
"perfect_foresight_with_expectation_errors_setup"
"perfect_foresight_with_expectation_errors_solver"
"compilation_setup" "model_remove" "model_options" "var_remove" "resid"
"std" "corr" "prior_function" "posterior_function" "end")
"perfect_foresight_with_expectation_errors_solver" "compilation_setup"
"model_remove" "model_options" "var_remove" "resid" "std" "corr"
"prior_function" "posterior_function" "end")
"Dynare statement keywords.")
;; Keywords that may appear in blocks, and that begin a statement which will be
;; closed by a semicolon
(defvar dynare-statements-like
'("stderr" "values" "periods" "scales" "restriction" "exclusion" "upper_cholesky" "lower_cholesky"
"bind" "relax" "error_bind" "error_relax" "add" "multiply")
'("stderr" "values" "periods" "scales" "restriction" "exclusion"
"upper_cholesky" "lower_cholesky" "equation" "bind" "relax" "error_bind"
"error_relax" "add" "multiply" "target" "auxname_target_nonstationary"
"component" "growth" "auxname" "kind")
"Dynare statements-like keywords.")
;; Those keywords that makes the lexer enter the DYNARE_BLOCK start condition
@ -96,13 +97,15 @@
;; referenced in another eval-when-compile statement in dynare-calculate-indentation
(eval-when-compile
(defvar dynare-blocks
'("model" "steady_state_model" "initval" "endval" "histval" "shocks" "heteroskedastic_shocks"
"shock_groups" "init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
"estimated_params_init" "estimated_params_bounds" "estimated_params_remove"
"osr_params_bounds" "observation_trends" "deterministic_trends" "optim_weights"
"homotopy_setup" "conditional_forecast_paths" "svar_identification"
"moment_calibration" "irf_calibration" "ramsey_constraints" "generate_irfs"
"matched_moments" "occbin_constraints" "model_replace" "verbatim" "pac_target_info")
'("model" "steady_state_model" "initval" "endval" "histval"
"filter_initial_state" "shocks" "heteroskedastic_shocks" "shock_groups"
"init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
"estimated_params_init" "estimated_params_bounds"
"estimated_params_remove" "osr_params_bounds" "observation_trends"
"deterministic_trends" "optim_weights" "homotopy_setup"
"conditional_forecast_paths" "svar_identification" "moment_calibration"
"irf_calibration" "ramsey_constraints" "generate_irfs" "matched_moments"
"occbin_constraints" "model_replace" "pac_target_info" "verbatim")
"Dynare block keywords."))
;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file)

View File

@ -0,0 +1,85 @@
// 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='dsmh');

View File

@ -0,0 +1,91 @@
// 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='hssmc',
posterior_sampler_options=('steps',10,
'lambda',2,
'particles', 20000,
'scale',.5,
'target', .25));

View File

@ -269,8 +269,6 @@ method_of_moments(
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
% , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
% , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

View File

@ -218,8 +218,6 @@ method_of_moments(
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
% , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
% , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

View File

@ -205,8 +205,6 @@ end
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
% , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
% , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

View File

@ -177,8 +177,6 @@ save('test_matrix.mat','weighting_matrix')
% , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
% , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver
% , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm
% , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
% , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver
% , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
% , schur_vec_tol=1e-11 % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

View File

@ -1,94 +1,94 @@
// 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.05;
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_.posterior_sampling_method = 'slice';
estimation(order=1,datafile='../fsdat_simul',nobs=192,silent_optimizer,loglinear,mh_replic=50,mh_nblocks=2,mh_drop=0.2, //mode_compute=0,cova_compute=0,
posterior_sampling_method='slice'
);
// continue with rotated slice
estimation(order=1,datafile='../fsdat_simul',silent_optimizer,nobs=192,loglinear,mh_replic=100,mh_nblocks=2,mh_drop=0.5,load_mh_file,//mode_compute=0,
posterior_sampling_method='slice',
posterior_sampler_options=('rotated',1,'use_mh_covariance_matrix',1)
);
options_.TeX=1;
generate_trace_plots(1:2);
options_.TeX=1;
// 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.05;
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_.posterior_sampling_method = 'slice';
estimation(order=1,datafile='../fsdat_simul',nobs=192,silent_optimizer,loglinear,mh_replic=50,mh_nblocks=2,mh_drop=0.2, //mode_compute=0,cova_compute=0,
posterior_sampling_method='slice'
);
// continue with rotated slice
estimation(order=1,datafile='../fsdat_simul',silent_optimizer,nobs=192,loglinear,mh_replic=100,mh_nblocks=2,mh_drop=0.5,load_mh_file,//mode_compute=0,
posterior_sampling_method='slice',
posterior_sampler_options=('rotated',1,'use_mh_covariance_matrix',1)
);
options_.TeX=1;
generate_trace_plots(1:2);
options_.TeX=1;

View File

@ -53,15 +53,15 @@ ln -s "$ROOT_DIRECTORY"/deps/mkoctfile64 /tmp/windeps/
# Go to source root directory
cd ..
common_meson_opts=(-Dbuildtype=release --cross-file scripts/windows-cross.ini)
common_meson_opts=(-Dbuildtype=release --cross-file windows/mingw-cross.ini)
# Create Windows 64-bit DLL binaries for MATLAB ≥ R2018b
meson setup --cross-file scripts/windows-cross-matlab.ini -Dmatlab_path=/tmp/windeps/matlab64/R2018b \
meson setup --cross-file windows/mingw-cross-matlab.ini -Dmatlab_path=/tmp/windeps/matlab64/R2018b \
"${common_meson_opts[@]}" build-win-matlab
meson compile -v -C build-win-matlab
# Create Windows DLL binaries for Octave/MinGW (64bit)
meson setup --cross-file scripts/windows-cross-octave.ini \
meson setup --cross-file windows/mingw-cross-octave.ini \
"${common_meson_opts[@]}" build-win-octave
meson compile -v -C build-win-octave

View File

@ -18,7 +18,7 @@ MATLAB64_VERSION = 20231122
## Build dependencies
# pacman -Ss mingw-w64-x86_64-boost
MINGW64_BOOST_VERSION = 1.83.0-1
MINGW64_BOOST_VERSION = 1.83.0-2
# pacman -Ss mingw-w64-x86_64-gsl
MINGW64_GSL_VERSION = 2.7.1-2
@ -40,11 +40,11 @@ MINGW64_LIBAEC_VERSION = 1.0.6-2
# Dependency of HDF5
# pacman -Ss mingw-w64-x86_64-openssl
MINGW64_OPENSSL_VERSION = 3.1.4-1
MINGW64_OPENSSL_VERSION = 3.2.0-1
# Dependency of HDF5
# pacman -Ss mingw-w64-x86_64-curl
MINGW64_CURL_VERSION = 8.4.0-1
MINGW64_CURL_VERSION = 8.5.0-1
# Dependency of curl (and of the MinGW compiler)
# pacman -Ss mingw-w64-x86_64-zstd
@ -81,22 +81,22 @@ MINGW64_LIBUNISTRING_VERSION = 1.1-1
## MinGW packages for the embedded compiler
# pacman -Ss mingw-w64-x86_64-gcc$
MINGW64_GCC_VERSION = 13.2.0-2
MINGW64_GCC_VERSION = 13.2.0-3
# pacman -Ss mingw-w64-x86_64-gmp
MINGW64_GMP_VERSION = 6.3.0-2
# pacman -Ss mingw-w64-x86_64-binutils
MINGW64_BINUTILS_VERSION = 2.41-2
MINGW64_BINUTILS_VERSION = 2.41-3
# pacman -Ss mingw-w64-x86_64-headers-git
MINGW64_HEADERS_VERSION = 11.0.0.r404.g3a137bd87-1
MINGW64_HEADERS_VERSION = 11.0.0.r442.ga27e7b27e-1
# pacman -Ss mingw-w64-x86_64-crt-git
MINGW64_CRT_VERSION = 11.0.0.r404.g3a137bd87-1
MINGW64_CRT_VERSION = 11.0.0.r442.ga27e7b27e-1
# pacman -Ss mingw-w64-x86_64-winpthreads-git
MINGW64_WINPTHREADS_VERSION = 11.0.0.r404.g3a137bd87-1
MINGW64_WINPTHREADS_VERSION = 11.0.0.r442.ga27e7b27e-1
# pacman -Ss mingw-w64-x86_64-isl
MINGW64_ISL_VERSION = 0.26-1

View File

@ -1,5 +1,5 @@
# Meson common cross file for targeting Windows from Linux
# To be included after windows-cross-{matlab,octave}.ini
# To be included after mingw-cross-{matlab,octave}.ini
[constants]
# For the architectural baseline, we follow MSYS2: