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

Merge branch 'dynare-irf_matching' 

Ref. !2159
kalman-mex
Sébastien Villemot 2023-09-08 09:45:51 +02:00
parent 6551b88fc8 2e742462b3
commit cb0f0e6701
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
42 changed files with 2083 additions and 1367 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
doc/parallel/parallel.tex
View File

@ -189,7 +189,7 @@ We have considered the following DYNARE components suitable to be parallelized u
\item the Random Walk- (and the analogous Independent-)-Metropolis-Hastings algorithm with multiple chains: the different chains are completely independent and do not require any communication between them, so it can be executed on different cores/CPUs/Computer Network easily; \item the Random Walk- (and the analogous Independent-)-Metropolis-Hastings algorithm with multiple chains: the different chains are completely independent and do not require any communication between them, so it can be executed on different cores/CPUs/Computer Network easily;
\item a number of procedures performed after the completion of Metropolis, that use the posterior MC sample: \item a number of procedures performed after the completion of Metropolis, that use the posterior MC sample:
\begin{enumerate} \begin{enumerate}
\item the diagnostic tests for the convergence of the Markov Chain \\(\texttt{McMCDiagnostics.m}); \item the diagnostic tests for the convergence of the Markov Chain \\(\texttt{mcmc\_diagnostics.m});
\item the function that computes posterior IRF's (\texttt{posteriorIRF.m}). \item the function that computes posterior IRF's (\texttt{posteriorIRF.m}).
\item the function that computes posterior statistics for filtered and smoothed variables, forecasts, smoothed shocks, etc.. \\ (\verb"prior_posterior_statistics.m"). \item the function that computes posterior statistics for filtered and smoothed variables, forecasts, smoothed shocks, etc.. \\ (\verb"prior_posterior_statistics.m").
\item the utility function that loads matrices of results and produces plots for posterior statistics (\texttt{pm3.m}). \item the utility function that loads matrices of results and produces plots for posterior statistics (\texttt{pm3.m}).
@ -725,8 +725,8 @@ So far, we have parallelized the following functions, by selecting the most comp
\verb"independent_metropolis_hastings.m", \\ \verb"independent_metropolis_hastings.m", \\
\verb"independent_metropolis_hastings_core.m"; \verb"independent_metropolis_hastings_core.m";
\item the cycle looping over estimated parameters computing univariate diagnostics:\\ \item the cycle looping over estimated parameters computing univariate diagnostics:\\
\verb"McMCDiagnostics.m", \\ \verb"mcmc_diagnostics.m", \\
\verb"McMCDiagnostics_core.m"; \verb"mcmc_diagnostics_core.m";
\item the Monte Carlo cycle looping over posterior parameter subdraws performing the IRF simulations (\verb"<*>_core1") and the cycle looping over exogenous shocks plotting IRF's charts (\verb"<*>_core2"):\\ \item the Monte Carlo cycle looping over posterior parameter subdraws performing the IRF simulations (\verb"<*>_core1") and the cycle looping over exogenous shocks plotting IRF's charts (\verb"<*>_core2"):\\
\verb"posteriorIRF.m", \\\verb"posteriorIRF_core1.m", \verb"posteriorIRF_core2.m"; \verb"posteriorIRF.m", \\\verb"posteriorIRF_core1.m", \verb"posteriorIRF_core2.m";
\item the Monte Carlo cycle looping over posterior parameter subdraws, that computes filtered, smoothed, forecasted variables and shocks:\\ \item the Monte Carlo cycle looping over posterior parameter subdraws, that computes filtered, smoothed, forecasted variables and shocks:\\

67
matlab/+mom/Jtest.m Normal file
View File

@ -0,0 +1,67 @@
function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
% function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
% -------------------------------------------------------------------------
% Computes the J-test statistic and p-value for a GMM/SMM estimation
% =========================================================================
% INPUTS
% xparam: [vector] estimated parameter vector
% objective_function: [function handle] objective function
% Woptflag: [logical] flag if optimal weighting matrix has already been computed
% oo_: [struct] results
% options_mom_: [struct] options
% bayestopt_: [struct] information on priors
% Bounds: [struct] bounds on parameters
% estim_params_: [struct] information on estimated parameters
% M_: [struct] information on the model
% nobs: [scalar] number of observations
% -------------------------------------------------------------------------
% OUTPUT
% oo_: [struct] updated results
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o mom.objective_function
% o mom.optimal_weighting_matrix
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if options_mom_.mom.mom_nbr > length(xparam)
% Get optimal weighting matrix for J test, if necessary
if ~Woptflag
W_opt = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
oo_J = oo_;
oo_J.mom.Sw = chol(W_opt);
fval = feval(objective_function, xparam, Bounds, oo_J, estim_params_, M_, options_mom_, bayestopt_);
else
fval = oo_.mom.Q;
end
% Compute J statistic
if strcmp(options_mom_.mom.mom_method,'SMM')
Variance_correction_factor = options_mom_.mom.variance_correction_factor;
elseif strcmp(options_mom_.mom.mom_method,'GMM')
Variance_correction_factor = 1;
end
oo_.mom.J_test.j_stat = nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam);
oo_.mom.J_test.p_val = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
fprintf('\nValue of J-test statistic: %f\n',oo_.mom.J_test.j_stat);
fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val);
end

270
matlab/+mom/default_option_mom_values.m Normal file
View File

@ -0,0 +1,270 @@
function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
% function options_mom_ = default_option_mom_values(options_mom_, options_, dname, doBayesianEstimation)
% Returns structure containing the options for method_of_moments command
% options_mom_ is local and contains default and user-specified values for
% all settings needed for the method of moments estimation. Some options,
% though, are set by the preprocessor into options_ and we copy these over.
% The idea is to be independent of options_ and have full control of the
% estimation instead of possibly having to deal with options chosen somewhere
% else in the mod file.
% =========================================================================
% INPUTS
% o options_mom_: [structure] information about all (user-specified and updated) settings used in estimation (options_mom_)
% o options_: [structure] information on global options
% o dname: [string] name of directory to store results
% o doBayesianEstimation [boolean] indicator whether we do Bayesian estimation
% -------------------------------------------------------------------------
% OUTPUTS
% o oo_: [structure] storage for results (oo_)
% o options_mom_: [structure] information about all (user-specified and updated) settings used in estimation (options_mom_)
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o set_default_option
% o user_has_matlab_license
% o user_has_octave_forge_package
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
mom_method = options_mom_.mom.mom_method; % this is a required option
% -------------------------------------------------------------------------
% LIMITATIONS
% -------------------------------------------------------------------------
if options_.logged_steady_state || options_.loglinear
error('method_of_moments: The loglinear option is not supported. Please append the required logged variables as auxiliary equations.')
else
options_mom_.logged_steady_state = 0;
options_mom_.loglinear = false;
end
options_mom_.hessian.use_penalized_objective = false; % penalized objective not yet
% options related to variable declarations
if isfield(options_,'trend_coeffs')
error('method_of_moments: %s does not allow for trend in data',mom_method)
end
% options related to endogenous prior restrictions are not supported
if ~isempty(options_.endogenous_prior_restrictions.irf) && ~isempty(options_.endogenous_prior_restrictions.moment)
fprintf('method_of_moments: Endogenous prior restrictions are not supported yet and will be skipped.\n')
end
options_mom_.endogenous_prior_restrictions.irf = {};
options_mom_.endogenous_prior_restrictions.moment = {};
options_mom_.mom.analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; % these are currently supported optimizers that are able to use the analytical_jacobian option
% -------------------------------------------------------------------------
% OPTIONS POSSIBLY SET BY THE USER
% -------------------------------------------------------------------------
% common settings
options_mom_ = set_default_option(options_mom_,'dirname',dname); % specify directory in which to store estimation output
options_mom_ = set_default_option(options_mom_,'graph_format','eps'); % specify the file format(s) for graphs saved to disk
options_mom_ = set_default_option(options_mom_,'nodisplay',false); % do not display the graphs, but still save them to disk
options_mom_ = set_default_option(options_mom_,'nograph',false); % do not create graphs (which implies that they are not saved to the disk nor displayed)
options_mom_ = set_default_option(options_mom_,'noprint',false); % do not print output to console
options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics
options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results
%options_mom_ = set_default_option(options_mom_,'verbosity',false); %
if doBayesianEstimation
options_mom_ = set_default_option(options_mom_,'plot_priors',true); % control plotting of priors
options_mom_ = set_default_option(options_mom_,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters
end
% specific method_of_moments settings
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
options_mom_.mom = set_default_option(options_mom_.mom,'bartlett_kernel_lag',20); % bandwith in optimal weighting matrix
options_mom_.mom = set_default_option(options_mom_.mom,'penalized_estimator',false); % include deviation from prior mean as additional moment restriction and use prior precision as weights
options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for numerical computation of standard errors
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix in objective function
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'OPTIMAL'}); % weighting matrix in moments distance objective function at each iteration of estimation;
% possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation.
end
if strcmp(mom_method,'SMM')
options_mom_.mom = set_default_option(options_mom_.mom,'burnin',500); % number of periods dropped at beginning of simulation
options_mom_.mom = set_default_option(options_mom_.mom,'bounded_shock_support',false); % trim shocks in simulation to +- 2 stdev
options_mom_.mom = set_default_option(options_mom_.mom,'seed',24051986); % seed used in simulations
options_mom_.mom = set_default_option(options_mom_.mom,'simulation_multiple',7); % multiple of the data length used for simulation
end
if strcmp(mom_method,'GMM')
options_mom_.mom = set_default_option(options_mom_.mom,'analytic_standard_errors',false); % compute standard errors numerically (0) or analytically (1). Analytical derivatives are only available for GMM.
end
% data related options
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
options_mom_ = set_default_option(options_mom_,'first_obs',1); % number of first observation
options_mom_ = set_default_option(options_mom_,'logdata',false); % if data is already in logs
options_mom_ = set_default_option(options_mom_,'nobs',NaN); % number of observations
options_mom_ = set_default_option(options_mom_,'prefilter',false); % demean each data series by its empirical mean and use centered moments
options_mom_ = set_default_option(options_mom_,'xls_sheet',1); % name of sheet with data in Excel, Octave does not support the empty string, rather use first sheet
options_mom_ = set_default_option(options_mom_,'xls_range',''); % range of data in Excel sheet
end
% optimization related
if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox'))
if strcmp(mom_method,'GMM') || strcmp(mom_method,'SMM')
options_mom_ = set_default_option(options_mom_,'mode_compute',13); % specifies lsqnonlin as default optimizer for minimization
end
else
options_mom_ = set_default_option(options_mom_,'mode_compute',4); % specifies csminwel as fallback default option for minimization
end
options_mom_ = set_default_option(options_mom_,'additional_optimizer_steps',[]); % vector of additional mode-finders run after mode_compute
options_mom_ = set_default_option(options_mom_,'optim_opt',[]); % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
options_mom_ = set_default_option(options_mom_,'silent_optimizer',false); % run minimization of moments distance silently without displaying results or saving files in between
options_mom_ = set_default_option(options_mom_,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
options_mom_.mom = set_default_option(options_mom_.mom,'analytic_jacobian',false); % use analytic Jacobian in optimization, only available for GMM and gradient-based optimizers
options_mom_.optimizer_vec = [options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)];
% perturbation related
options_mom_ = set_default_option(options_mom_,'order',1); % order of Taylor approximation in perturbation
options_mom_ = set_default_option(options_mom_,'pruning',false); % use pruned state space system at order>1
options_mom_ = set_default_option(options_mom_,'aim_solver',false); % use AIM algorithm to compute perturbation approximation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'k_order_solver',false); % use k_order_perturbation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false); % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false); % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6); % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
% if there are no unit roots one can use 1.0 (or slightly below) which we set as default; if they are possible, you may have have multiple unit roots and the accuracy decreases when computing the eigenvalues in lyapunov_symm
% Note that unit roots are only possible at first-order, at higher order we set it to 1 in pruned_state_space_system and focus only on stationary observables.
options_mom_ = set_default_option(options_mom_,'qz_zero_threshold',1e-6); % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
options_mom_ = set_default_option(options_mom_,'schur_vec_tol',1e-11); % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
% numerical algorithms
options_mom_ = set_default_option(options_mom_,'lyapunov_db',false); % doubling algorithm (disclyap_fast) to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_fp',false); % fixed-point algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false); % square-root-solver (dlyapchol) algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15); % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10); % convergence criterion used in the fixed point Lyapunov solver
options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm
options_mom_ = set_default_option(options_mom_,'sylvester_fp',false); % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12); % convergence criterion used in the fixed point Sylvester solver
% mode check plot
options_mom_.mode_check.nolik = false; % we don't do likelihood (also this initializes mode_check substructure)
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed minimum for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'neighbourhood_size',.5); % width of the window around the computed minimum to be displayed on the diagnostic plots. This width is expressed in percentage deviation. The Inf value is allowed, and will trigger a plot over the entire domain
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'symmetric_plots',true); % ensure that the check plots are symmetric around the minimum. A value of 0 allows to have asymmetric plots, which can be useful if the minimum is close to a domain boundary, or in conjunction with neighbourhood_size = Inf when the domain is not the entire real line
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'number_of_points',20); % number of points around the minimum where the target function is evaluated (for each parameter)
% -------------------------------------------------------------------------
% OPTIONS THAT NEED TO BE CARRIED OVER (E.G. SET BY THE PREPROCESSOR)
% -------------------------------------------------------------------------
% related to VAROBS block
options_mom_.varobs = options_.varobs; % observable variables in order they are declared in varobs
options_mom_.varobs_id = options_.varobs_id; % index for observable variables in M_.endo_names
options_mom_.obs_nbr = length(options_mom_.varobs); % number of observed variables
% related to call of dynare
options_mom_.console_mode = options_.console_mode;
options_mom_.parallel = options_.parallel;
options_mom_.parallel_info = options_.parallel_info;
% related to estimated_params and estimated_params_init blocks
options_mom_.use_calibration_initialization = options_.use_calibration_initialization;
% related to model block
options_mom_.linear = options_.linear;
options_mom_.use_dll = options_.use_dll;
options_mom_.block = options_.block;
options_mom_.bytecode = options_.bytecode;
% related to steady-state computations
options_mom_.homotopy_force_continue = options_.homotopy_force_continue;
options_mom_.homotopy_mode = options_.homotopy_mode;
options_mom_.homotopy_steps = options_.homotopy_steps;
options_mom_.markowitz = options_.markowitz;
options_mom_.solve_algo = options_.solve_algo;
options_mom_.solve_tolf = options_.solve_tolf;
options_mom_.solve_tolx = options_.solve_tolx;
options_mom_.steady = options_.steady;
options_mom_.steadystate = options_.steadystate;
options_mom_.steadystate_flag = options_.steadystate_flag;
%options_mom_.steadystate_partial
options_mom_.threads = options_.threads; % needed by resol
options_mom_.debug = options_.debug; % debug option needed by some functions, e.g. check_plot
% random numbers
options_mom_.DynareRandomStreams.seed = options_.DynareRandomStreams.seed;
options_mom_.DynareRandomStreams.algo = options_.DynareRandomStreams.algo;
% dataset_ related
options_mom_.dataset = options_.dataset;
options_mom_.initial_period = options_.initial_period;
% optimization related
if any(cellfun(@(x) isnumeric(x) && any(x == 2), options_mom_.optimizer_vec)) % simulated annealing (mode_compute=2)
options_mom_.saopt = options_.saopt;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 4), options_mom_.optimizer_vec)) % csminwel (mode_compute=4)
options_mom_.csminwel = options_.csminwel;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 5), options_mom_.optimizer_vec)) % newrat (mode_compute=5)
options_mom_.newrat = options_.newrat;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 6), options_mom_.optimizer_vec)) % gmhmaxlik (mode_compute=6)
options_mom_.gmhmaxlik = options_.gmhmaxlik;
options_mom_.mh_jscale = options_.mh_jscale;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 8), options_mom_.optimizer_vec)) % simplex variation on Nelder Mead algorithm (mode_compute=8)
options_mom_.simplex = options_.simplex;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 9), options_mom_.optimizer_vec)) % cmaes (mode_compute=9)
options_mom_.cmaes = options_.cmaes;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 10), options_mom_.optimizer_vec)) % simpsa (mode_compute=10)
options_mom_.simpsa = options_.simpsa;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 12), options_mom_.optimizer_vec)) % particleswarm (mode_compute=12)
options_mom_.particleswarm = options_.particleswarm;
end
if any(cellfun(@(x) isnumeric(x) && any(x == 101), options_mom_.optimizer_vec)) % solveopt (mode_compute=101)
options_mom_.solveopt = options_.solveopt;
end
if any(cellfun(@(x) isnumeric(x) && (any(x == 4) || any(x == 5)), options_mom_.optimizer_vec)) % used by csminwel and newrat
options_mom_.gradient_method = options_.gradient_method;
options_mom_.gradient_epsilon = options_.gradient_epsilon;
end
options_mom_.gstep = options_.gstep; % needed by hessian.m
options_mom_.trust_region_initial_step_bound_factor = options_.trust_region_initial_step_bound_factor; % used in dynare_solve for trust_region
% other
options_mom_.MaxNumberOfBytes = options_.MaxNumberOfBytes;
%options_mom_.MaximumNumberOfMegaBytes = options_.MaximumNumberOfMegaBytes;
% -------------------------------------------------------------------------
% DEFAULT VALUES
% -------------------------------------------------------------------------
options_mom_.analytic_derivation = 0;
options_mom_.analytic_derivation_mode = 0; % needed by get_perturbation_params_derivs.m, ie use efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m
options_mom_.figures = options_.figures; % needed by plot_priors.m
options_mom_.ramsey_policy = false; % needed by evaluate_steady_state
options_mom_.risky_steadystate = false; % needed by resol
options_mom_.jacobian_flag = true; % needed by dynare_solve

74
matlab/+mom/display_comparison_moments.m Normal file
View File

@ -0,0 +1,74 @@
function display_comparison_moments(M_, options_mom_, data_moments, model_moments)
% function display_comparison_moments(M_, options_mom_, data_moments, model_moments)
% -------------------------------------------------------------------------
% Displays and saves to disk the comparison of the data moments and the model moments
% =========================================================================
% INPUTS
% M_: [structure] model information
% options_mom_: [structure] method of moments options
% data_moments: [vector] data moments
% model_moments: [vector] model moments
% -------------------------------------------------------------------------
% OUTPUT
% No output, just displays and saves to disk the comparison of the data moments and the model moments
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o dyn_latex_table
% o dyntable
% o cellofchararraymaxlength
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
titl = ['Comparison of matched data moments and model moments (',options_mom_.mom.mom_method,')'];
headers = {'Moment','Data','Model'};
for jm = 1:size(M_.matched_moments,1)
lables_tmp = 'E[';
lables_tmp_tex = 'E \left[ ';
for jvar = 1:length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp M_.endo_names{M_.matched_moments{jm,1}(jvar)}];
lables_tmp_tex = [lables_tmp_tex, '{', M_.endo_names_tex{M_.matched_moments{jm,1}(jvar)}, '}'];
if M_.matched_moments{jm,2}(jvar) ~= 0
lables_tmp = [lables_tmp, '(', num2str(M_.matched_moments{jm,2}(jvar)), ')'];
lables_tmp_tex = [lables_tmp_tex, '_{t', num2str(M_.matched_moments{jm,2}(jvar)), '}'];
else
lables_tmp_tex = [lables_tmp_tex, '_{t}'];
end
if M_.matched_moments{jm,3}(jvar) > 1
lables_tmp = [lables_tmp, '^', num2str(M_.matched_moments{jm,3}(jvar))];
lables_tmp_tex = [lables_tmp_tex, '^{', num2str(M_.matched_moments{jm,3}(jvar)) '}'];
end
if jvar == length(M_.matched_moments{jm,1})
lables_tmp = [lables_tmp, ']'];
lables_tmp_tex = [lables_tmp_tex, ' \right]'];
else
lables_tmp = [lables_tmp, '*'];
lables_tmp_tex = [lables_tmp_tex, ' \times '];
end
end
labels{jm,1} = lables_tmp;
labels_TeX{jm,1} = lables_tmp_tex;
end
data_mat = [data_moments model_moments];
dyntable(options_mom_, titl, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
if options_mom_.TeX
dyn_latex_table(M_, options_mom_, titl, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7);
end

17
matlab/+mom/data_moments.m → matlab/+mom/get_data_moments.m
View File

@ -1,5 +1,5 @@
function [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_) function [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
% [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_) % [dataMoments, m_data] = get_data_moments(data, oo_, matched_moments_, options_mom_)
% This function computes the user-selected empirical moments from data % This function computes the user-selected empirical moments from data
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
@ -13,10 +13,10 @@ function [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, optio
% o m_data [T x numMom] selected empirical moments at each point in time % o m_data [T x numMom] selected empirical moments at each point in time
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o mom.run.m % o mom.run
% o mom.objective_function.m % o mom.objective_function
% ========================================================================= % =========================================================================
% Copyright © 2020-2021 Dynare Team % Copyright © 2020-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -49,7 +49,7 @@ for jm = 1:options_mom_.mom.mom_nbr
leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers
powers = matched_moments_{jm,3}; powers = matched_moments_{jm,3};
for jv = 1:length(vars) for jv = 1:length(vars)
jvar = (oo_.dr.obs_var == vars(jv)); jvar = (oo_.mom.obs_var == vars(jv));
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv); y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
if jv==1 if jv==1
@ -66,7 +66,4 @@ for jm = 1:options_mom_.mom.mom_nbr
end end
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1); m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
m_data(:,jm) = m_data_tmp; m_data(:,jm) = m_data_tmp;
end end
end %function end

80
matlab/+mom/matched_moments_block.m Normal file
View File

@ -0,0 +1,80 @@
function matched_moments = matched_moments_block(matched_moments, mom_method)
% function matched_moments = matched_moments_block(matched_moments, mom_method)
% -------------------------------------------------------------------------
% Checks and transforms matched_moments bock for further use in the estimation
% =========================================================================
% INPUTS
% matched_moments: [cell array] original matched_moments block
% mom_method: [string] method of moments method (GMM or SMM)
% -------------------------------------------------------------------------
% OUTPUT
% matched_moments: [cell array] transformed matched_moments block
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
matched_moments_orig = matched_moments;
% higher-order product moments not supported yet for GMM
if strcmp(mom_method, 'GMM') && any(cellfun(@sum,matched_moments(:,3))> 2)
error('method_of_moments: GMM does not yet support product moments higher than 2. Change your ''matched_moments'' block!');
end
% check for duplicate moment conditions
for jm = 1:size(matched_moments,1)
% expand powers to vector of ones
if any(matched_moments{jm,3}>1)
tmp1=[]; tmp2=[]; tmp3=[];
for jjm=1:length(matched_moments{jm,3})
tmp1 = [tmp1 repmat(matched_moments{jm,1}(jjm),[1 matched_moments{jm,3}(jjm)]) ];
tmp2 = [tmp2 repmat(matched_moments{jm,2}(jjm),[1 matched_moments{jm,3}(jjm)]) ];
tmp3 = [tmp3 repmat(1,[1 matched_moments{jm,3}(jjm)]) ];
end
matched_moments{jm,1} = tmp1;
matched_moments{jm,2} = tmp2;
matched_moments{jm,3} = tmp3;
end
% shift time structure to focus only on lags
matched_moments{jm,2} = matched_moments{jm,2} - max(matched_moments{jm,2});
% sort such that t=0 variable comes first
[matched_moments{jm,2},idx_sort] = sort(matched_moments{jm,2},'descend');
matched_moments{jm,1} = matched_moments{jm,1}(idx_sort);
matched_moments{jm,3} = matched_moments{jm,3}(idx_sort);
end
% find duplicate rows in cell array by making groups according to powers as we can then use cell2mat for the unique function
powers = cellfun(@sum,matched_moments(:,3))';
UniqueMomIdx = [];
for jpow = unique(powers)
idx1 = find(powers==jpow);
[~,idx2] = unique(cell2mat(matched_moments(idx1,:)),'rows');
UniqueMomIdx = [UniqueMomIdx idx1(idx2)];
end
% remove duplicate elements
DuplicateMoms = setdiff(1:size(matched_moments_orig,1),UniqueMomIdx);
if ~isempty(DuplicateMoms)
fprintf('Duplicate declared moments found and removed in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms))
fprintf('Dynare will continue with remaining moment conditions\n');
end
if strcmp(mom_method, 'SMM')
% for SMM: keep the original structure, but get rid of duplicate moments
matched_moments = matched_moments_orig(sort(UniqueMomIdx),:);
elseif strcmp(mom_method, 'GMM')
% for GMM we use the transformed matched_moments structure
matched_moments = matched_moments(sort(UniqueMomIdx),:);
end

129
matlab/+mom/mode_compute_gmm_smm.m Normal file
View File

@ -0,0 +1,129 @@
function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
% function [xparam1, oo_, Woptflag] = mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds)
% -------------------------------------------------------------------------
% Iterated method of moments for GMM and SMM, computes the minimum of the
% objective function (distance between data moments and model moments)
% for a sequence of optimizers and GMM/SMM iterations with different
% weighting matrices.
% =========================================================================
% INPUTS
% xparam0: [vector] vector of initialized parameters
% objective_function: [func handle] name of the objective function
% oo_: [structure] results
% M_: [structure] model information
% options_mom_: [structure] options
% estim_params_: [structure] information on estimated parameters
% bayestopt_: [structure] information on priors
% Bounds: [structure] bounds for optimization
% -------------------------------------------------------------------------
% OUTPUT
% xparam1: [vector] mode of objective function
% oo_: [structure] updated results
% Woptflag: [logical] true if optimal weighting matrix was computed
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o mom.optimal_weighting_matrix
% o mom.display_estimation_results_table
% o dynare_minimize_objective
% o mom.objective_function
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',options_mom_.mom.weighting_matrix)) || any(strcmpi('optimal',options_mom_.mom.weighting_matrix)))
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
end
for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
weighting_matrix = eye(options_mom_.mom.mom_nbr);
case 'diagonal'
fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) ));
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = diag(diag( mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) ));
end
case 'optimal'
fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
Woptflag = true;
end
otherwise % user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
try
load(options_mom_.mom.weighting_matrix{stage_iter},'weighting_matrix')
catch
error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',options_mom_.mom.weighting_matrix{stage_iter},'.mat !'])
end
[nrow, ncol] = size(weighting_matrix);
if ~isequal(nrow,ncol) || ~isequal(nrow,length(oo_.mom.data_moments)) %check if square and right size
error(['method_of_moments: ''weighting_matrix'' must be square and have ',num2str(length(oo_.mom.data_moments)),' rows and columns!'])
end
end
try % check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified ''weighting_matrix'' is not positive definite. Check whether your model implies stochastic singularity!')
end
for optim_iter = 1:length(options_mom_.optimizer_vec)
options_mom_.current_optimizer = options_mom_.optimizer_vec{optim_iter};
if options_mom_.optimizer_vec{optim_iter} == 0
xparam1 = xparam0; % no minimization, evaluate objective at current values
fval = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
else
if options_mom_.optimizer_vec{optim_iter} == 13
options_mom_.mom.vector_output = true;
else
options_mom_.mom.vector_output = false;
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(options_mom_.optimizer_vec{optim_iter},options_mom_.mom.analytic_jacobian_optimizers) %do this only for gradient-based optimizers
options_mom_.mom.compute_derivs = true;
else
options_mom_.mom.compute_derivs = false;
end
[xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_.name, bayestopt_, [],...
Bounds, oo_, estim_params_, M_, options_mom_);
if options_mom_.mom.vector_output
fval = fval'*fval;
end
end
fprintf('\nStage %d Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',stage_iter,optim_iter,fval)
if options_mom_.mom.verbose
oo_.mom = display_estimation_results_table(xparam1,NaN(size(xparam1)),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,sprintf('%s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter),sprintf('verbose_%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter));
end
xparam0 = xparam1;
end
options_mom_.vector_output = false;
[~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % get oo_.mom.model_moments for iterated GMM/SMM to compute optimal weighting matrix
end

257
matlab/+mom/objective_function.m
View File

@ -1,42 +1,43 @@
function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_) function [fval, info, exit_flag, df, junkHessian, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_) % [fval, info, exit_flag, df, junk1, oo_, M_] = objective_function(xparam, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function evaluates the objective function for GMM/SMM estimation % This function evaluates the objective function for method of moments estimation
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% o xparam1: current value of estimated parameters as returned by set_prior() % o xparam: [vector] current value of estimated parameters as returned by set_prior()
% o Bounds: structure containing parameter bounds % o Bounds: [structure] containing parameter bounds
% o oo_: structure for results % o oo_: [structure] for results
% o estim_params_: structure describing the estimated_parameters % o estim_params_: [structure] describing the estimated_parameters
% o M_ structure describing the model % o M_ [structure] describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_) % o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% o bayestopt_: [structure] information about the prior
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% OUTPUTS % OUTPUTS
% o fval: value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly) % o fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
% o info: vector storing error code and penalty % o info: [vector] information on error codes and penalties
% o exit_flag: 0 if error, 1 if no error % o exit_flag: [double] flag for exit status (0 if error, 1 if no error)
% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only) % o df: [matrix] analytical jacobian of the moment difference (wrt paramters), currently for GMM only
% o junk1: empty matrix required for optimizer interface (Hessian would go here) % o junkHessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here)
% o oo_: structure containing the results with the following updated fields: % o oo_: [structure] results with the following updated fields:
% - mom.model_moments [numMom x 1] vector with model moments % - oo_.mom.model_moments: [vector] model moments
% - mom.Q value of the quadratic form of the moment difference % - oo_.mom.Q: [double] value of the quadratic form of the moment difference
% - mom.model_moments_params_derivs % - oo_.mom.model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
% [numMom x numParams] Jacobian matrix of derivatives of model_moments with respect to estimated parameters % o M_: [structure] updated model structure
% (only for GMM with analytical derivatives)
% o M_: Matlab's structure describing the model
% o options_mom_: structure information about all settings (specified by the user, preprocessor, and taken from global options_)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% o mom.run.m % o mom.run
% o dynare_minimize_objective.m % o dynare_minimize_objective
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% o check_bounds_and_definiteness_estimation % o check_bounds_and_definiteness_estimation
% o pruned_state_space_system % o get_perturbation_params_derivs
% o resol % o mom.get_data_moments
% o set_all_parameters % o pruned_state_space_system
% o resol
% o set_all_parameters
% o simult_
% ========================================================================= % =========================================================================
% Copyright © 2020-2021 Dynare Team % Copyright © 2020-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -52,65 +53,65 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_f
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Author(s):
% o Willi Mutschler (willi@mutschler.eu)
% o Johannes Pfeifer (jpfeifer@uni-koeln.de)
% ========================================================================= % =========================================================================
%% TO DO
% check the info values and make use of meaningful penalties
% how do we do the penalty for the prior??
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
% 0. Initialization of the returned variables and others... % Initialization of the returned variables and others...
%------------------------------------------------------------------------------ %------------------------------------------------------------------------------
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian junkHessian = [];
if options_mom_.vector_output == 1 df = []; % required to be empty by e.g. newrat
if options_mom_.mom.penalized_estimator if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1)); if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
if options_mom_.mom.vector_output == 1
if options_mom_.mom.penalized_estimator
df = nan(size(oo_.mom.data_moments,1)+length(xparam),length(xparam));
else
df = nan(size(oo_.mom.data_moments,1),length(xparam));
end
else else
df = nan(size(oo_.mom.data_moments,1),length(xparam1)); df = nan(1,length(xparam));
end end
else
df = nan(1,length(xparam1));
end end
else
df=[]; %required to be empty by e.g. newrat
end end
junk1 = [];
junk2 = [];
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties % Get the structural parameters and define penalties
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
% Ensure that xparam1 is a column vector; particleswarm.m requires this. % Ensure that xparam1 is a column vector; particleswarm.m requires this.
xparam1 = xparam1(:); xparam = xparam(:);
M_ = set_all_parameters(xparam, estim_params_, M_);
M_ = set_all_parameters(xparam1, estim_params_, M_); [fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, Bounds);
[fval,info,exit_flag]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, Bounds);
if info(1) if info(1)
if options_mom_.vector_output == 1 % lsqnonlin requires vector output if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number; fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
end end
return return
end end
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
% 2. call resol to compute steady state and model solution % Call resol to compute steady state and model solution
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
% Compute linear approximation around the deterministic steady state % Compute linear approximation around the deterministic steady state
[dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_); [dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_);
% Return, with endogenous penalty when possible, if resol issues an error code % Return, with endogenous penalty when possible, if resol issues an error code
if info(1) if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||... if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ... info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86 info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
%meaningful second entry of output that can be used % meaningful second entry of output that can be used
fval = Inf; fval = Inf;
info(4) = info(2); info(4) = info(2);
exit_flag = 0; exit_flag = 0;
if options_mom_.vector_output == 1 % lsqnonlin requires vector output if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number; fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
end end
return return
@ -118,31 +119,32 @@ if info(1)
fval = Inf; fval = Inf;
info(4) = 0.1; info(4) = 0.1;
exit_flag = 0; exit_flag = 0;
if options_mom_.vector_output == 1 % lsqnonlin requires vector output if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number; fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
end end
return return
end end
end end
%--------------------------------------------------------------------------
% GMM: Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') if strcmp(options_mom_.mom.mom_method,'GMM')
%--------------------------------------------------------------------------
% 3. Set up pruned state-space system and compute model moments
%--------------------------------------------------------------------------
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
indpmodel = []; %initialize index for model parameters indpmodel = []; % initialize index for model parameters
if ~isempty(estim_params_.param_vals) if ~isempty(estim_params_.param_vals)
indpmodel = estim_params_.param_vals(:,1); %values correspond to parameters declaration order, row number corresponds to order in estimated_params indpmodel = estim_params_.param_vals(:,1); % values correspond to parameters declaration order, row number corresponds to order in estimated_params
end end
indpstderr=[]; %initialize index for stderr parameters indpstderr=[]; % initialize index for stderr parameters
if ~isempty(estim_params_.var_exo) if ~isempty(estim_params_.var_exo)
indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params indpstderr = estim_params_.var_exo(:,1); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
end end
indpcorr=[]; %initialize matrix for corr paramters indpcorr=[]; % initialize matrix for corr paramters
if ~isempty(estim_params_.corrx) if ~isempty(estim_params_.corrx)
indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params indpcorr = estim_params_.corrx(:,1:2); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
end end
if estim_params_.nvn || estim_params_.ncn %nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params if estim_params_.nvn || estim_params_.ncn % nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params
error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.') error('Analytic computation of standard errrors does not (yet) support measurement errors.\nInstead, define them explicitly as varexo and provide measurement equations in the model definition.\nAlternatively, use numerical standard errors.')
end end
modparam_nbr = estim_params_.np; % number of model parameters as declared in estimated_params modparam_nbr = estim_params_.np; % number of model parameters as declared in estimated_params
@ -151,27 +153,26 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr; totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr;
dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices dr.derivs = get_perturbation_params_derivs(M_, options_mom_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, 0); %analytic derivatives of perturbation matrices
oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr); oo_.mom.model_moments_params_derivs = NaN(options_mom_.mom.mom_nbr,totparam_nbr);
pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 1); pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.mom.obs_var, options_mom_.ar, 0, 1);
else else
pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0); pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.mom.obs_var, options_mom_.ar, 0, 0);
end end
oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1); oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1);
for jm = 1:size(M_.matched_moments,1) for jm = 1:size(M_.matched_moments,1)
% First moments % First moments
if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1) if ~options_mom_.prefilter && (sum(M_.matched_moments{jm,3}) == 1)
idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) ); idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}) );
oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1); oo_.mom.model_moments(jm,1) = pruned_state_space.E_y(idx1);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:); oo_.mom.model_moments_params_derivs(jm,:) = pruned_state_space.dE_y(idx1,:);
end end
end end
% Second moments % second moments
if (sum(M_.matched_moments{jm,3}) == 2) if (sum(M_.matched_moments{jm,3}) == 2)
idx1 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) ); idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) );
idx2 = (oo_.dr.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) ); idx2 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) );
if nnz(M_.matched_moments{jm,2}) == 0 if nnz(M_.matched_moments{jm,2}) == 0
% Covariance % covariance
if options_mom_.prefilter if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2); oo_.mom.model_moments(jm,1) = pruned_state_space.Var_y(idx1,idx2);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
@ -186,8 +187,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end end
end end
else else
% Autocovariance % autocovariance
lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in matched_moments are transformed such that first entry is always 0 and the second is a lag lag = -M_.matched_moments{jm,2}(2); %note that leads/lags in M_.matched_moments are transformed such that first entry is always 0 and the second is a lag
if options_mom_.prefilter if options_mom_.prefilter
oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag); oo_.mom.model_moments(jm,1) = pruned_state_space.Var_yi(idx1,idx2,lag);
if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian ) if options_mom_.mom.compute_derivs && ( options_mom_.mom.analytic_standard_errors || options_mom_.mom.analytic_jacobian )
@ -204,24 +205,23 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
end end
end end
end end
end
elseif strcmp(options_mom_.mom.mom_method,'SMM') %------------------------------------------------------------------------------
%------------------------------------------------------------------------------ % SMM: Compute Moments of the model solution for Gaussian innovations
% 3. Compute Moments of the model solution for normal innovations %------------------------------------------------------------------------------
%------------------------------------------------------------------------------ if strcmp(options_mom_.mom.mom_method,'SMM')
% create shock series with correct covariance matrix from iid standard normal shocks % create shock series with correct covariance matrix from iid standard normal shocks
i_exo_var = setdiff(1:M_.exo_nbr, find(diag(M_.Sigma_e) == 0 )); %find singular entries in covariance i_exo_var = setdiff(1:M_.exo_nbr, find(diag(M_.Sigma_e) == 0 )); % find singular entries in covariance
chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var)); chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
scaled_shock_series = zeros(size(options_mom_.mom.shock_series)); %initialize scaled_shock_series = zeros(size(options_mom_.mom.shock_series)); % initialize
scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; %set non-zero entries scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; % set non-zero entries
% simulate series % simulate series
y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order); y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf % provide meaningful penalty if data is nan or inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim))) if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.vector_output == 1 % lsqnonlin requires vector output if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = Inf(size(oo_.mom.Sw,1),1); fval = Inf(size(oo_.mom.Sw,1),1);
else else
fval = Inf; fval = Inf;
@ -229,73 +229,70 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
info(1)=180; info(1)=180;
info(4) = 0.1; info(4) = 0.1;
exit_flag = 0; exit_flag = 0;
if options_mom_.vector_output == 1 % lsqnonlin requires vector output if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number; fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number;
end end
return return
end end
% remove burn-in and focus on observables (note that y_sim is in declaration order)
% Remove burn-in and focus on observables (note that y_sim is in declaration order) y_sim = y_sim(oo_.dr.order_var(oo_.mom.obs_var) , end-options_mom_.mom.long+1:end)';
y_sim = y_sim(oo_.dr.order_var(oo_.dr.obs_var) , end-options_mom_.mom.long+1:end)';
if ~all(diag(M_.H)==0) if ~all(diag(M_.H)==0)
i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance
chol_S = chol(M_.H(i_ME,i_ME)); %decompose rest chol_S = chol(M_.H(i_ME,i_ME)); % decompose rest
shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); %initialize shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); % initialize
shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_ME)*chol_S; shock_mat(:,i_ME)=options_mom_.mom.ME_shock_series(:,i_ME)*chol_S;
y_sim = y_sim+shock_mat; y_sim = y_sim+shock_mat;
end end
% remove mean if centered moments
% Remove mean if centered moments
if options_mom_.prefilter if options_mom_.prefilter
y_sim = bsxfun(@minus, y_sim, mean(y_sim,1)); y_sim = bsxfun(@minus, y_sim, mean(y_sim,1));
end end
oo_.mom.model_moments = mom.data_moments(y_sim, oo_, M_.matched_moments, options_mom_); oo_.mom.model_moments = mom.get_data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
end end
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
% 4. Compute quadratic target function % Compute quadratic target function
%-------------------------------------------------------------------------- %--------------------------------------------------------------------------
moments_difference = oo_.mom.data_moments - oo_.mom.model_moments; moments_difference = oo_.mom.data_moments - oo_.mom.model_moments;
residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference;
oo_.mom.Q = residuals'*residuals;
if options_mom_.vector_output == 1 % lsqnonlin requires vector output
fval = residuals;
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam1-oo_.prior.mean)./sqrt(diag(oo_.prior.variance))];
end
else
fval = oo_.mom.Q;
if options_mom_.mom.penalized_estimator
fval=fval+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(xparam1-oo_.prior.mean);
end
end
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
if options_mom_.mom.penalized_estimator residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference;
dxparam1 = eye(length(xparam1)); oo_.mom.Q = residuals'*residuals;
if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
fval = residuals;
if options_mom_.mom.penalized_estimator
fval=[fval;(xparam-oo_.mom.prior.mean)./sqrt(diag(oo_.mom.prior.variance))];
end
else
fval = oo_.mom.Q;
if options_mom_.mom.penalized_estimator
fval=fval+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean);
end
end end
if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
for jp=1:length(xparam1) if options_mom_.mom.penalized_estimator
dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp); dxparam1 = eye(length(xparam));
dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference; end
for jp=1:length(xparam)
if options_mom_.vector_output == 1 % lsqnonlin requires vector output dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp);
if options_mom_.mom.penalized_estimator dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference;
df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.prior.variance))]; if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
if options_mom_.mom.penalized_estimator
df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.mom.prior.variance))];
else
df(:,jp) = dresiduals;
end
else else
df(:,jp) = dresiduals; df(:,jp) = dresiduals'*residuals + residuals'*dresiduals;
end if options_mom_.mom.penalized_estimator
else df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.mom.prior.variance*(xparam-oo_.mom.prior.mean)+(xparam-oo_.mom.prior.mean)'/oo_.mom.prior.variance*(dxparam1(:,jp));
df(:,jp) = dresiduals'*residuals + residuals'*dresiduals; end
if options_mom_.mom.penalized_estimator
df(:,jp)=df(:,jp)+(dxparam1(:,jp))'/oo_.prior.variance*(xparam1-oo_.prior.mean)+(xparam1-oo_.prior.mean)'/oo_.prior.variance*(dxparam1(:,jp));
end end
end end
end end
end end
end%main function end end % main function end

123
matlab/+mom/print_info_on_estimation_settings.m Normal file
View File

@ -0,0 +1,123 @@
function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters)
% function print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters)
% -------------------------------------------------------------------------
% Print information on the method of moments estimation settings to the console
% =========================================================================
% INPUTS
% options_mom_ [struct] Options for the method of moments estimation
% number_of_estimated_parameters [integer] Number of estimated parameters
% -------------------------------------------------------------------------
% OUTPUT
% No output, just displays the chosen settings
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% -------------------------------------------------------------------------
% This function calls
% o skipline
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
fprintf('\n---------------------------------------------------\n')
if strcmp(options_mom_.mom.mom_method,'SMM')
fprintf('Simulated method of moments with');
elseif strcmp(options_mom_.mom.mom_method,'GMM')
fprintf('General method of moments with');
end
if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
if options_mom_.prefilter
fprintf('\n - centered moments (prefilter=1)');
else
fprintf('\n - uncentered moments (prefilter=0)');
end
if options_mom_.mom.penalized_estimator
fprintf('\n - penalized estimation using deviation from prior mean and weighted with prior precision');
end
end
for i = 1:length(options_mom_.optimizer_vec)
if i == 1
str = '- optimizer (mode_compute';
else
str = ' (additional_optimizer_steps';
end
switch options_mom_.optimizer_vec{i}
case 0
fprintf('\n %s=0): no minimization',str);
case 1
fprintf('\n %s=1): fmincon',str);
case 2
fprintf('\n %s=2): continuous simulated annealing',str);
case 3
fprintf('\n %s=3): fminunc',str);
case 4
fprintf('\n %s=4): csminwel',str);
case 5
fprintf('\n %s=5): newrat',str);
case 6
fprintf('\n %s=6): gmhmaxlik',str);
case 7
fprintf('\n %s=7): fminsearch',str);
case 8
fprintf('\n %s=8): Dynare Nelder-Mead simplex',str);
case 9
fprintf('\n %s=9): CMA-ES',str);
case 10
fprintf('\n %s=10): simpsa',str);
case 11
skipline;
error('method_of_moments: online_auxiliary_filter (mode_compute=11) is only supported with likelihood-based estimation techniques!');
case 12
fprintf('\n %s=12): particleswarm',str);
case 101
fprintf('\n %s=101): SolveOpt',str);
case 102
fprintf('\n %s=102): simulannealbnd',str);
case 13
fprintf('\n %s=13): lsqnonlin',str);
otherwise
if ischar(options_mom_.optimizer_vec{i})
fprintf('\n %s=%s): user-defined',str,options_mom_.optimizer_vec{i});
else
skipline;
error('method_of_moments: Unknown optimizer!');
end
end
if options_mom_.silent_optimizer
fprintf(' (silent)');
end
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(options_mom_.optimizer_vec{i},options_mom_.mom.analytic_jacobian_optimizers)
fprintf(' (using analytical Jacobian)');
end
end
if options_mom_.order > 0
fprintf('\n - stochastic simulations with perturbation order: %d', options_mom_.order)
end
if options_mom_.order > 1 && options_mom_.pruning
fprintf(' (with pruning)')
end
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors
fprintf('\n - standard errors: analytic derivatives');
else
fprintf('\n - standard errors: numerical derivatives');
end
fprintf('\n - number of matched moments: %d', options_mom_.mom.mom_nbr);
end
fprintf('\n - number of parameters: %d', number_of_estimated_parameters);
fprintf('\n\n');

1158
matlab/+mom/run.m
View File

File diff suppressed because it is too large Load Diff

43
matlab/+mom/set_correct_bounds_for_stderr_corr.m Normal file
View File

@ -0,0 +1,43 @@
function Bounds = set_correct_bounds_for_stderr_corr(estim_params_,Bounds)
% function Bounds = set_correct_bounds_for_stderr_corr(estim_params_,Bounds)
% -------------------------------------------------------------------------
% Set correct bounds for standard deviation and corrrelation parameters
% =========================================================================
% INPUTS
% o estim_params_ [struct] information on estimated parameters
% o Bounds [struct] information on bounds
% -------------------------------------------------------------------------
% OUTPUT
% o Bounds [struct] updated bounds
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
number_of_estimated_parameters = estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np;
% set correct bounds for standard deviations and corrrelations
param_of_interest = (1:number_of_estimated_parameters)'<=estim_params_.nvx+estim_params_.nvn;
LB_below_0 = (Bounds.lb<0 & param_of_interest);
Bounds.lb(LB_below_0) = 0;
param_of_interest = (1:number_of_estimated_parameters)'> estim_params_.nvx+estim_params_.nvn & (1:number_of_estimated_parameters)'<estim_params_.nvx+estim_params_.nvn +estim_params_.ncx + estim_params_.ncn;
LB_below_minus_1 = (Bounds.lb<-1 & param_of_interest);
UB_above_1 = (Bounds.ub>1 & param_of_interest);
Bounds.lb(LB_below_minus_1) = -1;
Bounds.ub(UB_above_1) = 1;

58
matlab/+mom/transform_prior_to_laplace_prior.m Normal file
View File

@ -0,0 +1,58 @@
function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% function bayestopt_ = transform_prior_to_laplace_prior(bayestopt_)
% -------------------------------------------------------------------------
% Transforms the prior specification to a Laplace type of approximation:
% only the prior mean and standard deviation are relevant, all other shape
% information, except for the parameter bounds, is ignored.
% =========================================================================
% INPUTS
% bayestopt_ [structure] prior information
% -------------------------------------------------------------------------
% OUTPUT
% bayestopt_ [structure] Laplace prior information
% -------------------------------------------------------------------------
% This function is called by
% o mom.run
% =========================================================================
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
if any(setdiff([0;bayestopt_.pshape],[0,3]))
fprintf('\nNon-normal priors specified. Penalized estimation uses a Laplace type of approximation:');
fprintf('\nOnly the prior mean and standard deviation are relevant, all other shape information, except for the parameter bounds, is ignored.\n\n');
non_normal_priors = (bayestopt_.pshape~=3);
bayestopt_.pshape(non_normal_priors) = 3;
bayestopt_.p3(non_normal_priors) = -Inf*ones(sum(non_normal_priors),1);
bayestopt_.p4(non_normal_priors) = Inf*ones(sum(non_normal_priors),1);
bayestopt_.p6(non_normal_priors) = bayestopt_.p1(non_normal_priors);
bayestopt_.p7(non_normal_priors) = bayestopt_.p2(non_normal_priors);
bayestopt_.p5(non_normal_priors) = bayestopt_.p1(non_normal_priors);
end
if any(isinf(bayestopt_.p2)) % find infinite variance priors
inf_var_pars = bayestopt_.name(isinf(bayestopt_.p2));
disp_string = [inf_var_pars{1,:}];
for ii = 2:size(inf_var_pars,1)
disp_string = [disp_string,', ',inf_var_pars{ii,:}];
end
fprintf('The parameter(s) %s have infinite prior variance. This implies a flat prior.\n',disp_string);
fprintf('Dynare disables the matrix singularity warning\n');
if isoctave
warning('off','Octave:singular-matrix');
else
warning('off','MATLAB:singularMatrix');
end
end

4
matlab/CheckPath.m
View File

@ -12,7 +12,7 @@ function [DirectoryName, info] = CheckPath(type,dname)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2005-2017 Dynare Team % Copyright © 2005-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -31,7 +31,7 @@ function [DirectoryName, info] = CheckPath(type,dname)
info = 0; info = 0;
DirectoryName = [ dname '/' type ]; DirectoryName = [ dname filesep type ];
if ~isdir(dname) if ~isdir(dname)
% Make sure there isn't a file with the same name, see trac ticket #47 % Make sure there isn't a file with the same name, see trac ticket #47

14
matlab/CutSample.m
View File

@ -33,6 +33,8 @@ function CutSample(M_, options_, estim_params_)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
dispString = 'Estimation::mcmc';
% Get the path to the metropolis files. % Get the path to the metropolis files.
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
@ -47,7 +49,7 @@ npar=size(record.LastParameters,2);
mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]); mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]);
if ~length(mh_files) if ~length(mh_files)
error('Estimation::mcmc: I can''t find MH file to load here!') error('%s: I can''t find MH file to load here!',dispString)
end end
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
@ -71,9 +73,9 @@ end
% Save updated mh-history file. % Save updated mh-history file.
update_last_mh_history_file(MetropolisFolder, ModelName, record); update_last_mh_history_file(MetropolisFolder, ModelName, record);
fprintf('Estimation::mcmc: Total number of MH draws per chain: %d.\n',TotalNumberOfMhDraws); fprintf('%s: Total number of MH draws per chain: %d.\n',dispString,TotalNumberOfMhDraws);
fprintf('Estimation::mcmc: Total number of generated MH files: %d.\n',TotalNumberOfMhFiles); fprintf('%s: Total number of generated MH files: %d.\n',dispString,TotalNumberOfMhFiles);
fprintf('Estimation::mcmc: I''ll use mh-files %d to %d.\n',FirstMhFile,TotalNumberOfMhFiles); fprintf('%s: I''ll use mh-files %d to %d.\n',dispString,FirstMhFile,TotalNumberOfMhFiles);
fprintf('Estimation::mcmc: In MH-file number %d I''ll start at line %d.\n',FirstMhFile,FirstLine); fprintf('%s: In MH-file number %d I''ll start at line %d.\n',dispString,FirstMhFile,FirstLine);
fprintf('Estimation::mcmc: Finally I keep %d draws per chain.\n',TotalNumberOfMhDraws-FirstDraw+1); fprintf('%s: Finally I keep %d draws per chain.\n',dispString,TotalNumberOfMhDraws-FirstDraw+1);
skipline() skipline()

10
matlab/PosteriorIRF.m
View File

@ -16,7 +16,7 @@ function PosteriorIRF(type)
% functions associated with it(the _core1 and _core2). % functions associated with it(the _core1 and _core2).
% See also the comments posterior_sampler.m funtion. % See also the comments posterior_sampler.m funtion.
% Copyright © 2006-2018 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -36,6 +36,8 @@ function PosteriorIRF(type)
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
dispString = 'Estimation::mcmc';
% Set the number of periods % Set the number of periods
if isempty(options_.irf) || ~options_.irf if isempty(options_.irf) || ~options_.irf
options_.irf = 40; options_.irf = 40;
@ -287,7 +289,7 @@ if options_.TeX
end end
end end
fprintf('Estimation::mcmc: Posterior (dsge) IRFs...\n'); fprintf('%s: Posterior (dsge) IRFs...\n',dispString);
tit = M_.exo_names; tit = M_.exo_names;
kdx = 0; kdx = 0;
@ -327,7 +329,7 @@ if MAX_nirfs_dsgevar
VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr); VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr); DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr); HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
fprintf('Estimation::mcmc: Posterior (bvar-dsge) IRFs...\n'); fprintf('%s: Posterior (bvar-dsge) IRFs...\n',dispString);
tit = M_.exo_names; tit = M_.exo_names;
kdx = 0; kdx = 0;
for file = 1:NumberOfIRFfiles_dsgevar for file = 1:NumberOfIRFfiles_dsgevar
@ -457,4 +459,4 @@ if ~options_.nograph && ~options_.no_graph.posterior
end end
fprintf('Estimation::mcmc: Posterior IRFs, done!\n'); fprintf('%s: Posterior IRFs, done!\n',dispString);

13
matlab/check_posterior_sampler_options.m
View File

@ -1,10 +1,11 @@
function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_) function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_)
% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, fname, dname, options_, bounds, bayestopt_)
% function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_)
% initialization of posterior samplers % initialization of posterior samplers
% %
% INPUTS % INPUTS
% posterior_sampler_options: posterior sampler options % posterior_sampler_options: posterior sampler options
% fname: name of the mod-file
% dname: name of directory with metropolis folder
% options_: structure storing the options % options_: structure storing the options
% bounds: structure containing prior bounds % bounds: structure containing prior bounds
% bayestopt_: structure storing information about priors % bayestopt_: structure storing information about priors
@ -17,7 +18,7 @@ function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sam
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2015-2022 Dynare Team % Copyright © 2015-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -391,7 +392,7 @@ if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
end end
if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix if options_.load_mh_file && posterior_sampler_options.use_mh_covariance_matrix
[~, invhess] = compute_mh_covariance_matrix; [~, invhess] = compute_mh_covariance_matrix(bayestopt_,fname,dname);
posterior_sampler_options.invhess = invhess; posterior_sampler_options.invhess = invhess;
end end
@ -413,7 +414,7 @@ if strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
error('check_posterior_sampler_options:: This error should not occur, please contact developers.') error('check_posterior_sampler_options:: This error should not occur, please contact developers.')
end end
% % % if options_.load_mh_file && options_.use_mh_covariance_matrix, % % % if options_.load_mh_file && options_.use_mh_covariance_matrix,
% % % [~, invhess] = compute_mh_covariance_matrix; % % % [~, invhess] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname));
% % % posterior_sampler_options.invhess = invhess; % % % posterior_sampler_options.invhess = invhess;
% % % end % % % end
[V1, D]=eig(invhess); [V1, D]=eig(invhess);

21
matlab/compute_mh_covariance_matrix.m
View File

@ -1,10 +1,13 @@
function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix() function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname)
% function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,fname,dname)
% Estimation of the posterior covariance matrix, posterior mean, posterior mode and evaluation of the posterior kernel at the % 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 % estimated mode, using the draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in
% a file <M_.fname>_mh_mode.mat. % a file <fname>_mh_mode.mat.
% %
% INPUTS % INPUTS
% None. % o bayestopt_ [struct] characterizing priors
% o fname [string] name of model
% o dname [string] name of directory with metropolis folder
% %
% OUTPUTS % OUTPUTS
% o posterior_mean [double] (n*1) vector, posterior expectation of the parameters. % o posterior_mean [double] (n*1) vector, posterior expectation of the parameters.
@ -31,14 +34,10 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
MetropolisFolder = CheckPath('metropolis',dname);
BaseName = [MetropolisFolder filesep fname];
global M_ bayestopt_ record=load_last_mh_history_file(MetropolisFolder, fname);
MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
record=load_last_mh_history_file(MetropolisFolder, ModelName);
FirstMhFile = record.KeepedDraws.FirstMhFile; FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine; FirstLine = record.KeepedDraws.FirstLine;
@ -71,4 +70,4 @@ hh = inv(posterior_covariance);
fval = posterior_kernel_at_the_mode; fval = posterior_kernel_at_the_mode;
parameter_names = bayestopt_.name; parameter_names = bayestopt_.name;
save([M_.dname filesep 'Output' filesep M_.fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names'); save([dname filesep 'Output' filesep fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names');

24
matlab/convergence_diagnostics/McMCDiagnostics.m → matlab/convergence_diagnostics/mcmc_diagnostics.m
View File

@ -1,5 +1,5 @@
function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) function oo_ = mcmc_diagnostics(options_, estim_params_, M_, oo_)
% function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) % function oo_ = mcmc_diagnostics(options_, estim_params_, M_, oo_)
% Computes convergence tests % Computes convergence tests
% %
% INPUTS % INPUTS
@ -51,11 +51,11 @@ for b = 1:nblck
nfiles = length(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck' num2str(b) '.mat'])); nfiles = length(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck' num2str(b) '.mat']));
if ~isequal(NumberOfMcFilesPerBlock,nfiles) if ~isequal(NumberOfMcFilesPerBlock,nfiles)
issue_an_error_message = 1; issue_an_error_message = 1;
disp(['Estimation::mcmc::diagnostics: The number of MCMC files in chain ' num2str(b) ' is ' num2str(nfiles) ' while the mh-history files indicate that we should have ' num2str(NumberOfMcFilesPerBlock) ' MCMC files per chain!']) fprintf('The number of MCMC files in chain %u is %u while the mh-history files indicate that we should have %u MCMC files per chain!\n',b, nfiles, NumberOfMcFilesPerBlock);
end end
end end
if issue_an_error_message if issue_an_error_message
error('Estimation::mcmc::diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...') error('mcmc_diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...!');
end end
% compute inefficiency factor % compute inefficiency factor
@ -111,7 +111,7 @@ PastDraws = sum(record.MhDraws,1);
NumberOfDraws = PastDraws(1); NumberOfDraws = PastDraws(1);
if NumberOfDraws<=2000 if NumberOfDraws<=2000
warning(['estimation:: MCMC convergence diagnostics are not computed because the total number of iterations is not bigger than 2000!']) warning('MCMC convergence diagnostics are not computed because the total number of iterations is not bigger than 2000!');
return return
end end
@ -185,7 +185,7 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
if options_.convergence.rafterylewis.indicator if options_.convergence.rafterylewis.indicator
if any(options_.convergence.rafterylewis.qrs<0) || any(options_.convergence.rafterylewis.qrs>1) || length(options_.convergence.rafterylewis.qrs)~=3 ... if any(options_.convergence.rafterylewis.qrs<0) || any(options_.convergence.rafterylewis.qrs>1) || length(options_.convergence.rafterylewis.qrs)~=3 ...
|| (options_.convergence.rafterylewis.qrs(1)-options_.convergence.rafterylewis.qrs(2)<=0) || (options_.convergence.rafterylewis.qrs(1)-options_.convergence.rafterylewis.qrs(2)<=0)
fprintf('\nCONVERGENCE DIAGNOSTICS: Invalid option for raftery_lewis_qrs. Using the default of [0.025 0.005 0.95].\n') fprintf('\nInvalid option for raftery_lewis_qrs. Using the default of [0.025 0.005 0.95].\n');
options_.convergence.rafterylewis.qrs=[0.025 0.005 0.95]; options_.convergence.rafterylewis.qrs=[0.025 0.005 0.95];
end end
Raftery_Lewis_q=options_.convergence.rafterylewis.qrs(1); Raftery_Lewis_q=options_.convergence.rafterylewis.qrs(1);
@ -218,18 +218,18 @@ xx = Origin:StepSize:NumberOfDraws;
NumberOfLines = length(xx); NumberOfLines = length(xx);
if NumberOfDraws < Origin if NumberOfDraws < Origin
disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.') fprintf('The number of simulations is too small to compute the MCMC convergence diagnostics.\n');
return return
end end
if TeX && any(strcmp('eps',cellstr(options_.graph_format))) if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([OutputFolder '/' ModelName '_UnivariateDiagnostics.tex'],'w'); fidTeX = fopen([OutputFolder '/' ModelName '_UnivariateDiagnostics.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by mcmc_diagnostics.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
end end
disp('Estimation::mcmc::diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):') fprintf('Univariate convergence diagnostic, Brooks and Gelman (1998):\n');
% The mandatory variables for local/remote parallel % The mandatory variables for local/remote parallel
% computing are stored in localVars struct. % computing are stored in localVars struct.
@ -248,7 +248,7 @@ localVars.M_ = M_;
% Like sequential execution! % Like sequential execution!
if isnumeric(options_.parallel) if isnumeric(options_.parallel)
fout = McMCDiagnostics_core(localVars,1,npar,0); fout = mcmc_diagnostics_core(localVars,1,npar,0);
UDIAG = fout.UDIAG; UDIAG = fout.UDIAG;
clear fout clear fout
% Parallel execution! % Parallel execution!
@ -258,7 +258,7 @@ else
end end
NamFileInput={[M_.dname '/metropolis/'],[ModelName '_mh*_blck*.mat']}; NamFileInput={[M_.dname '/metropolis/'],[ModelName '_mh*_blck*.mat']};
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,NamFileInput,'McMCDiagnostics_core', localVars, [], options_.parallel_info); [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,NamFileInput,'mcmc_diagnostics_core', localVars, [], options_.parallel_info);
UDIAG = fout(1).UDIAG; UDIAG = fout(1).UDIAG;
for j=2:totCPU for j=2:totCPU
UDIAG = cat(3,UDIAG ,fout(j).UDIAG); UDIAG = cat(3,UDIAG ,fout(j).UDIAG);
@ -397,7 +397,7 @@ clear UDIAG;
% %
if TeX && any(strcmp('eps',cellstr(options_.graph_format))) if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([OutputFolder '/' ModelName '_MultivariateDiagnostics.tex'],'w'); fidTeX = fopen([OutputFolder '/' ModelName '_MultivariateDiagnostics.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); fprintf(fidTeX,'%% TeX eps-loader file generated by mcmc_diagnostics.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
end end

13
matlab/convergence_diagnostics/McMCDiagnostics_core.m → matlab/convergence_diagnostics/mcmc_diagnostics_core.m
View File

@ -1,5 +1,5 @@
function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab) function myoutput = mcmc_diagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab) % function myoutput = mcmc_diagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% Computes the Brooks/Gelman (1998) convergence diagnostics, both the % Computes the Brooks/Gelman (1998) convergence diagnostics, both the
% parameteric and the non-parameteric versions % parameteric and the non-parameteric versions
% %
@ -20,11 +20,10 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% 4nd column: sum of within sequence variances; used to compute mean within sequence variances % 4nd column: sum of within sequence variances; used to compute mean within sequence variances
% 5nd column: within sequence kurtosis % 5nd column: within sequence kurtosis
% 6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses % 6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses
% Averaging to compute mean moments is done in % Averaging to compute mean moments is done in mcmc_diagnostics
% McMCDiagnostics
% %
% ALGORITHM % ALGORITHM
% Computes part of the convergence diagnostics, the rest is computed in McMCDiagnostics.m . % Computes part of the convergence diagnostics, the rest is computed in mcmc_diagnostics.m.
% The methodology and terminology is based on: Brooks/Gelman (1998): General % The methodology and terminology is based on: Brooks/Gelman (1998): General
% Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational % Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational
% and Graphical Statistics, Volume 7, Number 4, Pages 434-455 % and Graphical Statistics, Volume 7, Number 4, Pages 434-455
@ -33,7 +32,7 @@ function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
% SPECIAL REQUIREMENTS. % SPECIAL REQUIREMENTS.
% None. % None.
% Copyright © 2006-2017 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -80,7 +79,7 @@ tmp = zeros(NumberOfDraws*nblck,3);
UDIAG = zeros(NumberOfLines,6,npar-fpar+1); UDIAG = zeros(NumberOfLines,6,npar-fpar+1);
if whoiam if whoiam
waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...']; waitbarString = ['Please wait... MCMC diagnostics (' int2str(fpar) 'of' int2str(npar) ')...'];
if Parallel(ThisMatlab).Local if Parallel(ThisMatlab).Local
waitbarTitle=['Local ']; waitbarTitle=['Local '];
else else

2
matlab/display_estimation_results_table.m
View File

@ -174,7 +174,7 @@ if ncn
end end
if any(xparam1(1:nvx+nvn)<0) if any(xparam1(1:nvx+nvn)<0)
warning('Some estimated standard deviations are negative. Dynare internally works with variances so that the sign does not matter. Nevertheless, it is recommended to impose either prior restrictions (Bayesian Estimation) or a lower bound (ML) to assure positive values.') warning(sprintf('Some estimated standard deviations are negative.\n Dynare internally works with variances so that the sign does not matter.\n Nevertheless, it is recommended to impose either prior restrictions (Bayesian Estimation)\n or a lower bound (ML) to assure positive values.'))
end end
OutputDirectoryName = CheckPath('Output',M_.dname); OutputDirectoryName = CheckPath('Output',M_.dname);

118
matlab/dynare_estimation_1.m
View File

@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright © 2003-2022 Dynare Team % Copyright © 2003-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -31,6 +31,8 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
dispString = 'Estimation::mcmc';
if ~exist([M_.dname filesep 'Output'],'dir') if ~exist([M_.dname filesep 'Output'],'dir')
if isoctave && octave_ver_less_than('7') && ~exist(M_.dname) if isoctave && octave_ver_less_than('7') && ~exist(M_.dname)
% See https://savannah.gnu.org/bugs/index.php?61166 % See https://savannah.gnu.org/bugs/index.php?61166
@ -293,37 +295,7 @@ 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
try check_hessian_at_the_mode(hh, xparam1, M_, estim_params_, options_, bounds);
chol(hh);
catch
skipline()
disp('POSTERIOR KERNEL OPTIMIZATION PROBLEM!')
disp(' (minus) the hessian matrix at the "mode" is not positive definite!')
disp('=> posterior variance of the estimated parameters are not positive.')
disp('You should try to change the initial values of the parameters using')
disp('the estimated_params_init block, or use another optimization routine.')
params_at_bound=find(abs(xparam1-bounds.ub)<1.e-10 | abs(xparam1-bounds.lb)<1.e-10);
if ~isempty(params_at_bound)
for ii=1:length(params_at_bound)
params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_);
end
disp_string=[params_at_bound_name{1,:}];
for ii=2:size(params_at_bound_name,1)
disp_string=[disp_string,', ',params_at_bound_name{ii,:}];
end
fprintf('\nThe following parameters are at the prior bound: %s\n', disp_string)
fprintf('Some potential solutions are:\n')
fprintf(' - Check your model for mistakes.\n')
fprintf(' - Check whether model and data are consistent (correct observation equation).\n')
fprintf(' - Shut off prior_trunc.\n')
fprintf(' - Change the optimization bounds.\n')
fprintf(' - Use a different mode_compute like 6 or 9.\n')
fprintf(' - Check whether the parameters estimated are identified.\n')
fprintf(' - Check prior shape (e.g. Inf density at bound(s)).\n')
fprintf(' - Increase the informativeness of the prior.\n')
end
warning('The results below are most likely wrong!');
end
end end
if options_.mode_check.status && ~options_.mh_posterior_mode_estimation if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
@ -404,73 +376,19 @@ if np > 0
save([M_.dname filesep 'Output' filesep M_.fname '_params.mat'],'pindx'); save([M_.dname filesep 'Output' filesep M_.fname '_params.mat'],'pindx');
end end
switch options_.MCMC_jumping_covariance invhess = set_mcmc_jumping_covariance(invhess, nx, options_.MCMC_jumping_covariance, bayestopt_, 'dynare_estimation_1');
case 'hessian' %Baseline
%do nothing and use hessian from mode_compute
case 'prior_variance' %Use prior variance
if any(isinf(bayestopt_.p2))
error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.')
else
hh = diag(1./(bayestopt_.p2.^2));
end
hsd = sqrt(diag(hh));
invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
case 'identity_matrix' %Use identity
invhess = eye(nx);
otherwise %user specified matrix in file
try
load(options_.MCMC_jumping_covariance,'jumping_covariance')
hh=jumping_covariance;
catch
error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat'])
end
[nrow, ncol]=size(hh);
if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size
error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns'])
end
try %check for positive definiteness
chol(hh);
hsd = sqrt(diag(hh));
invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
catch
error(['Specified jumping_covariance is not positive definite'])
end
end
if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
(any(bayestopt_.pshape >0 ) && options_.load_mh_file) %% not ML estimation (any(bayestopt_.pshape >0 ) && options_.load_mh_file) %% not ML estimation
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding %reset bounds as lb and ub must only be operational during mode-finding
outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub); bounds = set_mcmc_prior_bounds(xparam1, bayestopt_, options_, 'dynare_estimation_1');
if ~isempty(outside_bound_pars)
for ii=1:length(outside_bound_pars)
outside_bound_par_names{ii,1}=get_the_name(ii,0,M_,estim_params_,options_);
end
disp_string=[outside_bound_par_names{1,:}];
for ii=2:size(outside_bound_par_names,1)
disp_string=[disp_string,', ',outside_bound_par_names{ii,:}];
end
if options_.prior_trunc>0
error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0.'])
else
error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds.'])
end
end
% Tunes the jumping distribution's scale parameter % Tunes the jumping distribution's scale parameter
if options_.mh_tune_jscale.status if options_.mh_tune_jscale.status
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings') if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
%get invhess in case of use_mh_covariance_matrix options_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds,...
posterior_sampler_options_temp = options_.posterior_sampler_options.current_options; dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_);
posterior_sampler_options_temp.invhess = invhess;
posterior_sampler_options_temp = check_posterior_sampler_options(posterior_sampler_options_temp, options_);
options = options_.mh_tune_jscale;
options.rwmh = options_.posterior_sampler_options.rwmh;
options_.mh_jscale = calibrate_mh_scale_parameter(objective_function, ...
posterior_sampler_options_temp.invhess, xparam1, [bounds.lb,bounds.ub], ...
options, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_);
clear('posterior_sampler_options_temp','options')
bayestopt_.jscale(:) = options_.mh_jscale; bayestopt_.jscale(:) = options_.mh_jscale;
fprintf('mh_jscale has been set equal to %s\n', num2str(options_.mh_jscale)) fprintf('mh_jscale has been set equal to %s\n', num2str(options_.mh_jscale));
else else
warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!') warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!')
end end
@ -479,7 +397,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if options_.mh_replic || options_.load_mh_file if options_.mh_replic || options_.load_mh_file
posterior_sampler_options = options_.posterior_sampler_options.current_options; posterior_sampler_options = options_.posterior_sampler_options.current_options;
posterior_sampler_options.invhess = invhess; posterior_sampler_options.invhess = invhess;
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, options_, bounds, bayestopt_); [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, M_.fname, M_.dname, options_, bounds, bayestopt_);
% store current options in global % store current options in global
options_.posterior_sampler_options.current_options = posterior_sampler_options; options_.posterior_sampler_options.current_options = posterior_sampler_options;
if options_.mh_replic if options_.mh_replic
@ -492,7 +410,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
%% Here I discard first mh_drop percent of the draws: %% Here I discard first mh_drop percent of the draws:
CutSample(M_, options_, estim_params_); CutSample(M_, options_, estim_params_);
if options_.mh_posterior_mode_estimation if options_.mh_posterior_mode_estimation
[~,~,posterior_mode,~] = compute_mh_covariance_matrix(); [~,~,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'); oo_=fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_,estim_params_,bayestopt_,oo_,'posterior');
%reset qz_criterium %reset qz_criterium
options_.qz_criterium=qz_criterium_old; options_.qz_criterium=qz_criterium_old;
@ -504,7 +422,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end end
if ~options_.nodiagnostic if ~options_.nodiagnostic
if (options_.mh_replic>0 || (options_.load_mh_file && ~options_.load_results_after_load_mh)) if (options_.mh_replic>0 || (options_.load_mh_file && ~options_.load_results_after_load_mh))
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_); oo_= mcmc_diagnostics(options_, estim_params_, M_,oo_);
elseif options_.load_mh_file && options_.load_results_after_load_mh elseif options_.load_mh_file && options_.load_results_after_load_mh
if isfield(oo_load_mh.oo_,'convergence') if isfield(oo_load_mh.oo_,'convergence')
oo_.convergence=oo_load_mh.oo_.convergence; oo_.convergence=oo_load_mh.oo_.convergence;
@ -547,13 +465,13 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if ~(~isempty(options_.sub_draws) && options_.sub_draws==0) if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
if options_.bayesian_irf if options_.bayesian_irf
if error_flag if error_flag
error('Estimation::mcmc: I cannot compute the posterior IRFs!') error('%s: I cannot compute the posterior IRFs!',dispString)
end end
PosteriorIRF('posterior'); PosteriorIRF('posterior');
end end
if options_.moments_varendo if options_.moments_varendo
if error_flag if error_flag
error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!') error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
end end
if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results
[MetropolisFolder, info] = CheckPath('metropolis',M_.dname); [MetropolisFolder, info] = CheckPath('metropolis',M_.dname);
@ -578,16 +496,16 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
end end
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
if error_flag if error_flag
error('Estimation::mcmc: I cannot compute the posterior statistics!') error('%s: I cannot compute the posterior statistics!',dispString)
end end
if options_.order==1 && ~options_.particle.status if options_.order==1 && ~options_.particle.status
prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts
else else
error('Estimation::mcmc: Particle Smoothers are not yet implemented.') error('%s: Particle Smoothers are not yet implemented.',dispString)
end end
end end
else else
fprintf('Estimation:mcmc: sub_draws was set to 0. Skipping posterior computations.') fprintf('%s: sub_draws was set to 0. Skipping posterior computations.',dispString);
end end
xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
M_ = set_all_parameters(xparam1,estim_params_,M_); M_ = set_all_parameters(xparam1,estim_params_,M_);

149
matlab/dynare_estimation_init.m
View File

@ -80,26 +80,11 @@ if ~isfield(options_,'varobs')
error('VAROBS statement is missing!') error('VAROBS statement is missing!')
end end
% Checks on VAROBS
check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names);
% Set the number of observed variables. % Set the number of observed variables.
options_.number_of_observed_variables = length(options_.varobs); options_.number_of_observed_variables = length(options_.varobs);
% Check that each declared observed variable is also an endogenous variable.
for i = 1:options_.number_of_observed_variables
id = strmatch(options_.varobs{i}, M_.endo_names, 'exact');
if isempty(id)
error(['Unknown variable (' options_.varobs{i} ')!'])
end
end
% Check that a variable is not declared as observed more than once.
if length(unique(options_.varobs))<length(options_.varobs)
for i = 1:options_.number_of_observed_variables
if length(strmatch(options_.varobs{i},options_.varobs,'exact'))>1
error(['A variable cannot be declared as observed more than once (' options_.varobs{i} ')!'])
end
end
end
if options_.discretionary_policy if options_.discretionary_policy
if options_.order>1 if options_.order>1
error('discretionary_policy does not support order>1'); error('discretionary_policy does not support order>1');
@ -173,133 +158,7 @@ end
% Check that the provided mode_file is compatible with the current estimation settings. % Check that the provided mode_file is compatible with the current estimation settings.
if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0) && ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && sum(estim_params_.nvx+estim_params_.nvn+estim_params_.ncx+estim_params_.ncn+estim_params_.np)==0) && ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
number_of_estimated_parameters = length(xparam1); [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_);
mode_file = load(options_.mode_file);
if number_of_estimated_parameters>length(mode_file.xparam1)
% More estimated parameters than parameters in the mode file.
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your mode file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate ' int2str(number_of_estimated_parameters) ' parameters:'])
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded mode file (prior mean will be used, if possible).'])
else
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
if isempty(id)
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
if ~options_.mode_compute
% The posterior mode is not estimated.
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
error('Please remove the mode_file option.')
end
end
elseif number_of_estimated_parameters<length(mode_file.xparam1)
% Less estimated parameters than parameters in the mode file.
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your mode file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate only ' int2str(number_of_estimated_parameters) ' parameters:'])
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' deblank(bayestopt_.name(i,:)) ' is not present in the loaded mode file (prior mean will be used, if possible).'])
else
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
if isempty(id)
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
if ~options_.mode_compute
% The posterior mode is not estimated. If possible, fix the mode_file.
if isequal(length(xd),number_of_estimated_parameters)
disp('==> Fix mode file (remove unused parameters).')
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh = mode_file.hh(md,md);
end
else
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
end
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the mode_file option.')
end
end
else
% The number of declared estimated parameters match the number of parameters in the mode file.
% Check that the parameters in the mode file and according to the current mod file are identical.
if ~isfield(mode_file,'parameter_names')
disp(['The posterior mode file ' options_.mode_file ' has been generated using an older version of Dynare. It cannot be verified if it matches the present model. Proceed at your own risk.'])
mode_file.parameter_names=deblank(bayestopt_.name); %set names
end
if isequal(mode_file.parameter_names, bayestopt_.name)
xparam1 = mode_file.xparam1;
if isfield(mode_file,'hh')
hh = mode_file.hh;
end
else
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
% Check if this only an ordering issue or if the missing parameters can be initialized with the prior mean.
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)), mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded mode file.'])
else
xd = [xd; i];
md = [md; id];
end
end
if ~options_.mode_compute
% The posterior mode is not estimated
if isequal(length(xd), number_of_estimated_parameters)
% This is an ordering issue.
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh = mode_file.hh(md,md);
end
else
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
end
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh(xd,xd) = mode_file.hh(md,md);
end
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the mode_file option.')
end
end
end
end
skipline()
end end
%check for calibrated covariances before updating parameters %check for calibrated covariances before updating parameters
@ -627,7 +486,7 @@ if options_.load_results_after_load_mh
end end
if options_.mh_replic || options_.load_mh_file if options_.mh_replic || options_.load_mh_file
[current_options, options_, bayestopt_] = check_posterior_sampler_options([], options_, bounds, bayestopt_); [current_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = current_options; options_.posterior_sampler_options.current_options = current_options;
end end

2
matlab/get_the_name.m
View File

@ -32,7 +32,7 @@ function [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, options_)
%! @sp 2 %! @sp 2
%! @strong{This function is called by:} %! @strong{This function is called by:}
%! @sp 1 %! @sp 1
%! @ref{get_prior_info}, @ref{McMCDiagnostics}, @ref{mode_check}, @ref{PlotPosteriorDistributions}, @ref{plot_priors} %! @ref{get_prior_info}, @ref{mcmc_diagnostics}, @ref{mode_check}, @ref{PlotPosteriorDistributions}, @ref{plot_priors}
%! @sp 2 %! @sp 2
%! @strong{This function calls:} %! @strong{This function calls:}
%! @sp 1 %! @sp 1

47
matlab/initial_estimation_checks.m
View File

@ -161,49 +161,12 @@ if (any(BayesInfo.pshape >0 ) && DynareOptions.mh_replic) && DynareOptions.mh_n
error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.']) error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.'])
end end
old_steady_params=Model.params; %save initial parameters for check if steady state changes param values % check and display warnings if steady-state solves static model (except if diffuse_filter == 1) and if steady-state changes estimated parameters
[DynareResults.steady_state] = check_steady_state_changes_parameters(Model,EstimatedParameters,DynareResults,DynareOptions, [DynareOptions.diffuse_filter==0 DynareOptions.diffuse_filter==0] );
% % check if steady state solves static model (except if diffuse_filter == 1) % check and display warning if negative values of stderr or corr params are outside unit circle for Bayesian estimation
[DynareResults.steady_state, new_steady_params] = evaluate_steady_state(DynareResults.steady_state,[DynareResults.exo_steady_state; DynareResults.exo_det_steady_state],Model,DynareOptions,DynareOptions.diffuse_filter==0); if any(BayesInfo.pshape)
check_prior_stderr_corr(EstimatedParameters,BayesInfo);
if isfield(EstimatedParameters,'param_vals') && ~isempty(EstimatedParameters.param_vals)
%check whether steady state file changes estimated parameters
Model_par_varied=Model; %store Model structure
Model_par_varied.params(EstimatedParameters.param_vals(:,1))=Model_par_varied.params(EstimatedParameters.param_vals(:,1))*1.01; %vary parameters
[~, new_steady_params_2] = evaluate_steady_state(DynareResults.steady_state,[DynareResults.exo_steady_state; DynareResults.exo_det_steady_state],Model_par_varied,DynareOptions,DynareOptions.diffuse_filter==0);
changed_par_indices=find((old_steady_params(EstimatedParameters.param_vals(:,1))-new_steady_params(EstimatedParameters.param_vals(:,1))) ...
| (Model_par_varied.params(EstimatedParameters.param_vals(:,1))-new_steady_params_2(EstimatedParameters.param_vals(:,1))));
if ~isempty(changed_par_indices)
fprintf('\nThe steady state file internally changed the values of the following estimated parameters:\n')
disp(char(Model.param_names(EstimatedParameters.param_vals(changed_par_indices,1))))
fprintf('This will override the parameter values drawn from the proposal density and may lead to wrong results.\n')
fprintf('Check whether this is really intended.\n')
warning('The steady state file internally changes the values of the estimated parameters.')
end
end
if any(BayesInfo.pshape) % if Bayesian estimation
nvx=EstimatedParameters.nvx;
if nvx && any(BayesInfo.p3(1:nvx)<0)
warning('Your prior allows for negative standard deviations for structural shocks. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset=nvx;
nvn=EstimatedParameters.nvn;
if nvn && any(BayesInfo.p3(1+offset:offset+nvn)<0)
warning('Your prior allows for negative standard deviations for measurement error. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset = nvx+nvn;
ncx=EstimatedParameters.ncx;
if ncx && (any(BayesInfo.p3(1+offset:offset+ncx)<-1) || any(BayesInfo.p4(1+offset:offset+ncx)>1))
warning('Your prior allows for correlations between structural shocks larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end
offset = nvx+nvn+ncx;
ncn=EstimatedParameters.ncn;
if ncn && (any(BayesInfo.p3(1+offset:offset+ncn)<-1) || any(BayesInfo.p4(1+offset:offset+ncn)>1))
warning('Your prior allows for correlations between measurement errors larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end
end end
% display warning if some parameters are still NaN % display warning if some parameters are still NaN

2
matlab/marginal_density.m
View File

@ -48,7 +48,7 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws); TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws);
fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... '); fprintf('Estimation::marginal density: I''m computing the posterior mean and covariance... ');
[posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(); [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at_the_mode] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname);
MU = transpose(posterior_mean); MU = transpose(posterior_mean);
SIGMA = posterior_covariance; SIGMA = posterior_covariance;

4
matlab/metropolis_draw.m
View File

@ -72,8 +72,8 @@ if init
else else
if options_.sub_draws>NumberOfDraws*mh_nblck if options_.sub_draws>NumberOfDraws*mh_nblck
skipline() skipline()
disp(['Estimation::mcmc: The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!']) disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!'])
disp('Estimation::mcmc: You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).') disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
skipline() skipline()
xparams = 1; % xparams is interpreted as an error flag xparams = 1; % xparams is interpreted as an error flag
end end

14
matlab/pm3.m
View File

@ -24,7 +24,7 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
% See also the comment in posterior_sampler.m funtion. % See also the comment in posterior_sampler.m funtion.
% Copyright © 2007-2018 Dynare Team % Copyright © 2007-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -43,6 +43,8 @@ function pm3(n1,n2,ifil,B,tit1,tit2,tit3,tit_tex,names1,names2,name3,DirectoryNa
global options_ M_ oo_ global options_ M_ oo_
dispString = 'Estimation::mcmc';
nn = 3; nn = 3;
MaxNumberOfPlotsPerFigure = nn^2; % must be square MaxNumberOfPlotsPerFigure = nn^2; % must be square
varlist = names2; varlist = names2;
@ -76,7 +78,7 @@ HPD = zeros(2,n2,nvar);
if options_.estimation.moments_posterior_density.indicator if options_.estimation.moments_posterior_density.indicator
Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar); Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar);
end end
fprintf(['Estimation::mcmc: ' tit1 '\n']); fprintf(['%s: ' tit1 '\n'],dispString);
k = 0; k = 0;
filter_step_ahead_indicator=0; filter_step_ahead_indicator=0;
filter_covar_indicator=0; filter_covar_indicator=0;
@ -161,7 +163,7 @@ elseif filter_covar_indicator
oo_.FilterCovariance.post_deciles=post_deciles; oo_.FilterCovariance.post_deciles=post_deciles;
oo_.FilterCovariance.HPDinf=squeeze(hpd_interval(:,:,:,1)); oo_.FilterCovariance.HPDinf=squeeze(hpd_interval(:,:,:,1));
oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2)); oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2));
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); fprintf(['%s: ' tit1 ', done!\n'],dispString);
return return
elseif state_uncert_indicator elseif state_uncert_indicator
draw_dimension=4; draw_dimension=4;
@ -183,7 +185,7 @@ elseif state_uncert_indicator
oo_.Smoother.State_uncertainty.post_deciles=post_deciles; oo_.Smoother.State_uncertainty.post_deciles=post_deciles;
oo_.Smoother.State_uncertainty.HPDinf=squeeze(hpd_interval(:,:,:,1)); oo_.Smoother.State_uncertainty.HPDinf=squeeze(hpd_interval(:,:,:,1));
oo_.Smoother.State_uncertainty.HPDsup=squeeze(hpd_interval(:,:,:,2)); oo_.Smoother.State_uncertainty.HPDsup=squeeze(hpd_interval(:,:,:,2));
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); fprintf(['%s: ' tit1 ', done!\n'],dispString);
return return
end end
@ -280,7 +282,7 @@ else
end end
if strcmp(var_type,'_trend_coeff') || max(max(abs(Mean(:,:))))<=10^(-6) || all(all(isnan(Mean))) if strcmp(var_type,'_trend_coeff') || max(max(abs(Mean(:,:))))<=10^(-6) || all(all(isnan(Mean)))
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); fprintf(['%s: ' tit1 ', done!\n'],dispString);
return %not do plots return %not do plots
end end
%% %%
@ -378,4 +380,4 @@ if ~options_.nograph && ~options_.no_graph.posterior
end end
end end
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); fprintf(['%s: ' tit1 ', done!\n'],dispString);

12
matlab/posterior_sampler.m
View File

@ -53,6 +53,8 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
dispString = 'Estimation::mcmc';
vv = sampler_options.invhess; vv = sampler_options.invhess;
% Initialization of the sampler % Initialization of the sampler
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ... [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
@ -173,10 +175,10 @@ update_last_mh_history_file(MetropolisFolder, ModelName, record);
% Provide diagnostic output % Provide diagnostic output
skipline() skipline()
disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) fprintf('%s: Number of mh files: %u per block.\n', dispString, NewFile(1));
disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) fprintf('%s: Total number of generated files: %u.\n', dispString, NewFile(1)*nblck);
disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) fprintf('%s: Total number of iterations: %u.\n', dispString, (NewFile(1)-1)*MAX_nruns+irun-1);
disp(['Estimation::mcmc: Current acceptance ratio per chain: ']) fprintf('%s: Current acceptance ratio per chain:\n', dispString);
for i=1:nblck for i=1:nblck
if i<10 if i<10
disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%'])
@ -185,7 +187,7 @@ for i=1:nblck
end end
end end
if max(record.FunctionEvalPerIteration)>1 if max(record.FunctionEvalPerIteration)>1
disp(['Estimation::mcmc: Current function evaluations per iteration: ']) fprintf('%s: Current function evaluations per iteration:\n', dispString);
for i=1:nblck for i=1:nblck
if i<10 if i<10
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])

144
matlab/posterior_sampler_initialization.m
View File

@ -55,6 +55,8 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
dispString = 'Estimation::mcmc';
%Initialize outputs %Initialize outputs
ix2 = []; ix2 = [];
ilogpo2 = []; ilogpo2 = [];
@ -79,22 +81,22 @@ d = chol(vv);
if ~options_.load_mh_file && ~options_.mh_recover if ~options_.load_mh_file && ~options_.mh_recover
% Here we start a new Metropolis-Hastings, previous draws are discarded. % Here we start a new Metropolis-Hastings, previous draws are discarded.
if NumberOfBlocks > 1 if NumberOfBlocks > 1
disp('Estimation::mcmc: Multiple chains mode.') fprintf('%s: Multiple chains mode.\n',dispString);
else else
disp('Estimation::mcmc: One Chain mode.') fprintf('%s: One Chain mode.\n',dispString);
end end
% Delete old mh files if any... % Delete old mh files if any...
files = dir([BaseName '_mh*_blck*.mat']); files = dir([BaseName '_mh*_blck*.mat']);
if length(files) if length(files)
delete([BaseName '_mh*_blck*.mat']); delete([BaseName '_mh*_blck*.mat']);
disp('Estimation::mcmc: Old mh-files successfully erased!') fprintf('%s: Old mh-files successfully erased!\n',dispString);
end end
% Delete old Metropolis log file. % Delete old Metropolis log file.
file = dir([ MetropolisFolder '/metropolis.log']); file = dir([ MetropolisFolder '/metropolis.log']);
if length(file) if length(file)
delete([ MetropolisFolder '/metropolis.log']); delete([ MetropolisFolder '/metropolis.log']);
disp('Estimation::mcmc: Old metropolis.log file successfully erased!') fprintf('%s: Old metropolis.log file successfully erased!\n',dispString)
disp('Estimation::mcmc: Creation of a new metropolis.log file.') fprintf('%s: Creation of a new metropolis.log file.\n',dispString)
end end
fidlog = fopen([MetropolisFolder '/metropolis.log'],'w'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'w');
fprintf(fidlog,'%% MH log file (Dynare).\n'); fprintf(fidlog,'%% MH log file (Dynare).\n');
@ -116,8 +118,8 @@ if ~options_.load_mh_file && ~options_.mh_recover
%% check for proper filesep char in user defined paths %% check for proper filesep char in user defined paths
RecordFile0=strrep(RecordFile0,'\',filesep); RecordFile0=strrep(RecordFile0,'\',filesep);
if isempty(dir(RecordFile0)) if isempty(dir(RecordFile0))
disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_record option') fprintf('%s: Wrong value for mh_initialize_from_previous_mcmc_record option.\n',dispString);
error('Estimation::mcmc: path to record file is not found') error('%s: Path to record file is not found!',dispString)
else else
record0=load(RecordFile0); record0=load(RecordFile0);
end end
@ -133,31 +135,31 @@ if ~options_.load_mh_file && ~options_.mh_recover
record0=load_last_mh_history_file(MetropolisFolder0, ModelName0); record0=load_last_mh_history_file(MetropolisFolder0, ModelName0);
end end
if ~isnan(record0.MCMCConcludedSuccessfully) && ~record0.MCMCConcludedSuccessfully if ~isnan(record0.MCMCConcludedSuccessfully) && ~record0.MCMCConcludedSuccessfully
error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.') error('%s: You are trying to load an MCMC that did not finish successfully. Please use ''mh_recover''!',dispString);
end end
% mh_files = dir([ MetropolisFolder0 filesep ModelName0 '_mh*.mat']); % mh_files = dir([ MetropolisFolder0 filesep ModelName0 '_mh*.mat']);
% if ~length(mh_files) % if ~length(mh_files)
% error('Estimation::mcmc: I cannot find any MH file to load here!') % error('%s: I cannot find any MH file to load here!',dispString)
% end % end
disp('Estimation::mcmc: Initializing from past Metropolis-Hastings simulations...') fprintf('%s: Initializing from past Metropolis-Hastings simulations...\n',dispString);
disp(['Estimation::mcmc: Past MH path ' MetropolisFolder0 ]) fprintf('%s: Past MH path %s\n',dispString,MetropolisFolder0);
disp(['Estimation::mcmc: Past model name ' ModelName0 ]) fprintf('%s: Past model name %s\n', dispString, ModelName0);
fprintf(fidlog,' Loading initial values from previous MH\n'); fprintf(fidlog,' Loading initial values from previous MH\n');
fprintf(fidlog,' Past MH path: %s\n', MetropolisFolder0 ); fprintf(fidlog,' Past MH path: %s\n', MetropolisFolder0 );
fprintf(fidlog,' Past model name: %s\n', ModelName0); fprintf(fidlog,' Past model name: %s\n', ModelName0);
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
past_number_of_blocks = record0.Nblck; past_number_of_blocks = record0.Nblck;
if past_number_of_blocks ~= NumberOfBlocks if past_number_of_blocks ~= NumberOfBlocks
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') fprintf('%s: The specified number of blocks doesn''t match with the previous number of blocks!\n', dispString);
disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) fprintf('%s: You declared %u blocks, but the previous number of blocks was %u.\n', dispString, NumberOfBlocks, past_number_of_blocks);
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) fprintf('%s: I will run the Metropolis-Hastings with %u block.\n', dispString, past_number_of_blocks);
NumberOfBlocks = past_number_of_blocks; NumberOfBlocks = past_number_of_blocks;
options_.mh_nblck = NumberOfBlocks; options_.mh_nblck = NumberOfBlocks;
end end
if ~isempty(PriorFile0) if ~isempty(PriorFile0)
if isempty(dir(PriorFile0)) if isempty(dir(PriorFile0))
disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_prior option') fprintf('%s: Wrong value for mh_initialize_from_previous_mcmc_prior option.', dispString);
error('Estimation::mcmc: path to prior file is not found') error('%s: Path to prior file is not found!',dispString);
else else
bayestopt0 = load(PriorFile0); bayestopt0 = load(PriorFile0);
end end
@ -177,7 +179,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains
set_dynare_seed('default'); set_dynare_seed('default');
fprintf(fidlog,[' Initial values of the parameters:\n']); fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...') fprintf('%s: Searching for initial values...\n', dispString);
if ~options_.mh_initialize_from_previous_mcmc.status if ~options_.mh_initialize_from_previous_mcmc.status
ix2 = zeros(NumberOfBlocks,npar); ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(NumberOfBlocks,1); ilogpo2 = zeros(NumberOfBlocks,1);
@ -217,56 +219,54 @@ if ~options_.load_mh_file && ~options_.mh_recover
end end
init_iter = init_iter + 1; init_iter = init_iter + 1;
if init_iter > 100 && validate == 0 if init_iter > 100 && validate == 0
disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.']) fprintf('%s: I couldn''t get a valid initial value in 100 trials.\n', dispString);
if options_.nointeractive if options_.nointeractive
disp(['Estimation::mcmc: I reduce mh_init_scale_factor by 10 percent:']) fprintf('%s: I reduce ''mh_init_scale_factor'' by 10 percent:\n', dispString);
if isfield(options_,'mh_init_scale') if isfield(options_,'mh_init_scale')
options_.mh_init_scale = .9*options_.mh_init_scale; options_.mh_init_scale = .9*options_.mh_init_scale;
fprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.\n',options_.mh_init_scale) fprintf('%s: Parameter ''mh_init_scale'' is now equal to %f.\n',dispString,options_.mh_init_scale);
else else
options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor; options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor;
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is now equal to %f.\n',options_.mh_init_scale_factor) fprintf('%s: Parameter ''mh_init_scale_factor'' is now equal to %f.\n', dispString,options_.mh_init_scale_factor);
end end
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is now equal to %f.\n',options_.mh_init_scale_factor) fprintf('%s: Parameter ''mh_init_scale_factor'' is now equal to %f.\n', dispString,options_.mh_init_scale_factor);
else else
if isfield(options_,'mh_init_scale') if isfield(options_,'mh_init_scale')
disp(['Estimation::mcmc: You should reduce mh_init_scale...']) fprintf('%s: You should reduce mh_init_scale...\n',dispString);
fprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.\n',options_.mh_init_scale) fprintf('%s: Parameter ''mh_init_scale'' is equal to %f.\n',dispString,options_.mh_init_scale);
options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... '); options_.mh_init_scale_factor = input(sprintf('%s: Enter a new value... ',dispString));
else else
disp(['Estimation::mcmc: You should reduce mh_init_scale_factor...']) fprintf('%s: You should reduce ''mh_init_scale_factor''...\n',dispString);
fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is equal to %f.\n',options_.mh_init_scale_factor) fprintf('%s: Parameter ''mh_init_scale_factor'' is equal to %f.\n',dispString,options_.mh_init_scale_factor);
options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... '); options_.mh_init_scale_factor = input(sprintf('%s: Enter a new value... ',dispString));
end end
end end
trial = trial+1; trial = trial+1;
end end
end end
if trial > 10 && ~validate if trial > 10 && ~validate
disp(['Estimation::mcmc: I''m unable to find a starting value for block ' int2str(j)]) fprintf('%s: I''m unable to find a starting value for block %u.', dispString,j);
fclose(fidlog); fclose(fidlog);
return return
end end
end end
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
disp('Estimation::mcmc: Initial values found!') fprintf('%s: Initial values found!\n\n',dispString);
skipline()
else% Case 2: one chain (we start from the posterior mode) else% Case 2: one chain (we start from the posterior mode)
fprintf(fidlog,[' Initial values of the parameters:\n']); fprintf(fidlog,[' Initial values of the parameters:\n']);
candidate = transpose(xparam1(:));% candidate = transpose(xparam1(:));%
if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub) if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
ix2 = candidate; ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
disp('Estimation::mcmc: Initialization at the posterior mode.') fprintf('%s: Initialization at the posterior mode.\n\n',dispString);
skipline()
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']); fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
for i=1:length(ix2(1,:)) for i=1:length(ix2(1,:))
fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']); fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']);
end end
fprintf(fidlog,[' Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']); fprintf(fidlog,[' Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']);
else else
disp('Estimation::mcmc: Initialization failed...') fprintf('%s: Initialization failed...\n',dispString);
disp('Estimation::mcmc: The posterior mode lies outside the prior bounds.') fprintf('%s: The posterior mode lies outside the prior bounds.\n',dispString);
fclose(fidlog); fclose(fidlog);
return return
end end
@ -279,7 +279,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
% Delete the mh-history files % Delete the mh-history files
delete_mh_history_files(MetropolisFolder, ModelName); delete_mh_history_files(MetropolisFolder, ModelName);
% Create a new record structure % Create a new record structure
fprintf(['Estimation::mcmc: Write details about the MCMC... ']); fprintf('%s: Write details about the MCMC... ', dispString);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns); AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns; AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
record.Sampler = options_.posterior_sampler_options.posterior_sampling_method; record.Sampler = options_.posterior_sampler_options.posterior_sampling_method;
@ -312,8 +312,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
record.ProposalCovariance=d; record.ProposalCovariance=d;
fprintf('Ok!\n'); fprintf('Ok!\n');
id = write_mh_history_file(MetropolisFolder, ModelName, record); id = write_mh_history_file(MetropolisFolder, ModelName, record);
disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat']) fprintf('%s: Details about the MCMC are available in %s_mh_history_%u.mat\n\n', dispString,BaseName,id);
skipline()
fprintf(fidlog,[' CREATION OF THE MH HISTORY FILE!\n\n']); fprintf(fidlog,[' CREATION OF THE MH HISTORY FILE!\n\n']);
fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']); fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']); fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
@ -332,15 +331,15 @@ if ~options_.load_mh_file && ~options_.mh_recover
fclose(fidlog); fclose(fidlog);
elseif options_.load_mh_file && ~options_.mh_recover elseif options_.load_mh_file && ~options_.mh_recover
% Here we consider previous mh files (previous mh did not crash). % Here we consider previous mh files (previous mh did not crash).
disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...') fprintf('%s: I am loading past Metropolis-Hastings simulations...\n',dispString);
record=load_last_mh_history_file(MetropolisFolder, ModelName); record=load_last_mh_history_file(MetropolisFolder, ModelName);
if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully
error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.') error('%s: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.',dispString);
end end
record.MCMCConcludedSuccessfully=0; %reset indicator for this run record.MCMCConcludedSuccessfully=0; %reset indicator for this run
mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']); mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']);
if ~length(mh_files) if ~length(mh_files)
error('Estimation::mcmc: I cannot find any MH file to load here!') error('%s: I cannot find any MH file to load here!',dispString);
end end
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,'\n'); fprintf(fidlog,'\n');
@ -351,9 +350,9 @@ elseif options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
past_number_of_blocks = record.Nblck; past_number_of_blocks = record.Nblck;
if past_number_of_blocks ~= NumberOfBlocks if past_number_of_blocks ~= NumberOfBlocks
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') fprintf('%s: The specified number of blocks doesn''t match with the previous number of blocks!\n',dispString);
disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) fprintf('%s: You declared %u blocks, but the previous number of blocks was %u.\n', dispString,NumberOfBlocks,past_number_of_blocks);
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) fprintf('%s: I will run the Metropolis-Hastings with %u blocks.\n', dispString,past_number_of_blocks);
NumberOfBlocks = past_number_of_blocks; NumberOfBlocks = past_number_of_blocks;
options_.mh_nblck = NumberOfBlocks; options_.mh_nblck = NumberOfBlocks;
end end
@ -369,10 +368,10 @@ elseif options_.load_mh_file && ~options_.mh_recover
end end
ilogpo2 = record.LastLogPost; ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters; ix2 = record.LastParameters;
[d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d); [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d,dispString);
FirstBlock = 1; FirstBlock = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1); NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
fprintf('Estimation::mcmc: I am writing a new mh-history file... '); fprintf('%s: I am writing a new mh-history file... ',dispString);
record.MhDraws = [record.MhDraws;zeros(1,3)]; record.MhDraws = [record.MhDraws;zeros(1,3)];
NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber; NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile; NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
@ -387,28 +386,27 @@ elseif options_.load_mh_file && ~options_.mh_recover
record.InitialSeeds = record.LastSeeds; record.InitialSeeds = record.LastSeeds;
write_mh_history_file(MetropolisFolder, ModelName, record); write_mh_history_file(MetropolisFolder, ModelName, record);
fprintf('Done.\n') fprintf('Done.\n')
disp(['Estimation::mcmc: Ok. I have loaded ' int2str(NumberOfPreviousSimulations) ' simulations.']) fprintf('%s: Ok. I have loaded %u simulations.\n\n', dispString,NumberOfPreviousSimulations);
skipline()
fclose(fidlog); fclose(fidlog);
elseif options_.mh_recover elseif options_.mh_recover
% The previous metropolis-hastings crashed before the end! I try to recover the saved draws... % The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
disp('Estimation::mcmc: Recover mode!') fprintf('%s: Recover mode!\n',dispString);
record=load_last_mh_history_file(MetropolisFolder, ModelName); record=load_last_mh_history_file(MetropolisFolder, ModelName);
NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains. NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
options_.mh_nblck = NumberOfBlocks; options_.mh_nblck = NumberOfBlocks;
%% check consistency of options %% check consistency of options
if record.MhDraws(end,1)~=options_.mh_replic if record.MhDraws(end,1)~=options_.mh_replic
fprintf('\nEstimation::mcmc: You cannot specify a different mh_replic than in the chain you are trying to recover\n') fprintf('\n%s: You cannot specify a different mh_replic than in the chain you are trying to recover\n',dispString);
fprintf('Estimation::mcmc: I am resetting mh_replic to %u\n',record.MhDraws(end,1)) fprintf('%s: I am resetting mh_replic to %u\n',dispString,record.MhDraws(end,1));
options_.mh_replic=record.MhDraws(end,1); options_.mh_replic = record.MhDraws(end,1);
nruns = ones(NumberOfBlocks,1)*options_.mh_replic; nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
end end
if ~isnan(record.MAX_nruns(end,1)) %field exists if ~isnan(record.MAX_nruns(end,1)) %field exists
if record.MAX_nruns(end,1)~=MAX_nruns if record.MAX_nruns(end,1)~=MAX_nruns
fprintf('\nEstimation::mcmc: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n') fprintf('\n%s: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n',dispString);
fprintf('Estimation::mcmc: I am resetting MAX_nruns to %u\n',record.MAX_nruns(end,1)) fprintf('%s: I am resetting MAX_nruns to %u\n',dispString,record.MAX_nruns(end,1));
MAX_nruns=record.MAX_nruns(end,1); MAX_nruns=record.MAX_nruns(end,1);
end end
end end
@ -451,7 +449,7 @@ elseif options_.mh_recover
LastFileFullIndicator=1; LastFileFullIndicator=1;
end end
if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice') if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice')
[d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d); [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d,dispString);
end end
%% Now find out what exactly needs to be redone %% Now find out what exactly needs to be redone
% 1. Check if really something needs to be done % 1. Check if really something needs to be done
@ -464,10 +462,10 @@ elseif options_.mh_recover
% Quit if no crashed mcmc chain can be found as there are as many files as expected % Quit if no crashed mcmc chain can be found as there are as many files as expected
if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
if isnumeric(options_.parallel) if isnumeric(options_.parallel)
disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!') fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString);
disp(' You have to edit the mod file and remove the mh_recover option') fprintf(' You have to edit the mod file and remove the mh_recover option\n');
disp(' in the estimation command') fprintf(' in the estimation command.\n');
error('Estimation::mcmc: mh_recover option not required!') error('%s: mh_recover option not required!',dispString)
end end
end end
% 2. Something needs to be done; find out what % 2. Something needs to be done; find out what
@ -482,10 +480,10 @@ elseif options_.mh_recover
FBlock = zeros(NumberOfBlocks,1); FBlock = zeros(NumberOfBlocks,1);
while FirstBlock <= NumberOfBlocks while FirstBlock <= NumberOfBlocks
if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is not complete!']) fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock);
FBlock(FirstBlock)=1; FBlock(FirstBlock)=1;
else else
disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!']) fprintf('%s: Chain %u is complete!\n', dispString,FirstBlock);
end end
FirstBlock = FirstBlock+1; FirstBlock = FirstBlock+1;
end end
@ -537,8 +535,8 @@ elseif options_.mh_recover
record.InitialSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Unifor; record.InitialSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Unifor;
record.InitialSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Normal; record.InitialSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Normal;
else else
fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n') fprintf('%s: You are trying to recover a chain generated with an older Dynare version.\n',dispString);
fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n') fprintf('%s: I am using the default seeds to continue the chain.\n',dispString);
end end
if update_record if update_record
update_last_mh_history_file(MetropolisFolder, ModelName, record); update_last_mh_history_file(MetropolisFolder, ModelName, record);
@ -557,8 +555,8 @@ elseif options_.mh_recover
record.LastSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(ExpectedNumberOfMhFilesPerBlock)]).Unifor; record.LastSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(ExpectedNumberOfMhFilesPerBlock)]).Unifor;
record.LastSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(ExpectedNumberOfMhFilesPerBlock)]).Normal; record.LastSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(ExpectedNumberOfMhFilesPerBlock)]).Normal;
else else
fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n') fprintf('%s: You are trying to recover a chain generated with an older Dynare version.\n',dispString);
fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n') fprintf('%s: I am using the default seeds to continue the chain.\n',dispString);
end end
if isfield(loaded_results,'accepted_draws_this_chain') if isfield(loaded_results,'accepted_draws_this_chain')
record.AcceptanceRatio(FirstBlock)=loaded_results.accepted_draws_this_chain/record.MhDraws(end,1); record.AcceptanceRatio(FirstBlock)=loaded_results.accepted_draws_this_chain/record.MhDraws(end,1);
@ -572,32 +570,32 @@ elseif options_.mh_recover
FirstBlock = find(FBlock==1,1); FirstBlock = find(FBlock==1,1);
end end
function [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d) function [d,bayestopt_,record]=set_proposal_density_to_previous_value(record,options_,bayestopt_,d,dispString)
if ~options_.use_mh_covariance_matrix if ~options_.use_mh_covariance_matrix
if isfield(record,'ProposalCovariance') && isfield(record,'ProposalScaleVec') if isfield(record,'ProposalCovariance') && isfield(record,'ProposalScaleVec')
fprintf('Estimation::mcmc: Recovering the previous proposal density.\n') fprintf('%s: Recovering the previous proposal density.\n',dispString);
d=record.ProposalCovariance; d=record.ProposalCovariance;
bayestopt_.jscale=record.ProposalScaleVec; bayestopt_.jscale=record.ProposalScaleVec;
else else
if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice') if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice')
% pass through input d unaltered % pass through input d unaltered
if options_.mode_compute~=0 if options_.mode_compute~=0
fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by mode_compute\n.'); fprintf('%s: No stored previous proposal density found, continuing with the one implied by mode_compute.\n',dispString);
elseif ~isempty(options_.mode_file) elseif ~isempty(options_.mode_file)
fprintf('Estimation::mcmc: No stored previous proposal density found, continuing with the one implied by the mode_file\n.'); fprintf('%s: No stored previous proposal density found, continuing with the one implied by the mode_file.\n',dispString);
else else
error('Estimation::mcmc: No stored previous proposal density found, no mode-finding conducted, and no mode-file provided. I don''t know how to continue') error('%s: No stored previous proposal density found, no mode-finding conducted, and no mode-file provided. I don''t know how to continue!',dispString);
end end
end end
end end
else else
% pass through input d unaltered % pass through input d unaltered
fprintf('Estimation::mcmc: use_mh_covariance_matrix specified, continuing with proposal density implied by the previous MCMC run.\n.'); fprintf('%s: ''use_mh_covariance_matrix'' specified, continuing with proposal density implied by the previous MCMC run.\n',dispString);
end end
if isfield(record,'Sampler') if isfield(record,'Sampler')
if ~strcmp(record.Sampler,options_.posterior_sampler_options.posterior_sampling_method) if ~strcmp(record.Sampler,options_.posterior_sampler_options.posterior_sampling_method)
warning('Estimation::mcmc: The posterior_sampling_method %s selected differs from the %s of the original chain. This may create problems with the convergence diagnostics.',options_.posterior_sampler_options.posterior_sampling_method,record.Sampler) warning('%s: The posterior_sampling_method %s selected differs from the %s of the original chain. This may create problems with the convergence diagnostics.',dispString,options_.posterior_sampler_options.posterior_sampling_method,record.Sampler)
record.Sampler=options_.posterior_sampler_options.posterior_sampling_method; %update sampler used record.Sampler=options_.posterior_sampler_options.posterior_sampling_method; %update sampler used
end end
end end

4
matlab/set_prior.m
View File

@ -5,7 +5,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% INPUTS % INPUTS
% o estim_params_ [structure] characterizing parameters to be estimated. % o estim_params_ [structure] characterizing parameters to be estimated.
% o M_ [structure] characterizing the model. % o M_ [structure] characterizing the model.
% o options_ [structure] % o options_ [structure] characterizing the options.
% %
% OUTPUTS % OUTPUTS
% o xparam1 [double] vector of parameters to be estimated (initial values) % o xparam1 [double] vector of parameters to be estimated (initial values)
@ -18,7 +18,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None % None
% Copyright © 2003-2018 Dynare Team % Copyright © 2003-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %

69
matlab/utilities/estimation/check_hessian_at_the_mode.m Normal file
View File

@ -0,0 +1,69 @@
function check_hessian_at_the_mode(hessian_xparam1, xparam1, M_, estim_params_, options_, bounds)
% check_hessian_at_the_mode(hessian_xparam1, xparam1, M_, estim_params_, options_, bounds)
% -------------------------------------------------------------------------
% This function checks whether the hessian matrix at the mode is positive definite.
% =========================================================================
% INPUTS
% o hessian_xparam1: [matrix] hessian matrix at the mode
% o xparam1: [vector] vector of parameter values at the mode
% o M_: [structure] information about model
% o estim_params_: [structure] information about estimated parameters
% o options_: [structure] information about options
% o bounds: [structure] information about bounds
% -------------------------------------------------------------------------
% OUTPUTS
% none, displays a warning message if the hessian matrix is not positive definite
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_1.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
try
chol(hessian_xparam1);
catch
tol_bounds = 1.e-10;
skipline()
disp('OPTIMIZATION PROBLEM!')
disp(' (minus) the hessian matrix at the "mode" is not positive definite!')
disp('=> variance of the estimated parameters are not positive.')
disp('You should try to change the initial values of the parameters using')
disp('the estimated_params_init block, or use another optimization routine.')
params_at_bound = find(abs(xparam1-bounds.ub)<tol_bounds | abs(xparam1-bounds.lb)<tol_bounds);
if ~isempty(params_at_bound)
for ii=1:length(params_at_bound)
params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_);
end
disp_string=[params_at_bound_name{1,:}];
for ii=2:size(params_at_bound_name,1)
disp_string=[disp_string,', ',params_at_bound_name{ii,:}];
end
fprintf('\nThe following parameters are at the bound: %s\n', disp_string)
fprintf('Some potential solutions are:\n')
fprintf(' - Check your model for mistakes.\n')
fprintf(' - Check whether model and data are consistent (correct observation equation).\n')
fprintf(' - Shut off prior_trunc.\n')
fprintf(' - Change the optimization bounds.\n')
fprintf(' - Use a different mode_compute like 6 or 9.\n')
fprintf(' - Check whether the parameters estimated are identified.\n')
fprintf(' - Check prior shape (e.g. Inf density at bound(s)).\n')
fprintf(' - Increase the informativeness of the prior.\n')
end
warning('The results below are most likely wrong!');
end

163
matlab/utilities/estimation/check_mode_file.m Normal file
View File

@ -0,0 +1,163 @@
function [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_)
% function [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_)
% -------------------------------------------------------------------------
% Check that the provided mode_file is compatible with the current estimation settings.
% =========================================================================
% INPUTS
% o xparam1: [vector] current vector of parameter values at the mode
% o hh: [matrix] current hessian matrix at the mode
% o options_: [structure] information about options
% o bayestopt_: [structure] information about priors
% -------------------------------------------------------------------------
% OUTPUTS
% o xparam1: [vector] updated vector of parameter values at the mode
% o hh: [matrix] updated hessian matrix at the mode
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_init.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
number_of_estimated_parameters = length(xparam1);
mode_file = load(options_.mode_file);
if number_of_estimated_parameters>length(mode_file.xparam1)
% More estimated parameters than parameters in the mode_file.
skipline()
disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate ' int2str(number_of_estimated_parameters) ' parameters:'])
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).'])
else
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
if isempty(id)
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
if ~options_.mode_compute
% The mode is not estimated.
error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean or initialized values.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
error('Please remove the ''mode_file'' option.')
end
end
elseif number_of_estimated_parameters<length(mode_file.xparam1)
% Less estimated parameters than parameters in the mode_file.
skipline()
disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate only ' int2str(number_of_estimated_parameters) ' parameters:'])
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' deblank(bayestopt_.name(i,:)) ' is not present in the loaded ''mode_file'' (prior mean or initialized values will be used, if possible).'])
else
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
id = strmatch(mode_file.parameter_names{i},bayestopt_.name,'exact');
if isempty(id)
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
if ~options_.mode_compute
% The posterior mode is not estimated. If possible, fix the mode_file.
if isequal(length(xd),number_of_estimated_parameters)
disp('==> Fix mode_file (remove unused parameters).')
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh = mode_file.hh(md,md);
end
else
error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
end
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode_file using the prior mean or initialized values.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the ''mode_file'' option.')
end
end
else
% The number of declared estimated parameters match the number of parameters in the mode file.
% Check that the parameters in the mode file and according to the current mod file are identical.
if ~isfield(mode_file,'parameter_names')
disp(['The ''mode_file'' ' options_.mode_file ' has been generated using an older version of Dynare. It cannot be verified if it matches the present model. Proceed at your own risk.'])
mode_file.parameter_names=deblank(bayestopt_.name); %set names
end
if isequal(mode_file.parameter_names, bayestopt_.name)
xparam1 = mode_file.xparam1;
if isfield(mode_file,'hh')
hh = mode_file.hh;
end
else
skipline()
disp(['The ''mode_file'' ' options_.mode_file ' has been generated using another specification of the model or another model!'])
% Check if this is only an ordering issue or if the missing parameters can be initialized with the prior mean.
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)), mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded ''mode_file''.'])
else
xd = [xd; i];
md = [md; id];
end
end
if ~options_.mode_compute
% The mode is not estimated
if isequal(length(xd), number_of_estimated_parameters)
% This is an ordering issue.
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh = mode_file.hh(md,md);
end
else
error('Please change the ''mode_file'' option, the list of estimated parameters or set ''mode_compute''>0.')
end
else
% The mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode_file using the prior mean or initialized values.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh(xd,xd) = mode_file.hh(md,md);
end
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the ''mode_file'' option.')
end
end
end
end
skipline()

53
matlab/utilities/estimation/check_prior_stderr_corr.m Normal file
View File

@ -0,0 +1,53 @@
function check_prior_stderr_corr(estim_params_,bayestopt_)
% function check_prior_stderr_corr(estim_params_,bayestopt_)
% -------------------------------------------------------------------------
% Check if the prior allows for negative standard deviations and
% correlations larger than +-1. If so, issue a warning.
% =========================================================================
% INPUTS
% o estim_params_: [struct] information on estimated parameters
% o bayestopt_: [struct] information on priors
% -------------------------------------------------------------------------
% OUTPUTS
% none, but issues a warning if the prior allows for negative standard
% -------------------------------------------------------------------------
% This function is called by
% o initial_estimation_checks.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
nvx = estim_params_.nvx; % number of stderr parameters for structural shocks
if nvx && any(bayestopt_.p3(1:nvx)<0)
warning('Your prior allows for negative standard deviations for structural shocks. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset = nvx;
nvn = estim_params_.nvn; % number of stderr parameters for measurement errors
if nvn && any(bayestopt_.p3(1+offset:offset+nvn)<0)
warning('Your prior allows for negative standard deviations for measurement error. Due to working with variances, Dynare will be able to continue, but it is recommended to change your prior.')
end
offset = nvx+nvn;
ncx = estim_params_.ncx; % number of corr parameters for structural shocks
if ncx && (any(bayestopt_.p3(1+offset:offset+ncx)<-1) || any(bayestopt_.p4(1+offset:offset+ncx)>1))
warning('Your prior allows for correlations between structural shocks larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end
offset = nvx+nvn+ncx;
ncn = estim_params_.ncn; % number of corr parameters for measurement errors
if ncn && (any(bayestopt_.p3(1+offset:offset+ncn)<-1) || any(bayestopt_.p4(1+offset:offset+ncn)>1))
warning('Your prior allows for correlations between measurement errors larger than +-1 and will not integrate to 1 due to truncation. Please change your prior')
end

61
matlab/utilities/estimation/check_steady_state_changes_parameters.m Normal file
View File

@ -0,0 +1,61 @@
function [steady_state, info, steady_state_changes_parameters] = check_steady_state_changes_parameters(M_,estim_params_,oo_,options_,steadystate_check_flag_vec)
% function [steady_state, info, steady_state_changes_parameters] = check_steady_state_changes_parameters(M_,estim_params_,oo_,options_,steadystate_check_flag_vec)
% -------------------------------------------------------------------------
% Check if steady-state solves static model and if it changes estimated parameters
% =========================================================================
% INPUTS
% o M_: [struct] information on the model
% o estim_params_: [struct] information on estimated parameters
% o oo_: [struct] results
% o options_: [struct] information on options
% o steadystate_check_flag_vec: [vector] steadystate_check_flag for both checks (might be different in case of diffuse_filter)
% -------------------------------------------------------------------------
% OUTPUTS
% o steady_state: [vector] steady state
% o info: [scalar] 0 if steady state solves static model
% o steady_state_changes_parameters: [logical] true if steady state file changes estimated parameters
% -------------------------------------------------------------------------
% This function is called by
% o initial_estimation_checks.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
value_parameter_change = 1.01; % value with which parameters are slightly changed.
steady_state_changes_parameters = false; % initialize
% check if steady state solves static model
[steady_state, new_steady_params, info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,steadystate_check_flag_vec(1));
% check whether steady state file changes estimated parameters
if isfield(estim_params_,'param_vals') && ~isempty(estim_params_.param_vals)
old_steady_params = M_.params; % save initial parameters
M_par_varied = M_; % store Model structure
M_par_varied.params(estim_params_.param_vals(:,1)) = M_par_varied.params(estim_params_.param_vals(:,1))*value_parameter_change; % vary parameters
[~, new_steady_params_2] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_par_varied,options_,steadystate_check_flag_vec(2));
changed_par_indices = find((old_steady_params(estim_params_.param_vals(:,1))-new_steady_params(estim_params_.param_vals(:,1))) ...
| (M_par_varied.params(estim_params_.param_vals(:,1))-new_steady_params_2(estim_params_.param_vals(:,1))));
if ~isempty(changed_par_indices)
fprintf('\nThe steady state file internally changed the values of the following estimated parameters:\n')
disp(char(M_.param_names(estim_params_.param_vals(changed_par_indices,1))))
fprintf('This will override the parameter values and may lead to wrong results.\n')
fprintf('Check whether this is really intended.\n')
warning('The steady state file internally changes the values of the estimated parameters.')
steady_state_changes_parameters = true;
end
end

50
matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m Normal file
View File

@ -0,0 +1,50 @@
function check_varobs_are_endo_and_declared_once(varobs,endo_names)
% function check_varobs_are_endo_and_declared_once(varobs,endo_names)
% -------------------------------------------------------------------------
% Check that each declared observed variable:
% - is also an endogenous variable
% - is declared only once
% =========================================================================
% INPUTS
% o varobs: [cell] list of observed variables
% o endo_names: [cell] list of endogenous variables
% -------------------------------------------------------------------------
% OUTPUTS
% none, display an error message something is wrong with VAROBS
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_init.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
number_of_observed_variables = length(varobs);
for i = 1:number_of_observed_variables
if ~any(strcmp(varobs{i},endo_names))
error(['VAROBS: unknown variable (' varobs{i} ')!'])
end
end
% Check that a variable is not declared as observed more than once.
if length(unique(varobs))<length(varobs)
for i = 1:number_of_observed_variables
if sum(strcmp(varobs{i},varobs)) > 1
error(['VAROBS: a variable cannot be declared as observed more than once (' varobs{i} ')!'])
end
end
end

67
matlab/utilities/estimation/set_mcmc_jumping_covariance.m Normal file
View File

@ -0,0 +1,67 @@
function invhess = set_mcmc_jumping_covariance(invhess, xparam_nbr, MCMC_jumping_covariance, bayestopt_, stringForErrors)
% function invhess = set_mcmc_jumping_covariance(invhess, xparam_nbr, MCMC_jumping_covariance, bayestopt_, stringForErrors)
% -------------------------------------------------------------------------
% sets the jumping covariance matrix for the MCMC algorithm
% =========================================================================
% INPUTS
% o invhess: [matrix] already computed inverse of the hessian matrix
% o xparam_nbr: [integer] number of estimated parameters
% o MCMC_jumping_covariance: [string] name of option or file setting the jumping covariance matrix
% o bayestopt_: [struct] information on priors
% o stringForErrors: [string] string to be used in error messages
% -------------------------------------------------------------------------
% OUTPUTS
% o invhess: [matrix] jumping covariance matrix
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_1.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
switch MCMC_jumping_covariance
case 'hessian' % do nothing and use hessian from previous mode optimization
case 'prior_variance' % use prior variance
if any(isinf(bayestopt_.p2))
error('%s: Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.',stringForErrors);
else
hess = diag(1./(bayestopt_.p2.^2));
end
hsd = sqrt(diag(hess));
invhess = inv(hess./(hsd*hsd'))./(hsd*hsd');
case 'identity_matrix' % use identity
invhess = eye(xparam_nbr);
otherwise % user specified matrix in file
try
load(MCMC_jumping_covariance,'jumping_covariance')
hess = jumping_covariance;
catch
error(['%s: No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat!'],stringForErrors);
end
[nrow, ncol] = size(hess);
if ~isequal(nrow,ncol) && ~isequal(nrow,xparam_nbr) % check if square and right size
error(['%s: ''jumping_covariance'' matrix (loaded from ',options_.MCMC_jumping_covariance,'.mat) must be square and have ',num2str(xparam_nbr),' rows and columns!'],stringForErrors);
end
try % check for positive definiteness
chol(hess);
hsd = sqrt(diag(hess));
invhess = inv(hess./(hsd*hsd'))./(hsd*hsd');
catch
error('%s: Specified ''MCMC_jumping_covariance'' is not positive definite!',stringForErrors);
end
end

51
matlab/utilities/estimation/set_mcmc_prior_bounds.m Normal file
View File

@ -0,0 +1,51 @@
function bounds = set_mcmc_prior_bounds(xparam, bayestopt_, options_, stringForErrors)
% function bounds = set_mcmc_prior_bounds(xparam, bayestopt_, options_, stringForErrors)
% -------------------------------------------------------------------------
% Reset bounds as lower and upper bounds must only be operational during mode-finding
% =========================================================================
% INPUTS
% o xparam: [vector] vector of parameters
% o bayestopt_: [struct] information on priors
% o options_: [struct] Dynare options
% o stringForErrors: [string] string to be used in error messages
% -------------------------------------------------------------------------
% OUTPUTS
% o bounds: [struct] structure with fields lb and ub
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_1.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
bounds = prior_bounds(bayestopt_, options_.prior_trunc);
outside_bound_pars = find(xparam < bounds.lb | xparam > bounds.ub);
if ~isempty(outside_bound_pars)
for ii = 1:length(outside_bound_pars)
outside_bound_par_names{ii,1} = get_the_name(ii,0,M_,estim_params_,options_);
end
disp_string = [outside_bound_par_names{1,:}];
for ii = 2:size(outside_bound_par_names,1)
disp_string = [disp_string,', ',outside_bound_par_names{ii,:}];
end
if options_.prior_trunc > 0
error(['%s: Mode value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0!'],stringForErrors);
else
error(['%s: Mode value(s) of ', disp_string ,' are outside parameter bounds!'],stringForErrors);
end
end

48
matlab/utilities/estimation/tune_mcmc_mh_jscale_wrapper.m Normal file
View File

@ -0,0 +1,48 @@
function mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds, varargin)
% function mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds, varargin)
% -------------------------------------------------------------------------
% Wrapper to call the algorithm to tune the jumping scale parameter for the
% Metropolis-Hastings algorithm; currently only supports RW-MH algorithm.
% =========================================================================
% INPUTS
% o invhess: [matrix] jumping covariance matrix
% o options_: [struct] Dynare options
% o M_: [struct] Dynare model structure
% o objective_function: [function handle] objective function
% o xparam1: [vector] vector of estimated parameters at the mode
% o bounds: [struct] structure containing information on bounds
% o varargin: [cell] additional arguments to be passed to the objective function
% -------------------------------------------------------------------------
% OUTPUTS
% o mh_jscale: [double] tuned jumping scale parameter
% -------------------------------------------------------------------------
% This function is called by
% o dynare_estimation_1.m
% -------------------------------------------------------------------------
% Copyright © 2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================
% get invhess in case of use_mh_covariance_matrix
posterior_sampler_options_temp = options_.posterior_sampler_options.current_options;
posterior_sampler_options_temp.invhess = invhess;
posterior_sampler_options_temp = check_posterior_sampler_options(posterior_sampler_options_temp, M_.fname, M_.dname, options_);
opt = options_.mh_tune_jscale;
opt.rwmh = options_.posterior_sampler_options.rwmh;
mh_jscale = calibrate_mh_scale_parameter(objective_function, ...
posterior_sampler_options_temp.invhess, xparam1, [bounds.lb,bounds.ub], ...
opt, varargin{:});

3
tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc
View File

@ -193,9 +193,10 @@ INT(0)^2*INFL(1)^4; %redundant
@#endif @#endif
end; end;
M_.matched_moments_orig = M_.matched_moments;
method_of_moments( method_of_moments(
% Necessery options % Necessary options
mom_method = @{MoM_Method} % method of moments method; possible values: GMM|SMM mom_method = @{MoM_Method} % method of moments method; possible values: GMM|SMM
, datafile = 'AnScho_MoM_data_2.mat' % name of filename with data , datafile = 'AnScho_MoM_data_2.mat' % name of filename with data

2
tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod
View File

@ -152,7 +152,7 @@ end
method_of_moments( method_of_moments(
% Necessery options % Necessary options
mom_method = GMM % method of moments method; possible values: GMM|SMM mom_method = GMM % method of moments method; possible values: GMM|SMM
, datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data , datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data

2
tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod
View File

@ -139,7 +139,7 @@ end
@#for mommethod in ["SMM"] @#for mommethod in ["SMM"]
method_of_moments( method_of_moments(
% Necessery options % Necessary options
mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
, datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data

2
tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod
View File

@ -110,7 +110,7 @@ save('test_matrix.mat','weighting_matrix')
@#for mommethod in ["GMM", "SMM"] @#for mommethod in ["GMM", "SMM"]
method_of_moments( method_of_moments(
% Necessery options % Necessary options
mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
, datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data