.gitlab-ci.yml

@@ -107,7 +107,7 @@ pkg_macOS_x86_64:
   script:
     # Enforce the arm64 meson for rewrite, as a workaround to https://github.com/mesonbuild/meson/issues/12282
     - env PATH="/opt/homebrew/bin:$PATH" meson rewrite kwargs set project / version "$VERSION"
-    - ln -s ~/tarballs macOS/deps/x86_64
+    - ln -s ~/tarballs macOS/deps/
     - make -C macOS build-x86_64
   cache:
     key: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG"
@@ -127,7 +127,7 @@ pkg_macOS_arm64:
   script:
     # Enforce the arm64 meson for rewrite, as a workaround to https://github.com/mesonbuild/meson/issues/12282
     - env PATH="/opt/homebrew/bin:$PATH" meson rewrite kwargs set project / version "$VERSION"
-    - ln -s ~/tarballs macOS/deps/arm64
+    - ln -s ~/tarballs macOS/deps/
     - make -C macOS build-arm64
   cache:
     key: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG"

NEWS.md

+Announcement for Dynare 6.3 (on 2025-02-19)
+===========================================
+
+We are pleased to announce the release of Dynare 6.3.
+
+This maintenance release fixes various bugs.
+
+The Windows, macOS, MATLAB online and source packages are available for
+download at [the Dynare website](https://www.dynare.org/download/).
+
+This release is compatible with MATLAB versions ranging from 9.5 (R2018b) to
+24.2 (R2024b), and with GNU Octave versions ranging from 7.1.0 to 9.4.0 (NB:
+the Windows package requires version 9.4.0 specifically).
+
+Here is a list of the problems identified in version 6.2 and that have been
+fixed in version 6.3:
+
+* OccBin with option `smoother_inversion_filter` would crash the MCMC
+  estimation if the `filtered_variables` or `filter_step_ahead` options were
+  used
+* OccBin with option `smoother_inversion_filter` would use the PKF if
+  `mh_replic>0`
+* OccBin with option `smoothed_state_uncertainty` would crash Dynare if the
+  MCMC smoother was run after the classical one
+* OccBin's smoother would crash when encountering an internal error due to the
+  output of the linear smoother not having been computed
+* Calling `model_info` with `differentiate_forward_vars` would crash
+* The `identification` command would compute the asymptotic Hessian via
+  simulation instead of via moments as intended
+* The `identification` command would crash during prior sampling if the initial
+  draw did not solve the model
+* The `gsa_sample_file` was broken
+* Optimization algorithm `mode_compute=5` (`newrat`) would crash with
+  `analytic_derivation`
+* The `discretionary_policy` command would crash if the model was purely
+  forward-looking
+* The `dsample` command would crash
+* The `conditional_forecast_paths` block did not accept vector inputs
+* For MCMC chains with fewer than 6000 draws, the default number of `sub_draws`
+  used to compute posterior moments was incorrect
+* Bi-annual dates (e.g. `2024S1` or `2024H1`) were not accepted within Dynare
+  statements
+* Plotting `dseries` did not correctly show dates on the x-axis
+
+
 Announcement for Dynare 6.2 (on 2024-09-25)
 ===========================================
 

README.md

@@ -362,7 +362,7 @@ cd dynare
 ```
 - Configure Dynare from the source directory:
 ```sh
-meson setup -Dmatlab_path=<…> --buildtype=debugoptimized --prefer-static -Dfortran_args="['-B','/usr/local/lib']" build-matlab
+meson setup -Dmatlab_path=<…> --buildtype=debugoptimized --prefer-static -Dfortran_args="['-B','C:/msys64/usr/local/lib']" build-matlab
 ```
 where the path of MATLAB is specified. Note that you should use
 the MSYS2 notation and not put spaces in the MATLAB path, so you probably want

doc/manual/source/installation-and-configuration.rst

@@ -15,7 +15,7 @@ compilation steps are necessary in that case.
 In order to run Dynare, you need one of the following:
 
 * MATLAB, any version ranging from 9.5 (R2018b) to 24.2 (R2024b);
-* GNU Octave, any version ranging from 7.1.0 to 9.3.0, with the ``statistics`` package
+* GNU Octave, any version ranging from 7.1.0 to 9.4.0, with the ``statistics`` package
   from `Octave-Forge`_. Note however that the Dynare installer for Windows
   requires a more specific version of Octave, as indicated on the download
   page.
@@ -100,6 +100,18 @@ Debian, Ubuntu and Linux Mint).
 On macOS
 --------
 
+.. warning::
+
+   Installing into ``/Applications/dynare`` might fail if you have older versions of Dynare already installed in ``/Applications/Dynare``.
+   To fix this, modify the ownership by executing the following command in Terminal.app::
+
+     sudo chown -R "$USER":staff /Applications/Dynare
+
+   Alternatively, you can modify the installation path in the automated installed using *Customize* and *Location*.
+   After installation, the folder will contain several sub-directories, among which are ``matlab``, ``mex``, and ``doc``.
+   Several versions of Dynare can coexist (by default in ``/Applications/Dynare``),
+   as long as you correctly adjust your path settings (see :ref:`words-warning`).
+
 With MATLAB
 ^^^^^^^^^^^
 
@@ -109,15 +121,6 @@ and follow the instructions.
 This installation does not require administrative privileges.
 If for some reason admin rights are requested, use *Change Install Location* and select *Install for me only*.
 The default installation directory is ``/Applications/Dynare/x.y-arch``.
-Installing into ``/Applications/dynare`` might fail if you have older versions of Dynare already installed in ``/Applications/Dynare``.
-To fix this, modify the ownership by executing the following command in Terminal.app::
-
-  sudo chown -R "$USER":staff /Applications/Dynare
-
-Alternatively, you can modify the installation path in the automated installed using *Customize* and *Location*.
-After installation, the folder will contain several sub-directories, among which are ``matlab``, ``mex``, and ``doc``.
-Several versions of Dynare can coexist (by default in ``/Applications/Dynare``),
-as long as you correctly adjust your path settings (see :ref:`words-warning`).
 
 It is recommended to install the Xcode Command Line Tools (this is an Apple product)
 and GCC via Homebrew_ (see :ref:`prerequisites-macos`).
@@ -142,6 +145,8 @@ once)::
 
   octave:1> pkg install -forge io statistics control struct optim
 
+If you want to use the `x13` functionality of `dseries`, you also need to build the `x13as` binary. [#fx13]_
+
 On FreeBSD
 ----------
 
@@ -334,3 +339,8 @@ Dynare unusable.
 .. _Dynare wiki: https://git.dynare.org/Dynare/dynare/wikis
 .. _Octave-Forge: https://octave.sourceforge.io/
 .. _Homebrew: https://brew.sh
+
+
+.. rubric:: Footnotes
+
+.. [#fx13] See the instructions at `<https://forum.dynare.org/t/missing-installation-package/27350/4>`__.
\ No newline at end of file

doc/manual/source/the-model-file.rst

@@ -3037,14 +3037,16 @@ Finding the steady state with Dynare nonlinear solver
                 ``perfect_foresight_solver`` command with purely backward,
                 forward or static models, or with routines for semi-structural
                 models, and it must *not* be combined with option ``block`` of
-                the ``model`` block. Also note that for those models, the block
-                decomposition is performed as if ``mfs=3`` had been passed to
-                the ``model`` block, and the decomposition is slightly
-                different because it is computed in a time-recursive fashion
-                (*i.e.* in such a way that the simulation is meant to be done
-                with the outer loop on periods and the inner loop on blocks;
-                while for models with both leads and lags, the outer loop is on
-                blocks and the inner loop is on periods).
+                the :bck:`model` block or :comm:`model_options` command. Also
+                note that for those models, the block decomposition is
+                performed as if ``mfs=3`` had been passed to the :bck:`model`
+                block or :comm:`model_options` command, and the decomposition
+                is slightly different because it is computed in a
+                time-recursive fashion (*i.e.* in such a way that the
+                simulation is meant to be done with the outer loop on periods
+                and the inner loop on blocks; while for models with both leads
+                and lags, the outer loop is on blocks and the inner loop is on
+                periods).
            ``14``
@@ -3560,14 +3562,15 @@ Getting information about the model
     |br| Prints the equations and the Jacobian matrix of the dynamic
     model stored in the bytecode binary format file. Can only be used
-    in conjunction with the ``bytecode`` option of the ``model``
-    block.
+    in conjunction with the ``bytecode`` option of the :bck:`model`
+    block or :comm:`model_options` command.
 .. command:: print_bytecode_static_model ;
     |br| Prints the equations and the Jacobian matrix of the static model
     stored in the bytecode binary format file. Can only be used in
-    conjunction with the ``bytecode`` option of the ``model`` block.
+    conjunction with the ``bytecode`` option of the :bck:`model`
+    block or :comm:`model_options` command.
 .. _det-simul:
@@ -3825,9 +3828,9 @@ speed-up on large models.
        See :ref:`solve_algo <solvalg>`. Allows selecting the solver
        used with ``stack_solve_algo=7``. Also used for purely backward, forward
-       and static models (when neither the ``block`` nor the ``bytecode`` option
-       of the ``model`` block is specified); for those models, the values
-       ``12`` and ``14`` are especially relevant.
+       and static models (when neither the ``block`` nor the ``bytecode``
+       option of the :bck:`model` block or :comm:`model_options` command is specified);
+       for those models, the values ``12`` and ``14`` are especially relevant.
     .. option:: no_homotopy
@@ -4015,7 +4018,7 @@ speed-up on large models.
        shooting and relaxation approaches. Note that round off errors
        are more important with this mixed strategy (user should check
        the reported value of the maximum absolute error). Only
-       available with option ``stack_solve_algo==0``.
+       available with ``stack_solve_algo`` option equal to ``0``.
     .. option:: linear_approximation
@@ -4024,8 +4027,8 @@ speed-up on large models.
        model. The model must be stationary and a steady state 
        needs to be provided. Linearization is conducted about the 
        last defined steady state, which can derive from ``initval``,
-       ``endval`` or a subsequent ``steady``. Only available with option
-       ``stack_solve_algo==0`` or ``stack_solve_algo==7``.
+       ``endval`` or a subsequent ``steady``. Only available with
+       ``stack_solve_algo`` option equal to ``0`` or ``7``.
     .. option:: steady_solve_algo = INTEGER
@@ -4616,6 +4619,144 @@ and ``endval`` blocks which are given a special ``learnt_in`` option.
     the terminal condition for exogenous indexed ``k``, as anticipated from
     period ``s``, is stored in ``oo_.pfwee.terminal_info(k,s)``.
+
+Controlling the path of endogenous variables
+--------------------------------------------
+
+In the usual perfect foresight problem, the user controls the path of exogenous
+variables for the simulation periods and the initial and terminal
+conditions for endogenous variables, while Dynare solves for the path of
+endogenous variables for the simulation periods.
+
+However, Dynare offers the possibility of controlling the value of some
+endogenous variables for some simulation periods (in which case some exogenous
+variables must be left free and are thus solved for by Dynare, to avoid
+over-determination of the problem). This exercise is called “conditional
+forecasting” in some contexts (even though one may argue that this is not
+really forecasting, since perfect foresight by the agents is assumed; for the
+stochastic case, see the :comm:`conditional_forecast` command).
+
+The description of controlled endogenous variables is done using the
+``perfect_foresight_controlled_paths`` block. The information given therein is
+then processed by the :comm:`perfect_foresight_setup` (or
+:comm:`perfect_foresight_with_expectation_errors_setup`) command, so
+that the next :comm:`perfect_foresight_solver` (or
+:comm:`perfect_foresight_with_expectation_errors_solver`) command
+computes the simulation with controlled paths. In particular,
+:mvar:`oo_.exo_simul` will contain the computed value of exogenous variables
+that have been left free.
+
+.. block:: perfect_foresight_controlled_paths ;
+           perfect_foresight_controlled_paths(OPTIONS...);
+
+    |br| This block is used to tell the perfect foresight solver that the value
+    of some endogenous variables will be controlled (in other words, they will
+    be exogenized). It also gives the period(s) for which this control applies,
+    the value(s) imposed to the endogenous variable(s), and the exogenous
+    variable(s) that are left free at the same period(s) (in other words,
+    those exogenous are endogenized).
+
+    The block should contain one or more occurrences of the following
+    group of four lines::
+
+      exogenize ENDOGENOUS_NAME;
+      periods INTEGER[:INTEGER] | DATE[:DATE] [[,] INTEGER[:INTEGER] | DATE[:DATE]]...;
+      values DOUBLE | (EXPRESSION)  [[,] DOUBLE | (EXPRESSION) ]...;
+      endogenize EXOGENOUS_NAME;
+
+    Note that it is possible to have both
+    ``perfect_foresight_controlled_paths`` and regular :bck:`shocks` blocks in
+    the same ``.mod`` file (assuming of course that taken together they do not
+    impose inconsistent constraints).
+
+    The ``perfect_foresight_controlled_paths`` block requires that the
+    :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` option be equal to
+    either ``0``, ``1``, ``2``, ``3``, ``6`` or ``7``, and is incompatible with
+    the :opt:`block` and :opt:`bytecode` options of the :bck:`model` block and
+    :comm:`model_options` command.
+
+    *Options*
+
+    .. option:: learnt_in = INTEGER | DATE
+
+       Used in conjunction with
+       :comm:`perfect_foresight_with_expectation_errors_setup` and
+       :comm:`perfect_foresight_with_expectation_errors_solver` commands,
+       specifies the period or date at which this controlled paths block is
+       learnt by agents.
+
+    *Example (perfect foresight)*
+
+        ::
+
+            var c k;
+            varexo x z;
+
+            ...
+
+            shocks;
+              var x;
+              periods 1;
+              values 1.2;
+            end;
+
+            perfect_foresight_controlled_paths;
+              exogenize c;
+              periods 2 4:5;
+              values 1.6 1.7;
+              endogenize x;
+
+              exogenize k;
+              periods 7:9;
+              values 13;
+              endogenize z;
+            end;
+
+	    perfect_foresight_setup(periods = 100);
+	    perfect_foresight_solver;
+
+        In this example, the exogenous variable ``x`` is equal to 1.2 in
+        period 1, but in periods 2, 4 and 5 it will be endogenized so that
+        endogenous variable `c` is equal to 1.6 in period 2 and then 1.7 in
+        periods 4 and 5. Similarly, the exogenous variable ``z`` will be
+        endogenized in periods 7 to 9 so that the endogenous variable ``k`` is
+        equal to 13 over the same periods.
+
+    *Example (perfect foresight with expectation errors)*
+
+        ::
+
+            var c;
+            varexo x;
+
+            ...
+
+            perfect_foresight_controlled_paths;
+              exogenize c;
+              periods 2002Y 2003Y:2005Y;
+              values 1.6 1.7;
+              endogenize x;
+            end;
+
+            perfect_foresight_controlled_paths(learnt_in=2004Y);
+              exogenize c;
+              periods 2004Y:2005Y;
+              values 1.8;
+              endogenize x;
+            end;
+
+            perfect_foresight_with_expectation_errors_setup(periods = 30,
+                first_simulation_period = 2001Y);
+            perfect_foresight_with_expectation_errors_solver;
+
+        In this example, agents in year 2001 (at beginning of the simulation)
+        compute their plan under the assumption that the endogenous variable
+        ``c`` will be equal to 1.6 in year 2002 and 1.7 from years 2003 to
+        2005, and that exogenous variable ``x`` will behave so as to fulfill
+        that constraint. Then, when 2004 arrives, they recompute their plan
+        under the assumption that ``c`` will be equal to 1.8 in years 2004 and
+        2005 (and again that ``x`` will be endogenized accordingly).
+
 .. _stoch-sol:
 Stochastic solution and simulation
@@ -5342,7 +5483,7 @@ which is described below.
        agents believe that there will no more shocks after period
        :math:`t+S`. This is an experimental feature and can be quite
        slow. A non-zero value is not compatible with the ``bytecode``
-       option of the ``model`` block.
+       option of the :bck:`model` block or :comm:`model_options` command.
        Default: ``0``.
     .. option:: hybrid
@@ -7580,7 +7721,7 @@ observed variables.
        Triggers the computation of the posterior distribution of
        IRFs. The length of the IRFs are controlled by the ``irf``
        option. Results are stored in ``oo_.PosteriorIRF.dsge`` (see
-       below for a description of this variable).
+       below for a description of this variable). Not compatible with OccBin.
     .. option:: relative_irf
@@ -7928,7 +8069,7 @@ observed variables.
        stored in ``oo_.PosteriorTheoreticalMoments`` (see
        :mvar:`oo_.PosteriorTheoreticalMoments`). The number of lags in
        the autocorrelation function is controlled by the ``ar``
-       option.
+       option. Not compatible with OccBin.
     .. option:: contemporaneous_correlation
@@ -8013,7 +8154,7 @@ observed variables.
        the posterior mode. If a Metropolis-Hastings is computed, the
        distribution of forecasts is stored in variables
        ``oo_.PointForecast`` and ``oo_.MeanForecast``. See
-       :ref:`fore`, for a description of these variables.
+       :ref:`fore`, for a description of these variables. Not compatible with OccBin.
     .. option:: tex
@@ -12334,7 +12475,7 @@ Optimal policy under discretion
     You must ensure that your objective is quadratic. Regarding the model, it must
     either be linear or solved at first order with an analytical steady state provided.
     In the first case, you should set the ``linear`` option of the
-    ``model`` block.
+    :bck:`model` block or :comm:`model_options` command.
     It is possible to use the :comm:`estimation` command after the
     ``discretionary_policy`` command, in order to estimate the model with

license.txt

@@ -97,7 +97,7 @@ License: GPL-3+
 
 Files: matlab/optimization/cmaes.m
 Copyright: 2001-2012 Nikolaus Hansen
-           2012-2017 Dynare Team
+           2012-2023 Dynare Team
 License: GPL-3+
 
 Files: matlab/optimization/solvopt.m
@@ -173,7 +173,7 @@ Comment: The author gave authorization to redistribute
 Files: matlab/+gsa/Morris_Measure_Groups.m
        matlab/+gsa/Sampling_Function_2.m
 Copyright: 2005 European Commission
-           2012-2013 Dynare Team
+           2012-2023 Dynare Team
 License: GPL-3+
 Comment: Written by Jessica Cariboni and Francesca Campolongo
  Joint Research Centre, The European Commission,
@@ -220,14 +220,14 @@ License: GPL-3+
 
 Files: matlab/missing/stats/quantile.m
 Copyright: 2014-2016 Christopher Hummersone
-           2016-2017 Dynare Team
+           2016-2023 Dynare Team
 License: GPL-3+
 
 Files: matlab/missing/stats/corr.m
 Copyright: 1993-1996 Kurt Hornik
            1996-2015 John W. Eaton
            2013-2015 Julien Bect
-           2016-2017 Dynare Team
+           2016-2023 Dynare Team
 License: GPL-3+
 
 Files: matlab/lmmcp/catstruct.m
@@ -260,7 +260,7 @@ Copyright: 2000-2022 Frank Schorfheide
 License: CC-BY-SA-4.0
 
 Files: doc/*.rst doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
-Copyright: 1996-2022 Dynare Team
+Copyright: 1996-2025 Dynare Team
 License: GFDL-NIV-1.3+
 
 Files: doc/dr.tex doc/dr.bib
@@ -278,14 +278,14 @@ License: GPL-3+
 
 Files: scripts/dynare.el
 Copyright: 2010 Yannick Kalantzis
-           2019-2023 Dynare Team
+           2019-2025 Dynare Team
 License: GPL-3+
 
 Files: mex/sources/gensylv/gensylv.cc
        mex/sources/libkorder/kord/* mex/sources/libkorder/sylv/*
        mex/sources/libkorder/tl/* mex/sources/libkorder/utils/*
 Copyright: 2004-2011 Ondra Kamenik
-           2019-2023 Dynare Team
+           2019-2025 Dynare Team
 License: GPL-3+
 
 Files: mex/sources/sobol/sobol.f08
@@ -294,11 +294,11 @@ Copyright: 2004-2009 John Burkardt
 License: LGPL-3+
 
 Files: preprocessor/doc/macroprocessor/*
-Copyright: 2008-2021 Dynare Team
+Copyright: 2008-2024 Dynare Team
 License: CC-BY-SA-4.0
 
 Files: preprocessor/doc/preprocessor/*
-Copyright: 2007-2019 Dynare Team
+Copyright: 2007-2023 Dynare Team
 License: CC-BY-SA-4.0
 
 Files: contrib/ms-sbvar/utilities_dw/*

matlab/+identification/run.m

@@ -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 = [];

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/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

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/extended_path.m

@@ -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.use_first_order_solution_as_initial_guess, 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)

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.
 %
@@ -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;
@@ -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.
 %
@@ -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.use_first_order_solution_as_initial_guess, 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.use_first_order_solution_as_initial_guess, 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);

matlab/ep/setup_stochastic_perfect_foresight_model_solver.m

@@ -5,7 +5,7 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_)
 %  o options_               [struct]    Dynare's options structure
 %  o oo_                    [struct]    Dynare's results structure
 
-% Copyright © 2013-2023 Dynare Team
+% Copyright © 2013-2025 Dynare Team
 %
 % This file is part of Dynare.
 %