diff --git a/matlab/method_of_moments/method_of_moments_check_plot.m b/matlab/+mom/check_plot.m
similarity index 94%
rename from matlab/method_of_moments/method_of_moments_check_plot.m
rename to matlab/+mom/check_plot.m
index 44ab05f465b7b6f489db50ae397c388fec581b94..fe4ddadb548f779ad995614782fa1f2acb30d0c3 100644
--- a/matlab/method_of_moments/method_of_moments_check_plot.m
+++ b/matlab/+mom/check_plot.m
@@ -1,4 +1,5 @@
-function method_of_moments_check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
+function check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
+%function check_plot(fun,xparam,SE_vec,options_,M_,estim_params_,Bounds,bayestopt_,varargin)
 % Checks the estimated local minimum of the moment's distance objective
 
 
@@ -44,7 +45,7 @@ if ~exist([M_.dname filesep 'graphs'],'dir')
 end
 if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fidTeX = fopen([M_.dname, '/graphs/', M_.fname '_MoMCheckPlots.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by method_of_moments_check_plot.m (Dynare).\n');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by mom.check_plot.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
     fprintf(fidTeX,' \n');
 end
@@ -119,7 +120,7 @@ for plt = 1:nbplt
             else
                 y(i,1) = NaN;
                 if options_.debug
-                    fprintf('method_of_moments_check_plot:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
+                    fprintf('mom.check_plot:: 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
diff --git a/matlab/method_of_moments/method_of_moments_data_moments.m b/matlab/+mom/data_moments.m
similarity index 91%
rename from matlab/method_of_moments/method_of_moments_data_moments.m
rename to matlab/+mom/data_moments.m
index cd145df2b74170cfd64a6899c90db505cdb372d8..328205f51839ae88ba399fed26e5718433ebe454 100644
--- a/matlab/method_of_moments/method_of_moments_data_moments.m
+++ b/matlab/+mom/data_moments.m
@@ -1,5 +1,5 @@
-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_)
+function [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_)
+% [dataMoments, m_data] = data_moments(data, oo_, matched_moments_, options_mom_)
 % This function computes the user-selected empirical moments from data
 % =========================================================================
 % INPUTS
@@ -13,8 +13,8 @@ function [dataMoments, m_data] = method_of_moments_data_moments(data, oo_, match
 %  o m_data                  [T x numMom]       selected empirical moments at each point in time
 % -------------------------------------------------------------------------
 % This function is called by
-%  o method_of_moments.m
-%  o method_of_moments_objective_function.m
+%  o mom.run.m
+%  o mom.objective_function.m
 % =========================================================================
 % Copyright (C) 2020-2021 Dynare Team
 %
@@ -69,7 +69,4 @@ for jm = 1:options_mom_.mom.mom_nbr
 end
 
 
-end %function end
-
-
-
+end %function end
\ No newline at end of file
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/+mom/objective_function.m
similarity index 97%
rename from matlab/method_of_moments/method_of_moments_objective_function.m
rename to matlab/+mom/objective_function.m
index 27af647a5c3a459de7817f8a8111660c525d7cb8..d0f9672b0b3c6db392ef61e6835796fa82039dc7 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/+mom/objective_function.m
@@ -1,5 +1,5 @@
-function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
 % -------------------------------------------------------------------------
 % This function evaluates the objective function for GMM/SMM estimation
 % =========================================================================
@@ -27,7 +27,7 @@ function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_m
 %   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
 % -------------------------------------------------------------------------
 % This function is called by
-%  o method_of_moments.m
+%  o mom.run.m
 %  o dynare_minimize_objective.m
 % -------------------------------------------------------------------------
 % This function calls
@@ -250,7 +250,7 @@ elseif strcmp(options_mom_.mom.mom_method,'SMM')
     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_, M_.matched_moments, options_mom_);
+    oo_.mom.model_moments = mom.data_moments(y_sim, oo_, M_.matched_moments, options_mom_);
     
 end
 
diff --git a/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m b/matlab/+mom/optimal_weighting_matrix.m
similarity index 94%
rename from matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
rename to matlab/+mom/optimal_weighting_matrix.m
index d722d9a8ade7e52469f92feaf87cff2e2e24612f..4708db8f8f84aeaa0ccc41716068d0f846abaa12 100644
--- a/matlab/method_of_moments/method_of_moments_optimal_weighting_matrix.m
+++ b/matlab/+mom/optimal_weighting_matrix.m
@@ -1,5 +1,5 @@
-function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
-% W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_lag)
+function W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
+% W_opt = optimal_weighting_matrix(m_data, moments, q_lag)
 % -------------------------------------------------------------------------
 % This function computes the optimal weigthing matrix by a Bartlett kernel with maximum lag q_lag
 % Adapted from replication codes of
@@ -14,7 +14,7 @@ function W_opt = method_of_moments_optimal_weighting_matrix(m_data, moments, q_l
 %   o W_opt                  [numMom x numMom]  optimal weighting matrix
 % -------------------------------------------------------------------------
 % This function is called by
-%  o method_of_moments.m
+%  o mom.run.m
 % -------------------------------------------------------------------------
 % This function calls:
 %  o CorrMatrix (embedded)
diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/+mom/run.m
similarity index 97%
rename from matlab/method_of_moments/method_of_moments.m
rename to matlab/+mom/run.m
index 30d47f7fe533c0f661fb8c65ae5aad6ee97b6204..6cc44539e2f48d3c7e974fc007a4487196db184a 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/+mom/run.m
@@ -1,5 +1,5 @@
-function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
-%function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
+function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
+%function [oo_, options_mom_, M_] = run(bayestopt_, options_, oo_, estim_params_, M_, options_mom_)
 % -------------------------------------------------------------------------
 % This function performs a method of moments estimation with the following steps:
 %   Step 0: Check if required structures and options exist
@@ -49,11 +49,11 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
 %  o get_all_parameters.m
 %  o get_matrix_entries_for_psd_check.m
 %  o makedataset.m
-%  o method_of_moments_check_plot.m
-%  o method_of_moments_data_moments.m
-%  o method_of_moments_objective_function.m
-%  o method_of_moments_optimal_weighting_matrix
-%  o method_of_moments_standard_errors
+%  o mom.check_plot.m
+%  o mom.data_moments.m
+%  o mom.objective_function.m
+%  o mom.optimal_weighting_matrix
+%  o mom-standard_errors
 %  o plot_priors.m
 %  o print_info.m
 %  o prior_bounds.m
@@ -628,7 +628,7 @@ end
 fprintf('Computing data moments. Note that NaN values in the moments (due to leads and lags or missing data) are replaced by the mean of the corresponding moment\n');
 
 % Get data moments for the method of moments
-[oo_.mom.data_moments, oo_.mom.m_data] = method_of_moments_data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
+[oo_.mom.data_moments, oo_.mom.m_data] = mom.data_moments(dataset_.data, oo_, M_.matched_moments, options_mom_);
 
 % Get shock series for SMM and set variance correction factor
 if strcmp(options_mom_.mom.mom_method,'SMM')
@@ -715,7 +715,7 @@ end
 % -------------------------------------------------------------------------
 % Step 6: checks for objective function at initial parameters
 % -------------------------------------------------------------------------
-objective_function = str2func('method_of_moments_objective_function');
+objective_function = str2func('mom.objective_function');
 
 try
     % Check for NaN or complex values of moment-distance-funtion evaluated
@@ -851,19 +851,19 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
             fprintf('  - diagonal of optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
             if stage_iter == 1
                 fprintf('    and using data-moments as initial estimate of model-moments\n');
-                weighting_matrix = diag(diag(  method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag)  ));
+                weighting_matrix = diag(diag(  mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag)  ));
             else
                 fprintf('    and using previous stage estimate of model-moments\n');
-                weighting_matrix = diag(diag(  method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag)  ));
+                weighting_matrix = diag(diag(  mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag)  ));
             end            
         case 'optimal'
             fprintf('  - optimal weighting matrix (Bartlett kernel with %d lags)\n', options_mom_.mom.bartlett_kernel_lag);
             if stage_iter == 1
                 fprintf('    and using data-moments as initial estimate of model-moments\n');
-                weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
+                weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.data_moments, options_mom_.mom.bartlett_kernel_lag);
             else
                 fprintf('    and using previous stage estimate of model-moments\n');
-                weighting_matrix = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
+                weighting_matrix = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
                 Woptflag = true;
             end            
         otherwise %user specified matrix in file
@@ -922,7 +922,7 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
     [fval, ~, ~,~,~, oo_] = feval(objective_function, xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
     options_mom_.mom.compute_derivs = false; % reset to not compute derivatives in objective function during optimization
     
-    SE = method_of_moments_standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
+    SE = mom.standard_errors(xparam1, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Woptflag);
     
     % Store results in output structure
     oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_laplace,oo_.mom,prior_dist_names,sprintf('%s (STAGE %u)',options_mom_.mom.mom_method,stage_iter),sprintf('%s_stage_%u',lower(options_mom_.mom.mom_method),stage_iter));
@@ -934,7 +934,7 @@ end
 if options_mom_.mom.mom_nbr > length(xparam1)
     %get optimal weighting matrix for J test, if necessary
     if ~Woptflag
-        W_opt = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
+        W_opt = mom.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_, M_, options_mom_);
@@ -992,7 +992,7 @@ if options_mom_.TeX
 end
 
 if options_mom_.mode_check.status
-    method_of_moments_check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,...
+    mom.check_plot(objective_function,xparam1,SE,options_mom_,M_,estim_params_,Bounds,bayestopt_laplace,...
         Bounds, oo_, estim_params_, M_, options_mom_)
 end
 
diff --git a/matlab/method_of_moments/method_of_moments_standard_errors.m b/matlab/+mom/standard_errors.m
similarity index 88%
rename from matlab/method_of_moments/method_of_moments_standard_errors.m
rename to matlab/+mom/standard_errors.m
index f5de86162a35b20b261a5564fa6f9b43eaaa4cd9..43681a6340aa4b99958d577245b0958083a9795a 100644
--- a/matlab/method_of_moments/method_of_moments_standard_errors.m
+++ b/matlab/+mom/standard_errors.m
@@ -1,5 +1,5 @@
-function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
-% [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
+function [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
+% [SE_values, Asympt_Var] = standard_errors(xparam, objective_function, Bounds, oo_, estim_params_, M_, options_mom_, Wopt_flag)
 % -------------------------------------------------------------------------
 % This function computes standard errors to the method of moments estimates
 % Adapted from replication codes of
@@ -7,7 +7,7 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
 % =========================================================================
 % 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 objective_function        string of objective function
 %   o Bounds:                   structure containing parameter bounds
 %   o oo_:                      structure for results
 %   o estim_params_:            structure describing the estimated_parameters
@@ -20,13 +20,13 @@ function [SE_values, Asympt_Var] = method_of_moments_standard_errors(xparam, obj
 %   o Asympt_Var                 [nparam x nparam] asymptotic covariance matrix
 % -------------------------------------------------------------------------
 % This function is called by
-%  o method_of_moments.m
+%  o mom.run.m
 % -------------------------------------------------------------------------
 % This function calls:
 %  o get_the_name
 %  o get_error_message
-%  o method_of_moments_objective_function
-%  o method_of_moments_optimal_weighting_matrix  
+%  o mom.objective_function
+%  o mom.optimal_weighting_matrix  
 % =========================================================================
 % Copyright (C) 2020-2021 Dynare Team
 %
@@ -109,7 +109,7 @@ if Wopt_flag
     Asympt_Var  = 1/T*((D'*WW*D)\eye(dim_params));
 else
     % We do not have the optimal weighting matrix yet    
-    WWopt      = method_of_moments_optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
+    WWopt      = mom.optimal_weighting_matrix(oo_.mom.m_data, oo_.mom.model_moments, options_mom_.mom.bartlett_kernel_lag);
     S          = WWopt\eye(size(WWopt,1));
     AA         = (D'*WW*D)\eye(dim_params);
     Asympt_Var = 1/T*AA*D'*WW*S*WW*D*AA;
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 5aede4182bc0cf0bbabd28731a96dc464e8008de..ad68e100e33f536636a32f04f71d77cb63fc5642 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -63,7 +63,6 @@ p = {'/distributions/' ; ...
      '/optimization/' ; ...
      '/ols/' ; ...
      '/pac-tools/' ; ...
-     '/method_of_moments/' ; ...
      '/discretionary_policy/' ; ...
      '/accessors/' ; ...
      '/modules/dseries/src/' ; ...