diff --git a/matlab/+estimate/nls.m b/matlab/+estimate/nls.m
index 9b1cec30d..63a6e2e77 100644
--- a/matlab/+estimate/nls.m
+++ b/matlab/+estimate/nls.m
@@ -355,49 +355,11 @@ end
 
 % Set options if provided as input arguments to nls routine.
 oldopt = options_.optim_opt;
-if nargin>5
-    if mod(nargin-5, 2)
-        error('Options must come by key/value pairs.')
-    end
-    i = 1;
-    while i<nargin-5
-        if isequal(varargin{i}, 'noprint')
-            noprint = varargin{i+1};
-            i = i+2;
-            continue
-        else
-            if ~exist('opt', 'var')
-                opt = sprintf('''%s''', varargin{i});
-            else
-                opt = sprintf('%s,''%s''', opt, varargin{i});
-            end
-            if isnumeric(varargin{i+1})
-                opt = sprintf('%s,%s', opt, num2str(varargin{i+1}));
-            else
-                opt = sprintf('%s,''%s''', opt, varargin{i+1});
-            end
-            i = i+2;
-        end
-        options_.optim_opt = opt;
-    end
-end
-
-if ~exist('opt', 'var')
-    options_.optim_opt = [];
-end
-
-if ~exist('noprint', 'var')
-    noprint = false;
-end
-
-if nargin<5
-    % If default optimization algorithm is used (csminwel), do not print
-    % iterations.
-    options_.optim_opt = '''verbosity'',0';
-end
+[noprint, opt] = opt4nls(varargin);
+options_.optim_opt = opt;
 
 %
-% Check that we are able to evaluate the Sum of Squared Residuals on the initial guess
+% Check that we are able to evaluate the Sum of uared Residuals on the initial guess
 %
 
 ssr0 = ssrfun(params0);
diff --git a/matlab/+pac/+estimate/nls.m b/matlab/+pac/+estimate/nls.m
index caac195cb..0b3d409de 100644
--- a/matlab/+pac/+estimate/nls.m
+++ b/matlab/+pac/+estimate/nls.m
@@ -75,6 +75,8 @@ end
 [pacmodl, lhs, rhs, pnames, enames, xnames, ~, pid, eid, xid, ~, ipnames_, params, data, islaggedvariables] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
+dLHS = M_.aux_vars(strmatch(lhs,M_.endo_names, 'exact')==[M_.aux_vars(:).endo_index]).orig_expr;
+
 % Check that the error correction term is correct.
 if M_.pac.(pacmodl).ec.istarget(2)
     error(['\nThe error correction term in PAC equation (%s) is not correct.\nThe ' ...
@@ -261,34 +263,8 @@ end
 
 % Set options if provided as input arguments to nls routine.
 oldopt = options_.optim_opt;
-if nargin>5
-    if mod(nargin-5, 2)
-        error('Options must come by key/value pairs.')
-    end
-    i = 1;
-    opt = '';
-    while i<nargin-5
-        if i==1
-            opt = sprintf('''%s''', varargin{i});
-        else
-            opt = sprintf('%s,''%s''', opt, varargin{i});
-        end
-        if isnumeric(varargin{i+1})
-            opt = sprintf('%s,%s', opt, num2str(varargin{i+1}));
-        else
-            opt = sprintf('%s,''%s''', opt, varargin{i+1});
-        end
-        i = i+2;
-    end
-    options_.optim_opt = opt;
-else
-    options_.optim_opt = [];
-end
-if nargin<5
-    % If default optimization algorithm is used (csminwel), do not print
-    % iterations.
-    options_.optim_opt = '''verbosity'',0';
-end
+[noprint, opt] = opt4nls(varargin);
+options_.optim_opt = opt;
 
 %
 % Check that we are able to evaluate the Sum of Squared Residuals on the initial guess
@@ -344,6 +320,7 @@ oo_.pac.(pacmodl).lhs = dseries(lhs, range(1), 'lhs');
 oo_.pac.(pacmodl).fit = dseries(lhs-r, range(1), 'fit');
 oo_.pac.(pacmodl).residual = dseries(r, range(1), 'residual');
 oo_.pac.(pacmodl).ssr = SSR;
+oo_.pac.(pacmodl).s2 = SSR/T;
 oo_.pac.(pacmodl).R2 = 1-var(r)/var(lhs);
 oo_.pac.(pacmodl).parnames = fieldnames(params);
 oo_.pac.(pacmodl).estimator = params1;
@@ -353,3 +330,16 @@ oo_.pac.(pacmodl).student = params1./(sqrt(diag(C)));
 % Also save estimated parameters in M_
 M_.params(ipnames_) = params1;
 M_ = pac.update.parameters(pacmodl, M_, oo_, false);
+
+if ~noprint
+    title = ['NLS Estimation of equation ''' eqname ''''];
+    preamble = {['Dependent Variable: ' dLHS], ...
+                sprintf('Observations: %d from %s to %s\n', (range(end)-range(1))+1, range(1).char, range(end).char)};
+
+    afterward = {sprintf('R^2: %f', oo_.pac.(pacmodl).R2), ...
+                 sprintf('s^2: %f', oo_.pac.(pacmodl).s2)}; ...
+
+    dyn_table(title, preamble, afterward, oo_.pac.(pacmodl).parnames, ...
+              {'Estimates','t-statistic','Std. Error'}, 4, ...
+              [oo_.pac.(pacmodl).estimator oo_.pac.(pacmodl).student sqrt(diag(C))]);
+end
\ No newline at end of file
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index ad68e100e..7e605be6a 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -71,6 +71,7 @@ p = {'/distributions/' ; ...
      '/utilities/dataset/' ; ...
      '/utilities/general/' ; ...
      '/utilities/graphics/' ; ...
+     '/utilities/estimation/' ; ...
      '/modules/reporting/src/' ; ...
      '/modules/reporting/macros/'};
 
diff --git a/matlab/utilities/estimation/opt4nls.m b/matlab/utilities/estimation/opt4nls.m
new file mode 100644
index 000000000..e6067e578
--- /dev/null
+++ b/matlab/utilities/estimation/opt4nls.m
@@ -0,0 +1,53 @@
+function [noprint, opt] = opt4nls(varargin)
+
+% Sets options for NLS routines.
+
+% Copyright © 021 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/>.
+
+varargin = varargin{1};
+nargin = length(varargin);
+
+noprint = false;
+
+if nargin
+    if mod(nargin, 2)
+        error('Options must come by key/value pairs.')
+    end
+    i = 1;
+    while i<nargin
+        if isequal(varargin{i}, 'noprint')
+            noprint = varargin{i+1};
+            i = i+2;
+            continue
+        else
+            if ~exist('opt', 'var')
+                opt = sprintf('''%s''', varargin{i});
+            else
+                opt = sprintf('%s,''%s''', opt, varargin{i});
+            end
+            if isnumeric(varargin{i+1})
+                opt = sprintf('%s,%s', opt, num2str(varargin{i+1}));
+            else
+                opt = sprintf('%s,''%s''', opt, varargin{i+1});
+            end
+            i = i+2;
+        end
+    end
+else
+    opt = '';
+end
\ No newline at end of file
diff --git a/tests/pac/var-10e/example1.mod b/tests/pac/var-10e/example1.mod
index b185d7345..c7e8fe0e6 100644
--- a/tests/pac/var-10e/example1.mod
+++ b/tests/pac/var-10e/example1.mod
@@ -96,8 +96,4 @@ e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
 c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
 c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
 
-disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
-disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
-disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
-
 save('example1.estimation.mat', 'e_c_m_nls', 'c_z_1_nls', 'c_z_2_nls')
diff --git a/tests/pac/var-10e/example2.mod b/tests/pac/var-10e/example2.mod
index ddc56ad8c..fdce3ec6a 100644
--- a/tests/pac/var-10e/example2.mod
+++ b/tests/pac/var-10e/example2.mod
@@ -80,10 +80,6 @@ e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
 c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
 c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
 
-disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
-disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
-disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
-
 % Check consistency with disaggregated target
 ts1 = dseries('example1.data.mat');
 ts2 = dseries('example2.data.mat');
diff --git a/tests/pac/var-11e/example1.mod b/tests/pac/var-11e/example1.mod
index 89a56773d..13245257b 100644
--- a/tests/pac/var-11e/example1.mod
+++ b/tests/pac/var-11e/example1.mod
@@ -90,11 +90,3 @@ eparams.c_z_2   =  .2;
 edata = TrueData;
 edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
 pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
-
-e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
-c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
-c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
-
-disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
-disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
-disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
diff --git a/tests/pac/var-8e/example1.mod b/tests/pac/var-8e/example1.mod
index 9fb1bbc9b..95bd61dc1 100644
--- a/tests/pac/var-8e/example1.mod
+++ b/tests/pac/var-8e/example1.mod
@@ -84,11 +84,3 @@ eparams.c_z_2   =  .2;
 edata = TrueData;
 edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
 pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
-
-e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
-c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
-c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
-
-disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
-disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
-disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))
diff --git a/tests/pac/var-9e/example1.mod b/tests/pac/var-9e/example1.mod
index 4d0169f33..6b3dbc6ac 100644
--- a/tests/pac/var-9e/example1.mod
+++ b/tests/pac/var-9e/example1.mod
@@ -88,11 +88,3 @@ eparams.c_z_2   =  .2;
 edata = TrueData;
 edata.ez = dseries(NaN(TrueData.nobs, 1), 2000Q1, 'ez');
 pac.estimate.nls('zpac', eparams, edata, 2005Q1:2005Q1+4000, 'fmincon');
-
-e_c_m_nls = M_.params(strmatch('e_c_m', M_.param_names, 'exact'));
-c_z_1_nls = M_.params(strmatch('c_z_1', M_.param_names, 'exact'));
-c_z_2_nls = M_.params(strmatch('c_z_2', M_.param_names, 'exact'));
-
-disp(sprintf('Estimate of e_c_m: %f', e_c_m_nls))
-disp(sprintf('Estimate of c_z_1: %f', c_z_1_nls))
-disp(sprintf('Estimate of c_z_2: %f', c_z_2_nls))