diff --git a/doc/manual/source/bibliography.rst b/doc/manual/source/bibliography.rst index 06f4faac0..1df9cb36d 100644 --- a/doc/manual/source/bibliography.rst +++ b/doc/manual/source/bibliography.rst @@ -47,6 +47,7 @@ Bibliography * Hansen, Lars P. (1982): “Large sample properties of generalized method of moments estimators,” Econometrica, 50(4), 1029–1054. * Hansen, Nikolaus and Stefan Kern (2004): “Evaluating the CMA Evolution Strategy on Multimodal Test Functions”. In: *Eighth International Conference on Parallel Problem Solving from Nature PPSN VIII*, Proceedings, Berlin: Springer, 282–291. * Harvey, Andrew C. and Garry D.A. Phillips (1979): “Maximum likelihood estimation of regression models with autoregressive-moving average disturbances,” *Biometrika*, 66(1), 49–58. +* Herbst, Edward and Schorfheide, Frank (2014): "Sequential monte-carlo sampling for DSGE models," *Journal of Applied Econometrics*, 29, 1073-1098. * Herbst, Edward (2015): “Using the “Chandrasekhar Recursions” for Likelihood Evaluation of DSGE Models,” *Computational Economics*, 45(4), 693–705. * Ireland, Peter (2004): “A Method for Taking Models to the Data,” *Journal of Economic Dynamics and Control*, 28, 1205–26. * Iskrev, Nikolay (2010): “Local identification in DSGE models,” *Journal of Monetary Economics*, 57(2), 189–202. diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst index 990fc9b22..a1445ec2d 100644 --- a/doc/manual/source/the-model-file.rst +++ b/doc/manual/source/the-model-file.rst @@ -7454,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 @@ -7466,141 +7478,173 @@ observed variables. Available options are: - .. _prop_distrib: + .. _prop_distrib: - ``'proposal_distribution'`` + ``'proposal_distribution'`` - Specifies the statistical distribution used for the - proposal density. + Specifies the statistical distribution used for the + proposal density. - ``'rand_multivariate_normal'`` + ``'rand_multivariate_normal'`` - Use a multivariate normal distribution. This is the default. + Use a multivariate normal distribution. This is the default. - ``'rand_multivariate_student'`` + ``'rand_multivariate_student'`` - Use a multivariate student distribution. + Use a multivariate student distribution. - ``'student_degrees_of_freedom'`` + ``'student_degrees_of_freedom'`` - Specifies the degrees of freedom to be used with the - multivariate student distribution. Default: ``3``. + Specifies the degrees of freedom to be used with the + multivariate student distribution. Default: ``3``. - .. _usemhcov: + .. _usemhcov: - ``'use_mh_covariance_matrix'`` + ``'use_mh_covariance_matrix'`` - Indicates to use the covariance matrix of the draws - from a previous MCMC run to define the covariance of - the proposal distribution. Requires the - :opt:`load_mh_file` option to be specified. Default: - ``0``. + Indicates to use the covariance matrix of the draws + from a previous MCMC run to define the covariance of + the proposal distribution. Requires the + :opt:`load_mh_file` option to be specified. Default: ``0``. - .. _scale-file: + .. _scale-file: - ``'scale_file'`` + ``'scale_file'`` - Provides the name of a ``_mh_scale.mat`` file storing - the tuned scale factor from a previous run of - ``mode_compute=6``. + Provides the name of a ``_mh_scale.mat`` file storing + the tuned scale factor from a previous run of + ``mode_compute=6``. - .. _savetmp: + .. _savetmp: - ``'save_tmp_file'`` + ``'save_tmp_file'`` - Save the MCMC draws into a ``_mh_tmp_blck`` file at the - refresh rate of the status bar instead of just saving - the draws when the current ``_mh*_blck`` file is - full. Default: ``0`` + Save the MCMC draws into a ``_mh_tmp_blck`` file at the + refresh rate of the status bar instead of just saving + the draws when the current ``_mh*_blck`` file is + full. Default: ``0`` ``'independent_metropolis_hastings'`` - Takes the same options as in the case of - ``random_walk_metropolis_hastings``. + Takes the same options as in the case of ``random_walk_metropolis_hastings``. ``'slice'`` - ``'rotated'`` + Available options are: - Triggers rotated slice iterations using a covariance - matrix from initial burn-in iterations. Requires either - ``use_mh_covariance_matrix`` or - ``slice_initialize_with_mode``. Default: ``0``. + ``'rotated'`` - ``'mode_files'`` + Triggers rotated slice iterations using a covariance + matrix from initial burn-in iterations. Requires either + ``use_mh_covariance_matrix`` or + ``slice_initialize_with_mode``. Default: ``0``. - For multimodal posteriors, provide the name of a file - containing a ``nparam`` by ``nmodes`` variable called - ``xparams`` storing the different modes. This array - must have one column vector per mode and the estimated - parameters along the row dimension. With this info, the - code will automatically trigger the ``rotated`` and - ``mode`` options. Default: ``[]``. + ``'mode_files'`` - ``'slice_initialize_with_mode'`` + For multimodal posteriors, provide the name of a file + containing a ``nparam`` by ``nmodes`` variable called + ``xparams`` storing the different modes. This array + must have one column vector per mode and the estimated + parameters along the row dimension. With this info, the + code will automatically trigger the ``rotated`` and + ``mode`` options. Default: ``[]``. - The default for slice is to set ``mode_compute=0`` and - start the chain(s) from a random location in the prior - space. This option first runs the mode-finder and then - starts the chain from the mode. Together with - ``rotated``, it will use the inverse Hessian from the - mode to perform rotated slice iterations. Default: - ``0``. + ``'slice_initialize_with_mode'`` - ``'initial_step_size'`` + The default for slice is to set ``mode_compute=0`` and + start the chain(s) from a random location in the prior + space. This option first runs the mode-finder and then + starts the chain from the mode. Together with + ``rotated``, it will use the inverse Hessian from the + mode to perform rotated slice iterations. Default: + ``0``. - Sets the initial size of the interval in the - stepping-out procedure as fraction of the prior - support, i.e. the size will be ``initial_step_size * - (UB-LB)``. ``initial_step_size`` must be a real number - in the interval ``[0,1]``. Default: ``0.8``. + ``'initial_step_size'`` - ``'use_mh_covariance_matrix'`` + Sets the initial size of the interval in the + stepping-out procedure as fraction of the prior + support, i.e. the size will be ``initial_step_size * + (UB-LB)``. ``initial_step_size`` must be a real number + in the interval ``[0,1]``. Default: ``0.8``. - See :ref:`use_mh_covariance_matrix `. Must be - used with ``'rotated'``. Default: ``0``. + ``'use_mh_covariance_matrix'`` - ``'save_tmp_file'`` + See :ref:`use_mh_covariance_matrix `. Must be + used with ``'rotated'``. Default: ``0``. - See :ref:`save_tmp_file `. Default: ``1``. + ``'save_tmp_file'`` + + See :ref:`save_tmp_file `. Default: ``1``. ``'tailored_random_block_metropolis_hastings'`` - ``'proposal_distribution'`` + Available options are: + + ``'proposal_distribution'`` - Specifies the statistical distribution used for the - proposal density. See :ref:`proposal_distribution `. + Specifies the statistical distribution used for the + proposal density. See :ref:`proposal_distribution `. - ``new_block_probability = DOUBLE`` + ``new_block_probability = DOUBLE`` - Specifies the probability of the next parameter - belonging to a new block when the random blocking in - the TaRB Metropolis-Hastings algorithm is - conducted. The higher this number, the smaller is the - average block size and the more random blocks are - formed during each parameter sweep. Default: ``0.25``. + Specifies the probability of the next parameter + belonging to a new block when the random blocking in + the TaRB Metropolis-Hastings algorithm is + conducted. The higher this number, the smaller is the + average block size and the more random blocks are + formed during each parameter sweep. Default: ``0.25``. - ``mode_compute = INTEGER`` + ``mode_compute = INTEGER`` - Specifies the mode-finder run in every iteration for - every block of the TaRB Metropolis-Hastings - algorithm. See :opt:`mode_compute `. Default: ``4``. + Specifies the mode-finder run in every iteration for + every block of the TaRB Metropolis-Hastings + algorithm. See :opt:`mode_compute `. Default: ``4``. - ``optim = (NAME, VALUE,...)`` + ``optim = (NAME, VALUE,...)`` - Specifies the options for the mode-finder used in the - TaRB Metropolis-Hastings algorithm. See :opt:`optim - `. + Specifies the options for the mode-finder used in the + TaRB Metropolis-Hastings algorithm. See :opt:`optim + `. - ``'scale_file'`` + ``'scale_file'`` - See :ref:`scale_file `.. + See :ref:`scale_file `.. - ``'save_tmp_file'`` + ``'save_tmp_file'`` - See :ref:`save_tmp_file `. Default: ``1``. + See :ref:`save_tmp_file `. Default: ``1``. + ``'hssmc'`` + + Available options are: + + ``'particles'`` + + Number of particles. Default value is: 20000. + + ``'steps'`` + + Number of weights :math:`\phi_i\in[0,1]` on the likelihood function used to define a sequence of tempered likelihoods. This parameter is denoted :math:`N_{\phi}` in *Herbst and Schorfheide (2014)*, and we have :math:`\phi_1=0` and :math:`\phi_{N_\phi}=1`. Default value is: 25. + + ``'lambda'`` + + Positive parameter controling the sequence of weights :math:`\phi_i`, Default value is: 2. Weights are defined by: + + .. math:: + + \phi_i = \left(\frac{i-1}{N_{\phi}-1}\right)^{\lambda} + + for :math:`i=1,\ldots,N_{\phi}`. Usually we set :math:`\lambda>1`, so that :math:`\Delta \phi_i = \phi_i-\phi_{i-1}` is increasing with :math:`i`. + + ``'target'`` + + Acceptance rate target. Default value is: .25. + + ``'scale'`` + + Scale parameter in the mutation step (on the proposal covariance matrix of the MH iteration). Default value is: .5. + .. option:: moments_varendo Triggers the computation of the posterior distribution of the diff --git a/matlab/@dprior/admissible.m b/matlab/@dprior/admissible.m new file mode 100644 index 000000000..c3a46c36b --- /dev/null +++ b/matlab/@dprior/admissible.m @@ -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 . + +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 diff --git a/matlab/tempered_likelihood.m b/matlab/@dprior/length.m similarity index 50% rename from matlab/tempered_likelihood.m rename to matlab/@dprior/length.m index 34df1f64f..6e25d6360 100644 --- a/matlab/tempered_likelihood.m +++ b/matlab/@dprior/length.m @@ -1,7 +1,17 @@ -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 n = length(o) -% Copyright © 2022 Dynare Team +% 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. % @@ -18,7 +28,4 @@ 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 . - 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 = length(o.lb); diff --git a/matlab/@dprior/subsref.m b/matlab/@dprior/subsref.m index 4e6cf3ea1..faea66d53 100644 --- a/matlab/@dprior/subsref.m +++ b/matlab/@dprior/subsref.m @@ -23,9 +23,9 @@ 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'}) + elseif ismember(S(1).subs, {'draw','length'}) p = feval(S(1).subs, o); - elseif ismember(S(1).subs, {'draws', 'density', 'densities', 'moments'}) + 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 diff --git a/matlab/GetAllPosteriorDraws.m b/matlab/GetAllPosteriorDraws.m index 23222b40d..dcb6e80dc 100644 --- a/matlab/GetAllPosteriorDraws.m +++ b/matlab/GetAllPosteriorDraws.m @@ -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 . -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 \ No newline at end of file +end diff --git a/matlab/GetPosteriorMeanVariance.m b/matlab/GetPosteriorMeanVariance.m index 5a2c44af2..37ef07411 100644 --- a/matlab/GetPosteriorMeanVariance.m +++ b/matlab/GetPosteriorMeanVariance.m @@ -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 . -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 diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m index bd29b1f61..e54964abd 100644 --- a/matlab/GetPosteriorParametersStatistics.m +++ b/matlab/GetPosteriorParametersStatistics.m @@ -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 . -%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); \ No newline at end of file +post_var = oo.posterior_variance.(type).(name); diff --git a/matlab/Herbst_Schorfheide_sampler.m b/matlab/Herbst_Schorfheide_sampler.m deleted file mode 100644 index f6881bc84..000000000 --- a/matlab/Herbst_Schorfheide_sampler.m +++ /dev/null @@ -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 . - - -% 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)= 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). + 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 \ No newline at end of file +end diff --git a/matlab/SMC_samplers_initialization.m b/matlab/SMC_samplers_initialization.m deleted file mode 100644 index f144ab66b..000000000 --- a/matlab/SMC_samplers_initialization.m +++ /dev/null @@ -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 . - -%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() - - diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m index a3f127bd9..3c81972df 100644 --- a/matlab/check_posterior_sampler_options.m +++ b/matlab/check_posterior_sampler_options.m @@ -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') diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m index e7f0c6a59..e6af987f6 100644 --- a/matlab/compute_mh_covariance_matrix.m +++ b/matlab/compute_mh_covariance_matrix.m @@ -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 _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 . -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'); \ No newline at end of file +mode = transpose(mode); diff --git a/matlab/compute_posterior_covariance_matrix.m b/matlab/compute_posterior_covariance_matrix.m new file mode 100644 index 000000000..7fe0ca235 --- /dev/null +++ b/matlab/compute_posterior_covariance_matrix.m @@ -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 _mh_mode.mat, hssmc_mode.mat or dsmh__mode.mat under //. +% +% 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 . + + +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 diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m index 25ca6b34b..1c1b1abb7 100644 --- a/matlab/convergence_diagnostics/mcmc_diagnostics.m +++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m @@ -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 diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m index cbaa775a2..b3a500c26 100644 --- a/matlab/default_option_values.m +++ b/matlab/default_option_values.m @@ -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 ; diff --git a/matlab/dprintf.m b/matlab/dprintf.m index c47ed1856..e122f72fa 100644 --- a/matlab/dprintf.m +++ b/matlab/dprintf.m @@ -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 . -disp(sprintf(str, varargin{:})); \ No newline at end of file +disp(sprintf(str, varargin{:})); diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index b0ee92251..13f0d5c84 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -51,6 +51,9 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ... '/discretionary_policy/' ; ... '/distributions/' ; ... '/ep/' ; ... + '/estimation/'; ... + '/estimation/smc/'; ... + '/estimation/resampler/'; ... '/gsa/' ; ... '/kalman/' ; ... '/kalman/likelihood' ; ... diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 990bd3775..40206b4d0 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -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 diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 4932e3042..76fa750ff 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -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. + +n= length(weights); + +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 diff --git a/matlab/DSMH_sampler.m b/matlab/estimation/smc/dsmh.m similarity index 93% rename from matlab/DSMH_sampler.m rename to matlab/estimation/smc/dsmh.m index 1c91c3a1e..5382639d5 100644 --- a/matlab/DSMH_sampler.m +++ b/matlab/estimation/smc/dsmh.m @@ -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 . +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 diff --git a/matlab/estimation/smc/hssmc.m b/matlab/estimation/smc/hssmc.m new file mode 100644 index 000000000..59a0bd704 --- /dev/null +++ b/matlab/estimation/smc/hssmc.m @@ -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 . + + 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. - 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 diff --git a/matlab/estimation/smc/ishssmc.m b/matlab/estimation/smc/ishssmc.m new file mode 100644 index 000000000..6d2820d44 --- /dev/null +++ b/matlab/estimation/smc/ishssmc.m @@ -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 . + +bool = false; +if isfield(options_, 'posterior_sampler_options') + if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'hssmc') + bool = true; + end +end diff --git a/matlab/estimation/smc/smc_samplers_initialization.m b/matlab/estimation/smc/smc_samplers_initialization.m new file mode 100644 index 000000000..fb27f251e --- /dev/null +++ b/matlab/estimation/smc/smc_samplers_initialization.m @@ -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 . + +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() diff --git a/matlab/estimation/smc/tempered_likelihood.m b/matlab/estimation/smc/tempered_likelihood.m new file mode 100644 index 000000000..bcb1e7bb6 --- /dev/null +++ b/matlab/estimation/smc/tempered_likelihood.m @@ -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 . + +logpostkernel = -postkernelfun(xparam); +logprior = Prior.density(xparam); +loglikelihood = logpostkernel-logprior; +tlogpostkernel = lambda*loglikelihood + logprior; diff --git a/matlab/fill_mh_mode.m b/matlab/fill_mh_mode.m index 22339a22d..0d2581ce1 100644 --- a/matlab/fill_mh_mode.m +++ b/matlab/fill_mh_mode.m @@ -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_..mode and oo_..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 \ No newline at end of file +end diff --git a/matlab/issmc.m b/matlab/issmc.m new file mode 100644 index 000000000..a79d444db --- /dev/null +++ b/matlab/issmc.m @@ -0,0 +1,20 @@ +function bool = issmc(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 . + +bool = ishssmc(options_) || isdsmh(options_); diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m index e9e203686..323857e8b 100644 --- a/matlab/marginal_density.m +++ b/matlab/marginal_density.m @@ -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 diff --git a/matlab/mh_autocorrelation_function.m b/matlab/mh_autocorrelation_function.m index cf10a8279..953b542ef 100644 --- a/matlab/mh_autocorrelation_function.m +++ b/matlab/mh_autocorrelation_function.m @@ -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); diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index a557cc652..595715a2b 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -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']); diff --git a/matlab/prior_bounds.m b/matlab/prior_bounds.m index 5db18c81f..78027c4c2 100644 --- a/matlab/prior_bounds.m +++ b/matlab/prior_bounds.m @@ -1,6 +1,5 @@ -function bounds = prior_bounds(bayestopt, priortrunc) +function bounds = prior_bounds(bayestopt_, priortrunc) -% function bounds = prior_bounds(bayestopt) % computes bounds for prior density. % % INPUTS @@ -31,11 +30,11 @@ if nargin<2, priortrunc = 0.0; end 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; +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)); @@ -59,12 +58,12 @@ for i=1:length(p6) bounds.ub(i) = gaminv(1.0-priortrunc, p6(i), p7(i))+p3(i); end case 3 - if prior_trunc == 0 + 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 priortrunc==0 @@ -99,6 +98,6 @@ for i=1:length(p6) 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 diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 58f428c5a..6971f9465 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -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 \ No newline at end of file +end diff --git a/matlab/trace_plot.m b/matlab/trace_plot.m index d4cf5e1f1..9de1860b9 100644 --- a/matlab/trace_plot.m +++ b/matlab/trace_plot.m @@ -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))]; diff --git a/meson.build b/meson.build index 42d7017be..bf00091ba 100644 --- a/meson.build +++ b/meson.build @@ -805,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', diff --git a/tests/estimation/dsmc/fs2000.mod b/tests/estimation/dsmc/fs2000.mod new file mode 100644 index 000000000..647e50f95 --- /dev/null +++ b/tests/estimation/dsmc/fs2000.mod @@ -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'); diff --git a/tests/estimation/hssmc/fs2000.mod b/tests/estimation/hssmc/fs2000.mod new file mode 100644 index 000000000..d843c9e98 --- /dev/null +++ b/tests/estimation/hssmc/fs2000.mod @@ -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)); diff --git a/tests/estimation/slice/fs2000_slice.mod b/tests/estimation/slice/fs2000_slice.mod index f50a1ec4e..ecaf831e3 100644 --- a/tests/estimation/slice/fs2000_slice.mod +++ b/tests/estimation/slice/fs2000_slice.mod @@ -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; \ No newline at end of file +// 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;