diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index c4a590d0cb749aa5da5c5e458293c271dd9f1653..1072698f313d217a4866ab217e9afbdacf0dd1f5 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -68,23 +68,25 @@ while (t <= samplesize)
         initialguess = [];
     end
     if t>2
-        [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(innovations.positive_var_indx, ...
+        [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, y] = extended_path_core(innovations.positive_var_indx, ...
                                                                                                            spfm_exo_simul, ...
                                                                                                            endogenous_variables_paths(:,t-1), ...
                                                                                                            pfm, ...
                                                                                                            M_, ...
                                                                                                            options_, ...
                                                                                                            oo_, ...
-                                                                                                           initialguess);
+                                                                                                           initialguess, ...
+                                                                                                           y);
     else
-        [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, pfm, options_] = extended_path_core(innovations.positive_var_indx, ...
+        [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, y, pfm, options_] = extended_path_core(innovations.positive_var_indx, ...
                                                                                                                           spfm_exo_simul, ...
                                                                                                                           endogenous_variables_paths(:,t-1), ...
                                                                                                                           pfm, ...
                                                                                                                           M_, ...
                                                                                                                           options_, ...
                                                                                                                           oo_, ...
-                                                                                                                          initialguess);
+                                                                                                                          initialguess, ...
+                                                                                                                          []);
     end
     if ~info_convergence
         msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s)!', int2str(t));
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index 237d99b801ab3b5fbcf96ccb4f6529b4e3d70786..ec0aa4736ba8b650644af672e6ab931aae343310 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -1,11 +1,12 @@
-function [y, info_convergence, endogenousvariablespaths, pfm, options_] = extended_path_core(positive_var_indx, ...
+function [y1, info_convergence, endogenousvariablespaths, y, pfm, options_] = extended_path_core(positive_var_indx, ...
                                                                                              exo_simul, ...
                                                                                              initial_conditions ,...
                                                                                              pfm, ...
                                                                                              M_, ...
                                                                                              options_, ...
                                                                                              oo_, ...
-                                                                                             initialguess)
+                                                                                             initialguess, ...
+                                                                                             y)
 
 % Copyright © 2016-2025 Dynare Team
 %
@@ -41,7 +42,7 @@ stack_solve_algo = ep.stack_solve_algo;
 if init% Compute first order solution (Perturbation)...
     endo_simul = simult_(M_,options_,initial_conditions,oo_.dr,exo_simul(2:end,:),1);
 else
-    if nargin==8 && ~isempty(initialguess)
+    if nargin>7 && ~isempty(initialguess)
         % Note that the first column of initialguess should be equal to initial_conditions.
         endo_simul = initialguess;
     else
@@ -49,6 +50,10 @@ else
     end
 end
 
+if nargin~=9
+    y = [];
+end
+
 oo_.endo_simul = endo_simul;
 
 if debug
@@ -77,17 +82,17 @@ else
     switch(algo)
       case 0
         % Full tree of future trajectories.
-        if nargout>3
-            [flag, endogenousvariablespaths, errorcode, ~, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, options_, M_, pfm);
+        if nargout>4
+            [flag, endogenousvariablespaths, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm);
         else
-            [flag, endogenousvariablespaths, errorcode] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, options_, M_, pfm);
+            [flag, endogenousvariablespaths, errorcode, y] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm);
         end
       case 1
         % Sparse tree of future histories.
-        if nargout>3
-            [flag, endogenousvariablespaths, errorcode, ~, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, M_, pfm);
+        if nargout>4
+            [flag, endogenousvariablespaths, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm);
         else
-            [flag, endogenousvariablespaths, errorcode] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, M_, pfm);
+            [flag, endogenousvariablespaths, errorcode, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm);
         end
     end
     info_convergence = ~flag;
@@ -98,7 +103,7 @@ if ~info_convergence && ~options_.no_homotopy
 end
 
 if info_convergence
-    y = endogenousvariablespaths(:,2);
+    y1 = endogenousvariablespaths(:,2);
 else
-    y = NaN(size(endo_nbr,1));
+    y1 = NaN(size(endo_nbr,1));
 end
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
index 19e77a8315a6b2afdb96b742f01c66be4c9ed515..efdaf624721374a5855595a6c003ee0882145621 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
@@ -1,4 +1,4 @@
-function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, options_, M_, pfm)
+function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm)
 
 % Copyright © 2025 Dynare Team
 %
@@ -148,7 +148,9 @@ end
 
 pfm.Y = repmat(endo_simul(:),1,pfm.world_nbr);
 
-y = repmat(pfm.steady_state, pfm.dimension/pfm.ny, 1);
+if isempty(y)
+    y = repmat(pfm.steady_state, pfm.dimension/pfm.ny, 1);
+end
 
 if update_options_struct
     % Set algorithm
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
index ccdc882f61ae42003c44eb19e3761bf80d05bebf..19eb0673804beb1101b3a5fc3c0a15112e8694de 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
@@ -1,4 +1,4 @@
-function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, M_, pfm)
+function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm)
 
 % Copyright © 2012-2025 Dynare Team
 %
@@ -113,7 +113,9 @@ if update_pfm_struct
 
 end
 
-y = repmat(pfm.steady_state, pfm.block_nbr, 1);
+if isempty(y)
+    y = repmat(pfm.steady_state, pfm.block_nbr, 1);
+end
 
 if update_options_struct
     % Set algorithm