From 2a8abf4968514465046a24c722dc843911c6db24 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 19 Feb 2025 16:27:51 +0100
Subject: [PATCH] =?UTF-8?q?New=20=E2=80=9Cperfect=5Fforesight=5Fcontrolled?=
 =?UTF-8?q?=5Fpaths=E2=80=9D=20block?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 doc/manual/source/the-model-file.rst          | 138 ++++++++++++++++++
 license.txt                                   |   2 +-
 matlab/ep/extended_path.m                     |   4 +
 matlab/ep/extended_path_core.m                |   2 +-
 matlab/ep/extended_path_homotopy.m            |   6 +-
 .../controlled_paths_by_period.m              | 117 +++++++++++++++
 ...rolled_paths_substitute_stacked_jacobian.m |  43 ++++++
 .../perfect_foresight_setup.m                 |  12 ++
 .../perfect_foresight_solver.m                |  65 +++++++--
 .../perfect_foresight_solver_core.m           |  28 +++-
 ..._foresight_with_expectation_errors_setup.m |  10 +-
 ...foresight_with_expectation_errors_solver.m |  12 +-
 matlab/perfect-foresight-models/sim1.m        |  46 +++++-
 matlab/perfect-foresight-models/sim1_lbj.m    |  41 +++++-
 meson.build                                   |  27 +++-
 preprocessor                                  |   2 +-
 scripts/dynare.el                             |   6 +-
 .../homotopy_endval_steady_linearization.mod  |  29 +---
 .../homotopy_linearization.inc                |  28 ++++
 .../homotopy_linearization_condpaths.mod      |  25 ++++
 .../homotopy_marginal_linearization.mod       |  29 +---
 .../homotopy_marglinearization_condpaths.mod  |  25 ++++
 .../pfwee_homotopy.inc                        |  37 +++++
 .../pfwee_homotopy_linearization.mod          |  31 +---
 ...pfwee_homotopy_linearization_condpaths.mod |  22 +++
 .../pfwee_homotopy_marginal_linearization.mod |  37 +----
 .../pfwee_homotopy_marglin_condpaths.mod      |  22 +++
 .../ramst_pf_controlled_paths.mod             |  18 +++
 .../ramst_pf_controlled_paths_check.mod       |  35 +++++
 .../ramst_pf_controlled_paths_common.inc      |  50 +++++++
 .../ramst_pf_controlled_paths_lbj.mod         |  17 +++
 .../ramst_pfwee_controlled_paths.mod          |  29 ++++
 32 files changed, 838 insertions(+), 157 deletions(-)
 create mode 100644 matlab/perfect-foresight-models/controlled_paths_by_period.m
 create mode 100644 matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
 create mode 100644 tests/deterministic_simulations/homotopy_linearization.inc
 create mode 100644 tests/deterministic_simulations/homotopy_linearization_condpaths.mod
 create mode 100644 tests/deterministic_simulations/homotopy_marglinearization_condpaths.mod
 create mode 100644 tests/deterministic_simulations/pfwee_homotopy.inc
 create mode 100644 tests/deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod
 create mode 100644 tests/deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod
 create mode 100644 tests/deterministic_simulations/ramst_pf_controlled_paths.mod
 create mode 100644 tests/deterministic_simulations/ramst_pf_controlled_paths_check.mod
 create mode 100644 tests/deterministic_simulations/ramst_pf_controlled_paths_common.inc
 create mode 100644 tests/deterministic_simulations/ramst_pf_controlled_paths_lbj.mod
 create mode 100644 tests/deterministic_simulations/ramst_pfwee_controlled_paths.mod

diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index c52076ca1d..7b1a62eba4 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -4616,6 +4616,144 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
     the terminal condition for exogenous indexed ``k``, as anticipated from
     period ``s``, is stored in ``oo_.pfwee.terminal_info(k,s)``.
 
+
+Controlling the path of endogenous variables
+--------------------------------------------
+
+In the usual perfect foresight problem, the user controls the path of exogenous
+variables for the simulation periods and the initial and terminal
+conditions for endogenous variables, while Dynare solves for the path of
+endogenous variables for the simulation periods.
+
+However, Dynare offers the possibility of controlling the value of some
+endogenous variables for some simulation periods (in which case some exogenous
+variables must be left free and are thus solved for by Dynare, to avoid
+over-determination of the problem). This exercise is called “conditional
+forecasting” in some contexts (even though one may argue that this is not
+really forecasting, since perfect foresight by the agents is assumed; for the
+stochastic case, see the :comm:`conditional_forecast` command).
+
+The description of controlled endogenous variables is done using the
+``perfect_foresight_controlled_paths`` block. The information given therein is
+then processed by the :comm:`perfect_foresight_setup` (or
+:comm:`perfect_foresight_with_expectation_errors_setup`) command, so
+that the next :comm:`perfect_foresight_solver` (or
+:comm:`perfect_foresight_with_expectation_errors_solver`) command
+computes the simulation with controlled paths. In particular,
+:mvar:`oo_.exo_simul` will contain the computed value of exogenous variables
+that have been left free.
+
+.. block:: perfect_foresight_controlled_paths ;
+           perfect_foresight_controlled_paths(OPTIONS...);
+
+    |br| This block is used to tell the perfect foresight solver that the value
+    of some endogenous variables will be controlled (in other words, they will
+    be exogenized). It also gives the period(s) for which this control applies,
+    the value(s) imposed to the endogenous variable(s), and the exogenous
+    variable(s) that are left free at the same period(s) (in other words,
+    those exogenous are endogenized).
+
+    The block should contain one or more occurrences of the following
+    group of four lines::
+
+      exogenize ENDOGENOUS_NAME;
+      periods INTEGER[:INTEGER] | DATE[:DATE] [[,] INTEGER[:INTEGER] | DATE[:DATE]]...;
+      values DOUBLE | (EXPRESSION)  [[,] DOUBLE | (EXPRESSION) ]...;
+      endogenize EXOGENOUS_NAME;
+
+    Note that it is possible to have both
+    ``perfect_foresight_controlled_paths`` and regular :bck:`shocks` blocks in
+    the same ``.mod`` file (assuming of course that taken together they do not
+    impose inconsistent constraints).
+
+    The ``perfect_foresight_controlled_paths`` block requires that the
+    :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` option be equal to
+    either ``0``, ``1``, ``2``, ``3`` or ``6``, and is incompatible with the
+    :opt:`block`, :opt:`linear` and :opt:`bytecode` options of the :bck:`model`
+    block.
+
+    *Options*
+
+    .. option:: learnt_in = INTEGER | DATE
+
+       Used in conjunction with
+       :comm:`perfect_foresight_with_expectation_errors_setup` and
+       :comm:`perfect_foresight_with_expectation_errors_solver` commands,
+       specifies the period or date at which this controlled paths block is
+       learnt by agents.
+
+    *Example (perfect foresight)*
+
+        ::
+
+            var c k;
+            varexo x z;
+
+            ...
+
+            shocks;
+              var x;
+              periods 1;
+              values 1.2;
+            end;
+
+            perfect_foresight_controlled_paths;
+              exogenize c;
+              periods 2 4:5;
+              values 1.6 1.7;
+              endogenize x;
+
+              exogenize k;
+              periods 7:9;
+              values 13;
+              endogenize z;
+            end;
+
+	    perfect_foresight_setup(periods = 100);
+	    perfect_foresight_solver;
+
+        In this example, the exogenous variable ``x`` is equal to 1.2 in
+        period 1, but in periods 2, 4 and 5 it will be endogenized so that
+        endogenous variable `c` is equal to 1.6 in period 2 and then 1.7 in
+        periods 4 and 5. Similarly, the exogenous variable ``z`` will be
+        endogenized in periods 7 to 9 so that the endogenous variable ``k`` is
+        equal to 13 over the same periods.
+
+    *Example (perfect foresight with expectation errors)*
+
+        ::
+
+            var c;
+            varexo x;
+
+            ...
+
+            perfect_foresight_controlled_paths;
+              exogenize c;
+              periods 2002Y 2003Y:2005Y;
+              values 1.6 1.7;
+              endogenize x;
+            end;
+
+            perfect_foresight_controlled_paths(learnt_in=2004Y);
+              exogenize c;
+              periods 2004Y:2005Y;
+              values 1.8;
+              endogenize x;
+            end;
+
+            perfect_foresight_with_expectation_errors_setup(periods = 30,
+                first_simulation_period = 2001Y);
+            perfect_foresight_with_expectation_errors_solver;
+
+        In this example, agents in year 2001 (at beginning of the simulation)
+        compute their plan under the assumption that the endogenous variable
+        ``c`` will be equal to 1.6 in year 2002 and 1.7 from years 2003 to
+        2005, and that exogenous variable ``x`` will behave so as to fulfill
+        that constraint. Then, when 2004 arrives, they recompute their plan
+        under the assumption that ``c`` will be equal to 1.8 in years 2004 and
+        2005 (and again that ``x`` will be endogenized accordingly).
+
 .. _stoch-sol:
 
 Stochastic solution and simulation
diff --git a/license.txt b/license.txt
index 5b65abe896..6f5ea57bce 100644
--- a/license.txt
+++ b/license.txt
@@ -278,7 +278,7 @@ License: GPL-3+
 
 Files: scripts/dynare.el
 Copyright: 2010 Yannick Kalantzis
-           2019-2024 Dynare Team
+           2019-2025 Dynare Team
 License: GPL-3+
 
 Files: mex/sources/gensylv/gensylv.cc
diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 1072698f31..078ef9b3bf 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -36,6 +36,10 @@ function [ts,oo_] = extended_path(initialconditions, samplesize, exogenousvariab
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+if ~isempty(M_.perfect_foresight_controlled_paths)
+    error('extended_path command is not compatible with perfect_foresight_controlled_paths block')
+end
+
 [initialconditions, innovations, pfm, options_, oo_] = ...
     extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
 
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index ec0aa4736b..174b2867fe 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -76,7 +76,7 @@ if order == 0
     oo_.steady_state = steady_state;
     options_.solve_algo = solve_algo;
     options_.stack_solve_algo = stack_solve_algo;
-    [endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+    [endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
 else
     % Stochastic Extended Path
     switch(algo)
diff --git a/matlab/ep/extended_path_homotopy.m b/matlab/ep/extended_path_homotopy.m
index c4df1c1592..e653ad9a14 100644
--- a/matlab/ep/extended_path_homotopy.m
+++ b/matlab/ep/extended_path_homotopy.m
@@ -1,6 +1,6 @@
 function [info_convergence, endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M_, options_, oo_, pfm, ep, order, algo, method, debug)
 
-% Copyright © 2016-2023 Dynare Team
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -35,7 +35,7 @@ if ismember(method, [1, 2])
         oo_.endo_simul(:,1) = oo_.steady_state + weight*(endo_simul0(:,1) - oo_.steady_state);
         oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
         if order==0
-            [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+            [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
         else
             switch(algo)
               case 0
@@ -102,7 +102,7 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
         oo_.endo_simul = endo_simul;
         oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
         if order==0
-            [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+            [endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
         else
             switch(algo)
               case 0
diff --git a/matlab/perfect-foresight-models/controlled_paths_by_period.m b/matlab/perfect-foresight-models/controlled_paths_by_period.m
new file mode 100644
index 0000000000..ee2e64c9d9
--- /dev/null
+++ b/matlab/perfect-foresight-models/controlled_paths_by_period.m
@@ -0,0 +1,117 @@
+function controlled_paths_by_period = controlled_paths_by_period(M_, options_, info_period)
+% Reorganize M_.perfect_foresight_controlled_paths period by period.
+% “info_period” should only be set in a perfect foresight with expectation errors context,
+% in which case only the information available at that period is extracted.
+
+% Copyright © 2025 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/>.
+
+[periods, first_simulation_period] = get_simulation_periods(options_);
+
+controlled_paths_by_period = struct('exogenize_id', cell(periods, 1), 'endogenize_id', cell(periods, 1), 'values', cell(periods, 1), 'learnt_in', cell(periods, 1));
+for i=1:length(M_.perfect_foresight_controlled_paths)
+    learnt_in = M_.perfect_foresight_controlled_paths(i).learnt_in;
+    if isa(learnt_in, 'dates')
+        learnt_in = learnt_in - first_simulation_period + 1;
+    end
+    if nargin >= 3 && learnt_in > info_period
+        continue
+    end
+
+    prange = M_.perfect_foresight_controlled_paths(i).periods;
+    if isa(prange, 'dates')
+        prange = transpose(prange - first_simulation_period + 1);
+    end
+    if nargin >= 3 && any(prange < learnt_in)
+        error('perfect_foresight_controlled_paths(learnt_in=%d): the periods are inconsistent with the learnt_in value', learnt_in)
+    end
+    for p = prange
+        % Determine whether this is an change in expectations, and if yes overwrite the previous expected value
+        idx = find(controlled_paths_by_period(p).exogenize_id == M_.perfect_foresight_controlled_paths(i).exogenize_id ...
+                   & controlled_paths_by_period(p).endogenize_id == M_.perfect_foresight_controlled_paths(i).endogenize_id ...
+                   & controlled_paths_by_period(p).learnt_in < learnt_in);
+        if isempty(idx)
+            controlled_paths_by_period(p).exogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).exogenize_id;
+            controlled_paths_by_period(p).endogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).endogenize_id;
+            controlled_paths_by_period(p).values(end+1) = M_.perfect_foresight_controlled_paths(i).value;
+            controlled_paths_by_period(p).learnt_in(end+1) = learnt_in;
+        else
+            controlled_paths_by_period(p).values(idx) = M_.perfect_foresight_controlled_paths(i).value;
+            controlled_paths_by_period(p).learnt_in(idx) = learnt_in;
+        end
+    end
+end
+
+
+%% Do various sanity checks
+
+for p = 1:periods
+    exogenize_id = controlled_paths_by_period(p).exogenize_id;
+    [~, idx] = unique(exogenize_id, 'stable');
+    duplicate_idx = setdiff(1:numel(exogenize_id), idx);
+    if ~isempty(duplicate_idx)
+        error('perfect_foresight_controlled_paths: variable %s is exogenized two times in period %d', M_.endo_names{exogenize_id(duplicate_idx(1))}, p);
+    end
+
+    endogenize_id = controlled_paths_by_period(p).endogenize_id;
+    [~, idx] = unique(endogenize_id, 'stable');
+    duplicate_idx = setdiff(1:numel(endogenize_id), idx);
+    if ~isempty(duplicate_idx)
+        error('perfect_foresight_controlled_paths: variable %s is endogenized two times in period %d', M_.exo_names{endogenize_id(duplicate_idx(1))}, p);
+    end
+end
+
+if isfield(M_, 'det_shocks')
+    for i = 1:length(M_.det_shocks)
+        if M_.det_shocks(i).exo_det
+            continue
+        end
+        prange = M_.det_shocks(i).periods;
+        if isa(prange, 'dates')
+            prange = transpose(prange - first_simulation_period + 1);
+        end
+        for p = prange
+            exo_id = M_.det_shocks(i).exo_id;
+            if ismember(exo_id, controlled_paths_by_period(p).endogenize_id)
+                error('Exogenous variable %s is present in both the shocks and the perfect_foresight_controlled_paths in period %d', M_.exo_names{exo_id}, p);
+            end
+        end
+    end
+end
+
+if nargin >= 3 && ~isempty(M_.learnt_shocks)
+    for i = 1:length(M_.learnt_shocks)
+        learnt_in = M_.learnt_shocks.learnt_in;
+        if isa(learnt_in, 'dates')
+            learnt_in = learnt_in - first_simulation_period + 1;
+        end
+        if learnt_in > info_period
+            continue
+        end
+
+        prange = M_.learnt_shocks(i).periods;
+        if isa(prange, 'dates')
+            prange = transpose(prange - first_simulation_period + 1);
+        end
+        for p = prange
+            exo_id = M_.learnt_shocks(i).exo_id;
+            if ismember(exo_id, controlled_paths_by_period(p).endogenize_id)
+                error('Exogenous variable %s is present in both the shocks(learnt_in=%d) and the perfect_foresight_controlled_paths in period %d', learnt_in, M_.exo_names{exo_id}, p);
+            end
+        end
+    end
+end
diff --git a/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m b/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
new file mode 100644
index 0000000000..1be87d245d
--- /dev/null
+++ b/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
@@ -0,0 +1,43 @@
+function A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, controlled_paths_by_period, M_)
+% Given a stacked Jacobian for the perfect foresight problem, modifies some columns
+% according to the perfect_foresight_controlled_paths block.
+
+% Copyright © 2025 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/>.
+
+dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
+
+periods = numel(controlled_paths_by_period);
+
+for p = 1:periods
+    nsubst = length(controlled_paths_by_period(p).exogenize_id);
+    if nsubst > 0
+        if p == 1
+            y3n = [ y0; y(1:2*M_.endo_nbr) ];
+        elseif p == periods
+            y3n = [ y(end-2*M_.endo_nbr+1:end); yT ];
+        else
+            y3n = y((p-2)*M_.endo_nbr+(1:3*M_.endo_nbr));
+        end
+        g1 = dynamic_g1(y3n, exogenousvariables(p+M_.maximum_lag,:), M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr);
+
+        subst = [ zeros((p-1)*M_.endo_nbr, nsubst);
+                  g1(:,3*M_.endo_nbr+controlled_paths_by_period(p).endogenize_id);
+                  zeros((periods-p)*M_.endo_nbr, nsubst) ];
+        A(:,(p-1)*M_.endo_nbr+controlled_paths_by_period(p).exogenize_id) = subst;
+    end
+end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index 287316777f..fda2a81f82 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -78,3 +78,15 @@ end
 
 oo_ = make_ex_(M_,options_,oo_);
 oo_ = make_y_(M_,options_,oo_);
+
+if isempty(M_.perfect_foresight_controlled_paths)
+    oo_.deterministic_simulation.controlled_paths_by_period = [];
+else
+    for i=1:length(M_.perfect_foresight_controlled_paths)
+        learnt_in = M_.perfect_foresight_controlled_paths(i).learnt_in;
+        if ~isa(learnt_in, 'numeric') || ~isequal(learnt_in, 1)
+            error('A perfect_foresight_controlled_paths(learnt_in=...) block is present. You want to call perfect_foresight_with_expectations_error_setup and perfect_foresight_with_expectations_error_solver.')
+        end
+    end
+    oo_.deterministic_simulation.controlled_paths_by_period = controlled_paths_by_period(M_, options_);
+end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index d351f2864b..14e12e715e 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -149,7 +149,10 @@ if completed_share == 1
     if options_.simul.endval_steady
         oo_.steady_state = steady_state;
     end
-    % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock
+    % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock, except if controlled paths
+    if ~isempty(oo_.deterministic_simulation.controlled_paths_by_period)
+        oo_.exo_simul = exo_simul;
+    end
 
     if ~options_.noprint
         fprintf('Perfect foresight solution found.\n\n')
@@ -169,7 +172,10 @@ elseif options_.simul.homotopy_linearization_fallback && completed_share > 0
         % This is not a true steady state, but it is the closest we can get to
         oo_.steady_state = oo_.endo_simul(:, end);
     end
-    % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock (although with an approximation)
+    % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock (although with an approximation), except if controlled paths
+    if ~isempty(oo_.deterministic_simulation.controlled_paths_by_period)
+        oo_.exo_simul = exobase + (exo_simul - exobase)/completed_share;
+    end
 
     maxerror = recompute_maxerror(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
 
@@ -201,9 +207,9 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_
         fprintf('%s\n\n', repmat('*', 1, 80))
     end
     extra_simul_time_counter = tic;
-    [extra_success, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = create_scenario(M_,options_,oo_,extra_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state);
+    [extra_success, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state, extra_controlled_paths_by_period] = create_scenario(M_,options_,oo_,extra_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, exo_simul, steady_state, exo_steady_state);
     if extra_success
-        [extra_endo_simul, extra_success] = perfect_foresight_solver_core(extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state, M_, options_);
+        [extra_endo_simul, extra_success, ~, ~, ~, extra_exo_simul] = perfect_foresight_solver_core(extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state, extra_controlled_paths_by_period, M_, options_);
     end
     if ~extra_success
         if ~options_.noprint
@@ -228,7 +234,10 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_
             % This is not a true steady state, but it is the closest we can get to
             oo_.steady_state = oo_.endo_simul(:, end);
         end
-        % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock (although with an approximation)
+        % NB: no need to modify oo_.exo_simul and oo_.exo_steady_state, since we simulated 100% of the shock (although with an approximation), except if controlled paths
+        if ~isempty(oo_.deterministic_simulation.controlled_paths_by_period)
+            oo_.exo_simul = exo_simul + (exo_simul - extra_exo_simul)*(1-completed_share)/options_.simul.homotopy_marginal_linearization_fallback;
+        end
 
         maxerror = recompute_maxerror(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
 
@@ -299,6 +308,7 @@ success_counter = 0;
 iteration = 0;
 
 endo_simul = endoorig;
+exo_simul = exoorig;
 periods = get_simulation_periods(options_);
 
 while step > options_.simul.homotopy_min_step_size
@@ -306,6 +316,7 @@ while step > options_.simul.homotopy_min_step_size
     iteration = iteration+1;
 
     saved_endo_simul = endo_simul;
+    saved_exo_simul = exo_simul;
 
     new_share = completed_share + step; % Try this share, and see if it succeeds
 
@@ -316,7 +327,7 @@ while step > options_.simul.homotopy_min_step_size
 
     iter_time_counter = tic;
 
-    [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(M_,options_,oo_,new_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state);
+    [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period] = create_scenario(M_,options_,oo_,new_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, exo_simul, steady_state, exo_steady_state);
 
     if steady_success
         % At the first iteration, use the initial guess given by
@@ -332,10 +343,21 @@ while step > options_.simul.homotopy_min_step_size
             else
                 endo_simul(:, simperiods) = endobase(:, simperiods);
             end
+
+            % Also handle initial guess in exo_simul when there are controlled paths:
+            % first try the initial guess given by perfect_foresight_setup or the user,
+            % afterwards use the base scenario
+            if ~isempty(oo_.deterministic_simulation.controlled_paths_by_period)
+                if iteration == 1 && new_share == shareorig
+                    % Nothing to do, at this point exo_simul(:, simperiods) == exoorig(:, simperiods)
+                else
+                    exo_simul(simperiods,:) = exobase(simperiods,:);
+                end
+            end
         end
 
         % Solve for the paths of the endogenous variables.
-        [endo_simul, success, maxerror, solver_iter, per_block_status] = perfect_foresight_solver_core(endo_simul, exo_simul, steady_state, exo_steady_state, M_, options_);
+        [endo_simul, success, maxerror, solver_iter, per_block_status, exo_simul] = perfect_foresight_solver_core(endo_simul, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period, M_, options_);
     else
         success = false;
         maxerror = NaN;
@@ -378,6 +400,7 @@ while step > options_.simul.homotopy_min_step_size
         end
     else
         endo_simul = saved_endo_simul;
+        exo_simul = saved_exo_simul;
         success_counter = 0;
         step = step / 2;
         if ~options_.noprint
@@ -411,7 +434,7 @@ end
 fprintf('\n')
 
 
-function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(M_,options_,oo_,share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state)
+function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period] = create_scenario(M_,options_,oo_,share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, exo_simul, steady_state, exo_steady_state)
 % For a given share, comutes the exogenous path and also the initial and
 % terminal conditions for the endogenous path (but do not modify the initial
 % guess for endogenous)
@@ -424,12 +447,13 @@ function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state]
 %   shareorig        [double]    the share to which endoorig and exoorig correspond (typically 100%, except for perfect_foresight_with_expectation_errors_solver with homotopy and marginal linearization)
 %   endoorig         [matrix]    path of endogenous corresponding to shareorig of the shock (only initial and terminal conditions are used)
 %   exoorig          [matrix]    path of exogenous corresponding to shareorig of the shock
-%   endobase         [matrix]    path of endogenous corresponding to 0% of the shock (only initial and terminal conditions are used)
+%   endobase         [matrix]    path of endogenous corresponding to 0% of the shock (only initial and terminal conditions are used, except if oo_.deterministic_simulation.controlled_paths_by_period is not empty)
 %   exobase          [matrix]    path of exogenous corresponding to 0% of the shock
 %   initperiods      [vector]    period indices of initial conditions
 %   lastperiods      [vector]    period indices of terminal conditions
 %   recompute_final_steady_state [boolean] self-explanatory
 %   endo_simul       [matrix]    path of endogenous, used to construct the guess values (initial and terminal conditions are not used)
+%   exo_simul        [matrix]    path of exogenous, used to construct the guess values (only if oo_.deterministic_simulation.controlled_paths_by_period is not empty)
 %   steady_state     [vector]    steady state of endogenous, only used if terminal steady state is *not* recomputed by the function
 %   exo_steady_state [vector]    steady state of exogenous, only used if terminal steady state is *not* recomputed by the function
 %
@@ -439,6 +463,9 @@ function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state]
 %   exo_simul        [matrix]    path of exogenous corresponding to the scenario
 %   steady_state     [vector]    steady state of endogenous corresponding to the scenario (equal to the input if terminal steady state not recomputed)
 %   exo_steady_state [vector]    steady state of exogenous corresponding to the scenario (equal to the input if terminal steady state not recomputed)
+%   controlled_paths_by_period [struct] flips between endos and exos corresponding to the scenario
+
+saved_exo_simul = exo_simul; % For guess value of controlled paths
 
 % Compute convex combination for the path of exogenous
 exo_simul = exoorig*share/shareorig + exobase*(1-share/shareorig);
@@ -509,6 +536,26 @@ else
     endo_simul(:, lastperiods) = share/shareorig*endoorig(:, lastperiods)+(1-share/shareorig)*endobase(:, lastperiods);
 end
 
+controlled_paths_by_period = oo_.deterministic_simulation.controlled_paths_by_period;
+if ~isempty(controlled_paths_by_period)
+    for p = 1:length(controlled_paths_by_period)
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        controlled_paths_by_period(p).values = controlled_paths_by_period(p).values*share + endobase(controlled_paths_by_period(p).exogenize_id,p)'*(1-share);
+
+        % Handle guess values for endogenized exos
+        exo_ids = controlled_paths_by_period(p).endogenize_id;
+        if ~isempty(options_.simul.homotopy_exclude_varexo)
+            [is_excluded, excluded_exo_ids] = ismember(options_.simul.homotopy_exclude_varexo, M_.exo_names{exo_ids});
+            if any(is_excluded)
+                error('Exogenous %s cannot be in the homotopy_exclude_varexo option and a perfect_foresight_controlled_paths block at the same time', M_.exo_names{excluded_exo_ids(1)})
+            end
+        end
+        exo_simul(length(initperiods)+p,exo_ids) = saved_exo_simul(length(initperiods)+p,exo_ids);
+    end
+end
+
 
 function maxerror = recompute_maxerror(endo_simul, exo_simul, steady_state, M_, options_)
     % Computes ∞-norm of residuals for a given path of endogenous,
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index c561e59be8..c7e376f29f 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -1,4 +1,4 @@
-function [y, success, maxerror, iter, per_block_status] = perfect_foresight_solver_core(y, exo_simul, steady_state, exo_steady_state, M_, options_)
+function [y, success, maxerror, iter, per_block_status, exo_simul] = perfect_foresight_solver_core(y, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period, M_, options_)
 
 % Core function calling solvers for perfect foresight model
 %
@@ -7,6 +7,9 @@ function [y, success, maxerror, iter, per_block_status] = perfect_foresight_solv
 % - exo_simul           [matrix] path of exogenous
 % - steady_state        [vector] steady state of endogenous variables
 % - exo_steady_state    [vector] steady state of exogenous variables
+% - controlled_paths_by_period [struct] describes flips between endos and exos, typically
+%                          extracted from perfect_foresight_controlled_paths block;
+%                          Can be an empty array if no such flip is requested
 % - M_                  [struct] contains a description of the model.
 % - options_            [struct] contains various options.
 %
@@ -73,6 +76,25 @@ if options_.endogenous_terminal_period
     end
 end
 
+if ~isempty(controlled_paths_by_period)
+    assert(nargout >= 6); % Ensure modified exos are used
+    if ~ismember(options_.stack_solve_algo, [0 1 2 3 6])
+        error('perfect_foresight_controlled_paths is only available with stack_solve_algo option equal to 0, 1, 2, 3 or 6')
+    end
+    if options_.bytecode
+        error('perfect_foresight_controlled_paths is not available with the bytecode option')
+    end
+    if options_.block
+        error('perfect_foresight_controlled_paths is not available with the block option')
+    end
+    if options_.linear
+        error('perfect_foresight_controlled_paths is not available with the linear option')
+    end
+    if M_.maximum_endo_lead == 0 || M_.maximum_endo_lag == 0
+        error('perfect_foresight_controlled_paths is not available with purely backward, purely forward or static models')
+    end
+end
+
 
 if options_.linear && ismember(options_.stack_solve_algo, [0 7]) && ~options_.block ...
         && ~options_.bytecode && M_.maximum_endo_lead > 0 && M_.maximum_endo_lag > 0
@@ -124,10 +146,10 @@ else
                 if options_.linear_approximation
                     [y, success, maxerror] = sim1_linear(y, exo_simul, steady_state, exo_steady_state, M_, options_);
                 else
-                    [y, success, maxerror, iter] = sim1(y, exo_simul, steady_state, M_, options_);
+                    [y, success, maxerror, iter, exo_simul] = sim1(y, exo_simul, steady_state, controlled_paths_by_period, M_, options_);
                 end
               case {1 6}
-                [y, success, maxerror, iter] = sim1_lbj(y, exo_simul, steady_state, M_, options_);
+                [y, success, maxerror, iter, exo_simul] = sim1_lbj(y, exo_simul, steady_state, controlled_paths_by_period, M_, options_);
               case 7
                 if options_.linear_approximation
                     if options_.solve_algo == 10 || options_.solve_algo == 11
diff --git a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
index fce1048a21..529c00f869 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_setup.m
@@ -68,8 +68,8 @@ if exist(options_.datafile, 'file')
     end
 else
     %% No datafile option given, use the contents of shocks and endval blocks
-    if isempty(M_.learnt_shocks) && isempty(M_.learnt_endval)
-        warning('perfect_foresight_with_expectation_errors_setup: there is no shocks(learnt_in=...) or endval(learnt_in=...) block, and you did not pass the datafile option, so there is no point in using this command')
+    if isempty(M_.learnt_shocks) && isempty(M_.learnt_endval) && (isempty(M_.perfect_foresight_controlled_paths) || all(cellfun(@(x) (isa(x, 'numeric') && x == 1), {M_.perfect_foresight_controlled_paths.learnt_in})))
+        warning('perfect_foresight_with_expectation_errors_setup: there is no shocks(learnt_in=...), endval(learnt_in=...), or perfect_foresight_controlled_paths(learnt_in=...) block, and you did not pass the datafile option, so there is no point in using this command')
     end
 
     %% Check that dates can be processed, if any
@@ -168,6 +168,12 @@ else
     end
 end
 
+% Handle controlled paths
+oo_.pfwee.controlled_paths_by_period = struct('exogenize_id', cell(periods, periods), 'endogenize_id', cell(periods, periods), 'values', cell(periods, periods), 'learnt_in', cell(periods, periods)); % 1st dimension of cells is real time, 2nd is informational time
+for p = 1:periods
+    oo_.pfwee.controlled_paths_by_period(:, p) = controlled_paths_by_period(M_, options_, p);
+end
+
 % Build initial paths for endos and exos (only initial conditions are set, the rest is NaN)
 if isempty(oo_.initial_steady_state)
     oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+periods+M_.maximum_lead);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
index 22996f1cd1..54f75a6aa8 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m
@@ -28,6 +28,10 @@ if options_.pfwee.constant_simulation_length && ~isempty(options_.simul.last_sim
     error('Options constant_simulation_length and last_simulation_period cannot be used together')
 end
 
+if options_.pfwee.constant_simulation_length && ~all(cellfun(@(x) isempty(x), {oo_.pfwee.controlled_paths_by_period.endogenize_id}))
+    error('Options constant_simulation_length is incompatible with perfect_foresight_controlled_paths block')
+end
+
 [periods, first_simulation_period] = get_simulation_periods(options_);
 
 % Retrieve initial paths built by pfwee_setup
@@ -69,6 +73,11 @@ while info_period <= periods
     oo_.exo_simul(M_.maximum_lag+(1:periods-info_period+1), :) = oo_.pfwee.shocks_info(:, info_period:end, info_period)';
     oo_.exo_simul(M_.maximum_lag+periods-info_period+2:end, :) = repmat(oo_.exo_steady_state', sim_length+M_.maximum_lead-(periods-info_period+1), 1);
 
+    oo_.deterministic_simulation.controlled_paths_by_period = oo_.pfwee.controlled_paths_by_period(info_period:end, info_period);
+    if all(cellfun(@(x) isempty(x), {oo_.deterministic_simulation.controlled_paths_by_period.endogenize_id}))
+        oo_.deterministic_simulation.controlled_paths_by_period = []; % Expected by pf_solver when there is no controlled var
+    end
+
     options_.periods = sim_length;
     % The following two options are reset to empty, so as to avoid an inconsistency with periods
     options_.simul.first_simulation_period = dates();
@@ -108,7 +117,8 @@ while info_period <= periods
     increment = 1;
     while info_period+increment <= periods && ...
           all(oo_.pfwee.terminal_info(:, info_period) == oo_.pfwee.terminal_info(:, info_period+increment)) && ...
-          all(all(oo_.pfwee.shocks_info(:, info_period+increment:end, info_period) == oo_.pfwee.shocks_info(:, info_period+increment:end, info_period+increment)))
+          all(all(oo_.pfwee.shocks_info(:, info_period+increment:end, info_period) == oo_.pfwee.shocks_info(:, info_period+increment:end, info_period+increment))) && ...
+          isequal(oo_.pfwee.controlled_paths_by_period(info_period+increment:end, info_period), oo_.pfwee.controlled_paths_by_period(info_period+increment:end, info_period+increment))
         increment = increment + 1;
     end
     info_period = info_period + increment;
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 55f5f6f63f..e360a38dc6 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
+function [endogenousvariables, success, err, iter, exogenousvariables] = sim1(endogenousvariables, exogenousvariables, steadystate, controlled_paths_by_period, M_, options_)
 % [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
 % Performs deterministic simulations with lead or lag of one period, using
 % a basic Newton solver on sparse matrices.
@@ -6,8 +6,9 @@ function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, e
 %
 % INPUTS
 %   - endogenousvariables [double] N*(T+M_.maximum_lag+M_.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
-%   - exogenousvariables  [double] T*M array, paths for the exogenous variables.
+%   - exogenousvariables  [double] (T+M_.maximum_lag+M_.maximum_lead)*M array, paths for the exogenous variables.
 %   - steadystate         [double] N*1 array, steady state for the endogenous variables.
+%   - controlled_paths_by_period [struct] data from perfect_foresight_controlled_paths block
 %   - M_                  [struct] contains a description of the model.
 %   - options_            [struct] contains various options.
 % OUTPUTS
@@ -15,8 +16,10 @@ function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, e
 %   - success             [logical] Whether a solution was found
 %   - err                 [double] ∞-norm of the residual
 %   - iter                [integer] Number of iterations
+%   - exogenousvariables  [double] (T+M_.maximum_lag+M_.maximum_lead)*M array, paths for the exogenous variables
+%                                  (may be modified if perfect_foresight_controlled_paths present)
 
-% Copyright © 1996-2024 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,6 +42,22 @@ ny = M_.endo_nbr;
 periods = get_simulation_periods(options_);
 vperiods = periods*ones(1,options_.simul.maxit);
 
+if ~isempty(controlled_paths_by_period)
+    for p = 1:periods
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        endogenousvariables(controlled_paths_by_period(p).exogenize_id, p+M_.maximum_lag) = controlled_paths_by_period(p).values;
+    end
+
+    if options_.debug
+        error('Debugging not available with perfect_foresight_controlled_paths')
+    end
+    if options_.endogenous_terminal_period
+        error('The endogenous_terminal_period option not available with perfect_foresight_controlled_paths')
+    end
+end
+
 if M_.maximum_lag > 0
     y0 = endogenousvariables(:, M_.maximum_lag);
 else
@@ -70,6 +89,11 @@ for iter = 1:options_.simul.maxit
     [res, A] = perfect_foresight_problem(y, y0, yT, exogenousvariables, M_.params, steadystate, periods, M_, options_);
     % A is the stacked Jacobian with period x equations alongs the rows and
     % periods times variables (in declaration order) along the columns
+
+    if ~isempty(controlled_paths_by_period)
+        A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, controlled_paths_by_period, M_);
+    end
+
     if options_.debug && iter==1
         [row,col]=find(A);
         row=setdiff(1:periods*ny,row);
@@ -132,10 +156,22 @@ for iter = 1:options_.simul.maxit
     end
     if any(isnan(dy)) || any(isinf(dy))
         if verbose
-            display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint);
+            display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint || ~isempty(controlled_paths_by_period));
         end
     end
     y = y + dy;
+
+    if ~isempty(controlled_paths_by_period)
+        for p = 1:periods
+            endogenize_id = controlled_paths_by_period(p).endogenize_id;
+            exogenize_id = controlled_paths_by_period(p).exogenize_id;
+            if isempty(endogenize_id)
+                continue
+            end
+            y(exogenize_id+(p-1)*M_.endo_nbr) = controlled_paths_by_period(p).values;
+            exogenousvariables(p+M_.maximum_lag,endogenize_id) = exogenousvariables(p+M_.maximum_lag,endogenize_id) + dy(exogenize_id+(p-1)*M_.endo_nbr)';
+        end
+    end
 end
 
 endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, ny, periods);
@@ -154,7 +190,7 @@ if stop
             skipline()
             fprintf('Total time of simulation: %g.\n', etime(clock,h1))
             disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
-            display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint);
+            display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint || ~isempty(controlled_paths_by_period));
             disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
             printline(105)
         end
diff --git a/matlab/perfect-foresight-models/sim1_lbj.m b/matlab/perfect-foresight-models/sim1_lbj.m
index 491f93cd62..5c7442797d 100644
--- a/matlab/perfect-foresight-models/sim1_lbj.m
+++ b/matlab/perfect-foresight-models/sim1_lbj.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M_, options_)
+function [endogenousvariables, success, err, iter, exogenousvariables] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, controlled_paths_by_period, M_, options_)
 
 % Performs deterministic simulations with lead or lag on one period using the historical LBJ algorithm
 %
@@ -10,6 +10,8 @@ function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariable
 %   success             [logical]       Whether a solution was found
 %   err                 [double]        ∞-norm of Δendogenousvariables
 %   iter                [integer]       Number of iterations
+%   exogenousvariables  [matrix]        All exogenous variables of the model
+%                                       (may be modified if perfect_foresight_controlled_paths present)
 %
 % ALGORITHM
 %   Laffargue, Boucekkine, Juillard (LBJ)
@@ -20,7 +22,7 @@ function [endogenousvariables, success, err, iter] = sim1_lbj(endogenousvariable
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 1996-2024 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -72,6 +74,19 @@ function [r, g1] = dynamicmodel(it_)
         [r, T_order, T] = dynamic_resid(y3n, x, M_.params, steadystate);
         g1_sparse = dynamic_g1(y3n, x, M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr, T_order, T);
         g1 = full(g1_sparse(:, find(lead_lag_incidence')));
+
+        if ~isempty(controlled_paths_by_period)
+            p = it_ - M_.maximum_lag;
+            if p > 1 && ~isempty(controlled_paths_by_period(p-1).exogenize_id)
+                g1(:,nonzeros(lead_lag_incidence(1,controlled_paths_by_period(p-1).exogenize_id))) = 0;
+            end
+            if ~isempty(controlled_paths_by_period(p).exogenize_id)
+                g1(:,nyp+controlled_paths_by_period(p).exogenize_id) = g1_sparse(:,3*ny+controlled_paths_by_period(p).endogenize_id);
+            end
+            if p < periods && ~isempty(controlled_paths_by_period(p+1).exogenize_id)
+                g1(:,nonzeros(lead_lag_incidence(3,controlled_paths_by_period(p+1).exogenize_id))) = 0;
+            end
+        end
     end
 end
 
@@ -83,6 +98,15 @@ if verbose
     fprintf('MODEL SIMULATION :\n')
 end
 
+if ~isempty(controlled_paths_by_period)
+    for p = 1:periods
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        endogenousvariables(controlled_paths_by_period(p).exogenize_id,p+M_.maximum_lag) = controlled_paths_by_period(p).values;
+    end
+end
+
 h1 = clock;
 
 for iter = 1:options_.simul.maxit
@@ -103,6 +127,19 @@ for iter = 1:options_.simul.maxit
     end
     c = back_subst_lbj(c, ny, iyf, periods);
     endogenousvariables(:,M_.maximum_lag+(1:periods)) = endogenousvariables(:,M_.maximum_lag+(1:periods))+c;
+
+    if ~isempty(controlled_paths_by_period)
+        for p = 1:periods
+            endogenize_id = controlled_paths_by_period(p).endogenize_id;
+            exogenize_id = controlled_paths_by_period(p).exogenize_id;
+            if isempty(endogenize_id)
+                continue
+            end
+            endogenousvariables(exogenize_id,p+M_.maximum_lag) = controlled_paths_by_period(p).values;
+            exogenousvariables(p+M_.maximum_lag,endogenize_id) = exogenousvariables(p+M_.maximum_lag,endogenize_id) + c(exogenize_id,p)';
+        end
+    end
+
     err = max(max(abs(c)));
     if verbose
         fprintf('Iter: %s,\t err. = %s, \t time = %s\n', num2str(iter), num2str(err), num2str(etime(clock, h2)));
diff --git a/meson.build b/meson.build
index f6b2a7c01a..05a4f6b6af 100644
--- a/meson.build
+++ b/meson.build
@@ -1399,8 +1399,15 @@ mod_and_m_tests = [
   { 'test' : [ 'deterministic_simulations/rbc_det6.mod' ] },
   { 'test' : [ 'deterministic_simulations/homotopy.mod' ] },
   { 'test' : [ 'deterministic_simulations/homotopy_histval.mod' ] },
-  { 'test' : [ 'deterministic_simulations/homotopy_endval_steady_linearization.mod' ] },
-  { 'test' : [ 'deterministic_simulations/homotopy_marginal_linearization.mod' ] },
+  { 'test' : [
+    'deterministic_simulations/homotopy_endval_steady_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_marginal_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_linearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
+  { 'test' : [ 'deterministic_simulations/homotopy_marglinearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/homotopy_linearization.inc' ] },
   { 'test' : [ 'deterministic_simulations/homotopy_exclude_varexo.mod' ] },
   { 'test' : [ 'deterministic_simulations/rbc_det_exo_lag_2a.mod',
                'deterministic_simulations/rbc_det_exo_lag_2b.mod',
@@ -1435,8 +1442,14 @@ mod_and_m_tests = [
   { 'test' : [ 'deterministic_simulations/pfwee_learnt_in_dseries_bare_1st_period.mod' ],
     'extra' : [ 'deterministic_simulations/pfwee_learnt_in.inc' ] },
   { 'test' : [ 'deterministic_simulations/pfwee_multiple_shocks.mod' ] },
-  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization.mod' ] },
-  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marginal_linearization.mod' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marginal_linearization.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
+  { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod' ],
+    'extra' : [ 'deterministic_simulations/pfwee_homotopy.inc' ] },
   { 'test' : [ 'deterministic_simulations/lbj/rbc.mod' ] },
   { 'test' : [ 'lmmcp/rbc.mod' ] },
   { 'test' : [ 'lmmcp/sw_lmmcp.mod',
@@ -1612,6 +1625,12 @@ mod_and_m_tests = [
                'deprecated/ramsey_ex.mod' ] },
   { 'test' : [ 'deterministic_simulations/ramst.mod',
                'deprecated/ramst.mod' ] },
+  { 'test' : [ 'deterministic_simulations/ramst_pf_controlled_paths.mod',
+               'deterministic_simulations/ramst_pf_controlled_paths_lbj.mod',
+               'deterministic_simulations/ramst_pf_controlled_paths_check.mod' ],
+    'extra' : [ 'deterministic_simulations/ramst_pf_controlled_paths_common.inc'] },
+  { 'test' : [ 'deterministic_simulations/ramst_pfwee_controlled_paths.mod' ],
+    'extra' : [ 'deterministic_simulations/ramst_pf_controlled_paths_common.inc'] },
 
   # ECB modfiles
   { 'test' : [ 'var-expectations/1/example1.mod' ] },
diff --git a/preprocessor b/preprocessor
index 9372796b05..5da05127ca 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 9372796b0552b442f2f363cdc987001f4824c779
+Subproject commit 5da05127ca47c7221d5f047a497f6ec8e5c6ce16
diff --git a/scripts/dynare.el b/scripts/dynare.el
index e723bf6b78..50d6e1d252 100644
--- a/scripts/dynare.el
+++ b/scripts/dynare.el
@@ -1,7 +1,7 @@
 ;;; dynare.el --- major mode for editing Dynare mod files
 
 ;; Copyright © 2010 Yannick Kalantzis
-;; Copyright © 2019-2024 Dynare Team
+;; Copyright © 2019-2025 Dynare Team
 ;;
 ;; This program is free software; you can redistribute it and/or modify
 ;; it under the terms of the GNU General Public License as published by
@@ -88,7 +88,7 @@
   '("stderr" "values" "periods" "scales" "restriction" "exclusion"
     "upper_cholesky" "lower_cholesky" "equation" "bind" "relax" "error_bind"
     "error_relax" "add" "multiply" "target" "auxname_target_nonstationary"
-    "component" "growth" "auxname" "kind" "weights")
+    "component" "growth" "auxname" "kind" "weights" "exogenize" "endogenize")
   "Dynare statements-like keywords.")
 
 ;; Those keywords that makes the lexer enter the DYNARE_BLOCK start condition
@@ -106,7 +106,7 @@
       "conditional_forecast_paths" "svar_identification" "moment_calibration"
       "irf_calibration" "ramsey_constraints" "generate_irfs" "matched_moments"
       "occbin_constraints" "model_replace" "pac_target_info" "matched_irfs"
-      "matched_irfs_weights" "verbatim")
+      "matched_irfs_weights" "perfect_foresight_controlled_paths" "verbatim")
     "Dynare block keywords."))
 
 ;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file)
diff --git a/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod b/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
index 97b8fbae01..2dbb998a70 100644
--- a/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
+++ b/tests/deterministic_simulations/homotopy_endval_steady_linearization.mod
@@ -1,34 +1,7 @@
 // Example that triggers homotopy in perfect foresight simulation.
 // Tests the endval_steady, homotopy_linearization_fallback options, and a few more.
 
-var Consumption, Capital, LoggedProductivity;
-
-varexo LoggedProductivityInnovation;
-
-parameters beta, alpha, delta, rho;
-
-beta = .985;
-alpha = 1/3;
-delta = alpha/10;
-rho = .9;
-
-model;
-  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
-  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
-  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
-end;
-
-initval;
-  LoggedProductivityInnovation = 0;
-end;
-
-steady;
-
-endval;
-  LoggedProductivityInnovation = 1;
-  Consumption = 0.1;
-  Capital = 1;
-end;
+@#include "homotopy_linearization.inc"
 
 perfect_foresight_setup(periods=200, endval_steady);
 
diff --git a/tests/deterministic_simulations/homotopy_linearization.inc b/tests/deterministic_simulations/homotopy_linearization.inc
new file mode 100644
index 0000000000..2ce5608d7e
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy_linearization.inc
@@ -0,0 +1,28 @@
+var Consumption, Capital, LoggedProductivity;
+
+varexo LoggedProductivityInnovation;
+
+parameters beta, alpha, delta, rho;
+
+beta = .985;
+alpha = 1/3;
+delta = alpha/10;
+rho = .9;
+
+model;
+  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
+  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
+  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
+end;
+
+initval;
+  LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+endval;
+  LoggedProductivityInnovation = 1;
+  Consumption = 0.1;
+  Capital = 1;
+end;
diff --git a/tests/deterministic_simulations/homotopy_linearization_condpaths.mod b/tests/deterministic_simulations/homotopy_linearization_condpaths.mod
new file mode 100644
index 0000000000..529d23b0be
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy_linearization_condpaths.mod
@@ -0,0 +1,25 @@
+// Example that triggers homotopy in perfect foresight simulation.
+// Tests homotopy_linearization_fallback with perfect_foresight_controlled_paths
+
+@#include "homotopy_linearization.inc"
+
+perfect_foresight_controlled_paths;
+  exogenize LoggedProductivity;
+  periods 8:9;
+  values 5;
+  endogenize LoggedProductivityInnovation;
+end;
+
+perfect_foresight_setup(periods=200, endval_steady);
+
+perfect_foresight_solver(homotopy_max_completion_share = 0.7,
+                         homotopy_linearization_fallback,
+                         steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ~(all(abs(oo_.endo_simul(3, 9:10) - 5) < 1e-15) && abs(oo_.exo_simul(10) - 0.5) < 1e-15)
+   error('Homotopy with linearization and controlled paths failed')
+end
diff --git a/tests/deterministic_simulations/homotopy_marginal_linearization.mod b/tests/deterministic_simulations/homotopy_marginal_linearization.mod
index 080e11b67c..5202c708b1 100644
--- a/tests/deterministic_simulations/homotopy_marginal_linearization.mod
+++ b/tests/deterministic_simulations/homotopy_marginal_linearization.mod
@@ -1,34 +1,7 @@
 // Example that triggers homotopy in perfect foresight simulation.
 // Tests the homotopy_marginal_linearization_fallback and homotopy_step_size_increase_success_count options
 
-var Consumption, Capital, LoggedProductivity;
-
-varexo LoggedProductivityInnovation;
-
-parameters beta, alpha, delta, rho;
-
-beta = .985;
-alpha = 1/3;
-delta = alpha/10;
-rho = .9;
-
-model;
-  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
-  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
-  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
-end;
-
-initval;
-  LoggedProductivityInnovation = 0;
-end;
-
-steady;
-
-endval;
-  LoggedProductivityInnovation = 1;
-  Consumption = 0.1;
-  Capital = 1;
-end;
+@#include "homotopy_linearization.inc"
 
 perfect_foresight_setup(periods=200, endval_steady);
 
diff --git a/tests/deterministic_simulations/homotopy_marglinearization_condpaths.mod b/tests/deterministic_simulations/homotopy_marglinearization_condpaths.mod
new file mode 100644
index 0000000000..19c069725d
--- /dev/null
+++ b/tests/deterministic_simulations/homotopy_marglinearization_condpaths.mod
@@ -0,0 +1,25 @@
+// Example that triggers homotopy in perfect foresight simulation.
+// Tests homotopy_marginal_linearization_fallback with perfect_foresight_controlled_paths
+
+@#include "homotopy_linearization.inc"
+
+perfect_foresight_controlled_paths;
+  exogenize LoggedProductivity;
+  periods 8:9;
+  values 5;
+  endogenize LoggedProductivityInnovation;
+end;
+
+perfect_foresight_setup(periods=200, endval_steady);
+
+perfect_foresight_solver(homotopy_max_completion_share = 0.7,
+                         homotopy_marginal_linearization_fallback,
+                         steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ~(all(abs(oo_.endo_simul(3, 9:10) - 5) < 1e-14) && abs(oo_.exo_simul(10) - 0.5) < 1e-13)
+   error('Homotopy with marginal linearization and controlled paths failed')
+end
diff --git a/tests/deterministic_simulations/pfwee_homotopy.inc b/tests/deterministic_simulations/pfwee_homotopy.inc
new file mode 100644
index 0000000000..7f75759455
--- /dev/null
+++ b/tests/deterministic_simulations/pfwee_homotopy.inc
@@ -0,0 +1,37 @@
+var Consumption, Capital, LoggedProductivity;
+
+varexo LoggedProductivityInnovation;
+
+parameters beta, alpha, delta, rho;
+
+beta = .985;
+alpha = 1/3;
+delta = alpha/10;
+rho = .9;
+
+model;
+  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
+  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
+  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
+end;
+
+initval;
+  LoggedProductivityInnovation = 0;
+end;
+
+steady;
+
+endval;
+  LoggedProductivityInnovation = 0.6;
+end;
+
+endval(learnt_in = 5);
+  LoggedProductivityInnovation = 1;
+end;
+
+shocks(learnt_in = 3);
+  var LoggedProductivityInnovation;
+  periods 4;
+  values 0.2;
+end;
+
diff --git a/tests/deterministic_simulations/pfwee_homotopy_linearization.mod b/tests/deterministic_simulations/pfwee_homotopy_linearization.mod
index 72471830af..d4c89fa4b4 100644
--- a/tests/deterministic_simulations/pfwee_homotopy_linearization.mod
+++ b/tests/deterministic_simulations/pfwee_homotopy_linearization.mod
@@ -1,36 +1,7 @@
 /* Example that triggers homotopy in perfect foresight simulation with
    expectation errors, and tests linearization. */
 
-var Consumption, Capital, LoggedProductivity;
-
-varexo LoggedProductivityInnovation;
-
-parameters beta, alpha, delta, rho;
-
-beta = .985;
-alpha = 1/3;
-delta = alpha/10;
-rho = .9;
-
-model;
-  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
-  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
-  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
-end;
-
-initval;
-  LoggedProductivityInnovation = 0;
-end;
-
-steady;
-
-endval;
-  LoggedProductivityInnovation = 0.6;
-end;
-
-endval(learnt_in = 5);
-  LoggedProductivityInnovation = 1;
-end;
+@#include "pfwee_homotopy.inc"
 
 perfect_foresight_with_expectation_errors_setup(periods=200);
 perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_linearization_fallback, steady_solve_algo = 13);
diff --git a/tests/deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod b/tests/deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod
new file mode 100644
index 0000000000..e9a9eedbb3
--- /dev/null
+++ b/tests/deterministic_simulations/pfwee_homotopy_linearization_condpaths.mod
@@ -0,0 +1,22 @@
+/* Example that triggers homotopy in perfect foresight simulation with
+   expectation errors + controlled paths, and tests linearization. */
+
+@#include "pfwee_homotopy.inc"
+
+perfect_foresight_controlled_paths(learnt_in = 7);
+  exogenize LoggedProductivity;
+  periods 8:9;
+  values 5;
+  endogenize LoggedProductivityInnovation;
+end;
+
+perfect_foresight_with_expectation_errors_setup(periods=200);
+perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_linearization_fallback, steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ~(all(abs(oo_.endo_simul(3, 9:10) - 5) < 1e-15) && abs(oo_.exo_simul(10) - 0.5) < 1e-15)
+   error('Homotopy with linearization and controlled paths failed')
+end
diff --git a/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod b/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod
index 1151349d3f..b3fdfe20c3 100644
--- a/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod
+++ b/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod
@@ -1,42 +1,7 @@
 /* Example that triggers homotopy in perfect foresight simulation with
    expectation errors, and tests marginal linearization. */
 
-var Consumption, Capital, LoggedProductivity;
-
-varexo LoggedProductivityInnovation;
-
-parameters beta, alpha, delta, rho;
-
-beta = .985;
-alpha = 1/3;
-delta = alpha/10;
-rho = .9;
-
-model;
-  1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta);
-  Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption;
-  LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation;
-end;
-
-initval;
-  LoggedProductivityInnovation = 0;
-end;
-
-steady;
-
-endval;
-  LoggedProductivityInnovation = 0.6;
-end;
-
-endval(learnt_in = 5);
-  LoggedProductivityInnovation = 1;
-end;
-
-shocks(learnt_in = 3);
-  var LoggedProductivityInnovation;
-  periods 4;
-  values 0.2;
-end;
+@#include "pfwee_homotopy.inc"
 
 perfect_foresight_with_expectation_errors_setup(periods=200);
 perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_marginal_linearization_fallback, steady_solve_algo = 13);
diff --git a/tests/deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod b/tests/deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod
new file mode 100644
index 0000000000..cb64de238e
--- /dev/null
+++ b/tests/deterministic_simulations/pfwee_homotopy_marglin_condpaths.mod
@@ -0,0 +1,22 @@
+/* Example that triggers homotopy in perfect foresight simulation with
+   expectation errors + controlled paths, and tests marginal linearization. */
+
+@#include "pfwee_homotopy.inc"
+
+perfect_foresight_controlled_paths(learnt_in = 7);
+  exogenize LoggedProductivity;
+  periods 8:9;
+  values 5;
+  endogenize LoggedProductivityInnovation;
+end;
+
+perfect_foresight_with_expectation_errors_setup(periods=200);
+perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_marginal_linearization_fallback, steady_solve_algo = 13);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ~(all(abs(oo_.endo_simul(3, 9:10) - 5) < 1e-14) && abs(oo_.exo_simul(10) - 0.5) < 1e-14)
+   error('Homotopy with marginal linearization and controlled paths failed')
+end
diff --git a/tests/deterministic_simulations/ramst_pf_controlled_paths.mod b/tests/deterministic_simulations/ramst_pf_controlled_paths.mod
new file mode 100644
index 0000000000..84d382de21
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pf_controlled_paths.mod
@@ -0,0 +1,18 @@
+// Performs a simulation with a perfect_foresight_controlled_paths block
+// with the default stack_solve_algo=0
+
+@#define shocks = true
+@#include "ramst_pf_controlled_paths_common.inc"
+
+perfect_foresight_setup(periods=30, first_simulation_period = 2001Y);
+perfect_foresight_solver(stack_solve_algo = 0);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if oo_.endo_simul(1, 3) ~= 1.6 || any(oo_.endo_simul(1, 5:6) ~= 1.7) || any(oo_.endo_simul(2, 8:10) ~= 13) ...
+   || any(oo_.endo_simul(3, 4:6) ~= 0.5) || oo_.exo_simul(2, 1) ~= 1.2
+   error('Perfect foresight with controlled paths failed')
+end
+
diff --git a/tests/deterministic_simulations/ramst_pf_controlled_paths_check.mod b/tests/deterministic_simulations/ramst_pf_controlled_paths_check.mod
new file mode 100644
index 0000000000..b370f6001d
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pf_controlled_paths_check.mod
@@ -0,0 +1,35 @@
+/* Checks that exogenous computed by ramst_pf_controlled_paths{,lbj}.mod actually give back the
+   expected endogenous variables */
+
+@#define shocks = false
+@#include "ramst_pf_controlled_paths_common.inc"
+
+
+// stack_solve_algo=0
+
+S0 = load('ramst_pf_controlled_paths/Output/ramst_pf_controlled_paths_results.mat');
+
+perfect_foresight_setup(periods=30);
+
+oo_.exo_simul = S0.oo_.exo_simul;
+
+perfect_foresight_solver;
+
+if max(max(abs(oo_.endo_simul - S0.oo_.endo_simul))) > options_.dynatol.x
+    error('perfect_foresight_controlled_paths gives wrong result with stack_solve_algo=0')
+end
+
+
+// stack_solve_algo=1
+
+S1 = load('ramst_pf_controlled_paths_lbj/Output/ramst_pf_controlled_paths_lbj_results.mat');
+
+perfect_foresight_setup(periods=30);
+
+oo_.exo_simul = S1.oo_.exo_simul;
+
+perfect_foresight_solver;
+
+if max(max(abs(oo_.endo_simul - S1.oo_.endo_simul))) > options_.dynatol.x
+    error('perfect_foresight_controlled_paths gives wrong result with stack_solve_algo=1')
+end
diff --git a/tests/deterministic_simulations/ramst_pf_controlled_paths_common.inc b/tests/deterministic_simulations/ramst_pf_controlled_paths_common.inc
new file mode 100644
index 0000000000..e7680a65fd
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pf_controlled_paths_common.inc
@@ -0,0 +1,50 @@
+var c k y;
+varexo x e_y;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+y = 0.98*y(-1)+e_y;
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+y = 0;
+end;
+
+steady;
+
+@#if shocks
+perfect_foresight_controlled_paths;
+  exogenize c;
+  periods 2 4:5;
+  values 1.6 1.7;
+  endogenize x;
+
+  exogenize k;
+  periods 7:9;
+  values 13;
+  endogenize x;
+
+  exogenize y;
+  periods 2003Y:2005Y;
+  values 0.5;
+  endogenize e_y;
+end;
+
+shocks;
+  var x;
+  periods 1;
+  values 1.2;
+end;
+@#endif
diff --git a/tests/deterministic_simulations/ramst_pf_controlled_paths_lbj.mod b/tests/deterministic_simulations/ramst_pf_controlled_paths_lbj.mod
new file mode 100644
index 0000000000..76086a3649
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pf_controlled_paths_lbj.mod
@@ -0,0 +1,17 @@
+// Performs a simulation with a perfect_foresight_controlled_paths block
+// with the stack_solve_algo=1 (LBJ)
+
+@#define shocks = true
+@#include "ramst_pf_controlled_paths_common.inc"
+
+perfect_foresight_setup(periods=30, first_simulation_period = 2001Y);
+perfect_foresight_solver(stack_solve_algo = 1);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if oo_.endo_simul(1, 3) ~= 1.6 || any(oo_.endo_simul(1, 5:6) ~= 1.7) || any(oo_.endo_simul(2, 8:10) ~= 13) ...
+   || any(oo_.endo_simul(3, 4:6) ~= 0.5) || oo_.exo_simul(2, 1) ~= 1.2
+   error('Perfect foresight with controlled paths failed')
+end
diff --git a/tests/deterministic_simulations/ramst_pfwee_controlled_paths.mod b/tests/deterministic_simulations/ramst_pfwee_controlled_paths.mod
new file mode 100644
index 0000000000..d867d44409
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pfwee_controlled_paths.mod
@@ -0,0 +1,29 @@
+// Performs a pf_with_expectation_errors simulation with a perfect_foresight_controlled_paths block
+// In particular, tests that it’s possible to change the expected value of a given variable (c and y in period 4)
+
+@#define shocks = true
+@#include "ramst_pf_controlled_paths_common.inc"
+
+perfect_foresight_controlled_paths(learnt_in=2004Y);
+  exogenize c;
+  periods 4:6;
+  values 1.8;
+  endogenize x;
+
+  exogenize y;
+  periods 2004Y:2006Y;
+  values 0.6;
+  endogenize e_y;
+end;
+
+perfect_foresight_with_expectation_errors_setup(periods=30, first_simulation_period = 2001Y);
+perfect_foresight_with_expectation_errors_solver;
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if oo_.endo_simul(1, 3) ~= 1.6 || any(oo_.endo_simul(1, 5:7) ~= 1.8) || any(oo_.endo_simul(2, 8:10) ~= 13) ...
+   || oo_.endo_simul(3, 4) ~= 0.5 || any(oo_.endo_simul(3, 5:7) ~= 0.6) || oo_.exo_simul(2, 1) ~= 1.2
+   error('Perfect foresight with expectation errors and controlled paths failed')
+end
-- 
GitLab