diff --git a/matlab/+estimate/nls.m b/matlab/+estimate/nls.m
index c4a3e2a3e29f1983ce8ea1327184be1a53c35273..e73da463e810cb78f016fb1d6af5e35df18f6ce7 100644
--- a/matlab/+estimate/nls.m
+++ b/matlab/+estimate/nls.m
@@ -168,96 +168,16 @@ end
 % Rewrite and print the equation.
 %
 
-% List of objects to be replaced
-objNames = [pnames; enames; xnames];
-objIndex = [pid; eid; xid];
-objTypes = [ones(length(pid), 1); 2*ones(length(eid)+length(xid), 1)];
-
-[~,I] = sort(cellfun(@length, objNames), 'descend');
-objNames = objNames(I);
-objIndex = objIndex(I);
-objTypes = objTypes(I);
-
-% Substitute parameters and variables. Note that in the transformed model (we consider rhs
-% instead of RHS) the maximum lag is 1.
-for i=1:length(objNames)
-    switch objTypes(i)
-      case 1
-        rhs = strrep(rhs, objNames{i}, sprintf('DynareModel.params(%u)', objIndex(i)));
-      case 2
-        k = find(strcmp(objNames{i}, data.name));
-        if isempty(k)
-            error('Variable %s is missing in the database.', objNames{i})
-        end
-        j = regexp(rhs, ['\<', objNames{i}, '\>']);
-        if islaggedvariables
-            jlag = regexp(rhs, ['\<', objNames{i}, '\(-1\)']);
-            if ~isempty(jlag)
-                rhs = regexprep(rhs, ['\<' objNames{i} '\(-1\)'], sprintf('data(1:end-1,%u)', k));
-            end
-            if ~isempty(setdiff(j, jlag))
-                rhs = regexprep(rhs, ['\<' objNames{i} '\>'], sprintf('data(2:end,%u)', k));
-            end
-        else
-            rhs = regexprep(rhs, ['\<' objNames{i} '\>'], sprintf('data(:,%u)', k));
-        end
-        if contains(lhs, objNames{i})
-            if islaggedvariables
-                lhs = strrep(lhs, objNames{i}, sprintf('data(2:end,%u)', k));
-            else
-                lhs = strrep(lhs, objNames{i}, sprintf('data(:,%u)', k));
-            end
-        end
-    end
-end
-
-% Allow elementwise operations
-rhs = strrep(rhs, '^', '.^');
-rhs = strrep(rhs, '/', './');
-rhs = strrep(rhs, '*', '.*');
+[rhs, lhs] = rewrite_equation_with_tables(rhs, lhs, islaggedvariables, pnames, enames, xnames, pid, eid, xid, data);
 
 % Get list and indices of estimated parameters.
-pnames_ = fieldnames(params);
-ipnames_ = zeros(size(pnames_));
-for i=1:length(ipnames_)
-    ipnames_(i) = find(strcmp(pnames_{i}, M_.param_names));
-end
+ipnames_ = get_estimated_parameters_indices(params, pnames, eqname, M_);
 
 % Create a routine for evaluating the residuals of the nonlinear model
-fun = ['r_' eqname];
-fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
-fprintf(fid, 'function r = %s(params, data, DynareModel, DynareOutput)\n', fun);
-fprintf(fid, '\n');
-fprintf(fid, '%% Evaluates the residuals for equation %s.\n', eqname);
-fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
-fprintf(fid, '\n');
-for i=1:length(ipnames_)
-    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames_(i), i);
-end
-fprintf(fid, '\n');
-fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
-fclose(fid);
+write_residuals_routine(lhs, rhs, eqname, ipnames_, M_);
 
 % Create a routine for evaluating the sum of squared residuals of the nonlinear model
-fun = ['ssr_' eqname];
-fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
-fprintf(fid, 'function [s, fake1, fake2, fake3, fake4] = %s(params, data, DynareModel, DynareOutput)\n', fun);
-fprintf(fid, '\n');
-fprintf(fid, '%% Evaluates the sum of square residuals for equation %s.\n', eqname);
-fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
-fprintf(fid, '\n');
-fprintf(fid, 'fake1 = 0;\n');
-fprintf(fid, 'fake2 = [];\n');
-fprintf(fid, 'fake3 = [];\n');
-fprintf(fid, 'fake4 = [];\n');
-fprintf(fid, '\n');
-for i=1:length(ipnames_)
-    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames_(i), i);
-end
-fprintf(fid, '\n');
-fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
-fprintf(fid, 's = r''*r;\n');
-fclose(fid);
+write_ssr_routine(lhs, rhs, eqname, ipnames_, M_);
 
 % Workaround for Octave bug https://savannah.gnu.org/bugs/?46282
 % Octave will randomly fail to read the ssr_* file generated in the +folder
diff --git a/matlab/+pac/+estimate/init.m b/matlab/+pac/+estimate/init.m
index e06e1c6ab867338b3ebccb7fb305d44a03f3047b..891770b5ab1afedce359306b8f33091fa6994e21 100644
--- a/matlab/+pac/+estimate/init.m
+++ b/matlab/+pac/+estimate/init.m
@@ -59,29 +59,8 @@ pacmodel = M_.pac.(pacmodl);
 % Get the parameters and variables in the PAC equation.
 [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_);
 
-% Get the list of estimated parameters
-pnames_ = fieldnames(params);
-
-% Check that the estimated parameters are used in the PAC equation.
-ParametersNotInPAC = setdiff(pnames_, pnames);
-if ~isempty(ParametersNotInPAC)
-    skipline()
-    if length(ParametersNotInPAC)>1
-        list = sprintf('  %s\n', ParametersNotInPAC{:});
-        remk = sprintf('  The following parameters:\n\n%s\n  do not appear in the PAC equation.', list);
-    else
-        remk = sprintf('  Parameter %s does not appear in the PAC equation.', ParametersNotInPAC{1});
-    end
-    disp(remk)
-    skipline()
-    error('The estimated parameters must be used in equation %s.', eqname)
-end
-
-% Get indices of estimated parameters.
-ipnames_ = zeros(size(pnames_));
-for i=1:length(ipnames_)
-    ipnames_(i) = find(strcmp(pnames_{i}, M_.param_names));
-end
+% Get list and indices of estimated parameters.
+ipnames_ = get_estimated_parameters_indices(params, pnames, eqname, M_);
 
 % If equation is estimated by recursive OLS, ensure that the error
 % correction parameter comes first, followed by the autoregressive
@@ -112,6 +91,7 @@ ipnames_(isnan(ipnames_)) = [];
 
 % Reorder params if needed.
 [~, permutation] = ismember(ipnames__, ipnames_);
+pnames_ = fieldnames(params);
 pnames_ = pnames_(permutation);
 params = orderfields(params, permutation);
 
diff --git a/matlab/+pac/+estimate/nls.m b/matlab/+pac/+estimate/nls.m
index 5c688ae8f0fd397ab5ab0b5b70de3aa9f3a57bf8..131aab9af8cac6e934e7a986e8711fc17c44a8a6 100644
--- a/matlab/+pac/+estimate/nls.m
+++ b/matlab/+pac/+estimate/nls.m
@@ -84,87 +84,17 @@ if M_.pac.(pacmodl).ec.istarget(2)
            'and the level of the endogenous variable.'], pacmodl);
 end
 
-% List of objects to be replaced
-objNames = [pnames; enames; xnames];
-objIndex = [pid; eid; xid];
-objTypes = [ones(length(pid), 1); 2*ones(length(eid)+length(xid), 1)];
-
-[~,I] = sort(cellfun(@length, objNames), 'descend');
-objNames = objNames(I);
-objIndex = objIndex(I);
-objTypes = objTypes(I);
+%
+% Rewrite and print the equation.
+%
 
-% Substitute parameters and variables.
-for i=1:length(objNames)
-    switch objTypes(i)
-      case 1
-        rhs = strrep(rhs, objNames{i}, sprintf('DynareModel.params(%u)', objIndex(i)));
-      case 2
-        k = find(strcmp(objNames{i}, data.name));
-        if isempty(k)
-            error('Variable %s is missing in the database.', objNames{i})
-        end
-        j = regexp(rhs, ['\<', objNames{i}, '\>']);
-        if islaggedvariables
-            jlag = regexp(rhs, ['\<', objNames{i}, '\(-1\)']);
-            if ~isempty(jlag)
-                rhs = regexprep(rhs, ['\<' objNames{i} '\(-1\)'], sprintf('data(1:end-1,%u)', k));
-            end
-            if ~isempty(setdiff(j, jlag))
-                rhs = regexprep(rhs, ['\<' objNames{i} '\>'], sprintf('data(2:end,%u)', k));
-            end
-        else
-            rhs = regexprep(rhs, ['\<' objNames{i} '\>'], sprintf('data(:,%u)', k));
-        end
-        if contains(lhs, objNames{i})
-            if islaggedvariables
-                lhs = strrep(lhs, objNames{i}, sprintf('data(2:end,%u)', k));
-            else
-                lhs = strrep(lhs, objNames{i}, sprintf('data(:,%u)', k));
-            end
-        end
-    end
-end
+[rhs, lhs] = rewrite_equation_with_tables(rhs, lhs, islaggedvariables, pnames, enames, xnames, pid, eid, xid, data);
 
 % Create a routine for evaluating the residuals of the nonlinear model
-fun = ['r_' eqname];
-fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
-fprintf(fid, 'function r = %s(params, data, DynareModel, DynareOutput)\n', fun);
-fprintf(fid, '\n');
-fprintf(fid, '%% Evaluates the residuals for equation %s.\n', eqname);
-fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
-fprintf(fid, '\n');
-for i=1:length(ipnames_)
-    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames_(i), i);
-end
-fprintf(fid, '\n');
-fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl);
-fprintf(fid, '\n');
-fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
-fclose(fid);
+write_residuals_routine(lhs, rhs, eqname, ipnames_, M_, pacmodl);
 
 % Create a routine for evaluating the sum of squared residuals of the nonlinear model
-fun = ['ssr_' eqname];
-fid = fopen(['+' M_.fname filesep() fun '.m'], 'w');
-fprintf(fid, 'function [s, fake1, fake2, fake3, fake4] = %s(params, data, DynareModel, DynareOutput)\n', fun);
-fprintf(fid, '\n');
-fprintf(fid, '%% Evaluates the sum of square residuals for equation %s.\n', eqname);
-fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
-fprintf(fid, '\n');
-fprintf(fid, 'fake1 = 0;\n');
-fprintf(fid, 'fake2 = [];\n');
-fprintf(fid, 'fake3 = [];\n');
-fprintf(fid, 'fake4 = [];\n');
-fprintf(fid, '\n');
-for i=1:length(ipnames_)
-    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames_(i), i);
-end
-fprintf(fid, '\n');
-fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl);
-fprintf(fid, '\n');
-fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
-fprintf(fid, 's = r''*r;\n');
-fclose(fid);
+write_ssr_routine(lhs, rhs, eqname, ipnames_, M_, pacmodl);
 
 % Workaround for Octave bug https://savannah.gnu.org/bugs/?46282
 % Octave will randomly fail to read the ssr_* file generated in the +folder
diff --git a/matlab/utilities/estimation/get_estimated_parameters_indices.m b/matlab/utilities/estimation/get_estimated_parameters_indices.m
new file mode 100644
index 0000000000000000000000000000000000000000..aea1110a7b15a6bc5d47b7231e0413f52cdaf419
--- /dev/null
+++ b/matlab/utilities/estimation/get_estimated_parameters_indices.m
@@ -0,0 +1,40 @@
+function ipnames = get_estimated_parameters_indices(params, pnames, eqname, DynareModel)
+
+% Copyright © 2021 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/>.
+
+list_of_parameters = fieldnames(params);
+
+% Check that the estimated parameters are used in the PAC equation.
+parameters_not_in_equation = setdiff(list_of_parameters, pnames);
+if ~isempty(parameters_not_in_equation)
+    skipline()
+    if length(parameters_not_in_equation)>1
+        list = sprintf('  %s\n', parameters_not_in_equation{:});
+        remk = sprintf('  The following parameters:\n\n%s\n  do not appear in equation %s.', list, eqname);
+    else
+        remk = sprintf('  Parameter %s does not appear in equation %s.', parameters_not_in_equation{1}, eqname);
+    end
+    disp(remk)
+    skipline()
+    error('The estimated parameters must be used in equation %s.', eqname)
+end
+
+ipnames = zeros(size(list_of_parameters));
+for i=1:length(ipnames)
+    ipnames(i) = find(strcmp(list_of_parameters{i}, DynareModel.param_names));
+end
\ No newline at end of file
diff --git a/matlab/utilities/estimation/rewrite_equation_with_tables.m b/matlab/utilities/estimation/rewrite_equation_with_tables.m
new file mode 100644
index 0000000000000000000000000000000000000000..21c08ed50424bf4339d731f127ed3128937eabb8
--- /dev/null
+++ b/matlab/utilities/estimation/rewrite_equation_with_tables.m
@@ -0,0 +1,70 @@
+function [expression, lhs] = rewrite_equation_with_tables(expression, lhs, islaggedvariables, pnames, enames, xnames, pid, eid, xid, data)
+
+% Takes and expression of parameters and variables (with max. lag equal to one)
+% and replace parameters by elements of a vector and variables by columns of a
+% data matrix.
+
+% Copyright © 2021 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/>.
+
+% List of objects to be replaced
+objNames = [pnames; enames; xnames];
+objIndex = [pid; eid; xid];
+objTypes = [ones(length(pid), 1); 2*ones(length(eid)+length(xid), 1)];
+
+[~,I] = sort(cellfun(@length, objNames), 'descend');
+objNames = objNames(I);
+objIndex = objIndex(I);
+objTypes = objTypes(I);
+
+% Substitute parameters and variables.  Note that in the transformed model (we consider rhs
+% instead of RHS) the maximum lag is 1.
+for i=1:length(objNames)
+    switch objTypes(i)
+      case 1
+        expression = strrep(expression, objNames{i}, sprintf('DynareModel.params(%u)', objIndex(i)));
+      case 2
+        k = find(strcmp(objNames{i}, data.name));
+        if isempty(k)
+            error('Variable %s is missing in the database.', objNames{i})
+        end
+        j = regexp(expression, ['\<', objNames{i}, '\>']);
+        if islaggedvariables
+            jlag = regexp(expression, ['\<', objNames{i}, '\(-1\)']);
+            if ~isempty(jlag)
+                expression = regexprep(expression, ['\<' objNames{i} '\(-1\)'], sprintf('data(1:end-1,%u)', k));
+            end
+            if ~isempty(setdiff(j, jlag))
+                expression = regexprep(expression, ['\<' objNames{i} '\>'], sprintf('data(2:end,%u)', k));
+            end
+        else
+            expression = regexprep(expression, ['\<' objNames{i} '\>'], sprintf('data(:,%u)', k));
+        end
+        if contains(lhs, objNames{i})
+            if islaggedvariables
+                lhs = strrep(lhs, objNames{i}, sprintf('data(2:end,%u)', k));
+            else
+                lhs = strrep(lhs, objNames{i}, sprintf('data(:,%u)', k));
+            end
+        end
+    end
+end
+
+% Allow elementwise operations
+expression = strrep(expression, '^', '.^');
+expression = strrep(expression, '/', './');
+expression = strrep(expression, '*', '.*');
\ No newline at end of file
diff --git a/matlab/utilities/estimation/write_residuals_routine.m b/matlab/utilities/estimation/write_residuals_routine.m
new file mode 100644
index 0000000000000000000000000000000000000000..f16a4924ed85ff02bfa3db838903bfc71bd253f9
--- /dev/null
+++ b/matlab/utilities/estimation/write_residuals_routine.m
@@ -0,0 +1,38 @@
+function write_residuals_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl)
+
+% Creates a routine for evaluating the residuals of the nonlinear equation.
+
+% Copyright © 2021 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/>.
+
+fun = sprintf('r_%s', eqname);
+fid = fopen(['+' DynareModel.fname filesep() fun '.m'], 'w');
+fprintf(fid, 'function r = %s(params, data, DynareModel, DynareOutput)\n', fun);
+fprintf(fid, '\n');
+fprintf(fid, '%% Evaluates the residuals for equation %s.\n', eqname);
+fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
+fprintf(fid, '\n');
+for i=1:length(ipnames)
+    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames(i), i);
+end
+fprintf(fid, '\n');
+if nargin>5
+    fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl);
+    fprintf(fid, '\n');
+end
+fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
+fclose(fid);
\ No newline at end of file
diff --git a/matlab/utilities/estimation/write_ssr_routine.m b/matlab/utilities/estimation/write_ssr_routine.m
new file mode 100644
index 0000000000000000000000000000000000000000..1e001f5d5878d739f86379afd309db1eb5a58f52
--- /dev/null
+++ b/matlab/utilities/estimation/write_ssr_routine.m
@@ -0,0 +1,44 @@
+function write_ssr_routine(lhs, rhs, eqname, ipnames, DynareModel, pacmodl)
+
+% Creates a routine for evaluating the sum of squared residuals of the nonlinear equation.
+
+% Copyright © 2021 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/>.
+
+fun = sprintf('ssr_%s', eqname);
+fid = fopen(['+' DynareModel.fname filesep() fun '.m'], 'w');
+fprintf(fid, 'function [s, fake1, fake2, fake3, fake4] = %s(params, data, DynareModel, DynareOutput)\n', fun);
+fprintf(fid, '\n');
+fprintf(fid, '%% Evaluates the sum of square residuals for equation %s.\n', eqname);
+fprintf(fid, '%% File created by Dynare (%s).\n', datetime);
+fprintf(fid, '\n');
+fprintf(fid, 'fake1 = 0;\n');
+fprintf(fid, 'fake2 = [];\n');
+fprintf(fid, 'fake3 = [];\n');
+fprintf(fid, 'fake4 = [];\n');
+fprintf(fid, '\n');
+for i=1:length(ipnames)
+    fprintf(fid, 'DynareModel.params(%u) = params(%u);\n', ipnames(i), i);
+end
+fprintf(fid, '\n');
+if nargin>5
+    fprintf(fid, 'DynareModel = pac.update.parameters(''%s'', DynareModel, DynareOutput, false);\n', pacmodl);
+    fprintf(fid, '\n');
+end
+fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
+fprintf(fid, 's = r''*r;\n');
+fclose(fid);