diff --git a/matlab/+mom/Jtest.m b/matlab/+mom/Jtest.m
new file mode 100644
index 0000000000000000000000000000000000000000..0cf41df5dd8b51189662d5e53fd78a21af20a933
--- /dev/null
+++ b/matlab/+mom/Jtest.m
@@ -0,0 +1,67 @@
+function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
+% function oo_ = Jtest(xparam, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, nobs)
+% -------------------------------------------------------------------------
+% Computes the J-test statistic and p-value for a GMM/SMM estimation
+% =========================================================================
+% INPUTS
+%  xparam:              [vector]           estimated parameter vector
+%  objective_function:  [function handle]  objective function
+%  Woptflag:            [logical]          flag if optimal weighting matrix has already been computed
+%  oo_:                 [struct]           results
+%  options_mom_:        [struct]           options
+%  bayestopt_:          [struct]           information on priors
+%  Bounds:              [struct]           bounds on parameters
+%  estim_params_:       [struct]           information on estimated parameters
+%  M_:                  [struct]           information on the model
+%  nobs:                [scalar]           number of observations
+% -------------------------------------------------------------------------
+% OUTPUT
+%  oo_:                 [struct]           updated results
+% -------------------------------------------------------------------------
+% This function is called by
+%  o mom.run
+% -------------------------------------------------------------------------
+% This function calls
+%  o mom.objective_function
+%  o mom.optimal_weighting_matrix
+% =========================================================================
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+% =========================================================================
+
+if options_mom_.mom.mom_nbr > length(xparam)
+    % Get optimal weighting matrix for J test, if necessary
+    if ~Woptflag
+        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, xparam, Bounds, oo_J, estim_params_, M_, options_mom_, bayestopt_);
+    else
+        fval = oo_.mom.Q;
+    end
+    % Compute J statistic
+    if strcmp(options_mom_.mom.mom_method,'SMM')
+        Variance_correction_factor = options_mom_.mom.variance_correction_factor;
+    elseif strcmp(options_mom_.mom.mom_method,'GMM')
+        Variance_correction_factor = 1;
+    end
+    oo_.mom.J_test.j_stat          = nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
+    oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam);
+    oo_.mom.J_test.p_val           = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
+    fprintf('\nValue of J-test statistic: %f\n',oo_.mom.J_test.j_stat);
+    fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val);
+end
\ No newline at end of file
diff --git a/matlab/+mom/run.m b/matlab/+mom/run.m
index 6e8e81994f2a20f418e4feff1032ddbf7ff85eda..846bcea13e32489365dff58023ded5e5519a3baa 100644
--- a/matlab/+mom/run.m
+++ b/matlab/+mom/run.m
@@ -569,39 +569,14 @@ if strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_meth
 end
 
 
-% -------------------------------------------------------------------------
-% Step 8: J test
-% -------------------------------------------------------------------------
-if options_mom_.mom.mom_nbr > length(xparam1)
-    %get optimal weighting matrix for J test, if necessary
-    if ~Woptflag
-        W_opt = 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_);
-    end
-
-    % Compute J statistic
-    if strcmp(options_mom_.mom.mom_method,'SMM')    
-        Variance_correction_factor = options_mom_.mom.variance_correction_factor;
-    elseif strcmp(options_mom_.mom.mom_method,'GMM')
-        Variance_correction_factor=1;
-    end
-    oo_.mom.J_test.j_stat          = dataset_.nobs*Variance_correction_factor*fval/options_mom_.mom.weighting_matrix_scaling_factor;
-    oo_.mom.J_test.degrees_freedom = length(oo_.mom.model_moments)-length(xparam1);
-    oo_.mom.J_test.p_val           = 1-chi2cdf(oo_.mom.J_test.j_stat, oo_.mom.J_test.degrees_freedom);
-    fprintf('\nvalue of J-test statistic: %f\n',oo_.mom.J_test.j_stat)
-    fprintf('p-value of J-test statistic: %f\n',oo_.mom.J_test.p_val)
-end
-
-
 % -------------------------------------------------------------------------
 % display final estimation results
 % -------------------------------------------------------------------------
 if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
     % Store results in output structure
     oo_.mom = display_estimation_results_table(xparam1,SE,M_,options_mom_,estim_params_,bayestopt_,oo_.mom,prior_dist_names,options_mom_.mom.mom_method,lower(options_mom_.mom.mom_method));
-
+    % J test
+    oo_ = mom.Jtest(xparam1, objective_function, Woptflag, oo_, options_mom_, bayestopt_, Bounds, estim_params_, M_, dataset_.nobs);
 title = ['Comparison of data moments and model moments (',options_mom_.mom.mom_method,')'];
 headers = {'Moment','Data','Model'};
 for jm = 1:size(M_.matched_moments,1)