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

Merge branch 'moment_estimation' into 'master' 

First implementation of moment estimation

See merge request Dynare/dynare!1750
time-shift
Sébastien Villemot 2020-07-16 15:41:15 +00:00
parent cda3253d2c 19b7d202d2
commit 1dbbef9f2e
21 changed files with 2633 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
matlab/dynare_config.m
View File

@ -61,6 +61,7 @@ p = {'/distributions/' ; ...
'/cli/' ; ...
'/lmmcp/' ; ...
'/optimization/' ; ...
'/method_of_moments/' ; ...
'/discretionary_policy/' ; ...
'/accessors/' ; ...
'/modules/dseries/src/' ; ...

23
matlab/dynare_estimation_init.m
View File

@ -626,28 +626,7 @@ end
%% get the non-zero rows and columns of Sigma_e and H
H_non_zero_rows=find(~all(M_.H==0,1));
H_non_zero_columns=find(~all(M_.H==0,2));
if ~isequal(H_non_zero_rows,H_non_zero_columns')
error('Measurement error matrix not symmetric')
end
if isfield(estim_params_,'nvn_observable_correspondence')
estim_params_.H_entries_to_check_for_positive_definiteness=union(H_non_zero_rows,estim_params_.nvn_observable_correspondence(:,1));
else
estim_params_.H_entries_to_check_for_positive_definiteness=H_non_zero_rows;
end
Sigma_e_non_zero_rows=find(~all(M_.Sigma_e==0,1));
Sigma_e_non_zero_columns=find(~all(M_.Sigma_e==0,2));
if ~isequal(Sigma_e_non_zero_rows,Sigma_e_non_zero_columns')
error('Structual error matrix not symmetric')
end
if isfield(estim_params_,'var_exo') && ~isempty(estim_params_.var_exo)
estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=union(Sigma_e_non_zero_rows,estim_params_.var_exo(:,1));
else
estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=Sigma_e_non_zero_rows;
end
estim_params_= get_matrix_entries_for_psd_check(M_,estim_params_);
if options_.load_results_after_load_mh
if ~exist([M_.fname '_results.mat'],'file')

55
matlab/get_matrix_entries_for_psd_check.m Normal file
View File

@ -0,0 +1,55 @@
function estim_params_= get_matrix_entries_for_psd_check(M_,estim_params_)
% function estim_params_= get_matrix_entries_for_psd_check(M_)
% Get entries of Sigma_e and H to check for positive definiteness
%
% INPUTS
% M_: structure storing the model information
% estim_params_: structure storing information about estimated
% parameters
% OUTPUTS
% estim_params_: structure storing information about estimated
% parameters
%
% SPECIAL REQUIREMENTS
% none
% 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/>.
%% get the non-zero rows and columns of Sigma_e and H
H_non_zero_rows=find(~all(M_.H==0,1));
H_non_zero_columns=find(~all(M_.H==0,2));
if ~isequal(H_non_zero_rows,H_non_zero_columns') || (any(any(M_.H-M_.H'>1e-10)))
error('Measurement error matrix not symmetric')
end
if isfield(estim_params_,'nvn_observable_correspondence')
estim_params_.H_entries_to_check_for_positive_definiteness=union(H_non_zero_rows,estim_params_.nvn_observable_correspondence(:,1));
else
estim_params_.H_entries_to_check_for_positive_definiteness=H_non_zero_rows;
end
Sigma_e_non_zero_rows=find(~all(M_.Sigma_e==0,1));
Sigma_e_non_zero_columns=find(~all(M_.Sigma_e==0,2));
if ~isequal(Sigma_e_non_zero_rows,Sigma_e_non_zero_columns') || (any(any(M_.Sigma_e-M_.Sigma_e'>1e-10)))
error('Structual error matrix not symmetric')
end
if isfield(estim_params_,'var_exo') && ~isempty(estim_params_.var_exo)
estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=union(Sigma_e_non_zero_rows,estim_params_.var_exo(:,1));
else
estim_params_.Sigma_e_entries_to_check_for_positive_definiteness=Sigma_e_non_zero_rows;
end

918
matlab/method_of_moments/method_of_moments.m Normal file
View File

@ -0,0 +1,918 @@
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
% - Initialize other options
% - Get variable orderings and state space representation
% Step 2: Checks and transformations for matched moments structure
% Step 3: Checks and transformations for estimated parameters, priors, and bounds
% Step 4: Checks and transformations for data
% Step 5: Checks for steady state at initial parameters
% Step 6: Checks for objective function at initial parameters
% Step 7: Iterated method of moments estimation
% Step 8: J-Test
% Step 9: Clean up
% -------------------------------------------------------------------------
% 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 (user-specified and updated) 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 get_matrix_entries_for_psd_check.m
% o makedataset.m
% o method_of_moments_data_moments.m
% o method_of_moments_mode_check.m
% o method_of_moments_objective_function.m
% o method_of_moments_optimal_weighting_matrix
% o method_of_moments_standard_errors
% o plot_priors.m
% o print_info.m
% o prior_bounds.m
% o set_default_option.m
% o set_prior.m
% o set_state_space.m
% o set_all_parameters.m
% o test_for_deep_parameters_calibration.m
% =========================================================================
% 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
% - [ ] do we need dirname?
% - [ ] decide on default weighting matrix scheme, I would propose 2 stage with Diagonal of optimal matrix
% - [ ] check smm with product moments greater than 2
% -------------------------------------------------------------------------
% 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')
else
options_mom_.logged_steady_state = 0;
options_mom_.loglinear = false;
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); % include deviation from prior mean as additional moment restriction and use prior precision as weight
options_mom_.mom = set_default_option(options_mom_.mom,'verbose',false); % display and store intermediate estimation results
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix',{'DIAGONAL'; 'DIAGONAL'}); % weighting matrix in moments distance objective function at each iteration of estimation; cell of strings with
% possible values are 'OPTIMAL', 'IDENTITY_MATRIX' ,'DIAGONAL' or a filename. Size of cell determines stages in iterated estimation.
options_mom_.mom = set_default_option(options_mom_.mom,'weighting_matrix_scaling_factor',1); % scaling of weighting matrix
options_mom_.mom = set_default_option(options_mom_.mom,'se_tolx',1e-5); % step size for numerical computation of standard errors
options_mom_ = set_default_option(options_mom_,'order',1); % order of Taylor approximation in perturbation
options_mom_ = set_default_option(options_mom_,'pruning',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 5.\n')
options_mom_.mom.simulation_multiple = 5;
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 data is already in logs
options_mom_ = set_default_option(options_mom_,'nobs',NaN); % number of observations
options_mom_ = set_default_option(options_mom_,'prefilter',false); % demean each data series by its empirical mean and use centered moments
options_mom_ = set_default_option(options_mom_,'xls_sheet',1); % name of sheet with data in Excel
options_mom_ = set_default_option(options_mom_,'xls_range',''); % range of data in Excel sheet
% Recursive estimation and forecast are not supported
if numel(options_mom_.nobs)>1
error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''nobs''.');
end
if numel(options_mom_.first_obs)>1
error('method_of_moments: Recursive estimation and forecast for samples is not supported. Please set an integer as ''first_obs''.');
end
% Optimization options that can be set by the user in the mod file, otherwise default values are provided
options_mom_ = set_default_option(options_mom_,'huge_number',1e7); % value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons
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_,'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
% Mode_check plot options that can be set by the user in the mod file, otherwise default values are provided
options_mom_.mode_check.nolik = false; % we don't do likelihood (also this initializes mode_check substructure)
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'status',false); % plot the target function for values around the computed mode for each estimated parameter in turn. This is helpful to diagnose problems with the optimizer.
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'neighbourhood_size',.5); % width of the window around the mode to be displayed on the diagnostic plots. This width is expressed in percentage deviation. The Inf value is allowed, and will trigger a plot over the entire domain
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'symmetric_plots',true); % ensure that the check plots are symmetric around the mode. A value of 0 allows to have asymmetric plots, which can be useful if the posterior mode is close to a domain boundary, or in conjunction with mode_check_neighbourhood_size = Inf when the domain is not the entire real line
options_mom_.mode_check = set_default_option(options_mom_.mode_check,'number_of_points',20); % number of points around the mode where the target function is evaluated (for each parameter)
% Numerical algorithms options that can be set by the user in the mod file, otherwise default values are provided
options_mom_ = set_default_option(options_mom_,'aim_solver',false); % use AIM algorithm to compute perturbation approximation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'k_order_solver',false); % use k_order_perturbation instead of mjdgges
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction',false); % use cycle reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_cycle_reduction_tol',1e-7); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction',false); % use logarithmic reduction algorithm to solve the polynomial equation for retrieving the coefficients associated to the endogenous variables in the decision rule
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_maxiter',100); % maximum number of iterations used in the logarithmic reduction algorithm
options_mom_ = set_default_option(options_mom_,'dr_logarithmic_reduction_tol',1e-12); % convergence criterion used in the cycle reduction algorithm
options_mom_ = set_default_option(options_mom_,'lyapunov_db',false); % doubling algorithm (disclyap_fast) to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_fp',false); % fixed-point algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false); % square-root-solver (dlyapchol) algorithm to solve Lyapunov equation to compute variance-covariance matrix of state variables
options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15); % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10); % convergence criterion used in the fixed point Lyapunov solver
options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16); % convergence criterion used in the doubling algorithm
options_mom_ = set_default_option(options_mom_,'sylvester_fp',false); % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12); % convergence criterion used in the fixed point Sylvester solver
options_mom_ = set_default_option(options_mom_,'qz_criterium',1-1e-6); % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems [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
if options_mom_.order > 2
fprintf('Dynare will use ''k_order_solver'' as the order>2\n');
options_mom_.k_order_solver = true;
end
% -------------------------------------------------------------------------
% Step 1b: Options that are set by the preprocessor and need to be carried over
% -------------------------------------------------------------------------
% 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_.solve_tolx = options_.solve_tolx;
options_mom_.steady = options_.steady;
options_mom_.steadystate = options_.steadystate;
options_mom_.steadystate_flag = options_.steadystate_flag;
% options related to dataset
options_mom_.dataset = options_.dataset;
options_mom_.initial_period = options_.initial_period;
% options related to endogenous prior restrictions are not supported
options_mom_.endogenous_prior_restrictions.irf = {};
options_mom_.endogenous_prior_restrictions.moment = {};
if ~isempty(options_.endogenous_prior_restrictions.irf) && ~isempty(options_.endogenous_prior_restrictions.moment)
fprintf('Endogenous prior restrictions are not supported yet and will be skipped.\n')
end
% -------------------------------------------------------------------------
% Step 1c: Options related to optimizers
% -------------------------------------------------------------------------
% mode_compute = 1, 3, 7, 11, 102, 11, 13
% nothing to be done
% mode_compute = 2
options_mom_.saopt = options_.saopt;
% mode_compute = 4
options_mom_.csminwel = options_.csminwel;
% mode_compute = 5
options_mom_.newrat = options_.newrat;
options_mom_.gstep = options_.gstep;
% mode_compute = 6
options_mom_.gmhmaxlik = options_.gmhmaxlik;
options_mom_.mh_jscale = options_.mh_jscale;
% mode_compute = 8
options_mom_.simplex = options_.simplex;
% mode_compute = 9
options_mom_.cmaes = options_.cmaes;
% mode_compute = 10
options_mom_.simpsa = options_.simpsa;
% mode_compute = 12
options_mom_.particleswarm = options_.particleswarm;
% mode_compute = 101
options_mom_.solveopt = options_.solveopt;
options_mom_.gradient_method = options_.gradient_method;
options_mom_.gradient_epsilon = options_.gradient_epsilon;
options_mom_.analytic_derivation = 0;
options_mom_.vector_output= false; % specifies whether the objective function returns a vector
% -------------------------------------------------------------------------
% 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_.mom.index.E_y = false(options_mom_.obs_nbr,1); %unconditional first order product moments
options_mom_.mom.index.E_yy = false(options_mom_.obs_nbr,options_mom_.obs_nbr); %unconditional second order product moments
options_mom_.mom.index.E_yyt = false(options_mom_.obs_nbr,options_mom_.obs_nbr,0); %unconditional temporal second order product moments
options_mom_.mom.index.E_y_pos = zeros(options_mom_.obs_nbr,1); %position in matched moments block
options_mom_.mom.index.E_yy_pos = zeros(options_mom_.obs_nbr,options_mom_.obs_nbr); %position in matched moments block
options_mom_.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_.mom.index.E_y(vpos,1) = true;
options_mom_.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_.mom.index.E_yy(idx1,idx2) = true;
options_mom_.mom.index.E_yy(idx2,idx1) = true;
options_mom_.mom.index.E_yy_pos(idx1,idx2) = jm;
options_mom_.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_.mom.index.E_yyt(idx1,idx2,-lag2) = true;
options_mom_.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_.mom.index.E_y_pos); nonzeros(tril(options_mom_.mom.index.E_yy_pos)); nonzeros(options_mom_.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_.mom.index
matched_moments_ = matched_moments_(UniqueMomIdx,:);
if strcmp(options_mom_.mom.mom_method,'SMM')
options_mom_.mom=rmfield(options_mom_.mom,'index');
end
% Check if both prefilter and first moments were specified
options_mom_.mom.first_moment_indicator = find(cellfun(@(x) sum(abs(x))==1,matched_moments_(:,3)))';
if options_mom_.prefilter && ~isempty(options_mom_.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_.mom.first_moment_indicator');
matched_moments_(options_mom_.mom.first_moment_indicator,:)=[]; %remove first moments entries
options_mom_.mom.first_moment_indicator = [];
end
options_mom_.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.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')
if isoctave
warning('off','Octave:singular-matrix');
else
warning('off','MATLAB:singularMatrix');
end
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 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;
if options_mom_.mom.penalized_estimator
fprintf('Penalized estimation turned off as you did not declare priors\n')
options_mom_.mom.penalized_estimator = 0;
end
end
% 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
% provide info on missing observations
if any(any(isnan(dataset_.data)))
fprintf('missing observations will be replaced by the sample mean of the corresponding moment')
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_.mom.long = round(options_mom_.mom.simulation_multiple*options_mom_.nobs);
options_mom_.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_.mom.long+options_mom_.mom.burnin,M_.exo_nbr);
temp_shocks_ME = randn(smmstream,options_mom_.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_.mom.shock_series = temp_shocks;
options_mom_.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
% -------------------------------------------------------------------------
objective_function = str2func('method_of_moments_objective_function');
try
% Check for NaN or complex values of moment-distance-funtion evaluated
% at initial parameters and identity weighting matrix
oo_.mom.Sw = eye(options_mom_.mom.mom_nbr);
tic_id = tic;
[fval, info, ~, ~, ~, oo_, M_] = feval(objective_function, xparam0, Bounds, oo_, estim_params_, 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 deviation from prior mean and weighted with prior precision');
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.mom_nbr);
fprintf('\n - number of parameters: %d\n\n', length(xparam0));
% -------------------------------------------------------------------------
% Step 7b: Iterated method of moments estimation
% -------------------------------------------------------------------------
if size(options_mom_.mom.weighting_matrix,1)>1 && ~(any(strcmpi('diagonal',options_mom_.mom.weighting_matrix)) || any(strcmpi('optimal',options_mom_.mom.weighting_matrix)))
fprintf('\nYou did not specify the use of an optimal or diagonal weighting matrix. There is no point in running an iterated method of moments.\n')
end
optimizer_vec=[options_mom_.mode_compute,options_mom_.additional_optimizer_steps]; % at each stage one can possibly use different optimizers sequentially
for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
fprintf('Estimation stage %u\n',stage_iter);
Woptflag = false;
switch lower(options_mom_.mom.weighting_matrix{stage_iter})
case 'identity_matrix'
fprintf(' - identity weighting matrix\n');
weighting_matrix = eye(options_mom_.mom.mom_nbr);
case 'diagonal'
fprintf(' - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = diag(diag( method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag) ));
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = diag(diag( method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag) ));
end
case 'optimal'
fprintf(' - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
if stage_iter == 1
fprintf(' and using data-moments as initial estimate of model-moments\n');
weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
else
fprintf(' and using previous stage estimate of model-moments\n');
weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
Woptflag = true;
end
otherwise %user specified matrix in file
fprintf(' - user-specified weighting matrix\n');
try
load(options_mom_.mom.weighting_matrix{stage_iter},'weighting_matrix')
catch
error(['method_of_moments: No matrix named ''weighting_matrix'' could be found in ',options_mom_.mom.weighting_matrix{stage_iter},'.mat'])
end
[nrow, ncol] = size(weighting_matrix);
if ~isequal(nrow,ncol) || ~isequal(nrow,length(oo_.mom.data_moments)) %check if square and right size
error(['method_of_moments: weighting_matrix must be square and have ',num2str(length(oo_.mom.data_moments)),' rows and columns'])
end
end
try %check for positive definiteness of weighting_matrix
oo_.mom.Sw = chol(weighting_matrix);
catch
error('method_of_moments: Specified weighting_matrix is not positive definite. Check whether your model implies stochastic singularity.')
end
for optim_iter= 1:length(optimizer_vec)
if optimizer_vec(optim_iter)==13
options_mom_.vector_output = true;
else
options_mom_.vector_output = false;
end
[xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, optimizer_vec(optim_iter), options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
Bounds, oo_, estim_params_, 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,optim_iter,fval)
if options_mom_.mom.verbose
oo_.mom=display_estimation_results_table(xparam1,NaN(size(xparam1)),M_,options_mom_,estim_params_,bayestopt_laplace,oo_.mom,prior_dist_names,sprintf('%s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter),sprintf('verbose_%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter));
end
xparam0=xparam1;
end
options_mom_.vector_output = false;
% Update M_ and DynareResults (in particular to get oo_.mom.model_moments)
M_ = set_all_parameters(xparam1,estim_params_,M_);
[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 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
% -------------------------------------------------------------------------
% Step 8: J test
% -------------------------------------------------------------------------
if options_mom_.mom.mom_nbr > length(xparam1)
%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_.mom.variance_correction_factor;
elseif strcmp(options_mom_.mom.mom_method,'GMM')
Variance_correction_factor=1;
end
oo_.mom.J_test.j_stat = dataset_.nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam1);
oo_.mom.J_test.p_val = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
fprintf('\nvalue of J-test statistic: %f\n',oo_.mom.J_test.j_stat)
fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val)
end
% -------------------------------------------------------------------------
% Step 9: Display estimation results
% -------------------------------------------------------------------------
title = ['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 9: Clean up
% -------------------------------------------------------------------------
%reset warning state
if isoctave
warning('on')
else
warning on
end

71
matlab/method_of_moments/method_of_moments_data_moments.m Normal file
View File

@ -0,0 +1,71 @@
function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_)
% [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, matched_moments_, options_mom_)
% This function computes the user-selected empirical moments from data
% =========================================================================
% INPUTS
% o data [T x varobs_nbr] data set
% o oo_: [structure] storage for results
% o matched_moments_: [structure] information about selected moments to match in estimation
% o options_mom_: [structure] information about all settings (specified by the user, preprocessor, and taken from global options_)
% -------------------------------------------------------------------------
% OUTPUTS
% o dataMoments [numMom x 1] mean of selected empirical moments
% o m_data [T x numMom] selected empirical moments at each point in time
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% o method_of_moments_objective_function.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)
% =========================================================================
% Initialization
T = size(data,1); % Number of observations (T)
dataMoments = NaN(options_mom_.mom.mom_nbr,1);
m_data = NaN(T,options_mom_.mom.mom_nbr);
% Product moment for each time period, i.e. each row t contains y_t1(l1)^p1*y_t2(l2)^p2*...
% note that here we already are able to treat leads and lags and any power product moments
for jm = 1:options_mom_.mom.mom_nbr
vars = oo_.dr.inv_order_var(matched_moments_{jm,1})';
leadlags = matched_moments_{jm,2}; % lags are negative numbers and leads are positive numbers
powers = matched_moments_{jm,3};
for jv = 1:length(vars)
jvar = (oo_.dr.obs_var == vars(jv));
y = NaN(T,1); %Take care of T_eff instead of T for lags and NaN via mean with 'omitnan' option below
y( (1-min(leadlags(jv),0)) : (T-max(leadlags(jv),0)), 1) = data( (1+max(leadlags(jv),0)) : (T+min(leadlags(jv),0)), jvar).^powers(jv);
if jv==1
m_data_tmp = y;
else
m_data_tmp = m_data_tmp.*y;
end
end
% We replace NaN (due to leads and lags and missing values) with the corresponding mean
dataMoments(jm,1) = mean(m_data_tmp,'omitnan');
m_data_tmp(isnan(m_data_tmp)) = dataMoments(jm,1);
m_data(:,jm) = m_data_tmp;
end
end %function end

185
matlab/method_of_moments/method_of_moments_mode_check.m Normal file
View File

@ -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');

213
matlab/method_of_moments/method_of_moments_objective_function.m Normal file
View File

@ -0,0 +1,213 @@
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, junk1, junk2, 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
% 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 error, 1 if no 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)
% =========================================================================
%------------------------------------------------------------------------------
% 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.mom_nbr,1);
offset = 0;
% First moments
if ~options_mom_.prefilter && isfield(options_mom_.mom.index,'E_y') && nnz(options_mom_.mom.index.E_y) > 0
E_y = pruned_state_space.E_y;
E_y_nbr = nnz(options_mom_.mom.index.E_y);
oo_.mom.model_moments(offset+1:E_y_nbr,1) = E_y(options_mom_.mom.index.E_y);
offset = offset + E_y_nbr;
end
% Second moments
% Contemporaneous covariance
if isfield(options_mom_.mom.index,'E_yy') && nnz(options_mom_.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(tril(options_mom_.mom.index.E_yy));
oo_.mom.model_moments(offset+(1:E_yy_nbr),1) = E_yy(tril(options_mom_.mom.index.E_yy));
offset = offset + E_yy_nbr;
end
% Lead/lags covariance
if isfield(options_mom_.mom.index,'E_yyt') && nnz(options_mom_.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_.mom.index.E_yyt);
oo_.mom.model_moments(offset+(1:E_yyt_nbr),1) = E_yyt(options_mom_.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_.mom.shock_series)); %initialize
scaled_shock_series(:,i_exo_var) = options_mom_.mom.shock_series(:,i_exo_var)*chol_S; %set non-zero entries
% simulate series
y_sim = simult_(M_, options_mom_, dr.ys, dr, scaled_shock_series, options_mom_.order);
% provide meaningful penalty if data is nan or inf
if any(any(isnan(y_sim))) || any(any(isinf(y_sim)))
if options_mom_.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_.mom.long+1:end)';
if ~all(diag(M_.H)==0)
i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance
chol_S = chol(M_.H(i_ME,i_ME)); %decompose rest
shock_mat=zeros(size(options_mom_.mom.ME_shock_series)); %initialize
shock_mat(:,i_ME)=options_mom_.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

79
matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m Normal file
View File

@ -0,0 +1,79 @@
function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
% W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
% -------------------------------------------------------------------------
% This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag
% Adapted from replication codes of
% o Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
% =========================================================================
% INPUTS
% o m_data [T x numMom] selected data moments at each point in time
% o moments [numMom x 1] selected estimated moments (either data_moments or estimated model_moments)
% o q_lag [integer] Bartlett kernel maximum lag order
% -------------------------------------------------------------------------
% OUTPUTS
% o W_opt [numMom x numMom] optimal weighting matrix
% -------------------------------------------------------------------------
% This function is called by
% o method_of_moments.m
% -------------------------------------------------------------------------
% This function calls:
% o CorrMatrix (embedded)
% =========================================================================
% 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)
% =========================================================================
% Initialize
[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 data_moments or model_moments)
h_Func = m_data - repmat(moments',T,1);
% The required correlation matrices
GAMA_array = zeros(num_Mom,num_Mom,q_lag);
GAMA0 = Corr_Matrix(h_Func,T,num_Mom,0);
if q_lag > 0
for ii=1:q_lag
GAMA_array(:,:,ii) = Corr_Matrix(h_Func,T,num_Mom,ii);
end
end
% The estimate of S
S = GAMA0;
if q_lag > 0
for ii=1:q_lag
S = S + (1-ii/(q_lag+1))*(GAMA_array(:,:,ii) + GAMA_array(:,:,ii)');
end
end
% The estimate of W
W_opt = S\eye(size(S,1));
end
% The correlation matrix
function GAMA_corr = Corr_Matrix(h_Func,T,num_Mom,v)
GAMA_corr = zeros(num_Mom,num_Mom);
for t = 1+v:T
GAMA_corr = GAMA_corr + h_Func(t-v,:)'*h_Func(t,:);
end
GAMA_corr = GAMA_corr/T;
end

104
matlab/method_of_moments/method_of_moments_standard_errors.m Normal file
View File

@ -0,0 +1,104 @@
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
WW = oo_.mom.Sw'*oo_.mom.Sw;
if Wopt_flag
% We have the optimal weighting matrix
Asympt_Var = 1/T*((D'*WW*D)\eye(dim_params));
else
% We do not have the optimal weighting matrix yet
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'*WW*D)\eye(dim_params);
Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA;
end
SE_values = sqrt(diag(Asympt_Var));

17
matlab/optimization/dynare_minimize_objective.m
View File

@ -428,7 +428,7 @@ switch minimizer_algorithm
case 12
if isoctave
error('Option mode_compute=12 is not available under Octave')
elseif ~user_has_matlab_license('global_optimization_toolbox')
elseif ~user_has_matlab_license('GADS_Toolbox')
error('Option mode_compute=12 requires the Global Optimization Toolbox')
end
[LB, UB] = set_bounds_to_finite_values(bounds, options_.huge_number);
@ -523,6 +523,21 @@ switch minimizer_algorithm
end
func = @(x)objective_function(x,varargin{:});
[opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options);
case 13
% Matlab's lsqnonlin (Optimization toolbox needed).
if isoctave && ~user_has_octave_forge_package('optim')
error('Option mode_compute=13 requires the optim package')
elseif ~isoctave && ~user_has_matlab_license('optimization_toolbox')
error('Option mode_compute=13 requires the Optimization Toolbox')
end
optim_options = optimset('display','iter','MaxFunEvals',5000,'MaxIter',5000,'TolFun',1e-6,'TolX',1e-6);
if ~isempty(options_.optim_opt)
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
end
if options_.silent_optimizer
optim_options = optimset(optim_options,'display','off');
end
[opt_par_values,Resnorm,fval,exitflag,OUTPUT,LAMBDA,JACOB] = lsqnonlin(objective_function,start_par_value,bounds(:,1),bounds(:,2),optim_options,varargin{:});
otherwise
if ischar(minimizer_algorithm)
if exist(minimizer_algorithm)

21
matlab/plot_priors.m
View File

@ -1,19 +1,21 @@
function plot_priors(bayestopt_,M_,estim_params_,options_)
function plot_priors(bayestopt_,M_,estim_params_,options_,optional_title)
% function plot_priors
% plots prior density
%
% INPUTS
% o bayestopt_ [structure]
% o M_ [structure]
% o options_ [structure]
%
% o bayestopt_ [structure]
% o M_ [structure]
% o estim_params_ [structure]
% o options_ [structure]
% o optional_title [string]
% OUTPUTS
% None
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2004-2017 Dynare Team
% Copyright (C) 2004-2020 Dynare Team
%
% This file is part of Dynare.
%
@ -31,8 +33,11 @@ function plot_priors(bayestopt_,M_,estim_params_,options_)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
TeX = options_.TeX;
figurename = 'Priors';
if nargin<5
figurename = 'Priors';
else
figurename = optional_title;
end
npar = length(bayestopt_.p1);
[nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);

2
tests/.gitignore vendored
View File

@ -50,6 +50,8 @@ wsOct
!/ep/mean_preserving_spread.m
!/ep/rbcii_steady_state.m
!/estimation/fsdat_simul.m
!/estimation/method_of_moments/RBC_MoM_steady_helper.m
!/estimation/method_of_moments/RBC_Andreasen_Data_2.mat
!/expectations/expectation_ss_old_steadystate.m
!/external_function/extFunDeriv.m
!/external_function/extFunNoDerivs.m

6
tests/Makefile.am
View File

@ -47,6 +47,10 @@ MODFILES = \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/t_proposal/fs2000_student.mod \
estimation/tune_mh_jscale/fs2000.mod \
estimation/method_of_moments/AnScho_MoM.mod \
estimation/method_of_moments/RBC_MoM_Andreasen.mod \
estimation/method_of_moments/RBC_MoM_SMM_ME.mod \
estimation/method_of_moments/RBC_MoM_prefilter.mod \
moments/example1_var_decomp.mod \
moments/example1_bp_test.mod \
moments/test_AR1_spectral_density.mod \
@ -959,6 +963,8 @@ EXTRA_DIST = \
lmmcp/sw-common-header.inc \
lmmcp/sw-common-footer.inc \
estimation/tune_mh_jscale/fs2000.inc \
estimation/method_of_moments/RBC_MoM_common.inc \
estimation/method_of_moments/RBC_MoM_steady_helper.m \
histval_initval_file_unit_tests.m \
histval_initval_file/my_assert.m \
histval_initval_file/ramst_data.xls \

253
tests/estimation/method_of_moments/AnScho_MoM.mod Normal file
View File

@ -0,0 +1,253 @@
% DSGE model used in replication files of
% An, Sungbae and Schorfheide, Frank, (2007), Bayesian Analysis of DSGE Models, Econometric Reviews, 26, issue 2-4, p. 113-172.
% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
% =========================================================================
% 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
@#define estimParams = 1
% Note that we set the numerical optimization tolerance levels very large to speed up the testsuite
@#define optimizer = 4
var c p R g y z INFL INT YGR;
varexo e_r e_g e_z;
parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst;
varobs INT YGR INFL;
tau = 2;
nu = 0.1;
kap = 0.33;
cyst = 0.85;
psi1 = 1.5;
psi2 = 0.125;
rhor = 0.75;
rhog = 0.95;
rhoz = 0.9;
rrst = 1;
pist = 3.2;
gamst = 0.55;
model;
#pist2 = exp(pist/400);
#rrst2 = exp(rrst/400);
#bet = 1/rrst2;
#phi = tau*(1-nu)/nu/kap/pist2^2;
#gst = 1/cyst;
#cst = (1-nu)^(1/tau);
#yst = cst*gst;
#dy = y-y(-1);
1 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1));
(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+y(+1)-y+p(+1));
exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2;
R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(dy+z) + e_r/100;
g = rhog*g(-1) + e_g/100;
z = rhoz*z(-1) + e_z/100;
YGR = gamst+100*(dy+z);
INFL = pist+400*p;
INT = pist+rrst+4*gamst+400*R;
end;
steady_state_model;
z = 0; p = 0; g = 0; r = 0; c = 0; y = 0;
YGR = gamst; INFL = pist; INT = pist + rrst + 4*gamst;
end;
shocks;
var e_r = 0.20^2;
var e_g = 0.80^2;
var e_z = 0.45^2;
corr e_r,e_g = 0.2;
end;
@#if estimParams == 0
% Define only initial values without bounds
estimated_params;
%tau, 1.50;
%kap, 0.15;
psi1, 1.20;
psi2, 0.50;
rhor, 0.50;
%rhog, 0.50;
%rhoz, 0.50;
%rrst, 1.20;
%pist, 3.00;
gamst, 0.75;
stderr e_r, 0.30;
stderr e_g, 0.30;
stderr e_z, 0.30;
corr e_r,e_g, 0.10;
end;
@#endif
@#if estimParams == 1
% Define initial values and bounds
estimated_params;
%tau, 1.50, 1e-5, 10;
%kap, 0.15, 1e-5, 10;
psi1, 1.20, 1e-5, 10;
psi2, 0.50, 1e-5, 10;
rhor, 0.50, 1e-5, 0.99999;
%rhog, 0.50, 1e-5, 0.99999;
%rhoz, 0.50, 1e-5, 0.99999;
%rrst, 1.20, 1e-5, 10;
%pist, 3.00, 1e-5, 20;
gamst, 0.75, -5, 5;
stderr e_r, 0.30, 1e-8, 5;
stderr e_g, 0.30, 1e-8, 5;
stderr e_z, 0.30, 1e-8, 5;
corr e_r,e_g, 0.10, -1, 1;
end;
@#endif
@#if estimParams == 2
% Define prior distribution
estimated_params;
%tau, 1.50, 1e-5, 10, gamma_pdf, 2.00, 0.50;
%kap, 0.15, 1e-5, 10, gamma_pdf, 0.33, 0.10;
psi1, 1.20, 1e-5, 10, gamma_pdf, 1.50, 0.25;
psi2, 0.50, 1e-5, 10, gamma_pdf, 0.125, 0.25;
rhor, 0.50, 1e-5, 0.99999, beta_pdf, 0.50, 0.20;
%rhog, 0.50, 1e-5, 0.99999, beta_pdf, 0.80, 0.10;
%rhoz, 0.50, 1e-5, 0.99999, beta_pdf, 0.66, 0.15;
%rrst, 1.20, 1e-5, 10, gamma_pdf, 0.50, 0.50;
%pist, 3.00, 1e-5, 20, gamma_pdf, 7.00, 2.00;
gamst, 0.75, -5, 5, normal_pdf, 0.40, 0.20;
stderr e_r, 0.30, 1e-8, 5, inv_gamma_pdf, 0.50, 0.26;
stderr e_g, 0.30, 1e-8, 5, inv_gamma_pdf, 1.25, 0.65;
stderr e_z, 0.30, 1e-8, 5, inv_gamma_pdf, 0.63, 0.33;
corr e_r,e_g, 0.10, -1, 1, uniform_pdf, , , -1, 1;
end;
@#endif
% Simulate data
stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=750,drop=500);
save('AnScho_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
iYGR = strmatch('YGR', M_.endo_names,'exact');
iINFL = strmatch('INFL', M_.endo_names,'exact');
iINT = strmatch('INT', M_.endo_names,'exact');
% first entry: number of variable in declaration order
% second entry: lag
% third entry: power
matched_moments_ = {
%first-order product moments
[iYGR ] [0 ], [1 ];
[iINFL ] [0 ], [1 ];
[iINT ] [0 ], [1 ];
%second-order contemporenous product moments
[iYGR iYGR ] [0 0], [1 1];
[iYGR iINFL] [0 0], [1 1];
[iYGR iINT ] [0 0], [1 1];
[iINFL iINFL] [0 0], [1 1];
[iINFL iINT ] [0 0], [1 1];
[iINT iINT ] [0 0], [1 1];
%second-order temporal product moments
[iYGR iYGR ] [0 -1], [1 1];
%[iINT iYGR ] [0 -1], [1 1];
%[iINFL iYGR ] [0 -1], [1 1];
%[iYGR iINT ] [0 -1], [1 1];
[iINT iINT ] [0 -1], [1 1];
%[iINFL iINT ] [0 -1], [1 1];
%[iYGR iINFL] [0 -1], [1 1];
%[iINT iINFL] [0 -1], [1 1];
[iINFL iINFL] [0 -1], [1 1];
};
@#for mommethod in ["GMM", "SMM"]
method_of_moments(
% Necessery options
mom_method = @{mommethod} % method of moments method; possible values: GMM|SMM
, datafile = 'AnScho_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 numbers for the iterations in the 2-step feasible method of moments
% , prefilter=0 % demean each data series by its empirical mean and use centered moments
%
% Options for SMM
% , 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 = 250 % 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-5
% ,'TolX', 1e-6
% ) % 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

BIN
tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat Normal file
View File

Binary file not shown.

202
tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod Normal file
View File

@ -0,0 +1,202 @@
% Tests SMM and GMM routines
%
% 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
@#define estimParams = 1
% 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 c iv n;
@#if estimParams == 0
estimated_params;
DELTA, 0.025;
BETTA, 0.984;
B, 0.5;
ETAc, 2;
ALFA, 0.667;
RHOA, 0.979;
stderr u_a, 0.0072;
end;
@#endif
@#if estimParams == 1
estimated_params;
DELTA, , 0, 1;
BETTA, , 0, 1;
B, , 0, 1;
ETAc, , 0, 10;
ALFA, , 0, 1;
RHOA, , 0, 1;
stderr u_a, , 0, 1;
end;
@#endif
@#if estimParams == 2
estimated_params;
DELTA, 0.025, 0, 1, normal_pdf, 0.02, 0.5;
BETTA, 0.98, 0, 1, beta_pdf, 0.90, 0.25;
B, 0.45, 0, 1, normal_pdf, 0.40, 0.5;
%ETAl, 1, 0, 10, normal_pdf, 0.25, 0.0.1;
ETAc, 1.8, 0, 10, normal_pdf, 1.80, 0.5;
ALFA, 0.65, 0, 1, normal_pdf, 0.60, 0.5;
RHOA, 0.95, 0, 1, normal_pdf, 0.90, 0.5;
stderr u_a, 0.01, 0, 1, normal_pdf, 0.01, 0.5;
%THETA, 3.48, 0, 10, normal_pdf, 0.25, 0.0.1;
end;
@#endif
% Simulate data
%stoch_simul(order=@{orderApp},pruning,nodisplay,nomoments,periods=500);
%save('RBC_MoM_data_@{orderApp}.mat', options_.varobs{:} );
%pause(1);
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 in ] [0 0], [1 1];
[iiv iiv] [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];
[ic ic ] [0 -3], [1 1];
[in in ] [0 -3], [1 1];
[iiv iiv] [0 -3], [1 1];
[ic ic ] [0 -5], [1 1];
[in in ] [0 -5], [1 1];
[iiv iiv] [0 -5], [1 1];
};
method_of_moments(
% Necessery options
mom_method = GMM % method of moments method; possible values: GMM|SMM
, datafile = 'RBC_Andreasen_Data_2.mat' % name of filename with data
% Options for both GMM and SMM
%, bartlett_kernel_lag = 20 % bandwith in optimal weighting matrix
, order = 2 % 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 = ['DIAGONAL','OPTIMAL'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
%, weighting_matrix_scaling_factor=1
, additional_optimizer_steps = [13] % 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
%, burnin = 200
%
% 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 = 50 % 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 = 13 % specifies the optimizer for minimization of moments distance, note that by default there is a new optimizer
, optim = ('TolFun', 1D-6
,'TolX', 1D-6
) % 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
, se_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
, mode_check
%, mode_check_neighbourhood_size=0.5
%, mode_check_symmetric_plots=0
%, mode_check_number_of_points=25
);

191
tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod Normal file
View File

@ -0,0 +1,191 @@
%
% 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 = 1
@#define estimParams = 0
% Note that we will set the numerical optimization tolerance levels very large to speed up the testsuite
@#define optimizer = 5
@#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.025;
BETTA, 0.984;
B, 0.5;
%ETAl, 1;
ETAc, 1;
ALFA, 0.667;
RHOA, 0.979;
stderr u_a, 0.0072;
%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=250);
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 = ['identity_matrix'] % weighting matrix in moments distance objective function; possible values: OPTIMAL|IDENTITY_MATRIX|DIAGONAL|filename
, weighting_matrix_scaling_factor = 10
, burnin=250
%, 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

80
tests/estimation/method_of_moments/RBC_MoM_common.inc Normal file
View File

@ -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;

161
tests/estimation/method_of_moments/RBC_MoM_prefilter.mod Normal file
View File

@ -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(8,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

8
tests/estimation/method_of_moments/RBC_MoM_steady_helper.m Normal file
View File

@ -0,0 +1,8 @@
function N = RBC_MoM_steady_helper(THETA,ETAl,ETAc,BETTA,B,C_O_N,W)
if ETAc == 1 && ETAl == 1
N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
else
% No closed-form solution use a fixed-point algorithm
N0 = 1/3;
N = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,optimset('Display','off','TolX',1e-12,'TolFun',1e-12));
end

74
tests/estimation/method_of_moments/RBC_MoM_steadystate.m Normal file
View File

@ -0,0 +1,74 @@
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
function [ys,params,check] = RBCmodel_steadystate(ys,exo,M_,options_)
%% Step 0: initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display','off','TolX',1e-12,'TolFun',1e-12);
params = M_.params;
%% Step 1: read out parameters to access them with their name
for ii = 1:M_.param_nbr
eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end
%% Step 2: Check parameter restrictions
if ETAc*ETAl<1 % parameter violates restriction (here it is artifical)
check=1; %set failure indicator
return; %return without updating steady states
end
%% Step 3: Enter model equations here
A = 1;
RK = 1/BETTA - (1-DELTA);
K_O_N = (RK/(A*(1-ALFA)))^(-1/ALFA);
if K_O_N <= 0
check = 1; % set failure indicator
return; % return without updating steady states
end
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;
if C_O_N <= 0
check = 1; % set failure indicator
return; % return without updating steady states
end
% The labor level
if ETAc == 1 && ETAl == 1
N = (1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA/(1+(1-BETTA*B)*(C_O_N*(1-B))^-1*W/THETA);
else
% No closed-form solution use a fixed-point algorithm
N0 = 1/3;
[N,~,exitflag] = fsolve(@(N) THETA*(1-N)^(-ETAl)*N^ETAc - (1-BETTA*B)*(C_O_N*(1-B))^(-ETAc)*W, N0,options);
if exitflag <= 0
check = 1; % set failure indicator
return % return without updating steady states
end
end
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);
%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end
for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end
end