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