matlab/+hank/private/check_missingredundant.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/>.
+%
+% Check for missing and redundant fields in a structure.
+%
+% INPUTS
+% - f [structure]: structure to check
+% - f_name [char]: name to display in error/warning messages
+% - symbs [cell array of char]: list of expected field names
+% - nowarningredundant [logical]: if true, suppresses warnings about redundant fields
+%
+% OUTPUTS
+% - (none) This function throws an error if any expected field is missing,
+%   and issues a warning if redundant fields are found (unless nowarningredundant is true).
+%
+% DESCRIPTION
+% Verifies that all fields listed in `symbs` are present in the structure `f`.
+% Throws an error if any expected field is missing.
+% If `nowarningredundant` is false, checks if `f` contains fields not listed
+% in `symbs`, and issues a warning if redundant fields are found.
+function [] = check_missingredundant(f, f_name, symbs, nowarningredundant)
+   fields = fieldnames(f);
+   symbs_in_fields = ismember(symbs, fields);
+   if ~all(symbs_in_fields)
+      error('Misspecified steady-state input `ss`. The following fields are missing in `%s`: %s', f_name, strjoin(symbs(~symbs_in_fields)));
+   end
+   if ~nowarningredundant
+      fields_in_symbs = ismember(fields, symbs);
+      if ~all(fields_in_symbs)
+         warning('Steady-state input `ss`. The following fields are redundant in `%s`: %s', f_name, strjoin(fields(~fields_in_symbs)));
+      end
+   end
+end
\ No newline at end of file

matlab/+hank/private/check_permutation.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/>.
+%
+% Check and returns a permutation order for variables.
+%
+% INPUTS
+% - s [structure]: structure containing the permutation information
+% - f_name [char]: name of the field within `s` that should contain the permutation list
+% - s_name [char]: full name of the structure `s` for error messages
+% - symbs [cell array of char]: list of expected variable names
+%
+% OUTPUTS
+% - out_order [cell array or string array]: permutation order of variables
+%
+% DESCRIPTION
+% Checks that the field `f_name` exists in structure `s` and contains a valid
+% permutation (i.e., a cell array of character vectors or a string array).
+% If the field is missing, returns the default order given by `symbs`.
+% Throws an error if the specified variables are not consistent with `symbs`.
+function [out_order] = check_permutation(s, f_name, s_name, symbs)
+   if ~isfield(s, f_name)
+      out_order = symbs;
+   else
+      f = s.(f_name);
+      if ~isstring(f) && ~iscellstr(f)
+         error('Misspecified steady-state input `ss`: the `%s.%s` field should be a cell array of character vectors or a string array.', s_name, f_name);
+      end
+      err_var = setdiff(f, symbs);
+      if ~isempty(err_var)
+         error('Misspecified steady-state input `ss`: the set of variables of the `%s.%s` field is not consistent with the information in M_ and the other fields in the steady-state structure `ss`. Problematic variables: %s', s_name, f_name, strjoin(err_var));
+      end
+      out_order = f;
+   end
+end
\ No newline at end of file

matlab/+identification/analysis.m

@@ -297,6 +297,15 @@ if info(1) == 0 %no errors in solution
                 options_.analytic_derivation     = -2; %this sets asy_Hess=1 in dsge_likelihood.m                
                 [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, options_.varobs);
                 dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
+                % set info on missing data
+                if dataset_info.missing.state
+                    [dataset_info.missing.aindex, dataset_info.missing.number_of_observations, dataset_info.missing.no_more_missing_observations, dataset_info.missing.vindex] = ...
+                        describe_missing_data(dataset_.data);
+                else
+                    dataset_info.missing.aindex = num2cell(transpose(repmat(1:dataset_.vobs,dataset_.nobs,1)),1);
+                    dataset_info.missing.no_more_missing_observations = 1;
+                end
+
                 derivatives_info.no_DLIK = 1;
                 bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding 
                 %note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.

matlab/+identification/run.m

@@ -240,7 +240,7 @@ end
 
 % check for external draws, i.e. set pdraws0 for a gsa analysis
 if options_ident.gsa_sample_file
-    GSAFolder = checkpath('gsa',dname);
+    GSAFolder = CheckPath('gsa',dname);
     if options_ident.gsa_sample_file==1
         load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable');
     elseif options_ident.gsa_sample_file==2
@@ -300,6 +300,11 @@ options_.analytic_derivation=options_ident.analytic_derivation;
     % 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
 options_ = set_default_option(options_,'datafile','');
 options_.mode_compute = 0;
+if strcmp('slice',options_.posterior_sampler_options.posterior_sampling_method)
+    if strfind(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1')
+        options_.posterior_sampler_options.sampling_opt=strrep(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1','slice_initialize_with_mode'',0') % reset option to prevent crash in check_posterior_sampler_options.m if mode_compute is set to 0, see https://git.dynare.org/Dynare/dynare/-/issues/1957
+    end
+end
 options_.plot_priors = 0;
 options_.smoother = 1;
 options_.options_ident = [];
@@ -485,7 +490,7 @@ if iload <=0
             kk=0;
             while kk<50 && info(1)
                 kk=kk+1;
-                params = Prior.draw();
+                params = Prior.draw()'; %column vector is expected
                 options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
                 % perform identification analysis
                 [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ...
@@ -504,7 +509,10 @@ if iload <=0
         else
             % found a (random) point that solves the model
             fprintf('Found a random draw from the priors that solves the model:\n');
-            disp(params);
+            labels = name;
+            headers = {'Name', 'Value'};
+            lh = cellofchararraymaxlength(labels)+2;
+            dyntable(options_, 'Feasible draw', headers, labels, params', lh, 10, 6);
             fprintf('Identification now continues for this draw.');
             parameters = 'Random_prior_params';
             parameters_TeX = 'Random prior parameter draw';

matlab/+occbin/DSGE_smoother.m

@@ -85,7 +85,7 @@ if  options_.occbin.smoother.linear_smoother && nargin==12
     oo_.occbin.linear_smoother.alphahat0=alphahat0;
     oo_.occbin.linear_smoother.state_uncertainty0=state_uncertainty0;
 
-    fprintf('\nOccbin: linear smoother done.\n')
+    disp_verbose('Occbin: linear smoother done.',options_.verbosity)
     options_.occbin.smoother.status=true;
 end
 % if init_mode
@@ -123,22 +123,20 @@ occbin_options.opts_simul = opts_simul; % this builds the opts_simul options fie
 occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_binding_indicator;
 occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;
 
-error_indicator=false;
 options_.noprint = true;
 
-try
-    %blanket try-catch should be replaced be proper error handling, see https://git.dynare.org/Dynare/dynare/-/merge_requests/2226#note_20318
-    [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);%     T1=TT;
-catch ME
-    error_indicator=true;
-    disp(ME.message)
-    for iter = 1:numel(ME.stack)
-        ME.stack(iter)
-    end
+[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,~,error_indicator] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);%     T1=TT;
+
+if error_indicator(1) || isempty(alphahat0)
+    if ~options_.occbin.smoother.linear_smoother || nargin~=12 %make sure linear smoother results are set before using them
+        options_.occbin.smoother.status=false;
+        [~,etahat,~,~,~,~,~,~,~,~,~,~,~,~,~,~,alphahat0] = ...
+            DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
+         options_.occbin.smoother.status=true;
+    else
+        etahat= oo_.occbin.linear_smoother.etahat;
+        alphahat0= oo_.occbin.linear_smoother.alphahat0;
     end
-if error_indicator || isempty(alphahat0)
-    etahat= oo_.occbin.linear_smoother.etahat;
-    alphahat0= oo_.occbin.linear_smoother.alphahat0;
     base_regime = struct();
     if M_.occbin.constraint_nbr==1
         base_regime.regime = 0;
@@ -173,7 +171,7 @@ occbin_options.first_period_occbin_update = inf;
 opts_regime.binding_indicator=[];
 regime_history0 = regime_history;
 
-fprintf('Occbin smoother iteration 1.\n')
+disp_verbose('Occbin smoother iteration 1.',options_.verbosity)
 opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
 opts_simul.exo_pos = 1:M_.exo_nbr;
 opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
@@ -182,8 +180,7 @@ opts_simul.periods = size(opts_simul.SHOCKS,1);
 options_.occbin.simul=opts_simul;
 [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
 if out.error_flag
-    fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
-    print_info(out.error_flag, options_.noprint, options_)
+    disp_verbose('Occbin smoother:: simulation within smoother did not converge.',options_.verbosity)    
     oo_.occbin.smoother.error_flag=321;
     return;
 end
@@ -226,7 +223,7 @@ end
 
 while is_changed && maxiter>iter && ~is_periodic
     iter=iter+1;
-    fprintf('Occbin smoother iteration %u.\n', iter)
+    disp_verbose(sprintf('Occbin smoother iteration %u.', iter),options_.verbosity)
     occbin_options.opts_regime.regime_history=regime_history;
     [alphahat,etahat,epsilonhat,~,SteadyState,trend_coeff,~,T0,R0,P,~,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0]...
         = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);
@@ -245,8 +242,7 @@ while is_changed && maxiter>iter && ~is_periodic
     options_.occbin.simul=opts_simul;
     [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
     if out.error_flag
-        fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
-        print_info(out.error_flag, false, options_)
+        disp_verbose('Occbin smoother:: simulation within smoother did not converge.',options_.verbosity)
         oo_.occbin.smoother.error_flag=321;
         return;
     end
@@ -300,13 +296,13 @@ while is_changed && maxiter>iter && ~is_periodic
             eee(:,k) = eig(TT(:,:,k));
         end
         if options_.debug
-        err_eig(iter-1) = max(max(abs(sort(eee)-sort(sto_eee))));
-        err_alphahat(iter-1) = max(max(max(abs(alphahat-sto_alphahat))));
-        err_etahat(iter-1) = max(max(max(abs(etahat-sto_etahat{iter-1}))));
-        err_CC(iter-1) = max(max(max(abs(CC-sto_CC))));
-        err_RR(iter-1) = max(max(max(abs(RR-sto_RR))));
-        err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
-    end
+            err_eig(iter-1) = max(max(abs(sort(eee)-sort(sto_eee))));
+            err_alphahat(iter-1) = max(max(max(abs(alphahat-sto_alphahat))));
+            err_etahat(iter-1) = max(max(max(abs(etahat-sto_etahat{iter-1}))));
+            err_CC(iter-1) = max(max(max(abs(CC-sto_CC))));
+            err_RR(iter-1) = max(max(max(abs(RR-sto_RR))));
+            err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
+        end
     end
 
     if occbin_smoother_debug || is_periodic
@@ -391,22 +387,22 @@ if occbin_smoother_debug
 end
 
 if (maxiter==iter && is_changed) || is_periodic
-    disp('occbin.DSGE_smoother: smoother did not converge.')
-    fprintf('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.\n')
+    disp_verbose('occbin.DSGE_smoother: smoother did not converge.',options_.verbosity)
+    disp_verbose('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.',options_.verbosity)
     if is_periodic
         oo_.occbin.smoother.error_flag=0;
-        fprintf('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.\n')
-        fprintf('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.\n')
+        disp_verbose('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.',options_.verbosity)
+        disp_verbose('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.',options_.verbosity)
     else
-        fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
+        disp_verbose('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.',options_.verbosity)
         oo_.occbin.smoother=[];
         oo_.occbin.smoother.error_flag=322;
     end
 else
-    disp('occbin.DSGE_smoother: smoother converged.')
+    disp_verbose('occbin.DSGE_smoother: smoother converged.',options_.verbosity)
     oo_.occbin.smoother.error_flag=0;
     if occbin_smoother_fast && is_changed_start
-        disp('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge')
+        disp_verbose('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge',options_.verbosity)
     end
 end
 if (~is_changed || occbin_smoother_debug) && nargin==12
@@ -484,9 +480,9 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
                         fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
                         fprintf(fidTeX,'\\end{figure}\n');
                         fprintf(fidTeX,' \n');
+                    end
                 end
             end
-            end
         end
 
         if mod(j1,9)~=0 && j==M_.exo_nbr

matlab/+occbin/kalman_update_engine.m

@@ -109,6 +109,7 @@ if info==0
         newstart = regx(1).regimestart(end);
         diffstart = newstart-oldstart;
         regname = 'regimestart';
+        reg_string = 'regime';
     else
         newstart1 = regx(1).regimestart1(end);
         newstart2 = regx(1).regimestart2(end);
@@ -118,11 +119,13 @@ if info==0
         switch diffregime
             case 1
                 regname = 'regimestart1';
+                reg_string = 'regime1';
             case 2
                 regname = 'regimestart2';
+                reg_string = 'regime2';
+        end
     end
 end
-end
 if options_.occbin.filter.use_relaxation && diffstart>options_.occbin.filter.use_relaxation
     guess_regime = [base_regime base_regime];
     options_.occbin.filter.guess_regime = true;
@@ -132,9 +135,15 @@ if options_.occbin.filter.use_relaxation && diffstart>options_.occbin.filter.use
         % we reduce length until the converged regime does not change
 
         guess_regime(1).(regname)(end) = regx2(1).(regname)(end)-1;
-        if guess_regime(1).(regname)(end-1)==guess_regime(1).(regname)(end)
-            guess_regime(1).(regname)(end-1) = guess_regime(1).(regname)(end-1)-1;
-        end
+        if guess_regime(1).(regname)(end)==1
+            % make sure we enforce base regime
+            guess_regime(1).(regname)=guess_regime(1).(regname)(end);
+            guess_regime(1).(reg_string)=0;
+        else
+            if guess_regime(1).(regname)(end-1)==guess_regime(1).(regname)(end)
+                guess_regime(1).(regname)(end-1) = guess_regime(1).(regname)(end-1)-1;
+            end
+        end        
         if is_multivariate
             [ax2, a1x2, Px2, P1x2, vx2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, etahat2, alphahat2, V2] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
         else

matlab/default_option_values.m

@@ -12,7 +12,7 @@ function options_ = default_option_values(M_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2018-2024 Dynare Team
+% Copyright © 2018-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -195,8 +195,7 @@ ep.use_bytecode = 0;
 % Initialization of the perfect foresight equilibrium paths
 % * init=0, previous solution is used.
 % * init=1, a path generated with the first order reduced form is used.
-% * init=2, mix of cases 0 and 1.
-ep.init = 0;
+ep.use_first_order_solution_as_initial_guess = false;
 % Maximum number of iterations for the deterministic solver.
 ep.maxit = 500;
 % Number of periods for the perfect foresight model.
@@ -359,6 +358,11 @@ options_.k_order_solver = false; % by default do not use k_order_perturbation bu
 options_.partial_information = false;
 options_.conditional_variance_decomposition = [];
 
+% Heterogeneous agents
+options_.hank.tol_check_sum = 1e-6;
+options_.hank.nowarningredundant = true;
+options_.hank.nowarningdgrids = true;
+
 % Ramsey policy
 options_.ramsey_policy = false;
 options_.instruments = {};
@@ -453,7 +457,7 @@ options_.mh_posterior_mode_estimation = false;
 options_.smc_posterior_mode_estimation = false;
 options_.prefilter = 0;
 options_.presample = 0;
-options_.prior_trunc = 1e-10;
+options_.prior_trunc = 0;
 options_.smoother = false;
 options_.smoother_redux = false;
 options_.posterior_max_subsample_draws = 1200;

matlab/discretionary_policy/discretionary_policy.m

@@ -35,6 +35,9 @@ origorder = options_.order;
 options_.order = 1;
 [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list);
 
+if info(1)
+    return;
+end
 oo_.steady_state = oo_.dr.ys;
 
 if ~options_.noprint

matlab/discretionary_policy/discretionary_policy_1.m

@@ -42,6 +42,7 @@ if options_.steadystate_flag
     [ys,M_.params,info] = evaluate_steady_state_file(endo_steady_state,[exo_steady_state; exo_det_steady_state],M_, ...
                                                     options_,~options_.steadystate.nocheck);
     if info(1)
+        params=M_.params;
         return;
     end
 else
@@ -113,5 +114,7 @@ T=H(dr.order_var,dr.order_var);
 dr.ghu=G(dr.order_var,:);
 if M_.maximum_endo_lag
     Selection=M_.lead_lag_incidence(1,dr.order_var)>0;%select state variables
+else
+    Selection=[];
 end
 dr.ghx=T(:,Selection);

matlab/dsample.m

+function dsample(s1,s2)
+% function dsample(s1,s2)
+% This optional command permits to reduce the number of periods considered in following output commands.
+% If only one argument is provided, output is from period 1 to the period specified in the DSAMPLE command.
+% If two arguments are present output is done for the interval between the two periods.
+% DSAMPLE without arguments reset the sample to the one specified by PERIODS
+%
+% INPUTS
+%    s1:      first period
+%    s2:      last period
+%
+% OUTPUTS
+%    none
+%
+% SPECIAL REQUIREMENTS
+%    none
+% Note: this is a standalone command documented in the manual
+
+% Copyright © 2001-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/>.
+
+global options_
+
+options_.smpl = zeros(2,1) ;
+
+if nargin == 0
+    options_.smpl(1) = 1 ;
+    options_.smpl(2) = options_.periods ;
+elseif nargin == 1
+    if s1 > options_.periods
+        error('DSAMPLE: argument greater than number of periods');
+    end
+    options_.smpl(1) = 1 ;
+    options_.smpl(2) = s1 ;
+else
+    if s1 > options_.periods || s2 > options_.periods
+        error('DSAMPLE: one of the arguments is greater than number of periods');
+    end
+    options_.smpl(1) = s1 ;
+    options_.smpl(2) = s2 ;
+end
+
+% 02/23/01 MJ added error checking
\ No newline at end of file

dseries @ 657f8b8e

-Subproject commit 5436f91fe709b7e9f31302714cc1f9a7271001c6
+Subproject commit 657f8b8eaa4968c6f28434822d0f063f9c955935

matlab/ep/ep_problem_0.m

+function [res, A, info] = ep_problem_0(y, x, pfm)
+
+% Evaluate the residuals and stacked jacobian of a stochastic perfect
+% foresight, considering sequences of future innovations in a perfect n-ary tree.
+%
+% INPUTS:
+% - y      [double]   m×1 vector (endogenous variables in all periods and future worlds).
+% - x      [double]   q×1 vector of exogenous variables.
+% - pfm    [struct]   Definition of the perfect foresight model to be solved.
+%
+% OUTPUTS:
+% - res    [double]   m×1 vector, residuals of the stacked equations.
+% - A      [double]   m×m sparse matrix, jacobian of the stacked equations.
+% - info   [logical]  scalar
+%
+% REMARKS:
+% [1] The structure pfm holds the given initial condition for the states (pfm.y0) and the terminal condition
+
+info = false;
+
+params = pfm.params;
+steady_state = pfm.steady_state;
+ny = pfm.ny;
+periods = pfm.periods;
+dynamic_model = pfm.dynamic_model;
+lead_lag_incidence = pfm.lead_lag_incidence;
+i_cols_1 = pfm.i_cols_1;
+i_cols_j = pfm.i_cols_j;
+i_cols_T = pfm.i_cols_T;
+order = pfm.order;
+hybrid_order = pfm.hybrid_order;
+h_correction = pfm.h_correction;
+nodes = pfm.nodes;
+weights = pfm.weights;
+nnodes = pfm.nnodes;
+
+i_cols_p = pfm.i_cols_p;
+i_cols_s = pfm.i_cols_s;
+i_cols_f = pfm.i_cols_f;
+i_rows = pfm.i_rows;
+
+i_cols_Ap = pfm.i_cols_Ap;
+i_cols_As = pfm.i_cols_As;
+i_cols_Af = pfm.i_cols_Af;
+i_hc = pfm.i_hc;
+
+Y = pfm.Y;
+Y(pfm.i_upd_y) = y;
+
+A1 = pfm.A1;
+res = pfm.res;
+
+for i = 1:order+1
+    i_w_p = 1;
+    for j = 1:nnodes^(i-1)
+        innovation = x;
+        if i > 1
+            innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
+        end
+        if i <= order
+            for k=1:nnodes
+                if hybrid_order && i==order
+                    z = [Y(i_cols_p,i_w_p);
+                         Y(i_cols_s,j);
+                         Y(i_cols_f,(j-1)*nnodes+k)+h_correction(i_hc)];
+                else
+                    z = [Y(i_cols_p,i_w_p);
+                         Y(i_cols_s,j);
+                         Y(i_cols_f,(j-1)*nnodes+k)];
+                end
+
+                [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+                if i == 1
+                    % in first period we don't keep track of
+                    % predetermined variables
+                    i_cols_A = [i_cols_As - ny; i_cols_Af];
+                    A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
+                else
+                    i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+                    A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
+                end
+                res(:,i,j) = res(:,i,j)+weights(k)*d1;
+                i_cols_Af = i_cols_Af + ny;
+            end
+        else
+            z = [Y(i_cols_p,i_w_p);
+                 Y(i_cols_s,j);
+                 Y(i_cols_f,j)];
+            [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+            if i == 1
+                % in first period we don't keep track of
+                % predetermined variables
+                i_cols_A = [i_cols_As - ny; i_cols_Af];
+                A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
+            else
+                i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+                A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
+            end
+            res(:,i,j) = d1;
+            i_cols_Af = i_cols_Af + ny;
+        end
+        i_rows = i_rows + ny;
+        if mod(j,nnodes) == 0
+            i_w_p = i_w_p + 1;
+        end
+        if i > 1
+            if mod(j,nnodes) == 0
+                i_cols_Ap = i_cols_Ap + ny;
+            end
+            i_cols_As = i_cols_As + ny;
+        end
+    end
+    i_cols_p = i_cols_p + ny;
+    i_cols_s = i_cols_s + ny;
+    i_cols_f = i_cols_f + ny;
+end
+nzA = cell(periods,pfm.world_nbr);
+for j=1:pfm.world_nbr
+    i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
+    offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
+    offset_r = (j-1)*ny;
+    for i=order+2:periods
+        [d1,jacobian] = dynamic_model(Y(i_rows_y,j), x, params, steady_state, i+1);
+        if i == periods
+            [ir,ic,v] = find(jacobian(:,i_cols_T));
+        else
+            [ir,ic,v] = find(jacobian(:,i_cols_j));
+        end
+        nzA{i,j} = [offset_r+ir,offset_c+pfm.icA(ic), v]';
+        res(:,i,j) = d1;
+        i_rows_y = i_rows_y + ny;
+        offset_c = offset_c + pfm.world_nbr*ny;
+        offset_r = offset_r + pfm.world_nbr*ny;
+    end
+end
+A2 = [nzA{:}]';
+A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*pfm.world_nbr,pfm.dimension)];
+res = res(pfm.i_upd_r);

matlab/ep/ep_problem_2.m → matlab/ep/ep_problem_1.m

-function [res,A,info] = ep_problem_2(y,x,pm)
+function [res, A, info] = ep_problem_1(y, x, pfm)
 
-info = 0;
-res = [];
+% Evaluate the residuals and stacked jacobian of a stochastic perfect
+% foresight, considering sequences of future innovations in a sparse tree.
+%
+% INPUTS:
+% - y      [double]   m×1 vector (endogenous variables in all periods and future worlds).
+% - x      [double]   q×1 vector of exogenous variables.
+% - pfm    [struct]   Definition of the perfect foresight model to be solved.
+%
+% OUTPUTS:
+% - res    [double]   m×1 vector, residuals of the stacked equations.
+% - A      [double]   m×m sparse matrix, jacobian of the stacked equations.
+% - info   [logical]  scalar
+%
+% REMARKS:
+% [1] The structure pfm holds the given initial condition for the states (pfm.y0) and the terminal condition
+
+info = false;
 A = [];
 
-dynamic_model = pm.dynamic_model;
-ny = pm.ny;
-params = pm.params;
-steady_state = pm.steady_state;
-order = pm.order;
-nodes = pm.nodes;
-nnodes = pm.nnodes;
-weights = pm.weights;
-h_correction = pm.h_correction;
-dimension = pm.dimension;
-world_nbr = pm.world_nbr;
-nnzA = pm.nnzA;
-periods = pm.periods;
-i_rows = pm.i_rows';
-i_cols = pm.i_cols;
-nyp = pm.nyp;
-nyf = pm.nyf;
-hybrid_order = pm.hybrid_order;
-i_cols_1 = pm.i_cols_1;
-i_cols_j = pm.i_cols_j;
-icA = pm.icA;
-i_cols_T = pm.i_cols_T;
-eq_index = pm.eq_index;
+dynamic_model = pfm.dynamic_model;
+ny = pfm.ny;
+params = pfm.params;
+steady_state = pfm.steady_state;
+order = pfm.stochastic_order;
+nodes = pfm.nodes;
+nnodes = pfm.nnodes;
+weights = pfm.weights;
+h_correction = pfm.h_correction;
+dimension = pfm.dimension;
+world_nbr = pfm.world_nbr;
+nnzA = pfm.nnzA;
+periods = pfm.periods;
+i_rows = pfm.i_rows';
+i_cols = pfm.i_cols;
+nyp = pfm.nyp;
+nyf = pfm.nyf;
+hybrid_order = pfm.hybrid_order;
+i_cols_1 = pfm.i_cols_1;
+i_cols_j = pfm.i_cols_j;
+icA = pfm.icA;
+i_cols_T = pfm.i_cols_T;
+eq_index = pfm.eq_index;
+
+
 
 i_cols_p = i_cols(1:nyp);
 i_cols_s = i_cols(nyp+(1:ny));
 i_cols_f = i_cols(nyp+ny+(1:nyf));
-i_cols_A = i_cols;
 i_cols_Ap0 = i_cols_p;
 i_cols_As = i_cols_s;
 i_cols_Af0 = i_cols_f - ny;
@@ -40,9 +56,10 @@ i_hc = i_cols_f - 2*ny;
 nzA = cell(periods,world_nbr);
 res = zeros(ny,periods,world_nbr);
 Y = zeros(ny*(periods+2),world_nbr);
-Y(1:ny,1) = pm.y0;
+Y(1:ny,1) = pfm.y0;
 Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr);
-Y(pm.i_upd_y) = y;
+Y(pfm.i_upd_y) = y;
+
 offset_r0 = 0;
 for i = 1:order+1
     i_w_p = 1;
@@ -78,10 +95,10 @@ for i = 1:order+1
                 else
                     k1 = (nnodes-1)*(i-1)+k;
                 end
-                if hybrid_order == 2 && (k > 1 || i == order)
+                if hybrid_order && (k > 1 || i == order)
                     z = [Y(i_cols_p,1);
                          Y(i_cols_s,1);
-                         Y(i_cols_f,k1)+h_correction(i_hc)];
+                             Y(i_cols_f,k1)+h_correction(i_hc)];
                 else
                     z = [Y(i_cols_p,1);
                          Y(i_cols_s,1);
@@ -192,4 +209,4 @@ if nargout > 1
     iA = [nzA{:}]';
     A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension);
 end
-res = res(pm.i_upd_r);
\ No newline at end of file
+res = res(pfm.i_upd_r);

matlab/ep/euler_equation_error.m

 function e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
 % e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
 
-% Copyright © 2016-2023 Dynare Team
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -23,7 +23,7 @@ ep = options_.ep;
 y1 = extended_path_core(ep.periods, ...
                         M_.endo_nbr, M_.exo_nbr, ...
                         innovations.positive_var_indx, ...
-                        x, ep.init, y0, oo_.steady_state, ...
+                        x, ep.use_first_order_solution_as_initial_guess, y0, oo_.steady_state, ...
                         0, ...
                         ep.stochastic.order, M_, ...
                         pfm, ep.stochastic.algo, ...
@@ -38,7 +38,7 @@ for i=1:length(nodes)
     x2 = x1;
     x2(2,:) = x2(2,:) + nodes(i,:);
     y2 = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, ...
-                            innovations.positive_var_indx, x2, ep.init, ...
+                            innovations.positive_var_indx, x2, ep.use_first_order_solution_as_initial_guess, ...
                             y1, oo_.steady_state, 0, ...
                             ep.stochastic.order, M_, pfm, ep.stochastic.algo, ...
                             ep.solve_algo, ep.stack_solve_algo, options_.lmmcp, ...

matlab/ep/extended_path.m

@@ -19,7 +19,7 @@ function [ts,oo_] = extended_path(initialconditions, samplesize, exogenousvariab
 %
 % SPECIAL REQUIREMENTS
 
-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,18 +36,22 @@ 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/>.
 
-[initialconditions, innovations, pfm, ep, verbosity, options_, oo_] = ...
+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_);
 
-[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, samplesize,M_,options_,oo_);
+[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables, samplesize, M_, options_, oo_);
 
 % Initialize the matrix for the paths of the endogenous variables.
-endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
+endogenous_variables_paths = NaN(M_.endo_nbr, samplesize+1);
 endogenous_variables_paths(:,1) = initialconditions;
 
 % Set waitbar (graphic or text  mode)
 hh_fig = dyn_waitbar(0,'Please wait. Extended Path simulations...');
-set(hh_fig,'Name','EP simulations.');
+set(hh_fig,'Name', 'EP simulations.' );
 
 % Initialize while-loop index.
 t = 1;
@@ -67,14 +71,27 @@ while (t <= samplesize)
     else
         initialguess = [];
     end
-    [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
-                                                      spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
-                                                      oo_.steady_state, ...
-                                                      verbosity, ep.stochastic.order, ...
-                                                      M_, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
-                                                      options_.lmmcp, ...
-                                                      options_, ...
-                                                      oo_, initialguess);
+    if t>2
+        [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, ...
+                                                                                                           y);
+    else
+        [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, ...
+                                                                                                                          []);
+    end
     if ~info_convergence
         msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s)!', int2str(t));
         warning(msg)
@@ -106,4 +123,4 @@ if any(isnan(endogenous_variables_paths(:)))
 end
 ts = dseries(transpose(endogenous_variables_paths), initial_period, M_.endo_names);
 
-oo_.endo_simul = transpose(ts.data);
\ No newline at end of file
+oo_.endo_simul = transpose(ts.data);

matlab/ep/extended_path_core.m

-function [y, info_convergence, endogenousvariablespaths] = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ...
-                                                  exo_simul,init,initial_conditions,...
-                                                  steady_state, ...
-                                                  debug,order,M_,pfm,algo,solve_algo,stack_solve_algo,...
-                                                  olmmcp,options_,oo_,initialguess)
+function [y1, info_convergence, endogenousvariablespaths, y, pfm, options_] = extended_path_core(positive_var_indx, ...
+                                                                                             exo_simul, ...
+                                                                                             initial_conditions ,...
+                                                                                             pfm, ...
+                                                                                             M_, ...
+                                                                                             options_, ...
+                                                                                             oo_, ...
+                                                                                             initialguess, ...
+                                                                                             y)
 
-% Copyright © 2016-2023 Dynare Team
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -23,10 +27,22 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
 
 ep = options_.ep;
 
+periods = ep.periods;
+endo_nbr = M_.endo_nbr;
+exo_nbr = M_.exo_nbr;
+init = options_.ep.use_first_order_solution_as_initial_guess;
+steady_state = oo_.steady_state;
+debug = options_.verbosity;
+order = options_.ep.stochastic.order;
+
+algo = ep.stochastic.algo;
+solve_algo = ep.solve_algo;
+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==19 && ~isempty(initialguess)
+    if nargin>7 && ~isempty(initialguess)
         % Note that the first column of initialguess should be equal to initial_conditions.
         endo_simul = initialguess;
     else
@@ -34,6 +50,10 @@ else
     end
 end
 
+if nargin~=9
+    y = [];
+end
+
 oo_.endo_simul = endo_simul;
 
 if debug
@@ -41,30 +61,39 @@ if debug
 end
 
 if options_.bytecode && order > 0
-    error('Option order > 0 of extended_path command is not compatible with bytecode option.')
+    error('Option order > 0 of extended_path command with order>0 is not compatible with bytecode option.')
 end
 if options_.block && order > 0
-    error('Option order > 0 of extended_path command is not compatible with block option.')
+    error('Option order > 0 of extended_path command with order>0 is not compatible with block option.')
 end
 
 if order == 0
+    % Extended Path
     options_.periods = periods;
     options_.block = pfm.block;
     oo_.endo_simul = endo_simul;
     oo_.exo_simul = exo_simul;
     oo_.steady_state = steady_state;
-    options_.lmmcp = olmmcp;
     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)
-        case 0
-            [flag, endogenousvariablespaths] = ...
-            solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
-        case 1
-            [flag, endogenousvariablespaths] = ...
-            solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
+      case 0
+        % Full tree of future trajectories.
+        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, 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>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, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm);
+        end
     end
     info_convergence = ~flag;
 end
@@ -74,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

matlab/ep/extended_path_homotopy.m

 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,12 +35,12 @@ 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
                 [flag, endo_simul_new] = ...
-                    solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
+                    solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
               case 1
                 [flag, endo_simul_new] = ...
                     solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
@@ -102,12 +102,12 @@ 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
                 [flag, endo_simul_new] = ...
-                    solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
+                    solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
               case 1
                 [flag, endo_simul_new] = ...
                     solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);

matlab/ep/extended_path_initialization.m

-function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
+function [initial_conditions, innovations, pfm, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
 % [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
 % Initialization of the extended path routines.
 %
@@ -16,7 +16,7 @@ function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] =
 %
 % SPECIAL REQUIREMENTS
 
-% Copyright © 2016-2024 Dynare Team
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,8 +36,7 @@ function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] =
 ep  = options_.ep;
 
 % Set verbosity levels.
-options_.verbosity = ep.verbosity;
-verbosity = ep.verbosity+ep.debug;
+options_.verbosity = ep.verbosity+ep.debug;
 
 % Set maximum number of iterations for the deterministic solver.
 options_.simul.maxit = ep.maxit;
@@ -71,7 +70,7 @@ options_.stack_solve_algo = ep.stack_solve_algo;
 
 % Compute the first order reduced form if needed.
 dr = struct();
-if ep.init
+if ep.use_first_order_solution_as_initial_guess
     options_.order = 1;
     oo_.dr=set_state_space(dr,M_);
     [oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
@@ -101,7 +100,11 @@ end
 
 % hybrid correction
 pfm.hybrid_order = ep.stochastic.hybrid_order;
-if pfm.hybrid_order
+if pfm.hybrid_order==1
+    warning('extended_path:: hybrid=1 is equivalent to hybrid=0 (option value must be an integer greater than 1 to be effective).')
+end
+
+if pfm.hybrid_order>1
     oo_.dr = set_state_space(oo_.dr, M_);
     options = options_;
     options.order = pfm.hybrid_order;
@@ -110,29 +113,37 @@ else
     pfm.dr = [];
 end
 
+% Deactivate homotopy with SEP
+if ep.stochastic.order>0
+    options_.no_homotopy = true;
+end
+
 % number of nonzero derivatives
 pfm.nnzA = M_.NNZDerivatives(1);
 
 % setting up integration nodes if order > 0
-if ep.stochastic.order > 0
+if ep.stochastic.order>0
     [nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
     pfm.nodes = nodes;
     pfm.weights = weights;
     pfm.nnodes = nnodes;
     % compute number of blocks
-    [block_nbr,pfm.world_nbr] = get_block_world_nbr(ep.stochastic.algo,nnodes,ep.stochastic.order,ep.periods);
+    [pfm.block_nbr, pfm.world_nbr] = get_block_world_nbr(ep.stochastic.algo,nnodes,ep.stochastic.order,ep.periods);
 else
     block_nbr = ep.periods;
 end
 
 % set boundaries if mcp
-[lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
-pfm.eq_index = M_.dynamic_mcp_equations_reordering;
-if options_.ep.solve_algo == 10
-    options_.lmmcp.lb = repmat(lb,block_nbr,1);
-    options_.lmmcp.ub = repmat(ub,block_nbr,1);
-elseif options_.ep.solve_algo == 11
-    options_.mcppath.lb = repmat(lb,block_nbr,1);
-    options_.mcppath.ub = repmat(ub,block_nbr,1);
+if ~ep.stochastic.order
+    [lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
+    pfm.eq_index = M_.dynamic_mcp_equations_reordering;
+    if options_.ep.solve_algo == 10
+        options_.lmmcp.lb = repmat(lb, ep.periods, 1);
+        options_.lmmcp.ub = repmat(ub, ep.periods, 1);
+    elseif options_.ep.solve_algo == 11
+        options_.mcppath.lb = repmat(lb, ep.periods, 1);
+        options_.mcppath.ub = repmat(ub, ep.periods, 1);
+    end
+else
+    % For SEP, boundaries are set in solve_stochastic_perfect_foresight_model_{0,1}.m
 end
-pfm.block_nbr = block_nbr;

matlab/ep/extended_path_mc.m

 function Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
-% Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
+
 % Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
 % series of size T  is obtained by solving T perfect foresight models.
 %
@@ -19,7 +19,7 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
 %
 % SPECIAL REQUIREMENTS
 
-% Copyright © 2016-2023 Dynare Team
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +36,7 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-[initialconditions, innovations, pfm, ep, ~, options_, oo_] = ...
+[initialconditions, innovations, pfm, options_, oo_] = ...
     extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
 
 % Check the dimension of the first input argument
@@ -64,12 +64,12 @@ data = NaN(size(initialconditions, 1), samplesize+1, replic);
 vexo = NaN(innovations.effective_number_of_shocks, samplesize+1, replic);
 info = NaN(replic, 1);
 
-if ep.parallel
+if options_.ep.parallel
     % Use the Parallel toolbox.
     parfor i=1:replic
         innovations_ = innovations;
         oo__ = oo_;
-        [shocks, spfm_exo_simul, innovations_, oo__] = extended_path_shocks(innovations_, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo__);
+        [shocks, spfm_exo_simul, innovations_, oo__] = extended_path_shocks(innovations_, exogenousvariables(:,:,i), samplesize, M_, options_, oo__);
         endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
         endogenous_variables_paths(:,1) = initialconditions(:,1);
         exogenous_variables_paths = NaN(innovations_.effective_number_of_shocks,samplesize+1);
@@ -80,12 +80,21 @@ if ep.parallel
             t = t+1;
             spfm_exo_simul(2,:) = shocks(t-1,:);
             exogenous_variables_paths(:,t) = shocks(t-1,:);
-            [endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations_.positive_var_indx, ...
-                                                              spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
-                                                              oo__.steady_state, ...
-                                                              ep.verbosity, ep.stochastic.order, ...
-                                                              M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
-                                                              options_.lmmcp, options_, oo__);
+            if t>2
+                % Set initial guess for the solver (using the solution of the
+                % previous period problem).
+                initialguess = [endogenousvariablespaths(:, 2:end), oo_.steady_state];
+            else
+                initialguess = [];
+            end
+            [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(innovations.positive_var_indx, ...
+                                                                                                               spfm_exo_simul, ...
+                                                                                                               endogenous_variables_paths(:,t-1), ...
+                                                                                                               pfm, ...
+                                                                                                               M_, ...
+                                                                                                               options_, ...
+                                                                                                               oo__, ...
+                                                                                                               initialguess);
             if ~info_convergence
                 msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
                 warning(msg)
@@ -99,7 +108,7 @@ if ep.parallel
 else
     % Sequential approach.
     for i=1:replic
-        [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo_);
+        [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables(:,:,i), samplesize, M_, options_, oo_);
         endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
         endogenous_variables_paths(:,1) = initialconditions(:,1);
         exogenous_variables_paths = NaN(innovations.effective_number_of_shocks,samplesize+1);
@@ -109,12 +118,21 @@ else
             t = t+1;
             spfm_exo_simul(2,:) = shocks(t-1,:);
             exogenous_variables_paths(:,t) = shocks(t-1,:);
-            [endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
-                                                              spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
-                                                              oo_.steady_state, ...
-                                                              ep.verbosity, ep.stochastic.order, ...
-                                                              M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
-                                                              options_.lmmcp, options_, oo_);
+            if t>2
+                % Set initial guess for the solver (using the solution of the
+                % previous period problem).
+                initialguess = [endogenousvariablespaths(:, 2:end), oo_.steady_state];
+            else
+                initialguess = [];
+            end
+            [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(innovations.positive_var_indx, ...
+                                                                                                               spfm_exo_simul, ...
+                                                                                                               endogenous_variables_paths(:,t-1), ...
+                                                                                                               pfm, ...
+                                                                                                               M_, ...
+                                                                                                               options_, ...
+                                                                                                               oo_, ...
+                                                                                                               initialguess);
             if ~info_convergence
                 msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
                 warning(msg)

matlab/ep/extended_path_shocks.m

-function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
-% [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
-% Copyright © 2016-2023 Dynare Team
+function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables, sample_size, M_, options_, oo_)
+
+% Copyright © 2016-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,14 +19,14 @@ function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innov
 
 % Simulate shocks.
 if isempty(exogenousvariables)
-    switch ep.innovation_distribution
+    switch options_.ep.innovation_distribution
       case 'gaussian'
         shocks = zeros(sample_size, M_.exo_nbr);
         shocks(:,innovations.positive_var_indx) = transpose(transpose(innovations.covariance_matrix_upper_cholesky)*randn(innovations.effective_number_of_shocks,sample_size));
       case 'calibrated'
         options = options_;
         options.periods = options.ep.periods;
-        oo_local = make_ex_(M_,options,oo_);
+        oo_local = make_ex_(M_, options, oo_);
         shocks = oo_local.exo_simul(2:end,:);
       otherwise
         error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
@@ -38,4 +38,4 @@ end
 
 % Copy the shocks in exo_simul
 oo_.exo_simul = shocks;
-spfm_exo_simul = repmat(oo_.exo_steady_state',ep.periods+2,1);
\ No newline at end of file
+spfm_exo_simul = repmat(oo_.exo_steady_state',options_.ep.periods+2,1);