Code Review of GMM routines
- fix prefilter option - Implement iterative GMMtime-shift
parent
a40807caa9
commit
540f0454d2
|
@ -61,6 +61,7 @@ p = {'/distributions/' ; ...
|
||||||
'/cli/' ; ...
|
'/cli/' ; ...
|
||||||
'/lmmcp/' ; ...
|
'/lmmcp/' ; ...
|
||||||
'/optimization/' ; ...
|
'/optimization/' ; ...
|
||||||
|
'/method_of_moments/' ; ...
|
||||||
'/discretionary_policy/' ; ...
|
'/discretionary_policy/' ; ...
|
||||||
'/accessors/' ; ...
|
'/accessors/' ; ...
|
||||||
'/modules/dseries/src/' ; ...
|
'/modules/dseries/src/' ; ...
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
|
||||||
% -------------------------------------------------------------------------
|
|
||||||
% 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))<length(OptionsMoM.varobs)
|
|
||||||
for i = 1:OptionsMoM.obs_nbr
|
|
||||||
if sum(strcmp(OptionsMoM.varobs{i},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');
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% 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))<length(options_mom_.varobs)
|
||||||
|
for i = 1:options_mom_.obs_nbr
|
||||||
|
if sum(strcmp(options_mom_.varobs{i},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))'<estim_params_.nvx+estim_params_.nvn +estim_params_.ncx + estim_params_.ncn;
|
||||||
|
LB_below_minus_1=(Bounds.lb<-1 & param_of_interest);
|
||||||
|
UB_above_1=(Bounds.ub>1 & param_of_interest);
|
||||||
|
Bounds.lb(LB_below_minus_1)=-1;
|
||||||
|
Bounds.ub(UB_above_1)=1;
|
||||||
|
|
||||||
|
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');
|
|
@ -1,12 +1,12 @@
|
||||||
function [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_datamoments(data, DynareResults, MatchedMoments, OptionsMoM)
|
% [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_)
|
||||||
% This function computes the user-selected empirical moments from data
|
% This function computes the user-selected empirical moments from data
|
||||||
% =========================================================================
|
% =========================================================================
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% o data [T x varobs_nbr] data set
|
% o data [T x varobs_nbr] data set
|
||||||
% o DynareResults: [structure] storage for results (oo_)
|
% o oo_: [structure] storage for results
|
||||||
% o MatchedMoments: [structure] information about selected moments to match in estimation (matched_moments_)
|
% o matched_moments_: [structure] information about selected moments to match in estimation
|
||||||
% o OptionsMoM: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
|
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
|
||||||
% -------------------------------------------------------------------------
|
% -------------------------------------------------------------------------
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% o dataMoments [numMom x 1] mean of selected empirical moments
|
% 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
|
% Initialization
|
||||||
T = size(data,1); % Number of observations (T) and number of observables (ny)
|
T = size(data,1); % Number of observations (T) and number of observables (ny)
|
||||||
mom_nbr = OptionsMoM.mom_nbr;
|
dataMoments = NaN(options_mom_.mom_nbr,1);
|
||||||
dataMoments = nan(mom_nbr,1);
|
m_data = NaN(T,options_mom_.mom_nbr);
|
||||||
m_data = nan(T,mom_nbr);
|
% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*...
|
||||||
% Product moment for each time period, i.e. each row t contains yt1(l1)^p1*yt2(l2)^p2*...
|
|
||||||
% note that here we already are able to treat leads and lags and any power product moments
|
% note that here we already are able to treat leads and lags and any power product moments
|
||||||
for jm = 1:mom_nbr
|
for jm = 1:options_mom_.mom_nbr
|
||||||
vars = DynareResults.dr.inv_order_var(MatchedMoments{jm,1})';
|
vars = oo_.dr.inv_order_var(matched_moments_{jm,1})';
|
||||||
leadlags = MatchedMoments{jm,2}; % note that lags are negative numbers and leads are positive numbers
|
leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers
|
||||||
powers = MatchedMoments{jm,3};
|
powers = matched_moments_{jm,3};
|
||||||
for jv = 1:length(vars)
|
for jv = 1:length(vars)
|
||||||
jvar = DynareResults.dr.obs_var == vars(jv);
|
jvar = (oo_.dr.obs_var == vars(jv));
|
||||||
y = nan(T,1);
|
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);
|
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
|
||||||
if jv==1
|
if jv==1
|
||||||
m_data_tmp = y;
|
m_data_tmp = y;
|
||||||
else
|
else
|
||||||
m_data_tmp = m_data_tmp.*y;
|
m_data_tmp = m_data_tmp.*y;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
dataMoments(jm,1) = sum(m_data_tmp,'omitnan')/(T-sum(abs(leadlags)));
|
% We replace NaN (due to leads and lags and missing values) with the corresponding mean
|
||||||
% We replace nan (due to leads and lags) with the corresponding mean
|
dataMoments(jm,1) = mean(m_data_tmp,'omitnan');
|
||||||
% @wmutschl: this should also work for missing values, right?
|
|
||||||
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
|
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
|
||||||
m_data(:,jm) = m_data_tmp;
|
m_data(:,jm) = m_data_tmp;
|
||||||
end
|
end
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
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<eps
|
||||||
|
fprintf('Most negative variance %f for parameter %d (%s = %f)', s_min, k , bayestopt_.name{k}, xparam(k))
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
[nbplt,nr,nc,lr,lc,nstar] = pltorg(length(xparam));
|
||||||
|
|
||||||
|
if ~exist([M_.fname filesep 'graphs'],'dir')
|
||||||
|
mkdir(M_.fname,'graphs');
|
||||||
|
end
|
||||||
|
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
|
||||||
|
fidTeX = fopen([M_.fname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w');
|
||||||
|
fprintf(fidTeX,'%% TeX eps-loader file generated by method_of_moments_mode_check.m (Dynare).\n');
|
||||||
|
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
|
||||||
|
fprintf(fidTeX,' \n');
|
||||||
|
end
|
||||||
|
|
||||||
|
ll = options_.mode_check.neighbourhood_size;
|
||||||
|
if isinf(ll)
|
||||||
|
options_.mode_check.symmetric_plots = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
mcheck = struct('cross',struct(),'emode',struct());
|
||||||
|
|
||||||
|
for plt = 1:nbplt
|
||||||
|
if TeX
|
||||||
|
NAMES = [];
|
||||||
|
TeXNAMES = [];
|
||||||
|
end
|
||||||
|
hh = dyn_figure(options_.nodisplay,'Name','Mode check plots');
|
||||||
|
for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
|
||||||
|
subplot(nr,nc,k)
|
||||||
|
kk = (plt-1)*nstar+k;
|
||||||
|
[name,texname] = get_the_name(kk,TeX,M_,estim_params_,options_);
|
||||||
|
xx = xparam;
|
||||||
|
if xparam(kk)~=0 || ~isinf(Bounds.lb(kk)) || ~isinf(Bounds.lb(kk))
|
||||||
|
l1 = max(Bounds.lb(kk),(1-sign(xparam(kk))*ll)*xparam(kk)); m1 = 0; %lower bound
|
||||||
|
l2 = min(Bounds.ub(kk),(1+sign(xparam(kk))*ll)*xparam(kk)); %upper bound
|
||||||
|
else
|
||||||
|
%size info for 0 parameter is missing, use prior standard
|
||||||
|
%deviation
|
||||||
|
upper_bound=Bounds.lb(kk);
|
||||||
|
if isinf(upper_bound)
|
||||||
|
upper_bound=-1e-6*options_.huge_number;
|
||||||
|
end
|
||||||
|
lower_bound=Bounds.ub(kk);
|
||||||
|
if isinf(lower_bound)
|
||||||
|
lower_bound=-1e-6*options_.huge_number;
|
||||||
|
end
|
||||||
|
l1 = max(lower_bound,-bayestopt_.p2(kk)); m1 = 0; %lower bound
|
||||||
|
l2 = min(upper_bound,bayestopt_.p2(kk)); %upper bound
|
||||||
|
end
|
||||||
|
binding_lower_bound=0;
|
||||||
|
binding_upper_bound=0;
|
||||||
|
if isequal(xparam(kk),Bounds.lb(kk))
|
||||||
|
binding_lower_bound=1;
|
||||||
|
bound_value=Bounds.lb(kk);
|
||||||
|
elseif isequal(xparam(kk),Bounds.ub(kk))
|
||||||
|
binding_upper_bound=1;
|
||||||
|
bound_value=Bounds.ub(kk);
|
||||||
|
end
|
||||||
|
if options_.mode_check.symmetric_plots && ~binding_lower_bound && ~binding_upper_bound
|
||||||
|
if l2<(1+ll)*xparam(kk) %test whether upper bound is too small due to prior binding
|
||||||
|
l1 = xparam(kk) - (l2-xparam(kk)); %adjust lower bound to become closer
|
||||||
|
m1 = 1;
|
||||||
|
end
|
||||||
|
if ~m1 && (l1>(1-ll)*xparam(kk)) && (xparam(kk)+(xparam(kk)-l1)<Bounds.ub(kk)) % if lower bound was truncated and using difference from lower bound does not violate upper bound
|
||||||
|
l2 = xparam(kk) + (xparam(kk)-l1); %set upper bound to same distance as lower bound
|
||||||
|
end
|
||||||
|
end
|
||||||
|
z1 = l1:((xparam(kk)-l1)/(options_.mode_check.number_of_points/2)):xparam(kk);
|
||||||
|
z2 = xparam(kk):((l2-xparam(kk))/(options_.mode_check.number_of_points/2)):l2;
|
||||||
|
z = union(z1,z2);
|
||||||
|
if options_.mom.penalized_estimator
|
||||||
|
y = zeros(length(z),2);
|
||||||
|
dy=(xx-bayestopt_.p1)'/diag(bayestopt_.p2.^2)*(xx-bayestopt_.p1);
|
||||||
|
else
|
||||||
|
y = zeros(length(z),1);
|
||||||
|
end
|
||||||
|
for i=1:length(z)
|
||||||
|
xx(kk) = z(i);
|
||||||
|
[fval, info, exit_flag] = feval(fun,xx,varargin{:});
|
||||||
|
if exit_flag
|
||||||
|
y(i,1) = fval;
|
||||||
|
else
|
||||||
|
y(i,1) = NaN;
|
||||||
|
if options_.debug
|
||||||
|
fprintf('mode_check:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if options_.mom.penalized_estimator
|
||||||
|
prior=(xx-bayestopt_.p1)'/diag(bayestopt_.p2.^2)*(xx-bayestopt_.p1);
|
||||||
|
y(i,2) = (y(i,1)+prior-dy);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
mcheck.cross = setfield(mcheck.cross, name, [transpose(z), -y]);
|
||||||
|
mcheck.emode = setfield(mcheck.emode, name, xparam(kk));
|
||||||
|
fighandle=plot(z,-y);
|
||||||
|
hold on
|
||||||
|
yl=get(gca,'ylim');
|
||||||
|
plot( [xparam(kk) xparam(kk)], yl, 'c', 'LineWidth', 1)
|
||||||
|
NaN_index = find(isnan(y(:,1)));
|
||||||
|
zNaN = z(NaN_index);
|
||||||
|
yNaN = yl(1)*ones(size(NaN_index));
|
||||||
|
plot(zNaN,yNaN,'o','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',6);
|
||||||
|
if TeX
|
||||||
|
title(texname,'interpreter','latex')
|
||||||
|
else
|
||||||
|
title(name,'interpreter','none')
|
||||||
|
end
|
||||||
|
|
||||||
|
axis tight
|
||||||
|
if binding_lower_bound || binding_upper_bound
|
||||||
|
xl=get(gca,'xlim');
|
||||||
|
plot( [bound_value bound_value], yl, 'r--', 'LineWidth', 1)
|
||||||
|
xlim([xl(1)-0.5*binding_lower_bound*(xl(2)-xl(1)) xl(2)+0.5*binding_upper_bound*(xl(2)-xl(1))])
|
||||||
|
end
|
||||||
|
hold off
|
||||||
|
drawnow
|
||||||
|
end
|
||||||
|
if options_.mom.penalized_estimator
|
||||||
|
if isoctave
|
||||||
|
axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'),
|
||||||
|
else
|
||||||
|
axes('position',[0.3 0.01 0.42 0.05],'box','on'),
|
||||||
|
end
|
||||||
|
line_color=get(fighandle,'color');
|
||||||
|
plot([0.48 0.68],[0.5 0.5],'color',line_color{2})
|
||||||
|
hold on, plot([0.04 0.24],[0.5 0.5],'color',line_color{1})
|
||||||
|
set(gca,'xlim',[0 1],'ylim',[0 1],'xtick',[],'ytick',[])
|
||||||
|
text(0.25,0.5,'log-post')
|
||||||
|
text(0.69,0.5,'log-lik kernel')
|
||||||
|
end
|
||||||
|
dyn_saveas(hh,[M_.fname, '/graphs/', M_.fname '_MoMCheckPlots' int2str(plt) ],options_.nodisplay,options_.graph_format);
|
||||||
|
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
|
||||||
|
% TeX eps loader file
|
||||||
|
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
||||||
|
fprintf(fidTeX,'\\centering \n');
|
||||||
|
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_MoMCheckPlots%s}\n',options_.figures.textwidth*min(k/nc,1),[M_.fname, '/graphs/',M_.fname],int2str(plt));
|
||||||
|
fprintf(fidTeX,'\\caption{Method of Moments check plots.}');
|
||||||
|
fprintf(fidTeX,'\\label{Fig:MoMCheckPlots:%s}\n',int2str(plt));
|
||||||
|
fprintf(fidTeX,'\\end{figure}\n');
|
||||||
|
fprintf(fidTeX,' \n');
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
|
||||||
|
fclose(fidTeX);
|
||||||
|
end
|
||||||
|
|
||||||
|
OutputDirectoryName = CheckPath('modecheck',M_.dname);
|
||||||
|
save([OutputDirectoryName '/MoM_check_plot_data.mat'],'mcheck');
|
|
@ -0,0 +1,214 @@
|
||||||
|
function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_)
|
||||||
|
% [fval, info, exit_flag, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, matched_moments_, M_, options_mom_)
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% This function evaluates the objective function for GMM/SMM estimation
|
||||||
|
% =========================================================================
|
||||||
|
% INPUTS
|
||||||
|
% o xparam1: current value of estimated parameters as returned by set_prior()
|
||||||
|
% o Bounds: structure containing parameter bounds
|
||||||
|
% o oo_: structure for results
|
||||||
|
% o estim_params_: structure describing the estimated_parameters
|
||||||
|
% o matched_moments_: structure containing information about selected moments to match in estimation (matched_moments_)
|
||||||
|
% o M_ structure describing the model
|
||||||
|
% o options_mom_: 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 junk1: empty matrix required for optimizer interface
|
||||||
|
% o junk2: empty matrix required for optimizer interface
|
||||||
|
% o oo_: structure containing the results with the following updated fields:
|
||||||
|
% - mom.model_moments [numMom x 1] vector with model moments
|
||||||
|
% - mom.Q value of the quadratic form of the moment difference
|
||||||
|
% o M_: Matlab's structure describing the model
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% This function is called by
|
||||||
|
% o method_of_moments.m
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% This function calls
|
||||||
|
% o check_bounds_and_definiteness_estimation
|
||||||
|
% 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 <http://www.gnu.org/licenses/>.
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% 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
|
||||||
|
|
|
@ -1,17 +1,17 @@
|
||||||
function 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)
|
||||||
% Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag)
|
% 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
|
% 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
|
% INPUTS
|
||||||
% o m_data [T x numMom] selected empirical or theoretical moments at each point in time
|
% 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 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
|
% OUTPUTS
|
||||||
% o Wopt [numMom x numMom] optimal weighting matrix
|
% o W_opt [numMom x numMom] optimal weighting matrix
|
||||||
% -------------------------------------------------------------------------
|
% -------------------------------------------------------------------------
|
||||||
% This function is called by
|
% This function is called by
|
||||||
% o method_of_moments.m
|
% o method_of_moments.m
|
||||||
|
@ -42,45 +42,45 @@ function Wopt = method_of_moments_optimal_weighting_matrix(m_data, moments, qLag
|
||||||
% =========================================================================
|
% =========================================================================
|
||||||
|
|
||||||
% Initialize
|
% 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)
|
% center around moments (could be either data_moments or model_moments)
|
||||||
hFunc = m_data - repmat(moments',T,1);
|
h_Func = m_data - repmat(moments',T,1);
|
||||||
|
|
||||||
% The required correlation matrices
|
% The required correlation matrices
|
||||||
GAMA_array = zeros(numMom,numMom,qLag);
|
GAMA_array = zeros(num_Mom,num_Mom,q_lag);
|
||||||
GAMA0 = CorrMatrix(hFunc,T,numMom,0);
|
GAMA0 = Corr_Matrix(h_Func,T,num_Mom,0);
|
||||||
if qLag > 0
|
if q_lag > 0
|
||||||
for ii=1:qLag
|
for ii=1:q_lag
|
||||||
GAMA_array(:,:,ii) = CorrMatrix(hFunc,T,numMom,ii);
|
GAMA_array(:,:,ii) = Corr_Matrix(h_Func,T,num_Mom,ii);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
% The estimate of S
|
% The estimate of S
|
||||||
S = GAMA0;
|
S = GAMA0;
|
||||||
if qLag > 0
|
if q_lag > 0
|
||||||
for ii=1:qLag
|
for ii=1:q_lag
|
||||||
S = S + (1-ii/(qLag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)');
|
S = S + (1-ii/(q_lag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)');
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
% The estimate of W
|
% The estimate of W
|
||||||
Wopt = S\eye(size(S,1));
|
W_opt = S\eye(size(S,1));
|
||||||
|
|
||||||
% Check positive definite W
|
% Check positive definite W
|
||||||
try
|
try
|
||||||
chol(Wopt);
|
chol(W_opt);
|
||||||
catch err
|
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
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
% The correlation matrix
|
% The correlation matrix
|
||||||
function GAMAcorr = CorrMatrix(hFunc,T,numMom,v)
|
function GAMA_corr = Corr_Matrix(h_Func,T,num_Mom,v)
|
||||||
GAMAcorr = zeros(numMom,numMom);
|
GAMA_corr = zeros(num_Mom,num_Mom);
|
||||||
for t = 1+v:T
|
for t = 1+v:T
|
||||||
GAMAcorr = GAMAcorr + hFunc(t-v,:)'*hFunc(t,:);
|
GAMA_corr = GAMA_corr + h_Func(t-v,:)'*h_Func(t,:);
|
||||||
end
|
end
|
||||||
GAMAcorr = GAMAcorr/T;
|
GAMA_corr = GAMA_corr/T;
|
||||||
end
|
end
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
% -------------------------------------------------------------------------
|
||||||
|
% 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));
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
|
||||||
% -------------------------------------------------------------------------
|
|
||||||
% 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(xparam1<Bounds.lb)
|
|
||||||
if ~isequal(OptionsMoM.mode_compute,1)
|
|
||||||
k = find(xparam1<Bounds.lb);
|
|
||||||
fval = Inf;
|
|
||||||
exit_flag = 0;
|
|
||||||
info(1) = 41;
|
|
||||||
info(4)= sum((Bounds.lb(k)-xparam1(k)).^2);
|
|
||||||
return
|
|
||||||
elseif OptionsMoM.mode_compute == 13
|
|
||||||
fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number;
|
|
||||||
return
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the parameters.
|
|
||||||
if any(xparam1>Bounds.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
|
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
|
||||||
% -------------------------------------------------------------------------
|
|
||||||
% 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(xparam1<Bounds.lb)
|
|
||||||
if ~isequal(OptionsMoM.mode_compute,1)
|
|
||||||
k = find(xparam1<Bounds.lb);
|
|
||||||
fval = Inf;
|
|
||||||
exit_flag = 0;
|
|
||||||
info(1) = 41;
|
|
||||||
info(4)= sum((Bounds.lb(k)-xparam1(k)).^2);
|
|
||||||
return
|
|
||||||
elseif OptionsMoM.mode_compute == 13
|
|
||||||
fval = ones(size(DynareResults.mom.dataMoments,1),1)*OptionsMoM.huge_number;
|
|
||||||
return
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the parameters.
|
|
||||||
if any(xparam1>Bounds.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
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
|
||||||
% -------------------------------------------------------------------------
|
|
||||||
% 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));
|
|
|
@ -25,6 +25,8 @@ MODFILES = \
|
||||||
measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \
|
measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \
|
||||||
TeX/fs2000_corr_ME.mod \
|
TeX/fs2000_corr_ME.mod \
|
||||||
estimation/MH_recover/fs2000_recover_tarb.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 \
|
estimation/fs2000.mod \
|
||||||
gsa/ls2003a.mod \
|
gsa/ls2003a.mod \
|
||||||
optimizers/fs2000_8.mod \
|
optimizers/fs2000_8.mod \
|
||||||
|
@ -961,6 +963,8 @@ EXTRA_DIST = \
|
||||||
lmmcp/sw-common-header.inc \
|
lmmcp/sw-common-header.inc \
|
||||||
lmmcp/sw-common-footer.inc \
|
lmmcp/sw-common-footer.inc \
|
||||||
estimation/tune_mh_jscale/fs2000.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_unit_tests.m \
|
||||||
histval_initval_file/my_assert.m \
|
histval_initval_file/my_assert.m \
|
||||||
histval_initval_file/ramst_data.xls \
|
histval_initval_file/ramst_data.xls \
|
||||||
|
|
|
@ -25,7 +25,7 @@
|
||||||
@#define estimParams = 1
|
@#define estimParams = 1
|
||||||
|
|
||||||
% Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
|
% 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;
|
var c p R g y z INFL INT YGR;
|
||||||
varexo e_r e_g e_z;
|
varexo e_r e_g e_z;
|
||||||
|
@ -195,7 +195,7 @@ matched_moments_ = {
|
||||||
, 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
|
% , verbose % display and store intermediate estimation results
|
||||||
, weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
|
, 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
|
, 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
|
% , prefilter=0 % demean each data series by its empirical mean and use centered moments
|
||||||
%
|
%
|
||||||
% Options for SMM
|
% Options for SMM
|
||||||
|
|
|
@ -1,7 +1,5 @@
|
||||||
% RBC model used in replication files of
|
% Tests SMM and GMM routines
|
||||||
% 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
|
% Copyright (C) 2020 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
|
@ -21,92 +19,72 @@
|
||||||
% =========================================================================
|
% =========================================================================
|
||||||
|
|
||||||
% Define testscenario
|
% Define testscenario
|
||||||
@#define orderApp = 3
|
@#define orderApp = 1
|
||||||
@#define estimParams = 1
|
@#define estimParams = 1
|
||||||
|
|
||||||
% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
|
% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
|
||||||
@#define optimizer = 13
|
@#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;
|
@#include "RBC_MoM_common.inc"
|
||||||
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;
|
|
||||||
|
|
||||||
shocks;
|
shocks;
|
||||||
var u_a = 1;
|
var u_a; stderr 0.0072;
|
||||||
end;
|
end;
|
||||||
|
|
||||||
|
varobs n c iv;
|
||||||
|
|
||||||
|
|
||||||
@#if estimParams == 0
|
@#if estimParams == 0
|
||||||
estimated_params;
|
estimated_params;
|
||||||
DELTA, 0.02;
|
DELTA, 0.025;
|
||||||
BETTA, 0.9;
|
BETTA, 0.98;
|
||||||
B, 0.4;
|
B, 0.45;
|
||||||
%ETAl, 1;
|
%ETAl, 1;
|
||||||
ETAc, 1.5;
|
ETAc, 1.8;
|
||||||
ALFA, 0.6;
|
ALFA, 0.65;
|
||||||
RHOA, 0.9;
|
RHOA, 0.95;
|
||||||
STDA, 0.01;
|
stderr u_a, 0.01;
|
||||||
%THETA, 3.48;
|
%THETA, 3.48;
|
||||||
end;
|
end;
|
||||||
@#endif
|
@#endif
|
||||||
|
|
||||||
@#if estimParams == 1
|
@#if estimParams == 1
|
||||||
estimated_params;
|
estimated_params;
|
||||||
DELTA, 0.02, 0, 1;
|
DELTA, , 0, 1;
|
||||||
BETTA, 0.90, 0, 1;
|
BETTA, , 0, 1;
|
||||||
B, 0.40, 0, 1;
|
B, , 0, 1;
|
||||||
%ETAl, 1, 0, 10;
|
%ETAl, 1, 0, 10;
|
||||||
ETAc, 1.80, 0, 10;
|
ETAc, , 0, 10;
|
||||||
ALFA, 0.60, 0, 1;
|
ALFA, , 0, 1;
|
||||||
RHOA, 0.90, 0, 1;
|
RHOA, , 0, 1;
|
||||||
STDA, 0.01, 0, 1;
|
stderr u_a, , 0, 1;
|
||||||
%THETA, 3.48, 0, 10;
|
%THETA, 3.48, 0, 10;
|
||||||
end;
|
end;
|
||||||
@#endif
|
@#endif
|
||||||
|
|
||||||
@#if estimParams == 2
|
@#if estimParams == 2
|
||||||
estimated_params;
|
estimated_params;
|
||||||
DELTA, 0.02, 0, 1, normal_pdf, 0.02, 0.1;
|
DELTA, 0.025, 0, 1, normal_pdf, 0.02, 0.5;
|
||||||
BETTA, 0.90, 0, 1, normal_pdf, 0.90, 0.1;
|
BETTA, 0.98, 0, 1, beta_pdf, 0.90, 0.25;
|
||||||
B, 0.40, 0, 1, normal_pdf, 0.40, 0.1;
|
B, 0.45, 0, 1, normal_pdf, 0.40, 0.5;
|
||||||
%ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1;
|
%ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1;
|
||||||
ETAc, 1.80, 0, 10, normal_pdf, 1.80, 0.1;
|
ETAc, 1.8, 0, 10, normal_pdf, 1.80, 0.5;
|
||||||
ALFA, 0.60, 0, 1, normal_pdf, 0.60, 0.1;
|
ALFA, 0.65, 0, 1, normal_pdf, 0.60, 0.5;
|
||||||
RHOA, 0.90, 0, 1, normal_pdf, 0.90, 0.1;
|
RHOA, 0.95, 0, 1, normal_pdf, 0.90, 0.5;
|
||||||
STDA, 0.01, 0, 1, normal_pdf, 0.01, 0.1;
|
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;
|
%THETA, 3.48, 0, 10, normal_pdf, 0.25, 0.0.1;
|
||||||
end;
|
end;
|
||||||
@#endif
|
@#endif
|
||||||
|
|
||||||
% Simulate data
|
% 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{:} );
|
save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
|
||||||
pause(1);
|
pause(1);
|
||||||
|
|
||||||
|
|
||||||
|
estimated_params_init(use_calibration);
|
||||||
|
end;
|
||||||
|
|
||||||
%--------------------------------------------------------------------------
|
%--------------------------------------------------------------------------
|
||||||
% Method of Moments Estimation
|
% Method of Moments Estimation
|
||||||
|
@ -128,7 +106,7 @@ matched_moments_ = {
|
||||||
[ic ic ] [0 0], [1 1];
|
[ic ic ] [0 0], [1 1];
|
||||||
[ic iiv] [0 0], [1 1];
|
[ic iiv] [0 0], [1 1];
|
||||||
[ic in ] [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 iiv] [0 0], [1 1];
|
||||||
[iiv in ] [0 0], [1 1];
|
[iiv in ] [0 0], [1 1];
|
||||||
% [in ic ] [0 0], [1 1];
|
% [in ic ] [0 0], [1 1];
|
||||||
|
@ -154,8 +132,9 @@ matched_moments_ = {
|
||||||
% , penalized_estimator % use penalized optimization
|
% , 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
|
% , verbose % display and store intermediate estimation results
|
||||||
, weighting_matrix = OPTIMAL % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
|
, weighting_matrix = ['optimal','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_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
|
% , prefilter=0 % demean each data series by its empirical mean and use centered moments
|
||||||
%
|
%
|
||||||
% Options for SMM
|
% Options for SMM
|
||||||
|
@ -189,7 +168,7 @@ matched_moments_ = {
|
||||||
%, optim = ('TolFun', 1e-3
|
%, optim = ('TolFun', 1e-3
|
||||||
% ,'TolX', 1e-5
|
% ,'TolX', 1e-5
|
||||||
% ) % a list of NAME and VALUE pairs to set options for the optimization routines. Available options depend on mode_compute
|
% ) % 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
|
% , tolf = 1e-5 % convergence criterion on function value for numerical differentiation
|
||||||
% , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
|
% , tolx = 1e-6 % convergence criterion on funciton input for numerical differentiation
|
||||||
%
|
%
|
||||||
|
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
% =========================================================================
|
||||||
|
|
||||||
|
% 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
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -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;
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||||
|
% =========================================================================
|
||||||
|
|
||||||
|
% 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
|
Loading…
Reference in New Issue