From d574705b4a944a17ea5a718ef194483d0537ae5b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 30 Nov 2022 14:42:54 +0100
Subject: [PATCH] Design and performance improvement to solve_algo={12,14}
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Use the new time-recursive block decomposition computed by the preprocessor
for:
- the simulation of backward models with “simul_backward”
- the perfect foresight simulation of purely backward/forward/static models

Also note that in this case, the preprocessor now defaults to “mfs=3” (i.e. it
minimizes the set of feedback variables and tries to renormalize equations).

This replaces the previous algorithm based on Dulmage-Mendelsohn (dmperm), plus
an ad hoc identification of some equations that can be evaluated (those with a
LHS equal to a variable, the log of a variable, or the diff-log of a variable).

By the way, the block_trust_region MEX has been modified so that it accepts a
boolean argument to decide whether it performs a Dulmage-Mendelsohn
decomposition (if not, then it performs a simple trust region on the whole
nonlinear system).

This provides a significant performance improvement (of almost an order of
magnitude for solve_algo=14 on a 700 equations model).
---
 doc/manual/source/the-model-file.rst          |  36 +--
 matlab/backward/simul_backward_model.m        |   6 +-
 .../simul_backward_nonlinear_model_.m         |  61 +++--
 matlab/dynare_solve.m                         | 142 ++++--------
 .../perfect_foresight_solver_core.m           |   3 +
 .../setup_time_recursive_block_simul.m        |  40 ++++
 .../sim1_purely_backward.m                    |  67 ++++--
 .../sim1_purely_forward.m                     |  65 ++++--
 .../sim1_purely_static.m                      |  63 +++--
 matlab/setup_solvers.m                        |  97 --------
 .../block_trust_region/mexFunction.f08        | 215 ++++++------------
 preprocessor                                  |   2 +-
 tests/Makefile.am                             |  49 +++-
 tests/nonlinearsolvers.m                      |   4 +-
 tests/solve_algo_12_14/backward_model.inc     |  28 +++
 tests/solve_algo_12_14/purely_backward_12.mod |  13 ++
 tests/solve_algo_12_14/purely_backward_14.mod |  13 ++
 .../purely_backward_common.inc                |  34 +++
 .../purely_backward_reference.mod             |   5 +
 tests/solve_algo_12_14/purely_forward_12.mod  |  13 ++
 tests/solve_algo_12_14/purely_forward_14.mod  |  13 ++
 .../purely_forward_common.inc                 |  60 +++++
 .../purely_forward_reference.mod              |   5 +
 tests/solve_algo_12_14/purely_static_12.mod   |  11 +
 tests/solve_algo_12_14/purely_static_14.mod   |  11 +
 .../solve_algo_12_14/purely_static_common.inc |  32 +++
 .../purely_static_reference.mod               |   5 +
 tests/solve_algo_12_14/simul_backward_12.mod  |  12 +
 tests/solve_algo_12_14/simul_backward_14.mod  |  12 +
 .../simul_backward_common.inc                 |  15 ++
 .../simul_backward_reference.mod              |   7 +
 31 files changed, 698 insertions(+), 441 deletions(-)
 create mode 100644 matlab/perfect-foresight-models/setup_time_recursive_block_simul.m
 delete mode 100644 matlab/setup_solvers.m
 create mode 100644 tests/solve_algo_12_14/backward_model.inc
 create mode 100644 tests/solve_algo_12_14/purely_backward_12.mod
 create mode 100644 tests/solve_algo_12_14/purely_backward_14.mod
 create mode 100644 tests/solve_algo_12_14/purely_backward_common.inc
 create mode 100644 tests/solve_algo_12_14/purely_backward_reference.mod
 create mode 100644 tests/solve_algo_12_14/purely_forward_12.mod
 create mode 100644 tests/solve_algo_12_14/purely_forward_14.mod
 create mode 100644 tests/solve_algo_12_14/purely_forward_common.inc
 create mode 100644 tests/solve_algo_12_14/purely_forward_reference.mod
 create mode 100644 tests/solve_algo_12_14/purely_static_12.mod
 create mode 100644 tests/solve_algo_12_14/purely_static_14.mod
 create mode 100644 tests/solve_algo_12_14/purely_static_common.inc
 create mode 100644 tests/solve_algo_12_14/purely_static_reference.mod
 create mode 100644 tests/solve_algo_12_14/simul_backward_12.mod
 create mode 100644 tests/solve_algo_12_14/simul_backward_14.mod
 create mode 100644 tests/solve_algo_12_14/simul_backward_common.inc
 create mode 100644 tests/solve_algo_12_14/simul_backward_reference.mod

diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 2ebc3b4dc..bebf95f34 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3011,27 +3011,27 @@ Finding the steady state with Dynare nonlinear solver
 
            ``12``
 
-                Specialized version of ``2`` for models where all the equations
-                have one endogenous variable on the left hand side and where
-                each equation determines a different endogenous variable. Only
-                expressions allowed on the left hand side are the natural
-                logarithm of an endogenous variable, the first difference of an
-                endogenous variable (with the ``diff`` operator), or the first
-                difference of the logarithm of an endogenous variable.
-                Univariate blocks are solved by evaluating the expression on the
-                right hand side.
+                Computes a block decomposition and then applies a Newton-type
+                solver on those smaller blocks rather than on the full
+                nonlinear system. This is similar to ``2``, but is typically
+                more efficient. The block decomposition is done at the
+                preprocessor level, which brings two benefits: it identifies
+                blocks that can be evaluated rather than solved; and evulations
+                of the residual and Jacobian of the model are more efficient
+                because only the relevant elements are recomputed at every
+                iteration.
 
            ``14``
 
-                Specialized version of ``4`` for models where all the equations
-                have one endogenous variable on the left hand side and where
-                each equation determines a different endogenous variable. Only
-                expressions allowed on the left hand side are the natural
-                logarithm of an endogenous variable, the first difference of an
-                endogenous variable (with the ``diff`` operator), or the first
-                difference of the logarithm of an endogenous variable..
-                Univariate blocks are solved by evaluating the expression on the
-                right hand side.
+                Computes a block decomposition and then applies a trust region
+                solver with autoscaling on those smaller blocks rather than on
+                the full nonlinear system. This is similar to ``4``, but is
+                typically more efficient. The block decomposition is done at
+                the preprocessor level, which brings two benefits: it
+                identifies blocks that can be evaluated rather than solved; and
+                evulations of the residual and Jacobian of the model are more
+                efficient because only the relevant elements are recomputed at
+                every iteration.
 
        |br| Default value is ``4``.
 
diff --git a/matlab/backward/simul_backward_model.m b/matlab/backward/simul_backward_model.m
index 103334624..eb4f13088 100644
--- a/matlab/backward/simul_backward_model.m
+++ b/matlab/backward/simul_backward_model.m
@@ -47,10 +47,6 @@ if ~M_.maximum_lag
     return
 end
 
-if ismember(options_.solve_algo, [12,14]) && ~M_.possible_to_use_solve_algo_12_14
-    error(M_.message_solve_algo_12_14)
-end
-
 if nargin<3
     Innovations = [];
 else
@@ -100,4 +96,4 @@ if options_.linear
     simulation = simul_backward_linear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
 else
     simulation = simul_backward_nonlinear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
-end
\ No newline at end of file
+end
diff --git a/matlab/backward/simul_backward_nonlinear_model_.m b/matlab/backward/simul_backward_nonlinear_model_.m
index 084ff0452..c369d5d16 100644
--- a/matlab/backward/simul_backward_nonlinear_model_.m
+++ b/matlab/backward/simul_backward_nonlinear_model_.m
@@ -40,36 +40,64 @@ function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, sampl
 
 debug = false;
 
-model_dynamic_s = str2func('dynamic_backward_model_for_simulation');
-
 if ~isempty(innovations)
     DynareOutput.exo_simul(initialconditions.nobs+(1:samplesize),:) = innovations;
 end
 
+if ismember(DynareOptions.solve_algo, [12,14])
+    [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(DynareModel);
+end
+
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+    % NB: do as few computations as possible inside this function, since it is
+    % called a very large number of times
+    y_dynamic(feedback_vars_idx) = z;
+    [~, ~, r, J] = feval(func, y_dynamic, x, DynareModel.params, DynareOutput.steady_state, ...
+                         sparse_rowval, sparse_colval, sparse_colptr, T);
+end
+
 % Simulations (call a Newton-like algorithm for each period).
 for it = initialconditions.nobs+(1:samplesize)
     if debug
         dprintf('Period t = %s.', num2str(it-initialconditions.nobs));
     end
     y_ = DynareOutput.endo_simul(:,it-1);
-    ylag = y_(iy1);                    % Set lagged variables.
     y = y_;                            % A good guess for the initial conditions is the previous values for the endogenous variables.
     try
         if ismember(DynareOptions.solve_algo, [12,14])
-            [DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = dynare_solve(model_dynamic_s, y, ...
-                                                                                       DynareOptions.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
-                                                                                       DynareOptions,  ...
-                                                                                       DynareModel.isloggedlhs, DynareModel.isauxdiffloggedrhs, DynareModel.endo_names, DynareModel.lhs, ...
-                                                                                       model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+            x = DynareOutput.exo_simul(it,:);
+            T = NaN(DynareModel.block_structure.dyn_tmp_nbr);
+            y_dynamic = [y_; y; NaN(DynareModel.endo_nbr, 1)];
+            for blk = 1:length(DynareModel.block_structure.block)
+                sparse_rowval = DynareModel.block_structure.block(blk).g1_sparse_rowval;
+                sparse_colval = DynareModel.block_structure.block(blk).g1_sparse_colval;
+                sparse_colptr = DynareModel.block_structure.block(blk).g1_sparse_colptr;
+                if DynareModel.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
+                    [z, errorflag, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+                                                                   DynareOptions.simul.maxit, DynareOptions.dynatol.f, ...
+                                                                   DynareOptions.dynatol.x, DynareOptions, ...
+                                                                   feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+                    if errorflag
+                        error('Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+                    end
+                    y_dynamic(feedback_vars_idxs{blk}) = z;
+                end
+                %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+                %% Also update the temporary terms vector.
+                [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, DynareModel.params, ...
+                                       DynareOutput.steady_state, sparse_rowval, sparse_colval, ...
+                                       sparse_colptr, T);
+            end
+            DynareOutput.endo_simul(:,it) = y_dynamic(DynareModel.endo_nbr+(1:DynareModel.endo_nbr));
         else
             [DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = ...
-                dynare_solve(model_dynamic_s, y, ...
+                dynare_solve(@dynamic_backward_model_for_simulation, y, ...
                              DynareOptions.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
                              DynareOptions, ...
-                             model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
-        end
-        if errorflag
-            error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+                             model_dynamic, y_(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+            if errorflag
+                error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+            end
         end
     catch
         DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
@@ -111,7 +139,7 @@ for it = initialconditions.nobs+(1:samplesize)
         %
         % Evaluate and check the residuals
         %
-        [r, J] = feval(model_dynamic_s, ytm, model_dynamic, ytm(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+        [r, J] = feval(@dynamic_backward_model_for_simulation, ytm, model_dynamic, ytm(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
         residuals_evaluating_to_nan = isnan(r);
         residuals_evaluating_to_inf = isinf(r);
         residuals_evaluating_to_complex = ~isreal(r);
@@ -142,6 +170,8 @@ end
 ysim = DynareOutput.endo_simul(1:DynareModel.orig_endo_nbr,:);
 xsim = DynareOutput.exo_simul;
 
+end
+
 function display_names_of_problematic_equations(DynareModel, eqtags, TruthTable)
 for i=1:DynareModel.orig_endo_nbr
     if TruthTable(i)
@@ -152,4 +182,5 @@ for i=DynareModel.orig_endo_nbr+1:DynareModel.endo_nbr
     if TruthTable(i)
         dprintf(' - Auxiliary equation for %s', DynareModel.endo_names{i})
     end
-end
\ No newline at end of file
+end
+end
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 17934f30e..04266e618 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -58,21 +58,10 @@ else
     in0 = nn;
 end
 
-% Get first element of varargin if solve_algo ∈ {12,14} and rename varargin.
-if ismember(options.solve_algo, [12, 14])
-    isloggedlhs = varargin{1};
-    isauxdiffloggedrhs = varargin{2};
-    endo_names = varargin{3};
-    lhs = varargin{4};
-    args = varargin(5:end);
-else
-    args = varargin;
-end
-
 % checking initial values
 % TODO We should have an option to deactivate the randomization.
 if jacobian_flag
-    [fvec, fjac] = feval(f, x, args{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
     wrong_initial_guess_flag = false;
     if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)))
         if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec))< tolf
@@ -88,7 +77,7 @@ if jacobian_flag
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = rand(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, args{:});
+            [fvec, fjac] = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
         end
         % If all previous attempts failed, try with real numbers.
@@ -96,7 +85,7 @@ if jacobian_flag
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = randn(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, args{:});
+            [fvec, fjac] = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
         end
         % Last tentative, ff all previous attempts failed, try with negative numbers.
@@ -104,12 +93,12 @@ if jacobian_flag
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = -rand(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, args{:});
+            [fvec, fjac] = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
         end
     end
 else
-    fvec = feval(f, x, args{:});
+    fvec = feval(f, x, varargin{:});
     fjac = zeros(nn, nn);
     if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec)) < tolf
         % return if initial value solves the problem except if a mixed complementarity problem is to be solved (complementarity conditions may not be satisfied)
@@ -125,7 +114,7 @@ else
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = rand(in0, 1)*10;
-            fvec = feval(f, x, args{:});
+            fvec = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec));
         end
         % If all previous attempts failed, try with real numbers.
@@ -133,7 +122,7 @@ else
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = randn(in0, 1)*10;
-            fvec = feval(f, x, args{:});
+            fvec = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec));
         end
         % Last tentative, ff all previous attempts failed, try with negative numbers.
@@ -141,7 +130,7 @@ else
         while wrong_initial_guess_flag && tentative_number<=in0*10
             tentative_number = tentative_number+1;
             x(idx) = -rand(in0, 1)*10;
-            fvec = feval(f, x, args{:});
+            fvec = feval(f, x, varargin{:});
             wrong_initial_guess_flag = ~all(isfinite(fvec));
         end
     end
@@ -199,7 +188,7 @@ if options.solve_algo == 0
     end
 
     if ~isoctave
-        [x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, args{:});
+        [x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, varargin{:});
     else
         % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
         if ischar(f)
@@ -207,7 +196,7 @@ if options.solve_algo == 0
         else
             f2 = f;
         end
-        [x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, args{:}), x, options4fsolve);
+        [x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, varargin{:}), x, options4fsolve);
     end
     if errorcode==1
         errorflag = false;
@@ -220,91 +209,42 @@ if options.solve_algo == 0
     else
         errorflag = true;
     end
-elseif options.solve_algo==1
-    [x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, args{:});
-    [fvec, fjac] = feval(f, x, args{:});
+elseif ismember(options.solve_algo, [1, 12])
+    %% NB: It is the responsibility of the caller to deal with the block decomposition if solve_algo=12
+    [x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, varargin{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
 elseif options.solve_algo==9
-    [x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, args{:});
-    [fvec, fjac] = feval(f, x, args{:});
-elseif ismember(options.solve_algo, [2, 12, 4])
-    if ismember(options.solve_algo, [2, 12])
+    [x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, varargin{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
+elseif ismember(options.solve_algo, [2, 4])
+    if options.solve_algo == 2
         solver = @solve1;
     else
         solver = @trust_region;
     end
-    specializedunivariateblocks = options.solve_algo == 12;
     if ~jacobian_flag
         fjac = zeros(nn,nn) ;
         dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
         for j = 1:nn
             xdh = x ;
             xdh(j) = xdh(j)+dh(j) ;
-            fjac(:,j) = (feval(f, xdh, args{:})-fvec)./dh(j) ;
+            fjac(:,j) = (feval(f, xdh, varargin{:})-fvec)./dh(j) ;
         end
     end
     [j1,j2,r,s] = dmperm(fjac);
-    JAC = abs(fjac(j1,j2))>0;
     if options.debug
-        disp(['DYNARE_SOLVE (solve_algo=2|4|12): number of blocks = ' num2str(length(r)-1)]);
+        disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r)-1)]);
     end
-    l = 0;
-    fre = false;
     for i=length(r)-1:-1:1
         blocklength = r(i+1)-r(i);
         if options.debug
-            dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u of size %u.', i, blocklength);
+            dprintf('DYNARE_SOLVE (solve_algo=2|4): solving block %u of size %u.', i, blocklength);
         end
         j = r(i):r(i+1)-1;
-        if specializedunivariateblocks
-            if options.debug
-                dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u by evaluating RHS.', i);
-            end
-            if isequal(blocklength, 1)
-                if i<length(r)-1
-                    if fre || any(JAC(r(i), s(i)+(1:l)))
-                        % Reevaluation of the residuals is required because the current RHS depends on
-                        % variables that potentially have been updated previously.
-                        z = feval(f, x, args{:});
-                        l = 0;
-                        fre = false;
-                    end
-                else
-                    % First iteration requires the evaluation of the residuals.
-                    z = feval(f, x, args{:});
-                end
-                l = l+1;
-                if isequal(lhs{j1(j)}, endo_names{j2(j)}) || isequal(lhs{j1(j)}, sprintf('log(%s)', endo_names{j2(j)}))
-                    if isloggedlhs(j1(j))
-                        x(j2(j)) = exp(log(x(j2(j)))-z(j1(j)));
-                    else
-                        x(j2(j)) = x(j2(j))-z(j1(j));
-                    end
-                else
-                    if options.debug
-                        dprintf('LHS variable is not determined by RHS expression (%u).', j1(j))
-                        dprintf('%s -> %s', lhs{j1(j)}, endo_names{j2(j)})
-                    end
-                    if ~isempty(regexp(lhs{j1(j)}, '\<AUX_DIFF_(\d*)\>', 'once'))
-                        if isauxdiffloggedrhs(j1(j))
-                            x(j2(j)) = exp(log(x(j2(j)))+z(j1(j)));
-                        else
-                            x(j2(j)) = x(j2(j))+z(j1(j));
-                        end
-                    else
-                        error('Algorithm solve_algo=%u cannot be used with this nonlinear problem.', options.solve_algo)
-                    end
-                end
-                continue
-            end
-        else
-            if options.debug
-                dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
-            end
-        end
         blockcolumns=s(i+1)-s(i);
         if blockcolumns ~= blocklength
             %non-square-block in DM; check whether initial value is solution
-            [fval_check, fjac] = feval(f, x, args{:});
+            [fval_check, fjac] = feval(f, x, varargin{:});
             if norm(fval_check(j1(j))) < tolf
                 errorflag = false;
                 errorcode = 0;
@@ -317,10 +257,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
                 options.gstep, ...
                 tolf, options.solve_tolx, maxit, ...
                 options.trust_region_initial_step_bound_factor, ...
-                options.debug, args{:});
-            fre = true;
+                options.debug, varargin{:});
         else
-            fprintf('\nDYNARE_SOLVE (solve_algo=2|4|12): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.\n')
+            fprintf('\nDYNARE_SOLVE (solve_algo=2|4): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.\n')
             %overdetermined block
             errorflag = true;
             errorcode = 0;
@@ -329,31 +268,31 @@ elseif ismember(options.solve_algo, [2, 12, 4])
             return
         end
     end
-    fvec = feval(f, x, args{:});
+    fvec = feval(f, x, varargin{:});
     if max(abs(fvec))>tolf
         disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
         [x, errorflag, errorcode] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
                                            options.gstep, tolf, options.solve_tolx, maxit, ...
                                            options.trust_region_initial_step_bound_factor, ...
-                                           options.debug, args{:});
+                                           options.debug, varargin{:});
     end
-    [fvec, fjac] = feval(f, x, args{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
 elseif options.solve_algo==3
     if jacobian_flag
-        [x, errorcode] = csolve(f, x, f, tolf, maxit, args{:});
+        [x, errorcode] = csolve(f, x, f, tolf, maxit, varargin{:});
     else
-        [x, errorcode] = csolve(f, x, [], tolf, maxit, args{:});
+        [x, errorcode] = csolve(f, x, [], tolf, maxit, varargin{:});
     end
     if errorcode==0
         errorflag = false;
     else
         errorflag = true;
     end
-    [fvec, fjac] = feval(f, x, args{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
 elseif options.solve_algo==10
     % LMMCP
     olmmcp = options.lmmcp;
-    [x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, args{:});
+    [x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, varargin{:});
     eq_to_check=find(isfinite(olmmcp.lb) | isfinite(olmmcp.ub));
     eq_to_ignore=eq_to_check(x(eq_to_check,:)<=olmmcp.lb(eq_to_check)+eps | x(eq_to_check,:)>=olmmcp.ub(eq_to_check)-eps);
     fvec(eq_to_ignore)=0;
@@ -373,7 +312,7 @@ elseif options.solve_algo == 11
     omcppath = options.mcppath;
     global mcp_data
     mcp_data.func = f;
-    mcp_data.args = args;
+    mcp_data.args = varargin;
     try
         [x, fval, jac, mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
     catch
@@ -384,18 +323,15 @@ elseif options.solve_algo == 11
     eq_to_ignore=eq_to_check(x(eq_to_check,:)<=omcppath.lb(eq_to_check)+eps | x(eq_to_check,:)>=omcppath.ub(eq_to_check)-eps);
     fvec(eq_to_ignore)=0;
 elseif ismember(options.solve_algo, [13, 14])
+    %% NB: It is the responsibility of the caller to deal with the block decomposition if solve_algo=14
     if ~jacobian_flag
-        error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')
+        error('DYNARE_SOLVE: option solve_algo=13 needs computed Jacobian')
     end
-    auxstruct = struct();
-    if options.solve_algo == 14
-        auxstruct.lhs = lhs;
-        auxstruct.endo_names = endo_names;
-        auxstruct.isloggedlhs = isloggedlhs;
-        auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
-    end
-    [x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, args{:});
-    [fvec, fjac] = feval(f, x, args{:});
+    [x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, ...
+                                                   options.trust_region_initial_step_bound_factor, ...
+                                                   options.solve_algo == 13, ... % Only block-decompose with Dulmage-Mendelsohn for 13, not for 14
+                                                   options.debug, varargin{:});
+    [fvec, fjac] = feval(f, x, varargin{:});
 else
     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
 end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 65120b5f8..7127000bf 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -50,6 +50,9 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
 end
 
 if options_.block
+    if M_.block_structure.time_recursive
+        error('Internal error: can''t perform stacked perfect foresight simulation with time-recursive block decomposition')
+    end
     if options_.bytecode
         try
             oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
diff --git a/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m b/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m
new file mode 100644
index 000000000..d9b671b41
--- /dev/null
+++ b/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m
@@ -0,0 +1,40 @@
+function [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_)
+%function [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_)
+%
+% For solve_algo={12,14}, precompute the function handles for per-block dynamic files
+% (it is surprisingly costly to recompute them within simulation loops).
+% Also precompute indices of feedback variables (also brings some performance gains).
+% By the way, do other sanity checks on block decomposition.
+
+% Copyright © 2022 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 ~isfield(M_, 'block_structure')
+    error('Can''t use solve_algo=12 nor solve_algo=14, because the block decomposition of the dynamic model failed')
+end
+if ~M_.block_structure.time_recursive
+    error('Can''t use solve_algo=12 nor solve_algo=14, because the model is not purely backward/static/forward or you gave the ''block'' option to the ''model'' block')
+end
+
+funcs = cell(length(M_.block_structure.block), 1);
+feedback_vars_idxs = cell(length(M_.block_structure.block), 1);
+for blk = 1:length(M_.block_structure.block)
+    funcs{blk} = str2func(sprintf('%s.sparse.block.dynamic_%d', M_.fname, blk));
+    feedback_vars_idxs{blk} = M_.endo_nbr+M_.block_structure.block(blk).variable((M_.block_structure.block(blk).endo_nbr-M_.block_structure.block(blk).mfs+1):end); % Indices of feedback variables in the dynamic y vector (of size 3n)
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index ca0dcd011..c22a2d7ac 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -20,38 +20,71 @@ function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables,
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 ny0 = nnz(M.lead_lag_incidence(2,:));    % Number of variables at current period
-iyb = M.lead_lag_incidence(1,:)>0;       % Logical vector (for lagged variables)
 
 if ny0 ~= M.endo_nbr
     error('All endogenous variables must appear at the current period!')
 end
 
-if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
-    error(M.message_solve_algo_12_14)
+if ismember(options.solve_algo, [12,14])
+    [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+    iyb = M.lead_lag_incidence(1,:)>0;       % Logical vector (for lagged variables)
+    dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+    dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
 end
 
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+    % NB: do as few computations as possible inside this function, since it is
+    % called a very large number of times
+    y_dynamic(feedback_vars_idx) = z;
+    [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+                         sparse_rowval, sparse_colval, sparse_colptr, T);
+end
 
 info.status = true;
 
 for it = M.maximum_lag + (1:options.periods)
     y = endogenousvariables(:,it-1);        % Values at previous period, also used as guess value for current period
-    ylag = y(iyb);
     if ismember(options.solve_algo, [12,14])
-        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
-                                                     options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
-                                                     options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
-                                                     dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+        x = exogenousvariables(it,:);
+        T = NaN(M.block_structure.dyn_tmp_nbr);
+        y_dynamic = [y; y; NaN(M.endo_nbr, 1)];
+        for blk = 1:length(M.block_structure.block)
+            sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+            sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+            sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+            if M.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
+                [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+                                                           options.simul.maxit, options.dynatol.f, ...
+                                                           options.dynatol.x, options, ...
+                                                           feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+                if check
+                    info.status = false;
+                    if options.debug
+                        dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+                    end
+                end
+                y_dynamic(feedback_vars_idxs{blk}) = z;
+            end
+            %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+            %% Also update the temporary terms vector.
+            [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+                                   steadystate, sparse_rowval, sparse_colval, ...
+                                   sparse_colptr, T);
+        end
+        endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
     else
+        ylag = y(iyb);
         [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
                                                      options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
                                                      options, dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+        if check
+            info.status = false;
+            dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)
+            break
+        end
+        endogenousvariables(:,it) = tmp;
     end
-    if check
-        info.status = false;
-        dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)
-        break
-    end
-    endogenousvariables(:,it) = tmp;
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index 715c544e8..a71715e2a 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -19,34 +19,71 @@ function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables,
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 ny0 = nnz(M.lead_lag_incidence(1,:));    % Number of variables at current period
-iyf = find(M.lead_lag_incidence(2,:)>0); % Indices of variables at next period
 
 if ny0 ~= M.endo_nbr
     error('All endogenous variables must appear at the current period!')
 end
 
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_forward_model_for_simulation');
+if ismember(options.solve_algo, [12,14])
+    [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+    iyf = find(M.lead_lag_incidence(2,:)>0); % Indices of variables at next period
+    dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+    dynamicmodel_s = str2func('dynamic_forward_model_for_simulation');
+end
+
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+    % NB: do as few computations as possible inside this function, since it is
+    % called a very large number of times
+    y_dynamic(feedback_vars_idx) = z;
+    [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+                         sparse_rowval, sparse_colval, sparse_colptr, T);
+end
 
 info.status = true;
 
 for it = options.periods:-1:1
     yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
     if ismember(options.solve_algo, [12,14])
-        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
-                                                     options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
-                                                     options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
-                                                     dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
+        x = exogenousvariables(it,:);
+        T = NaN(M.block_structure.dyn_tmp_nbr);
+        y_dynamic = [NaN(M.endo_nbr, 1); yf; yf];
+        for blk = 1:length(M.block_structure.block)
+            sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+            sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+            sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+            if M.block_structure.block(blk).Simulation_Type ~= 2 % Not an evaluate backward block
+                [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+                                                           options.simul.maxit, options.dynatol.f, ...
+                                                           options.dynatol.x, options, ...
+                                                           feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+                if check
+                    info.status = false;
+                    if options.debug
+                        dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+                    end
+                end
+                y_dynamic(feedback_vars_idxs{blk}) = z;
+            end
+            %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+            %% Also update the temporary terms vector.
+            [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+                                   steadystate, sparse_rowval, sparse_colval, ...
+                                   sparse_colptr, T);
+        end
+        endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
     else
         [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
                                                      options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
                                                      options, ...
                                                      dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
+        if check
+            info.status = false;
+            dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+            break
+        end
+        endogenousvariables(:,it) = tmp(1:M.endo_nbr);
     end
-    if check
-        info.status = false;
-        dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
-        break
-    end
-    endogenousvariables(:,it) = tmp(1:M.endo_nbr);
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m
index 8563d730d..9e0c6623b 100644
--- a/matlab/perfect-foresight-models/sim1_purely_static.m
+++ b/matlab/perfect-foresight-models/sim1_purely_static.m
@@ -23,12 +23,20 @@ if nnz(M.lead_lag_incidence(1,:)) ~= M.endo_nbr
     error('All endogenous variables must appear at the current period!')
 end
 
-if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
-    error(M.message_solve_algo_12_14)
+if ismember(options.solve_algo, [12,14])
+    [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+    dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+    dynamicmodel_s = str2func('dynamic_static_model_for_simulation');
 end
 
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_static_model_for_simulation');
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+    % NB: do as few computations as possible inside this function, since it is
+    % called a very large number of times
+    y_dynamic(feedback_vars_idx) = z;
+    [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+                         sparse_rowval, sparse_colval, sparse_colptr, T);
+end
 
 info.status = true;
 
@@ -36,21 +44,46 @@ y = endogenousvariables(:,1);
 
 for it = 1:options.periods
     if ismember(options.solve_algo, [12,14])
-        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
-                                                     options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
-                                                     options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
-                                                     dynamicmodel, exogenousvariables, M.params, steadystate, it);
+        x = exogenousvariables(it,:);
+        T = NaN(M.block_structure.dyn_tmp_nbr);
+        y_dynamic = [NaN(M.endo_nbr, 1); y; NaN(M.endo_nbr, 1)];
+        for blk = 1:length(M.block_structure.block)
+            sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+            sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+            sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+            if M.block_structure.block(blk).Simulation_Type ~= 1  % Not an evaluate forward block
+                [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+                                                           options.simul.maxit, options.dynatol.f, ...
+                                                           options.dynatol.x, options, ...
+                                                           feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+                if check
+                    info.status = false;
+                    if options.debug
+                        dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+                    end
+                end
+                y_dynamic(feedback_vars_idxs{blk}) = z;
+            end
+            %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+            %% Also update the temporary terms vector.
+            [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+                                   steadystate, sparse_rowval, sparse_colval, ...
+                                   sparse_colptr, T);
+        end
+        endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
     else
         [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
                                                      options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
                                                      options, dynamicmodel, exogenousvariables, M.params, steadystate, it);
-    end
-    if check
-        info.status = false;
-        if options.debug
-            dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+        if check
+            info.status = false;
+            if options.debug
+                dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+            end
         end
+        endogenousvariables(:,it) = tmp;
     end
-    endogenousvariables(:,it) = tmp;
     y = endogenousvariables(:,it);
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/setup_solvers.m b/matlab/setup_solvers.m
deleted file mode 100644
index 13aac7b53..000000000
--- a/matlab/setup_solvers.m
+++ /dev/null
@@ -1,97 +0,0 @@
-function Model = setup_solvers(Model)
-
-% Setup solve_algo={12,14} by identifying equations with a log on the left hand side.
-%
-% INPUTS
-% - Model     [struct]      Model description, aka M_.
-% - Options   [struct]      Dynare's options, aka options_.
-%
-% OUTPUTS
-% - Model     [struct]      Updated model description.
-
-% Copyright © 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 <https://www.gnu.org/licenses/>.
-
-cannot_use_solve_algo_12_14 = false;
-
-try
-    json = loadjson_(sprintf('%s/model/json/modfile.json', Model.fname));
-catch
-    cannot_use_solve_algo_12_14 = true;
-    message = 'Algorithms solve_algo={12,14} require json output of the model (use json=compute option)';
-end
-
-if ~cannot_use_solve_algo_12_14
-
-    lhs = cell(length(json.model),1);
-    isauxdiffloggedrhs = false(length(json.model), 1);
-
-    for i = 1:length(json.model)
-        if length(json.model)>1
-            lhs{i} = json.model{i}.lhs;
-        else
-            lhs{i} = json.model.lhs;
-        end
-        if isempty(regexp(lhs{i}, '^\w+$|^log\(\w+\)$'))
-            cannot_use_solve_algo_12_14 = true;
-            message = sprintf('With solve_algo={12,14}, each equation must have on the left hand side a single variable or logged variable (equation %d does not satisfy this condition).', i);
-            break
-        end
-        if length(json.model)>1
-            rhs = json.model{i}.rhs;
-        else
-            rhs = json.model.lhs;
-        end
-        if i>Model.orig_endo_nbr &&  ~isempty(regexp(lhs{i}, '\<AUX_DIFF_(\d*)\>', 'once')) && ismember(lhs{i}, lhs(1:i-1)) && ...
-                ~isempty(regexp(rhs, 'log\(\w*\)-log\(\w*\(-1\)\)', 'once'))
-            isauxdiffloggedrhs(i) = true;
-        end
-    end
-
-end
-
-if ~cannot_use_solve_algo_12_14
-
-    islog = @(x) ~isempty(regexp(x, 'log\(\w*\)', 'once'));
-    
-    lhs0 = lhs;
-    for i=1:length(json.model)
-        if islog(lhs{i})
-            lhs0{i} = strrep(strrep(lhs{i}, 'log(', ''), ')', '');
-        end
-    end
-    
-    if ~isequal(length(unique(lhs0(1:Model.orig_endo_nbr))), length(lhs0(1:Model.orig_endo_nbr)))
-        cannot_use_solve_algo_12_14 = true;
-        message = sprintf('With solve_algo={12,14}, each equation must determine a different endogenous variable.')
-    end
-
-end
-
-if cannot_use_solve_algo_12_14
-    Model.isloggedlhs = {};
-    Model.lhs = {};
-    Model.isauxdiffloggedrhs = [];
-    Model.possible_to_use_solve_algo_12_14 = false;
-    Model.message_solve_algo_12_14 = message;
-else
-    Model.isloggedlhs = cellfun(islog, lhs);
-    Model.lhs = lhs;
-    Model.isauxdiffloggedrhs = isauxdiffloggedrhs;
-    Model.possible_to_use_solve_algo_12_14 = true;
-    Model.message_solve_algo_12_14 = '';
-end
\ No newline at end of file
diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 59da65547..11f9fccf2 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/mex/sources/block_trust_region/mexFunction.f08
@@ -35,16 +35,11 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   integer :: maxiter
   real(real64), dimension(:), allocatable :: fvec
   real(real64), dimension(:,:), allocatable :: fjac
-  logical :: debug, specializedunivariateblocks
+  logical :: debug, block_decompose
   character(len=80) :: debug_msg
-  logical(mxLogical), dimension(:), pointer :: isloggedlhs => null(), &
-       isauxdiffloggedrhs => null()
-  type(c_ptr) :: endo_names, lhs
-  logical :: fre ! True if the last block has been solved (i.e. not evaluated), so that residuals must be updated
-  integer, dimension(:), allocatable :: evaled_cols ! If fre=.false., lists the columns that have been evaluated so far without updating the residuals
 
-  if (nrhs < 4 .or. nlhs /= 3) then
-     call mexErrMsgTxt("Must have at least 7 inputs and exactly 3 outputs")
+  if (nrhs < 8 .or. nlhs /= 3) then
+     call mexErrMsgTxt("Must have at least 8 inputs and exactly 3 outputs")
   end if
 
   if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
@@ -73,178 +68,96 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   end if
 
   if (.not. (mxIsLogicalScalar(prhs(7)))) then
-     call mexErrMsgTxt("Seventh argument (debug) should be a logical scalar")
+     call mexErrMsgTxt("Seventh argument (block_decompose) should be a logical scalar")
   end if
 
-  if (.not. (mxIsStruct(prhs(8)) .and. &
-       (mxGetNumberOfFields(prhs(8)) == 0 .or. mxGetNumberOfFields(prhs(8)) == 4))) then
-     call mexErrMsgTxt("Eighth argument should be a struct with either 0 or 4 fields")
+  if (.not. (mxIsLogicalScalar(prhs(8)))) then
+     call mexErrMsgTxt("Eigth argument (debug) should be a logical scalar")
   end if
 
-  specializedunivariateblocks = (mxGetNumberOfFields(prhs(8)) == 4)
-
   func => prhs(1)
   tolf = mxGetScalar(prhs(3))
   tolx = mxGetScalar(prhs(4))
   maxiter = int(mxGetScalar(prhs(5)))
   factor = mxGetScalar(prhs(6))
-  debug = mxGetScalar(prhs(7)) == 1._c_double
-  extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 9 and subsequent ones
+  block_decompose = mxGetScalar(prhs(7)) == 1._c_double
+  debug = mxGetScalar(prhs(8)) == 1._c_double
+  extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 8 and subsequent ones
   associate (x_mat => mxGetPr(prhs(2)))
     allocate(x(size(x_mat)))
     x = x_mat
   end associate
-  if (specializedunivariateblocks) then
-     block
-       type(c_ptr) :: tmp
-       tmp = mxGetField(prhs(8), 1_mwIndex, "isloggedlhs")
-       if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
-          call mexErrMsgTxt("Seventh argument must have a 'isloggedlhs' field of type logical, of same size as second argument")
-       end if
-       isloggedlhs => mxGetLogicals(tmp)
-
-       tmp = mxGetField(prhs(8), 1_mwIndex, "isauxdiffloggedrhs")
-       if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
-          call mexErrMsgTxt("Seventh argument must have a 'isauxdiffloggedrhs' field of type &
-               &logical, of same size as second argument")
-       end if
-       isauxdiffloggedrhs => mxGetLogicals(tmp)
-
-       lhs = mxGetField(prhs(8), 1_mwIndex, "lhs")
-       if (.not. (c_associated(lhs) .and. mxIsCell(lhs) .and. mxGetNumberOfElements(lhs) == size(x))) then
-          call mexErrMsgTxt("Seventh argument must have a 'lhs' field of type cell, of same size as second argument")
-       end if
-
-       endo_names = mxGetField(prhs(8), 1_mwIndex, "endo_names")
-       if (.not. (c_associated(endo_names) .and. mxIsCell(endo_names) .and. mxGetNumberOfElements(endo_names) == size(x))) then
-          call mexErrMsgTxt("Seventh argument must have a 'endo_names' field of type cell, of same size as second argument")
-       end if
-     end block
-
-     allocate(evaled_cols(0))
-     fre = .false.
-  end if
 
   allocate(fvec(size(x)))
   allocate(fjac(size(x), size(x)))
 
-  ! Compute block decomposition
-  nullify(x_indices, f_indices, x_all)
-  call matlab_fcn(x, fvec, fjac)
-  call dm_blocks(fjac, blocks)
+  if (block_decompose) then
+     ! Compute block decomposition
+     nullify(x_indices, f_indices, x_all)
+     call matlab_fcn(x, fvec, fjac)
+     call dm_blocks(fjac, blocks)
 
-  if (debug) then
-     write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): number of blocks = ', i0)") size(blocks)
-     call mexPrintf_trim_newline(debug_msg)
-  end if
-
-  ! Solve the system, starting from bottom-rightmost block
-  do i = size(blocks),1,-1
      if (debug) then
-        write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' of size ', i0)") &
-             i, size(blocks(i)%col_indices)
+        write (debug_msg, "('DYNARE_SOLVE (solve_algo=13): number of blocks = ', i0)") size(blocks)
         call mexPrintf_trim_newline(debug_msg)
      end if
 
-     if (specializedunivariateblocks .and. size(blocks(i)%col_indices) == 1) then
+     ! Solve the system, starting from bottom-rightmost block
+     do i = size(blocks),1,-1
         if (debug) then
-           write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' by evaluating RHS')") i
+           write (debug_msg, "('DYNARE_SOLVE (solve_algo=13): solving block ', &
+                &i0, ' of size ', i0)") &
+                i, size(blocks(i)%col_indices)
            call mexPrintf_trim_newline(debug_msg)
         end if
-        associate (eq => blocks(i)%row_indices(1), var => blocks(i)%col_indices(1))
-          if (fre .or. any(abs(fjac(eq, evaled_cols)) > 0._real64)) then
-             ! Reevaluation of the residuals is required because the current RHS depends on
-             ! variables that potentially have been updated previously.
-             nullify(x_indices, f_indices, x_all)
-             call matlab_fcn(x, fvec)
-             deallocate(evaled_cols) ! This shouldn’t be necessary, but it crashes otherwise with gfortran 8
-             allocate(evaled_cols(0))
-             fre = .false.
-          end if
-          evaled_cols = [ evaled_cols, var]
-          block
-            ! An associate() construct for lhs_eq and endo_name_var makes the
-            ! code crash (with double free) using gfortran 8. Hence use a block
-            character(kind=c_char, len=:), allocatable :: lhs_eq, endo_name_var
-            lhs_eq = mxArrayToString(mxGetCell(lhs, int(eq, mwIndex)))
-            endo_name_var = mxArrayToString(mxGetCell(endo_names, int(var, mwIndex)))
-            if (lhs_eq == endo_name_var .or. lhs_eq == "log(" // endo_name_var // ")") then
-               if (isloggedlhs(eq)) then
-                  x(var) = exp(log(x(var)) - fvec(eq))
-               else
-                  x(var) = x(var) - fvec(eq)
-               end if
-            else
-               if (debug) then
-                  write (debug_msg, "('LHS variable is not determined by RHS expression (', i0, ')')") eq
-                  call mexPrintf_trim_newline(debug_msg)
-                  write (debug_msg, "(a, ' -> ', a)") lhs_eq, endo_name_var
-                  call mexPrintf_trim_newline(debug_msg)
-               end if
-               if (lhs_eq(1:9) == "AUX_DIFF_" .or. lhs_eq(1:13) == "log(AUX_DIFF_") then
-                  if (isauxdiffloggedrhs(eq)) then
-                     x(var) = exp(log(x(var)) + fvec(eq))
-                  else
-                     x(var) = x(var) + fvec(eq)
-                  end if
-               else
-                  call mexErrMsgTxt("Algorithm solve_algo=14 cannot be used with this nonlinear problem")
-               end if
-            end if
-          end block
-        end associate
-        cycle
+
+        if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
+           ! Non-square block in DM decomposition
+           ! Before erroring out, check whether we are not already at the solution for this block
+           ! See also #1851
+           if (norm2(fvec(blocks(i)%row_indices)) < tolf) then
+              cycle
+           else
+              call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13): the Dulmage-Mendelsohn &
+                   &decomposition returned a non-square block. This means that the &
+                   &Jacobian is singular. You may want to try another value for solve_algo.")
+           end if
+        end if
+
+        block
+          real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
+          x_indices => blocks(i)%col_indices
+          f_indices => blocks(i)%row_indices
+          x_all => x
+          x_block = x(x_indices)
+          call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
+          x(x_indices) = x_block
+        end block
+     end do
+
+     ! Verify that we have a solution
+     ! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
+     ! this check would almost always fail (because the 2-norm of the full fvec is
+     ! larger than the 2-norm of its sub-vectors)
+     ! If the check fails, this normally means that the block decomposition was
+     ! incorrect (because some element of the Jacobian was numerically zero at the
+     ! guess value while not being symbolically zero)
+     nullify(x_indices, f_indices, x_all)
+     call matlab_fcn(x, fvec)
+     if (maxval(abs(fvec)) > tolf) then
+        if (debug) &
+             call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13): residuals&
+             & still too large, solving for the whole model")
+        call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
      else
-        if (debug) then
-           write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' with trust region routine')") i
-           call mexPrintf_trim_newline(debug_msg)
+        if (size(blocks) > 1) then
+           ! Note that the value of info may be different across blocks
+           info = 1
         end if
      end if
-
-     if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
-        ! Non-square block in DM decomposition
-        ! Before erroring out, check whether we are not already at the solution for this block
-        ! See also #1851
-        if (norm2(fvec(blocks(i)%row_indices)) < tolf) then
-           cycle
-        else
-           call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13|14): the Dulmage-Mendelsohn &
-                &decomposition returned a non-square block. This means that the &
-                &Jacobian is singular. You may want to try another value for solve_algo.")
-        end if
-     end if
-
-     block
-       real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
-       x_indices => blocks(i)%col_indices
-       f_indices => blocks(i)%row_indices
-       x_all => x
-       x_block = x(x_indices)
-       call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
-       x(x_indices) = x_block
-     end block
-
-     fre = .true.
-  end do
-
-  ! Verify that we have a solution
-  ! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
-  ! this check would almost always fail (because the 2-norm of the full fvec is
-  ! larger than the 2-norm of its sub-vectors)
-  ! If the check fails, this normally means that the block decomposition was
-  ! incorrect (because some element of the Jacobian was numerically zero at the
-  ! guess value while not being symbolically zero)
-  nullify(x_indices, f_indices, x_all)
-  call matlab_fcn(x, fvec)
-  if (maxval(abs(fvec)) > tolf) then
-     if (debug) &
-          call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13|14): residuals still too large, solving for the whole model")
+  else ! No block decomposition
+     nullify(x_indices, f_indices, x_all)
      call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
-  else
-     if (size(blocks) > 1) then
-        ! Note that the value of info may be different across blocks
-        info = 1
-     end if
   end if
 
   plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
diff --git a/preprocessor b/preprocessor
index f725c534e..c48248fc0 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit f725c534ef1344813421c2404de2cb4f3774d113
+Subproject commit c48248fc0d65ed70fb13dd508ed0254f113cbcb6
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b8b043553..e0f8aead5 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -498,7 +498,19 @@ MODFILES = \
 	log_transform/example1.mod \
 	log_transform/ramst.mod \
 	log_transform/fs2000_nonstationary.mod \
-	log_transform/nk_ramsey.mod
+	log_transform/nk_ramsey.mod \
+	solve_algo_12_14/simul_backward_reference.mod \
+	solve_algo_12_14/simul_backward_12.mod \
+	solve_algo_12_14/simul_backward_14.mod \
+	solve_algo_12_14/purely_backward_reference.mod \
+	solve_algo_12_14/purely_backward_12.mod \
+	solve_algo_12_14/purely_backward_14.mod \
+	solve_algo_12_14/purely_static_reference.mod \
+	solve_algo_12_14/purely_static_12.mod \
+	solve_algo_12_14/purely_static_14.mod \
+	solve_algo_12_14/purely_forward_reference.mod \
+	solve_algo_12_14/purely_forward_12.mod \
+	solve_algo_12_14/purely_forward_14.mod
 
 ECB_MODFILES = \
 	var-expectations/1/example1.mod \
@@ -994,6 +1006,27 @@ model-inversion/nk-2/z_check_inversion.o.trs: model-inversion/nk-2/invert.o.trs
 ep/rbcii_MCP.m.trs: ep/rbcii.m.trs
 ep/rbcii_MCP.o.trs: ep/rbcii.o.trs
 
+solve_algo_12_14/simul_backward_12.m.trs: solve_algo_12_14/simul_backward_reference.m.trs
+solve_algo_12_14/simul_backward_12.o.trs: solve_algo_12_14/simul_backward_reference.o.trs
+solve_algo_12_14/simul_backward_14.m.trs: solve_algo_12_14/simul_backward_reference.m.trs
+solve_algo_12_14/simul_backward_14.o.trs: solve_algo_12_14/simul_backward_reference.o.trs
+
+solve_algo_12_14/purely_backward_12.m.trs: solve_algo_12_14/purely_backward_reference.m.trs
+solve_algo_12_14/purely_backward_12.o.trs: solve_algo_12_14/purely_backward_reference.o.trs
+solve_algo_12_14/purely_backward_14.m.trs: solve_algo_12_14/purely_backward_reference.m.trs
+solve_algo_12_14/purely_backward_14.o.trs: solve_algo_12_14/purely_backward_reference.o.trs
+
+solve_algo_12_14/purely_static_12.m.trs: solve_algo_12_14/purely_static_reference.m.trs
+solve_algo_12_14/purely_static_12.o.trs: solve_algo_12_14/purely_static_reference.o.trs
+solve_algo_12_14/purely_static_14.m.trs: solve_algo_12_14/purely_static_reference.m.trs
+solve_algo_12_14/purely_static_14.o.trs: solve_algo_12_14/purely_static_reference.o.trs
+
+solve_algo_12_14/purely_forward_12.m.trs: solve_algo_12_14/purely_forward_reference.m.trs
+solve_algo_12_14/purely_forward_12.o.trs: solve_algo_12_14/purely_forward_reference.o.trs
+solve_algo_12_14/purely_forward_14.m.trs: solve_algo_12_14/purely_forward_reference.m.trs
+solve_algo_12_14/purely_forward_14.o.trs: solve_algo_12_14/purely_forward_reference.o.trs
+
+
 observation_trends_and_prefiltering/MCMC: m/observation_trends_and_prefiltering/MCMC o/observation_trends_and_prefiltering/MCMC
 m/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.m.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
 o/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.o.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
@@ -1210,6 +1243,11 @@ model-inversion: m/model-inversion o/model-inversion
 m/model-inversion: $(patsubst %.mod, %.m.trs, $(filter model-inversion/%.mod, $(MODFILES)))
 o/model-inversion: $(patsubst %.mod, %.o.trs, $(filter model-inversion/%.mod, $(MODFILES)))
 
+solve_algo_12_14: m/solve_algo_12_14 o/solve_algo_12_14
+m/solve_algo_12_14: $(patsubst %.mod, %.m.trs, $(filter solve_algo_12_14/%.mod, $(MODFILES)))
+o/solve_algo_12_14: $(patsubst %.mod, %.o.trs, $(filter solve_algo_12_14/%.mod, $(MODFILES)))
+
+
 # ECB files
 M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
 
@@ -1448,7 +1486,12 @@ EXTRA_DIST = \
 	solver-test-functions/variablydimensioned.m \
 	solver-test-functions/watson.m \
 	solver-test-functions/wood.m \
-	ramst_normcdf_and_friends.inc
+	ramst_normcdf_and_friends.inc \
+	solve_algo_12_14/backward_model.inc \
+	solve_algo_12_14/simul_backward_common.inc \
+	solve_algo_12_14/purely_backward_common.inc \
+	solve_algo_12_14/purely_static_common.inc \
+	solve_algo_12_14/purely_forward_common.inc
 
 if ENABLE_MATLAB
 check-local: check-matlab
@@ -1639,6 +1682,8 @@ clean-local:
 
 	rm -rf tests/pac/var-12/toto
 
+	rm -f solve_algo_12_14/simul_backward_ref.mat
+
 	find . -name "*.tex" -type f -delete
 	find . -name "*.aux" -type f -delete
 	find . -name "*.log" -type f -delete
diff --git a/tests/nonlinearsolvers.m b/tests/nonlinearsolvers.m
index aab3e2a7a..b8e07f3c4 100644
--- a/tests/nonlinearsolvers.m
+++ b/tests/nonlinearsolvers.m
@@ -31,8 +31,6 @@ tolx = 1e-6;
 maxit = 50;
 factor = 10;
 
-auxstruct = struct();
-
 % List of function handles
 objfun = { @rosenbrock,
            @powell1,
@@ -73,7 +71,7 @@ for i=1:length(objfun)
         x = objfun{i}();
     end
     try
-        [x, errorflag, exitflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
+        [x, errorflag, exitflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, true, false);
         if isequal(func2str(objfun{i}), 'powell2')
             if ~errorflag
                 testFailed = testFailed+1;
diff --git a/tests/solve_algo_12_14/backward_model.inc b/tests/solve_algo_12_14/backward_model.inc
new file mode 100644
index 000000000..43c097c68
--- /dev/null
+++ b/tests/solve_algo_12_14/backward_model.inc
@@ -0,0 +1,28 @@
+/* A simple purely backward model, used for testing both “simul_backward” and
+   perfect foresight simulations, consisting of:
+  - a VAR(3) with one variable in level, one in log and one log-differentiated;
+  - two other variables that depend on the contemporaneous values of the VAR(3),
+    and which together for a simultaneous block.
+  Those five equations are repeated an arbitrary number of times (n). */
+
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, y@{i}, z@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_y@{i}, e_z@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model(use_dll);
+@#for i in 1:n
+x@{i} = 0.1*x@{i}(-1) + 0.2*log(y@{i}(-1)) + 0.3*diff(log(z@{i}(-1))) + e_x@{i};
+[name = 'y']
+log(y@{i}) = 0.4*x@{i}(-1) + 0.5*log(y@{i}(-1)) + 0.6*diff(log(z@{i}(-1))) + e_y@{i};
+[name = 'z']
+/* NB: we choose 0.5 as steady state for z, since initial value is drawn from
+   [0,1] (in the “simul_backward” case) */
+diff(log(z@{i})) = -0.8*(log(z@{i}(-1)) - log(0.5)) + 0.1*x@{i}(-1) + 0.2*log(y@{i}(-1)) - 0.3*diff(log(z@{i}(-1))) + e_z@{i};
+
+t@{i} = x@{i} + 2*log(y@{i}) + 3*diff(log(z@{i})) + s@{i} + e_t@{i};
+s@{i} = 4*x@{i}(-1) + 3*t@{i} + e_s@{i};
+@#endfor
+end;
diff --git a/tests/solve_algo_12_14/purely_backward_12.mod b/tests/solve_algo_12_14/purely_backward_12.mod
new file mode 100644
index 000000000..ffcfbc9ed
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_12.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely backward model with “solve_algo=12”
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_backward_reference/Output/purely_backward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 5e-5
+    error('Incorrect results for perfect foresight with solve_algo=12')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_backward_14.mod b/tests/solve_algo_12_14/purely_backward_14.mod
new file mode 100644
index 000000000..238c18855
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_14.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely backward model with “solve_algo=14”
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_backward_reference/Output/purely_backward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 5e-5
+    error('Incorrect results for perfect foresight with solve_algo=14')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_backward_common.inc b/tests/solve_algo_12_14/purely_backward_common.inc
new file mode 100644
index 000000000..089bc87df
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_common.inc
@@ -0,0 +1,34 @@
+@#include "backward_model.inc"
+
+initval;
+@#for i in 1:n
+y@{i} = 1;
+z@{i} = 0.5;
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 1;
+values @{0.1/n};
+
+var e_y@{i};
+periods 2;
+values @{0.05+0.1/n};
+
+var e_z@{i};
+periods 3;
+values @{0.1+0.1/n};
+
+var e_t@{i};
+periods 4;
+values @{0.15+0.1/n};
+
+var e_s@{i};
+periods 5;
+values @{0.2+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_backward_reference.mod b/tests/solve_algo_12_14/purely_backward_reference.mod
new file mode 100644
index 000000000..397cc4ace
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely backward model with default “solve_algo” value
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/purely_forward_12.mod b/tests/solve_algo_12_14/purely_forward_12.mod
new file mode 100644
index 000000000..46334bf89
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_12.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely forward model with “solve_algo=12”
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_forward_reference/Output/purely_forward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 2e-4
+    error('Incorrect results for perfect foresight with solve_algo=12')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_forward_14.mod b/tests/solve_algo_12_14/purely_forward_14.mod
new file mode 100644
index 000000000..c3c9eda65
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_14.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely forward model with “solve_algo=14”
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_forward_reference/Output/purely_forward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 2e-4
+    error('Incorrect results for perfect foresight with solve_algo=14')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_forward_common.inc b/tests/solve_algo_12_14/purely_forward_common.inc
new file mode 100644
index 000000000..16df0d2dc
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_common.inc
@@ -0,0 +1,60 @@
+/* A simple purely forward model, consisting of:
+  - a forward VAR(3) with one variable in level, one in log and one log-forward-differentiated;
+  - two other variables that depend on the contemporaneous values of the VAR(3),
+    and which together for a simultaneous block.
+  Those five equations are repeated an arbitrary number of times (n). */
+
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, y@{i}, z@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_y@{i}, e_z@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model;
+@#for i in 1:n
+x@{i} = 0.1*x@{i}(+1) + 0.2*log(y@{i}(+1)) + 0.3*(log(z@{i})-log(z@{i}(+1))) + e_x@{i};
+[name = 'y']
+log(y@{i}) = 0.4*x@{i}(+1) + 0.5*log(y@{i}(+1)) + 0.6*(log(z@{i})-log(z@{i}(+1))) + e_y@{i};
+[name = 'z']
+/* NB: we choose 0.5 as steady state for z, since initial value is drawn from
+   [0,1] (in the “simul_backward” case) */
+log(z@{i})-log(z@{i}(+1)) = -0.8*(log(z@{i}(+1)) - log(0.5)) + 0.1*x@{i}(+1) + 0.2*log(y@{i}(+1)) - 0.3*(log(z@{i}) - log(z@{i}(+1))) + e_z@{i};
+
+t@{i} = x@{i} + 2*log(y@{i}) + 3*(log(z@{i})-log(z@{i}(+1))) + s@{i} + e_t@{i};
+s@{i} = 4*x@{i}(+1) + 3*t@{i} + e_s@{i};
+@#endfor
+end;
+
+initval;
+@#for i in 1:n
+y@{i} = 1;
+z@{i} = 0.5;
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 99;
+values @{0.1/n};
+
+var e_y@{i};
+periods 98;
+values @{0.05+0.1/n};
+
+var e_z@{i};
+periods 97;
+values @{0.1+0.1/n};
+
+var e_t@{i};
+periods 96;
+values @{0.15+0.1/n};
+
+var e_s@{i};
+periods 95;
+values @{0.2+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_forward_reference.mod b/tests/solve_algo_12_14/purely_forward_reference.mod
new file mode 100644
index 000000000..70ec7cfb4
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely forward model with default “solve_algo” value
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/purely_static_12.mod b/tests/solve_algo_12_14/purely_static_12.mod
new file mode 100644
index 000000000..9b4da090f
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_12.mod
@@ -0,0 +1,11 @@
+// Check the correctedness of perfect foresight simulation of a purely static model with “solve_algo=12”
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_static_reference/Output/purely_static_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 1e-16
+    error('Incorrect results for perfect foresight with solve_algo=12')
+end
diff --git a/tests/solve_algo_12_14/purely_static_14.mod b/tests/solve_algo_12_14/purely_static_14.mod
new file mode 100644
index 000000000..3044f4ebf
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_14.mod
@@ -0,0 +1,11 @@
+// Check the correctedness of perfect foresight simulation of a purely static model with “solve_algo=14”
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_static_reference/Output/purely_static_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 1e-16
+    error('Incorrect results for perfect foresight with solve_algo=14')
+end
diff --git a/tests/solve_algo_12_14/purely_static_common.inc b/tests/solve_algo_12_14/purely_static_common.inc
new file mode 100644
index 000000000..14189814f
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_common.inc
@@ -0,0 +1,32 @@
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model;
+@#for i in 1:n
+x@{i} = e_x@{i};
+t@{i} = x@{i} + s@{i} + e_t@{i};
+s@{i} = 3*t@{i} + e_s@{i};
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 1;
+values @{0.1/n};
+
+var e_t@{i};
+periods 2;
+values @{0.05+0.1/n};
+
+var e_s@{i};
+periods 3;
+values @{0.1+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_static_reference.mod b/tests/solve_algo_12_14/purely_static_reference.mod
new file mode 100644
index 000000000..0a619b41c
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely static model with default “solve_algo” value
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/simul_backward_12.mod b/tests/solve_algo_12_14/simul_backward_12.mod
new file mode 100644
index 000000000..108d8f603
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_12.mod
@@ -0,0 +1,12 @@
+// Check the correctedness of “simul_backward” with “solve_algo=12”
+
+@#include "simul_backward_common.inc"
+
+options_.solve_algo = 12;
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+ref = load('simul_backward_ref.mat');
+
+if max(max(abs(s.data - ref.s.data))) > 5e-4
+    error('Incorrect results for simul_backward with solve_algo=12')
+end
diff --git a/tests/solve_algo_12_14/simul_backward_14.mod b/tests/solve_algo_12_14/simul_backward_14.mod
new file mode 100644
index 000000000..8816f9360
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_14.mod
@@ -0,0 +1,12 @@
+// Check the correctedness of “simul_backward” with “solve_algo=14”
+
+@#include "simul_backward_common.inc"
+
+options_.solve_algo = 14;
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+ref = load('simul_backward_ref.mat');
+
+if max(max(abs(s.data - ref.s.data))) > 5e-4
+    error('Incorrect results for simul_backward with solve_algo=14')
+end
diff --git a/tests/solve_algo_12_14/simul_backward_common.inc b/tests/solve_algo_12_14/simul_backward_common.inc
new file mode 100644
index 000000000..86f7fabcc
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_common.inc
@@ -0,0 +1,15 @@
+@#include "backward_model.inc"
+
+@#define simlength = 100
+
+initialconditions = dseries(rand(2, @{5*n}), '2021Q3', { ...
+@#for i in 1:n
+  'x@{i}', 'y@{i}', 'z@{i}', 't@{i}', 's@{i}', ...
+@#endfor
+});
+
+innovations = dseries(randn(@{simlength}, @{5*n}), '2022Q1', { ...
+@#for i in 1:n
+  'e_x@{i}', 'e_y@{i}', 'e_z@{i}', 'e_t@{i}', 'e_s@{i}', ...
+@#endfor
+});
diff --git a/tests/solve_algo_12_14/simul_backward_reference.mod b/tests/solve_algo_12_14/simul_backward_reference.mod
new file mode 100644
index 000000000..cf5284140
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_reference.mod
@@ -0,0 +1,7 @@
+// Reference simulation with “simul_backward” and default “solve_algo” value
+
+@#include "simul_backward_common.inc"
+
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+save simul_backward_ref.mat s