diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index ee74ebf11a9303f36a1e7f59993061487a87f3a5..4109f08dc9fdfeaf25eed52780961edc12b4c2ea 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -61,6 +61,7 @@ p = {'/distributions/' ; ...
'/cli/' ; ...
'/lmmcp/' ; ...
'/optimization/' ; ...
+ '/method_of_moments/' ; ...
'/discretionary_policy/' ; ...
'/accessors/' ; ...
'/modules/dseries/src/' ; ...
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index c893be9c63ded305bbcf704a5a657a659ff91d21..27cc95b34de93c8f74799e04293dc5938a429e2b 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -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')
diff --git a/matlab/get_matrix_entries_for_psd_check.m b/matlab/get_matrix_entries_for_psd_check.m
new file mode 100644
index 0000000000000000000000000000000000000000..1ddd80d3c53f4f032d833b87724733c9d4a9cfba
--- /dev/null
+++ b/matlab/get_matrix_entries_for_psd_check.m
@@ -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
\ No newline at end of file
diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
new file mode 100644
index 0000000000000000000000000000000000000000..faa832f567ca50f47e59568add9f4c97d49149ce
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -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
diff --git a/matlab/method_of_moments/method_of_moments_data_moments.m b/matlab/method_of_moments/method_of_moments_data_moments.m
new file mode 100644
index 0000000000000000000000000000000000000000..38d205a272f7e6155e52b81e4c9aa3c7a067a24d
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments_data_moments.m
@@ -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
+
+
+
diff --git a/matlab/method_of_moments/method_of_moments_mode_check.m b/matlab/method_of_moments/method_of_moments_mode_check.m
new file mode 100644
index 0000000000000000000000000000000000000000..1df8af84f11105a0431bd222d5c70a0662539b4e
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments_mode_check.m
@@ -0,0 +1,185 @@
+function method_of_moments_mode_check(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
+% Checks the estimated ML mode or Posterior mode.
+
+
+% Copyright (C) 2020 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <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');
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/method_of_moments/method_of_moments_objective_function.m
new file mode 100644
index 0000000000000000000000000000000000000000..e4f6714aaf91e6ae6b9cbe1c22d0756453c7d850
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -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
+
diff --git a/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
new file mode 100644
index 0000000000000000000000000000000000000000..7dde93568d45b011e081e63b615c5ffbcb9749be
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
@@ -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
\ No newline at end of file
diff --git a/matlab/method_of_moments/method_of_moments_standard_errors.m b/matlab/method_of_moments/method_of_moments_standard_errors.m
new file mode 100644
index 0000000000000000000000000000000000000000..785f4e8a29ab75c227e711f2690295c4d11e99e9
--- /dev/null
+++ b/matlab/method_of_moments/method_of_moments_standard_errors.m
@@ -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));
\ No newline at end of file
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index a7e6a83f5dc7f4158d5dd5e5e0e896d2d97fd23a..475392f593fbd46a1a0800af847ebd056544f263 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -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)
diff --git a/matlab/plot_priors.m b/matlab/plot_priors.m
index 02d7826688b28ca9c0d071f121d652d49b9890d2..7d8aec8f7c3616807da5d229c4d1d38f92682dfd 100644
--- a/matlab/plot_priors.m
+++ b/matlab/plot_priors.m
@@ -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);
diff --git a/tests/.gitignore b/tests/.gitignore
index 501fbc71d03f0f5e8098101e43bfa26287473633..722a27d5d7d11d53c5cdc860106a4364cda2f7c0 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -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
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 3e5575d906f70e29ebfbe0664887831c50fb5867..88302799d47117c221df3fed009779808e20f563 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -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 \
diff --git a/tests/estimation/method_of_moments/AnScho_MoM.mod b/tests/estimation/method_of_moments/AnScho_MoM.mod
new file mode 100644
index 0000000000000000000000000000000000000000..674111120283a75c0f3ac4a07b12555248d32bca
--- /dev/null
+++ b/tests/estimation/method_of_moments/AnScho_MoM.mod
@@ -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
+
diff --git a/tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat b/tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0b2ba62defdaab77aa663f2907ae16801837f6b6
Binary files /dev/null and b/tests/estimation/method_of_moments/RBC_Andreasen_Data_2.mat differ
diff --git a/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod b/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod
new file mode 100644
index 0000000000000000000000000000000000000000..feb47cb3a7db966a98538d123d5eaa782161c06b
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_Andreasen.mod
@@ -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
+ );
+
+
+
diff --git a/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod b/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a0e4ea654cd699c3db19dffc5a75645ba02f3cab
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_SMM_ME.mod
@@ -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
+
+
+
diff --git a/tests/estimation/method_of_moments/RBC_MoM_common.inc b/tests/estimation/method_of_moments/RBC_MoM_common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..625a9c50b88616ac61628a67eec2f587a1e2e461
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_common.inc
@@ -0,0 +1,80 @@
+% RBC model used in replication files of
+% Andreasen, Fernández-Villaverde, Rubio-Ramírez (2018): "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", Review of Economic Studies, 85(1):1-49.
+% Adapted by Willi Mutschler (@wmutschl, willi@mutschler.eu)
+% =========================================================================
+% Copyright (C) 2020 Dynare Team
+
+var k $K$
+ c $C$
+ a $A$
+ iv $I$
+ y $Y$
+ la $\lambda$
+ n $N$
+ rk ${r^k}$
+ w $W$
+ ;
+
+predetermined_variables k;
+
+varexo u_a ${\varepsilon^{a}}$
+ ;
+
+parameters DELTA $\delta$
+ BETTA $\beta$
+ B $B$
+ ETAl $\eta_l$
+ ETAc $\eta_c$
+ THETA $\theta$
+ ALFA $\alpha$
+ RHOA $\rho_a$
+ ;
+
+DELTA = 0.025;
+BETTA = 0.984;
+B = 0.5;
+ETAl = 1;
+ETAc = 2;
+THETA = 3.48;
+ALFA = 0.667;
+RHOA = 0.979;
+
+model;
+0 = -exp(la) +(exp(c)-B*exp(c(-1)))^(-ETAc) - BETTA*B*(exp(c(+1))-B*exp(c))^(-ETAc);
+0 = -THETA*(1-exp(n))^-ETAl + exp(la)*exp(w);
+0 = -exp(la) + BETTA*exp(la(+1))*(exp(rk(+1)) + (1-DELTA));
+0 = -exp(a)*(1-ALFA)*exp(k)^(-ALFA)*exp(n)^(ALFA) + exp(rk);
+0 = -exp(a)*ALFA*exp(k)^(1-ALFA)*exp(n)^(ALFA-1) + exp(w);
+0 = -exp(c) - exp(iv) + exp(y);
+0 = -exp(y) + exp(a)*exp(k)^(1-ALFA)*exp(n)^(ALFA);
+0 = -exp(k(+1)) + (1-DELTA)*exp(k) + exp(iv);
+0 = -log(exp(a)) + RHOA*log(exp(a(-1))) + u_a;
+end;
+
+steady_state_model;
+A = 1;
+RK = 1/BETTA - (1-DELTA);
+K_O_N = (RK/(A*(1-ALFA)))^(-1/ALFA);
+W = A*ALFA*(K_O_N)^(1-ALFA);
+IV_O_N = DELTA*K_O_N;
+Y_O_N = A*K_O_N^(1-ALFA);
+C_O_N = Y_O_N - IV_O_N;
+
+N=RBC_MoM_steady_helper(THETA,ETAl,ETAc,BETTA,B,C_O_N,W);
+C=C_O_N*N;
+Y=Y_O_N*N;
+IV=IV_O_N*N;
+K=K_O_N*N;
+LA = (C-B*C)^(-ETAc)-BETTA*B*(C-B*C)^(-ETAc);
+
+k=log(K);
+c=log(C);
+a=log(A);
+iv=log(IV);
+y=log(Y);
+la=log(LA);
+n=log(N);
+rk=log(RK);
+w=log(W);
+
+end;
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod b/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod
new file mode 100644
index 0000000000000000000000000000000000000000..22924d0666800b38b523efa406562c68a75dc02b
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_prefilter.mod
@@ -0,0 +1,161 @@
+% Tests SMM and GMM routines with prefilter, explicit initialization, and estimated_params_init(use_calibration);
+%
+% Copyright (C) 2020 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <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
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m b/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m
new file mode 100644
index 0000000000000000000000000000000000000000..08185c1e16ef38dedd379c6e46d726df66b545fe
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_steady_helper.m
@@ -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
\ No newline at end of file
diff --git a/tests/estimation/method_of_moments/RBC_MoM_steadystate.m b/tests/estimation/method_of_moments/RBC_MoM_steadystate.m
new file mode 100644
index 0000000000000000000000000000000000000000..ba4ef9240b522ac5bbf83fc41380a50f3f1805f8
--- /dev/null
+++ b/tests/estimation/method_of_moments/RBC_MoM_steadystate.m
@@ -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