diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 05d20af691ee238380d258fadd288023224694f5..8d38504c3bbe71dad11095f20833e47d36be6906 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, M_, options_)
 
 % Core function calling solvers for perfect foresight model
 %
@@ -17,7 +17,7 @@ function [y, success, maxerror, iter, per_block_status] = perfect_foresight_solv
 % - iter                [integer] Number of iterations of the underlying nonlinear solver (empty for non-iterative methods)
 % - per_block_status    [struct] In the case of block decomposition, provides per-block solver status information (empty if no block decomposition)
 
-% Copyright © 2015-2024 Dynare Team
+% Copyright © 2015-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -55,6 +55,25 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
     options_.linear_approximation = true;
 end
 
+if ~isempty(M_.perfect_foresight_controlled_paths)
+    assert(nargout >= 6); % Ensure modified exos are used
+    if options_.bytecode
+        error('perfect_foresight_controlled_paths is not compatible with bytecode option')
+    end
+    if options_.block
+        error('perfect_foresight_controlled_paths is not compatible with block option')
+    end
+    if options_.linear
+        error('perfect_foresight_controlled_paths is not compatible with linear option')
+    end
+    if ~ismember(options_.stack_solve_algo, [0 2 3])
+        error('perfect_foresight_controlled_paths is only available with stack_solve_algo equal to 0, 2 or 3')
+    end
+    if M_.maximum_endo_lead == 0 || M_.maximum_endo_lag == 0
+        error('perfect_foresight_controlled_paths is not compatible with purely backward, purely forward or static models')
+    end
+end
+
 maxerror = [];
 iter = [];
 per_block_status = [];
@@ -103,7 +122,7 @@ else
                     end
                     [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, M_, options_);
                 end
               case {1 6}
                 if options_.linear_approximation
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 55f5f6f63fd658761d110a57119de0ae52dbe7bd..3cb5d7a0f939496d39b5c544d23bcad475445215 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, 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.
@@ -15,8 +15,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 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 +41,19 @@ ny = M_.endo_nbr;
 periods = get_simulation_periods(options_);
 vperiods = periods*ones(1,options_.simul.maxit);
 
+if ~isempty(M_.perfect_foresight_controlled_paths)
+    for i=1:length(M_.perfect_foresight_controlled_paths)
+        endogenousvariables(M_.perfect_foresight_controlled_paths(i).exogeneize_id, M_.perfect_foresight_controlled_paths(i).periods + M_.maximum_lag) = M_.perfect_foresight_controlled_paths(i).value;
+    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 +85,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(M_.perfect_foresight_controlled_paths)
+        A = stacked_jacobian_substitute_controlled_paths(A, y, y0, yT, exogenousvariables, steadystate, periods, M_);
+    end
+
     if options_.debug && iter==1
         [row,col]=find(A);
         row=setdiff(1:periods*ny,row);
@@ -135,7 +155,16 @@ for iter = 1:options_.simul.maxit
             display_critical_variables(reshape(dy,[ny periods])', M_, options_.noprint);
         end
     end
-    y = y + dy;
+        y = y + dy;
+    if ~isempty(M_.perfect_foresight_controlled_paths)
+        for i=1:length(M_.perfect_foresight_controlled_paths)
+            p = M_.perfect_foresight_controlled_paths(i).periods;
+            endogenize_id = M_.perfect_foresight_controlled_paths(i).endogenize_id;
+            exogenize_id = M_.perfect_foresight_controlled_paths(i).exogenize_id;
+            y(exogenize_id + (p-1)*M_.endo_nbr) = M_.perfect_foresight_controlled_paths(i).values;
+            exogenousvariables(p,endogenize_id)=exogenousvariables(p,endogenize_id)+dy(exogenize_id + (p-1)*M_.endo_nbr);
+        end
+    end
 end
 
 endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, ny, periods);
@@ -278,7 +307,7 @@ end
 
 function display_critical_variables(dyy, M_, noprint)
 
-if noprint
+if noprint || ~isempty(M_.perfect_foresight_controlled_paths)
     return
 end
 
diff --git a/matlab/perfect-foresight-models/stacked_jacobian_substitute_controlled_paths.m b/matlab/perfect-foresight-models/stacked_jacobian_substitute_controlled_paths.m
new file mode 100644
index 0000000000000000000000000000000000000000..60be7aebc58da8ae7af70026ecfdc3d866101554
--- /dev/null
+++ b/matlab/perfect-foresight-models/stacked_jacobian_substitute_controlled_paths.m
@@ -0,0 +1,48 @@
+function A = stacked_jacobian_substitute_controlled_paths(A, y, y0, yT, exogenousvariables, steadystate, periods, M_)
+
+% 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/>.
+
+% Reorganize M_.perfect_foresight_controlled_paths period by period
+flips_by_period = struct('exogeneize_id', cell(periods, 1), 'endogeneize_id', cell(periods, 1));
+for i=1:length(M_.perfect_foresight_controlled_paths)
+    for p=1:M_.perfect_foresight_controlled_paths(i).periods
+        flips_by_period(p).exogeneize_id(end+1) = M_.perfect_foresight_controlled_paths(i).exogeneize_id;
+        flips_by_period(p).endogeneize_id(end+1) = M_.perfect_foresight_controlled_paths(i).endogeneize_id;
+    end
+end
+
+dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
+
+for p = 1:periods
+    nsubst = length(flips_by_periods(p).exogeneize_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_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+flips_by_periods(p).endogeneize_id);
+                  zeros((periods-p)*M_.endo_nbr, nsubst) ];
+        A(:,(p-1)*M_.endo_nbr+flips_by_periods(p).exogeneize_id) = subst;
+    end
+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 0000000000000000000000000000000000000000..b3a229e7842a878d5002f0bc1f1f5341187e9954
--- /dev/null
+++ b/tests/deterministic_simulations/ramst_pf_controlled_paths.mod
@@ -0,0 +1,51 @@
+var c k;
+varexo x;
+
+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);
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady;
+
+check;
+
+perfect_foresight_controlled_paths;
+  exogenize c;
+  periods 2 3:4;
+  values 1.6 1.7;
+  endogenize x;
+ 
+  exogenize k;
+  periods 5:7;
+  values 13;
+  endogenize x;
+end;
+
+shocks;
+var x;
+periods 1;
+values 1.2;
+end;
+
+
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+rplot c;
+rplot k;