diff --git a/doc/parallel/parallel.tex b/doc/parallel/parallel.tex index 849efe4c9..dd3ea6864 100644 --- a/doc/parallel/parallel.tex +++ b/doc/parallel/parallel.tex @@ -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 a number of procedures performed after the completion of Metropolis, that use the posterior MC sample: \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 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}). @@ -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_core.m"; \item the cycle looping over estimated parameters computing univariate diagnostics:\\ -\verb"McMCDiagnostics.m", \\ -\verb"McMCDiagnostics_core.m"; +\verb"mcmc_diagnostics.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"):\\ \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:\\ diff --git a/matlab/+mom/Jtest.m b/matlab/+mom/Jtest.m new file mode 100644 index 000000000..0cf41df5d --- /dev/null +++ b/matlab/+mom/Jtest.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/+mom/default_option_mom_values.m b/matlab/+mom/default_option_mom_values.m new file mode 100644 index 000000000..c0bc017d2 --- /dev/null +++ b/matlab/+mom/default_option_mom_values.m @@ -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 . +% ========================================================================= + + +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 \ No newline at end of file diff --git a/matlab/+mom/display_comparison_moments.m b/matlab/+mom/display_comparison_moments.m new file mode 100644 index 000000000..bcd1584ed --- /dev/null +++ b/matlab/+mom/display_comparison_moments.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/+mom/data_moments.m b/matlab/+mom/get_data_moments.m similarity index 90% rename from matlab/+mom/data_moments.m rename to matlab/+mom/get_data_moments.m index 8845f5508..0a2ced1c4 100644 --- a/matlab/+mom/data_moments.m +++ b/matlab/+mom/get_data_moments.m @@ -1,5 +1,5 @@ -function [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_) -% [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] = get_data_moments(data, oo_, matched_moments_, options_mom_) % This function computes the user-selected empirical moments from data % ========================================================================= % 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 % ------------------------------------------------------------------------- % This function is called by -% o mom.run.m -% o mom.objective_function.m +% o mom.run +% o mom.objective_function % ========================================================================= -% Copyright © 2020-2021 Dynare Team +% Copyright © 2020-2023 Dynare Team % % 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 powers = matched_moments_{jm,3}; 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( (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 @@ -66,7 +66,4 @@ for jm = 1:options_mom_.mom.mom_nbr end m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1); m_data(:,jm) = m_data_tmp; -end - - -end %function end \ No newline at end of file +end \ No newline at end of file diff --git a/matlab/+mom/matched_moments_block.m b/matlab/+mom/matched_moments_block.m new file mode 100644 index 000000000..a0b44c70b --- /dev/null +++ b/matlab/+mom/matched_moments_block.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/+mom/mode_compute_gmm_smm.m b/matlab/+mom/mode_compute_gmm_smm.m new file mode 100644 index 000000000..4b934bc51 --- /dev/null +++ b/matlab/+mom/mode_compute_gmm_smm.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/+mom/objective_function.m b/matlab/+mom/objective_function.m index 13a301d7d..419c7112b 100644 --- a/matlab/+mom/objective_function.m +++ b/matlab/+mom/objective_function.m @@ -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_) -% [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_] = 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 -% o xparam1: current value of estimated parameters as returned by set_prior() -% o Bounds: structure containing parameter bounds -% o oo_: structure for results -% o estim_params_: structure describing the estimated_parameters -% 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 xparam: [vector] current value of estimated parameters as returned by set_prior() +% o Bounds: [structure] containing parameter bounds +% o oo_: [structure] for results +% o estim_params_: [structure] describing the estimated_parameters +% 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 bayestopt_: [structure] information about the prior % ------------------------------------------------------------------------- % OUTPUTS -% o fval: 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 exit_flag: 0 if error, 1 if no error -% o df: analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only) -% o junk1: empty matrix required for optimizer interface (Hessian would go here) -% o oo_: structure containing the results with the following updated fields: -% - mom.model_moments [numMom x 1] vector with model moments -% - mom.Q value of the quadratic form of the moment difference -% - mom.model_moments_params_derivs -% [numMom x numParams] Jacobian matrix of derivatives of model_moments with respect to estimated parameters -% (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_) +% o fval: [double] value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly) +% o info: [vector] information on error codes and penalties +% o exit_flag: [double] flag for exit status (0 if error, 1 if no error) +% o df: [matrix] analytical jacobian of the moment difference (wrt paramters), currently for GMM only +% o junkHessian: [matrix] empty matrix required for optimizer interface (Hessian would typically go here) +% o oo_: [structure] results with the following updated fields: +% - oo_.mom.model_moments: [vector] model moments +% - oo_.mom.Q: [double] value of the quadratic form of the moment difference +% - oo_.mom.model_moments_params_derivs: [matrix] analytical jacobian of the model moments wrt estimated parameters (currently for GMM only) +% o M_: [structure] updated model structure % ------------------------------------------------------------------------- % This function is called by -% o mom.run.m -% o dynare_minimize_objective.m +% o mom.run +% o dynare_minimize_objective % ------------------------------------------------------------------------- % This function calls -% o check_bounds_and_definiteness_estimation -% o pruned_state_space_system -% o resol -% o set_all_parameters +% o check_bounds_and_definiteness_estimation +% o get_perturbation_params_derivs +% o mom.get_data_moments +% 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. % @@ -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 % along with Dynare. If not, see . -% ------------------------------------------------------------------------- -% 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 - if options_mom_.vector_output == 1 - if options_mom_.mom.penalized_estimator - df = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1)); +junkHessian = []; +df = []; % required to be empty by e.g. newrat +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + 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 - df = nan(size(oo_.mom.data_moments,1),length(xparam1)); + df = nan(1,length(xparam)); end - else - df = nan(1,length(xparam1)); end -else - df=[]; %required to be empty by e.g. newrat 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. -xparam1 = xparam1(:); - -M_ = set_all_parameters(xparam1, estim_params_, M_); - -[fval,info,exit_flag]=check_bounds_and_definiteness_estimation(xparam1, M_, estim_params_, Bounds); +xparam = xparam(:); +M_ = set_all_parameters(xparam, estim_params_, M_); +[fval,info,exit_flag] = check_bounds_and_definiteness_estimation(xparam, M_, estim_params_, Bounds); 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; end return 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 [dr, info, M_, oo_] = resol(0, M_, options_mom_, oo_); - % Return, with endogenous penalty when possible, if resol issues an error code if info(1) 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) == 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; info(4) = info(2); 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; end return @@ -118,31 +119,32 @@ if info(1) fval = Inf; info(4) = 0.1; 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; end return end end + +%-------------------------------------------------------------------------- +% GMM: Set up pruned state-space system and compute model moments +%-------------------------------------------------------------------------- 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 ) - indpmodel = []; %initialize index for model parameters + indpmodel = []; % initialize index for model parameters 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 - indpstderr=[]; %initialize index for stderr parameters + indpstderr=[]; % initialize index for stderr parameters 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 - indpcorr=[]; %initialize matrix for corr paramters + indpcorr=[]; % initialize matrix for corr paramters 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 - 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.') end 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; 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); - 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 - 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 - oo_.mom.model_moments = NaN(options_mom_.mom.mom_nbr,1); for jm = 1:size(M_.matched_moments,1) % First moments 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); 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,:); end end - % Second moments + % second moments if (sum(M_.matched_moments{jm,3}) == 2) - idx1 = (oo_.dr.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)) ); + idx1 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(1)) ); + idx2 = (oo_.mom.obs_var == find(oo_.dr.order_var==M_.matched_moments{jm,1}(2)) ); if nnz(M_.matched_moments{jm,2}) == 0 - % Covariance + % covariance if options_mom_.prefilter 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 ) @@ -186,8 +187,8 @@ if strcmp(options_mom_.mom.mom_method,'GMM') end end else - % 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 + % autocovariance + 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 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 ) @@ -204,24 +205,23 @@ if strcmp(options_mom_.mom.mom_method,'GMM') end end end +end -elseif strcmp(options_mom_.mom.mom_method,'SMM') - %------------------------------------------------------------------------------ - % 3. Compute Moments of the model solution for normal innovations - %------------------------------------------------------------------------------ - +%------------------------------------------------------------------------------ +% SMM: Compute Moments of the model solution for Gaussian innovations +%------------------------------------------------------------------------------ +if strcmp(options_mom_.mom.mom_method,'SMM') % 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)); - 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 = 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 % simulate series y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order); % provide meaningful penalty if data is nan or inf 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); else fval = Inf; @@ -229,73 +229,70 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM') info(1)=180; info(4) = 0.1; 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; end return end - - % Remove burn-in and focus on observables (note that y_sim is in declaration order) - y_sim = y_sim(oo_.dr.order_var(oo_.dr.obs_var) , end-options_mom_.mom.long+1:end)'; - + % 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)'; if ~all(diag(M_.H)==0) 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 - shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); %initialize + chol_S = chol(M_.H(i_ME,i_ME)); % decompose rest + 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; y_sim = y_sim+shock_mat; end - - % Remove mean if centered moments + % remove mean if centered moments if options_mom_.prefilter y_sim = bsxfun(@minus, y_sim, mean(y_sim,1)); 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 + %-------------------------------------------------------------------------- -% 4. Compute quadratic target function +% Compute quadratic target function %-------------------------------------------------------------------------- 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 options_mom_.mom.penalized_estimator - dxparam1 = eye(length(xparam1)); +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*moments_difference; + 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 - - for jp=1:length(xparam1) - dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp); - dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference; - - if options_mom_.vector_output == 1 % lsqnonlin requires vector output - if options_mom_.mom.penalized_estimator - df(:,jp)=[dresiduals;dxparam1(:,jp)./sqrt(diag(oo_.prior.variance))]; + if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian + if options_mom_.mom.penalized_estimator + dxparam1 = eye(length(xparam)); + end + for jp=1:length(xparam) + dmoments_difference = - oo_.mom.model_moments_params_derivs(:,jp); + dresiduals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*oo_.mom.Sw*dmoments_difference; + 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 - df(:,jp) = dresiduals; - end - else - df(:,jp) = dresiduals'*residuals + residuals'*dresiduals; - 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)); + df(:,jp) = dresiduals'*residuals + residuals'*dresiduals; + if options_mom_.mom.penalized_estimator + 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)); + end end end end end -end%main function end +end % main function end diff --git a/matlab/+mom/print_info_on_estimation_settings.m b/matlab/+mom/print_info_on_estimation_settings.m new file mode 100644 index 000000000..45dfe832a --- /dev/null +++ b/matlab/+mom/print_info_on_estimation_settings.m @@ -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 . +% ========================================================================= +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'); \ No newline at end of file diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m index 720006edb..1a57a9be4 100644 --- a/matlab/+mom/run.m +++ b/matlab/+mom/run.m @@ -1,68 +1,84 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_) -%function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_) +% function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_) % ------------------------------------------------------------------------- % This function performs a method of moments estimation with the following steps: -% Step 0: Check if required structures and options exist -% Step 1: - Prepare options_mom_ structure -% - Carry over options from the preprocessor -% - Initialize other options -% - Get variable orderings and state space representation -% Step 2: Checks and transformations for matched moments structure -% Step 3: Checks and transformations for estimated parameters, priors, and bounds -% Step 4: Checks and transformations for data -% Step 5: Checks for steady state at initial parameters -% Step 6: Checks for objective function at initial parameters -% Step 7: Iterated method of moments estimation -% Step 8: J-Test -% Step 9: Clean up +% o Checking if required structures and options exist +% o Preparing local options_mom_ structure +% o Checking the options and the compatibility of the settings +% o Initializations of variables, orderings and state space representation +% o Checks and transformations for matched moments structure +% o Checks and transformations for estimated parameters, priors, and bounds +% o Checks and transformations for data +% o Checks for objective function at initial parameters +% o GMM/SMM: iterated method of moments estimation +% o GMM/SMM: J-Test and fit of moments% +% o Display of results +% o Clean up % ------------------------------------------------------------------------- % This function is inspired by replication codes accompanied to the following papers: +% GMM/SMM: % o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49. % o Born, Pfeifer (2014): "Risk Matters: Comment", American Economic Review, 104(12):4231-4239. % o Mutschler (2018): "Higher-order statistics for DSGE models", Econometrics and Statistics, 6:44-56. % ========================================================================= % INPUTS -% o bayestopt_: [structure] information about priors -% o options_: [structure] information about global options -% o oo_: [structure] storage for results -% o estim_params_: [structure] information about estimated parameters -% o M_: [structure] information about model with -% o matched_moments: [cell] information about selected moments to match in estimation -% vars: matched_moments{:,1}); -% lead/lags: matched_moments{:,2}; -% powers: matched_moments{:,3}; -% o options_mom_: [structure] information about settings specified by the user +% o bayestopt_: [structure] information about priors +% o options_: [structure] information about global options +% o oo_: [structure] storage for results +% o estim_params_: [structure] information about estimated parameters +% o M_: [structure] information about model with +% o matched_moments: [cell] information about selected moments to match in GMM/SMM estimation +% vars: matched_moments{:,1}); +% lead/lags: matched_moments{:,2}; +% powers: matched_moments{:,3}; +% o options_mom_: [structure] information about settings specified by the user % ------------------------------------------------------------------------- % 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_) +% o oo_: [structure] storage for results (oo_) +% o options_mom_: [structure] information about all (user-specified and updated) settings used in estimation (options_mom_) +% o M_: [structure] updated information about model % ------------------------------------------------------------------------- % This function is called by % o driver.m % ------------------------------------------------------------------------- % This function calls -% o check_for_calibrated_covariances.m -% o check_prior_bounds.m -% o do_parameter_initialization.m -% o dynare_minimize_objective.m -% o evaluate_steady_state -% o get_all_parameters.m -% o get_matrix_entries_for_psd_check.m -% o makedataset.m -% o mom.check_plot.m -% o mom.data_moments.m -% o mom.objective_function.m +% o cellofchararraymaxlength +% o check_for_calibrated_covariances +% o check_prior_bounds +% o check_prior_stderr_corr +% o check_steady_state_changes_parameters +% o check_varobs_are_endo_and_declared_once +% o display_estimation_results_table +% o do_parameter_initialization +% o dyn_latex_table +% o dynare_minimize_objective +% o dyntable +% o get_all_parameters +% o get_dynare_random_generator_state +% o get_matrix_entries_for_psd_check +% o M_.fname '_prior_restrictions' +% o makedataset +% o mom.check_plot +% o mom.default_option_mom_values +% o mom.get_data_moments +% o mom.matched_moments_block +% o mom.objective_function % o mom.optimal_weighting_matrix -% o mom-standard_errors -% o plot_priors.m -% o print_info.m -% o prior_bounds.m -% o set_default_option.m -% o set_prior.m -% o set_state_space.m -% o set_all_parameters.m -% o test_for_deep_parameters_calibration.m - +% o mom.print_info_on_estimation_settings +% o mom.set_correct_bounds_for_stderr_corr +% o mom.standard_errors +% o plot_priors +% o prior_bounds +% o priordens +% o print_info +% o set_all_parameters +% o set_dynare_random_generator_state +% o set_prior +% o set_state_space +% o skipline +% o test_for_deep_parameters_calibration +% o warning_config +% ========================================================================= % Copyright © 2020-2023 Dynare Team % % This file is part of Dynare. @@ -80,20 +96,82 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . % ------------------------------------------------------------------------- -% Author(s): +% Maintaining Author(s): % o Willi Mutschler (willi@mutschler.eu) -% o Johannes Pfeifer (jpfeifer@uni-koeln.de) +% o Johannes Pfeifer (johannes.pfeifer@unibw.de) % ========================================================================= -%% TO DO LIST -% - [ ] add IRF matching -% - [ ] speed up pruned_state_space_system (by using doubling with old initial values, hardcoding zeros, other "tricks" used in e.g. nlma) -% - [ ] add option to use autocorrelations (we have useautocorr in identification toolbox already) -% - [ ] SMM with extended path -% - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox) -% - [ ] dirname option to save output to different directory not yet implemented -% - [ ] display scaled moments +% ------------------------------------------------------------------------- +% TO DO LISTS +% ------------------------------------------------------------------------- +% GENERAL +% - document all options in manual +% - document analytic_jacobian better +% - make endogenous_prior_restrictions work +% - dirname option to save output to different directory not yet implemented +% - create test for prior restrictions file +% - add mode_file option +% - implement penalty objective +% - test optimizers +% GMM/SMM +% - speed up pruned_state_space_system (by using doubling with old initial values, hardcoding zeros, other "tricks" used in e.g. nlma) +% - add option to use autocorrelations (we have useautocorr in identification toolbox already) +% - SMM with extended path +% - deal with measurement errors (once @wmutschl has implemented this in identification toolbox) +% - display scaled moments +% - enable first moments despite prefilter +% - do "true" Bayesian GMM and SMM not only penalized +fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method) + + +% ------------------------------------------------------------------------- +% checks if required structures exist +% ------------------------------------------------------------------------- +if isempty(estim_params_) % structure storing the info about estimated parameters in the estimated_params block + if ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0) + error('method_of_moments: You need to provide an ''estimated_params'' block!') + else + error('method_of_moments: The ''estimated_params'' block must not be empty!') + end +end +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + if ~isfield(M_,'matched_moments') || isempty(M_.matched_moments) % structure storing the moments used for GMM and SMM estimation + error('method_of_moments: You need to provide a ''matched_moments'' block for ''mom_method=%s''!',options_mom_.mom.mom_method) + end +end +if (~isempty(estim_params_.var_endo) || ~isempty(estim_params_.corrn)) && strcmp(options_mom_.mom.mom_method, 'GMM') + error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specify them as a structural shock!') +end +doBayesianEstimation = [estim_params_.var_exo(:,5); estim_params_.var_endo(:,5); estim_params_.corrx(:,6); estim_params_.corrn(:,6); estim_params_.param_vals(:,5)]; +if all(doBayesianEstimation~=0) + doBayesianEstimation = true; +elseif all(doBayesianEstimation==0) + doBayesianEstimation = false; +else + error('method_of_moments: Estimation must be either fully Frequentist or fully Bayesian. Maybe you forgot to specify a prior distribution!') +end +if ~isfield(options_,'varobs') + error('method_of_moments: VAROBS statement is missing!') +end +check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names); + + +% ------------------------------------------------------------------------- +% options_mom_ structure +% ------------------------------------------------------------------------- +% 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. +options_mom_ = mom.default_option_mom_values(options_mom_, options_, M_.dname, doBayesianEstimation); + + +% ------------------------------------------------------------------------- +% workarounds +% ------------------------------------------------------------------------- % The TeX option crashes MATLAB R2014a run with "-nodisplay" option % (as is done from the testsuite). % Since we can’t directly test whether "-nodisplay" has been passed, @@ -110,543 +188,263 @@ if isfield(options_mom_, 'TeX') && options_mom_.TeX && ~isoctave && matlab_ver_l warning('Disabling TeX option due to a bug in MATLAB R2014a with -nodisplay') options_mom_.TeX = false; end - -% ------------------------------------------------------------------------- -% Step 0: Check if required structures and options exist -% ------------------------------------------------------------------------- -if isempty(estim_params_) % structure storing the info about estimated parameters in the estimated_params block - if ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0) - error('method_of_moments: You need to provide an ''estimated_params'' block') - else - error('method_of_moments: The ''estimated_params'' block must not be empty') - end -end -if ~isfield(M_,'matched_moments') || isempty(M_.matched_moments) % structure storing the moments used for the method of moments estimation - error('method_of_moments: You need to provide a ''matched_moments'' block') -end -if ~isempty(bayestopt_) && any(bayestopt_.pshape==0) && any(bayestopt_.pshape~=0) - error('method_of_moments: Estimation must be either fully classical or fully Bayesian. Maybe you forgot to specify a prior distribution.') -end - -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.\n') -else - options_mom_.logged_steady_state = 0; - options_mom_.loglinear = false; -end - -fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method) - -% ------------------------------------------------------------------------- -% Step 1a: Prepare options_mom_ structure -% ------------------------------------------------------------------------- -% 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. - -% Method of Moments estimation options that can be set by the user in the mod file, otherwise default values are provided if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.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 weight - options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results - 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. - 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,'se_tolx',1e-5); % step size for numerical computation of standard errors - 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 higher-order - % Checks for perturbation order - if options_mom_.order < 1 - error('method_of_moments: The order of the Taylor approximation cannot be 0!') +% temporary workaround for https://git.dynare.org/Dynare/dseries/-/issues/51 + if options_mom_.xls_sheet~=1 + evalin('base','options_.xls_sheet=options_mom_.xls_sheet'); end + if ~isempty(options_mom_.xls_range) + evalin('base','options_.xls_range=options_mom_.xls_range'); + end +end + + +% ------------------------------------------------------------------------- +% checks on settings +% ------------------------------------------------------------------------- +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + if numel(options_mom_.nobs) > 1 + error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''!'); + end + if numel(options_mom_.first_obs) > 1 + error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''first_obs''!'); + end +end +if options_mom_.order < 1 + error('method_of_moments: The order of the Taylor approximation cannot be 0!') +end +if options_mom_.order > 2 + fprintf('Dynare will use ''k_order_solver'' as the order>2\n'); + options_mom_.k_order_solver = true; end if strcmp(options_mom_.mom.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 if options_mom_.mom.simulation_multiple < 1 - fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 5.\n') + fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 7.\n') options_mom_.mom.simulation_multiple = 7; end end if strcmp(options_mom_.mom.mom_method,'GMM') - % Check for pruning with GMM at higher order + % require pruning with GMM at higher order if options_mom_.order > 1 && ~options_mom_.pruning fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n'); options_mom_.pruning = true; end if options_mom_.order > 3 - error('method_of_moments: perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM.\n'); + error('method_of_moments: Perturbation orders higher than 3 are not implemented for GMM estimation, try using SMM!'); end end -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. -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 -% initialize flag to compute derivs in objective function (needed for GMM with either analytic_standard_errors or analytic_jacobian ) -options_mom_.mom.compute_derivs = false; - -% General options that can be set by the user in the mod file, otherwise default values are provided -options_mom_ = set_default_option(options_mom_,'dirname',M_.dname); % specify directory in which to store estimation output [not yet working] -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_,'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 -options_mom_ = set_default_option(options_mom_,'TeX',false); % print TeX tables and graphics -options_mom_ = set_default_option(options_mom_,'verbosity',false); % print TeX tables and graphics - -% Data and model options that can be set by the user in the mod file, otherwise default values are provided -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 -options_mom_ = set_default_option(options_mom_,'xls_range',''); % range of data in Excel sheet -% temporary workaround for https://git.dynare.org/Dynare/dseries/-/issues/51 -if options_mom_.xls_sheet~=1 - evalin('base','options_.xls_sheet=options_mom_.xls_sheet'); -end -if ~isempty(options_mom_.xls_range) - evalin('base','options_.xls_range=options_mom_.xls_range'); +if options_mom_.mom.analytic_jacobian && ~strcmp(options_mom_.mom.mom_method,'GMM') + options_mom_.mom.analytic_jacobian = false; + fprintf('\n''analytic_jacobian'' option will be dismissed as it only works with GMM.\n'); end -% Recursive estimation and forecast are not supported -if numel(options_mom_.nobs)>1 - error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''.'); -end -if numel(options_mom_.first_obs)>1 - error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''first_obs''.'); -end - -% Optimization options that can be set by the user in the mod file, otherwise default values are provided -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 -if (isoctave && user_has_octave_forge_package('optim')) || (~isoctave && user_has_matlab_license('optimization_toolbox')) - options_mom_ = set_default_option(options_mom_,'mode_compute',13); % specifies lsqnonlin as default optimizer for minimization of moments distance -else - options_mom_ = set_default_option(options_mom_,'mode_compute',4); % specifies csminwel as fallback default option for minimization of moments distance -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 -% Check plot options that can be set by the user in the mod file, otherwise default values are provided -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) - -% Numerical algorithms options that can be set by the user in the mod file, otherwise default values are provided -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_,'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 -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. -options_mom_ = set_default_option(options_mom_,'trust_region_initial_step_bound_factor',1); % used in dynare_solve for trust_region -if options_mom_.order > 2 - fprintf('Dynare will use ''k_order_solver'' as the order>2\n'); - options_mom_.k_order_solver = true; -end % ------------------------------------------------------------------------- -% Step 1b: Options that are set by the preprocessor and need to be carried over +% initializations % ------------------------------------------------------------------------- +% create output directories to store results +CheckPath('method_of_moments',M_.dname); +CheckPath('graphs',options_mom_.dirname); +% initialize options that might change +options_mom_.mom.compute_derivs = false; % flag to compute derivs in objective function (might change for GMM with either analytic_standard_errors or analytic_jacobian (dependent on optimizer)) +options_mom_.mom.vector_output = false; % specifies whether the objective function returns a vector +% decision rule +oo_.dr = set_state_space(oo_.dr,M_,options_mom_); % get state-space representation +oo_.mom.obs_var = []; % create index of observed variables in DR order +for i = 1:options_mom_.obs_nbr + oo_.mom.obs_var = [oo_.mom.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))]; +end -% options related to VAROBS -if ~isfield(options_,'varobs') - error('method_of_moments: VAROBS statement is missing!') -else - options_mom_.varobs = options_.varobs; % observable variables in declaration order - options_mom_.obs_nbr = length(options_mom_.varobs); % number of observed variables - % Check that each declared observed variable is also an endogenous variable - for i = 1:options_mom_.obs_nbr - if ~any(strcmp(options_mom_.varobs{i},M_.endo_names)) - error(['method_of_moments: Unknown variable (' options_mom_.varobs{i} ')!']) - end + +% ------------------------------------------------------------------------- +% matched_moments: checks and transformations +% ------------------------------------------------------------------------- +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + M_.matched_moments = mom.matched_moments_block(M_.matched_moments, options_mom_.mom.mom_method); + % Check if both prefilter and first moments were specified + first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))'; + if options_mom_.prefilter && ~isempty(first_moment_indicator) + fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block.\n'); + M_.matched_moments(first_moment_indicator,:)=[]; %remove first moments entries end - - % Check that a variable is not declared as observed more than once - if length(unique(options_mom_.varobs))1 - error(['method_of_moments: A variable cannot be declared as observed more than once (' options_mom_.varobs{i} ')!']) - end - end + options_mom_.mom.mom_nbr = size(M_.matched_moments,1); + % Get maximum lag number for autocovariances/autocorrelations + options_mom_.ar = max(cellfun(@max,M_.matched_moments(:,2))) - min(cellfun(@min,M_.matched_moments(:,2))); + % Check that only observed variables are involved in moments + not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),oo_.mom.obs_var); + if ~isempty(not_observed_variables) + skipline; + error('method_of_moments: You specified moments involving %s, but it is not a varobs!',M_.endo_names{oo_.dr.order_var(not_observed_variables)}) end end -% options related to variable declarations -if isfield(options_,'trend_coeffs') - error('method_of_moments: %s does not allow for trend in data',options_mom_.mom.mom_method) -end - -% options related to estimated_params and estimated_params_init -options_mom_.use_calibration_initialization = options_.use_calibration_initialization; - -% options 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; - -% options related to steady command -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 related to dataset -options_mom_.dataset = options_.dataset; -options_mom_.initial_period = options_.initial_period; - -% options related to endogenous prior restrictions are not supported -options_mom_.endogenous_prior_restrictions.irf = {}; -options_mom_.endogenous_prior_restrictions.moment = {}; -if ~isempty(options_.endogenous_prior_restrictions.irf) && ~isempty(options_.endogenous_prior_restrictions.moment) - fprintf('Endogenous prior restrictions are not supported yet and will be skipped.\n') -end % ------------------------------------------------------------------------- -% Step 1c: Options related to optimizers +% estimated parameters: checks and transformations on values, priors, bounds % ------------------------------------------------------------------------- -% mode_compute = 1, 3, 7, 11, 102, 11, 13 -% nothing to be done -% mode_compute = 2 -options_mom_.saopt = options_.saopt; -% mode_compute = 4 -options_mom_.csminwel = options_.csminwel; -% mode_compute = 5 -options_mom_.newrat = options_.newrat; -options_mom_.gstep = options_.gstep; -% mode_compute = 6 -options_mom_.gmhmaxlik = options_.gmhmaxlik; -options_mom_.mh_jscale = options_.mh_jscale; -% mode_compute = 8 -options_mom_.simplex = options_.simplex; -% mode_compute = 9 -options_mom_.cmaes = options_.cmaes; -% mode_compute = 10 -options_mom_.simpsa = options_.simpsa; -% mode_compute = 12 -options_mom_.particleswarm = options_.particleswarm; -% mode_compute = 101 -options_mom_.solveopt = options_.solveopt; - -options_mom_.gradient_method = options_.gradient_method; -options_mom_.gradient_epsilon = options_.gradient_epsilon; -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_.vector_output= false; % specifies whether the objective function returns a vector - -optimizer_vec=[options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)]; % at each stage one can possibly use different optimizers sequentially - -analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; %these are currently supported, see to-do list - -% ------------------------------------------------------------------------- -% Step 1d: Other options that need to be initialized -% ------------------------------------------------------------------------- -options_mom_.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m -options_mom_.figures.textwidth = 0.8; %needed by plot_priors.m -options_mom_.ramsey_policy = 0; % needed by evaluate_steady_state -options_mom_ = set_default_option(options_mom_,'debug',false); %neeeded by e.g. check_plot -options_mom_.risky_steadystate = false; %needed by resol -options_mom_.threads = options_.threads; %needed by resol -options_mom_.jacobian_flag = true; -options_mom_.gstep = options_.gstep; - -% options_mom.dsge_var = false; %needed by check_list_of_variables -% options_mom.bayesian_irf = false; %needed by check_list_of_variables -% options_mom.moments_varendo = false; %needed by check_list_of_variables -% options_mom.smoother = false; %needed by check_list_of_variables -% options_mom.filter_step_ahead = []; %needed by check_list_of_variables -% options_mom.forecast = 0; -%options_mom_ = set_default_option(options_mom_,'endo_vars_for_moment_computations_in_estimation',[]); - -% ------------------------------------------------------------------------- -% Step 1e: Get variable orderings and state space representation -% ------------------------------------------------------------------------- -oo_.dr = set_state_space(oo_.dr,M_,options_mom_); -% Get index of observed variables in DR order -oo_.dr.obs_var = []; -for i=1:options_mom_.obs_nbr - oo_.dr.obs_var = [oo_.dr.obs_var; find(strcmp(options_mom_.varobs{i}, M_.endo_names(oo_.dr.order_var)))]; -end - -% ------------------------------------------------------------------------- -% Step 2: Checks and transformations for matched_moments structure -% ------------------------------------------------------------------------- -M_.matched_moments_orig = M_.matched_moments; %save original structure - -% higher-order product moments not supported yet for GMM -if strcmp(options_mom_.mom.mom_method, 'GMM') && any(cellfun(@sum,M_.matched_moments(:,3))> 2) - error('method_of_moments: GMM does not yet support product moments higher than 2. Change row %d in ''matched_moments'' block.',jm); -end - -% check for duplicate moment conditions -for jm=1:size(M_.matched_moments,1) - % expand powers to vector of ones - if any(M_.matched_moments{jm,3}>1) - tmp1=[]; tmp2=[]; tmp3=[]; - for jjm=1:length(M_.matched_moments{jm,3}) - tmp1 = [tmp1 repmat(M_.matched_moments{jm,1}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ]; - tmp2 = [tmp2 repmat(M_.matched_moments{jm,2}(jjm),[1 M_.matched_moments{jm,3}(jjm)]) ]; - tmp3 = [tmp3 repmat(1,[1 M_.matched_moments{jm,3}(jjm)]) ]; - end - M_.matched_moments{jm,1} = tmp1; - M_.matched_moments{jm,2} = tmp2; - M_.matched_moments{jm,3} = tmp3; - end - % shift time structure to focus only on lags - M_.matched_moments{jm,2} = M_.matched_moments{jm,2} - max(M_.matched_moments{jm,2}); - % sort such that t=0 variable comes first - [M_.matched_moments{jm,2},idx_sort] = sort(M_.matched_moments{jm,2},'descend'); - M_.matched_moments{jm,1} = M_.matched_moments{jm,1}(idx_sort); - M_.matched_moments{jm,3} = M_.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,M_.matched_moments(:,3))'; -UniqueMomIdx = []; -for jpow = unique(powers) - idx1 = find(powers==jpow); - [~,idx2] = unique(cell2mat(M_.matched_moments(idx1,:)),'rows'); - UniqueMomIdx = [UniqueMomIdx idx1(idx2)]; -end - -% remove duplicate elements -DuplicateMoms = setdiff(1:size(M_.matched_moments_orig,1),UniqueMomIdx); -if ~isempty(DuplicateMoms) - fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows:\n %s.\n',num2str(DuplicateMoms)) - fprintf('Dynare will continue with remaining moment conditions\n'); -end - -if strcmp(options_mom_.mom.mom_method, 'SMM') - % for SMM we can keep the original structure but get rid of duplicate moments - M_.matched_moments = M_.matched_moments_orig(sort(UniqueMomIdx),:); -elseif strcmp(options_mom_.mom.mom_method, 'GMM') - % for GMM we use the transformed matched_moments structure - M_.matched_moments = M_.matched_moments(sort(UniqueMomIdx),:); -end - -% Check if both prefilter and first moments were specified -first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,M_.matched_moments(:,3)))'; -if options_mom_.prefilter && ~isempty(first_moment_indicator) - fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block.\n'); - M_.matched_moments(first_moment_indicator,:)=[]; %remove first moments entries -end -options_mom_.mom.mom_nbr = size(M_.matched_moments,1); - -% Get maximum lag number for autocovariances/autocorrelations -options_mom_.ar = max(cellfun(@max,M_.matched_moments(:,2))) - min(cellfun(@min,M_.matched_moments(:,2))); - -%check that only observed variables are involved in moments -not_observed_variables=setdiff(oo_.dr.inv_order_var([M_.matched_moments{:,1}]),oo_.dr.obs_var); -if ~isempty(not_observed_variables) - error('\nmethod_of_moments: You specified moments involving %s, but it is not a varobs.',M_.endo_names{oo_.dr.order_var(not_observed_variables)}) -end - -% ------------------------------------------------------------------------- -% Step 3: Checks and transformations for estimated parameters, priors, and bounds -% ------------------------------------------------------------------------- - -% Set priors and bounds over the estimated parameters +% set priors and bounds over the estimated parameters [xparam0, estim_params_, bayestopt_, lb, ub, M_] = set_prior(estim_params_, M_, options_mom_); +number_of_estimated_parameters = length(xparam0); +hessian_xparam0 = []; % initialize hessian -% Check measurement errors -if (estim_params_.nvn || estim_params_.ncn) && strcmp(options_mom_.mom.mom_method, 'GMM') - error('method_of_moments: GMM estimation does not support measurement error(s) yet. Please specifiy them as a structural shock.') +% check if enough moments for estimation +if strcmp(options_mom_.mom.mom_method, 'GMM') || strcmp(options_mom_.mom.mom_method, 'SMM') + if options_mom_.mom.mom_nbr < length(xparam0) + skipline; + error('method_of_moments: There must be at least as many moments as parameters for a %s estimation!',options_mom_.mom.mom_method); + end + skipline(2); end -% Check if enough moments for estimation -if options_mom_.mom.mom_nbr < length(xparam0) - fprintf('\n'); - error('method_of_moments: We must have at least as many moments as parameters for a method of moments estimation.') -end -fprintf('\n\n') - -% Check if a _prior_restrictions.m file exists +% check if a _prior_restrictions.m file exists if exist([M_.fname '_prior_restrictions.m'],'file') options_mom_.prior_restrictions.status = 1; options_mom_.prior_restrictions.routine = str2func([M_.fname '_prior_restrictions']); end -bayestopt_laplace=bayestopt_; - -% Check on specified priors and penalized estimation -if any(bayestopt_laplace.pshape > 0) % prior specified - if ~options_mom_.mom.penalized_estimator - fprintf('\nPriors were specified, but the penalized_estimator-option was not set.\n') - fprintf('Dynare sets penalized_estimator to 1. Conducting %s with penalty.\n',options_mom_.mom.mom_method) - options_mom_.mom.penalized_estimator=1; - end - if any(setdiff([0;bayestopt_laplace.pshape],[0,3])) - fprintf('\nNon-normal priors specified. %s with penalty uses a Laplace type of approximation.\n',options_mom_.mom.mom_method) - fprintf('Only 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_laplace.pshape~=3; - bayestopt_laplace.pshape(non_normal_priors) = 3; - bayestopt_laplace.p3(non_normal_priors) = -Inf*ones(sum(non_normal_priors),1); - bayestopt_laplace.p4(non_normal_priors) = Inf*ones(sum(non_normal_priors),1); - bayestopt_laplace.p6(non_normal_priors) = bayestopt_laplace.p1(non_normal_priors); - bayestopt_laplace.p7(non_normal_priors) = bayestopt_laplace.p2(non_normal_priors); - bayestopt_laplace.p5(non_normal_priors) = bayestopt_laplace.p1(non_normal_priors); - end - if any(isinf(bayestopt_laplace.p2)) %find infinite variance priors - inf_var_pars=bayestopt_laplace.name(isinf(bayestopt_laplace.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'); +% check on specified priors and penalized estimation (which uses Laplace approximated priors) +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + bayestopt_orig = bayestopt_; + if any(bayestopt_.pshape > 0) % prior specified + if ~options_mom_.mom.penalized_estimator + fprintf('\nPriors were specified, but the penalized_estimator-option was not set.'); + fprintf('\nDynare sets penalized_estimator to 1. Conducting %s with penalty.\n',options_mom_.mom.mom_method); + options_mom_.mom.penalized_estimator = 1; end + bayestopt_ = mom.transform_prior_to_laplace_prior(bayestopt_); end end -% Check for calibrated covariances before updating parameters +% check for calibrated covariances before updating parameters estim_params_ = check_for_calibrated_covariances(xparam0,estim_params_,M_); -% Checks on parameter calibration and initialization -xparam1_calib = get_all_parameters(estim_params_,M_); %get calibrated parameters -if ~any(isnan(xparam1_calib)) %all estimated parameters are calibrated - estim_params_.full_calibration_detected=1; +% checks on parameter calibration and initialization +xparam_calib = get_all_parameters(estim_params_,M_); % get calibrated parameters +if ~any(isnan(xparam_calib)) % all estimated parameters are calibrated + estim_params_.full_calibration_detected = true; else - estim_params_.full_calibration_detected=0; + estim_params_.full_calibration_detected = false; end -if options_mom_.use_calibration_initialization %set calibration as starting values - if ~isempty(bayestopt_laplace) && all(bayestopt_laplace.pshape==0) && any(all(isnan([xparam1_calib xparam0]),2)) - error('method_of_moments: When using the use_calibration option with %s without prior, the parameters must be explicitly initialized.',options_mom_.mom.mom_method) +if options_mom_.use_calibration_initialization % set calibration as starting values + if ~isempty(bayestopt_) && ~doBayesianEstimation && any(all(isnan([xparam_calib xparam0]),2)) + error('method_of_moments: When using the use_calibration option with %s without prior, the parameters must be explicitly initialized!',options_mom_.mom.mom_method); else - [xparam0,estim_params_]=do_parameter_initialization(estim_params_,xparam1_calib,xparam0); %get explicitly initialized parameters that have precedence over calibrated values + [xparam0,estim_params_] = do_parameter_initialization(estim_params_,xparam_calib,xparam0); % get explicitly initialized parameters that have precedence over calibrated values end end -% Check initialization -if ~isempty(bayestopt_laplace) && all(bayestopt_laplace.pshape==0) && any(isnan(xparam0)) - error('method_of_moments: %s without penalty requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block ',options_mom_.mom.mom_method) +% check initialization +if ~isempty(bayestopt_) && ~doBayesianEstimation && any(isnan(xparam0)) + error('method_of_moments: Frequentist %s requires all estimated parameters to be initialized, either in an estimated_params or estimated_params_init-block!',options_mom_.mom.mom_method); end -% Set and check parameter bounds -if ~isempty(bayestopt_laplace) && any(bayestopt_laplace.pshape > 0) - % Plot prior densities +% set and check parameter bounds +if ~isempty(bayestopt_) && doBayesianEstimation + % plot prior densities if ~options_mom_.nograph && options_mom_.plot_priors - plot_priors(bayestopt_,M_,estim_params_,options_mom_) - plot_priors(bayestopt_laplace,M_,estim_params_,options_mom_,'Laplace approximated priors') + if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + plot_priors(bayestopt_orig,M_,estim_params_,options_mom_,'Original priors'); % only for visual inspection (not saved to disk, because overwritten in next call to plot_priors) + plot_priors(bayestopt_,M_,estim_params_,options_mom_,'Laplace approximated priors'); + clear('bayestopt_orig'); % make sure stale structure cannot be used + end end - % Set prior bounds - Bounds = prior_bounds(bayestopt_laplace, options_mom_.prior_trunc); + % set prior bounds + Bounds = prior_bounds(bayestopt_, options_mom_.prior_trunc); Bounds.lb = max(Bounds.lb,lb); Bounds.ub = min(Bounds.ub,ub); -else % estimated parameters but no declared priors - % No priors are declared so Dynare will estimate the parameters - % with inequality constraints for the parameters. +else + % no priors are declared so Dynare will estimate the parameters with + % Frequentist methods using inequality constraints for the parameters Bounds.lb = lb; Bounds.ub = ub; if options_mom_.mom.penalized_estimator fprintf('Penalized estimation turned off as you did not declare priors\n') options_mom_.mom.penalized_estimator = 0; - end + end end -% Set correct bounds for standard deviations and corrrelations -param_of_interest=(1:length(xparam0))'<=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:length(xparam0))'> estim_params_.nvx+estim_params_.nvn & (1:length(xparam0))'1 & param_of_interest); -Bounds.lb(LB_below_minus_1)=-1; -Bounds.ub(UB_above_1)=1; -clear('bayestopt_','LB_below_0','LB_below_minus_1','UB_above_1','param_of_interest');%make sure stale structure cannot be used +% set correct bounds for standard deviations and correlations +Bounds = mom.set_correct_bounds_for_stderr_corr(estim_params_,Bounds); -% Test if initial values of the estimated parameters are all between the prior lower and upper bounds +% test if initial values of the estimated parameters are all between the prior lower and upper bounds if options_mom_.use_calibration_initialization try - check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_laplace) + check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_); catch last_error fprintf('Cannot use parameter values from calibration as they violate the prior bounds.') rethrow(last_error); end else - check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_laplace) + check_prior_bounds(xparam0,Bounds,M_,estim_params_,options_mom_,bayestopt_); end -estim_params_= get_matrix_entries_for_psd_check(M_,estim_params_); +% check for positive definiteness +estim_params_ = get_matrix_entries_for_psd_check(M_,estim_params_); -% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file). +% set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file) M_.sigma_e_is_diagonal = true; if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances') M_.sigma_e_is_diagonal = false; end -% storing prior parameters in MoM info structure for penalized minimization -oo_.prior.mean = bayestopt_laplace.p1; -oo_.prior.variance = diag(bayestopt_laplace.p2.^2); +% storing prior parameters in results +oo_.mom.prior.mean = bayestopt_.p1; +oo_.mom.prior.mode = bayestopt_.p5; +oo_.mom.prior.variance = diag(bayestopt_.p2.^2); +oo_.mom.prior.hyperparameters.first = bayestopt_.p6; +oo_.mom.prior.hyperparameters.second = bayestopt_.p7; -% Set all parameters +% set all parameters M_ = set_all_parameters(xparam0,estim_params_,M_); -%provide warning if there is NaN in parameters +% provide warning if there is NaN in parameters test_for_deep_parameters_calibration(M_); -% ------------------------------------------------------------------------- -% Step 4: Checks and transformations for data -% ------------------------------------------------------------------------- +if doBayesianEstimation + % warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle + check_prior_stderr_corr(estim_params_,bayestopt_); -% Check if datafile has same name as mod file -[~,name,~] = fileparts(options_mom_.datafile); -if strcmp(name,M_.fname) - error('method_of_moments: Data-file and mod-file are not allowed to have the same name. Please change the name of the data file.') + % check value of prior density + [~,~,~,info]= priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); + if any(info) + fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',bayestopt_.name{info,1}) + error('The initial value of the prior is -Inf!') + end end + +% ------------------------------------------------------------------------- +% datafile: checks and transformations +% ------------------------------------------------------------------------- % Build dataset -dataset_ = makedataset(options_mom_); - -% set options for old interface from the ones for new interface -if ~isempty(dataset_) - options_mom_.nobs = dataset_.nobs; +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + % Check if datafile has same name as mod file + [~,name,~] = fileparts(options_mom_.datafile); + if strcmp(name,M_.fname) + error('method_of_moments: ''datafile'' and mod file are not allowed to have the same name; change the name of the ''datafile''!') + end + dataset_ = makedataset(options_mom_); + % set options for old interface from the ones for new interface + if ~isempty(dataset_) + options_mom_.nobs = dataset_.nobs; + end + % Check length of data for estimation of second moments + if options_mom_.ar > options_mom_.nobs+1 + error('method_of_moments: Dataset is too short to compute higher than first moments!'); + end + % Provide info on data moments handling + fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment.\n'); + % Get data moments for the method of moments + [oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_); + if ~isreal(dataset_.data) + error('method_of_moments: The data moments contain complex values!') + end end -% Check length of data for estimation of second moments -if options_mom_.ar > options_mom_.nobs+1 - error('method_of_moments: Data set is too short to compute second moments'); -end -% Provide info on data moments handling -fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment\n'); - -% Get data moments for the method of moments -[oo_.mom.data_moments, oo_.mom.m_data] = mom.data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_); - -% Get shock series for SMM and set variance correction factor +% ------------------------------------------------------------------------- +% SMM: Get shock series fand set variance correction factor +% ------------------------------------------------------------------------- if strcmp(options_mom_.mom.mom_method,'SMM') options_mom_.mom.long = round(options_mom_.mom.simulation_multiple*options_mom_.nobs); options_mom_.mom.variance_correction_factor = (1+1/options_mom_.mom.simulation_multiple); @@ -659,7 +457,7 @@ if strcmp(options_mom_.mom.mom_method,'SMM') [state_u,state_n] = get_dynare_random_generator_state; %get state for later resetting set_dynare_random_generator_state(options_mom_.mom.seed,options_mom_.mom.seed); temp_shocks = randn(options_mom_.mom.long+options_mom_.mom.burnin,M_.exo_nbr); - temp_shocks_ME = randn(options_mom_.mom.long,length(M_.H)); + temp_shocks_ME = randn(options_mom_.mom.long,length(M_.H)); set_dynare_random_generator_state(state_u,state_n); %reset state for later resetting end if options_mom_.mom.bounded_shock_support == 1 @@ -675,84 +473,60 @@ if strcmp(options_mom_.mom.mom_method,'SMM') end end + % ------------------------------------------------------------------------- -% Step 5: checks for steady state at initial parameters +% checks for steady state at initial parameters % ------------------------------------------------------------------------- -% setting steadystate_check_flag option +% check if steady state solves static model and if steady-state changes estimated parameters if options_mom_.steadystate.nocheck - steadystate_check_flag = 0; + steadystate_check_flag_vec = [0 1]; else - steadystate_check_flag = 1; + steadystate_check_flag_vec = [1 1]; end - -old_steady_params=M_.params; %save initial parameters for check if steady state changes param values -% Check steady state at initial model parameter values -[oo_.steady_state, new_steady_params, info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_mom_,steadystate_check_flag); +[oo_.steady_state, info, steady_state_changes_parameters] = check_steady_state_changes_parameters(M_, estim_params_, oo_, options_mom_, steadystate_check_flag_vec); if info(1) - fprintf('\nmethod_of_moments: The steady state at the initial parameters cannot be computed.\n') + fprintf('\nThe steady state at the initial parameters cannot be computed.\n') print_info(info, 0, options_mom_); end - -% check whether steady state file changes estimated parameters -if isfield(estim_params_,'param_vals') && ~isempty(estim_params_.param_vals) - Model_par_varied=M_; %store M_ structure - - Model_par_varied.params(estim_params_.param_vals(:,1))=Model_par_varied.params(estim_params_.param_vals(:,1))*1.01; %vary parameters - [~, new_steady_params_2] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],Model_par_varied,options_mom_,true); - - changed_par_indices=find((old_steady_params(estim_params_.param_vals(:,1))-new_steady_params(estim_params_.param_vals(:,1))) ... - | (Model_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 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.') - if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors - fprintf('For analytical standard errors, the parameter-Jacobians of the dynamic model and of the steady-state will be computed numerically\n'), - fprintf('(re-set options_mom_.analytic_derivation_mode= -2)'), - options_mom_.analytic_derivation_mode= -2; - end - end +if steady_state_changes_parameters && strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors + fprintf('For analytical standard errors, the parameter-Jacobians of the dynamic model and of the steady-state will be computed numerically,\n'); + fprintf('because the steady-state changes estimated parameters. Option ''analytic_derivation_mode'' reset to -2.'); + options_mom_.analytic_derivation_mode = -2; end % display warning if some parameters are still NaN test_for_deep_parameters_calibration(M_); -% If steady state of observed variables is non zero, set noconstant equal 0 -if all(abs(oo_.steady_state(oo_.dr.order_var(oo_.dr.obs_var)))<1e-9) - options_mom_.noconstant = 0; %identifying the constant based on just the initial parameter value is not feasible -else - options_mom_.noconstant = 0; -end % ------------------------------------------------------------------------- -% Step 6: checks for objective function at initial parameters +% checks for objective function at initial parameters % ------------------------------------------------------------------------- objective_function = str2func('mom.objective_function'); - try - % Check for NaN or complex values of moment-distance-funtion evaluated - % at initial parameters and identity weighting matrix - oo_.mom.Sw = eye(options_mom_.mom.mom_nbr); - tic_id = tic; - [fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, M_, options_mom_); - elapsed_time = toc(tic_id); + % Check for NaN or complex values of moment-distance-funtion evaluated at initial parameters + if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM') + oo_.mom.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix + end + tic_id = tic; + [fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, M_, options_mom_, bayestopt_); + elapsed_time = toc(tic_id); if isnan(fval) - error('method_of_moments: The initial value of the objective function is NaN') + error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!') elseif imag(fval) - error('method_of_moments: The initial value of the objective function is complex') + error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!') end if info(1) > 0 disp('method_of_moments: Error in computing the objective function for initial parameter values') print_info(info, options_mom_.noprint, options_mom_) end - fprintf('Initial value of the moment objective function with %4.1f times identity weighting matrix: %6.4f \n\n', options_mom_.mom.weighting_matrix_scaling_factor, fval); + fprintf('Initial value of the moment objective function'); + if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM') + fprintf(' with %4.1f times identity weighting matrix', options_mom_.mom.weighting_matrix_scaling_factor); + end + fprintf(': %6.4f \n\n', fval); fprintf('Time required to compute objective function once: %5.4f seconds \n', elapsed_time); - -catch last_error% if check fails, provide info on using calibration if present +catch last_error % if check fails, provide info on using calibration if present if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values skipline(1); fprintf('There was an error in computing the moments for initial parameter values.\n') @@ -765,264 +539,60 @@ catch last_error% if check fails, provide info on using calibration if present rethrow(last_error); end -% ------------------------------------------------------------------------- -% Step 7a: Method of moments estimation: print some info -% ------------------------------------------------------------------------- -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 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 - -for i = 1:length(optimizer_vec) - if i == 1 - str = '- optimizer (mode_compute'; - else - str = ' (additional_optimizer_steps'; - end - switch 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 - error('\nmethod_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(optimizer_vec{i}) - fprintf('\n %s=%s): user-defined',str,optimizer_vec{i}); - else - error('method_of_moments: Unknown optimizer, please contact the developers ') - end - end - if options_mom_.silent_optimizer - fprintf(' (silent)'); - end - if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{i},analytic_jacobian_optimizers) - fprintf(' (using analytical Jacobian)'); - end -end -fprintf('\n - perturbation order: %d', options_mom_.order) -if options_mom_.order > 1 && options_mom_.pruning - fprintf(' (with pruning)') -end -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); -fprintf('\n - number of parameters: %d', length(xparam0)); -fprintf('\n\n'); % ------------------------------------------------------------------------- -% Step 7b: Iterated method of moments estimation +% print some info to console % ------------------------------------------------------------------------- -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 +mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters); -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(optimizer_vec) - options_mom_.current_optimizer = optimizer_vec{optim_iter}; - if 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 optimizer_vec{optim_iter}==13 - options_mom_.vector_output = true; - else - options_mom_.vector_output = false; - end - if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{optim_iter},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, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],... - Bounds, oo_, estim_params_, M_, options_mom_); - if options_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_laplace,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; - % Update M_ and DynareResults (in particular to get oo_.mom.model_moments) - M_ = set_all_parameters(xparam1,estim_params_,M_); +% ------------------------------------------------------------------------- +% GMM/SMM: iterated estimation +% ------------------------------------------------------------------------- +if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM') + % compute mode + [xparam1, oo_, Woptflag] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_, M_, options_mom_, estim_params_, bayestopt_, Bounds); + % compute standard errors at mode + options_mom_.mom.vector_output = false; % make sure flag is reset + M_ = set_all_parameters(xparam1,estim_params_,M_); % update M_ and DynareResults (in particular to get oo_.mom.model_moments) if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_standard_errors - options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag + options_mom_.mom.compute_derivs = true; % for GMM we compute derivatives analytically in the objective function with this flag end - [fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); + [~, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_); % compute model moments and oo_.mom.model_moments_params_derivs options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization - SE = mom.standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag); - + % mode check plots + if options_mom_.mode_check.status + mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_,... + Bounds, oo_, estim_params_, M_, options_mom_); + end +end + + +% ------------------------------------------------------------------------- +% display final estimation results +% ------------------------------------------------------------------------- +if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM') % Store results in output structure - oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_laplace,oo_.mom,prior_dist_names,sprintf('%s (STAGE %u)',options_mom_.mom.mom_method,stage_iter),sprintf('%s_stage_%u',lower(options_mom_.mom.mom_method),stage_iter)); + oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,options_mom_.mom.mom_method,lower(options_mom_.mom.mom_method)); + % J test + oo_ = mom.Jtest(xparam1, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, dataset_.nobs); + % display comparison of model moments and data moments + mom.display_comparison_moments(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments); end + % ------------------------------------------------------------------------- -% Step 8: J test +% clean up % ------------------------------------------------------------------------- -if options_mom_.mom.mom_nbr > length(xparam1) - %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, xparam1, Bounds, oo_j, estim_params_, M_, options_mom_); - 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 = dataset_.nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor; - oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam1); - 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 - -% ------------------------------------------------------------------------- -% Step 9: Display estimation results -% ------------------------------------------------------------------------- -title = ['Comparison of 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=[oo_.mom.data_moments oo_.mom.model_moments ]; -dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7); -if options_mom_.TeX - dyn_latex_table(M_, options_mom_, title, ['comparison_moments_', options_mom_.mom.mom_method], headers, labels_TeX, data_mat, cellofchararraymaxlength(labels)+2, 10, 7); -end - -if options_mom_.mode_check.status - mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,... - Bounds, oo_, estim_params_, M_, options_mom_) -end - fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method) -% ------------------------------------------------------------------------- -% Step 9: Clean up -% ------------------------------------------------------------------------- %reset warning state warning_config; if isoctave && isfield(options_, 'prior_restrictions') && ... - isfield(options_.prior_restrictions, 'routine') + isfield(options_mom_.prior_restrictions, 'routine') % Octave crashes if it tries to save function handles (to the _results.mat file) % See https://savannah.gnu.org/bugs/?43215 - options_.prior_restrictions.routine = []; -end + options_mom_.prior_restrictions.routine = []; +end \ No newline at end of file diff --git a/matlab/+mom/set_correct_bounds_for_stderr_corr.m b/matlab/+mom/set_correct_bounds_for_stderr_corr.m new file mode 100644 index 000000000..d65e9473b --- /dev/null +++ b/matlab/+mom/set_correct_bounds_for_stderr_corr.m @@ -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 . +% ========================================================================= + +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)'1 & param_of_interest); +Bounds.lb(LB_below_minus_1) = -1; +Bounds.ub(UB_above_1) = 1; \ No newline at end of file diff --git a/matlab/+mom/transform_prior_to_laplace_prior.m b/matlab/+mom/transform_prior_to_laplace_prior.m new file mode 100644 index 000000000..f2cdf7cd2 --- /dev/null +++ b/matlab/+mom/transform_prior_to_laplace_prior.m @@ -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 . +% ========================================================================= +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 \ No newline at end of file diff --git a/matlab/CheckPath.m b/matlab/CheckPath.m index 794c2005b..f32f8fa22 100644 --- a/matlab/CheckPath.m +++ b/matlab/CheckPath.m @@ -12,7 +12,7 @@ function [DirectoryName, info] = CheckPath(type,dname) % SPECIAL REQUIREMENTS % none -% Copyright © 2005-2017 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -31,7 +31,7 @@ function [DirectoryName, info] = CheckPath(type,dname) info = 0; -DirectoryName = [ dname '/' type ]; +DirectoryName = [ dname filesep type ]; if ~isdir(dname) % Make sure there isn't a file with the same name, see trac ticket #47 diff --git a/matlab/CutSample.m b/matlab/CutSample.m index e29f010d5..937ff413e 100644 --- a/matlab/CutSample.m +++ b/matlab/CutSample.m @@ -33,6 +33,8 @@ function CutSample(M_, options_, estim_params_) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +dispString = 'Estimation::mcmc'; + % Get the path to the metropolis files. MetropolisFolder = CheckPath('metropolis',M_.dname); @@ -47,7 +49,7 @@ npar=size(record.LastParameters,2); mh_files = dir([ MetropolisFolder ,filesep, M_.fname '_mh*.mat' ]); 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 TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); @@ -71,9 +73,9 @@ end % Save updated mh-history file. update_last_mh_history_file(MetropolisFolder, ModelName, record); -fprintf('Estimation::mcmc: Total number of MH draws per chain: %d.\n',TotalNumberOfMhDraws); -fprintf('Estimation::mcmc: Total number of generated MH files: %d.\n',TotalNumberOfMhFiles); -fprintf('Estimation::mcmc: I''ll use mh-files %d to %d.\n',FirstMhFile,TotalNumberOfMhFiles); -fprintf('Estimation::mcmc: In MH-file number %d I''ll start at line %d.\n',FirstMhFile,FirstLine); -fprintf('Estimation::mcmc: Finally I keep %d draws per chain.\n',TotalNumberOfMhDraws-FirstDraw+1); +fprintf('%s: Total number of MH draws per chain: %d.\n',dispString,TotalNumberOfMhDraws); +fprintf('%s: Total number of generated MH files: %d.\n',dispString,TotalNumberOfMhFiles); +fprintf('%s: I''ll use mh-files %d to %d.\n',dispString,FirstMhFile,TotalNumberOfMhFiles); +fprintf('%s: In MH-file number %d I''ll start at line %d.\n',dispString,FirstMhFile,FirstLine); +fprintf('%s: Finally I keep %d draws per chain.\n',dispString,TotalNumberOfMhDraws-FirstDraw+1); skipline() \ No newline at end of file diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index 007fb89f5..b5356fd7b 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -16,7 +16,7 @@ function PosteriorIRF(type) % functions associated with it(the _core1 and _core2). % See also the comments posterior_sampler.m funtion. -% Copyright © 2006-2018 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -36,6 +36,8 @@ function PosteriorIRF(type) global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info +dispString = 'Estimation::mcmc'; + % Set the number of periods if isempty(options_.irf) || ~options_.irf options_.irf = 40; @@ -287,7 +289,7 @@ if options_.TeX end end -fprintf('Estimation::mcmc: Posterior (dsge) IRFs...\n'); +fprintf('%s: Posterior (dsge) IRFs...\n',dispString); tit = M_.exo_names; kdx = 0; @@ -327,7 +329,7 @@ if MAX_nirfs_dsgevar VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr); DistribIRFdsgevar = zeros(options_.irf,9,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; kdx = 0; for file = 1:NumberOfIRFfiles_dsgevar @@ -457,4 +459,4 @@ if ~options_.nograph && ~options_.no_graph.posterior end -fprintf('Estimation::mcmc: Posterior IRFs, done!\n'); +fprintf('%s: Posterior IRFs, done!\n',dispString); diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m index d518abedb..d00d1913d 100644 --- a/matlab/check_posterior_sampler_options.m +++ b/matlab/check_posterior_sampler_options.m @@ -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, 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_) % initialization of posterior samplers % % INPUTS % posterior_sampler_options: posterior sampler options +% fname: name of the mod-file +% dname: name of directory with metropolis folder % options_: structure storing the options % bounds: structure containing prior bounds % bayestopt_: structure storing information about priors @@ -17,7 +18,7 @@ function [posterior_sampler_options, options_, bayestopt_] = check_posterior_sam % SPECIAL REQUIREMENTS % none -% Copyright © 2015-2022 Dynare Team +% Copyright © 2015-2023 Dynare Team % % This file is part of Dynare. % @@ -391,7 +392,7 @@ if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice') end 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; 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.') end % % % 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; % % % end [V1, D]=eig(invhess); diff --git a/matlab/compute_mh_covariance_matrix.m b/matlab/compute_mh_covariance_matrix.m index c8f96e270..3e1abab15 100644 --- a/matlab/compute_mh_covariance_matrix.m +++ b/matlab/compute_mh_covariance_matrix.m @@ -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 % estimated mode, using the draws from a metropolis-hastings. The estimated posterior mode and covariance matrix are saved in -% a file _mh_mode.mat. +% a file _mh_mode.mat. % % INPUTS -% None. +% o bayestopt_ [struct] characterizing priors +% o fname [string] name of model +% o dname [string] name of directory with metropolis folder % % OUTPUTS % 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 % along with Dynare. If not, see . +MetropolisFolder = CheckPath('metropolis',dname); +BaseName = [MetropolisFolder filesep fname]; -global M_ bayestopt_ - -MetropolisFolder = CheckPath('metropolis',M_.dname); -ModelName = M_.fname; -BaseName = [MetropolisFolder filesep ModelName]; - -record=load_last_mh_history_file(MetropolisFolder, ModelName); +record=load_last_mh_history_file(MetropolisFolder, fname); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; @@ -71,4 +70,4 @@ hh = inv(posterior_covariance); fval = posterior_kernel_at_the_mode; parameter_names = bayestopt_.name; -save([M_.dname filesep 'Output' filesep M_.fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names'); \ No newline at end of file +save([dname filesep 'Output' filesep fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names'); \ No newline at end of file diff --git a/matlab/convergence_diagnostics/McMCDiagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m similarity index 93% rename from matlab/convergence_diagnostics/McMCDiagnostics.m rename to matlab/convergence_diagnostics/mcmc_diagnostics.m index cbe02ce8d..8ad512a68 100644 --- a/matlab/convergence_diagnostics/McMCDiagnostics.m +++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m @@ -1,5 +1,5 @@ -function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) -% function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) +function oo_ = mcmc_diagnostics(options_, estim_params_, M_, oo_) +% function oo_ = mcmc_diagnostics(options_, estim_params_, M_, oo_) % Computes convergence tests % % INPUTS @@ -51,11 +51,11 @@ for b = 1:nblck nfiles = length(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck' num2str(b) '.mat'])); if ~isequal(NumberOfMcFilesPerBlock,nfiles) 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 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 % compute inefficiency factor @@ -111,7 +111,7 @@ PastDraws = sum(record.MhDraws,1); NumberOfDraws = PastDraws(1); 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 end @@ -185,7 +185,7 @@ if nblck == 1 % Brooks and Gelman tests need more than one block if options_.convergence.rafterylewis.indicator 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) - 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]; end Raftery_Lewis_q=options_.convergence.rafterylewis.qrs(1); @@ -218,18 +218,18 @@ xx = Origin:StepSize:NumberOfDraws; NumberOfLines = length(xx); 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 end if TeX && any(strcmp('eps',cellstr(options_.graph_format))) 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,' \n'); 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 % computing are stored in localVars struct. @@ -248,7 +248,7 @@ localVars.M_ = M_; % Like sequential execution! if isnumeric(options_.parallel) - fout = McMCDiagnostics_core(localVars,1,npar,0); + fout = mcmc_diagnostics_core(localVars,1,npar,0); UDIAG = fout.UDIAG; clear fout % Parallel execution! @@ -258,7 +258,7 @@ else end 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; for j=2:totCPU UDIAG = cat(3,UDIAG ,fout(j).UDIAG); @@ -397,7 +397,7 @@ clear UDIAG; % if TeX && any(strcmp('eps',cellstr(options_.graph_format))) 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,' \n'); end diff --git a/matlab/convergence_diagnostics/McMCDiagnostics_core.m b/matlab/convergence_diagnostics/mcmc_diagnostics_core.m similarity index 93% rename from matlab/convergence_diagnostics/McMCDiagnostics_core.m rename to matlab/convergence_diagnostics/mcmc_diagnostics_core.m index 2cb710f64..a9ab4581e 100644 --- a/matlab/convergence_diagnostics/McMCDiagnostics_core.m +++ b/matlab/convergence_diagnostics/mcmc_diagnostics_core.m @@ -1,5 +1,5 @@ -function myoutput = McMCDiagnostics_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) +% function myoutput = mcmc_diagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab) % Computes the Brooks/Gelman (1998) convergence diagnostics, both the % 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 % 5nd column: within sequence kurtosis % 6nd column: sum of within sequence kurtoses; used to compute mean within sequence kurtoses -% Averaging to compute mean moments is done in -% McMCDiagnostics +% Averaging to compute mean moments is done in mcmc_diagnostics % % 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 % Methods for Monitoring Convergence of Iterative Simulations, Journal of Computational % 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. % None. -% Copyright © 2006-2017 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -80,7 +79,7 @@ tmp = zeros(NumberOfDraws*nblck,3); UDIAG = zeros(NumberOfLines,6,npar-fpar+1); 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 waitbarTitle=['Local ']; else diff --git a/matlab/display_estimation_results_table.m b/matlab/display_estimation_results_table.m index 92c8a5348..4d0b7031a 100644 --- a/matlab/display_estimation_results_table.m +++ b/matlab/display_estimation_results_table.m @@ -174,7 +174,7 @@ if ncn end 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 OutputDirectoryName = CheckPath('Output',M_.dname); diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 8bb1337ae..29aa7842d 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname) % SPECIAL REQUIREMENTS % none -% Copyright © 2003-2022 Dynare Team +% Copyright © 2003-2023 Dynare Team % % 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 +dispString = 'Estimation::mcmc'; + if ~exist([M_.dname filesep 'Output'],'dir') if isoctave && octave_ver_less_than('7') && ~exist(M_.dname) % 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 if ~options_.mh_posterior_mode_estimation && options_.cova_compute - try - 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 + check_hessian_at_the_mode(hh, xparam1, M_, estim_params_, options_, bounds); end 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'); end -switch options_.MCMC_jumping_covariance - 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 +invhess = set_mcmc_jumping_covariance(invhess, nx, options_.MCMC_jumping_covariance, bayestopt_, 'dynare_estimation_1'); if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... (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 - outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > 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(['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 + %reset bounds as lb and ub must only be operational during mode-finding + bounds = set_mcmc_prior_bounds(xparam1, bayestopt_, options_, 'dynare_estimation_1'); % Tunes the jumping distribution's scale parameter if options_.mh_tune_jscale.status if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings') - %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, 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') + options_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_, M_, objective_function, xparam1, bounds,... + dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_); 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 warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!') end @@ -479,7 +397,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... if options_.mh_replic || options_.load_mh_file posterior_sampler_options = options_.posterior_sampler_options.current_options; posterior_sampler_options.invhess = invhess; - [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, 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 options_.posterior_sampler_options.current_options = posterior_sampler_options; 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: CutSample(M_, options_, estim_params_); 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'); %reset qz_criterium options_.qz_criterium=qz_criterium_old; @@ -504,7 +422,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... end if ~options_.nodiagnostic 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 if isfield(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 options_.bayesian_irf if error_flag - error('Estimation::mcmc: I cannot compute the posterior IRFs!') + error('%s: I cannot compute the posterior IRFs!',dispString) end PosteriorIRF('posterior'); end if options_.moments_varendo 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 if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results [MetropolisFolder, info] = CheckPath('metropolis',M_.dname); @@ -578,16 +496,16 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... end if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast if error_flag - error('Estimation::mcmc: I cannot compute the posterior statistics!') + error('%s: I cannot compute the posterior statistics!',dispString) end if options_.order==1 && ~options_.particle.status prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts else - error('Estimation::mcmc: Particle Smoothers are not yet implemented.') + error('%s: Particle Smoothers are not yet implemented.',dispString) end end 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 xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); M_ = set_all_parameters(xparam1,estim_params_,M_); diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index a81627bc2..179d39f50 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -80,26 +80,11 @@ if ~isfield(options_,'varobs') error('VAROBS statement is missing!') end +% Checks on VAROBS +check_varobs_are_endo_and_declared_once(options_.varobs,M_.endo_names); % Set the number of observed variables. 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))1 - error(['A variable cannot be declared as observed more than once (' options_.varobs{i} ')!']) - end - end -end - if options_.discretionary_policy if options_.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. 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); - 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 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() + [xparam1, hh] = check_mode_file(xparam1, hh, options_, bayestopt_); end %check for calibrated covariances before updating parameters @@ -627,7 +486,7 @@ if options_.load_results_after_load_mh end 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; end diff --git a/matlab/get_the_name.m b/matlab/get_the_name.m index dfeb7c909..11de078cf 100644 --- a/matlab/get_the_name.m +++ b/matlab/get_the_name.m @@ -32,7 +32,7 @@ function [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, options_) %! @sp 2 %! @strong{This function is called by:} %! @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 %! @strong{This function calls:} %! @sp 1 diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m index 47adc09aa..1340f7e5e 100644 --- a/matlab/initial_estimation_checks.m +++ b/matlab/initial_estimation_checks.m @@ -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.']) 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) -[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 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 +% check and display warning if negative values of stderr or corr params are outside unit circle for Bayesian estimation +if any(BayesInfo.pshape) + check_prior_stderr_corr(EstimatedParameters,BayesInfo); end % display warning if some parameters are still NaN diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m index 48014ba61..2f8c45bd7 100644 --- a/matlab/marginal_density.m +++ b/matlab/marginal_density.m @@ -48,7 +48,7 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); TODROP = floor(options_.mh_drop*TotalNumberOfMhDraws); 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); SIGMA = posterior_covariance; diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m index 690dcac5d..505d1dbeb 100644 --- a/matlab/metropolis_draw.m +++ b/matlab/metropolis_draw.m @@ -72,8 +72,8 @@ if init else if options_.sub_draws>NumberOfDraws*mh_nblck 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('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(['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('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() xparams = 1; % xparams is interpreted as an error flag end diff --git a/matlab/pm3.m b/matlab/pm3.m index c2d19f092..658a58eb7 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -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. -% Copyright © 2007-2018 Dynare Team +% Copyright © 2007-2023 Dynare Team % % 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_ +dispString = 'Estimation::mcmc'; + nn = 3; MaxNumberOfPlotsPerFigure = nn^2; % must be square varlist = names2; @@ -76,7 +78,7 @@ HPD = zeros(2,n2,nvar); if options_.estimation.moments_posterior_density.indicator Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar); end -fprintf(['Estimation::mcmc: ' tit1 '\n']); +fprintf(['%s: ' tit1 '\n'],dispString); k = 0; filter_step_ahead_indicator=0; filter_covar_indicator=0; @@ -161,7 +163,7 @@ elseif filter_covar_indicator oo_.FilterCovariance.post_deciles=post_deciles; oo_.FilterCovariance.HPDinf=squeeze(hpd_interval(:,:,:,1)); oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2)); - fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); + fprintf(['%s: ' tit1 ', done!\n'],dispString); return elseif state_uncert_indicator draw_dimension=4; @@ -183,7 +185,7 @@ elseif state_uncert_indicator oo_.Smoother.State_uncertainty.post_deciles=post_deciles; oo_.Smoother.State_uncertainty.HPDinf=squeeze(hpd_interval(:,:,:,1)); oo_.Smoother.State_uncertainty.HPDsup=squeeze(hpd_interval(:,:,:,2)); - fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); + fprintf(['%s: ' tit1 ', done!\n'],dispString); return end @@ -280,7 +282,7 @@ else end 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 end %% @@ -378,4 +380,4 @@ if ~options_.nograph && ~options_.no_graph.posterior end end -fprintf(['Estimation::mcmc: ' tit1 ', done!\n']); +fprintf(['%s: ' tit1 ', done!\n'],dispString); diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m index fbc395895..8a98bfc24 100644 --- a/matlab/posterior_sampler.m +++ b/matlab/posterior_sampler.m @@ -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 % along with Dynare. If not, see . +dispString = 'Estimation::mcmc'; + vv = sampler_options.invhess; % Initialization of the sampler [ 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 skipline() -disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) -disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) -disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) -disp(['Estimation::mcmc: Current acceptance ratio per chain: ']) +fprintf('%s: Number of mh files: %u per block.\n', dispString, NewFile(1)); +fprintf('%s: Total number of generated files: %u.\n', dispString, NewFile(1)*nblck); +fprintf('%s: Total number of iterations: %u.\n', dispString, (NewFile(1)-1)*MAX_nruns+irun-1); +fprintf('%s: Current acceptance ratio per chain:\n', dispString); for i=1:nblck if i<10 disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) @@ -185,7 +187,7 @@ for i=1:nblck end end 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 if i<10 disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index 2ac7f6963..36c308777 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -55,6 +55,8 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +dispString = 'Estimation::mcmc'; + %Initialize outputs ix2 = []; ilogpo2 = []; @@ -79,22 +81,22 @@ d = chol(vv); if ~options_.load_mh_file && ~options_.mh_recover % Here we start a new Metropolis-Hastings, previous draws are discarded. if NumberOfBlocks > 1 - disp('Estimation::mcmc: Multiple chains mode.') + fprintf('%s: Multiple chains mode.\n',dispString); else - disp('Estimation::mcmc: One Chain mode.') + fprintf('%s: One Chain mode.\n',dispString); end % Delete old mh files if any... files = dir([BaseName '_mh*_blck*.mat']); if length(files) delete([BaseName '_mh*_blck*.mat']); - disp('Estimation::mcmc: Old mh-files successfully erased!') + fprintf('%s: Old mh-files successfully erased!\n',dispString); end % Delete old Metropolis log file. file = dir([ MetropolisFolder '/metropolis.log']); if length(file) delete([ MetropolisFolder '/metropolis.log']); - disp('Estimation::mcmc: Old metropolis.log file successfully erased!') - disp('Estimation::mcmc: Creation of a new metropolis.log file.') + fprintf('%s: Old metropolis.log file successfully erased!\n',dispString) + fprintf('%s: Creation of a new metropolis.log file.\n',dispString) end fidlog = fopen([MetropolisFolder '/metropolis.log'],'w'); 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 RecordFile0=strrep(RecordFile0,'\',filesep); if isempty(dir(RecordFile0)) - disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_record option') - error('Estimation::mcmc: path to record file is not found') + fprintf('%s: Wrong value for mh_initialize_from_previous_mcmc_record option.\n',dispString); + error('%s: Path to record file is not found!',dispString) else record0=load(RecordFile0); end @@ -133,31 +135,31 @@ if ~options_.load_mh_file && ~options_.mh_recover record0=load_last_mh_history_file(MetropolisFolder0, ModelName0); end 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 % mh_files = dir([ MetropolisFolder0 filesep ModelName0 '_mh*.mat']); % 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 - disp('Estimation::mcmc: Initializing from past Metropolis-Hastings simulations...') - disp(['Estimation::mcmc: Past MH path ' MetropolisFolder0 ]) - disp(['Estimation::mcmc: Past model name ' ModelName0 ]) + fprintf('%s: Initializing from past Metropolis-Hastings simulations...\n',dispString); + fprintf('%s: Past MH path %s\n',dispString,MetropolisFolder0); + fprintf('%s: Past model name %s\n', dispString, ModelName0); fprintf(fidlog,' Loading initial values from previous MH\n'); fprintf(fidlog,' Past MH path: %s\n', MetropolisFolder0 ); fprintf(fidlog,' Past model name: %s\n', ModelName0); fprintf(fidlog,' \n'); past_number_of_blocks = record0.Nblck; if past_number_of_blocks ~= NumberOfBlocks - disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') - disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) - disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) + fprintf('%s: The specified number of blocks doesn''t match with the previous number of blocks!\n', dispString); + fprintf('%s: You declared %u blocks, but the previous number of blocks was %u.\n', dispString, NumberOfBlocks, past_number_of_blocks); + fprintf('%s: I will run the Metropolis-Hastings with %u block.\n', dispString, past_number_of_blocks); NumberOfBlocks = past_number_of_blocks; options_.mh_nblck = NumberOfBlocks; end if ~isempty(PriorFile0) if isempty(dir(PriorFile0)) - disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_prior option') - error('Estimation::mcmc: path to prior file is not found') + fprintf('%s: Wrong value for mh_initialize_from_previous_mcmc_prior option.', dispString); + error('%s: Path to prior file is not found!',dispString); else bayestopt0 = load(PriorFile0); 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 set_dynare_seed('default'); 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 ix2 = zeros(NumberOfBlocks,npar); ilogpo2 = zeros(NumberOfBlocks,1); @@ -217,56 +219,54 @@ if ~options_.load_mh_file && ~options_.mh_recover end init_iter = init_iter + 1; 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 - 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') 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 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 - 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 if isfield(options_,'mh_init_scale') - disp(['Estimation::mcmc: You should reduce mh_init_scale...']) - fprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.\n',options_.mh_init_scale) - options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... '); + fprintf('%s: You should reduce mh_init_scale...\n',dispString); + fprintf('%s: Parameter ''mh_init_scale'' is equal to %f.\n',dispString,options_.mh_init_scale); + options_.mh_init_scale_factor = input(sprintf('%s: Enter a new value... ',dispString)); else - disp(['Estimation::mcmc: You should reduce mh_init_scale_factor...']) - fprintf('Estimation::mcmc: Parameter mh_init_scale_factor is equal to %f.\n',options_.mh_init_scale_factor) - options_.mh_init_scale_factor = input('Estimation::mcmc: Enter a new value... '); + fprintf('%s: You should reduce ''mh_init_scale_factor''...\n',dispString); + fprintf('%s: Parameter ''mh_init_scale_factor'' is equal to %f.\n',dispString,options_.mh_init_scale_factor); + options_.mh_init_scale_factor = input(sprintf('%s: Enter a new value... ',dispString)); end end trial = trial+1; end end 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); return end end fprintf(fidlog,' \n'); - disp('Estimation::mcmc: Initial values found!') - skipline() + fprintf('%s: Initial values found!\n\n',dispString); else% Case 2: one chain (we start from the posterior mode) fprintf(fidlog,[' Initial values of the parameters:\n']); candidate = transpose(xparam1(:));% if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub) ix2 = candidate; ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); - disp('Estimation::mcmc: Initialization at the posterior mode.') - skipline() + fprintf('%s: Initialization at the posterior mode.\n\n',dispString); fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']); for i=1:length(ix2(1,:)) fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']); end fprintf(fidlog,[' Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']); else - disp('Estimation::mcmc: Initialization failed...') - disp('Estimation::mcmc: The posterior mode lies outside the prior bounds.') + fprintf('%s: Initialization failed...\n',dispString); + fprintf('%s: The posterior mode lies outside the prior bounds.\n',dispString); fclose(fidlog); return end @@ -279,7 +279,7 @@ if ~options_.load_mh_file && ~options_.mh_recover % Delete the mh-history files delete_mh_history_files(MetropolisFolder, ModelName); % 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); AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns; record.Sampler = options_.posterior_sampler_options.posterior_sampling_method; @@ -312,8 +312,7 @@ if ~options_.load_mh_file && ~options_.mh_recover record.ProposalCovariance=d; fprintf('Ok!\n'); id = write_mh_history_file(MetropolisFolder, ModelName, record); - disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat']) - skipline() + fprintf('%s: Details about the MCMC are available in %s_mh_history_%u.mat\n\n', dispString,BaseName,id); 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 lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']); @@ -332,15 +331,15 @@ if ~options_.load_mh_file && ~options_.mh_recover fclose(fidlog); elseif options_.load_mh_file && ~options_.mh_recover % 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); 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 record.MCMCConcludedSuccessfully=0; %reset indicator for this run mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']); 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 fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fprintf(fidlog,'\n'); @@ -351,9 +350,9 @@ elseif options_.load_mh_file && ~options_.mh_recover fprintf(fidlog,' \n'); past_number_of_blocks = record.Nblck; if past_number_of_blocks ~= NumberOfBlocks - disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') - disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) - disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) + fprintf('%s: The specified number of blocks doesn''t match with the previous number of blocks!\n',dispString); + fprintf('%s: You declared %u blocks, but the previous number of blocks was %u.\n', dispString,NumberOfBlocks,past_number_of_blocks); + fprintf('%s: I will run the Metropolis-Hastings with %u blocks.\n', dispString,past_number_of_blocks); NumberOfBlocks = past_number_of_blocks; options_.mh_nblck = NumberOfBlocks; end @@ -369,10 +368,10 @@ elseif options_.load_mh_file && ~options_.mh_recover end ilogpo2 = record.LastLogPost; 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; 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)]; NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber; NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile; @@ -387,28 +386,27 @@ elseif options_.load_mh_file && ~options_.mh_recover record.InitialSeeds = record.LastSeeds; write_mh_history_file(MetropolisFolder, ModelName, record); fprintf('Done.\n') - disp(['Estimation::mcmc: Ok. I have loaded ' int2str(NumberOfPreviousSimulations) ' simulations.']) - skipline() + fprintf('%s: Ok. I have loaded %u simulations.\n\n', dispString,NumberOfPreviousSimulations); fclose(fidlog); elseif options_.mh_recover % 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); NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains. options_.mh_nblck = NumberOfBlocks; %% check consistency of options 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('Estimation::mcmc: I am resetting mh_replic to %u\n',record.MhDraws(end,1)) - options_.mh_replic=record.MhDraws(end,1); + fprintf('\n%s: You cannot specify a different mh_replic than in the chain you are trying to recover\n',dispString); + fprintf('%s: I am resetting mh_replic to %u\n',dispString,record.MhDraws(end,1)); + options_.mh_replic = record.MhDraws(end,1); nruns = ones(NumberOfBlocks,1)*options_.mh_replic; end if ~isnan(record.MAX_nruns(end,1)) %field exists 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('Estimation::mcmc: I am resetting MAX_nruns to %u\n',record.MAX_nruns(end,1)) + fprintf('\n%s: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n',dispString); + fprintf('%s: I am resetting MAX_nruns to %u\n',dispString,record.MAX_nruns(end,1)); MAX_nruns=record.MAX_nruns(end,1); end end @@ -451,7 +449,7 @@ elseif options_.mh_recover LastFileFullIndicator=1; end 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 %% Now find out what exactly needs to be redone % 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 if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) if isnumeric(options_.parallel) - disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!') - disp(' You have to edit the mod file and remove the mh_recover option') - disp(' in the estimation command') - error('Estimation::mcmc: mh_recover option not required!') + fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString); + fprintf(' You have to edit the mod file and remove the mh_recover option\n'); + fprintf(' in the estimation command.\n'); + error('%s: mh_recover option not required!',dispString) end end % 2. Something needs to be done; find out what @@ -482,10 +480,10 @@ elseif options_.mh_recover FBlock = zeros(NumberOfBlocks,1); while FirstBlock <= NumberOfBlocks 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; else - disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!']) + fprintf('%s: Chain %u is complete!\n', dispString,FirstBlock); end FirstBlock = FirstBlock+1; end @@ -537,8 +535,8 @@ elseif options_.mh_recover record.InitialSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Unifor; record.InitialSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Normal; else - fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n') - fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n') + fprintf('%s: You are trying to recover a chain generated with an older Dynare version.\n',dispString); + fprintf('%s: I am using the default seeds to continue the chain.\n',dispString); end if update_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).Normal=loaded_results.LastSeeds.(['file' int2str(ExpectedNumberOfMhFilesPerBlock)]).Normal; else - fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n') - fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n') + fprintf('%s: You are trying to recover a chain generated with an older Dynare version.\n',dispString); + fprintf('%s: I am using the default seeds to continue the chain.\n',dispString); end if isfield(loaded_results,'accepted_draws_this_chain') 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); 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 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; bayestopt_.jscale=record.ProposalScaleVec; else if ~isequal(options_.posterior_sampler_options.posterior_sampling_method,'slice') % pass through input d unaltered 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) - 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 - 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 else % 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 if isfield(record,'Sampler') 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 end end diff --git a/matlab/set_prior.m b/matlab/set_prior.m index bc4fff302..635568dc8 100644 --- a/matlab/set_prior.m +++ b/matlab/set_prior.m @@ -5,7 +5,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params % INPUTS % o estim_params_ [structure] characterizing parameters to be estimated. % o M_ [structure] characterizing the model. -% o options_ [structure] +% o options_ [structure] characterizing the options. % % OUTPUTS % 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 % None -% Copyright © 2003-2018 Dynare Team +% Copyright © 2003-2023 Dynare Team % % This file is part of Dynare. % diff --git a/matlab/utilities/estimation/check_hessian_at_the_mode.m b/matlab/utilities/estimation/check_hessian_at_the_mode.m new file mode 100644 index 000000000..8fa880253 --- /dev/null +++ b/matlab/utilities/estimation/check_hessian_at_the_mode.m @@ -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 . +% ========================================================================= + +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). +% ========================================================================= + +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 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() \ No newline at end of file diff --git a/matlab/utilities/estimation/check_prior_stderr_corr.m b/matlab/utilities/estimation/check_prior_stderr_corr.m new file mode 100644 index 000000000..0cdc7b71c --- /dev/null +++ b/matlab/utilities/estimation/check_prior_stderr_corr.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/utilities/estimation/check_steady_state_changes_parameters.m b/matlab/utilities/estimation/check_steady_state_changes_parameters.m new file mode 100644 index 000000000..1cd40b293 --- /dev/null +++ b/matlab/utilities/estimation/check_steady_state_changes_parameters.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m b/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m new file mode 100644 index 000000000..75daffc91 --- /dev/null +++ b/matlab/utilities/estimation/check_varobs_are_endo_and_declared_once.m @@ -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 . +% ========================================================================= + +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)) 1 + error(['VAROBS: a variable cannot be declared as observed more than once (' varobs{i} ')!']) + end + end +end \ No newline at end of file diff --git a/matlab/utilities/estimation/set_mcmc_jumping_covariance.m b/matlab/utilities/estimation/set_mcmc_jumping_covariance.m new file mode 100644 index 000000000..b0cb99d49 --- /dev/null +++ b/matlab/utilities/estimation/set_mcmc_jumping_covariance.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/utilities/estimation/set_mcmc_prior_bounds.m b/matlab/utilities/estimation/set_mcmc_prior_bounds.m new file mode 100644 index 000000000..76bbf09ca --- /dev/null +++ b/matlab/utilities/estimation/set_mcmc_prior_bounds.m @@ -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 . +% ========================================================================= + +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 \ No newline at end of file diff --git a/matlab/utilities/estimation/tune_mcmc_mh_jscale_wrapper.m b/matlab/utilities/estimation/tune_mcmc_mh_jscale_wrapper.m new file mode 100644 index 000000000..d0be1e6e8 --- /dev/null +++ b/matlab/utilities/estimation/tune_mcmc_mh_jscale_wrapper.m @@ -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 . +% ========================================================================= + +% 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{:}); \ No newline at end of file diff --git a/tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc b/tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc index 6ba8a8f8a..dc2eb4971 100644 --- a/tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc +++ b/tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc @@ -193,9 +193,10 @@ INT(0)^2*INFL(1)^4; %redundant @#endif end; +M_.matched_moments_orig = M_.matched_moments; method_of_moments( -% Necessery options +% Necessary options mom_method = @{MoM_Method} % method of moments method; possible values: GMM|SMM , datafile = 'AnScho_MoM_data_2.mat' % name of filename with data diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod index 41dcdd8c9..e8af3f66a 100644 --- a/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod +++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod @@ -152,7 +152,7 @@ end method_of_moments( - % Necessery options + % Necessary options mom_method = GMM % method of moments method; possible values: GMM|SMM , datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod index 60920a15c..d8e6077de 100644 --- a/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod +++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod @@ -139,7 +139,7 @@ end @#for mommethod in ["SMM"] method_of_moments( - % Necessery options + % Necessary options mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data diff --git a/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod b/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod index 00b89a5a2..6bc36f3d4 100644 --- a/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod +++ b/tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod @@ -110,7 +110,7 @@ save('test_matrix.mat','weighting_matrix') @#for mommethod in ["GMM", "SMM"] method_of_moments( - % Necessery options + % Necessary options mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data