diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index edcc52d1b5eac606a56b394721be281a414de16f..e762e63c487e5d0b65eef9eeb27156b345ebf7aa 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -6,31 +6,47 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% 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 matched_moments structure
+% o Checks and transformations for matched_irfs and matched_irfs_weights 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 Mode computation: optimization
+% - GMM/SMM: iterated optimization
+% - IRF_MATCHING: optimization
+% o Bayesian MCMC estimation
% o Display of results
+% - GMM/SMM: J-Test and fit of moments
+% - IRF_MATCHING: fit of irfs
% o Clean up
% -------------------------------------------------------------------------
+% Note that we call a "mode" the minimum of the objective function, i.e.
+% the parameter vector that minimizes the distance between the moments/irfs
+% computed from the model and the moments/irfs computed from the data.
+% -------------------------------------------------------------------------
% 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.
-% =========================================================================
+% o Ruge-Murcia (2007): "Methods to Estimate Dynamic Stochastic General Equilibrium Models", Journal of Economic Dynamics and Control, 31(8):2599-2636.
+% IRF MATCHING:
+% o Christiano, Trabandt, Walentin (2010): "DSGE Models for Monetary Policy Analysis." In Handbook of Monetary Economics, 3:285–367.
+% o Christiano, Eichenbaum, Trabandt (2016): "Unemployment and Business Cycles." Econometrica, 84: 1523-1569.
+% o Ruge-Murcia (2020): "Estimating Nonlinear Dynamic Equilibrium Models by Matching Impulse Responses", Economics Letters, 197.
+% -------------------------------------------------------------------------
% INPUTS
% o bayestopt_: [structure] information about priors
% o options_: [structure] information about global options
-% o oo_: [structure] storage for results
+% o oo_: [structure] 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 matched_irfs: [cell] information about selected irfs to match in IRF_MATCHING estimation
+% o matched_irfs_weights: [cell] information about entries in weight matrix for an IRF_MATCHING estimation
% o options_mom_: [structure] information about settings specified by the user
% -------------------------------------------------------------------------
% OUTPUTS
@@ -42,25 +58,26 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o driver.m
% -------------------------------------------------------------------------
% This function calls
-% o cellofchararraymaxlength
% o check_for_calibrated_covariances
+% o check_mode_file
+% o check_posterior_sampler_options
% o check_prior_bounds
% o check_prior_stderr_corr
% o check_steady_state_changes_parameters
% o check_varobs_are_endo_and_declared_once
+% o check_hessian_at_the_mode
% 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 mode_check
+% o mom.check_irf_matching_file
% o mom.default_option_mom_values
% o mom.get_data_moments
+% o mom.matched_irfs_blocks
% o mom.matched_moments_block
% o mom.objective_function
% o mom.optimal_weighting_matrix
@@ -77,8 +94,14 @@ function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_,
% o set_state_space
% o skipline
% o test_for_deep_parameters_calibration
+% o transform_prior_to_laplace_prior
% o warning_config
-% =========================================================================
+% -------------------------------------------------------------------------
+% Maintaining Author(s):
+% o Willi Mutschler (willi@mutschler.eu)
+% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
+% -------------------------------------------------------------------------
+
% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
@@ -95,32 +118,47 @@ 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 <https://www.gnu.org/licenses/>.
-% -------------------------------------------------------------------------
-% Maintaining Author(s):
-% o Willi Mutschler (willi@mutschler.eu)
-% o Johannes Pfeifer (johannes.pfeifer@unibw.de)
-% =========================================================================
+
% -------------------------------------------------------------------------
-% TO DO LISTS
+% TO DO LISTS AND IDEAS
% -------------------------------------------------------------------------
% GENERAL
+% - PREPROCESSOR CHANGE CALL TO FUNCTION (SEE ABOVE)
% - 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
+% - do endogenous_prior_restrictions work?, create test for prior restrictions file
+% - implement penalty objective function for optimization
+% - test optimizers (what about analytic Jacobians? Also check if mode compute is a string and additional optimizers also a string)
+% - factorize mode_compute codes
+% - analytic_jacobian with mode_compute 5
+% - decide on default mode_compute
+% - use same names for variables (e.g. for weighting matrix)
+% - mom.objective_function: check the info values and make use of meaningful penalties (which numbers do we use??)
+% - add Approximate Bayesian Computation (ABC) option
+% - merge mode_compute functions for GMM/SMM and IRF_MATCHING
% GMM/SMM
+% - do true Bayesian MCMC sampling and not just penalized
% - 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
+% - enable to use first moments even if prefilter option is set
+% IRF_MATCHING/SMM
+% - add option to do simulations with extended path
+% - add option to do simulations with perfect_foresight and perfect_foresight_with_expectation_errors
+% - add option to do simulations with occbin
+% - factorize stoch_simul irf codes that are the same as in mom.objective_function
+% IRF_MATCHING
+% - add analytic_standard_errors and analytic_jacobian (at least for order=1)
+% - What about iterating over mode_compute and updating the weighting matrix? Is that also useful for IRF_MATCHING?
+% Can we also use optimal_weighting_matrix for IRF_MATCHING or is this done outside of our codes?
+% - check all mode compute options (and also optim options)
+% - use_penalized_objective_for_hessian
+% - do we need bayesian_irf? If not remove as option
+% - check order > 1
+% - print more info to console
+
fprintf('\n==== Method of Moments Estimation (%s) ====\n\n',options_mom_.mom.mom_method);
@@ -139,6 +177,10 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
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
+elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
+ if ~isfield(M_,'matched_irfs') || isempty(M_.matched_irfs) % structure storing the irfs used for matching
+ error('method_of_moments: You need to provide a ''matched_irfs'' 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!');
@@ -231,28 +273,50 @@ end
% -------------------------------------------------------------------------
-% estimated parameters: checks and transformations on values, priors, bounds
+% matched_irfs: checks and transformations
+% -------------------------------------------------------------------------
+if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
+ [oo_.mom.data_moments, oo_.mom.weighting_info.W, options_mom_.mom.irfIndex, options_mom_.irf] = mom.matched_irfs_blocks(M_.matched_irfs, M_.matched_irfs_weights, options_mom_.varobs_id, options_mom_.obs_nbr, M_.exo_nbr, M_.endo_names);
+ oo_.mom.weighting_info.Winv = inv(oo_.mom.weighting_info.W);
+ oo_.mom.weighting_info.Winv_logdet = 2*sum(log(diag(chol(oo_.mom.weighting_info.Winv)))); % use this robust option to avoid inf/nan
+ options_mom_.mom.mom_nbr = length(options_mom_.mom.irfIndex);
+end
+
+% -------------------------------------------------------------------------
+% irf_matching_file: checks and transformations
+% -------------------------------------------------------------------------
+if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
+ [options_mom_.mom.irf_matching_file.name, options_mom_.mom.irf_matching_file.path] = mom.check_irf_matching_file(options_mom_.mom.irf_matching_file.name);
+ % check for irf_matching_file
+ if ~( isempty(options_mom_.mom.irf_matching_file.path) || strcmp(options_mom_.mom.irf_matching_file.path,'.') )
+ fprintf('\nAdding %s to MATLAB''s path.\n',options_mom_.mom.irf_matching_file.path);
+ addpath(options_mom_.mom.irf_matching_file.path);
+ end
+end
+
+
+% -------------------------------------------------------------------------
+% estimated parameters: checks and transformations on values, priors, bounds, posterior options
% -------------------------------------------------------------------------
% 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 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);
+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);
% 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
-
+% check that the provided mode_file is compatible with the current estimation settings
+if ~isempty(options_mom_.mode_file) && ( ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation) )
+ [xparam0, hessian_xparam0] = check_mode_file(xparam0, hessian_xparam0, options_mom_, bayestopt_);
+end
% 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_;
@@ -265,7 +329,6 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
bayestopt_ = mom.transform_prior_to_laplace_prior(bayestopt_);
end
end
-
% check for calibrated covariances before updating parameters
estim_params_ = check_for_calibrated_covariances(estim_params_,M_);
@@ -283,12 +346,10 @@ if options_mom_.use_calibration_initialization % set calibration as starting val
[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_) && ~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_) && doBayesianEstimation
% plot prior densities
@@ -297,6 +358,8 @@ if ~isempty(bayestopt_) && doBayesianEstimation
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
+ else
+ plot_priors(bayestopt_,M_,estim_params_,options_mom_,'Priors');
end
end
% set prior bounds
@@ -304,19 +367,22 @@ if ~isempty(bayestopt_) && doBayesianEstimation
BoundsInfo.lb = max(BoundsInfo.lb,lb);
BoundsInfo.ub = min(BoundsInfo.ub,ub);
else
- % no priors are declared so Dynare will estimate the parameters with
- % Frequentist methods using inequality constraints for the parameters
+ % no priors are declared so Dynare will estimate the parameters with Frequentist methods using inequality constraints for the parameters
BoundsInfo.lb = lb;
BoundsInfo.ub = ub;
- if options_mom_.mom.penalized_estimator
+ if (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator
fprintf('Penalized estimation turned off as you did not declare priors\n');
options_mom_.mom.penalized_estimator = 0;
+ else
+ if isfield(options_mom_,'mh_replic') && options_mom_.mh_replic > 0
+ fprintf('Setting ''mh_replic=0'' as you did not declare priors.\n');
+ options_mom_.mh_replic = 0;
+ end
end
end
% set correct bounds for standard deviations and correlations
BoundsInfo = mom.set_correct_bounds_for_stderr_corr(estim_params_,BoundsInfo);
-
% test if initial values of the estimated parameters are all between the prior lower and upper bounds
if options_mom_.use_calibration_initialization
try
@@ -328,35 +394,47 @@ if options_mom_.use_calibration_initialization
else
check_prior_bounds(xparam0,BoundsInfo,M_,estim_params_,options_mom_,bayestopt_);
end
-
% 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)
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 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;
-
+% storing prior parameters in results structure
+if doBayesianEstimation || ( (strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) && options_mom_.mom.penalized_estimator)
+ 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;
+end
% set all parameters
M_ = set_all_parameters(xparam0,estim_params_,M_);
-
% provide warning if there is NaN in parameters
test_for_deep_parameters_calibration(M_);
-
+% set jscale
+if doBayesianEstimationMCMC
+ if ~strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method,'slice')
+ if isempty(options_mom_.mh_jscale)
+ options_mom_.mh_jscale = 2.38/sqrt(number_of_estimated_parameters); % use optimal value for univariate normal distribution (check_posterior_sampler_options and mode_compute=6 may overwrite this setting)
+ end
+ bayestopt_.jscale(find(isnan(bayestopt_.jscale))) = options_mom_.mh_jscale;
+ end
+end
+% initialization of posterior sampler options
+if doBayesianEstimationMCMC
+ [current_options, options_mom_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_mom_, BoundsInfo, bayestopt_);
+ options_mom_.posterior_sampler_options.current_options = current_options;
+ if strcmp(current_options.posterior_sampling_method,'slice') && current_options.use_mh_covariance_matrix && ~current_options.rotated
+ error('method_of_moments: Using the slice sampler with the ''use_mh_covariance_matrix'' option requires also setting the ''rotated'' option!');
+ end
+end
+% warning if prior allows that stderr parameters are negative or corr parameters are outside the unit circle
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 value of prior density
- [~,~,~,info]= priordens(xparam0,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
+ [~,~,~,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!');
@@ -367,9 +445,9 @@ end
% -------------------------------------------------------------------------
% datafile: checks and transformations
% -------------------------------------------------------------------------
-% Build dataset
+% build dataset
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
- % Check if datafile has same name as mod file
+ % 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''!');
@@ -379,13 +457,13 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
if ~isempty(dataset_)
options_mom_.nobs = dataset_.nobs;
end
- % Check length of data for estimation of second moments
+ % 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
+ % 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
+ % get data moments for the method of moments
[oo_.mom.data_moments, oo_.mom.m_data] = mom.get_data_moments(dataset_.data, options_mom_.mom.obs_var, oo_.dr.inv_order_var, M_.matched_moments, options_mom_);
if ~isreal(dataset_.data)
error('method_of_moments: The data moments contain complex values!');
@@ -394,7 +472,7 @@ end
% -------------------------------------------------------------------------
-% SMM: Get shock series fand set variance correction factor
+% SMM: Get shock series and 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);
@@ -428,7 +506,6 @@ end
% -------------------------------------------------------------------------
% checks for steady state at initial parameters
% -------------------------------------------------------------------------
-
% check if steady state solves static model and if steady-state changes estimated parameters
if options_mom_.steadystate.nocheck
steadystate_check_flag_vec = [0 1];
@@ -445,7 +522,6 @@ if steady_state_changes_parameters && strcmp(options_mom_.mom.mom_method,'GMM')
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_);
@@ -455,17 +531,25 @@ test_for_deep_parameters_calibration(M_);
% -------------------------------------------------------------------------
objective_function = str2func('mom.objective_function');
try
- % Check for NaN or complex values of moment-distance-funtion evaluated at initial parameters
+ % 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')
- weighting_info.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
+ oo_.mom.weighting_info.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
end
tic_id = tic;
- [fval, info] = feval(objective_function, xparam0, oo_.mom.data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ [fval, info] = feval(objective_function, xparam0, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
elapsed_time = toc(tic_id);
if isnan(fval)
- error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!');
+ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
+ error('method_of_moments: The initial value of the objective function with identity weighting matrix is NaN!');
+ else
+ error('method_of_moments: The initial value of the objective function is NaN!');
+ end
elseif imag(fval)
- error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!');
+ if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
+ error('method_of_moments: The initial value of the objective function with identity weighting matrix is complex!');
+ else
+ error('method_of_moments: The initial value of the objective function is complex!');
+ end
end
if info(1) > 0
disp('method_of_moments: Error in computing the objective function for initial parameter values')
@@ -494,15 +578,33 @@ end
% -------------------------------------------------------------------------
% print some info to console
% -------------------------------------------------------------------------
-mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters);
+mom.print_info_on_estimation_settings(options_mom_, number_of_estimated_parameters, doBayesianEstimation);
% -------------------------------------------------------------------------
-% GMM/SMM: iterated estimation
+% compute mode for GMM/SMM
% -------------------------------------------------------------------------
if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
- % compute mode
[xparam1, oo_.mom.weighting_info, oo_.mom.verbose] = mom.mode_compute_gmm_smm(xparam0, objective_function, oo_.mom.m_data, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+end
+
+% -------------------------------------------------------------------------
+% compute mode for IRF matching
+% -------------------------------------------------------------------------
+if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
+ if ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation)
+ [xparam1, hessian_xparam1, fval, oo_.mom.verbose] = mom.mode_compute_irf_matching(xparam0, hessian_xparam0, objective_function, doBayesianEstimation, oo_.mom.weighting_info, oo_.mom.data_moments, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ else
+ xparam1 = xparam0;
+ hessian_xparam1 = hessian_xparam0;
+ end
+end
+
+
+% -------------------------------------------------------------------------
+% compute standard errors and initialize covariance of the proposal distribution
+% -------------------------------------------------------------------------
+if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
% 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 oo_ (in particular to get oo_.mom.model_moments)
@@ -511,31 +613,229 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
end
[~, ~, ~, ~, ~, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs] = feval(objective_function, xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
- [stdh, hessian_xparam1] = mom.standard_errors(xparam1, objective_function, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.m_data, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
- hessian_xparam1 = inv(hessian_xparam1); % mom.standard_errors returns the asymptotic covariance matrix
+ [stdh, invhess] = mom.standard_errors(xparam1, objective_function, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.m_data, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ if options_mom_.cova_compute
+ hessian_xparam1 = inv(invhess);
+ end
+else
+ if ~doBayesianEstimation || ~options_mom_.mh_posterior_mode_estimation
+ if doBayesianEstimation
+ oo_.mom.posterior.optimization.mode = xparam1;
+ if exist('fval','var')
+ oo_.mom.posterior.optimization.log_density = -fval;
+ end
+ end
+ if options_mom_.cova_compute
+ hsd = sqrt(diag(hessian_xparam1)); % represent the curvature (or second derivatives) of the likelihood with respect to each parameter being estimated.
+ invhess = inv(hessian_xparam1./(hsd*hsd'))./(hsd*hsd'); % before taking the inverse scale the Hessian matrix by dividing each of its elements by the outer product of hsd such that the diagonal of the resulting matrix is approximately 1. This kind of scaling can help in regularizing the matrix and potentially improves its condition number, which in turn can make the matrix inversion more stable.
+ stdh = sqrt(diag(invhess));
+ if doBayesianEstimation
+ oo_.mom.posterior.optimization.Variance = invhess;
+ end
+ end
+ else
+ variances = bayestopt_.p2.*bayestopt_.p2;
+ idInf = isinf(variances);
+ variances(idInf) = 1;
+ invhess = options_mom_.mh_posterior_mode_estimation*diag(variances);
+ xparam1 = bayestopt_.p5;
+ idNaN = isnan(xparam1);
+ xparam1(idNaN) = bayestopt_.p1(idNaN);
+ outside_bound_pars=find(xparam1 < BoundsInfo.lb | xparam1 > BoundsInfo.ub);
+ xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars);
+ end
+ if ~options_mom_.cova_compute
+ stdh = NaN(length(xparam1),1);
+ end
+end
+
+
+% -------------------------------------------------------------------------
+% display estimation results at mode
+% -------------------------------------------------------------------------
+if doBayesianEstimation && ~options_mom_.mom.penalized_estimator && ~options_mom_.mh_posterior_mode_estimation
+ % display table with Bayesian mode estimation results and store parameter estimates and standard errors in oo_
+ oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, 'Posterior', 'posterior');
+ % Laplace approximation to the marginal log density
+ if options_mom_.cova_compute
+ estim_params_nbr = size(xparam1,1);
+ if ispd(invhess)
+ log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
+ likelihood = feval(objective_function, xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ oo_.mom.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
+ else
+ oo_.mom.MarginalDensity.LaplaceApproximation = NaN;
+ end
+ fprintf('\nLog data density [Laplace approximation] is %f.\n',oo_.mom.MarginalDensity.LaplaceApproximation);
+ end
+elseif ~doBayesianEstimation || (doBayesianEstimation && options_mom_.mom.penalized_estimator)
+ % display table with Frequentist estimation results and store parameter estimates and standard errors in oo_
+ oo_.mom = display_estimation_results_table(xparam1, stdh, M_, options_mom_, estim_params_, bayestopt_, oo_.mom, prior_dist_names, options_mom_.mom.mom_method, lower(options_mom_.mom.mom_method));
end
% -------------------------------------------------------------------------
% checks for mode and hessian at the mode
% -------------------------------------------------------------------------
+if (~doBayesianEstimation && options_mom_.cova_compute) || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation && options_mom_.cova_compute)
+ check_hessian_at_the_mode(hessian_xparam1, xparam1, M_, estim_params_, options_, BoundsInfo);
+end
if options_mom_.mode_check.status
- mode_check(objective_function, xparam1, hessian_xparam1, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, true,...
- oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ if ~doBayesianEstimation || (doBayesianEstimation && ~options_mom_.mh_posterior_mode_estimation)
+ mode_check(objective_function, xparam1, diag(stdh), options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, true,... % use diag(stdh) instead of hessian_xparam1 as mode_check uses diagonal elements
+ oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ end
+end
+
+% -------------------------------------------------------------------------
+% Bayesian MCMC estimation
+% -------------------------------------------------------------------------
+if doBayesianEstimationMCMC
+ invhess = set_mcmc_jumping_covariance(invhess, length(xparam1), options_mom_.MCMC_jumping_covariance, bayestopt_, 'method_of_moments');
+ % reset bounds as lb and ub must only be operational during mode-finding
+ BoundsInfo = set_mcmc_prior_bounds(xparam1, bayestopt_, options_mom_, 'method_of_moments');
+ % tunes the jumping distribution's scale parameter
+ if isfield(options_mom_,'mh_tune_jscale') && options_mom_.mh_tune_jscale.status
+ if strcmp(options_mom_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
+ options_mom_.mh_jscale = tune_mcmc_mh_jscale_wrapper(invhess, options_mom_, M_, objective_function, xparam1, BoundsInfo,...
+ oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+ bayestopt_.jscale(:) = options_mom_.mh_jscale;
+ fprintf('mh_tune_jscale: mh_jscale has been set equal to %s.\n', num2str(options_mom_.mh_jscale));
+ else
+ warning('mh_tune_jscale is only available with ''random_walk_metropolis_hastings''!')
+ end
+ end
+ % run MCMC sampling
+ posterior_sampler_options = options_mom_.posterior_sampler_options.current_options;
+ posterior_sampler_options.invhess = invhess;
+ [posterior_sampler_options, options_mom_, bayestopt_] = check_posterior_sampler_options(posterior_sampler_options, M_.fname, M_.dname, options_mom_, BoundsInfo, bayestopt_,'method_of_moments');
+ options_mom_.posterior_sampler_options.current_options = posterior_sampler_options; % store current options
+ if options_mom_.mh_replic>0
+ posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,BoundsInfo,oo_.mom.data_moments,oo_.mom.weighting_info,options_mom_,M_,estim_params_,bayestopt_,oo_,'method_of_moments::mcmc');
+ end
+ CutSample(M_, options_mom_, 'method_of_moments::mcmc'); % discard first mh_drop percent of the draws
+ if options_mom_.mh_posterior_mode_estimation
+ % skip optimizer-based mode-finding and instead compute the mode based on a run of a MCMC
+ [~,~,posterior_mode,~] = compute_mh_covariance_matrix(bayestopt_,M_.fname,M_.dname,'method_of_moments');
+ oo_.mom = fill_mh_mode(posterior_mode',NaN(length(posterior_mode),1),M_,options_mom_,estim_params_,bayestopt_,oo_.mom,'posterior');
+ return
+ else
+ % get stored results if required
+ if options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
+ oo_load_mh = load([M_.dname filesep 'method_of_moments' filesep M_.fname '_mom_results'],'oo_');
+ end
+ % convergence diagnostics
+ if ~options_mom_.nodiagnostic
+ if (options_mom_.mh_replic>0 || (options_mom_.load_mh_file && ~options_mom_.load_results_after_load_mh))
+ oo_.mom = mcmc_diagnostics(options_mom_, estim_params_, M_, oo_.mom);
+ elseif options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
+ if isfield(oo_load_mh.oo_.mom,'convergence')
+ oo_.mom.convergence = oo_load_mh.oo_.mom.convergence;
+ end
+ end
+ end
+ % statistics and plots for posterior draws
+ if options_mom_.mh_replic || (options_mom_.load_mh_file && ~options_mom_.load_results_after_load_mh)
+ [~,oo_.mom] = marginal_density(M_, options_mom_, estim_params_, oo_.mom, bayestopt_, 'method_of_moments');
+ oo_.mom = GetPosteriorParametersStatistics(estim_params_, M_, options_mom_, bayestopt_, oo_.mom, prior_dist_names);
+ if ~options_mom_.nograph
+ oo_.mom = PlotPosteriorDistributions(estim_params_, M_, options_mom_, bayestopt_, oo_.mom);
+ end
+ [oo_.mom.posterior.metropolis.mean,oo_.mom.posterior.metropolis.Variance] = GetPosteriorMeanVariance(M_,options_mom_.mh_drop);
+ elseif options_mom_.load_mh_file && options_mom_.load_results_after_load_mh
+ % load fields from previous MCMC run stored in results-file
+ field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
+ 'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
+ 'prior_density',...% fields set by PlotPosteriorDistributions
+ };
+ for field_iter=1:size(field_names,2)
+ if isfield(oo_load_mh.oo_.mom,field_names{1,field_iter})
+ oo_.mom.(field_names{1,field_iter}) = oo_load_mh.oo_.mom.(field_names{1,field_iter});
+ end
+ end
+ if isfield(oo_load_mh.oo_.mom,'MarginalDensity') && isfield(oo_load_mh.oo_.mom.MarginalDensity,'ModifiedHarmonicMean') % field set by marginal_density
+ oo_.mom.MarginalDensity.ModifiedHarmonicMean = oo_load_mh.oo_.mom.MarginalDensity.ModifiedHarmonicMean;
+ end
+ if isfield(oo_load_mh.oo_.mom,'posterior') && isfield(oo_load_mh.oo_.mom.posterior,'metropolis') % field set by GetPosteriorMeanVariance
+ oo_.mom.posterior.metropolis = oo_load_mh.oo_.mom.posterior.metropolis;
+ end
+ end
+ [error_flag,~,options_mom_]= metropolis_draw(1,options_mom_,estim_params_,M_);
+ if ~(~isempty(options_mom_.sub_draws) && options_mom_.sub_draws==0)
+ % THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
+ %if options_mom_.bayesian_irf
+ % if error_flag
+ % error('method_of_moments: Cannot compute the posterior IRFs!');
+ % end
+ % PosteriorIRF('posterior','method_of_moments::mcmc');
+ %end
+ % if options_mom_.moments_varendo
+ % if error_flag
+ % error('method_of_moments: Cannot compute the posterior moments for the endogenous variables!');
+ % end
+ % if options_mom_.load_mh_file && options_mom_.mh_replic==0 %user wants to recompute results
+ % [MetropolisFolder, info] = CheckPath('metropolis',options_mom_.dirname);
+ % if ~info
+ % generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'};
+ % for ii=1:length(generic_post_data_file_name)
+ % delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']);
+ % end
+ % % restore compatibility for loading pre-4.6.2 runs where estim_params_ was not saved; see 6e06acc7 and !1944
+ % NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_posterior_draws*' ]));
+ % if NumberOfDrawsFiles>0
+ % temp=load([M_.dname '/metropolis/' M_.fname '_posterior_draws1']);
+ % if ~isfield(temp,'estim_params_')
+ % for file_iter=1:NumberOfDrawsFiles
+ % save([M_.dname '/metropolis/' M_.fname '_posterior_draws' num2str(file_iter)],'estim_params_','-append')
+ % end
+ % end
+ % end
+ % end
+ % end
+ % oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
+ % end
+ else
+ fprintf('''sub_draws'' was set to 0. Skipping posterior computations.');
+ end
+ xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_.mom,options_);
+ end
+ % MAYBE USEFUL????
+ % % Posterior correlations
+ % ExtremeCorrBound=0.7;
+ % if ~isnan(ExtremeCorrBound)
+ % tril_para_correlation_matrix=tril(para_correlation_matrix,-1);
+ % [RowIndex,ColIndex]=find(abs(tril_para_correlation_matrix)>ExtremeCorrBound);
+ % ExtremeCorrParams=cell(length(RowIndex),3);
+ % for i=1:length(RowIndex)
+ % ExtremeCorrParams{i,1}=char(parameter_names(RowIndex(i),:));
+ % ExtremeCorrParams{i,2}=char(parameter_names(ColIndex(i),:));
+ % ExtremeCorrParams{i,3}=tril_para_correlation_matrix(RowIndex(i),ColIndex(i));
+ % end
+ % end
+ % disp(' ');
+ % disp(['Correlations of Parameters (at Posterior Mode) > ',num2str(ExtremeCorrBound)]);
+ % disp(ExtremeCorrParams)
end
% -------------------------------------------------------------------------
% display final estimation results
% -------------------------------------------------------------------------
+M_ = set_all_parameters(xparam1,estim_params_,M_); % update parameters
+[~, ~, ~, ~, ~, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.model_moments_params_derivs, oo_.mom.irf_model_varobs] = objective_function(xparam1, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state); % store final results in oo_.mom
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,stdh,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.J_test = mom.Jtest(xparam1, objective_function, oo_.mom.Q, oo_.mom.model_moments, oo_.mom.m_data, oo_.mom.data_moments, oo_.mom.weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
- % display comparison of model moments and data moments
- mom.display_comparison_moments(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
+elseif strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
+ if ~options_mom_.nograph
+ mom.graph_comparison_irfs(M_.matched_irfs,oo_.mom.irf_model_varobs,options_mom_.varobs_id,options_mom_.irf,options_mom_.relative_irf,M_.endo_names,M_.exo_names,M_.exo_names_tex,M_.dname,M_.fname,options_mom_.graph_format,options_mom_.TeX,options_mom_.nodisplay,options_mom_.figures.textwidth)
+ end
end
+% display comparison of model moments/irfs and data moments/irfs
+mom.display_comparison_moments_irfs(M_, options_mom_, oo_.mom.data_moments, oo_.mom.model_moments);
+% save results to _mom_results.mat
+save([M_.dname filesep 'method_of_moments' filesep M_.fname '_mom_results.mat'], 'oo_', 'options_mom_', 'M_', 'estim_params_', 'bayestopt_');
+
fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method);
% -------------------------------------------------------------------------
@@ -548,3 +848,6 @@ if isoctave && isfield(options_mom_, 'prior_restrictions') && ...
% See https://savannah.gnu.org/bugs/?43215
options_mom_.prior_restrictions.routine = [];
end
+if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && ~isempty(options_mom_.mom.irf_matching_file.path) && ~strcmp(options_mom_.mom.irf_matching_file.path,'.')
+ rmpath(options_mom_.irf_matching_file.path); % remove path to irf_matching_file
+end
\ No newline at end of file