From 540f0454d2a6c96de5a19c739040b654ba33826f Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 29 Jun 2020 07:40:34 +0200 Subject: [PATCH] Code Review of GMM routines - fix prefilter option - Implement iterative GMM --- matlab/dynare_config.m | 1 + matlab/method_of_moments.m | 994 ------------------ matlab/method_of_moments/method_of_moments.m | 906 ++++++++++++++++ .../method_of_moments_data_moments.m} | 36 +- .../method_of_moments_mode_check.m | 185 ++++ .../method_of_moments_objective_function.m | 214 ++++ ...thod_of_moments_optimal_weighting_matrix.m | 48 +- .../method_of_moments_standard_errors.m | 105 ++ matlab/method_of_moments_GMM.m | 224 ---- matlab/method_of_moments_SMM.m | 223 ---- matlab/method_of_moments_standard_errors.m | 105 -- tests/Makefile.am | 4 + .../method_of_moments/AnScho_MoM.mod | 10 +- .../estimation/method_of_moments/RBC_MoM.mod | 101 +- .../method_of_moments/RBC_MoM_SMM_ME.mod | 189 ++++ .../method_of_moments/RBC_MoM_common.inc | 80 ++ .../method_of_moments/RBC_MoM_prefilter.mod | 161 +++ 17 files changed, 1931 insertions(+), 1655 deletions(-) delete mode 100644 matlab/method_of_moments.m create mode 100644 matlab/method_of_moments/method_of_moments.m rename matlab/{method_of_moments_datamoments.m => method_of_moments/method_of_moments_data_moments.m} (64%) create mode 100644 matlab/method_of_moments/method_of_moments_mode_check.m create mode 100644 matlab/method_of_moments/method_of_moments_objective_function.m rename matlab/{ => method_of_moments}/method_of_moments_optimal_weighting_matrix.m (62%) create mode 100644 matlab/method_of_moments/method_of_moments_standard_errors.m delete mode 100644 matlab/method_of_moments_GMM.m delete mode 100644 matlab/method_of_moments_SMM.m delete mode 100644 matlab/method_of_moments_standard_errors.m create mode 100644 tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod create mode 100644 tests/estimation/method_of_moments/RBC_MoM_common.inc create mode 100644 tests/estimation/method_of_moments/RBC_MoM_prefilter.mod diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index ee74ebf11..4109f08dc 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -61,6 +61,7 @@ p = {'/distributions/' ; ... '/cli/' ; ... '/lmmcp/' ; ... '/optimization/' ; ... + '/method_of_moments/' ; ... '/discretionary_policy/' ; ... '/accessors/' ; ... '/modules/dseries/src/' ; ... diff --git a/matlab/method_of_moments.m b/matlab/method_of_moments.m deleted file mode 100644 index f933b4b93..000000000 --- a/matlab/method_of_moments.m +++ /dev/null @@ -1,994 +0,0 @@ -function [DynareResults, OptionsMoM] = method_of_moments(BayesInfo, DynareOptions, DynareResults, EstimatedParameters, Model, MatchedMoments, OptionsMoM) -%function [oo_, options_mom] = method_of_moments(M, options, oo, bayestopt, estim_params, matched_moments, 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 OptionsMoM structure -% - Carry over Options from the preprocessor -% - Other options that need to be initialized -% - Get variable orderings and state space representation -% Step 2: Checks and transformations for matched moments structure (preliminary) -% 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: Method of moments estimation: print some info, first-stage, and second-stage -% Step 8: Clean up -% ------------------------------------------------------------------------- -% This function is inspired by replication codes accompanied to the following papers: -% 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 BayesInfo: [structure] information about priors (bayestopt_) -% o DynareOptions: [structure] information about global options (options_) -% o DynareResults: [structure] storage for results (oo_) -% o EstimatedParameters: [structure] information about estimated parameters (estim_params_) -% o Model: [structure] information about model (M_) -% o MatchedMoments: [structure] information about selected moments to match in estimation (matched_moments_) -% o OptionsMoM: [structure] information about settings specified by the user (options_mom_) -% ------------------------------------------------------------------------- -% OUTPUTS -% o DynareResults: [structure] storage for results (oo_) -% o OptionsMoM: [structure] information about all used settings used in estimation (options_mom_) -% ------------------------------------------------------------------------- -% 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.m -% o evaluate_steady_state -% o get_all_parameters.m -% o makedataset.m -% o method_of_moments_datamoments.m -% 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 -% o list_of_parameters_calibrated_as_Inf -% o list_of_parameters_calibrated_as_NaN -% ========================================================================= -% Copyright (C) 2020 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 . -% ------------------------------------------------------------------------- -% Author(s): -% o Willi Mutschler (willi@mutschler.eu) -% o Johannes Pfeifer (jpfeifer@uni-koeln.de) -% ========================================================================= - -%% TO DO LIST -% - [ ] penalized estimation: how does penalized_estimator work? What is -% penalized_estimator? Not all optimizer make use of this...what is special -% about mode_compute=1 in objective functions. Do we need global objective_function_penalty_base in objective function -% - [ ] make csminwel work -% - [ ] why does lsqnonlin take less time in Andreasen toolbox? -% - [ ] recheck different optimizers if they are useful -% - [ ] test prefilter option (i.e. centered moments) -% - [ ] how to deal with logged_steady_state -% - [ ] mode_check plots? -% - [ ] test prior restrictions file -% - [ ] test user-specified weightning matrix -% - [ ] test non-symetric Sigma_e -% - [ ] which qz_criterium value? -% - [ ] are missing observations a problem? in method_of_moments_datamoments.m nan are replaced by mean of moment -% - [ ] check negative priors on std errors and above 1 for correlations -% - [ ] add measurement errors -% - [ ] add IRF matching -% - [ ] what are trend_coeffs and how to deal with them? -% - [ ] once interface is established: provide code to remove duplicate declared moments in matched moments block -% - [ ] test estimated_params_init(use calibration) -% - [ ] test estimated_params_bounds block -% - [ ] test what happens if all parameters will be estimated but some/all are not calibrated -% - [ ] speed up lyapunov equation by using doubling with old initial values -% - [ ] check what happens if parameters are set to INF or NAN -% - [ ] provide a table with dataMoments and final modelMoments -% - [ ] check smm at order > 3 without pruning -% - [ ] provide option to use analytical derivatives to compute std errors (similar to what we already do in identification) -% - [ ] add Bayesian GMM/SMM estimation -% - [ ] useautocorr -% - [ ] instead of mom_steps, iterate over optimizers using additional_optimizer_steps - -% ------------------------------------------------------------------------- -% Step 0: Check if required structures and options exist -% ------------------------------------------------------------------------- -if isempty(EstimatedParameters) % structure storing the info about estimated parameters in the estimated_params block - error('method_of_moments: You need to provide an ''estimated_params'' block') -end -if isempty(MatchedMoments) % 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(BayesInfo) && any(BayesInfo.pshape==0) && any(BayesInfo.pshape~=0) - error('method_of_moments: Estimation must be either fully classical or fully Bayesian. Maybe you forgot to specify a prior distribution.') -end -fprintf('\n==== Method of Moments Estimation ====\n\n') - -% ------------------------------------------------------------------------- -% Step 1a: Prepare OptionsMoM structure -% ------------------------------------------------------------------------- -% OptionsMoM 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(OptionsMoM.mom_method,'GMM') || strcmp(OptionsMoM.mom_method,'SMM') - OptionsMoM = set_default_option(OptionsMoM,'bartlett_kernel_lag',20); % bandwith in optimal weighting matrix - OptionsMoM = set_default_option(OptionsMoM,'order',1); % order of Taylor approximation in perturbation - OptionsMoM = set_default_option(OptionsMoM,'penalized_estimator',false); % @wmutschl: provide description - OptionsMoM = set_default_option(OptionsMoM,'pruning',true); % use pruned state space system at higher-order - OptionsMoM = set_default_option(OptionsMoM,'verbose',false); % display and store intermediate estimation results - OptionsMoM = set_default_option(OptionsMoM,'weighting_matrix','identity_matrix'); % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename - OptionsMoM = set_default_option(OptionsMoM,'additional_optimizer_steps',[]); % Number of iterations in the steps of the 2-step feasible method of moments - % Checks for perturbation order - if OptionsMoM.order < 1 - error('method_of_moments:: The order of the Taylor approximation cannot be 0!') - elseif OptionsMoM.order > 2 && (~isfield(OptionsMoM,'k_order_solver') || ~OptionsMoM.k_order_solver) - fprintf('method_of_moments: For perturbation order k>2, we add the k_order_solver option.\n'); - OptionsMoM.k_order_solver = 1; - end -end - -if strcmp(OptionsMoM.mom_method,'SMM') - objective_function = str2func('method_of_moments_SMM'); - OptionsMoM = set_default_option(OptionsMoM,'bounded_shock_support',false); % trim shocks in simulation to +- 2 stdev - OptionsMoM = set_default_option(OptionsMoM,'drop',500); % number of periods dropped at beginning of simulation - OptionsMoM = set_default_option(OptionsMoM,'seed',24051986); % seed used in simulations - OptionsMoM = set_default_option(OptionsMoM,'simulation_multiple',5); % multiple of the data length used for simulation - if OptionsMoM.simulation_multiple < 1 - fprintf('The simulation horizon is shorter than the data. Dynare resets the simulation_multiple to 2.\n') - OptionsMoM.smm.simulation_multiple = 2; - end -end -if strcmp(OptionsMoM.mom_method,'GMM') - objective_function = str2func('method_of_moments_GMM'); - % Check for pruning with GMM at higher order - if OptionsMoM.order > 1 && ~OptionsMoM.pruning - fprintf('GMM at higher order only works with pruning, so we set pruning option to 1.\n'); - OptionsMoM.pruning = true; - end -end - -% General options that can be set by the user in the mod file, otherwise default values are provided -OptionsMoM = set_default_option(OptionsMoM,'dirname',Model.fname); % directory in which to store estimation output -OptionsMoM = set_default_option(OptionsMoM,'graph_format','eps'); % specify the file format(s) for graphs saved to disk -OptionsMoM = set_default_option(OptionsMoM,'nodisplay',false); % do not display the graphs, but still save them to disk -OptionsMoM = set_default_option(OptionsMoM,'nograph',false); % do not create graphs (which implies that they are not saved to the disk nor displayed) -OptionsMoM = set_default_option(OptionsMoM,'noprint',false); % do not print output to console -OptionsMoM = set_default_option(OptionsMoM,'plot_priors',true); % control plotting of priors -OptionsMoM = set_default_option(OptionsMoM,'prior_trunc',1e-10); % probability of extreme values of the prior density that is ignored when computing bounds for the parameters -OptionsMoM = set_default_option(OptionsMoM,'TeX',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 -OptionsMoM = set_default_option(OptionsMoM,'first_obs',1); % number of first observation -OptionsMoM = set_default_option(OptionsMoM,'logdata',false); % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data) -OptionsMoM = set_default_option(OptionsMoM,'loglinear',false); % computes a log-linear approximation of the model instead of a linear approximation -OptionsMoM = set_default_option(OptionsMoM,'nobs',NaN); % number of observations -OptionsMoM = set_default_option(OptionsMoM,'prefilter',false); % demean each data series by its empirical mean and use centered moments -OptionsMoM = set_default_option(OptionsMoM,'xls_sheet',1); % name of sheet with data in Excel -OptionsMoM = set_default_option(OptionsMoM,'xls_range',''); % range of data in Excel sheet -% Recursive estimation and forecast are not supported -if numel(OptionsMoM.nobs)>1 - error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''.'); -end -if numel(OptionsMoM.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 -if strcmp(OptionsMoM.mom_method, 'GMM') - OptionsMoM = set_default_option(OptionsMoM,'analytic_derivation',0); % use analytic derivatives to compute standard errors for GMM -elseif isfield(OptionsMoM,'analytic_derivation') - fprintf('Only GMM supports analytic derivation to compute standard errors, we reset ''analytic_derivation'' to 0.\n') - OptionsMoM.analytic_derivation = 0; -else - OptionsMoM.analytic_derivation = 0; -end -OptionsMoM = set_default_option(OptionsMoM,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons -OptionsMoM = set_default_option(OptionsMoM,'mode_compute',13); % specifies the optimizer for minimization of moments distance -OptionsMoM = set_default_option(OptionsMoM,'optim_opt',[]); % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute -OptionsMoM = set_default_option(OptionsMoM,'silent_optimizer',false); % run minimization of moments distance silently without displaying results or saving files in between -if ~isfield(OptionsMoM,'dynatol') - OptionsMoM.dynatol = {}; -end -OptionsMoM.dynatol = set_default_option(OptionsMoM.dynatol,'f', 1e-5);% convergence criterion on function value for numerical differentiation -OptionsMoM.dynatol = set_default_option(OptionsMoM.dynatol,'x', 1e-5);% convergence criterion on funciton input for numerical differentiation - -% Numerical algorithms options that can be set by the user in the mod file, otherwise default values are provided -OptionsMoM = set_default_option(OptionsMoM,'aim_solver',false); % Use AIM algorithm to compute perturbation approximation -OptionsMoM = set_default_option(OptionsMoM,'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 -OptionsMoM = set_default_option(OptionsMoM,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm -OptionsMoM = set_default_option(OptionsMoM,'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 -OptionsMoM = set_default_option(OptionsMoM,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm -OptionsMoM = set_default_option(OptionsMoM,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the cycle reduction algorithm -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_db',false); % doubling algorithm (disclyap_fast) to solve Lyapunov equation to compute variance-covariance matrix of state variables -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_fp',false); % fixed-point algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_srs',false); % square-root-solver (dlyapchol) algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_complex_threshold',1e-15); % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_fixed_point_tol',1e-10); % convergence criterion used in the fixed point Lyapunov solver -OptionsMoM = set_default_option(OptionsMoM,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm -OptionsMoM = set_default_option(OptionsMoM,'sylvester_fp',false); % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models -OptionsMoM = set_default_option(OptionsMoM,'sylvester_fixed_point_tol',1e-12); % convergence criterion used in the fixed point Sylvester solver -OptionsMoM = set_default_option(OptionsMoM,'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 [IS THIS CORRET @wmutschl] -OptionsMoM = set_default_option(OptionsMoM,'qz_zero_threshold',1e-6); % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition - -% ------------------------------------------------------------------------- -% Step 1b: Options that are set by the preprocessor and (probably) need to be carried over -% ------------------------------------------------------------------------- - -% options related to VAROBS -if ~isfield(DynareOptions,'varobs') - error('method_of_moments: VAROBS statement is missing!') -else - OptionsMoM.varobs = DynareOptions.varobs; % observable variables in declaration order - OptionsMoM.obs_nbr = length(OptionsMoM.varobs);% number of observed variables - % Check that each declared observed variable is also an endogenous variable - for i = 1:OptionsMoM.obs_nbr - if ~any(strcmp(OptionsMoM.varobs{i},Model.endo_names)) - error(['method_of_moments: Unknown variable (' OptionsMoM.varobs{i} ')!']) - end - end - - % Check that a variable is not declared as observed more than once. - if length(unique(OptionsMoM.varobs))1 - error(['method_of_moments: A variable cannot be declared as observed more than once (' OptionsMoM.varobs{i} ')!']) - end - end - end -end - -% options related to variable declarations -if ~isfield(DynareOptions,'trend_coeffs') - %BayesInfo.with_trend = 0; -else - error('method_of_moments: %s does not allow for trend in data',OptionsMoM.mom_method) -end - -% options related to estimated_params and estimated_params_init -OptionsMoM.use_calibration_initialization = DynareOptions.use_calibration_initialization; - -% options related to model; block -OptionsMoM.linear = DynareOptions.linear; -OptionsMoM.use_dll = DynareOptions.use_dll; -OptionsMoM.block = DynareOptions.block; -OptionsMoM.bytecode = DynareOptions.bytecode; - -% options related to steady; command -OptionsMoM.homotopy_force_continue = DynareOptions.homotopy_force_continue; -OptionsMoM.homotopy_mode = DynareOptions.homotopy_mode; -OptionsMoM.homotopy_steps = DynareOptions.homotopy_steps; -OptionsMoM.logged_steady_state = DynareOptions.logged_steady_state; % @wmutschl: when and how does this get set? -OptionsMoM.markowitz = DynareOptions.markowitz; -OptionsMoM.solve_algo = DynareOptions.solve_algo; -OptionsMoM.solve_tolf = DynareOptions.solve_tolf; -OptionsMoM.steady = DynareOptions.steady; -OptionsMoM.steadystate = DynareOptions.steadystate; -OptionsMoM.steadystate_flag = DynareOptions.steadystate_flag; - -% options related to dataset -OptionsMoM.dataset = DynareOptions.dataset; -OptionsMoM.initial_period = DynareOptions.initial_period; - -% options related to endogenous prior restrictions -OptionsMoM.endogenous_prior_restrictions.irf = {}; -OptionsMoM.endogenous_prior_restrictions.moment = {}; -if ~isempty(DynareOptions.endogenous_prior_restrictions.irf) && ~isempty(DynareOptions.endogenous_prior_restrictions.moment) - fprintf('Endogenous prior restrictions are not supported yet and will be skipped.\n') -end - - -% ------------------------------------------------------------------------- -% Step 1c: Options related to optimizers -% ------------------------------------------------------------------------- -% mode_compute = 2 -OptionsMoM.saopt = DynareOptions.saopt; -% mode_compute = 4 -OptionsMoM.csminwel = DynareOptions.csminwel; -if OptionsMoM.mode_compute == 4 - error('method_of_moments:optimizer','method_of_moments: csminwel optimizer (mode_compute=4) is not yet supported (due to penalized_objective handling).\n Choose a different optimizer, e.g. lsqnonlin (mode_compute=13), fminsearch (mode_compute=7), SolveOpt (mode_compute=101).') -end -% mode_compute = 5 (not yet) -if OptionsMoM.mode_compute == 5 - error('method_of_moments:optimizer','method_of_moments: newrat optimizer (mode_compute=5) is not yet supported.\n Choose a different optimizer, e.g. lsqnonlin (mode_compute=13), fminsearch (mode_compute=7), SolveOpt (mode_compute=101).') -end -% mode_compute = 6 -if OptionsMoM.mode_compute == 6 - error('method_of_moments:optimizer','method_of_moments: mode_compute=6 is not compatible with a method of moments estimation.\n Choose a different optimizer, e.g. lsqnonlin (mode_compute=13), fminsearch (mode_compute=7), SolveOpt (mode_compute=101).') -end -% mode_compute = 8 -OptionsMoM.simplex = DynareOptions.simplex; -% mode_compute = 9 -OptionsMoM.cmaes = DynareOptions.cmaes; -% mode_compute = 10 -OptionsMoM.simpsa = DynareOptions.simpsa; -% mode_compute = 11 -if OptionsMoM.mode_compute == 11 - error('method_of_moments:optimizer','method_of_moments: mode_compute=11 is not compatible with a method of moments estimation.\n Choose a different optimizer, e.g. lsqnonlin (mode_compute=13), fminsearch (mode_compute=7), SolveOpt (mode_compute=101).') -end -% mode_compute = 12 -OptionsMoM.particleswarm = DynareOptions.particleswarm; -if OptionsMoM.mode_compute == 12 - error('method_of_moments:optimizer','method_of_moments: mode_compute=12 is not yet supported (due to penalized_objective handling).\n Choose a different optimizer, e.g. lsqnonlin (mode_compute=13), fminsearch (mode_compute=7), SolveOpt (mode_compute=101).') -end -% mode_compute = 101 -OptionsMoM.solveopt = DynareOptions.solveopt; - -OptionsMoM.gradient_method = DynareOptions.gradient_method; -OptionsMoM.gradient_epsilon = DynareOptions.gradient_epsilon; - -% ------------------------------------------------------------------------- -% Step 1d: Other options that need to be initialized -% ------------------------------------------------------------------------- -OptionsMoM.initialize_estimated_parameters_with_the_prior_mode = 0; % needed by set_prior.m -OptionsMoM.figures.textwidth = 0.8; %needed by plot_priors.m -OptionsMoM.ramsey_policy = 0; % needed by evaluate_steady_state -OptionsMoM.debug = false; %needed by resol.m -OptionsMoM.risky_steadystate = false; %needed by resol -OptionsMoM.threads = DynareOptions.threads; %needed by resol -OptionsMoM.jacobian_flag = true; -OptionsMoM.gstep = DynareOptions.gstep; -OptionsMoM.solve_tolf = DynareOptions.solve_tolf; -OptionsMoM.solve_tolx = DynareOptions.solve_tolx; - - -% 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; -%OptionsMoM = set_default_option(OptionsMoM,'endo_vars_for_moment_computations_in_estimation',[]); - -% ------------------------------------------------------------------------- -% Step 1e: Get variable orderings and state space representation -% ------------------------------------------------------------------------- -DynareResults.dr = set_state_space(DynareResults.dr,Model,OptionsMoM); -% Get index of observed variables in DR order -DynareResults.dr.obs_var = []; -for i=1:OptionsMoM.obs_nbr - DynareResults.dr.obs_var = [DynareResults.dr.obs_var; find(strcmp(OptionsMoM.varobs{i}, Model.endo_names(DynareResults.dr.order_var)))]; -end - -% ------------------------------------------------------------------------- -% Step 2: Checks and transformations for matched moments structure (preliminary) -% ------------------------------------------------------------------------- -% Note that we do not have a preprocessor interface yet for this, so this -% will need much improvement later on. @wmutschl -if strcmp(OptionsMoM.mom_method, 'GMM') - % Initialize indices - OptionsMoM.index.E_y = false(OptionsMoM.obs_nbr,1); %unconditional first order product moments - OptionsMoM.index.E_yy = false(OptionsMoM.obs_nbr,OptionsMoM.obs_nbr); %unconditional second order product moments - OptionsMoM.index.E_yyt = false(OptionsMoM.obs_nbr,OptionsMoM.obs_nbr,0); %unconditional temporal second order product moments - OptionsMoM.index.E_y_pos = zeros(OptionsMoM.obs_nbr,1); %position in matched moments block - OptionsMoM.index.E_yy_pos = zeros(OptionsMoM.obs_nbr,OptionsMoM.obs_nbr); %position in matched moments block - OptionsMoM.index.E_yyt_pos = zeros(OptionsMoM.obs_nbr,OptionsMoM.obs_nbr,0); %position in matched moments block -end - -for jm=1:size(MatchedMoments,1) - % higher-order product moments not supported yet for GMM - if strcmp(OptionsMoM.mom_method, 'GMM') && sum(MatchedMoments{jm,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 if declared variables are also observed (needed as otherwise the dataset variables won't coincide) - if any(~ismember(DynareResults.dr.inv_order_var(MatchedMoments{jm,1})', DynareResults.dr.obs_var)) - error('method_of_moments: Variables in row %d in ''matched_moments'' block need to be declared as VAROBS.', jm) - end - - if strcmp(OptionsMoM.mom_method, 'GMM') - % Check (for now) that only lags are declared - if any(MatchedMoments{jm,2}>0) - error('method_of_moments: Leads in row %d in the ''matched_moments'' block are not supported for GMM, shift the moments and declare only lags.', jm) - end - % Check (for now) that first declared variable has zero lag - if MatchedMoments{jm,2}(1)~=0 - error('method_of_moments: The first variable declared in row %d in the ''matched_moments'' block is not allowed to have a lead or lag for GMM;\n reorder the variables in the row such that the first variable has zero lag!',jm) - end - vars = DynareResults.dr.inv_order_var(MatchedMoments{jm,1})'; - if sum(MatchedMoments{jm,3}) == 1 - % First-order product moment - vpos = (DynareResults.dr.obs_var == vars); - OptionsMoM.index.E_y(vpos,1) = true; - OptionsMoM.index.E_y_pos(vpos,1) = jm; - elseif sum(MatchedMoments{jm,3}) == 2 - % Second-order product moment - idx1 = (DynareResults.dr.obs_var == vars(1)); - idx2 = (DynareResults.dr.obs_var == vars(2)); - lag1 = MatchedMoments{jm,2}(1); - lag2 = MatchedMoments{jm,2}(2); - if lag1==0 && lag2==0 % contemporenous covariance matrix - OptionsMoM.index.E_yy(idx1,idx2) = true; - OptionsMoM.index.E_yy(idx2,idx1) = true; - OptionsMoM.index.E_yy_pos(idx1,idx2) = jm; - OptionsMoM.index.E_yy_pos(idx2,idx1) = jm; - elseif lag1==0 && lag2 < 0 - OptionsMoM.index.E_yyt(idx1,idx2,-lag2) = true; - OptionsMoM.index.E_yyt_pos(idx1,idx2,-lag2) = jm; - end - end - end -end - -% @wmutschl: add check for duplicate moments by using the cellfun and unique functions -if strcmp(OptionsMoM.mom_method,'GMM') - %Remove duplicate elements - UniqueMomIdx = [nonzeros(OptionsMoM.index.E_y_pos); nonzeros(triu(OptionsMoM.index.E_yy_pos)); nonzeros(OptionsMoM.index.E_yyt_pos)]; - DuplicateMoms = setdiff(1:size(MatchedMoments,1),UniqueMomIdx); - if ~isempty(DuplicateMoms) - fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows: %s.\n',num2str(DuplicateMoms)) - end - %reorder MatchedMoments to be compatible with OptionsMoM.index - MatchedMoments = MatchedMoments(UniqueMomIdx,:); -else - fprintf('For SMM we do not check yet for duplicate moment declarations in ''matched_moments'' block. You need to check this manually.\n\n') -end - -% Check if both prefilter and first moments were specified -OptionsMoM.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,MatchedMoments(:,3)))'; -if OptionsMoM.prefilter && ~isempty(OptionsMoM.first_moment_indicator) - fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block in rows: %s.\n',num2str(OptionsMoM.index.E_y_pos')); - MatchedMoments = MatchedMoments(OptionsMoM.first_moment_indicator,:); %remove first moments entries - OptionsMoM.first_moment_indicator = []; -end -OptionsMoM.mom_nbr = size(MatchedMoments,1); - - - -% Get maximum lag number for autocovariances/autocorrelations -OptionsMoM.ar = max(cellfun(@max,MatchedMoments(:,2))) - min(cellfun(@min,MatchedMoments(:,2))); - -% ------------------------------------------------------------------------- -% Step 3: Checks and transformations for estimated parameters, priors, and bounds -% ------------------------------------------------------------------------- - -% Set priors over the estimated parameters -if ~isempty(EstimatedParameters) && ~(isfield(EstimatedParameters,'nvx') && (size(EstimatedParameters.var_exo,1)+size(EstimatedParameters.var_endo,1)+size(EstimatedParameters.corrx,1)+size(EstimatedParameters.corrn,1)+size(EstimatedParameters.param_vals,1))==0) - [xparam0, EstimatedParameters, BayesInfo, lb, ub, Model] = set_prior(EstimatedParameters, Model, OptionsMoM); -end - -% Check measurement errors -if EstimatedParameters.nvn || EstimatedParameters.ncn - error('method_of_moments: moment estimation does not support measurement error(s) yet. Please specifiy them as a structural shock.') -end - -% Check if a _prior_restrictions.m file exists -if exist([Model.fname '_prior_restrictions.m']) - OptionsMoM.prior_restrictions.status = 1; - OptionsMoM.prior_restrictions.routine = str2func([Model.fname '_prior_restrictions']); -end - -% Check on specified priors and penalized estimation -if ~isempty(EstimatedParameters) && ~(isfield(EstimatedParameters,'nvx') && (size(EstimatedParameters.var_exo,1)+size(EstimatedParameters.var_endo,1)+size(EstimatedParameters.corrx,1)+size(EstimatedParameters.corrn,1)+size(EstimatedParameters.param_vals,1))==0) - if any(BayesInfo.pshape > 0) % prior specified - if any(setdiff([0;BayesInfo.pshape],[0,3])) - if ~OptionsMoM.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',OptionsMoM.mom_method) - OptionsMoM.penalized_estimator=1; - end - fprintf('\nNon-normal priors specified. %s with penalty uses a Laplace type of approximation.\n',OptionsMoM.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') - end - if any(isinf(BayesInfo.p2)) - inf_var_pars=BayesInfo.name(isinf(BayesInfo.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') - warning('off','MATLAB:singularMatrix'); - end - end -end -if OptionsMoM.penalized_estimator && OptionsMoM.mode_compute==13 - error('method_of_moments: Penalized estimator option is not compatible with mode_compute=13. Deactivate penalized_estimator, don''t use priors, or choose a different optimizer.') -end - - -% Check for calibrated covariances before updating parameters -if ~isempty(EstimatedParameters) && ~(isfield(EstimatedParameters,'nvx') && sum(EstimatedParameters.nvx+EstimatedParameters.nvn+EstimatedParameters.ncx+EstimatedParameters.ncn+EstimatedParameters.np)==0) - EstimatedParameters = check_for_calibrated_covariances(xparam0,EstimatedParameters,Model); -end - -% Checks on parameter calibration and initialization -xparam1_calib = get_all_parameters(EstimatedParameters,Model); %get calibrated parameters -if ~any(isnan(xparam1_calib)) %all estimated parameters are calibrated - EstimatedParameters.full_calibration_detected=1; -else - EstimatedParameters.full_calibration_detected=0; -end -if OptionsMoM.use_calibration_initialization %set calibration as starting values - if ~isempty(BayesInfo) && all(BayesInfo.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 properly initialized.',OptionsMoM.mom_method) - else - [xparam0,EstimatedParameters]=do_parameter_initialization(EstimatedParameters,xparam1_calib,xparam0); %get explicitly initialized parameters that have precedence to calibrated values - end -end - -% Check on initialization @wmutschl: why without penalty? -if ~isempty(BayesInfo) && all(BayesInfo.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 ',OptionsMoM.mom_method) -end - -% Set and check parameter bounds -if ~isempty(EstimatedParameters) && ~(all(strcmp(fieldnames(EstimatedParameters),'full_calibration_detected')) || (isfield(EstimatedParameters,'nvx') && sum(EstimatedParameters.nvx+EstimatedParameters.nvn+EstimatedParameters.ncx+EstimatedParameters.ncn+EstimatedParameters.np)==0)) - if ~isempty(BayesInfo) && any(BayesInfo.pshape > 0) - % Plot prior densities - if ~OptionsMoM.nograph && OptionsMoM.plot_priors - plot_priors(BayesInfo,Model,EstimatedParameters,OptionsMoM) - end - % Set prior bounds - Bounds = prior_bounds(BayesInfo, OptionsMoM.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. - Bounds.lb = lb; - Bounds.ub = ub; - end - % Test if initial values of the estimated parameters are all between the prior lower and upper bounds - if OptionsMoM.use_calibration_initialization - try - check_prior_bounds(xparam0,Bounds,Model,EstimatedParameters,OptionsMoM,BayesInfo) - catch - e = lasterror(); - fprintf('Cannot use parameter values from calibration as they violate the prior bounds.') - rethrow(e); - end - else - check_prior_bounds(xparam0,Bounds,Model,EstimatedParameters,OptionsMoM,BayesInfo) - end -end - -% Check symmetric Sigma_e -Sigma_e_non_zero_rows = find(~all(Model.Sigma_e==0,1)); -Sigma_e_non_zero_columns = find(~all(Model.Sigma_e==0,2)); -if ~isequal(Sigma_e_non_zero_rows,Sigma_e_non_zero_columns') - error('method_of_moments: Structual error matrix not symmetric') -end -if isfield(EstimatedParameters,'var_exo') && ~isempty(EstimatedParameters.var_exo) - EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness = union(Sigma_e_non_zero_rows,EstimatedParameters.var_exo(:,1)); -else - EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness = Sigma_e_non_zero_rows; -end - -% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file). -Model.sigma_e_is_diagonal = true; -if EstimatedParameters.ncx || any(nnz(tril(Model.Correlation_matrix,-1))) || isfield(EstimatedParameters,'calibrated_covariances') - Model.sigma_e_is_diagonal = false; -end - -% Provide some warnings on standard errors and correlations of shocks -if any(BayesInfo.pshape) - if EstimatedParameters.nvx && any(BayesInfo.p3(1:EstimatedParameters.nvx)<0) - warning('Your prior allows for negative standard deviations for structural shocks. It is recommended to change your prior.') - end - offset=EstimatedParameters.nvx; - if EstimatedParameters.nvn && any(BayesInfo.p3(1+offset:offset+EstimatedParameters.nvn)<0) - warning('Your prior allows for negative standard deviations for measurement error. It is recommended to change your prior.') - end - offset = EstimatedParameters.nvx+EstimatedParameters.nvn; - if EstimatedParameters.ncx && (any(BayesInfo.p3(1+offset:offset+EstimatedParameters.ncx)<-1) || any(BayesInfo.p4(1+offset:offset+EstimatedParameters.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 = EstimatedParameters.nvx+EstimatedParameters.nvn+EstimatedParameters.ncx; - if EstimatedParameters.ncn && (any(BayesInfo.p3(1+offset:offset+EstimatedParameters.ncn)<-1) || any(BayesInfo.p4(1+offset:offset+EstimatedParameters.ncn)>1)) - warning('Your prior allows for correlations between measurement errors larger than +-1 and will not integrate to 1 due to truncation. Please change your prior') - end -end - -% storing prior parameters in MoM info structure for penalized minimization -DynareResults.prior.pnames = {''; 'beta'; 'gamm'; 'norm'; 'invg'; 'unif'; 'invg2'; ''; 'weibl'}; -DynareResults.prior.p1 = BayesInfo.p1; -DynareResults.prior.p2 = BayesInfo.p2; -% DynareResults.prior.mode = BayesInfo.p5; -% DynareResults.prior.variance = diag(BayesInfo.p2.^2); -% DynareResults.prior.hyperparameters.first = BayesInfo.p6; -% DynareResults.prior.hyperparameters.second = BayesInfo.p7; - - -% Set all parameters -Model = set_all_parameters(xparam0,EstimatedParameters,Model); - - -% ------------------------------------------------------------------------- -% Step 4: Checks and transformations for data -% ------------------------------------------------------------------------- - -% Check if datafile has same name as mod file -if ~isempty(OptionsMoM.datafile) - [~,name,~] = fileparts(OptionsMoM.datafile); - if strcmp(name,Model.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.') - end -end - -% Build dataset -[DynareDataset, ~, ~] = makedataset(OptionsMoM); - -% set options for old interface from the ones for new interface -if ~isempty(DynareDataset) - OptionsMoM.nobs = DynareDataset.nobs; -end - -% missing observations are not supported yet -if any(any(isnan(DynareDataset.data))) - error('method_of_moments: missing observations are not supported') -end - -% Check length of data for estimation of second moments -if OptionsMoM.ar > OptionsMoM.nobs+1 - error('method_of_moments: Data set is too short to compute second moments'); -end - -% Get data moments for the method of moments -[DynareResults.mom.dataMoments, DynareResults.mom.m_data] = method_of_moments_datamoments(DynareDataset.data, DynareResults, MatchedMoments, OptionsMoM); - -if OptionsMoM.prefilter - if sum(abs(DynareDataset.mean))/DynareDataset.nobs >1e-9 - fprintf('The mean of the data is:\n') - disp(DynareDataset.mean); - error('method_of_moments: You are trying to perform a method of moments estimation with centered moments (prefilter=1) using uncentered data.') - end -elseif ~isempty(OptionsMoM.first_moment_indicator) - if sum(abs(DynareResults.mom.dataMoments(OptionsMoM.first_moment_indicator)))/sum(OptionsMoM.first_moment_indicator) <1e-2 - fprintf('The mean of the data for which Dynare tries to match first moments is:\n') - disp(DynareResults.mom.dataMoments(OptionsMoM.first_moment_indicator)'); - warning('method_of_moments:data_mean_zero','method_of_moments: You are trying to perform a method of moments estimation with uncentered moments (prefilter=0),\n but the data are (almost) mean 0. Check if this is desired.') - end -end - -% Get shock series for SMM and set variance correction factor -if strcmp(OptionsMoM.mom_method,'SMM') - OptionsMoM.long = round(OptionsMoM.simulation_multiple*OptionsMoM.nobs); - OptionsMoM.variance_correction_factor = (1+1/OptionsMoM.simulation_multiple); - % draw shocks for SMM - smmstream = RandStream('mt19937ar','Seed',OptionsMoM.seed); - temp_shocks = randn(smmstream,OptionsMoM.long+OptionsMoM.drop,Model.exo_nbr); - if OptionsMoM.bounded_shock_support == 1 - temp_shocks(temp_shocks>2) = 2; - temp_shocks(temp_shocks<-2) = -2; - end - OptionsMoM.shock_series = temp_shocks; -end - -% ------------------------------------------------------------------------- -% Step 5: checks for steady state at initial parameters -% ------------------------------------------------------------------------- - -% Check for logged steady state ...is this necessary @wmutschl -if OptionsMoM.logged_steady_state - DynareResults.dr.ys=exp(DynareResults.dr.ys); - DynareResults.steady_state=exp(DynareResults.steady_state); - OptionsMoM.logged_steady_state=0; -end - -% setting steadystate_check_flag option -if OptionsMoM.steadystate.nocheck - steadystate_check_flag = 0; -else - steadystate_check_flag = 1; -end - -% Check steady state at initial model parameter values -[DynareResults.steady_state, params, info] = evaluate_steady_state(DynareResults.steady_state,Model,OptionsMoM,DynareResults,steadystate_check_flag); -if info(1) - fprintf('\nmethod_of_moments: The steady state at the initial parameters cannot be computed.\n') - print_info(info, 0, OptionsMoM); -end - -try - % check if steady state solves static model - [DynareResults.steady_state, new_steady_params] = evaluate_steady_state(DynareResults.steady_state,Model,OptionsMoM,DynareResults,1); - - % check whether steady state file changes estimated parameters - if isfield(EstimatedParameters,'param_vals') && ~isempty(EstimatedParameters.param_vals) - Model0=Model; %store Model structure - old_steady_params=Model.params; %save initial parameters for check if steady state changes param values - - Model0.params(EstimatedParameters.param_vals(:,1))=Model0.params(EstimatedParameters.param_vals(:,1))*1.01; %vary parameters - [~, new_steady_params_2] = evaluate_steady_state(DynareResults.steady_state,Model0,OptionsMoM,DynareResults,1); - - changed_par_indices=find((old_steady_params(EstimatedParameters.param_vals(:,1))-new_steady_params(EstimatedParameters.param_vals(:,1))) ... - | (Model0.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 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.') - end - end - - % display warning if some parameters are still NaN - test_for_deep_parameters_calibration(Model); - -catch % if check fails, provide info on using calibration if present - e = lasterror(); - if EstimatedParameters.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') - fprintf('If this is not a problem with the setting of options (check the error message below),\n') - fprintf('you should try using the calibrated version of the model as starting values. To do\n') - fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n') - fprintf('command (and after the estimated_params-block so that it does not get overwritten):\n'); - skipline(2); - end - rethrow(e); -end - -% If steady state of observed variables is non zero, set noconstant equal 0 -if (~OptionsMoM.loglinear && all(abs(DynareResults.steady_state(DynareResults.dr.order_var(DynareResults.dr.obs_var)))<1e-9)) || (OptionsMoM.loglinear && all(abs(log(DynareResults.steady_state(DynareResults.dr.order_var(DynareResults.dr.obs_var))))<1e-9)) - OptionsMoM.noconstant = 1; -else - OptionsMoM.noconstant = 0; - % If the data are prefiltered then there must not be constants in the - % measurement equation of the DSGE model - if OptionsMoM.prefilter - skipline() - disp('You have specified the option "prefilter" to demean your data but the') - disp('steady state of of the observed variables is non zero.') - disp('Either change the measurement equations, by centering the observed') - disp('variables in the model block, or drop the prefiltering.') - error('method_of_moments: The option ''prefilter'' is inconsistent with the non-zero mean measurement equations in the model.') - end -end - - -% ------------------------------------------------------------------------- -% Step 6: checks for objective function at initial parameters -% ------------------------------------------------------------------------- -try - % Check for NaN or complex values of moment-distance-funtion evaluated - % at initial parameters and identity weighting matrix - DynareResults.mom.Sw = eye(OptionsMoM.mom_nbr); - tic_id = tic; - [fval, info, exit_flag, DynareResults, Model, OptionsMoM] = feval(objective_function, xparam0, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); - if OptionsMoM.mode_compute == 13 - fval = fval'*fval; - end - elapsed_time = toc(tic_id); - if isnan(fval) - error('method_of_moments: The initial value of the target function is NaN') - elseif imag(fval) - error('method_of_moments: The initial value of the target function is complex') - end - if info(1) > 0 - disp('method_of_moments: Error in computing the target function for initial parameter values') - print_info(info, OptionsMoM.noprint, OptionsMoM) - end - if OptionsMoM.prefilter==1 - if (~OptionsMoM.loglinear && any(abs(DynareResults.steady_state(DynareResults.dr.order_var(DynareResults.dr.obs_var)))>1e-9)) || (OptionsMoM.loglinear && any(abs(log(DynareResults.steady_state(DynareResults.dr.order_var(DynareResults.dr.obs_var))))>1e-9)) - disp(['You are trying to estimate a model with a non zero steady state for the observed endogenous']) - disp(['variables using demeaned data!']) - error('method_of_moments: You should change something in your mod file...') - end - end - fprintf('Initial value of the moment objective function with identity weighting matrix: %6.4f \n\n', fval); - fprintf('Time required to compute target function once: %5.4f seconds \n', elapsed_time); - -catch % if check fails, provide info on using calibration if present - e = lasterror(); - if EstimatedParameters.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') - fprintf('If this is not a problem with the setting of options (check the error message below),\n') - fprintf('you should try using the calibrated version of the model as starting values. To do\n') - fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n') - fprintf('command (and after the estimated_params-block so that it does not get overwritten):\n'); - skipline(2); - end - rethrow(e); -end - -if OptionsMoM.mode_compute == 0 %We only report value of moments distance at initial value of the parameters - fprintf('no minimization of moments distance due to ''mode_compute=0''\n') - return -end - -% ------------------------------------------------------------------------- -% Step 7a: Method of moments estimation: print some info -% ------------------------------------------------------------------------- -fprintf('\n---------------------------------------------------\n') -if strcmp(OptionsMoM.mom_method,'SMM') - fprintf('Simulated method of moments with'); -elseif strcmp(OptionsMoM.mom_method,'GMM') - fprintf('General method of moments with'); -end -if OptionsMoM.prefilter - fprintf('\n - centered moments (prefilter=1)'); -else - fprintf('\n - uncentered moments (prefilter=0)'); -end -if OptionsMoM.penalized_estimator - fprintf('\n - penalized estimation using priors'); -end -if OptionsMoM.mode_compute == 1; fprintf('\n - optimizer (mode_compute=1): fmincon'); -elseif OptionsMoM.mode_compute == 2; fprintf('\n - optimizer (mode_compute=2): continuous simulated annealing'); -elseif OptionsMoM.mode_compute == 3; fprintf('\n - optimizer (mode_compute=3): fminunc'); -elseif OptionsMoM.mode_compute == 4; fprintf('\n - optimizer (mode_compute=4): csminwel'); -elseif OptionsMoM.mode_compute == 7; fprintf('\n - optimizer (mode_compute=7): fminsearch'); -elseif OptionsMoM.mode_compute == 8; fprintf('\n - optimizer (mode_compute=8): Dynare Nelder-Mead simplex'); -elseif OptionsMoM.mode_compute == 9; fprintf('\n - optimizer (mode_compute=9): CMA-ES'); -elseif OptionsMoM.mode_compute == 10; fprintf('\n - optimizer (mode_compute=10): simpsa'); -elseif OptionsMoM.mode_compute == 12; fprintf('\n - optimizer (mode_compute=12): particleswarm'); -elseif OptionsMoM.mode_compute == 101; fprintf('\n - optimizer (mode_compute=101): SolveOpt'); -elseif OptionsMoM.mode_compute == 102; fprintf('\n - optimizer (mode_compute=102): simulannealbnd'); -elseif OptionsMoM.mode_compute == 13; fprintf('\n - optimizer (mode_compute=13): lsqnonlin'); -end -if OptionsMoM.silent_optimizer - fprintf(' (silent)'); -end -fprintf('\n - perturbation order: %d', OptionsMoM.order) -if OptionsMoM.order > 1 && OptionsMoM.pruning -fprintf(' (with pruning)') -end -fprintf('\n - number of matched moments: %d', OptionsMoM.mom_nbr); -fprintf('\n - number of parameters: %d', length(xparam0)); -% Check if enough moments for estimation -if OptionsMoM.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') - -% ------------------------------------------------------------------------- -% Step 7b: Method of moments estimation: First-stage -% ------------------------------------------------------------------------- -fprintf('First-stage estimation\n'); -switch lower(OptionsMoM.weighting_matrix) - case 'identity_matrix' - fprintf(' - identity weighting matrix\n'); - DynareResults.mom.Sw = eye(length(DynareResults.mom.dataMoments)); - case 'diagonal' - %@wmutschl: better description in fprintf - fprintf(' - diagonal weighting matrix: diagonal of Newey-West estimator with lag order %d\n', OptionsMoM.bartlett_kernel_lag); - fprintf(' and data moments as estimate of unconditional moments\n'); - Wopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.dataMoments, OptionsMoM.bartlett_kernel_lag); - DynareResults.mom.Sw = chol(diag(diag(Wopt))); - case 'optimal' - %@wmutschl: better description in fprintf - fprintf(' - weighting matrix: optimal. At first-stage we use diagonal of Newey-West estimator with lag order %d\n', OptionsMoM.bartlett_kernel_lag); - fprintf(' and the data moments as initial estimate of unconditional moments\n'); - Wopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.dataMoments, OptionsMoM.bartlett_kernel_lag); - DynareResults.mom.Sw = chol(diag(diag(Wopt))); - otherwise %user specified matrix in file - fprintf(' - weighting matrix: user-specified\n'); - try - load(OptionsMoM.weighting_matrix,'weighting_matrix') - catch - error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',OptionsMoM.weighting_matrix,'.mat']) - end - [nrow, ncol] = size(weighting_matrix); - if ~isequal(nrow,ncol) && ~isequal(nrow,length(DynareResults.mom.dataMoments)) %check if square and right size - error(['method_of_moments: weighting_matrix must be square and have ',num2str(length(DynareResults.mom.dataMoments)),' rows and columns']) - end - try %check for positive definiteness - DynareResults.Sw = chol(weighting_matrix); - hsd = sqrt(diag(weighting_matrix)); - inv(weighting_matrix./(hsd*hsd'))./(hsd*hsd'); - catch - error('method_of_moments: Specified weighting_matrix is not positive definite') - end -end -Woptflag = 0; -xparam1 = xparam0; -for istep1 = 1:2 - [xparam1, fval, exitflag, hessian_mat, OptionsMoM] = dynare_minimize_objective(objective_function, xparam1, OptionsMoM.mode_compute, OptionsMoM, [Bounds.lb Bounds.ub], BayesInfo.name, BayesInfo, [],... - Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); - if OptionsMoM.mode_compute == 13 - fval = fval'*fval; - end - fprintf('\nIteration %d value of minimized moment''s distance target function: %f.\n',istep1,fval) - if OptionsMoM.verbose - DynareResults.mom=display_estimation_results_table(xparam1,NaN(size(xparam1)),Model,OptionsMoM,EstimatedParameters,BayesInfo,DynareResults.mom,DynareResults.prior.pnames,sprintf('%s (FIRST-STAGE ITERATION %d) verbose',OptionsMoM.mom_method,istep1),sprintf('verbose_%s_1st_stage_iter_%d',lower(OptionsMoM.mom_method),istep1)); - end -end - -% Update Model and DynareResults (in particular DynareResults.mom.modelMoments) -Model = set_all_parameters(xparam1,EstimatedParameters,Model); -[fval, ~, ~, DynareResults, ~, ~] = feval(objective_function, xparam1, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); -if OptionsMoM.mode_compute == 13 - fval = fval'*fval; -end - -% Compute Standard errors -SE_1 = method_of_moments_standard_errors(xparam1, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Woptflag); - -% Store first-stage results in output structure -DynareResults.mom = display_estimation_results_table(xparam1,SE_1,Model,OptionsMoM,EstimatedParameters,BayesInfo,DynareResults.mom,DynareResults.prior.pnames,sprintf('%s (FIRST-STAGE)',OptionsMoM.mom_method),sprintf('%s_1st_stage',lower(OptionsMoM.mom_method))); - -% ------------------------------------------------------------------------- -% Step 7c: Method of moments estimation: Second-stage -% ------------------------------------------------------------------------- -fprintf('Second-stage estimation\n'); -switch lower(OptionsMoM.weighting_matrix) - case 'identity_matrix' - fprintf(' - weighting matrix: identity\n'); - DynareResults.mom.Sw = eye(length(DynareResults.mom.dataMoments)); - case 'diagonal' - fprintf(' - weighting matrix: diagonal of Newey-West estimator with lag order %d\n', OptionsMoM.bartlett_kernel_lag); - fprintf(' and based on first-stage estimate of unconditional model moments\n'); - Wopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.modelMoments, OptionsMoM.bartlett_kernel_lag); - DynareResults.mom.Sw = chol(diag(diag(Wopt))); - case 'optimal' - fprintf(' - weighting matrix: Newey-West estimator with lag order %d\n', OptionsMoM.bartlett_kernel_lag); - fprintf(' and based on first-stage estimate of unconditional model moments\n'); - Wopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.modelMoments, OptionsMoM.bartlett_kernel_lag); - DynareResults.mom.Sw = chol(Wopt); - Woptflag = 1; - fprintf(' rank of optimal weighting matrix: %d\n',rank(Wopt)); - otherwise %keep user specified matrix in file - fprintf(' - weighting matrix: user-specified\n'); -end - -xparam2 = xparam1; -for istep2 = 1:2 - [xparam2, fval, exitflag, hessian_mat, OptionsMoM] = dynare_minimize_objective(objective_function, xparam2, OptionsMoM.mode_compute, OptionsMoM, [Bounds.lb Bounds.ub], BayesInfo.name, BayesInfo, [],... - Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); - if OptionsMoM.mode_compute == 13 - fval = fval'*fval; - end - fprintf('\n - iteration %d value of minimized moment''s distance target function: %f.\n',istep2,fval) - if OptionsMoM.verbose - DynareResults.mom=display_estimation_results_table(xparam2,NaN(size(xparam2)),Model,OptionsMoM,EstimatedParameters,BayesInfo,DynareResults.mom,DynareResults.prior.pnames,sprintf('%s (SECOND-STAGE ITERATION %d) verbose',OptionsMoM.mom_method,istep2),sprintf('verbose_%s_2nd_stage_iter_%d',lower(OptionsMoM.mom_method),istep2)); - end -end - -% Update Model and DynareResults (in particular DynareResults.mom.modelMoments) -Model = set_all_parameters(xparam2,EstimatedParameters,Model); -[fval, ~, ~, DynareResults, ~, ~] = feval(objective_function, xparam2, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); -if OptionsMoM.mode_compute == 13 - fval = fval'*fval; -end - -% Compute Standard errors -SE_2 = method_of_moments_standard_errors(xparam2, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Woptflag); - -% Store second-stage results in output structure -DynareResults.mom = display_estimation_results_table(xparam2,SE_2,Model,OptionsMoM,EstimatedParameters,BayesInfo,DynareResults.mom,DynareResults.prior.pnames,sprintf('%s (SECOND-STAGE)',OptionsMoM.mom_method),sprintf('%s_2nd_stage',lower(OptionsMoM.mom_method))); - -% Compute J statistic -if strcmp(OptionsMoM.mom_method,'SMM') - Variance_correction_factor = OptionsMoM.variance_correction_factor; -elseif strcmp(OptionsMoM.mom_method,'GMM') - Variance_correction_factor=1; -end -DynareResults.mom.J_test.j_stat = DynareDataset.nobs*Variance_correction_factor*fval; -DynareResults.mom.J_test.degrees_freedom = length(DynareResults.mom.modelMoments)-length(xparam2); -DynareResults.mom.J_test.p_val = 1-chi2cdf(DynareResults.mom.J_test.j_stat, DynareResults.mom.J_test.degrees_freedom); -fprintf('\n p-value of J-test: %f\n',DynareResults.mom.J_test.p_val) - -fprintf('\n==== Method of Moments Estimation Completed ====\n\n') - -% ------------------------------------------------------------------------- -% Step 8: Clean up -% ------------------------------------------------------------------------- -% restore warnings -warning('on','MATLAB:singularMatrix'); diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m new file mode 100644 index 000000000..684ed8632 --- /dev/null +++ b/matlab/method_of_moments/method_of_moments.m @@ -0,0 +1,906 @@ +function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, matched_moments_, options_mom_) +%function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, matched_moments_, 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 +% - Other options that need to be initialized +% - Get variable orderings and state space representation +% Step 2: Checks and transformations for matched moments structure (preliminary) +% 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: Method of moments estimation: print some info, first-stage, and second-stage +% Step 8: Clean up +% ------------------------------------------------------------------------- +% This function is inspired by replication codes accompanied to the following papers: +% 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 +% 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 +% ------------------------------------------------------------------------- +% OUTPUTS +% o oo_: [structure] storage for results (oo_) +% o options_mom_: [structure] information about all used settings used in estimation (options_mom_) +% ------------------------------------------------------------------------- +% 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 makedataset.m +% o method_of_moments_data_moments.m +% 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 +% ========================================================================= +% Copyright (C) 2020 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 . +% ------------------------------------------------------------------------- +% Author(s): +% o Willi Mutschler (willi@mutschler.eu) +% o Johannes Pfeifer (jpfeifer@uni-koeln.de) +% ========================================================================= + +%% TO DO LIST +% - [ ] why does lsqnonlin take less time in Andreasen toolbox? +% - [ ] test user-specified weightning matrix +% - [ ] which qz_criterium value? +% - [ ] document that in method_of_moments_data_moments.m NaN are replaced by mean of moment +% - [ ] add IRF matching +% - [ ] test estimated_params_bounds block +% - [ ] test what happens if all parameters will be estimated but some/all are not calibrated +% - [ ] speed up lyapunov equation by using doubling with old initial values +% - [ ] check smm at order > 3 without pruning +% - [ ] provide option to use analytical derivatives to compute std errors (similar to what we already do in identification) +% - [ ] add Bayesian GMM/SMM estimation +% - [ ] useautocorr + +% ------------------------------------------------------------------------- +% Stuff that needs to be taken care by the preprocessor +% ------------------------------------------------------------------------- + +if ~isfield(options_mom_.mom,'mom_method') % Estimation method, required + error('method_of_moments: You need to provide a ''mom_method''. Possible values are GMM or SMM.'); +else + if ~(strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')) + error('method_of_moments: The provided ''mom_method'' needs to be GMM or SMM.'); + end + objective_function = str2func('method_of_moments_objective_function'); +end +if ~isfield(options_mom_,'datafile') || isempty(options_mom_.datafile) % Filename of data, required + error('method_of_moments: You need to supply a ''datafile''.'); +end +% if order > 2 then we need to make sure that k_order_solver is selected +options_mom_.k_order_solver = options_.k_order_solver; +if isfield(options_mom_,'order') && options_mom_.order > 2 + if ~options_.k_order_solver + error('method_of_moments: For perturbation order k>2 the k_order_solver option needs to be added. Workaround: run stoch_simul(order=k) before method_of_moments.'); + end +end +% preprocessor needs to create all files as in stoch_simul(order=1|2|3) + +% ------------------------------------------------------------------------- +% 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 isempty(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') +end + +fprintf('\n==== Method of Moments (%s) Estimation ====\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); % @wmutschl: provide description + 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','identity_matrix'); % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1000); % scaling of weighting matrix + options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for standard error + options_mom_ = set_default_option(options_mom_,'order',1); % order of Taylor approximation in perturbation + options_mom_ = set_default_option(options_mom_,'pruning',true); % 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!') + end +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',5); % 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 2.\n') + options_mom_.mom.simulation_multiple = 2; + end +end +if strcmp(options_mom_.mom.mom_method,'GMM') + % Check for 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 +end + +% 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_.fname); % 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_,'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 + +% 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 loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data) +options_mom_ = set_default_option(options_mom_,'loglinear',false); % we do not allow it here, but it needs to be set for makedataset +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 +% 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 +if strcmp(options_mom_.mom.mom_method, 'GMM') + options_mom_ = set_default_option(options_mom_,'analytic_derivation',0); % use analytic derivatives to compute standard errors for GMM +elseif isfield(options_mom_,'analytic_derivation') + fprintf('Only GMM supports analytic derivation to compute standard errors, we reset ''analytic_derivation'' to 0.\n') + options_mom_.analytic_derivation = 0; +else + options_mom_.analytic_derivation = 0; +end +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_ = set_default_option(options_mom_,'mode_compute',13); % specifies the optimizer for minimization of moments distance +options_mom_ = set_default_option(options_mom_,'vector_output',false); % specifies the whether the objective function returns a vector +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_.solve_tolf = set_default_option(options_mom_,'solve_tolf', eps^(1/3));% convergence criterion on function value for steady state finding +options_mom_.solve_tolx = set_default_option(options_mom_,'solve_tolx', eps^(2/3));% convergence criterion on function input for steady state finding + +% 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 +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 [IS THIS CORRET @wmutschl] +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 + +% ------------------------------------------------------------------------- +% Step 1b: Options that are set by the preprocessor and (probably) need to be carried over +% ------------------------------------------------------------------------- + +% 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 + 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 + 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_.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 +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 + +options_mom_.mode_check = options_.mode_check; + +% ------------------------------------------------------------------------- +% Step 1c: Options related to optimizers +% ------------------------------------------------------------------------- +% 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; + +% ------------------------------------------------------------------------- +% 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_.debug = false; %needed by resol.m +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 (preliminary) +% ------------------------------------------------------------------------- +% Note that we do not have a preprocessor interface yet for this, so this +% will need much improvement later on. @wmutschl + +% Initialize indices +options_mom_.index.E_y = false(options_mom_.obs_nbr,1); %unconditional first order product moments +options_mom_.index.E_yy = false(options_mom_.obs_nbr,options_mom_.obs_nbr); %unconditional second order product moments +options_mom_.index.E_yyt = false(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %unconditional temporal second order product moments +options_mom_.index.E_y_pos = zeros(options_mom_.obs_nbr,1); %position in matched moments block +options_mom_.index.E_yy_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr); %position in matched moments block +options_mom_.index.E_yyt_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %position in matched moments block + +for jm=1:size(matched_moments_,1) + % higher-order product moments not supported yet for GMM + if strcmp(options_mom_.mom.mom_method, 'GMM') && sum(matched_moments_{jm,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 if declared variables are also observed (needed as otherwise the dataset variables won't coincide) + if any(~ismember(oo_.dr.inv_order_var(matched_moments_{jm,1})', oo_.dr.obs_var)) + error('method_of_moments: Variables in row %d in ''matched_moments'' block need to be declared as VAROBS.', jm) + end + + if strcmp(options_mom_.mom.mom_method, 'GMM') + % Check (for now) that only lags are declared + if any(matched_moments_{jm,2}>0) + error('method_of_moments: Leads in row %d in the ''matched_moments'' block are not supported for GMM, shift the moments and declare only lags.', jm) + end + % Check (for now) that first declared variable has zero lag + if matched_moments_{jm,2}(1)~=0 + error('method_of_moments: The first variable declared in row %d in the ''matched_moments'' block is not allowed to have a lead or lag for GMM;\n reorder the variables in the row such that the first variable has zero lag!',jm) + end + end + vars = oo_.dr.inv_order_var(matched_moments_{jm,1})'; + if sum(matched_moments_{jm,3}) == 1 + % First-order product moment + vpos = (oo_.dr.obs_var == vars); + options_mom_.index.E_y(vpos,1) = true; + options_mom_.index.E_y_pos(vpos,1) = jm; + matched_moments_{jm,4}=['E(',M_.endo_names{matched_moments_{jm,1}},')']; + matched_moments_{jm,5}=['$E(',M_.endo_names_tex{matched_moments_{jm,1}},')$']; + elseif sum(matched_moments_{jm,3}) == 2 + % Second-order product moment + idx1 = (oo_.dr.obs_var == vars(1)); + idx2 = (oo_.dr.obs_var == vars(2)); + lag1 = matched_moments_{jm,2}(1); + lag2 = matched_moments_{jm,2}(2); + if lag1==0 && lag2==0 % contemporaneous covariance matrix + options_mom_.index.E_yy(idx1,idx2) = true; + options_mom_.index.E_yy(idx2,idx1) = true; + options_mom_.index.E_yy_pos(idx1,idx2) = jm; + options_mom_.index.E_yy_pos(idx2,idx1) = jm; + matched_moments_{jm,4}=['E(',M_.endo_names{matched_moments_{jm,1}(1)},',',M_.endo_names{matched_moments_{jm,1}(2)},')']; + matched_moments_{jm,5}=['$E({',M_.endo_names_tex{matched_moments_{jm,1}(1)},'}_t,{',M_.endo_names_tex{matched_moments_{jm,1}(1)},'}_t)$']; + elseif lag1==0 && lag2 < 0 + options_mom_.index.E_yyt(idx1,idx2,-lag2) = true; + options_mom_.index.E_yyt_pos(idx1,idx2,-lag2) = jm; + matched_moments_{jm,4}=['E(',M_.endo_names{matched_moments_{jm,1}(1)},',',M_.endo_names{matched_moments_{jm,1}(2)},'(',num2str(lag2),'))']; + matched_moments_{jm,5}=['$E({',M_.endo_names_tex{matched_moments_{jm,1}(1)},'}_t\times{',M_.endo_names_tex{matched_moments_{jm,1}(1)},'_{t',num2str(lag2) ,'})$']; + end + end +end + + +% @wmutschl: add check for duplicate moments by using the cellfun and unique functions +%Remove duplicate elements +UniqueMomIdx = [nonzeros(options_mom_.index.E_y_pos); nonzeros(triu(options_mom_.index.E_yy_pos)); nonzeros(options_mom_.index.E_yyt_pos)]; +DuplicateMoms = setdiff(1:size(matched_moments_,1),UniqueMomIdx); +if ~isempty(DuplicateMoms) + fprintf('Found and removed duplicate declared moments in ''matched_moments'' block in rows: %s.\n',num2str(DuplicateMoms)) +end +%reorder matched_moments_ to be compatible with options_mom_.index +matched_moments_ = matched_moments_(UniqueMomIdx,:); +if strcmp(options_mom_.mom.mom_method,'SMM') + options_mom_=rmfield(options_mom_,'index'); +end + +% Check if both prefilter and first moments were specified +options_mom_.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,matched_moments_(:,3)))'; +if options_mom_.prefilter && ~isempty(options_mom_.first_moment_indicator) + fprintf('Centered moments requested (prefilter option is set); therefore, ignore declared first moments in ''matched_moments'' block in rows: %u.\n',options_mom_.first_moment_indicator'); + matched_moments_(options_mom_.first_moment_indicator,:)=[]; %remove first moments entries + options_mom_.first_moment_indicator = []; +end +options_mom_.mom_nbr = size(matched_moments_,1); + +% Get maximum lag number for autocovariances/autocorrelations +options_mom_.ar = max(cellfun(@max,matched_moments_(:,2))) - min(cellfun(@min,matched_moments_(:,2))); + +% ------------------------------------------------------------------------- +% Step 3: Checks and transformations for estimated parameters, priors, and bounds +% ------------------------------------------------------------------------- + +% Set priors and bounds over the estimated parameters +[xparam0, estim_params_, bayestopt_, lb, ub, M_] = set_prior(estim_params_, M_, options_mom_); + +% 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.') +end + +% Check if enough moments for estimation +if options_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 +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, not ML + 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') + warning('off','MATLAB:singularMatrix'); + end +end + +% 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; +else + estim_params_.full_calibration_detected=0; +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) + else + [xparam0,estim_params_]=do_parameter_initialization(estim_params_,xparam1_calib,xparam0); %get explicitly initialized parameters that have precedence over calibrated values + end +end + +% Check whether on 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) +end + +% Set and check parameter bounds +if ~isempty(bayestopt_laplace) && any(bayestopt_laplace.pshape > 0) + % 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') + end + % Set prior bounds + Bounds = prior_bounds(bayestopt_laplace, 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. + Bounds.lb = lb; + Bounds.ub = ub; +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 + +% 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) + 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) +end + +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 MoM info structure for penalized minimization +oo_.prior.mean = bayestopt_laplace.p1; +oo_.prior.variance = diag(bayestopt_laplace.p2.^2); + +% 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_); + +% ------------------------------------------------------------------------- +% Step 4: Checks and transformations for data +% ------------------------------------------------------------------------- + +% 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.') +end + +% 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; +end + +% missing observations are not supported yet +if any(any(isnan(dataset_.data))) + error('method_of_moments: missing observations are not supported') +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 + +% Get data moments for the method of moments +[oo_.mom.data_moments, oo_.mom.m_data] = method_of_moments_data_moments(dataset_.data, oo_, matched_moments_, options_mom_); + +% Get shock series for SMM and set variance correction factor +if strcmp(options_mom_.mom.mom_method,'SMM') + options_mom_.long = round(options_mom_.mom.simulation_multiple*options_mom_.nobs); + options_mom_.variance_correction_factor = (1+1/options_mom_.mom.simulation_multiple); + % draw shocks for SMM + smmstream = RandStream('mt19937ar','Seed',options_mom_.mom.seed); + temp_shocks = randn(smmstream,options_mom_.long+options_mom_.mom.burnin,M_.exo_nbr); + temp_shocks_ME = randn(smmstream,options_mom_.long,length(M_.H)); + if options_mom_.mom.bounded_shock_support == 1 + temp_shocks(temp_shocks>2) = 2; + temp_shocks(temp_shocks<-2) = -2; + temp_shocks_ME(temp_shocks_ME<-2) = -2; + temp_shocks_ME(temp_shocks_ME<-2) = -2; + end + options_mom_.shock_series = temp_shocks; + options_mom_.ME_shock_series = temp_shocks_ME; +end + +% ------------------------------------------------------------------------- +% Step 5: checks for steady state at initial parameters +% ------------------------------------------------------------------------- + +% setting steadystate_check_flag option +if options_mom_.steadystate.nocheck + steadystate_check_flag = 0; +else + steadystate_check_flag = 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,M_,options_mom_,oo_,steadystate_check_flag); +if info(1) + fprintf('\nmethod_of_moments: The 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,Model_par_varied,options_mom_,oo_,1); + + 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.') + end +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 = 1; +else + options_mom_.noconstant = 0; +end + +% ------------------------------------------------------------------------- +% Step 6: checks for objective function at initial parameters +% ------------------------------------------------------------------------- +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_nbr); + tic_id = tic; + [fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_); + elapsed_time = toc(tic_id); + if isnan(fval) + error('method_of_moments: The initial value of the objective function is NaN') + elseif imag(fval) + error('method_of_moments: The initial value of the objective function 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('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 + 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') + fprintf('If this is not a problem with the setting of options (check the error message below),\n') + fprintf('you should try using the calibrated version of the model as starting values. To do\n') + fprintf('this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n') + fprintf('command (and after the estimated_params-block so that it does not get overwritten):\n'); + skipline(2); + end + rethrow(last_error); +end + +if options_mom_.mode_compute == 0 %We only report value of moments distance at initial value of the parameters + fprintf('No minimization of moments distance due to ''mode_compute=0''\n') + return +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 priors'); +end +if options_mom_.mode_compute == 1; fprintf('\n - optimizer (mode_compute=1): fmincon'); +elseif options_mom_.mode_compute == 2; fprintf('\n - optimizer (mode_compute=2): continuous simulated annealing'); +elseif options_mom_.mode_compute == 3; fprintf('\n - optimizer (mode_compute=3): fminunc'); +elseif options_mom_.mode_compute == 4; fprintf('\n - optimizer (mode_compute=4): csminwel'); +elseif options_mom_.mode_compute == 5; fprintf('\n - optimizer (mode_compute=5): newrat'); +elseif options_mom_.mode_compute == 6; fprintf('\n - optimizer (mode_compute=6): gmhmaxlik'); +elseif options_mom_.mode_compute == 7; fprintf('\n - optimizer (mode_compute=7): fminsearch'); +elseif options_mom_.mode_compute == 8; fprintf('\n - optimizer (mode_compute=8): Dynare Nelder-Mead simplex'); +elseif options_mom_.mode_compute == 9; fprintf('\n - optimizer (mode_compute=9): CMA-ES'); +elseif options_mom_.mode_compute == 10; fprintf('\n - optimizer (mode_compute=10): simpsa'); +elseif options_mom_.mode_compute == 11; fprintf('\n - optimizer (mode_compute=11): online_auxiliary_filter'); +elseif options_mom_.mode_compute == 12; fprintf('\n - optimizer (mode_compute=12): particleswarm'); +elseif options_mom_.mode_compute == 101; fprintf('\n - optimizer (mode_compute=101): SolveOpt'); +elseif options_mom_.mode_compute == 102; fprintf('\n - optimizer (mode_compute=102): simulannealbnd'); +elseif options_mom_.mode_compute == 13; fprintf('\n - optimizer (mode_compute=13): lsqnonlin'); +elseif ischar(minimizer_algorithm); fprintf(['\n - user-defined optimizer: ' minimizer_algorithm]); +else + error('method_of_moments: Unknown optimizer, please contact the developers ') +end +if options_mom_.silent_optimizer + fprintf(' (silent)'); +end +fprintf('\n - perturbation order: %d', options_mom_.order) +if options_mom_.order > 1 && options_mom_.pruning + fprintf(' (with pruning)') +end +fprintf('\n - number of matched moments: %d', options_mom_.mom_nbr); +fprintf('\n - number of parameters: %d\n\n', length(xparam0)); + +% ------------------------------------------------------------------------- +% Step 7b: Method of moments estimation: First-stage +% ------------------------------------------------------------------------- +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 covariance 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) + Woptflag = 0; + fprintf('Estimation stage %u\n',stage_iter); + switch lower(options_mom_.mom.weighting_matrix{stage_iter}) + case 'identity_matrix' + fprintf(' - identity weighting matrix\n'); + oo_.mom.Sw = eye(options_mom_.mom_nbr); + case 'diagonal' + %@wmutschl: better description in fprintf + fprintf(' - diagonal weighting matrix: diagonal of Newey-West estimator with lag order %d\n', options_mom_.mom.bartlett_kernel_lag); + fprintf(' and data moments as estimate of unconditional moments\n'); + W_opt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag); + oo_.mom.Sw = chol(diag(diag(W_opt))); + case 'optimal' + %@wmutschl: better description in fprintf + fprintf(' - weighting matrix: optimal. At first-stage we use diagonal of Newey-West estimator with lag order %d\n', options_mom_.mom.bartlett_kernel_lag); + fprintf(' and the data moments as initial estimate of unconditional moments\n'); + W_opt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag); + oo_.mom.Sw = chol(diag(diag(W_opt))); + Woptflag = 1; + otherwise %user specified matrix in file + fprintf(' - weighting matrix: user-specified\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 + try %check for positive definiteness + oo_.Sw = chol(weighting_matrix); + catch + error('method_of_moments: Specified weighting_matrix is not positive definite') + end + end + + optimizer_vec=[options_mom_.mode_compute,options_mom_.additional_optimizer_steps]; + + for istep= 1:length(optimizer_vec) + if optimizer_vec(istep)==13 + options_mom_.vector_output = true; + else + options_mom_.vector_output = false; + end + [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, optimizer_vec(istep), options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],... + Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_); + if options_mom_.vector_output + fval = fval'*fval; + end + fprintf('\nStage %d Iteration %d: value of minimized moment distance objective function: %12.10f.\n',stage_iter,istep,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 (FIRST-STAGE ITERATION %d) verbose',options_mom_.mom.mom_method,istep),sprintf('verbose_%s_1st_stage_iter_%d',lower(options_mom_.mom.mom_method),istep)); + end + xparam0=xparam1; + end + options_mom_.vector_output = false; + + % Update M_ and DynareResults (in particular oo_.mom.model_moments) + M_ = set_all_parameters(xparam1,estim_params_,M_); + [fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_); + + % Compute Standard errors + SE = method_of_moments_standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_, Woptflag); + + % Store first-stage 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)); + + +end + +%get optimal weighting matrix for J test, if necessary +if ~Woptflag + W_opt = method_of_moments_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_, matched_moments_, M_, options_mom_); +end + +% Compute J statistic +if strcmp(options_mom_.mom.mom_method,'SMM') + Variance_correction_factor = options_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) + +title = ['Data moments and model moments (',options_mom_.mom.mom_method,')']; +headers = {'Moment','Data','Model','% dev. target'}; +labels= matched_moments_(:,4); +data_mat=[oo_.mom.data_moments oo_.mom.model_moments 100*abs((oo_.mom.model_moments-oo_.mom.data_moments)./oo_.mom.data_moments)]; +dyntable(options_mom_, title, headers, labels, data_mat, cellofchararraymaxlength(labels)+2, 10, 7); +if options_mom_.TeX + lh = cellofchararraymaxlength(labels)+2; + labels_TeX = matched_moments_(:,5); + dyn_latex_table(M_, options_mom_, title, 'sim_corr_matrix', headers, labels_TeX, data_mat, lh, 10, 7); +end + +if options_mom_.mode_check.status + method_of_moments_mode_check(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,... + Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_) +end + +fprintf('\n==== Method of Moments Estimation (%s) Completed ====\n\n',options_mom_.mom.mom_method) + +% ------------------------------------------------------------------------- +% Step 8: Clean up +% ------------------------------------------------------------------------- +% restore warnings +warning('on','MATLAB:singularMatrix'); diff --git a/matlab/method_of_moments_datamoments.m b/matlab/method_of_moments/method_of_moments_data_moments.m similarity index 64% rename from matlab/method_of_moments_datamoments.m rename to matlab/method_of_moments/method_of_moments_data_moments.m index 4fd48c99d..c447ea1a9 100644 --- a/matlab/method_of_moments_datamoments.m +++ b/matlab/method_of_moments/method_of_moments_data_moments.m @@ -1,12 +1,12 @@ -function [dataMoments, m_data] = method_of_moments_datamoments(data, DynareResults, MatchedMoments, OptionsMoM) -% [dataMoments, m_data] = method_of_moments_datamoments(data, DynareResults, MatchedMoments, OptionsMoM) +function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_) +% [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_) % This function computes the user-selected empirical moments from data % ========================================================================= % INPUTS % o data [T x varobs_nbr] data set -% o DynareResults: [structure] storage for results (oo_) -% o MatchedMoments: [structure] information about selected moments to match in estimation (matched_moments_) -% o OptionsMoM: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_) +% o oo_: [structure] storage for results +% o matched_moments_: [structure] information about selected moments to match in estimation +% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_) % ------------------------------------------------------------------------- % OUTPUTS % o dataMoments [numMom x 1] mean of selected empirical moments @@ -40,28 +40,26 @@ function [dataMoments, m_data] = method_of_moments_datamoments(data, DynareResul % Initialization T = size(data,1); % Number of observations (T) and number of observables (ny) -mom_nbr = OptionsMoM.mom_nbr; -dataMoments = nan(mom_nbr,1); -m_data = nan(T,mom_nbr); -% Product moment for each time period, i.e. each row t contains yt1(l1)^p1*yt2(l2)^p2*... +dataMoments = NaN(options_mom_.mom_nbr,1); +m_data = NaN(T,options_mom_.mom_nbr); +% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*... % note that here we already are able to treat leads and lags and any power product moments -for jm = 1:mom_nbr - vars = DynareResults.dr.inv_order_var(MatchedMoments{jm,1})'; - leadlags = MatchedMoments{jm,2}; % note that lags are negative numbers and leads are positive numbers - powers = MatchedMoments{jm,3}; +for jm = 1:options_mom_.mom_nbr + vars = oo_.dr.inv_order_var(matched_moments_{jm,1})'; + 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 = DynareResults.dr.obs_var == vars(jv); - y = nan(T,1); - 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); + jvar = (oo_.dr.obs_var == vars(jv)); + y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via nanmean 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 m_data_tmp = y; else m_data_tmp = m_data_tmp.*y; end end - dataMoments(jm,1) = sum(m_data_tmp,'omitnan')/(T-sum(abs(leadlags))); - % We replace nan (due to leads and lags) with the corresponding mean - % @wmutschl: this should also work for missing values, right? + % We replace NaN (due to leads and lags and missing values) with the corresponding mean + dataMoments(jm,1) = mean(m_data_tmp,'omitnan'); m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1); m_data(:,jm) = m_data_tmp; end diff --git a/matlab/method_of_moments/method_of_moments_mode_check.m b/matlab/method_of_moments/method_of_moments_mode_check.m new file mode 100644 index 000000000..1df8af84f --- /dev/null +++ b/matlab/method_of_moments/method_of_moments_mode_check.m @@ -0,0 +1,185 @@ +function method_of_moments_mode_check(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin) +% Checks the estimated ML mode or Posterior mode. + + +% Copyright (C) 2020 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 . + +TeX = options_.TeX; +if ~isempty(SE_vec) + [ s_min, k ] = min(SE_vec); +end + +fval = feval(fun,xparam,varargin{:}); + +if ~isempty(SE_vec) + skipline() + disp('MODE CHECK') + skipline() + fprintf('Fval obtained by the minimization routine: %f', fval); + skipline() + if s_min(1-ll)*xparam(kk)) && (xparam(kk)+(xparam(kk)-l1). +% ------------------------------------------------------------------------- +% Author(s): +% o Willi Mutschler (willi@mutschler.eu) +% o Johannes Pfeifer (jpfeifer@uni-koeln.de) +% ========================================================================= +% To Do: check penalized estimation for different optimizers, what is special about mode_compute=1 [@wmutschl] + +%------------------------------------------------------------------------------ +% 0. Initialization of the returned variables and others... +%------------------------------------------------------------------------------ + +junk1 = []; +junk2 = []; + +%-------------------------------------------------------------------------- +% 1. Get the structural parameters & define penalties +%-------------------------------------------------------------------------- + +[fval,info,exit_flag,M_]=check_bounds_and_definiteness_estimation(xparam1, M_, options_mom_, estim_params_, Bounds); +if info(1) + if options_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 +%-------------------------------------------------------------------------- + +% Compute linear approximation around the deterministic steady state +[dr, info, M_, options_mom_, 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 + fval = Inf; + info(4) = info(2); + exit_flag = 0; + if options_mom_.vector_output == 1 % lsqnonlin requires vector output + fval = ones(size(oo_.mom.data_moments,1),1)*options_mom_.huge_number; + end + return + else + fval = Inf; + info(4) = 0.1; + exit_flag = 0; + if options_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 + +if strcmp(options_mom_.mom.mom_method,'GMM') + %-------------------------------------------------------------------------- + % 3. Set up pruned state-space system and compute model moments + %-------------------------------------------------------------------------- + pruned_state_space = pruned_state_space_system(M_, options_mom_, dr, oo_.dr.obs_var, options_mom_.ar, 0, 0); + + oo_.mom.model_moments = NaN(options_mom_.mom_nbr,1); + offset = 0; + % First moments + if ~options_mom_.prefilter && isfield(options_mom_.index,'E_y') && nnz(options_mom_.index.E_y) > 0 + E_y = pruned_state_space.E_y; + E_y_nbr = nnz(options_mom_.index.E_y); + oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.index.E_y); + offset = offset + E_y_nbr; + end + % Second moments + % Contemporaneous covariance + if isfield(options_mom_.index,'E_yy') && nnz(options_mom_.index.E_yy) > 0 + if options_mom_.prefilter + E_yy = pruned_state_space.Var_y; + else + E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y'; + end + E_yy_nbr = nnz(triu(options_mom_.index.E_yy)); + oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(triu(options_mom_.index.E_yy)); + offset = offset + E_yy_nbr; + end + % Lead/lags covariance + if isfield(options_mom_.index,'E_yyt') && nnz(options_mom_.index.E_yyt) > 0 + if options_mom_.prefilter + E_yyt = pruned_state_space.Var_yi; + else + E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]); + end + E_yyt_nbr = nnz(options_mom_.index.E_yyt); + oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.index.E_yyt); + end + +elseif strcmp(options_mom_.mom.mom_method,'SMM') + %------------------------------------------------------------------------------ + % 3. Compute Moments of the model solution for normal innovations + %------------------------------------------------------------------------------ + + % 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 + chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var)); + scaled_shock_series = zeros(size(options_mom_.shock_series)); %initialize + scaled_shock_series(:,i_exo_var) = options_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_.mode_compute==13 + fval = Inf(size(oo_.mom.Sw,1),1); + else + fval = Inf; + end + info(1)=180; + info(4) = 0.1; + exit_flag = 0; + if options_mom_.mode_compute == 13 + fval = ones(size(oo_.mom.dataMoments,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_.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_.ME_shock_series)); %initialize + shock_mat(:,i_ME)=options_mom_.ME_shock_series(:,i_exo_var)*chol_S; + y_sim = y_sim+shock_mat; + end + + % Remove mean if centered moments + if options_mom_.prefilter + y_sim = bsxfun(@minus, y_sim, mean(y_sim,1)); + end + oo_.mom.model_moments = method_of_moments_data_moments(y_sim, oo_, matched_moments_, options_mom_); + +end + +%-------------------------------------------------------------------------- +% 4. 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 + + +end%main function end + diff --git a/matlab/method_of_moments_optimal_weighting_matrix.m b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m similarity index 62% rename from matlab/method_of_moments_optimal_weighting_matrix.m rename to matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m index 8860a7f6b..e234fbf5b 100644 --- a/matlab/method_of_moments_optimal_weighting_matrix.m +++ b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m @@ -1,17 +1,17 @@ -function Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag) -% Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag) +function [W_opt, normalization_factor]= method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag) +% W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag) % ------------------------------------------------------------------------- -% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag qlag +% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag % Adapted from replication codes of -% 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 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. % ========================================================================= % INPUTS % o m_data [T x numMom] selected empirical or theoretical moments at each point in time % o moments [numMom x 1] mean of selected empirical or theoretical moments -% o qlag [integer] Bartlett kernel maximum lag order +% o q_lag [integer] Bartlett kernel maximum lag order % ------------------------------------------------------------------------- % OUTPUTS -% o Wopt [numMom x numMom] optimal weighting matrix +% o W_opt [numMom x numMom] optimal weighting matrix % ------------------------------------------------------------------------- % This function is called by % o method_of_moments.m @@ -42,45 +42,45 @@ function Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag % ========================================================================= % Initialize -[T,numMom] = size(m_data); %note that in m_data nan values (due to leads or lags in matchedmoments) are removed so T is the effective sample size +[T,num_Mom] = size(m_data); %note that in m_data NaN values (due to leads or lags in matched_moments and missing data) were replaced by the mean -% center around moments (could be either datamoments or modelmoments) -hFunc = m_data - repmat(moments',T,1); +% center around moments (could be either data_moments or model_moments) +h_Func = m_data - repmat(moments',T,1); % The required correlation matrices -GAMA_array = zeros(numMom,numMom,qLag); -GAMA0 = CorrMatrix(hFunc,T,numMom,0); -if qLag > 0 - for ii=1:qLag - GAMA_array(:,:,ii) = CorrMatrix(hFunc,T,numMom,ii); +GAMA_array = zeros(num_Mom,num_Mom,q_lag); +GAMA0 = Corr_Matrix(h_Func,T,num_Mom,0); +if q_lag > 0 + for ii=1:q_lag + GAMA_array(:,:,ii) = Corr_Matrix(h_Func,T,num_Mom,ii); end end % The estimate of S S = GAMA0; -if qLag > 0 - for ii=1:qLag - S = S + (1-ii/(qLag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)'); +if q_lag > 0 + for ii=1:q_lag + S = S + (1-ii/(q_lag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)'); end end % The estimate of W -Wopt = S\eye(size(S,1)); +W_opt = S\eye(size(S,1)); % Check positive definite W try - chol(Wopt); + chol(W_opt); catch err - error('method_of_moments: The optimal weighting matrix is not positive definite. Check whether your model implies stochastic singularity\n') + error('method_of_moments: The optimal weighting matrix is not positive definite. Check whether your model implies stochastic singularity.\n') end end % The correlation matrix -function GAMAcorr = CorrMatrix(hFunc,T,numMom,v) - GAMAcorr = zeros(numMom,numMom); +function GAMA_corr = Corr_Matrix(h_Func,T,num_Mom,v) + GAMA_corr = zeros(num_Mom,num_Mom); for t = 1+v:T - GAMAcorr = GAMAcorr + hFunc(t-v,:)'*hFunc(t,:); + GAMA_corr = GAMA_corr + h_Func(t-v,:)'*h_Func(t,:); end - GAMAcorr = GAMAcorr/T; + GAMA_corr = GAMA_corr/T; end \ No newline at end of file diff --git a/matlab/method_of_moments/method_of_moments_standard_errors.m b/matlab/method_of_moments/method_of_moments_standard_errors.m new file mode 100644 index 000000000..54553d737 --- /dev/null +++ b/matlab/method_of_moments/method_of_moments_standard_errors.m @@ -0,0 +1,105 @@ +function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_, Wopt_flag) +% [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_, Wopt_flag) +% ------------------------------------------------------------------------- +% This function computes standard errors to the method of moments estimates +% Adapted from replication codes of +% 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. +% ========================================================================= +% INPUTS +% o xparam: value of estimated parameters as returned by set_prior() +% o objective_function string of objective function, either method_of_moments_GMM.m or method_of_moments_SMM.m +% o Bounds: structure containing parameter bounds +% o oo_: structure for results +% o estim_params_: structure describing the estimated_parameters +% o matched_moments_: structure containing information about selected moments to match in estimation +% 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 Wopt_flag: indicator whether the optimal weighting is actually used +% ------------------------------------------------------------------------- +% OUTPUTS +% o SE_values [nparam x 1] vector of standard errors +% o Asympt_Var [nparam x nparam] asymptotic covariance matrix +% ------------------------------------------------------------------------- +% This function is called by +% o method_of_moments.m +% ------------------------------------------------------------------------- +% This function calls: +% o get_the_name +% o get_error_message +% o GMM_objective_function +% o SMM_objective_function.m +% o method_of_moments_optimal_weighting_matrix +% ========================================================================= +% Copyright (C) 2020 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 . +% ------------------------------------------------------------------------- +% Author(s): +% o Willi Mutschler (willi@mutschler.eu) +% o Johannes Pfeifer (jpfeifer@uni-koeln.de) +% ========================================================================= + +% Some dimensions +num_mom = size(oo_.mom.model_moments,1); +dim_params = size(xparam,1); +D = zeros(num_mom,dim_params); +eps_value = options_mom_.mom.se_tolx; + +for i=1:dim_params + %Positive step + xparam_eps_p = xparam; + xparam_eps_p(i,1) = xparam_eps_p(i) + eps_value; + [~, info_p, ~, ~,~, oo__p] = feval(objective_function, xparam_eps_p, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_); + + % Negative step + xparam_eps_m = xparam; + xparam_eps_m(i,1) = xparam_eps_m(i) - eps_value; + [~, info_m, ~, ~,~, oo__m] = feval(objective_function, xparam_eps_m, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_); + + % The Jacobian: + if nnz(info_p)==0 && nnz(info_m)==0 + D(:,i) = (oo__p.mom.model_moments - oo__m.mom.model_moments)/(2*eps_value); + else + problpar = get_the_name(i,options_mom_.TeX, M_, estim_params_, options_mom_); + message_p = get_error_message(info_p, options_mom_); + message_m = get_error_message(info_m, options_mom_); + + warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m) + Asympt_Var = NaN(length(xparam),length(xparam)); + SE_values = NaN(length(xparam),1); + return + end +end + +T = options_mom_.nobs; %Number of observations +if isfield(options_mom_,'variance_correction_factor') + T = T*options_mom_.variance_correction_factor; +end + +if Wopt_flag + % We have the optimal weighting matrix + WW = oo_.mom.Sw'*oo_.mom.Sw; + Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params)); +else + % We do not have the optimal weighting matrix yet + WWused = oo_.mom.Sw'*oo_.mom.Sw; + WWopt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag); + S = WWopt\eye(size(WWopt,1)); + AA = (D'*WWused*D)\eye(dim_params); + Asympt_Var = 1/T*AA*D'*WWused*S*WWused*D*AA; +end + +SE_values = sqrt(diag(Asympt_Var)); \ No newline at end of file diff --git a/matlab/method_of_moments_GMM.m b/matlab/method_of_moments_GMM.m deleted file mode 100644 index d70756a63..000000000 --- a/matlab/method_of_moments_GMM.m +++ /dev/null @@ -1,224 +0,0 @@ -function [fval, info, exit_flag, DynareResults, Model, OptionsMoM] = method_of_moments_GMM(xparam1, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM) -% [fval, info, exit_flag, DynareResults, Model, OptionsMoM] = method_of_moments_GMM(xparam1, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM) -% ------------------------------------------------------------------------- -% This function evaluates the objective function for GMM estimation -% ========================================================================= -% INPUTS -% o xparam1: current value of estimated parameters as returned by set_prior() -% o Bounds: structure containing parameter bounds -% o DynareResults: structure for results (oo_) -% o EstimatedParameters: structure describing the estimated_parameters (estim_params_) -% o MatchedMoments: structure containing information about selected moments to match in estimation (matched_moments_) -% o Model structure describing the Model -% o OptionsMoM: structure information about all settings (specified by the user, preprocessor, and taken from global options_) -% ------------------------------------------------------------------------- -% 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 no error, 1 of error -% o DynareResults: structure containing the results with the following updated fields: -% - mom.modelMoments [numMom x 1] vector with model moments -% - mom.Q value of the quadratic form of the moment difference -% o Model: Matlab's structure describing the Model -% ------------------------------------------------------------------------- -% This function is called by -% o driver.m -% ------------------------------------------------------------------------- -% This function calls -% o ispd -% o pruned_state_space_system -% o resol -% o set_all_parameters -% ========================================================================= -% Copyright (C) 2020 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 . -% ------------------------------------------------------------------------- -% Author(s): -% o Willi Mutschler (willi@mutschler.eu) -% o Johannes Pfeifer (jpfeifer@uni-koeln.de) -% ========================================================================= -% To Do: check penalized estimation for different optimizers, what is special about mode_compute=1 [@wmutschl] - -%-------------------------------------------------------------------------- -% 0. Initialization of the returned variables and others... -%-------------------------------------------------------------------------- -exit_flag = 1; -info = zeros(4,1); -xparam1 = xparam1(:); % Ensure that xparam1 is a column vector - -%-------------------------------------------------------------------------- -% 1. Get the structural parameters & define penalties -%-------------------------------------------------------------------------- - -% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the parameters. -if any(xparam1Bounds.ub) - if ~isequal(OptionsMoM.mode_compute,1) - k = find(xparam1>Bounds.ub); - fval = Inf; - exit_flag = 0; - info(1) = 42; - info(4)= sum((xparam1(k)-Bounds.ub(k)).^2); - return - elseif OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - return - end -end - -% Set all parameters -Model = set_all_parameters(xparam1,EstimatedParameters,Model); - -% Test if Q is positive definite. -if ~issquare(Model.Sigma_e) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances') - [Q_is_positive_definite, penalty] = ispd(Model.Sigma_e(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness)); - if ~Q_is_positive_definite - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - else - fval = Inf; - exit_flag = 0; - info(1) = 43; - info(4) = penalty; - end - return - end - if isfield(EstimatedParameters,'calibrated_covariances') - correct_flag=check_consistency_covariances(Model.Sigma_e); - if ~correct_flag - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - else - penalty = sum(Model.Sigma_e(EstimatedParameters.calibrated_covariances.position).^2); - fval = Inf; - exit_flag = 0; - info(1) = 71; - info(4) = penalty; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - end - return - end - end -end - -%-------------------------------------------------------------------------- -% 2. call resol to compute steady state and model solution -%-------------------------------------------------------------------------- - -% Compute linear approximation around the deterministic steady state -[dr, info, Model, OptionsMoM, DynareResults] = resol(0, Model, OptionsMoM, DynareResults); - -% Return, with endogenous penalty when possible, if resol issues an error code -if info(1) - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - return - else - 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 - fval = Inf; - info(4) = info(2); - exit_flag = 0; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - return - else - fval = Inf; - info(4) = 0.1; - exit_flag = 0; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - return - end - end -end - -%-------------------------------------------------------------------------- -% 3. Set up pruned state-space system and compute model moments -%-------------------------------------------------------------------------- -pruned_state_space = pruned_state_space_system(Model, OptionsMoM, dr, DynareResults.dr.obs_var, OptionsMoM.ar, 0, 0); - -DynareResults.mom.modelMoments = nan(OptionsMoM.mom_nbr,1); -offset = 0; -% First moments -if isfield(OptionsMoM.index,'E_y') && nnz(OptionsMoM.index.E_y) > 0 && ~OptionsMoM.prefilter - E_y = pruned_state_space.E_y; - E_y_nbr = nnz(OptionsMoM.index.E_y); - DynareResults.mom.modelMoments(offset+1:E_y_nbr,1) = E_y(OptionsMoM.index.E_y); - offset = offset + E_y_nbr; -end -% Second moments -if isfield(OptionsMoM.index,'E_yy') && nnz(OptionsMoM.index.E_yy) > 0 - if OptionsMoM.prefilter - E_yy = pruned_state_space.Var_y; - else - E_yy = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y'; - end - E_yy_nbr = nnz(triu(OptionsMoM.index.E_yy)); - DynareResults.mom.modelMoments(offset+(1:E_yy_nbr),1) = E_yy(triu(OptionsMoM.index.E_yy)); - offset = offset + E_yy_nbr; -end - -if isfield(OptionsMoM.index,'E_yyt') && nnz(OptionsMoM.index.E_yyt) > 0 - if OptionsMoM.prefilter - E_yyt = pruned_state_space.Var_yi; - else - E_yyt = pruned_state_space.Var_yi + repmat(pruned_state_space.E_y*pruned_state_space.E_y',[1 1 size(pruned_state_space.Var_yi,3)]); - end - E_yyt_nbr = nnz(OptionsMoM.index.E_yyt); - DynareResults.mom.modelMoments(offset+(1:E_yyt_nbr),1) = E_yyt(OptionsMoM.index.E_yyt); -end - - -%-------------------------------------------------------------------------- -% 4. Compute quadratic target function -%-------------------------------------------------------------------------- -moments_difference = DynareResults.mom.dataMoments - DynareResults.mom.modelMoments; -residuals = DynareResults.mom.Sw*moments_difference; -DynareResults.mom.Q = residuals'*residuals; -if OptionsMoM.mode_compute == 13 % lsqnonlin - fval = residuals; -else - fval = DynareResults.mom.Q; - if OptionsMoM.penalized_estimator - fval=fval+(xparam1-DynareResults.prior.p1)'/diag(DynareResults.prior.p2)*(xparam1-DynareResults.prior.p1); - end -end - - -end%main function end - diff --git a/matlab/method_of_moments_SMM.m b/matlab/method_of_moments_SMM.m deleted file mode 100644 index ea1f09155..000000000 --- a/matlab/method_of_moments_SMM.m +++ /dev/null @@ -1,223 +0,0 @@ -function [fval, info, exit_flag, DynareResults, Model, OptionsMoM] = method_of_moments_SMM(xparam1, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM) -% [fval, info, exit_flag, DynareResults, Model, OptionsMoM] = method_of_moments_SMM(xparam1, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM) -% ------------------------------------------------------------------------- -% This function evaluates the objective function for SMM estimation -% ========================================================================= -% INPUTS -% o xparam1: current value of estimated parameters as returned by set_prior() -% o Bounds: structure containing parameter bounds -% o DynareResults: structure for results (oo_) -% o EstimatedParameters: structure describing the estimated_parameters (estim_params_) -% o MatchedMoments: structure containing information about selected moments to match in estimation (matched_moments_) -% o Model structure describing the Model -% o OptionsMoM: structure information about all settings (specified by the user, preprocessor, and taken from global options_) -% ------------------------------------------------------------------------- -% 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 no error, 1 of error -% o DynareResults: structure containing the results with the following updated fields: -% - mom.modelMoments [numMom x 1] vector with model moments -% - mom.Q value of the quadratic form of the moment difference -% o Model: Matlab's structure describing the Model -% ------------------------------------------------------------------------- -% This function is called by -% o driver.m -% ------------------------------------------------------------------------- -% This function calls -% o bsxfun -% o ispd -% o method_of_moments_datamoments -% o resol -% o set_all_parameters -% o simult_ -% ========================================================================= -% Copyright (C) 2020 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 . -% ------------------------------------------------------------------------- -% Author(s): -% o Willi Mutschler (willi@mutschler.eu) -% o Johannes Pfeifer (jpfeifer@uni-koeln.de) -% ========================================================================= -% To Do: check penalized estimation for different optimizers, what is special about mode_compute=1 [@wmutschl] - -%------------------------------------------------------------------------------ -% 0. Initialization of the returned variables and others... -%------------------------------------------------------------------------------ -exit_flag = 1; -info = zeros(4,1); -xparam1 = xparam1(:); % Ensure that xparam1 is a column vector - -%------------------------------------------------------------------------------ -% 1. Get the structural parameters & define penalties -%------------------------------------------------------------------------------ - -% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the parameters. -if any(xparam1Bounds.ub) - if ~isequal(OptionsMoM.mode_compute,1) - k = find(xparam1>Bounds.ub); - fval = Inf; - exit_flag = 0; - info(1) = 42; - info(4)= sum((xparam1(k)-Bounds.ub(k)).^2); - return - elseif OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - return - end -end - -% Set all parameters -Model = set_all_parameters(xparam1,EstimatedParameters,Model); - -% Test if Q is positive definite. -if ~issquare(Model.Sigma_e) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances') - [Q_is_positive_definite, penalty] = ispd(Model.Sigma_e(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness)); - if ~Q_is_positive_definite - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - else - fval = Inf; - exit_flag = 0; - info(1) = 43; - info(4) = penalty; - end - return - end - if isfield(EstimatedParameters,'calibrated_covariances') - correct_flag=check_consistency_covariances(Model.Sigma_e); - if ~correct_flag - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - else - penalty = sum(Model.Sigma_e(EstimatedParameters.calibrated_covariances.position).^2); - fval = Inf; - exit_flag = 0; - info(1) = 71; - info(4) = penalty; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - end - return - end - end -end - -%------------------------------------------------------------------------------ -% 2. call resol to compute steady state and model solution -%------------------------------------------------------------------------------ - -% Compute linear approximation around the deterministic steady state -[dr, info, Model, OptionsMoM, DynareResults] = resol(0, Model, OptionsMoM, DynareResults); - -% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol). -if info(1) - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - return - else - 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 - fval = Inf; - info(4) = info(2); - exit_flag = 0; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - return - else - fval = Inf; - info(4) = 0.1; - exit_flag = 0; - if OptionsMoM.mode_compute == 13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - end - return - end - end -end - - -%------------------------------------------------------------------------------ -% 3. Compute Moments of the model solution for normal innovations -%------------------------------------------------------------------------------ - -% create shock series with correct covariance matrix from iid standard normal shocks -i_exo_var = setdiff(1:Model.exo_nbr, find(diag(Model.Sigma_e) == 0 )); %find singular entries in covariance -chol_S = chol(Model.Sigma_e(i_exo_var,i_exo_var)); -scaled_shock_series = zeros(size(OptionsMoM.shock_series)); %initialize -scaled_shock_series(:,i_exo_var) = OptionsMoM.shock_series(:,i_exo_var)*chol_S; %set non-zero entries - -% simulate series -y_sim = simult_(Model, OptionsMoM, dr.ys, dr, scaled_shock_series, OptionsMoM.order); -% provide meaningful penalty if data is nan or inf -if any(any(isnan(y_sim))) || any(any(isinf(y_sim))) - if OptionsMoM.mode_compute==13 - fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number; - else - fval = Inf; - end - info(1)=180; - info(4) = 0.1; - exit_flag = 0; - return -end - -% Remove burning and focus on observables (note that y_sim is in declaration order) -y_sim = y_sim(DynareResults.dr.order_var(DynareResults.dr.obs_var) , end-OptionsMoM.long:end)'; - -% Remove mean if centered moments -if OptionsMoM.prefilter - y_sim = bsxfun(@minus, y_sim, mean(y_sim,1)); -end -DynareResults.mom.modelMoments = method_of_moments_datamoments(y_sim, DynareResults, MatchedMoments, OptionsMoM); - - -%------------------------------------------------------------------------------ -% 4. Compute quadratic target function -%------------------------------------------------------------------------------ -moments_difference = DynareResults.mom.dataMoments - DynareResults.mom.modelMoments; -residuals = DynareResults.mom.Sw*moments_difference; -DynareResults.mom.Q = residuals'*residuals; -if OptionsMoM.mode_compute == 13 % lsqnonlin - fval = residuals; -else - fval = DynareResults.mom.Q; - if OptionsMoM.penalized_estimator - fval=fval+(xparam1-DynareResults.prior.p1)'/diag(DynareResults.mom.p2)*(xparam1-DynareResults.mom.p1); - end -end - -end %main function end \ No newline at end of file diff --git a/matlab/method_of_moments_standard_errors.m b/matlab/method_of_moments_standard_errors.m deleted file mode 100644 index 27396446f..000000000 --- a/matlab/method_of_moments_standard_errors.m +++ /dev/null @@ -1,105 +0,0 @@ -function [SEvalues, AVar] = method_of_moments_standard_errors(xparam, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Wopt_flag) -% [SEvalues, AVar] = method_of_moments_standard_errors(xparam, objective_function, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM, Wopt_flag) -% ------------------------------------------------------------------------- -% This function computes standard errors to the method of moments estimates -% Adapted from replication codes of -% 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. -% ========================================================================= -% INPUTS -% o xparam: value of estimated parameters as returned by set_prior() -% o objective_function string of objective function, either method_of_moments_GMM.m or method_of_moments_SMM.m -% o Bounds: structure containing parameter bounds -% o DynareResults: structure for results (oo_) -% o EstimatedParameters: structure describing the estimated_parameters (estim_params_) -% o MatchedMoments: structure containing information about selected moments to match in estimation (matched_moments_) -% o Model structure describing the Model -% o OptionsMoM: structure information about all settings (specified by the user, preprocessor, and taken from global options_) -% o Wopt_flag: indicator whether the optimal weighting is actually used -% ------------------------------------------------------------------------- -% OUTPUTS -% o SEvalues [nparam x 1] vector of standard errors -% o AVar [nparam x nparam] asymptotic covariance matrix -% ------------------------------------------------------------------------- -% This function is called by -% o method_of_moments.m -% ------------------------------------------------------------------------- -% This function calls: -% o get_the_name -% o get_error_message -% o method_of_moments_GMM.m (objective function) -% o method_of_moments_SMM.m (objective function) -% o method_of_moments_optimal_weighting_matrix -% ========================================================================= -% Copyright (C) 2020 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 . -% ------------------------------------------------------------------------- -% Author(s): -% o Willi Mutschler (willi@mutschler.eu) -% o Johannes Pfeifer (jpfeifer@uni-koeln.de) -% ========================================================================= - -% Some dimensions -numMom = size(DynareResults.mom.modelMoments,1); -dimParams = size(xparam,1); -D = zeros(numMom,dimParams); -epsValue = OptionsMoM.dynatol.x; - -for i=1:dimParams - %Positive step - xparam_eps_p = xparam; - xparam_eps_p(i,1) = xparam_eps_p(i) + epsValue; - [~, info_p, exit_flag_p, DynareResults_p, ~, ~] = feval(objective_function, xparam_eps_p, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); - - % Negative step - xparam_eps_m = xparam; - xparam_eps_m(i,1) = xparam_eps_m(i) - epsValue; - [~, info_m, exit_flag_m, DynareResults_m, ~, ~] = feval(objective_function, xparam_eps_m, Bounds, DynareResults, EstimatedParameters, MatchedMoments, Model, OptionsMoM); - - % The Jacobian: - if nnz(info_p)==0 && nnz(info_m)==0 - D(:,i) = (DynareResults_p.mom.modelMoments - DynareResults_m.mom.modelMoments)/(2*epsValue); - else - problpar = get_the_name(i,OptionsMoM.TeX, Model, EstimatedParameters, OptionsMoM); - message_p = get_error_message(info_p, OptionsMoM); - message_m = get_error_message(info_m, OptionsMoM); - - warning('method_of_moments:info','Cannot compute the Jacobian for parameter %s - no standard errors available\n %s %s\nCheck your bounds and/or priors, or use a different optimizer.\n',problpar, message_p, message_m) - AVar = NaN(length(xparam),length(xparam)); - SEvalues = NaN(length(xparam),1); - return - end -end - -T = OptionsMoM.nobs; %Number of observations -if isfield(OptionsMoM,'variance_correction_factor') - T = T*OptionsMoM.variance_correction_factor; -end - -if Wopt_flag - % We have the optimal weighting matrix - WW = DynareResults.mom.Sw'*DynareResults.mom.Sw; - AVar = 1/T*((D'*WW*D)\eye(dimParams)); -else - % We do not have the optimal weighting matrix yet - WWused = DynareResults.mom.Sw'*DynareResults.mom.Sw; - WWopt = method_of_moments_optimal_weighting_matrix(DynareResults.mom.m_data, DynareResults.mom.modelMoments, OptionsMoM.bartlett_kernel_lag); - S = WWopt\eye(size(WWopt,1)); - AA = (D'*WWused*D)\eye(dimParams); - AVar = 1/T*AA*D'*WWused*S*WWused*D*AA; -end - -SEvalues = sqrt(diag(AVar)); \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index 1904e08c6..cfaf424ff 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -25,6 +25,8 @@ MODFILES = \ measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \ TeX/fs2000_corr_ME.mod \ estimation/MH_recover/fs2000_recover_tarb.mod \ + estimation/method_of_moments/RBC_MoM.mod \ + estimation/method_of_moments/RBC_MoM_SMM_ME.mod \ estimation/fs2000.mod \ gsa/ls2003a.mod \ optimizers/fs2000_8.mod \ @@ -961,6 +963,8 @@ EXTRA_DIST = \ lmmcp/sw-common-header.inc \ lmmcp/sw-common-footer.inc \ estimation/tune_mh_jscale/fs2000.inc \ + estimation/method_of_moments/RBC_MoM_common.inc \ + estimation/method_of_moments/RBC_MoM_steady_helper.m \ histval_initval_file_unit_tests.m \ histval_initval_file/my_assert.m \ histval_initval_file/ramst_data.xls \ diff --git a/tests/estimation/method_of_moments/AnScho_MoM.mod b/tests/estimation/method_of_moments/AnScho_MoM.mod index 96d2be001..7b482897c 100644 --- a/tests/estimation/method_of_moments/AnScho_MoM.mod +++ b/tests/estimation/method_of_moments/AnScho_MoM.mod @@ -25,7 +25,7 @@ @#define estimParams = 1 % Note that we set the numerical optimization tolerance levels very large to speed up the testsuite -@#define optimizer = 1 +@#define optimizer = 4 var c p R g y z INFL INT YGR; varexo e_r e_g e_z; @@ -190,12 +190,12 @@ matched_moments_ = { % Options for both GMM and SMM % , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix - , order = @{orderApp} % order of Taylor approximation in perturbation + , order = @{orderApp} % order of Taylor approximation in perturbation % , penalized_estimator % use penalized optimization - , pruning % use pruned state space system at higher-order + , pruning % use pruned state space system at higher-order % , verbose % display and store intermediate estimation results - , weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename - , mom_steps = [2 2] % vector of numbers for the iterations in the 2-step feasible method of moments + , weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + , additional_optimizer_steps = [4] % vector of numbers for the iterations in the 2-step feasible method of moments % , prefilter=0 % demean each data series by its empirical mean and use centered moments % % Options for SMM diff --git a/tests/estimation/method_of_moments/RBC_MoM.mod b/tests/estimation/method_of_moments/RBC_MoM.mod index 5e3f50d6a..7714ac5f2 100644 --- a/tests/estimation/method_of_moments/RBC_MoM.mod +++ b/tests/estimation/method_of_moments/RBC_MoM.mod @@ -1,7 +1,5 @@ -% RBC model used in replication files of -% 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. -% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu) -% ========================================================================= +% Tests SMM and GMM routines +% % Copyright (C) 2020 Dynare Team % % This file is part of Dynare. @@ -21,92 +19,72 @@ % ========================================================================= % Define testscenario -@#define orderApp = 3 +@#define orderApp = 1 @#define estimParams = 1 % Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite @#define optimizer = 13 -var k c a iv y la n rk w; -predetermined_variables k; -varexo u_a; -varobs n c iv; -parameters DELTA BETTA B ETAl ETAc THETA ALFA RHOA STDA; -DELTA = 0.025; -BETTA = 0.984; -B = 0.5; -ETAl = 1; -ETAc = 2; -THETA = 3.48; -ALFA = 0.667; -RHOA = 0.979; -STDA = 0.0072; - -model; -0 = -exp(la) +(exp(c)-B*exp(c(-1)))^(-ETAc) - BETTA*B*(exp(c(+1))-B*exp(c))^(-ETAc); -0 = -THETA*(1-exp(n))^-ETAl + exp(la)*exp(w); -0 = -exp(la) + BETTA*exp(la(+1))*(exp(rk(+1)) + (1-DELTA)); -0 = -exp(a)*(1-ALFA)*exp(k)^(-ALFA)*exp(n)^(ALFA) + exp(rk); -0 = -exp(a)*ALFA*exp(k)^(1-ALFA)*exp(n)^(ALFA-1) + exp(w); -0 = -exp(c) - exp(iv) + exp(y); -0 = -exp(y) + exp(a)*exp(k)^(1-ALFA)*exp(n)^(ALFA); -0 = -exp(k(+1)) + (1-DELTA)*exp(k) + exp(iv); -0 = -log(exp(a)) + RHOA*log(exp(a(-1))) + STDA*u_a; -end; +@#include "RBC_MoM_common.inc" shocks; -var u_a = 1; -end; +var u_a; stderr 0.0072; +end; + +varobs n c iv; + @#if estimParams == 0 estimated_params; - DELTA, 0.02; - BETTA, 0.9; - B, 0.4; + DELTA, 0.025; + BETTA, 0.98; + B, 0.45; %ETAl, 1; - ETAc, 1.5; - ALFA, 0.6; - RHOA, 0.9; - STDA, 0.01; + ETAc, 1.8; + ALFA, 0.65; + RHOA, 0.95; + stderr u_a, 0.01; %THETA, 3.48; end; @#endif @#if estimParams == 1 estimated_params; - DELTA, 0.02, 0, 1; - BETTA, 0.90, 0, 1; - B, 0.40, 0, 1; + DELTA, , 0, 1; + BETTA, , 0, 1; + B, , 0, 1; %ETAl, 1, 0, 10; - ETAc, 1.80, 0, 10; - ALFA, 0.60, 0, 1; - RHOA, 0.90, 0, 1; - STDA, 0.01, 0, 1; + ETAc, , 0, 10; + ALFA, , 0, 1; + RHOA, , 0, 1; + stderr u_a, , 0, 1; %THETA, 3.48, 0, 10; end; @#endif @#if estimParams == 2 estimated_params; - DELTA, 0.02, 0, 1, normal_pdf, 0.02, 0.1; - BETTA, 0.90, 0, 1, normal_pdf, 0.90, 0.1; - B, 0.40, 0, 1, normal_pdf, 0.40, 0.1; + DELTA, 0.025, 0, 1, normal_pdf, 0.02, 0.5; + BETTA, 0.98, 0, 1, beta_pdf, 0.90, 0.25; + B, 0.45, 0, 1, normal_pdf, 0.40, 0.5; %ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1; - ETAc, 1.80, 0, 10, normal_pdf, 1.80, 0.1; - ALFA, 0.60, 0, 1, normal_pdf, 0.60, 0.1; - RHOA, 0.90, 0, 1, normal_pdf, 0.90, 0.1; - STDA, 0.01, 0, 1, normal_pdf, 0.01, 0.1; + ETAc, 1.8, 0, 10, normal_pdf, 1.80, 0.5; + ALFA, 0.65, 0, 1, normal_pdf, 0.60, 0.5; + RHOA, 0.95, 0, 1, normal_pdf, 0.90, 0.5; + stderr u_a, 0.01, 0, 1, normal_pdf, 0.01, 0.5; %THETA, 3.48, 0, 10, normal_pdf, 0.25, 0.0.1; end; @#endif % Simulate data -stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500); +stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=500); save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} ); pause(1); +estimated_params_init(use_calibration); +end; %-------------------------------------------------------------------------- % Method of Moments Estimation @@ -128,7 +106,7 @@ matched_moments_ = { [ic ic ] [0 0], [1 1]; [ic iiv] [0 0], [1 1]; [ic in ] [0 0], [1 1]; -% [iiv ic ] [0 0], [1 1]; + [iiv ic ] [0 0], [1 1]; [iiv iiv] [0 0], [1 1]; [iiv in ] [0 0], [1 1]; % [in ic ] [0 0], [1 1]; @@ -150,12 +128,13 @@ matched_moments_ = { % Options for both GMM and SMM % , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix - , order = @{orderApp} % order of Taylor approximation in perturbation + , order = @{orderApp} % order of Taylor approximation in perturbation % , penalized_estimator % use penalized optimization - , pruning % use pruned state space system at higher-order + , pruning % use pruned state space system at higher-order % , verbose % display and store intermediate estimation results - , weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename - , mom_steps = [2 2] % vector of numbers for the iterations in the 2-step feasible method of moments + , weighting_matrix = ['optimal','optimal'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + , weighting_matrix_scaling_factor=1 + %, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute % , prefilter=0 % demean each data series by its empirical mean and use centered moments % % Options for SMM @@ -189,7 +168,7 @@ matched_moments_ = { %, optim = ('TolFun', 1e-3 % ,'TolX', 1e-5 % ) % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute - , silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between + %, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between % , tolf = 1e-5 % convergence criterion on function value for numerical differentiation % , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation % diff --git a/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod b/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod new file mode 100644 index 000000000..a5199c1c0 --- /dev/null +++ b/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod @@ -0,0 +1,189 @@ +% +% 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 . +% ========================================================================= + +% Define testscenario +@#define orderApp = 2 +@#define estimParams = 0 + +% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite +@#define optimizer = 13 + +@#include "RBC_MoM_common.inc" + +shocks; +var u_a; stderr 0.0072; +var n; stderr 0.01; +end; + +varobs n c iv; + +@#if estimParams == 0 +estimated_params; + DELTA, 0.02; + BETTA, 0.9; + B, 0.4; + %ETAl, 1; + ETAc, 1.5; + ALFA, 0.6; + RHOA, 0.9; + stderr u_a, 0.010; + %THETA, 3.48; + stderr n, 0.01; + +end; +@#endif + +@#if estimParams == 1 +estimated_params; + DELTA, 0.02, 0, 1; + BETTA, 0.90, 0, 1; + B, 0.40, 0, 1; + %ETAl, 1, 0, 10; + ETAc, 1.80, 0, 10; + ALFA, 0.60, 0, 1; + RHOA, 0.90, 0, 1; + stderr u_a, 0.01, 0, 1; + stderr n, 0.01, 0, 1; +end; +@#endif + +@#if estimParams == 2 +estimated_params; + DELTA, 0.02, 0, 1, normal_pdf, 0.02, 0.5; + BETTA, 0.90, 0, 1, beta_pdf, 0.90, 0.25; + B, 0.40, 0, 1, normal_pdf, 0.40, 0.5; + %ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1; + ETAc, 1.80, 0, 10, normal_pdf, 1.80, 0.5; + ALFA, 0.60, 0, 1, normal_pdf, 0.60, 0.5; + RHOA, 0.90, 0, 1, normal_pdf, 0.90, 0.5; + stderr u_a, 0.01, 0, 1, normal_pdf, 0.01, 0.5; + stderr n, 0.001, 0, 1, normal_pdf, 0.01, 0.5; +end; +@#endif + +% Simulate data +stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500); +save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} ); +pause(1); + + + +%-------------------------------------------------------------------------- +% Method of Moments Estimation +%-------------------------------------------------------------------------- +% matched_moments blocks : We don't have an interface yet + +% get indices in declaration order +ic = strmatch('c', M_.endo_names,'exact'); +iiv = strmatch('iv', M_.endo_names,'exact'); +in = strmatch('n', M_.endo_names,'exact'); +% first entry: number of variable in declaration order +% second entry: lag +% third entry: power + +matched_moments_ = { + [ic ] [0 ], [1 ]; + [in ] [0 ], [1 ]; + [iiv ] [0 ], [1 ]; + [ic ic ] [0 0], [1 1]; + [ic iiv] [0 0], [1 1]; + [ic in ] [0 0], [1 1]; + [iiv ic ] [0 0], [1 1]; + [iiv iiv] [0 0], [1 1]; + [iiv in ] [0 0], [1 1]; +% [in ic ] [0 0], [1 1]; +% [in iiv] [0 0], [1 1]; + [in in ] [0 0], [1 1]; + [ic ic ] [0 -1], [1 1]; + [in in ] [0 -1], [1 1]; + [iiv iiv] [0 -1], [1 1]; +% [iiv iiv] [0 -1], [1 1]; +}; + + + +@#for mommethod in ["SMM"] + method_of_moments( + % Necessery options + mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM + , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data + + % Options for both GMM and SMM + % , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix + , order = @{orderApp} % order of Taylor approximation in perturbation + % , penalized_estimator % use penalized optimization + , pruning % use pruned state space system at higher-order + % , verbose % display and store intermediate estimation results + , weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + , additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute + % , prefilter=0 % demean each data series by its empirical mean and use centered moments + % + % Options for SMM + % , bounded_shock_support % trim shocks in simulation to +- 2 stdev + % , drop = 500 % number of periods dropped at beginning of simulation + % , seed = 24051986 % seed used in simulations + % , simulation_multiple = 5 % multiple of the data length used for simulation + % + % General options + %, dirname = 'MM' % directory in which to store estimation output + % , graph_format = EPS % specify the file format(s) for graphs saved to disk + % , nodisplay % do not display the graphs, but still save them to disk + % , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed) + % , noprint % do not print stuff to console + % , plot_priors = 1 % control plotting of priors + % , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters + % , TeX % print TeX tables and graphics + % + % Data and model options + %, first_obs = 501 % number of first observation + % , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data) + % , loglinear % computes a log-linear approximation of the model instead of a linear approximation + %, nobs = 500 % number of observations + % , xls_sheet = willi % name of sheet with data in Excel + % , xls_range = B2:D200 % range of data in Excel sheet + % + % Optimization options that can be set by the user in the mod file, otherwise default values are provided + % , analytic_derivation % uses analytic derivatives to compute standard errors for GMM + %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons + , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer + %, optim = ('TolFun', 1e-3 + % ,'TolX', 1e-5 + % ) % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute + %, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between + % , tolf = 1e-5 % convergence criterion on function value for numerical differentiation + % , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation + % + % % Numerical algorithms options + % , aim_solver % Use AIM algorithm to compute perturbation approximation + % , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION + % , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm + % , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm + % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm + % , k_order_solver % use k_order_solver in higher order perturbation approximations + % , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER + % , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver + % , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver + % , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm + % , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT + % , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver + % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl] + % , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition + ); +@#endfor + + + diff --git a/tests/estimation/method_of_moments/RBC_MoM_common.inc b/tests/estimation/method_of_moments/RBC_MoM_common.inc new file mode 100644 index 000000000..625a9c50b --- /dev/null +++ b/tests/estimation/method_of_moments/RBC_MoM_common.inc @@ -0,0 +1,80 @@ +% RBC model used in replication files of +% 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. +% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu) +% ========================================================================= +% Copyright (C) 2020 Dynare Team + +var k $K$ + c $C$ + a $A$ + iv $I$ + y $Y$ + la $\lambda$ + n $N$ + rk ${r^k}$ + w $W$ + ; + +predetermined_variables k; + +varexo u_a ${\varepsilon^{a}}$ + ; + +parameters DELTA $\delta$ + BETTA $\beta$ + B $B$ + ETAl $\eta_l$ + ETAc $\eta_c$ + THETA $\theta$ + ALFA $\alpha$ + RHOA $\rho_a$ + ; + +DELTA = 0.025; +BETTA = 0.984; +B = 0.5; +ETAl = 1; +ETAc = 2; +THETA = 3.48; +ALFA = 0.667; +RHOA = 0.979; + +model; +0 = -exp(la) +(exp(c)-B*exp(c(-1)))^(-ETAc) - BETTA*B*(exp(c(+1))-B*exp(c))^(-ETAc); +0 = -THETA*(1-exp(n))^-ETAl + exp(la)*exp(w); +0 = -exp(la) + BETTA*exp(la(+1))*(exp(rk(+1)) + (1-DELTA)); +0 = -exp(a)*(1-ALFA)*exp(k)^(-ALFA)*exp(n)^(ALFA) + exp(rk); +0 = -exp(a)*ALFA*exp(k)^(1-ALFA)*exp(n)^(ALFA-1) + exp(w); +0 = -exp(c) - exp(iv) + exp(y); +0 = -exp(y) + exp(a)*exp(k)^(1-ALFA)*exp(n)^(ALFA); +0 = -exp(k(+1)) + (1-DELTA)*exp(k) + exp(iv); +0 = -log(exp(a)) + RHOA*log(exp(a(-1))) + u_a; +end; + +steady_state_model; +A = 1; +RK = 1/BETTA - (1-DELTA); +K_O_N = (RK/(A*(1-ALFA)))^(-1/ALFA); +W = A*ALFA*(K_O_N)^(1-ALFA); +IV_O_N = DELTA*K_O_N; +Y_O_N = A*K_O_N^(1-ALFA); +C_O_N = Y_O_N - IV_O_N; + +N=RBC_MoM_steady_helper(THETA,ETAl,ETAc,BETTA,B,C_O_N,W); +C=C_O_N*N; +Y=Y_O_N*N; +IV=IV_O_N*N; +K=K_O_N*N; +LA = (C-B*C)^(-ETAc)-BETTA*B*(C-B*C)^(-ETAc); + +k=log(K); +c=log(C); +a=log(A); +iv=log(IV); +y=log(Y); +la=log(LA); +n=log(N); +rk=log(RK); +w=log(W); + +end; \ No newline at end of file diff --git a/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod b/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod new file mode 100644 index 000000000..de5d12d5f --- /dev/null +++ b/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod @@ -0,0 +1,161 @@ +% Tests SMM and GMM routines with prefilter, explicit initialization, and estimated_params_init(use_calibration); +% +% Copyright (C) 2020 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 . +% ========================================================================= + +% Define testscenario +@#define orderApp = 2 + +% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite +@#define optimizer = 13 + + +@#include "RBC_MoM_common.inc" + +shocks; +var u_a; stderr 0.0072; +end; + +varobs n c iv; + +% Simulate data +stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=250,TeX); +save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} ); +pause(1); + +set_param_value('DELTA',NaN); + +estimated_params; + DELTA, 0.025, 0, 1; + BETTA, , 0, 1; + B, , 0, 1; + %ETAl, 1, 0, 10; + ETAc, , 0, 10; + ALFA, , 0, 1; + RHOA, , 0, 1; + stderr u_a, , 0, 1; + %THETA, 3.48, 0, 10; +end; + +estimated_params_init(use_calibration); +end; + +%-------------------------------------------------------------------------- +% Method of Moments Estimation +%-------------------------------------------------------------------------- +% matched_moments blocks : We don't have an interface yet + +% get indices in declaration order +ic = strmatch('c', M_.endo_names,'exact'); +iiv = strmatch('iv', M_.endo_names,'exact'); +in = strmatch('n', M_.endo_names,'exact'); +% first entry: number of variable in declaration order +% second entry: lag +% third entry: power + +matched_moments_ = { + [ic ] [0 ], [1 ]; + [in ] [0 ], [1 ]; + [iiv ] [0 ], [1 ]; + [ic ic ] [0 0], [1 1]; + [ic iiv] [0 0], [1 1]; + [ic in ] [0 0], [1 1]; + [iiv ic ] [0 0], [1 1]; + [iiv iiv] [0 0], [1 1]; + [iiv in ] [0 0], [1 1]; +% [in ic ] [0 0], [1 1]; +% [in iiv] [0 0], [1 1]; + [in in ] [0 0], [1 1]; + [ic ic ] [0 -1], [1 1]; + [in in ] [0 -1], [1 1]; + [iiv iiv] [0 -1], [1 1]; +% [iiv iiv] [0 -1], [1 1]; +}; + +weighting_matrix=diag([1000;ones(6,1)]); +save('test_matrix.mat','weighting_matrix') + +@#for mommethod in ["GMM", "SMM"] + method_of_moments( + % Necessery options + mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM + , datafile = 'RBC_MoM_data_@{orderApp}.mat' % name of filename with data + + % Options for both GMM and SMM + % , bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix + , order = @{orderApp} % order of Taylor approximation in perturbation + % , penalized_estimator % use penalized optimization + , pruning % use pruned state space system at higher-order + % , verbose % display and store intermediate estimation results +% , weighting_matrix = 'test_matrix.mat' % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + , weighting_matrix =['test_matrix.mat','optimal'] + %, weighting_matrix = optimal % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename + %, additional_optimizer_steps = [4] % vector of additional mode-finders run after mode_compute + , prefilter=1 % demean each data series by its empirical mean and use centered moments + , se_tolx=1e-5 + % + % Options for SMM + % , bounded_shock_support % trim shocks in simulation to +- 2 stdev + , burnin = 500 % number of periods dropped at beginning of simulation + % , seed = 24051986 % seed used in simulations + % , simulation_multiple = 5 % multiple of the data length used for simulation + % + % General options + %, dirname = 'MM' % directory in which to store estimation output + % , graph_format = EPS % specify the file format(s) for graphs saved to disk + % , nodisplay % do not display the graphs, but still save them to disk + % , nograph % do not create graphs (which implies that they are not saved to the disk nor displayed) + % , noprint % do not print stuff to console + % , plot_priors = 1 % control plotting of priors + % , prior_trunc = 1e-10 % probability of extreme values of the prior density that is ignored when computing bounds for the parameters + % , TeX % print TeX tables and graphics + % + % Data and model options + %, first_obs = 501 % number of first observation + % , logdata % if loglinear is set, this option is necessary if the user provides data already in logs, otherwise the log transformation will be applied twice (this may result in complex data) + % , loglinear % computes a log-linear approximation of the model instead of a linear approximation + %, nobs = 500 % number of observations + % , xls_sheet = willi % name of sheet with data in Excel + % , xls_range = B2:D200 % range of data in Excel sheet + % + % Optimization options that can be set by the user in the mod file, otherwise default values are provided + % , analytic_derivation % uses analytic derivatives to compute standard errors for GMM + %, huge_number=1D10 % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons + , mode_compute = @{optimizer} % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer + %, optim = ('TolFun', 1e-3 + % ,'TolX', 1e-5 + % ) % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute + %, silent_optimizer % run minimization of moments distance silently without displaying results or saving files in between + % + % % Numerical algorithms options + % , aim_solver % Use AIM algorithm to compute perturbation approximation + % , dr=default % method used to compute the decision rule; possible values are DEFAULT, CYCLE_REDUCTION, LOGARITHMIC_REDUCTION + % , dr_cycle_reduction_tol = 1e-7 % convergence criterion used in the cycle reduction algorithm + % , dr_logarithmic_reduction_maxiter = 100 % maximum number of iterations used in the logarithmic reduction algorithm + % , dr_logarithmic_reduction_tol = 1e-12 % convergence criterion used in the cycle reduction algorithm + % , k_order_solver % use k_order_solver in higher order perturbation approximations + % , lyapunov = DEFAULT % algorithm used to solve lyapunov equations; possible values are DEFAULT, FIXED_POINT, DOUBLING, SQUARE_ROOT_SOLVER + % , lyapunov_complex_threshold = 1e-15 % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver + % , lyapunov_fixed_point_tol = 1e-10 % convergence criterion used in the fixed point Lyapunov solver + % , lyapunov_doubling_tol = 1e-16 % convergence criterion used in the doubling algorithm + % , sylvester = default % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT + % , sylvester_fixed_point_tol = 1e-12 % convergence criterion used in the fixed point Sylvester solver + % , qz_criterium = 0.999999 % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [IS THIS CORRET @wmutschl] + % , qz_zero_threshold = 1e-6 % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition + ); +@#endfor \ No newline at end of file