diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 7b0ef485c8438bcba2f2e7f191438c5c61867d2e..9422d030165c1b6e5d7714cc1fdd7db7f2c9ef65 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -329,6 +329,8 @@ options_.minimal_solving_periods = 1;
 options_.endogenous_terminal_period = false;
 options_.no_homotopy = false;
 options_.simul.endval_steady = false;
+options_.simul.first_simulation_period = dates();
+options_.simul.last_simulation_period = dates();
 
 options_.simul.homotopy_max_completion_share = 1;
 options_.simul.homotopy_min_step_size = 1e-3;
diff --git a/matlab/perfect-foresight-models/get_simulation_periods.m b/matlab/perfect-foresight-models/get_simulation_periods.m
new file mode 100644
index 0000000000000000000000000000000000000000..20eadd5d93df90819be16a57928e4cc885f9088d
--- /dev/null
+++ b/matlab/perfect-foresight-models/get_simulation_periods.m
@@ -0,0 +1,42 @@
+function [periods, first_simulation_period, last_simulation_period] = get_simulation_periods(options_)
+% Given the options_ structure, returns the number of simulation periods, and
+% also possibly the first and last simulation periods (as dates objects).
+% Also verify the consistency of the three options if all are given by the user.
+
+% Copyright © 2024 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/>.
+
+first_simulation_period = options_.simul.first_simulation_period;
+last_simulation_period = options_.simul.last_simulation_period;
+
+if options_.periods > 0
+    periods = options_.periods;
+    if ~isempty(first_simulation_period) && ~isempty(last_simulation_period)
+        if periods ~= last_simulation_period - first_simulation_period + 1
+            error('Options first_simulation_period, last_simulation_period and periods are inconsistent')
+        end
+    elseif ~isempty(first_simulation_period) && isempty(last_simulation_period)
+        last_simulation_period = first_simulation_period + periods - 1;
+    elseif isempty(first_simulation_period) && ~isempty(last_simulation_period)
+        first_simulation_period = last_simulation_period - periods + 1;
+    end
+else
+    if isempty(options_.simul.first_simulation_period) || isempty(options_.simul.last_simulation_period)
+        error('Dynare:periodsNotSet', 'Options first_simulation_period and last_simulation_period must be set when periods is not provided')
+    end
+    periods = last_simulation_period - first_simulation_period + 1;
+end
diff --git a/matlab/perfect-foresight-models/make_ex_.m b/matlab/perfect-foresight-models/make_ex_.m
index f3a9f483fb3c0a5b6ff8a182b8112c07bc1bd38d..e98f4c4b5dffb48d5644c97a66e1256ef0ebe20a 100644
--- a/matlab/perfect-foresight-models/make_ex_.m
+++ b/matlab/perfect-foresight-models/make_ex_.m
@@ -10,7 +10,7 @@ function oo_ = make_ex_(M_, options_, oo_)
 % OUTPUTS
 % - oo_          [struct]   Updated dynare results structure
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -27,6 +27,17 @@ function oo_ = make_ex_(M_, options_, oo_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+try
+    periods = get_simulation_periods(options_);
+catch ME
+    if strcmp(ME.identifier, 'Dynare:periodsNotSet')
+        % This function is called from dyn_forecast in some contexts where periods is not set
+        periods = 0;
+    else
+        rethrow(ME);
+    end
+end
+
 if isempty(oo_.exo_steady_state)
     oo_.exo_steady_state = zeros(M_.exo_nbr,1);
 end
@@ -38,13 +49,13 @@ end
 if isempty(oo_.initval_series)
     if isempty(M_.exo_histval)
         if isempty(oo_.initial_exo_steady_state)
-            oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
+            oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
         else
-            oo_.exo_simul = [ repmat(oo_.initial_exo_steady_state',M_.maximum_lag,1) ; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1) ];
+            oo_.exo_simul = [ repmat(oo_.initial_exo_steady_state',M_.maximum_lag,1) ; repmat(oo_.exo_steady_state',periods+M_.maximum_lead,1) ];
         end
     else
         if isempty(oo_.initial_exo_steady_state)
-            oo_.exo_simul = [M_.exo_histval'; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1)];
+            oo_.exo_simul = [M_.exo_histval'; repmat(oo_.exo_steady_state',periods+M_.maximum_lead,1)];
         else
             error('histval and endval cannot be used simultaneously')
         end
@@ -52,17 +63,17 @@ if isempty(oo_.initval_series)
 else
     if M_.exo_nbr > 0
         x = oo_.initval_series{M_.exo_names{:}}.data;
-        oo_.exo_simul = x(M_.orig_maximum_lag-M_.maximum_lag+1:M_.orig_maximum_lag + options_.periods + M_.maximum_lead,:);
+        oo_.exo_simul = x(M_.orig_maximum_lag-M_.maximum_lag+1:M_.orig_maximum_lag + periods + M_.maximum_lead,:);
         if ~isempty(M_.exo_histval)
             oo_.exo_simul(1:M_.maximum_lag, :) ...
                 = M_.exo_histval(:, 1:M_.maximum_lag)';
         end
     else
-        oo_.exo_simul=zeros(M_.maximum_lag + options_.periods + M_.maximum_lead,M_.exo_nbr);
+        oo_.exo_simul=zeros(M_.maximum_lag + periods + M_.maximum_lead,M_.exo_nbr);
     end
     if M_.exo_det_nbr > 0
         x_det = oo_.initval_series{M_.exo_det_names{:}}.data;
-        oo_.exo_det_simul = x_det(M_.orig_maximum_lag-M_.maximum_lag+1:M_.orig_maximum_lag + options_.periods + M_.maximum_lead,:);
+        oo_.exo_det_simul = x_det(M_.orig_maximum_lag-M_.maximum_lag+1:M_.orig_maximum_lag + periods + M_.maximum_lead,:);
         if ~isempty(M_.exo_det_histval)
             oo_.exo_det_simul(1:M_.maximum_lag, :) ...
                 = M_.exo_det_histval(:, 1:M_.maximum_lag)';
@@ -72,9 +83,9 @@ end
 % Initialize oo_.exo_det_simul
 if M_.exo_det_nbr > 0
     if isempty(M_.exo_det_histval)
-        oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
+        oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
     else
-        oo_.exo_det_simul = [M_.exo_det_histval'; repmat(oo_.exo_det_steady_state',options_.periods+M_.maximum_lead,1)];
+        oo_.exo_det_simul = [M_.exo_det_histval'; repmat(oo_.exo_det_steady_state',periods+M_.maximum_lead,1)];
     end
 end
 
diff --git a/matlab/perfect-foresight-models/make_y_.m b/matlab/perfect-foresight-models/make_y_.m
index c421421032d8bf6cc6c583a41b3f0a853db71b13..d653010b6e7b60d392379787071b96d30bc72a2e 100644
--- a/matlab/perfect-foresight-models/make_y_.m
+++ b/matlab/perfect-foresight-models/make_y_.m
@@ -10,7 +10,7 @@ function oo_=make_y_(M_, options_, oo_)
 % OUTPUTS
 % - oo_         [struct]   Updated dynare results structure
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -27,6 +27,8 @@ function oo_=make_y_(M_, options_, oo_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if isempty(oo_.steady_state)
     oo_.steady_state = zeros(M_.endo_nbr,1);
 end
@@ -34,9 +36,9 @@ end
 if isempty(oo_.initval_series)
     if isempty(M_.endo_histval)
         if isempty(oo_.initial_steady_state)
-            oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+options_.periods+M_.maximum_lead);
+            oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+periods+M_.maximum_lead);
         else
-            oo_.endo_simul = [repmat(oo_.initial_steady_state, 1, M_.maximum_lag) repmat(oo_.steady_state, 1,options_.periods+M_.maximum_lead)];
+            oo_.endo_simul = [repmat(oo_.initial_steady_state, 1, M_.maximum_lag) repmat(oo_.steady_state, 1,periods+M_.maximum_lead)];
         end
     else
         if ~isempty(oo_.initial_steady_state)
@@ -45,11 +47,11 @@ if isempty(oo_.initval_series)
         % the first NaNs take care of the case where there are lags > 1 on
         % exogenous variables
         oo_.endo_simul = [M_.endo_histval ...
-                          repmat(oo_.steady_state, 1, options_.periods+M_.maximum_lead)];
+                          repmat(oo_.steady_state, 1, periods+M_.maximum_lead)];
     end
 else
     y = oo_.initval_series{M_.endo_names{:}}.data;
-    oo_.endo_simul = y(M_.orig_maximum_lag - M_.maximum_lag + 1:M_.orig_maximum_lag + options_.periods + ...
+    oo_.endo_simul = y(M_.orig_maximum_lag - M_.maximum_lag + 1:M_.orig_maximum_lag + periods + ...
                        M_.maximum_lead, :)';
     if ~isempty(M_.endo_histval)
         if ~isempty(oo_.initial_steady_state)
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index cf943d001e3fd25404c686c06576bd91858037ef..bcbc17d9e90af24ad49c3f15aa3a50ea5ffcb79e 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -14,7 +14,7 @@ function oo_=perfect_foresight_setup(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -47,21 +47,19 @@ if size(M_.lead_lag_incidence,2)-nnz(M_.lead_lag_incidence(M_.maximum_endo_lag+1
     error(mess)
 end
 
-if options_.periods == 0
-    error('PERFECT_FORESIGHT_SETUP: number of periods for the simulation isn''t specified')
-end
+periods = get_simulation_periods(options_);
 
-if ~isempty(M_.det_shocks) && options_.periods<max([M_.det_shocks.periods])
+if ~isempty(M_.det_shocks) && periods < max([M_.det_shocks.periods])
     % Some expected shocks happen after the terminal period.
     mess = sprintf('\nPERFECT_FORESIGHT_SETUP: Problem with the declaration of the expected shocks:\n');
     for i=1:length(M_.det_shocks)
-        if any(M_.det_shocks(i).periods>options_.periods)
+        if any(M_.det_shocks(i).periods > periods)
             mess = sprintf('%s\n   At least one expected value for %s has been declared after the terminal period.', mess, M_.exo_names{M_.det_shocks(i).exo_id});
         end
     end
     disp(mess)
     skipline()
-    error('PERFECT_FORESIGHT_SETUP: Please check the declaration of the shocks or increase the value of the periods option.')
+    error('PERFECT_FORESIGHT_SETUP: Please check the declaration of the shocks or increase the number of periods in the simulation.')
 end
 
 if options_.simul.endval_steady && M_.maximum_lead == 0
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 63cbe0f2b4f652ec724ae11dd49d02a21dfe850d..5b5c71d2d1869ea28264ae2955d074e14f342ee6 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -52,7 +52,7 @@ if (~isempty(M_.learnt_shocks) || ~isempty(M_.learnt_endval)) && ~no_error_if_le
     error('A shocks(learnt_in=...) or endval(learnt_in=...) block is present. You want to call perfect_foresight_with_expectations_error_setup and perfect_foresight_with_expectations_error_solver.')
 end
 
-periods = options_.periods;
+[periods, first_simulation_period, last_simulation_period] = get_simulation_periods(options_);
 
 if options_.debug
     model_static = str2func([M_.fname,'.sparse.static_resid']);
@@ -274,15 +274,15 @@ if ~isempty(per_block_status)
     oo_.deterministic_simulation.block = per_block_status;
 end
 
-if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
-    initial_period = oo_.initval_series.dates(1)+(M_.orig_maximum_lag-1);
-elseif ~isdates(options_.initial_period) && isnan(options_.initial_period)
-    initial_period = dates(1,1);
-else
-    initial_period = options_.initial_period;
+if isempty(first_simulation_period)
+    if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
+        first_simulation_period = oo_.initval_series.dates(1)+(M_.orig_maximum_lag-1);
+    else
+        first_simulation_period = dates(1,1);
+    end
 end
 
-ts = dseries([transpose(oo_.endo_simul(1:M_.orig_endo_nbr,:)), oo_.exo_simul], initial_period, [M_.endo_names(1:M_.orig_endo_nbr); M_.exo_names]);
+ts = dseries([transpose(oo_.endo_simul(1:M_.orig_endo_nbr,:)), oo_.exo_simul], first_simulation_period - M_.maximum_lag, [M_.endo_names(1:M_.orig_endo_nbr); M_.exo_names]);
 
 if isfield(oo_, 'initval_series') && ~isempty(oo_.initval_series)
     names = ts.name;
@@ -319,6 +319,7 @@ success_counter = 0;
 iteration = 0;
 
 endo_simul = endoorig;
+periods = get_simulation_periods(options_);
 
 while step > options_.simul.homotopy_min_step_size
 
@@ -347,7 +348,7 @@ while step > options_.simul.homotopy_min_step_size
             if iteration == 1 && new_share == shareorig
                 % Nothing to do, at this point endo_simul(:, simperiods) == endoorig(:, simperiods)
             elseif M_.maximum_lead > 0
-                endo_simul(:, simperiods) = repmat(endo_simul(:, lastperiods(1)), 1, options_.periods);
+                endo_simul(:, simperiods) = repmat(endo_simul(:, lastperiods(1)), 1, periods);
             else
                 endo_simul(:, simperiods) = endobase(:, simperiods);
             end
@@ -555,6 +556,8 @@ function maxerror = recompute_maxerror(endo_simul, exo_simul, steady_state, M_,
 
 function check_input_arguments(options_, M_, oo_)
 
+periods = get_simulation_periods(options_);
+
 if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 7
     error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 7')
 end
@@ -567,10 +570,10 @@ if options_.block && ~options_.bytecode && options_.stack_solve_algo == 5
     error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 5 without bytecode option')
 end
 
-if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_lag+options_.periods+M_.maximum_lead ])
+if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_lag+periods+M_.maximum_lead ])
 
     if options_.initval_file
-        fprintf('PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Check whether your initval-file provides %d periods.',M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead)
+        fprintf('PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Check whether your initval-file provides %d periods.',M_.maximum_endo_lag+periods+M_.maximum_endo_lead)
         error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
     else
         error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
@@ -578,9 +581,9 @@ if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maxi
 end
 
 if (M_.exo_nbr > 0) && ...
-        (isempty(oo_.exo_simul) || any(size(oo_.exo_simul) ~= [ M_.maximum_lag+options_.periods+M_.maximum_lead, M_.exo_nbr ]))
+        (isempty(oo_.exo_simul) || any(size(oo_.exo_simul) ~= [ M_.maximum_lag+periods+M_.maximum_lead, M_.exo_nbr ]))
     if options_.initval_file
-        fprintf('PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size. Check whether your initval-file provides %d periods.',M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead)
+        fprintf('PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size. Check whether your initval-file provides %d periods.',M_.maximum_endo_lag+periods+M_.maximum_endo_lead)
         error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size.')
     else
         error('perfect_foresight_solver:ArgCheck','PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index fdeaa60e4b9a35dd2ef7f0fc93c414fc17619702..05d20af691ee238380d258fadd288023224694f5 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -41,7 +41,7 @@ elseif options_.stack_solve_algo==7 && options_.solve_algo == 11
     options_.lmmcp.status = 1; %Path solver
 end
 
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 
 if options_.linear_approximation && ~(isequal(options_.stack_solve_algo,0) || isequal(options_.stack_solve_algo,7))
     error('perfect_foresight_solver: Option linear_approximation is only available with option stack_solve_algo equal to 0 or 7.')
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 8074e7bebab1355306ffa3ecf061df148f9df214..545f566c5089cca9ca182a9d29c81f34c33d5734 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
@@ -7,7 +7,7 @@ function oo_=perfect_foresight_with_expectation_errors_setup(M_, options_, oo_)
 % OUTPUTS
 %   oo_                 [structure] storing the results
 
-% Copyright © 2021-2023 Dynare Team
+% Copyright © 2021-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -31,7 +31,7 @@ if ~isempty(oo_.initval_series)
     error('perfect_foresight_with_expectation_errors_setup: cannot be used in conjunction with histval_file/initval_file')
 end
 
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 
 %% Initialize informational structures
 oo_.pfwee.terminal_info = NaN(M_.exo_nbr, periods); % 2nd dimension is informational time
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 a7b63bfd5ae61741a9fddedc6e6d1a3aee1ff5c1..c5d9a10e39ee5d84277f69337b6d17576e3682cd 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
@@ -7,7 +7,7 @@ function oo_=perfect_foresight_with_expectation_errors_solver(M_, options_, oo_)
 % OUTPUTS
 %   oo_                 [structure] storing the results
 
-% Copyright © 2021-2023 Dynare Team
+% Copyright © 2021-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -24,8 +24,11 @@ function oo_=perfect_foresight_with_expectation_errors_solver(M_, options_, oo_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-% Same for periods (it will be modified before calling perfect_foresight_solver if constants_simulation_length option is false)
-periods = options_.periods;
+if options_.pfwee.constant_simulation_length && ~isempty(options_.simul.last_simulation_period)
+    error('Options constant_simulation_length and last_simulation_period cannot be used together')
+end
+
+periods = get_simulation_periods(options_);
 
 % Retrieve initial paths built by pfwee_setup
 % (the versions in oo_ will be truncated before calling perfect_foresight_solver)
@@ -67,6 +70,9 @@ while info_period <= periods
     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);
 
     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();
+    options_.simul.last_simulation_period = dates();
 
     if info_period > 1 && homotopy_completion_share < 1 && options_.simul.homotopy_marginal_linearization_fallback > 0
         marginal_linearization_previous_raw_sims.sim1.endo_simul = oo_.deterministic_simulation.sim1.endo_simul(:, info_period:end);
@@ -107,4 +113,4 @@ end
 
 % Set final paths
 oo_.endo_simul = endo_simul;
-oo_.exo_simul = exo_simul;
\ No newline at end of file
+oo_.exo_simul = exo_simul;
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index b2f189c31646475edc3a94ef2ef99026ffb37efd..55f5f6f63fd658761d110a57119de0ae52dbe7bd 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -36,7 +36,7 @@ function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, e
 verbose = options_.verbosity && ~options_.noprint;
 
 ny = M_.endo_nbr;
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 vperiods = periods*ones(1,options_.simul.maxit);
 
 if M_.maximum_lag > 0
@@ -141,7 +141,7 @@ end
 endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, ny, periods);
 
 if options_.endogenous_terminal_period
-    periods = options_.periods;
+    periods = get_simulation_periods(options_);
     err = evaluate_max_dynamic_residual(str2func([M_.fname,'.sparse.dynamic_resid']), endogenousvariables, exogenousvariables, M_.params, steadystate, periods, M_.maximum_lag);
 end
 
diff --git a/matlab/perfect-foresight-models/sim1_lbj.m b/matlab/perfect-foresight-models/sim1_lbj.m
index 4af6af8ba51e17080a48b00616aa4ae1bd5a9d4b..491f93cd6239814a68473cd0a864ff4451a8c3de 100644
--- a/matlab/perfect-foresight-models/sim1_lbj.m
+++ b/matlab/perfect-foresight-models/sim1_lbj.m
@@ -76,6 +76,7 @@ function [r, g1] = dynamicmodel(it_)
 end
 
 verbose = options_.verbosity;
+periods = get_simulation_periods(options_);
 
 if verbose
     printline(56)
@@ -86,13 +87,13 @@ h1 = clock;
 
 for iter = 1:options_.simul.maxit
     h2 = clock;
-    c = zeros(ny*options_.periods, nrc);
+    c = zeros(ny*periods, nrc);
     [d1, jacobian] = dynamicmodel(M_.maximum_lag+1);
     jacobian = [jacobian(:,iz), -d1];
     ic = 1:ny;
     icp = iyp;
     c (ic,:) = jacobian(:,is)\jacobian(:,isf1);
-    for it_ = M_.maximum_lag+(2:options_.periods)
+    for it_ = M_.maximum_lag+(2:periods)
         [d1, jacobian] = dynamicmodel(it_);
         jacobian = [jacobian(:,iz), -d1];
         jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:);
@@ -100,8 +101,8 @@ for iter = 1:options_.simul.maxit
         icp = icp + ny;
         c (ic,:) = jacobian(:,is)\jacobian(:,isf1);
     end
-    c = back_subst_lbj(c, ny, iyf, options_.periods);
-    endogenousvariables(:,M_.maximum_lag+(1:options_.periods)) = endogenousvariables(:,M_.maximum_lag+(1:options_.periods))+c;
+    c = back_subst_lbj(c, ny, iyf, periods);
+    endogenousvariables(:,M_.maximum_lag+(1:periods)) = endogenousvariables(:,M_.maximum_lag+(1:periods))+c;
     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/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m
index 26c52184cb3330231c53b0bc1700f3d377a9e30b..30c11700cf8e66d83fea7824d9fc3f4552c76c0a 100644
--- a/matlab/perfect-foresight-models/sim1_linear.m
+++ b/matlab/perfect-foresight-models/sim1_linear.m
@@ -63,7 +63,7 @@ nx = M_.exo_nbr;
 
 maximum_lag = M_.maximum_lag;
 
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 
 params = M_.params;
 
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index 2e2f50b325eed1c578cb50e2e12ba49730a14746..e005db254dc830314728ff784a47fd332ad2e6a5 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -2,7 +2,7 @@ function [endogenousvariables, success] = sim1_purely_backward(endogenousvariabl
 
 % Performs deterministic simulation of a purely backward model
 
-% Copyright © 2012-2023 Dynare Team
+% Copyright © 2012-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,6 +19,8 @@ function [endogenousvariables, success] = sim1_purely_backward(endogenousvariabl
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if ismember(options_.solve_algo, [12,14])
     [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
 else
@@ -36,7 +38,7 @@ end
 
 success = true;
 
-for it = M_.maximum_lag + (1:options_.periods)
+for it = M_.maximum_lag + (1:periods)
     y = endogenousvariables(:,it-1);        % Values at previous period, also used as guess value for current period
     x = exogenousvariables(it,:);
     if ismember(options_.solve_algo, [12,14])
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index b49d4180da695899a229714894e03cac9ff2ac9f..3a828506245fea85cffbc2060567d53c4426a193 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -1,7 +1,7 @@
 function [endogenousvariables, success] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M_, options_)
 % Performs deterministic simulation of a purely forward model
 
-% Copyright © 2012-2023 Dynare Team
+% Copyright © 2012-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -18,6 +18,8 @@ function [endogenousvariables, success] = sim1_purely_forward(endogenousvariable
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if ismember(options_.solve_algo, [12,14])
     [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
 else
@@ -35,7 +37,7 @@ end
 
 success = true;
 
-for it = options_.periods:-1:1
+for it = periods:-1:1
     yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
     x = exogenousvariables(it,:);
     if ismember(options_.solve_algo, [12,14])
diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m
index 1500485daa005dfeff878dbb075fc6c6a71c8152..a22c83b785efa0dc7703a1c6f385620f52391752 100644
--- a/matlab/perfect-foresight-models/sim1_purely_static.m
+++ b/matlab/perfect-foresight-models/sim1_purely_static.m
@@ -2,7 +2,7 @@ function [endogenousvariables, success] = sim1_purely_static(endogenousvariables
 
 % Performs deterministic simulation of a purely static model
 
-% Copyright © 2021-2023 Dynare Team
+% Copyright © 2021-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,6 +19,8 @@ function [endogenousvariables, success] = sim1_purely_static(endogenousvariables
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if ismember(options_.solve_algo, [12,14])
     [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_);
 else
@@ -38,7 +40,7 @@ success = true;
 
 y = endogenousvariables(:,1);
 
-for it = 1:options_.periods
+for it = 1:periods
     x = exogenousvariables(it,:);
     if ismember(options_.solve_algo, [12,14])
         T = NaN(M_.block_structure.dyn_tmp_nbr);
diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
index 3a48588cf302fe8add07cb3334bbacd0d2149e82..14db6dee693ee24108b535248d5c5b5147395361 100644
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
@@ -31,6 +31,7 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
 cutoff = 1e-15;
 
 switch options_.stack_solve_algo
@@ -55,7 +56,7 @@ if options_.verbosity
     skipline()
 end
 
-T=NaN(M_.block_structure.dyn_tmp_nbr, options_.periods+M_.maximum_lag+M_.maximum_lead);
+T=NaN(M_.block_structure.dyn_tmp_nbr, periods+M_.maximum_lag+M_.maximum_lead);
 
 maxerror = 0;
 nblocks = length(M_.block_structure.block);
@@ -71,9 +72,9 @@ for blk = 1:nblocks
     switch M_.block_structure.block(blk).Simulation_Type
         case {1, 2} % evaluate{Forward,Backward}
             if M_.block_structure.block(blk).Simulation_Type == 1
-                range = M_.maximum_lag+1:M_.maximum_lag+options_.periods;
+                range = M_.maximum_lag+1:M_.maximum_lag+periods;
             else
-                range = M_.maximum_lag+options_.periods:-1:M_.maximum_lag+1;
+                range = M_.maximum_lag+periods:-1:M_.maximum_lag+1;
             end
             for it_ = range
                 if it_ > 1 && it_ < size(y, 2)
@@ -97,7 +98,7 @@ for blk = 1:nblocks
         case {3, 4, 6, 7} % solve{Forward,Backward}{Simple,Complete}
             is_forward = M_.block_structure.block(blk).Simulation_Type == 3 || M_.block_structure.block(blk).Simulation_Type == 6;
             y_index = M_.block_structure.block(blk).variable(end-M_.block_structure.block(blk).mfs+1:end);
-            [y, T, success, maxblkerror, iter] = solve_one_boundary(fh_dynamic, y, exo_simul, M_.params, steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, options_.periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_);
+            [y, T, success, maxblkerror, iter] = solve_one_boundary(fh_dynamic, y, exo_simul, M_.params, steady_state, T, y_index, M_.block_structure.block(blk).NNZDerivatives, periods, M_.block_structure.block(blk).is_linear, blk, M_.maximum_lag, options_.simul.maxit, options_.dynatol.f, cutoff, options_.stack_solve_algo, is_forward, true, false, M_, options_);
         case {5, 8} % solveTwoBoundaries{Simple,Complete}
             if ismember(options_.stack_solve_algo, [1 6])
                 [y, T, success, maxblkerror, iter] = solve_two_boundaries_lbj(fh_dynamic, y, exo_simul, steady_state, T, blk, options_, M_);
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
index 2e919a53e9091686d4d2cab74a65bb7674995508..10ab54c114a678a45497c8d78d753a34b7319ec1 100644
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -17,17 +17,19 @@ function [endogenousvariables, success] = solve_stacked_linear_problem(endogenou
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if M_.maximum_lag > 0
     y0 = endogenousvariables(:, M_.maximum_lag);
 else
     y0 = NaN(M_.endo_nbr, 1);
 end
 if M_.maximum_lead > 0
-    yT = endogenousvariables(:, M_.maximum_lag+options_.periods+1);
+    yT = endogenousvariables(:, M_.maximum_lag+periods+1);
 else
     yT = NaN(M_.endo_nbr, 1);
 end
-z = endogenousvariables(:,M_.maximum_lag+(1:options_.periods));
+z = endogenousvariables(:,M_.maximum_lag+(1:periods));
 
 % Evaluate the residuals and Jacobian of the dynamic model at the deterministic steady state.
 y3n = repmat(steadystate_y, 3, 1);
@@ -50,7 +52,7 @@ x = bsxfun(@minus, exogenousvariables, steadystate_x');
                                            options_, ...
                                            jacobian, y0-steadystate_y, yT-steadystate_y, ...
                                            x, M_.params, steadystate_y, ...
-                                           M_.maximum_lag, options_.periods, M_.endo_nbr);
+                                           M_.maximum_lag, periods, M_.endo_nbr);
 
 if all(imag(y)<.1*options_.dynatol.x)
     if ~isreal(y)
@@ -60,7 +62,7 @@ else
     check = 1;
 end
 
-endogenousvariables = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,options_.periods), steadystate_y) yT];
+endogenousvariables = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,periods), steadystate_y) yT];
 
 success = ~check;
 
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index 4ec67beacf096d449248ba8d0998cf593898f238..61a7b3294136de34150b71205749cdee62da3974 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -31,17 +31,19 @@ function [endogenousvariables, success, maxerror] = solve_stacked_problem(endoge
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+periods = get_simulation_periods(options_);
+
 if M_.maximum_lag > 0
     y0 = endogenousvariables(:, M_.maximum_lag);
 else
     y0 = NaN(M_.endo_nbr, 1);
 end
 if M_.maximum_lead > 0
-    yT = endogenousvariables(:, M_.maximum_lag+options_.periods+1);
+    yT = endogenousvariables(:, M_.maximum_lag+periods+1);
 else
     yT = NaN(M_.endo_nbr, 1);
 end
-z = endogenousvariables(:,M_.maximum_lag+(1:options_.periods));
+z = endogenousvariables(:,M_.maximum_lag+(1:periods));
 
 if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
     [lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
@@ -50,11 +52,11 @@ if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementari
         ub = ub - steadystate_y;
     end
     if options_.solve_algo == 10
-        options_.lmmcp.lb = repmat(lb,options_.periods,1);
-        options_.lmmcp.ub = repmat(ub,options_.periods,1);
+        options_.lmmcp.lb = repmat(lb,periods,1);
+        options_.lmmcp.ub = repmat(ub,periods,1);
     elseif options_.solve_algo == 11
-        options_.mcppath.lb = repmat(lb,options_.periods,1);
-        options_.mcppath.ub = repmat(ub,options_.periods,1);
+        options_.mcppath.lb = repmat(lb,periods,1);
+        options_.mcppath.ub = repmat(ub,periods,1);
     end
     dynamic_resid_function = str2func([M_.fname,'.sparse.dynamic_resid']);
     dynamic_g1_function = str2func([M_.fname,'.sparse.dynamic_g1']);
@@ -64,14 +66,14 @@ if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementari
                                                  dynamic_resid_function, dynamic_g1_function, y0, yT, ...
                                                  exogenousvariables, M_.params, steadystate, ...
                                                  M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr, ...
-                                                 M_.maximum_lag, options_.periods, M_.endo_nbr, ...
+                                                 M_.maximum_lag, periods, M_.endo_nbr, ...
                                                  M_.dynamic_mcp_equations_reordering);
     eq_to_ignore=find(isfinite(lb) | isfinite(ub));
 
 else
     [y, check, res, ~, errorcode] = dynare_solve(@perfect_foresight_problem, z(:), ...
                                                options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, ...
-                                               options_, y0, yT, exogenousvariables, M_.params, steadystate, options_.periods, M_, options_);
+                                               options_, y0, yT, exogenousvariables, M_.params, steadystate, periods, M_, options_);
 end
 
 if all(imag(y)<.1*options_.dynatol.x)
@@ -82,9 +84,9 @@ else
     check = 1;
 end
 
-endogenousvariables(:, M_.maximum_lag+(1:options_.periods)) = reshape(y, M_.endo_nbr, options_.periods);
+endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, M_.endo_nbr, periods);
 residuals=zeros(size(endogenousvariables));
-residuals(:, M_.maximum_lag+(1:options_.periods)) = reshape(res, M_.endo_nbr, options_.periods);
+residuals(:, M_.maximum_lag+(1:periods)) = reshape(res, M_.endo_nbr, periods);
 if (options_.solve_algo == 10 || options_.solve_algo == 11)% mixed complementarity problem
     residuals(eq_to_ignore,bsxfun(@le, endogenousvariables(eq_to_ignore,:), lb(eq_to_ignore)+eps) | bsxfun(@ge,endogenousvariables(eq_to_ignore,:),ub(eq_to_ignore)-eps))=0;
 end
diff --git a/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m b/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m
index 752cc8665925505047616f99d2e5ed7907e38685..8ea326e22ac000c1e81634568cc868a65a6b2c2b 100644
--- a/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m
+++ b/matlab/perfect-foresight-models/solve_two_boundaries_lbj.m
@@ -25,7 +25,7 @@ function [y, T, success, err, iter] = solve_two_boundaries_lbj(fh, y, x, steady_
 %   simulation of dynamic models with forward variables through the use
 %   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.
 
-% Copyright © 2023 Dynare Team
+% Copyright © 2023-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -46,7 +46,7 @@ sparse_rowval = M_.block_structure.block(blk).g1_sparse_rowval;
 sparse_colval = M_.block_structure.block(blk).g1_sparse_colval;
 sparse_colptr = M_.block_structure.block(blk).g1_sparse_colptr;
 
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 
 % NB: notations are deliberately similar to those of sim1_lbj.m
 
diff --git a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
index 3fdee24facedafaf7844ddfb627da4a3864120d3..e3934e2579cefe5092d4b7f6f0bf59559cf56302 100644
--- a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
+++ b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
@@ -44,7 +44,7 @@ function [y, T, success, max_res, iter] = solve_two_boundaries_stacked(fh, y, x,
 
 Blck_size = M_.block_structure.block(Block_Num).mfs;
 y_index = M_.block_structure.block(Block_Num).variable(end-Blck_size+1:end);
-periods = options_.periods;
+periods = get_simulation_periods(options_);
 y_kmin = M_.maximum_lag;
 stack_solve_algo = options_.stack_solve_algo;
 
diff --git a/preprocessor b/preprocessor
index a992fb1b9c730583edef312f08177bc705d3dddc..35b68b1411ae18c1403882ab784897983c3193fd 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit a992fb1b9c730583edef312f08177bc705d3dddc
+Subproject commit 35b68b1411ae18c1403882ab784897983c3193fd
diff --git a/tests/deterministic_simulations/homotopy.mod b/tests/deterministic_simulations/homotopy.mod
index a1043102e83e815452538496cf165af529771fe4..e1e2dc10757f685f3a9c3fd7d27fc5259ff54e36 100644
--- a/tests/deterministic_simulations/homotopy.mod
+++ b/tests/deterministic_simulations/homotopy.mod
@@ -31,8 +31,6 @@ steady_state_model;
   Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
 end;
 
-set_time(1Q1);
-
 initval;
   LoggedProductivityInnovation = 0;
 end;
@@ -45,7 +43,7 @@ end;
 
 steady;
 
-perfect_foresight_setup(periods=200);
+perfect_foresight_setup(periods=200, first_simulation_period = 1Q2);
 perfect_foresight_solver;
 
 if ~oo_.deterministic_simulation.status
diff --git a/tests/deterministic_simulations/homotopy_histval.mod b/tests/deterministic_simulations/homotopy_histval.mod
index 551bab78ffa04d70214cea615317a2b50bf5b5a1..554672b661853a45ba2504a2ace14b163c058e72 100644
--- a/tests/deterministic_simulations/homotopy_histval.mod
+++ b/tests/deterministic_simulations/homotopy_histval.mod
@@ -31,8 +31,6 @@ steady_state_model;
   Consumption = exp(LoggedProductivity)*Capital^alpha-delta*Capital;
 end;
 
-set_time(1Q1);
-
 initval;
   LoggedProductivityInnovation = 0;
   LoggedProductivity = 10;
@@ -46,7 +44,7 @@ histval;
  LoggedProductivity(0)=10;
 end;
 
-perfect_foresight_setup(periods=200);
+perfect_foresight_setup(periods=200, first_simulation_period = 1Q2);
 perfect_foresight_solver;
 
 if ~oo_.deterministic_simulation.status