Compare revisions

353a21ca · 87abdff3 · c62d9db0 · 334eaaa0 · d2667f6e · 152fb491
--- a/.clang-format
+++ b/.clang-format
@@ -8,6 +8,12 @@
 #  https://clang.llvm.org/docs/ClangFormatStyleOptions.html
 # Effective configuration can be obtained with:
 #  clang-format --dump-config
+# The RemoveParentheses and RemoveSemicolon option are not permanently set,
+# because the clang-format manual (as of version 19) states that these
+# options can lead to incorrect formatting and thus their result should be
+# carefully reviewed.
 Language: Cpp
 Standard: c++20
 ColumnLimit: 100
@@ -23,7 +29,8 @@ InsertNewlineAtEOF: true
 PackConstructorInitializers: NextLine
 PPIndentWidth: 1
 PointerAlignment: Left
-RemoveSemicolon: true
+# RemoveParentheses: ReturnStatement
+# RemoveSemicolon: true
 SpaceAfterTemplateKeyword: false
 SpaceBeforeParens: ControlStatements
 SpaceBeforeCpp11BracedList: true
--- a/.gitlab-ci.yml
+++ b/.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"
@@ -142,6 +142,20 @@ pkg_macOS_arm64:
    expire_in: 3 days
  needs: [ "build_doc" ]
+pkg_matlab_online:
+  stage: pkg
+  script:
+    - meson rewrite kwargs set project / version "$VERSION"
+    - cd scripts/matlab-online && ./packageDynare.sh
+  tags:
+    - matlab-online
+  artifacts:
+    paths:
+      - scripts/matlab-online/*.mltbx
+    expire_in: 3 days
+  needs: []
+  when: manual
 test_matlab:
  stage: test
  script:
@@ -162,6 +176,7 @@ test_old_matlab:
    paths:
      - build-old-matlab/meson-logs/testlog.txt
    when: always
+  needs: []
  when: manual
 test_octave:
@@ -188,7 +203,7 @@ test_clang_tidy:
  stage: test
  script:
    # Hack needed for meson < 1.6.0 which only looks for unversioned clang-tidy
-    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-16 ~/.local/bin/clang-tidy
+    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-19 ~/.local/bin/clang-tidy
    - export PATH="$HOME/.local/bin:$PATH"
    - meson setup -Dbuild_for=octave build-clang-tidy
    - ninja -C build-clang-tidy clang-tidy

--- a/NEWS.md
+++ b/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)
 ===========================================

--- a/README.md
+++ b/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

--- a/doc/manual/source/bibliography.rst
+++ b/doc/manual/source/bibliography.rst
@@ -6,6 +6,7 @@ Bibliography
 * Abramowitz, Milton and Irene A. Stegun (1964): “Handbook of Mathematical Functions”, Courier Dover Publications.
 * Adjemian, Stéphane, Matthieu Darracq Parriès and Stéphane Moyen (2008): “Towards a monetary policy evaluation framework”, *European Central Bank Working Paper*, 942.
+* Adjemian, Stéphane and Michel Juillard (2025): “Stochastic Extended Path”, *Dynare Working Papers*, 84, CEPREMAP.
 * Aguiar, Mark and Gopinath, Gita (2004): “Emerging Market Business Cycles: The Cycle is the Trend,” *NBER* Working Paper, 10734.
 * Amisano, Gianni and Tristani, Oreste (2010): “Euro area inflation persistence in an estimated nonlinear DSGE model”, *Journal of Economic Dynamics and Control*, 34(10), 1837–1858.
 * Andreasen, Martin M., Jesús Fernández-Villaverde, and Juan Rubio-Ramírez (2018): “The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications,” *Review of Economic Studies*, 85(1), 1-49.
@@ -65,7 +66,7 @@ Bibliography
 * Koopman, S. J. and J. Durbin (2003): “Filtering and Smoothing of State Vector for Diffuse State Space Models,” *Journal of Time Series Analysis*, 24(1), 85–98.
 * Kuntsevich, Alexei V. and Franz Kappel (1997): “SolvOpt - The solver for local nonlinear optimization problems (version 1.1, Matlab, C, FORTRAN)”, University of Graz, Graz, Austria.
 * Laffargue, Jean-Pierre (1990): “Résolution d’un modèle macroéconomique avec anticipations rationnelles”, *Annales d’Économie et Statistique*, 17, 97–119.
-* Liu, Jane and Mike West (2001): “Combined parameter and state estimation in simulation-based filtering”, in *Sequential Monte Carlo Methods in Practice*, Eds. Doucet, Freitas and Gordon, Springer Verlag.
+* Liu, Jane and Mike West (2001): “Combined parameter and state estimation in simulation-based filtering”, in *Sequential Monte Carlo Methods in Practice*, Eds. Doucet, Freitas and Gordon, Springer Verlag, Chapter 10, 197-223.
 * Murray, Lawrence M., Emlyn M. Jones and John Parslow (2013): “On Disturbance State-Space Models and the Particle Marginal Metropolis-Hastings Sampler”, *SIAM/ASA Journal on Uncertainty Quantification*, 1, 494–521.
 * Mutschler, Willi (2015): “Identification of DSGE models - The effect of higher-order approximation and pruning“, *Journal of Economic Dynamics & Control*, 56, 34-54.
 * Mutschler, Willi (2018): “Higher-order statistics for DSGE models”, *Econometrics and Statistics*, 6(C), 44-56.

--- a/doc/manual/source/conf.py
+++ b/doc/manual/source/conf.py
 # -*- coding: utf-8 -*-
-# Copyright © 2018-2024 Dynare Team
+# Copyright © 2018-2025 Dynare Team
 #
 # This file is part of Dynare.
 #
@@ -34,7 +34,7 @@ html_static_path = ['_static']
 master_doc = 'index'
 project = u'Dynare'
-copyright = u'1996–2024 Dynare Team'
+copyright = u'1996–2025 Dynare Team'
 author = u'Dynare Team'
 add_function_parentheses = False

--- a/doc/manual/source/index.rst
+++ b/doc/manual/source/index.rst
@@ -26,7 +26,7 @@ The following people used to be members of the team:
 * Ferhat Mihoubi
 * George Perendia
-Copyright © 1996-2024, Dynare Team.
+Copyright © 1996-2025, Dynare Team.
 Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.

--- a/doc/manual/source/installation-and-configuration.rst
+++ b/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.2.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.
@@ -24,6 +24,7 @@ The following optional extensions are also useful to benefit from
 extra features, but are in no way required:
 * If under MATLAB: the 
  * Optimization Toolbox (providing various optimizers 
    like ``fminsearch``, ``fmincon``, or ``fminunc``, used in e.g. ``mode_compute``,
    ``opt_algo`` or ``ALGO``), 
@@ -99,6 +100,19 @@ 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
 ^^^^^^^^^^^
@@ -108,19 +122,14 @@ 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`).
+To deinstall Dynare, simply delete the folder where you installed the program. The package installer does 
+not put any files anywhere else in the system.
 With Octave
 ^^^^^^^^^^^
@@ -141,6 +150,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
 ----------
@@ -333,3 +344,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
--- a/doc/manual/source/running-dynare.rst
+++ b/doc/manual/source/running-dynare.rst
@@ -188,10 +188,13 @@ by the ``dynare`` command.
        Instructs Dynare to no create a logfile of this run in
        ``FILENAME.log.`` The default is to create the logfile.
-    .. option:: output=second|third
+    .. option:: output=first|second|third
       Instructs the preprocessor to output derivatives of the dynamic model at
-       least up to the given order.
+       least up to the given order. The `first` option is useful in
+       larger models when debugging steady state computation, because it allows
+       overriding the default computation and output of dynamic second order derivatives 
+       in case of the mod-file not containing commands for further computations.
    .. option:: language=matlab|julia

--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -54,7 +54,7 @@ are introduced by ``/*`` and terminated by ``*/``.
    */
 Note that these comment marks should not be used in native MATLAB code regions
-where the `%` should be preferred instead to introduce a comment. In a
+where the ``%`` should be preferred instead to introduce a comment. In a
 ``verbatim`` block, see :ref:`verbatim`, this would result in a crash since
 ``//`` is not a valid MATLAB statement).
@@ -86,10 +86,10 @@ observed:
  (see :ref:`macro-exp`);
 * VARIABLE_NAME (sometimes VAR_NAME) indicates a variable name
  starting with an alphabetical character and can’t contain:
-  ‘()+-\*/^=!;:@#.’ or accentuated characters;
+  ``()+-\*/^=!;:@#.`` or accentuated characters;
 * PARAMETER_NAME (sometimes PARAM_NAME) indicates a parameter name
  starting with an alphabetical character and can’t contain:
-  ‘()+-\*/^=!;:@#.’ or accentuated characters;
+  ``()+-\*/^=!;:@#.`` or accentuated characters;
 * LATEX_NAME (sometimes TEX_NAME) indicates a valid
  LaTeX expression in math mode (not including the
  dollar signs);
@@ -99,12 +99,12 @@ observed:
  extension or if the filename contains a non-alphanumeric character;
 * QUOTED_STRING indicates an arbitrary string enclosed between (single)
  quotes;
-* DATE indicates a time period which can be either a year (e.g. `2024Y` or
+* DATE indicates a time period which can be either a year (e.g. ``2024Y`` or
-  `2024A`), a half-year (`2024S1` or `2024H1`), a quarter (`2024Q2`) or a month
+  ``2024A``), a half-year (``2024S1`` or ``2024H1``), a quarter (``2024Q2``) or a month
-  (`2024M3`) (see :ref:`dates in a mod file`). Optionally, the time period can
+  (``2024M3``) (see :ref:`dates in a mod file`). Optionally, the time period can
  be followed by a plus sign and a number of periods, in which case the date is
-  shifted accordingly (e.g. `2023Q1+6` is accepted and is equivalent to
+  shifted accordingly (e.g. ``2023Q1+6`` is accepted and is equivalent to
-  `2024Q3`).
+  ``2024Q3``).
 .. _var-decl:
@@ -959,6 +959,15 @@ The model is declared inside a ``model`` block:
        MODEL_EXPRESSION;
+    .. warning::
+        In Dynare, only equality signs can delineate the left and right-hand side of an 
+        equation. If Dynare encounters an expression like ``a>=b``, this will therefore not 
+        define an inequality constraint. Rather, it is interpreted as the homogenous equation
+        `(a>=b)=0;`, i.e., the Boolean `(a>=b)` must evaluate to 0. Inequality constraints 
+        in Dynare instead need to be set up either via OccBin or as mixed complementarity problems.
    |br| Inside the model block, Dynare allows the creation of
    *model-local variables*, which constitute a simple way to share a
    common expression between several equations. The syntax consists
@@ -2981,8 +2990,9 @@ Finding the steady state with Dynare nonlinear solver
           ``5``
                Newton algorithm with a sparse Gaussian elimination (SPE)
-                solver at each iteration (requires ``bytecode`` option, see
+                solver at each iteration. This algorithm requires the
-                :ref:`model-decl`).
+                :opt:`bytecode` option. The :ref:`markowitz <steady_markowitz>`
+                option can be used to control the behaviour of the algorithm.
           ``6``
@@ -3037,14 +3047,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
+                the :bck:`model` block or :comm:`model_options` command. Also
-                decomposition is performed as if ``mfs=3`` had been passed to
+                note that for those models, the block decomposition is
-                the ``model`` block, and the decomposition is slightly
+                performed as if ``mfs=3`` had been passed to the :bck:`model`
-                different because it is computed in a time-recursive fashion
+                block or :comm:`model_options` command, and the decomposition
-                (*i.e.* in such a way that the simulation is meant to be done
+                is slightly different because it is computed in a
-                with the outer loop on periods and the inner loop on blocks;
+                time-recursive fashion (*i.e.* in such a way that the
-                while for models with both leads and lags, the outer loop is on
+                simulation is meant to be done with the outer loop on periods
-                blocks and the inner loop is 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 +3572,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``
+    in conjunction with the ``bytecode`` option of the :bck:`model`
-    block.
+    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:
@@ -3772,13 +3785,20 @@ speed-up on large models.
               Use a Newton algorithm with a Generalized Minimal Residual
               (GMRES) solver at each iteration, applied on the stacked system
-               of all equations in all periods.
+               of all equations in all periods. The following options can be
+               used to control the behaviour of the algorithm:
+               :opt:`preconditioner <preconditioner = OPTION>`, :opt:`iter_tol
+               <iter_tol = DOUBLE>`, :opt:`iter_maxit <iter_maxit = INTEGER>`,
+               :opt:`gmres_restart <gmres_restart = INTEGER>`.
           ``3``
               Use a Newton algorithm with a Stabilized Bi-Conjugate Gradient
               (BiCGStab) solver at each iteration, applied on the stacked
-               system of all equations in all periods.
+               system of all equations in all periods. The following options
+               can be used to control the behaviour of the algorithm:
+               :opt:`preconditioner <preconditioner = OPTION>`, :opt:`iter_tol
+               <iter_tol = DOUBLE>`, :opt:`iter_maxit <iter_maxit = INTEGER>`.
           ``4``
@@ -3792,8 +3812,11 @@ speed-up on large models.
               Use the Laffargue-Boucekkine-Juillard (LBJ) algorithm proposed
               in *Juillard (1996)* on top of a sparse Gaussian elimination
               (SPE) solver. The latter takes advantage of the similarity of
-               the Jacobian across periods when searching for the pivots
+               the Jacobian across periods when searching for the pivots. This
-               (requires ``bytecode`` option, see :ref:`model-decl`).
+               algorithm requires the :opt:`bytecode` option. The following
+               options can be used to control the behaviour of the algorithm:
+               :ref:`markowitz <dynamic_markowitz>`,
+               :opt:`minimal_solving_periods <minimal_solving_periods = INTEGER>`.
           ``6``
@@ -3816,6 +3839,109 @@ speed-up on large models.
               trigger the computation of the solution with a trust
               region algorithm.
+    .. option:: preconditioner = OPTION
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option specifies which preconditioner will be used
+       in combination with the iterative sparse linear solver (either GMRES or
+       BiCGStab). Possible values for OPTION are:
+           ``umfiter``
+               At the first iteration of the nonlinear Newton solver, compute
+               the full LU decomposition with complete pivoting of the linear
+               system (and use it to solve that first iteration rather than
+               using the iterative solver). This LU decomposition is then used
+               as the preconditioner in further Newton iterations. Inspired
+               from TROLL’s option with the same name.
+           ``iterstack``
+               Compute the LU decomposition with complete pivoting for only a
+               few simulation periods within the stacked Jacobian (which is a
+               block tridiagonal matrix). Then repeat that LU decomposition
+               over the block diagonal to construct a preconditioner for the
+               linear system with all simulation periods. If the total number
+               of simulations periods is not a multiple of the number of
+               periods used for the small LU, then an additional LU is computed
+               for the remainder. The following options can be used to control
+               the construction of this preconditioner:
+               :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>`,
+               :opt:`iterstack_nperiods <iterstack_nperiods = INTEGER>`,
+               :opt:`iterstack_nlu <iterstack_nlu = INTEGER>`,
+               :opt:`iterstack_relu <iterstack_relu = DOUBLE>`.
+               Inspired from TROLL’s solver with the same name.
+           ``ilu``
+               Use an incomple LU decomposition as the preconditioner,
+               recomputed at every iteration of the nonlinear Newton solver.
+       |br| Default value is ``umfiter``.
+    .. option:: iter_tol = DOUBLE
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option controls the relative tolerance of the
+       iterative linear solver (either GMRES or BiCGStab). It corresponds to
+       the ``tol`` option of the ``gmres`` and ``bicgstab`` MATLAB/Octave
+       functions. Note that the perfect foresight solver uses an *absolute*
+       tolerance for determining convergence, so this option should be used
+       with care, and the default is meant to suit most situations.
+       Default: the value of the :opt:``tolf <tolf = DOUBLE>`` option, divided
+       by 10 times the infinite norm of the right-hand side of the linear system.
+    .. option:: iter_maxit = INTEGER
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2`` or ``3``, this option controls the maximum number of iterations of
+       the iterative linear solver (either GMRES or BiCGStab). It corresponds
+       to the ``maxit`` option of the ``gmres`` and ``bicgstab`` MATLAB/Octave
+       functions. It should not be confused with the :opt:`maxit <maxit = INTEGER>`
+       option, which controls the (outer) nonlinear Newton loop, while
+       ``iter_maxit`` controls the (inner) linear loop. Default: ``500``.
+    .. option:: gmres_restart = INTEGER
+       When :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` is equal to
+       ``2``, this option controls the number of iterations before restart of
+       the GMRES algorithm. It corresponds to the ``restart`` option of the
+       ``gmres`` MATLAB/Octave function. Default: ``100``.
+    .. option:: iterstack_maxlu = INTEGER
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the maximum size of the matrix for which the
+       small LU will computed. The actual size of the matrix will be determined
+       by the largest number of periods that, multiplied by the number of
+       equations, is less than the value of the option. Note that when combined
+       with block decomposition (see :opt:`block`), blocks for which the
+       whole stacked system is less than this option will be solved using a
+       regular LU decomposition instead of the iterative linear solver.
+       Default: ``20000``.
+    .. option:: iterstack_nperiods = INTEGER
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the number of periods used for the small LU.
+       If nonzero, this option overrides the :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>`
+       and :opt:`iterstack_nlu <iterstack_nlu = INTEGER>` options. Default: ``0``.
+    .. option:: iterstack_nlu = INTEGER
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, specifies the number of times that the small LU should be
+       repeated in the large preconditioner. If nonzero, this option overrides
+       the :opt:`iterstack_maxlu <iterstack_maxlu = INTEGER>` option. Default:
+       ``0``.
+    .. option:: iterstack_relu = DOUBLE
+       When :opt:`preconditioner <preconditioner = OPTION>` is equal to
+       ``iterstack``, controls the relative position of the small LU within the
+       whole stacked system. Must be a number between ``0`` and ``1``. Default:
+       ``0.5``.
    .. option:: robust_lin_solve
       Triggers the use of a robust linear solver for the default
@@ -3825,9 +3951,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
+       and static models (when neither the ``block`` nor the ``bytecode``
-       of the ``model`` block is specified); for those models, the values
+       option of the :bck:`model` block or :comm:`model_options` command is specified);
-       ``12`` and ``14`` are especially relevant.
+       for those models, the values ``12`` and ``14`` are especially relevant.
    .. option:: no_homotopy
@@ -3868,6 +3994,11 @@ speed-up on large models.
       the homotopy procedure has been able to find a solution, then the approximate
       solution returned is :math:`\frac{y(s^*)-y(0)}{s^*}`.
+       If linearization is triggered, the variable
+       :mvar:`oo_.deterministic_simulation.homotopy_linearization` is set, and
+       the simulation corresponding to share :math:`s^*` is stored in
+       :mvar:`oo_.deterministic_simulation.sim1`.
    .. option:: homotopy_marginal_linearization_fallback [= DOUBLE]
       Whenever the homotopy procedure is not able to find a solution for 100%
@@ -3887,6 +4018,13 @@ speed-up on large models.
       of :math:`\epsilon` is ``0.01`` by default, but can be modified by
       passing some other value to the option.
+       If marginal linearization is triggered, the variable
+       :mvar:`oo_.deterministic_simulation.homotopy_marginal_linearization` is
+       set. Moreover, the simulation corresponding to share :math:`s^*` is
+       stored in :mvar:`oo_.deterministic_simulation.sim1`, and the one
+       corresponding to share :math:`s^*-\epsilon` is stored in
+       :mvar:`oo_.deterministic_simulation.sim2`.
    .. option:: homotopy_max_completion_share = DOUBLE
       Instructs Dynare, within the homotopy procedure, to not try to compute
@@ -3896,7 +4034,7 @@ speed-up on large models.
       ``homotopy_marginal_linearization_fallback`` option. It is typically
       used in situations where it is known that homotopy will fail to go
       beyond a certain point, so as to save computing time, while at the same
-       time getting an approximate solution.
+       time getting an approximate solution. Default: ``1``.
    .. option:: homotopy_exclude_varexo = (VARIABLE_NAME...)
@@ -3904,11 +4042,13 @@ speed-up on large models.
       procedure, *i.e.* which must be kept at their value corresponding to
       100% of the shock during all homotopy iterations.
+    .. _dynamic_markowitz:
    .. option:: markowitz = DOUBLE
-       Value of the Markowitz criterion, used to select the
+       Value of the Markowitz criterion, used to select the pivot (see
-       pivot. Only used when ``stack_solve_algo = 5``. Default:
+       :ref:`markowitz <steady_markowitz>` for more details). Only used when
-       ``0.5``.
+       ``stack_solve_algo = 5``. Default: ``0.5``.
    .. option:: minimal_solving_periods = INTEGER
@@ -4003,7 +4143,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
@@ -4012,8 +4152,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
+       ``endval`` or a subsequent ``steady``. Only available with
-       ``stack_solve_algo==0`` or ``stack_solve_algo==7``.
+       ``stack_solve_algo`` option equal to ``0`` or ``7``.
    .. option:: steady_solve_algo = INTEGER
@@ -4048,13 +4188,16 @@ speed-up on large models.
    *Output*
    The simulated endogenous variables are available in global matrix
-    ``oo_.endo_simul``.
+    :mvar:`oo_.endo_simul`.
    If any of the ``first_simulation_period`` or ``last_simulation_period``
    option was passed to the preceding :comm:`perfect_foresight_setup` command,
    a time series object containing both endogenous and exogenous variables is
    stored in the workspace variable :mvar:`Simulated_time_series`.
+    The variable :mvar:`oo_.deterministic_simulation.status` indicates whether
+    the simulation was successful or not.
 .. command:: simul ;
             simul (OPTIONS...);
@@ -4132,6 +4275,42 @@ speed-up on large models.
    variables, so in practice it is either 1 or 0 (the latter value corresponds
    to a purely backward or static model).
+.. matvar:: oo_.deterministic_simulation.status
+    |br| Set to ``true`` by the :comm:`perfect_foresight_solver` command if the
+    simulation succeeded, otherwise set to ``false``.
+.. matvar:: oo_.deterministic_simulation.homotopy_linearization
+    |br| Set to ``true`` by the :comm:`perfect_foresight_solver` command if
+    linearization has been used to compute an approximate solution.
+.. matvar:: oo_.deterministic_simulation.homotopy_marginal_linearization
+    |br| Set to ``true`` by the :comm:`perfect_foresight_solver` command if
+    marginal linearization has been used to compute an approximate solution.
+.. matvar:: oo_.deterministic_simulation.sim1
+    |br| Set by the :comm:`perfect_foresight_solver` command if either
+    linearization or marginal linearization has been used to compute an
+    approximate solution. This structure contains the simulation corresponding
+    to the greatest share of the shocks for which an exact solution could be
+    computed. The subfield ``homotopy_completion_share`` contains that share.
+    The subfields ``endo_simul``, ``exo_simul``, ``steady_state`` and
+    ``exo_steady_state`` respectively contain the path of endogenous, the path
+    of exogenous, the steady state of endogenous and the steady state of
+    exogenous for that simulation (with the same conventions as the fields of
+    of the same name in ``oo_``).
+.. matvar:: oo_.deterministic_simulation.sim2
+    |br| Set by the :comm:`perfect_foresight_solver` command if marginal
+    linearization has been used to compute an approximate solution. This
+    structure contains the simulation corresponding to a share marginally
+    smaller that the one in :mvar:`oo_.deterministic_simulation.sim1`. The
+    subfields are the same as in :mvar:`oo_.deterministic_simulation.sim1`.
 Perfect foresight with expectation errors
 -----------------------------------------
@@ -4565,6 +4744,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
@@ -5286,12 +5603,13 @@ which is described below.
       If order is greater than ``0`` Dynare uses a gaussian
       quadrature to take into account the effects of future
-       uncertainty. If ``order`` :math:`=S` then the time series for
+       uncertainty; this is called *stochastic* extended path, see *Adjemian
+       and Juillard (2025)*. If ``order`` :math:`=S` then the time series for
       the endogenous variables are generated by assuming that the
       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
@@ -5309,6 +5627,13 @@ which is described below.
       rate or a model with irreversible investment). For specifying the
       necessary complementarity conditions, see :opt:`lmmcp`.
+    .. option:: use_first_order_solution
+       Utilize the model simulation based on a first-order local
+       approximation as the initial guess for the nonlinear solver in
+       each period. If this is not applied, solution in previous
+       period is used.
 Typology and ordering of variables
 ----------------------------------
@@ -6934,13 +7259,10 @@ observed variables.
           ``11``
-                This is not strictly speaking an optimization
+                Currently not in use. The Liu and West (2020) filter that 
-                algorithm. The (estimated) parameters are treated as
+                used to be available under this option value is now triggered with 
-                state variables and estimated jointly with the
+                ``posterior_sampling_method='online'``.
-                original state variables of the model using a
-                nonlinear filter. The algorithm implemented in Dynare
-                is described in *Liu and West (2001)*, and works with
-                ``k`` order local approximations of the model.
           ``12``
@@ -7055,9 +7377,9 @@ observed variables.
    .. option:: prior_trunc = DOUBLE
-       Probability of extreme values of the prior density that is
+       Probability of extreme values of the prior density in each tail that is
       ignored when computing bounds for the parameters. Default:
-       ``1e-32``.
+       ``1e-10`` for ``posterior_sampling_method=slice`` and ``0`` otherwise .
    .. option:: huge_number = DOUBLE
@@ -7525,7 +7847,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
@@ -7642,6 +7964,15 @@ observed variables.
               sampler proposed by *Waggoner, Wu and Zha (2016)* instead of the
               standard Random-Walk Metropolis-Hastings.
+           ``'online'``
+                Instructs Dynare to treat the (estimated) parameters as
+                state variables and estimate them jointly with the
+                original state variables of the model using a
+                nonlinear filter. The algorithm implemented in Dynare
+                is described in *Liu and West (2001)*, and works with
+                ``k`` order local approximations of the model.
    .. option:: posterior_sampler_options = (NAME, VALUE, ...)
       A list of NAME and VALUE pairs. Can be used to set options for
@@ -7864,7 +8195,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
@@ -7949,7 +8280,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
@@ -8503,7 +8834,8 @@ observed variables.
           ``'liu_west_delta'``
-               Set the value for delta for the Liu/West online filter. Default: ``0.99``.
+               Set the value for delta for the Liu/West online filter (``posterior_sampling_method='online'``). 
+               Default: ``0.99``.
           ``'unscented_alpha'``
@@ -9534,7 +9866,7 @@ Method of moments specific blocks
    * the first column contains the names of the endogenous variables
    * the second column contains the names of the exogenous variables
    * the third column contains a nested cell array that contains
-    the list of horizons, values and weights.
+      the list of horizons, values and weights.
 .. block:: matched_irfs_weights ;
@@ -11666,9 +11998,9 @@ the :comm:`bvar_forecast` command.
            oo_.conditional_forecast.uncond.FORECAST_MOMENT.VARIABLE_NAME
-    .. matvar:: forecasts.instruments
+    .. matvar:: oo_.conditional_forecast.instruments
-        Variable set by the ``conditional_forecast command``. Stores
+        Variable set by the ``conditional_forecast`` command. Stores
        the names of the exogenous instruments.
@@ -12269,7 +12601,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
@@ -12616,8 +12948,9 @@ Performing sensitivity analysis
        If equal to ``0``, ANOVA mapping (Type I error) If equal to
        ``1``, Screening analysis (Type II error). If equal to ``2``,
        Analytic derivatives (similar to Type II error, only valid
-        when identification=1). Default: ``1`` when
+        when identification=1). The ANOVA mapping requires the SS-ANOVA-R MATLAB 
-        ``identification=1``, ``0`` otherwise.
+        Toolbox available at https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip 
+        Default: ``1`` when ``identification=1``, ``0`` otherwise.
    .. option:: morris_nliv = INTEGER

--- a/doc/manual/source/time-series.rst
+++ b/doc/manual/source/time-series.rst
@@ -268,6 +268,39 @@ The dates class
            do4 = dates('Q',1950, 1);
            do5 = dates('D',1973, 1, 25);
+    |br|
+    A ``dates`` object with multiple elements can be considered a one-dimensional array of dates. Standard array operations can be applied to a ``dates`` object:
+     - square brackets can be used to concatenate dates objects::
+         >> A = dates('1938Q4');
+         >> B = dates('1945Q3');
+         >> C = [A, B];
+     - semicolons can be used to create ranges of dates::
+         >> A = dates('2009Q2');
+         >> B = A:A+2;
+         >> B
+         B = <dates: 2009Q2, 2009Q3, 2009Q4>
+     - objects can be indexed by an integer or a vector of integer::
+         >> B(1)
+         ans = <dates: 2009Q2>
+         >> B(end)
+         ans = <dates: 2009Q4>
+         >> B(1:2)
+         ans = <dates: 2009Q2, 2009Q3>
+    |br|
    A list of the available methods, by alphabetical order, is given
    below. Note that by default the methods do not allow in place
@@ -306,8 +339,10 @@ The dates class
    particularly useful if the methods are called in loops, since this
    saves the object instantiation overhead.
+    |br|
    .. datesmethod:: C = append (A, B)
-                     C = append_ (A, B)
+                     append_ (B)
        |br| Appends ``dates`` object ``B``, or a string that can be
        interpreted as a date, to the ``dates`` object ``A``. If ``B``
@@ -791,8 +826,8 @@ The dates class
    .. datesmethod:: C = pop (A)
                     C = pop (A, B)
-                     C = pop_ (A)
+                     pop_ ()
-                     C = pop_ (A, B)
+                     pop_ (B)
        |br| Pop method for ``dates`` class. If only one input is
        provided, the method removes the last element of a ``dates``
@@ -815,7 +850,7 @@ The dates class
    .. datesmethod:: C = remove (A, B)
-                     C = remove_ (A, B)
+                     remove_ (B)
        |br| Remove method for ``dates`` class. Both inputs have to be ``dates`` objects, removes dates in ``B`` from ``A``.
@@ -855,7 +890,7 @@ The dates class
    .. datesmethod:: B = sort (A)
-                     B = sort_ (A)
+                     sort_ ()
        |br| Sort method for ``dates`` objects. Returns a ``dates`` object
        with elements sorted by increasing order.
@@ -941,7 +976,7 @@ The dates class
    .. datesmethod:: B = unique (A)
-                     B = unique_ (A)
+                     unique_ ()
        |br| Overloads the MATLAB/Octave ``unique`` function. Returns
        a ``dates`` object with repetitions removed (only the last
@@ -1124,7 +1159,7 @@ The dseries class
    .. dseriesmethod:: A = abs (B)
-                       abs_ (B)
+                       abs_ ()
        |br| Overloads the ``abs()`` function for ``dseries``
        objects. Returns the absolute value of the variables in
@@ -1154,9 +1189,32 @@ The dseries class
                1973Q2 | 0.51222 | 0.4948
                1973Q3 | 0.99791 | 0.22677
+        *Example (in-place modification version)*
+            ::
+                >> ts0 = dseries(randn(3,2),'1973Q1',{'A1'; 'A2'},{'A_1'; 'A_2'});
+                >> ts0
+                ts0 is a dseries object:
+                       | A1       | A2
+                1973Q1 | -0.67284 | 1.4367
+                1973Q2 | -0.51222 | -0.4948
+                1973Q3 | 0.99791  | 0.22677
+                >> ts0.abs_();
+                >> ts0
+                ts0 is a dseries object:
+                       | abs(A1) | abs(A2)
+                1973Q1 | 0.67284 | 1.4367
+                1973Q2 | 0.51222 | 0.4948
+                1973Q3 | 0.99791 | 0.22677
    .. dseriesmethod:: [A, B] = align (A, B)
-                       align_ (A, B)
+                       align_ (B)
        If ``dseries`` objects ``A`` and ``B`` are defined on
        different time ranges, this function extends ``A`` and/or
@@ -1211,8 +1269,60 @@ The dseries class
                2001Q2 | 0.12801
+        *Example (in-place modification version)*
+            ::
+                >> ts0 = dseries(rand(5,1),dates('2000Q1')); % 2000Q1 -> 2001Q1
+                >> ts1 = dseries(rand(3,1),dates('2000Q4')); % 2000Q4 -> 2001Q2
+                >> ts0
+                ts0 is a dseries object:
+                       | Variable_1
+                2000Q1 | 0.80028
+                2000Q2 | 0.14189
+                2000Q3 | 0.42176
+                2000Q4 | 0.91574
+                2001Q1 | 0.79221
+                >> ts1
+                ts1 is a dseries object:
+                       | Variable_1
+                2000Q4 | 0.95949
+                2001Q1 | 0.65574
+                2001Q2 | 0.035712
+                >> align_(ts0, ts1);                         % 2000Q1 -> 2001Q2
+                >> ts0
+                ts0 is a dseries object:
+                       | Variable_1
+                2000Q1 | 0.80028
+                2000Q2 | 0.14189
+                2000Q3 | 0.42176
+                2000Q4 | 0.91574
+                2001Q1 | 0.79221
+                2001Q2 | NaN
+                >> ts1
+                ts1 is a dseries object:
+                       | Variable_1
+                2000Q1 | NaN
+                2000Q2 | NaN
+                2000Q3 | NaN
+                2000Q4 | 0.95949
+                2001Q1 | 0.65574
+                2001Q2 | 0.035712
    .. dseriesmethod:: C = backcast (A, B[, diff])
-                       backcast_ (A, B[, diff])
+                       backcast_ (B[, diff])
        Backcasts ``dseries`` object ``A`` with ``dseries`` object B's
        growth rates (except if the last optional argument, ``diff``,
@@ -1220,8 +1330,8 @@ The dseries class
        ``dseries`` objects must have the same frequency.
-    .. dseriesmethod:: B = baxter_king_filter (A, hf, lf, K)
+    .. dseriesmethod:: B = baxter_king_filter (A[, hf[, lf[, K]]])
-                       baxter_king_filter_ (A, hf, lf, K)
+                       baxter_king_filter_ ([hf[, lf[, K]]])
        |br| Implementation of the *Baxter and King* (1999) band pass
        filter for ``dseries`` objects. This filter isolates business
@@ -1267,7 +1377,7 @@ The dseries class
    .. dseriesmethod:: B = center (A[, geometric])
-                       center_ (A[, geometric])
+                       center_ ([geometric])
       |br| Centers variables in ``dseries`` object ``A`` around their
       arithmetic means, except if the optional argument ``geometric``
@@ -1276,7 +1386,7 @@ The dseries class
    .. dseriesmethod:: C = chain (A, B)
-                       chain_ (A, B)
+                       chain_ (B)
        |br| Merge two ``dseries`` objects along the time
        dimension. The two objects must have the same number of
@@ -1289,7 +1399,7 @@ The dseries class
            ::
-                >> ts = dseries([1; 2; 3; 4],dates(`1950Q1'))
+                >> ts = dseries([1; 2; 3; 4],dates('1950Q1'))
                ts is a dseries object:
@@ -1299,7 +1409,7 @@ The dseries class
                1950Q3 | 3
                1950Q4 | 4
-                >> us = dseries([3; 4; 5; 6],dates(`1950Q3'))
+                >> us = dseries([3; 4; 5; 6],dates('1950Q3'))
                us is a dseries object:
@@ -1321,6 +1431,25 @@ The dseries class
                1951Q1 | 5
                1951Q2 | 6
+        *Example (in-place modification version)*
+            ::
+                >> ts = dseries([1; 2; 3; 4],dates('1950Q1'))
+                >> us = dseries([3; 4; 5; 6],dates('1950Q3'))
+                >> ts.chain_(us);
+                >> ts
+                ts is a dseries object:
+                       | Variable_1
+                1950Q1 | 1
+                1950Q2 | 2
+                1950Q3 | 3
+                1950Q4 | 4
+                1951Q1 | 5
+                1951Q2 | 6
    .. dseriesmethod:: [error_flag, message ] = check (A)
@@ -1391,7 +1520,7 @@ The dseries class
    .. dseriesmethod:: B = cumprod (A[, d[, v]])
-                      cumprod_ (A[, d[, v]])
+                      cumprod_ ([d[, v]])
        |br| Overloads the MATLAB/Octave ``cumprod`` function for
        ``dseries`` objects. The cumulated product cannot be computed
@@ -1454,7 +1583,7 @@ The dseries class
    .. dseriesmethod:: B = cumsum (A[, d[, v]])
-                       cumsum (A[, d[, v]])
+                       cumsum_ ([d[, v]])
        |br| Overloads the MATLAB/Octave ``cumsum`` function for
        ``dseries`` objects. The cumulated sum cannot be computed if
@@ -1525,22 +1654,23 @@ The dseries class
                10Y | 10.1416
-    .. dseriesmethod:: B = detrend (A, m)
+    .. dseriesmethod:: B = detrend (A[, m])
-                       detrend_ (A, m)
+                       detrend_ ([m])
        |br| Detrends ``dseries`` object ``A`` with a fitted
-        polynomial of order ``m``. Note that each variable is
+        polynomial of order ``m``. Default value fir ``m`` is 0 (time
-        detrended with a different polynomial.
+        series are detrended by removing the average). Note that each
+        variable is detrended with a different polynomial.
    .. dseriesmethod:: B = dgrowth (A)
-                       dgrowth_ (A)
+                       dgrowth_ ()
        |br| Computes daily growth rates.
    .. dseriesmethod:: B = diff (A)
-                       diff_ (A)
+                       diff_ ()
        |br| Returns the first difference of ``dseries`` object ``A``.
@@ -1630,6 +1760,29 @@ The dseries class
                >> ts0 = dseries(rand(10,1));
                >> ts1 = ts0.exp();
+        *Exemple (in-place modification version)*
+            ::
+                >> ts0 = dseries(rand(3,1))
+                ts0 is a dseries object:
+                   | Variable_1
+                1Y | 0.82953
+                2Y | 0.84909
+                3Y | 0.37253
+                >> ts0.exp_();
+                >> ts0
+                ts0 is a dseries object:
+                   | Variable_1
+                1Y | 2.2922
+                2Y | 2.3375
+                3Y | 1.4514
    .. dseriesmethod:: C = extract (A, B[, ...])
@@ -1678,9 +1831,40 @@ The dseries class
                1973Q4 | -0.03705 | -0.35899  | 0.85838   | -1.4675  | -2.1666  | -0.62032
+    .. dseriesmethod:: fill_(name, v)
+       |br| Assign the value ``v`` to the variable ``name`` in a
+       dseries object. If ``name`` is a character row array, it should
+       correspond to an existing variable within the dseries
+       object. When ``v`` is a scalar, its value will be applied to
+       all periods uniformly. If ``v`` is a vector, its length must
+       match the number of observations in the dseries object.You can
+       invoke this method for a batch of variables by providing a 1 by
+       n cell array of character row arrays as the first argument. When
+       "v" is a row vector with n elements, the method will be applied
+       uniformly across all periods. If "v" is a matrix, it must have
+       n columns, and the number of rows should correspond to the
+       number of periods.
+       *Example*
+            ::
+                >> ts = dseries(rand(3,3));
+                >> ts.fill_({'Variable_1', 'Variable_3'}, [1 3]);
+                >> ts
+                ts is a dseries object:
+                   | Variable_1 | Variable_2 | Variable_3
+                1Y | 1          | 0.91338    | 3
+                2Y | 1          | 0.63236    | 3
+                3Y | 1          | 0.09754    | 3
    .. dseriesmethod:: f = firstdate (A)
-       |br| Returns the first period in ``dseries`` object ``A``.
+       |br| Returns the initial period in the ``dseries`` object ``A``.
    .. dseriesmethod:: f = firstobservedperiod (A)
@@ -1711,7 +1895,67 @@ The dseries class
                     4
-    .. dseriesmethod:: D = horzcat (A, B[, ...])
+    .. dseriesmethod:: l = ge (A, B)
+                       l = gt (A, B)
+       |br| Overloads the ``gt`` (>) and ``ge`` (>=) binary operators. Returns a logical array.
+       *Example*
+           ::
+                >> ts = dseries(randn(3,1))
+                ts is a dseries object:
+                   | Variable_1
+                1Y | -1.2075
+                2Y | 0.71724
+                3Y | 1.6302
+                >> ts>1
+                ans =
+                  3x1 logical array
+                   0
+                   0
+                   1
+                >> ds = dseries(randn(3,1))
+                ds is a dseries object:
+                   | Variable_1
+                1Y | 0.48889
+                2Y | 1.0347
+                3Y | 0.72689
+                >> ds>ts
+                ans =
+                  3x1 logical array
+                   1
+                   1
+                   0
+    .. dseriesmethod:: B = hdiff (A)
+                       hdiff_ ()
+       |br| Computes bi-annual differences.
+    .. dseriesmethod:: B = hgrowth (A)
+                       hgrowth_ ()
+       |br| Computes bi-annual growth rates.
+   .. dseriesmethod:: D = horzcat (A, B[, ...])
        |br| Overloads the ``horzcat`` MATLAB/Octave’s method for
        ``dseries`` objects. Returns a ``dseries`` object ``D``
@@ -1747,7 +1991,7 @@ The dseries class
    .. dseriesmethod:: B = hpcycle (A[, lambda])
-                       hpcycle_ (A[, lambda])
+                       hpcycle_ ([lambda])
        |br| Extracts the cycle component from a ``dseries`` ``A``
        object using the *Hodrick and Prescott (1997)* filter and
@@ -1788,7 +2032,7 @@ The dseries class
    .. dseriesmethod:: B = hptrend (A[, lambda])
-                       hptrend_ (A[, lambda])
+                       hptrend_ ([lambda])
        |br| Extracts the trend component from a ``dseries`` A object
        using the *Hodrick and Prescott (1997)* filter and returns a
@@ -1885,10 +2129,12 @@ The dseries class
    .. dseriesmethod:: B = lag (A[, p])
-                       lag_ (A[, p])
+                       lag_ ([p])
-        |br| Returns lagged time series. Default value of integer scalar ``p``, the number
+        |br| Returns lagged time series. Default value of integer
-        of lags, is ``1``.
+        scalar ``p``, the number of lags, is ``1``. The `dseries`
+        class overloads the parentheses, so that `ts.lag(p)` is
+        equivalent to `ts(-p)`.
        *Example*
@@ -1953,7 +2199,7 @@ The dseries class
    .. dseriesmethod:: l = lastdate (B)
-        |br| Returns the last period in ``dseries`` object ``B``.
+        |br| Retrieves the final period from the ``dseries`` object ``B``.
        *Example*
@@ -1967,23 +2213,76 @@ The dseries class
    .. dseriesmethod:: f = lastobservedperiod (A)
-       |br| Returns the last period where all the variables in ``dseries`` object ``A`` are observed (non NaN).
+       |br| Returns the last period in which all variables of the
+       ``dseries`` object ``A`` are fully observed (i.e., contain no
+       NaN values).
    .. dseriesmethod:: f = lastobservedperiods (A)
-       |br| Returns for each variable the last period without missing
+       |br| Returns the last period without missing observations for
-       observations in ``dseries`` object ``A``. Output argument ``f`` is a
+       each variable in the ``dseries`` object ``A``. The output
-       structure, each field name is the name of a variable in ``A``, each field
+       argument ``f`` is a structure where each field name corresponds
-       content is a singleton ``date`` object.
+       to a variable in ``A``, and the content of each field is a
+       singleton ``date`` object.
+    .. dseriesmethod:: l = le (A, B)
+                       l = lt (A, B)
+       |br| Overloads the ``gt`` (<) and ``ge`` (<=) binary operators. Returns a logical array.
+       *Example*
+           ::
+                >> ts = dseries(randn(3,1))
+                ts is a dseries object:
+                   | Variable_1
+                1Y | -1.2075
+                2Y | 0.71724
+                3Y | 1.6302
+                >> ts<1
+                ans =
+                  3x1 logical array
+                   1
+                   1
+                   0
+                >> ds = dseries(randn(3,1))
+                ds is a dseries object:
+                   | Variable_1
+                1Y | 0.48889
+                2Y | 1.0347
+                3Y | 0.72689
+                >> ds<ts
+                ans =
+                  3x1 logical array
+                   0
+                   0
+                   1
    .. dseriesmethod:: B = lead (A[, p])
-                       lead_ (A[, p])
+                       lead_ ([p])
-        |br| Returns lead time series. Default value of integer scalar
+        |br| Returns a lead time series. The default value for the
-        ``p``, the number of leads, is ``1``. As in the ``lag``
+        integer scalar ``p``, which represents the number of leads, is
-        method, the ``dseries`` class overloads the parenthesis so
+        ``1``. Similar to the ``lag`` method, the ``dseries`` class
-        that ``ts.lead(p)`` is equivalent to ``ts(p)``.
+        overloads the parentheses, making ``ts.lead(p)`` equivalent to
+        ``ts(p)``.
        *Example*
@@ -2012,11 +2311,11 @@ The dseries class
        *Remark*
-        The overloading of the parenthesis for ``dseries`` objects,
+        The overload of parentheses for ``dseries`` objects simplifies
-        allows to easily create new ``dseries`` objects by
+        the creation of new ``dseries`` instances by enabling the
-        copying/pasting equations declared in the ``model`` block. For
+        direct copying and pasting of equations defined within the
-        instance, if an Euler equation is defined in the ``model``
+        ``model`` block. For example, if an Euler equation is
-        block::
+        specified in the ``model`` block,::
            model;
            ...
@@ -2057,7 +2356,7 @@ The dseries class
    .. dseriesmethod:: B = log (A)
-                       log_ (A)
+                       log_ ()
        |br| Overloads the MATLAB/Octave ``log`` function for
        ``dseries`` objects.
@@ -2070,37 +2369,36 @@ The dseries class
                >> ts1 = ts0.log();
    .. dseriesmethod:: B = mdiff (A)
-                       mdiff_ (A)
+                       mdiff_ ()
                       B = mgrowth (A)
-                       mgrowth_ (A)
+                       mgrowth_ ()
-       |br| Computes monthly differences or growth rates of variables in
+       |br| Calculates the monthly differences or growth rates of
-       ``dseries`` object ``A``.
+       variables in the ``dseries`` object ``A``.
    .. dseriesmethod:: B = mean (A[, geometric])
-        |br| Overloads the MATLAB/Octave ``mean`` function for
+        |br| This function overloads the MATLAB/Octave ``mean``
-        ``dseries`` objects. Returns the mean of each variable in
+        function specifically for ``dseries`` objects. It calculates
-        ``dseries`` object ``A``. If the second argument is ``true``
+        the mean for each variable within the ``dseries`` object
-        the geometric mean is computed, otherwise (default) the
+        ``A``. If the second argument is set to ``true``, the
-        arithmetic mean is reported.
+        geometric mean is calculated; otherwise, the arithmetic mean
+        is computed by default.
    .. dseriesmethod:: C = merge (A, B[, legacy])
-        |br| Merges two ``dseries`` objects ``A`` and ``B`` in
+        |br| Merges two ``dseries`` objects, ``A`` and ``B``, into a new
-        ``dseries`` object ``C``. Objects ``A`` and ``B`` need to have
+        ``dseries`` object ``C``. The objects ``A`` and ``B`` must share a
-        common frequency but can be defined on different time
+        common frequency, although they can cover different time
-        ranges. If a variable, say ``x``, is defined both in
+        ranges. If a variable, such as ``x``, exists in both ``dseries``
-        ``dseries`` objects ``A`` and ``B``, then the ``merge`` will
+        objects, the ``merge`` function will prioritize the definition
-        select the variable ``x`` as defined in the second input
+        from the second input, ``B``, while retaining the values from
-        argument, ``B``, except for the NaN elements in ``B`` if
+        ``A`` for any corresponding periods where ``B`` has NaN
-        corresponding elements in ``A`` (ie same periods) are well
+        values. This behavior can be altered by setting the optional
-        defined numbers. This behaviour can be changed by setting the
+        argument ``legacy`` to true, in which case the second variable
-        optional argument ``legacy`` equal to true, in which case the
+        will replace the first, even if it contains NaN values.
-        second variable overwrites the first one even if the second
-        variable has NaNs.
        *Example*
@@ -2148,30 +2446,29 @@ The dseries class
    .. dseriesmethod:: C = minus (A, B)
        |br| Overloads the MATLAB/Octave ``minus`` (``-``) operator
-        for ``dseries`` objects, element by element subtraction. If
+        for ``dseries`` objects, allowing for element-by-element
-        both ``A`` and ``B`` are ``dseries`` objects, they do not need
+        subtraction. When both ``A`` and ``B`` are ``dseries``
-        to be defined over the same time ranges. If ``A`` and ``B``
+        objects, they do not need to be defined over the same time
-        are ``dseries`` objects with :math:`T_A` and :math:`T_B`
+        ranges. If ``A`` and ``B`` have :math:`T_A` and :math:`T_B`
        observations and :math:`N_A` and :math:`N_B` variables, then
-        :math:`N_A` must be equal to :math:`N_B` or :math:`1` and
+        :math:`N_A` must equal :math:`N_B` or :math:`1`, and
-        :math:`N_B` must be equal to :math:`N_A` or :math:`1`. If
+        :math:`N_B` must equal :math:`N_A` or :math:`1`. If
-        :math:`T_A=T_B`, ``isequal(A.init,B.init)`` returns ``1`` and
+        :math:`T_A=T_B`, ``isequal(A.init,B.init)`` returns ``1``,
-        :math:`N_A=N_B`, then the ``minus`` operator will compute for
+        and :math:`N_A=N_B`, then the ``minus`` operator will compute
-        each couple :math:`(t,n)`, with :math:`1\le t\le T_A` and
+        for each pair :math:`(t,n)`, where :math:`1\le t\le T_A` and
-        :math:`1\le n\le N_A`,
+        :math:`1\le n\le N_A`, the operation
-        ``C.data(t,n)=A.data(t,n)-B.data(t,n)``. If :math:`N_B` is
+        ``C.data(t,n)=A.data(t,n)-B.data(t,n)``. If :math:`N_B` equals
-        equal to :math:`1` and :math:`N_A>1`, the smaller ``dseries``
+        :math:`1` and :math:`N_A>1`, the smaller ``dseries`` object
-        object (``B``) is “broadcast” across the larger ``dseries``
+        (``B``) is “broadcasted” across the larger ``dseries`` (``A``),
-        (``A``) so that they have compatible shapes, the ``minus``
+        ensuring compatible shapes for the subtraction of the variable
-        operator will subtract the variable defined in ``B`` from each
+        defined in ``B`` from each variable in ``A``. If ``B`` is a
-        variable in ``A``. If ``B`` is a double scalar, then the
+        double scalar, the ``minus`` method will subtract ``B`` from
-        method ``minus`` will subtract ``B`` from all the
+        all observations and variables in ``A``. If ``B`` is a row
-        observations/variables in ``A``. If ``B`` is a row vector of
+        vector of length :math:`N_A`, the ``minus`` method will
-        length :math:`N_A`, then the ``minus`` method will subtract
+        subtract ``B(i)`` from all observations of variable ``i``, for
-        ``B(i)`` from all the observations of variable ``i``, for
        :math:`i=1,...,N_A`. If ``B`` is a column vector of length
-        :math:`T_A`, then the ``minus`` method will subtract ``B``
+        :math:`T_A`, the ``minus`` method will subtract ``B`` from all
-        from all the variables.
+        the variables.
        *Example*
@@ -2218,15 +2515,16 @@ The dseries class
    .. dseriesmethod:: C = mpower (A, B)
-        |br| Overloads the MATLAB/Octave ``mpower`` (``^``) operator for ``dseries``
+        |br| Overloads the MATLAB/Octave ``mpower`` (``^``) operator
-        objects and computes element-by-element power. ``A`` is a
+        for ``dseries`` objects, performing element-wise
+        exponentiation. Given a ``dseries`` object ``A`` with ``N``
+        variables and ``T`` observations, if ``B`` is a real scalar,
+        then ``mpower(A,B)`` yields a ``dseries`` object ``C`` where
+        ``C.data(t,n) = A.data(t,n)^B``. If ``B`` is also a
        ``dseries`` object with ``N`` variables and ``T``
-        observations. If ``B`` is a real scalar, then ``mpower(A,B)``
+        observations, then ``mpower(A,B)`` produces a ``dseries``
-        returns a ``dseries`` object ``C`` with
+        object ``C`` such that ``C.data(t,n) =
-        ``C.data(t,n)=A.data(t,n)^C``. If ``B`` is a ``dseries``
+        A.data(t,n)^{C.data(t,n)}``.
-        object with ``N`` variables and ``T`` observations then
-        ``mpower(A,B)`` returns a ``dseries`` object ``C`` with
-        ``C.data(t,n)=A.data(t,n)^C.data(t,n)``.
        *Example*
@@ -2254,31 +2552,31 @@ The dseries class
    .. dseriesmethod:: C = mrdivide (A, B)
-        |br| Overloads the MATLAB/Octave ``mrdivide`` (``/``) operator for
+        |br| Overloads the MATLAB/Octave ``mrdivide`` (``/``) operator
-        ``dseries`` objects, element by element division (like the
+        for ``dseries`` objects, enabling element-wise division
-        ``./`` MATLAB/Octave operator). If both ``A`` and ``B`` are
+        similar to the ``./`` operator in MATLAB/Octave. When both
-        ``dseries`` objects, they do not need to be defined over the
+        ``A`` and ``B`` are ``dseries`` objects, they can have
-        same time ranges. If ``A`` and ``B`` are ``dseries`` objects
+        different time ranges. If ``A`` contains :math:`T_A`
-        with :math:`T_A` and :math:`T_B` observations and :math:`N_A`
+        observations and :math:`N_A` variables, and ``B`` has
-        and :math:`N_B` variables, then :math:`N_A` must be equal to
+        :math:`T_B` observations and :math:`N_B` variables, then
-        :math:`N_B` or :math:`1` and :math:`N_B` must be equal to
+        :math:`N_A` must equal :math:`N_B` or :math:`1`, and vice
-        :math:`N_A` or :math:`1`. If :math:`T_A=T_B`,
+        versa. If :math:`T_A=T_B` and ``isequal(A.init,B.init)``
-        ``isequal(A.init,B.init)`` returns ``1`` and :math:`N_A=N_B`,
+        returns ``1``, along with :math:`N_A=N_B`, the ``mrdivide``
-        then the ``mrdivide`` operator will compute for each couple
+        operator calculates for each pair :math:`(t,n)`, where
-        :math:`(t,n)`, with :math:`1\le t\le T_A` and :math:`1\le n\le
+        :math:`1\le t\le T_A` and :math:`1\le n\le N_A`, the value of
-        N_A`, ``C.data(t,n)=A.data(t,n)/B.data(t,n)``. If :math:`N_B`
+        ``C.data(t,n)=A.data(t,n)/B.data(t,n)``. If :math:`N_B` equals
-        is equal to :math:`1` and :math:`N_A>1`, the smaller
+        :math:`1` and :math:`N_A>1`, the smaller ``dseries`` object
-        ``dseries`` object (``B``) is “broadcast” across the larger
+        (``B``) is “broadcast” across the larger one (``A``) to ensure
-        ``dseries`` (``A``) so that they have compatible shapes. In
+        compatible shapes. In this case, the ``mrdivide`` operator
-        this case the ``mrdivide`` operator will divide each variable
+        divides each variable in ``A`` by the variable in ``B``,
-        defined in A by the variable in B, observation per
+        observation by observation. If ``B`` is a double scalar, then
-        observation. If B is a double scalar, then ``mrdivide`` will
+        ``mrdivide`` will divide all observations and variables in
-        divide all the observations/variables in ``A`` by ``B``. If
+        ``A`` by ``B``. If ``B`` is a row vector of length
-        ``B`` is a row vector of length :math:`N_A`, then ``mrdivide``
+        :math:`N_A`, then ``mrdivide`` will divide each observation of
-        will divide all the observations of variable ``i`` by
+        variable ``i`` by ``B(i)``, for :math:`i=1,...,N_A`. If ``B``
-        ``B(i)``, for :math:`i=1,...,N_A`. If ``B`` is a column vector
+        is a column vector of length :math:`T_A`, then ``mrdivide``
-        of length :math:`T_A`, then ``mrdivide`` will perform a
+        will perform an element-wise division of all variables by
-        division of all the variables by ``B``, element by element.
+        ``B``.
        *Example*
@@ -2307,62 +2605,62 @@ The dseries class
    .. dseriesmethod:: C = mtimes (A, B)
        |br| Overloads the MATLAB/Octave ``mtimes`` (``*``) operator
-        for ``dseries`` objects and the Hadammard product (the .*
+        for ``dseries`` objects, enabling element-wise multiplication
-        MATLAB/Octave operator). If both ``A`` and ``B`` are
+        similar to the ``.*`` operator in MATLAB/Octave. When both
-        ``dseries`` objects, they do not need to be defined over the
+        ``A`` and ``B`` are ``dseries`` objects, they can have
-        same time ranges. If ``A`` and ``B`` are ``dseries`` objects
+        different time ranges. If ``A`` contains :math:`T_A`
-        with :math:`T_A` and :math:`_B` observations and :math:`N_A`
+        observations and :math:`N_A` variables, and ``B`` has
-        and :math:`N_B` variables, then :math:`N_A` must be equal to
+        :math:`T_B` observations and :math:`N_B` variables, then
-        :math:`N_B` or :math:`1` and :math:`N_B` must be equal to
+        :math:`N_A` must equal :math:`N_B` or :math:`1`, and vice
-        :math:`N_A` or :math:`1`. If :math:`T_A=T_B`,
+        versa. If :math:`T_A=T_B` and ``isequal(A.init,B.init)``
-        ``isequal(A.init,B.init)`` returns ``1`` and :math:`N_A=N_B`,
+        returns ``1``, along with :math:`N_A=N_B`, the ``mtimes``
-        then the ``mtimes`` operator will compute for each couple
+        operator calculates for each pair :math:`(t,n)`, where
-        :math:`(t,n)`, with :math:`1\le t\le T_A` and :math:`1\le n\le
+        :math:`1\le t\le T_A` and :math:`1\le n\le N_A`, the value of
-        N_A`, ``C.data(t,n)=A.data(t,n)*B.data(t,n)``. If :math:`N_B`
+        ``C.data(t,n)=A.data(t,n)*B.data(t,n)``. If :math:`N_B` equals
-        is equal to :math:`1` and :math:`N_A>1`, the smaller
+        :math:`1` and :math:`N_A>1`, the smaller ``dseries`` object
-        ``dseries`` object (``B``) is “broadcast” across the larger
+        (``B``) is “broadcasted” across the larger one (``A``) to ensure
-        ``dseries`` (``A``) so that they have compatible shapes,
+        compatible shapes. In this case, the ``mtimes`` operator
-        ``mtimes`` operator will multiply each variable defined in
+        multiply each variable in ``A`` by the variable in ``B``,
-        ``A`` by the variable in ``B``, observation per
+        observation by observation. If ``B`` is a double scalar, then
-        observation. If ``B`` is a double scalar, then the method
+        ``mtimes`` will multiply all observations and variables in
-        ``mtimes`` will multiply all the observations/variables in
        ``A`` by ``B``. If ``B`` is a row vector of length
-        :math:`N_A`, then the ``mtimes`` method will multiply all the
+        :math:`N_A`, then ``mtimes`` will multiply each observation of
-        observations of variable ``i`` by ``B(i)``, for
+        variable ``i`` by ``B(i)``, for :math:`i=1,...,N_A`. If ``B``
-        :math:`i=1,...,N_A`. If ``B`` is a column vector of length
+        is a column vector of length :math:`T_A`, then ``mtimes``
-        :math:`T_A`, then the ``mtimes`` method will perform a
+        will perform an element-wise multiplication of all variables by
-        multiplication of all the variables by ``B``, element by
+        ``B``.
-        element.
    .. dseriesmethod:: B = nanmean (A[, geometric])
        |br| Overloads the MATLAB/Octave ``nanmean`` function for
-        ``dseries`` objects. Returns the mean of each variable in
+        ``dseries`` objects. Computes the mean of each variable in the
-        ``dseries`` object ``A`` ignoring the NaN values. If the
+        ``dseries`` object ``A``, excluding NaN values. If the second
-        second argument is ``true`` the geometric mean is computed,
+        argument is ``true``, the geometric mean is calculated;
-        otherwise (default) the arithmetic mean is reported.
+        otherwise, the default is to report the arithmetic mean.
    .. dseriesmethod:: B = nanstd (A[, geometric])
        |br| Overloads the MATLAB/Octave ``nanstd`` function for
-        ``dseries`` objects. Returns the standard deviation of each
+        ``dseries`` objects. This function calculates the standard
-        variable in ``dseries`` object ``A`` ignoring the NaN
+        deviation for each variable within the ``dseries`` object
-        values. If the second argument is ``true`` the geometric std
+        ``A``, while disregarding any NaN values. If the second
-        is computed, default value of the second argument is
+        argument is set to ``true``, the geometric standard deviation
+        will be computed; the default value for the second argument is
        ``false``.
    .. dseriesmethod:: C = ne (A, B)
        |br| Overloads the MATLAB/Octave ``ne`` (not equal, ``~=``)
-        operator. ``dseries`` objects ``A`` and ``B`` must have the
+        operator. The ``dseries`` objects ``A`` and ``B`` must contain
-        same number of observations (say, :math:`T`) and variables
+        the same number of observations (denoted as :math:`T`) and
-        (:math:`N`). The returned argument is a :math:`T` by :math:`N`
+        variables (denoted as :math:`N`). The output is a :math:`T` by
-        matrix of zeros and ones. Element :math:`(i,j)` of ``C`` is
+        :math:`N` matrix consisting of zeros and ones. The element
-        equal to ``1`` if and only if observation :math:`i` for
+        :math:`(i,j)` of the matrix ``C`` is equal to ``1`` if and
-        variable :math:`j` in ``A`` and ``B`` are not equal.
+        only if observation :math:`i` for variable :math:`j` in ``A``
+        and ``B`` are not equal.
        *Example*
@@ -2399,25 +2697,25 @@ The dseries class
    .. dseriesmethod:: B = onesidedhpcycle (A[, lambda[, init]])
-                       onesidedhpcycle_ (A[, lambda[, init]])
+                       onesidedhpcycle_ ([lambda[, init]])
        |br| Extracts the cycle component from a ``dseries`` ``A``
-        object using a one sided HP filter (with a Kalman filter) and
+        object using a one-sided HP filter (implemented with a Kalman
-        returns a ``dseries`` object, ``B``. The default value for
+        filter) and returns a ``dseries`` object, ``B``. The default
-        ``lambda``, the smoothing parameter, is ``1600``. By default,
+        value for ``lambda``, the smoothing parameter, is set to
-        if ``ìnit`` is not provided, the initial value is based on the
+        ``1600``. By default, if ``init`` is not provided, the initial
-        first two observations.
+        value is determined from the first two observations.
    .. dseriesmethod:: B = onesidedhptrend (A[, lambda[, init]])
                       onesidedhptrend_ (A[, lambda[, init]])
        |br| Extracts the trend component from a ``dseries`` ``A``
-        object using a one sided HP filter (with a Kalman filter) and
+        object using a one-sided HP filter (implemented with a Kalman
-        returns a ``dseries`` object, ``B``. The default value for
+        filter) and returns a ``dseries`` object, ``B``. The default
-        ``lambda``, the smoothing parameter, is ``1600``. By default,
+        value for ``lambda``, the smoothing parameter, is set to
-        if ``ìnit`` is not provided, the initial value is based on the
+        ``1600``. By default, if ``init`` is not provided, the initial
-        first two observations.
+        value is derived from the first two observations.
    .. dseriesmethod:: h = plot (A)
@@ -2425,27 +2723,28 @@ The dseries class
                       h = plot (A[, ...])
                       h = plot (A, B[, ...])
-        |br| Overloads MATLAB/Octave’s ``plot`` function for
+        |br| Overloads the MATLAB/Octave ``plot`` function for
-        ``dseries`` objects. Returns a MATLAB/Octave plot handle, that
+        ``dseries`` objects. This function returns a MATLAB/Octave
-        can be used to modify the properties of the plotted time
+        plot handle, which can be utilized to modify the properties of
-        series. If only one ``dseries`` object, ``A``, is passed as
+        the plotted time series. If a single ``dseries`` object,
-        argument, then the plot function will put the associated dates
+        ``A``, is provided as an argument, the plot function will
-        on the x-abscissa. If this ``dseries`` object contains only
+        place the corresponding dates on the x-axis. If this
-        one variable, additional arguments can be passed to modify the
+        ``dseries`` object contains only one variable, additional
-        properties of the plot (as one would do with the
+        arguments can be included to adjust the plot properties,
-        MATLAB/Octave’s version of the plot function). If ``dseries``
+        similar to how one would with MATLAB/Octave's original plot
-        object ``A`` contains more than one variable, it is not
+        function. However, if the ``dseries`` object ``A`` has more
-        possible to pass these additional arguments and the properties
+        than one variable, additional arguments cannot be passed, and
-        of the plotted time series must be modified using the returned
+        modifications to the plotted time series properties must be
-        plot handle and the MATLAB/Octave ``set`` function (see
+        done using the returned plot handle alongside the
-        example below). If two ``dseries`` objects, ``A`` and ``B``,
+        MATLAB/Octave ``set`` function (refer to the example
-        are passed as input arguments, the plot function will plot the
+        below). When two ``dseries`` objects, ``A`` and ``B``, are
-        variables in ``A`` against the variables in ``B`` (the number
+        passed as input arguments, the plot function will display the
-        of variables in each object must be the same otherwise an
+        variables in ``A`` against those in ``B`` (it is essential
-        error is issued). Again, if each object contains only one
+        that both objects contain the same number of variables;
-        variable, additional arguments can be passed to modify the
+        otherwise, an error will occur). Once more, if each object
-        properties of the plotted time series, otherwise the
+        includes only one variable, additional arguments can be
-        MATLAB/Octave ``set`` command has to be used.
+        utilized to alter the plotted time series properties;
+        otherwise, the MATLAB/Octave ``set`` command must be employed.
        *Example*
@@ -2486,7 +2785,7 @@ The dseries class
    .. dseriesmethod:: C = plus (A, B)
        |br| Overloads the MATLAB/Octave ``plus`` (``+``) operator for
-        ``dseries`` objects, element by element addition. If both
+        ``dseries`` objects, allowing for element-wise addition. When both
        ``A`` and ``B`` are ``dseries`` objects, they do not need to
        be defined over the same time ranges. If ``A`` and ``B`` are
        ``dseries`` objects with :math:`T_A` and :math:`T_B`
@@ -2495,28 +2794,28 @@ The dseries class
        :math:`N_B` must be equal to :math:`N_A` or :math:`1`. If
        :math:`T_A=T_B`, ``isequal(A.init,B.init)`` returns ``1`` and
        :math:`N_A=N_B`, then the ``plus`` operator will compute for
-        each couple :math:`(t,n)`, with :math:`1\le t\le T_A` and
+        each pair :math:`(t,n)`, with :math:`1\le t\le T_A` and
        :math:`1\le n\le N_A`,
        ``C.data(t,n)=A.data(t,n)+B.data(t,n)``. If :math:`N_B` is
        equal to :math:`1` and :math:`N_A>1`, the smaller ``dseries``
-        object (``B``) is “broadcast” across the larger ``dseries``
+        object (``B``) is “broadcasted” across the larger ``dseries``
-        (``A``) so that they have compatible shapes, the plus operator
+        (``A``) to ensure compatible shapes, the plus operator
        will add the variable defined in ``B`` to each variable in
        ``A``. If ``B`` is a double scalar, then the method ``plus``
        will add ``B`` to all the observations/variables in ``A``. If
        ``B`` is a row vector of length :math:`N_A`, then the ``plus``
        method will add ``B(i)`` to all the observations of variable
-        ``i``, for :math:`i=1,...,N_A`. If ``B`` is a column vector of
+        ``i``, for :math:`i=1,\ldots,N_A`. If ``B`` is a column vector of
        length :math:`T_A`, then the ``plus`` method will add ``B`` to
        all the variables.
    .. dseriesmethod:: C = pop (A[, B])
-                       pop_ (A[, B])
+                       pop_ ([B])
-        |br| Removes variable ``B`` from ``dseries`` object ``A``. By
+        |br| Removes the variable ``B`` from the ``dseries`` object
-        default, if the second argument is not provided, the last
+        ``A``. By default, if the second argument is not specified, the
-        variable is removed.
+        last variable is removed.
        *Example*
@@ -2535,26 +2834,25 @@ The dseries class
    .. dseriesmethod:: A = projection (A, info, periods)
-        |br| Projects variables in dseries object ``A``. ``info`` is
+        |br| Projects variables in the dseries object ``A``. The
-        is a :math:`n \times 3` cell array. Each row provides
+        ``info`` variable is a :math:`n \times 3` cell array, where
-        informations necessary to project a variable. The first column
+        each row contains essential information for projecting a
-        contains the name of variable (row char array). the second
+        variable. The first column holds the variable name (as a
-        column contains the name of the method used to project the
+        character array), while the second column indicates the
-        associated variable (row char array), possible values are
+        projection method used (also a character array). The possible
-        ``'Trend'``, ``'Constant'``, and ``'AR'``. Last column
+        values for this column are ``'Trend'``, ``'Constant'``, and
-        provides quantitative information about the projection. If the
+        ``'AR'``. The third column provides quantitative details
-        second column value is ``'Trend'``, the third column value is
+        related to the projection: if the second column is
-        the growth factor of the (exponential) trend. If the second
+        ``'Trend'``, the third column specifies the growth factor of
-        column value is ``'Constant'``, the third column value is the
+        the (exponential) trend; if ``'Constant'``, it indicates the
-        level of the variable. If the second column value is ``'AR'``,
+        variable's level; and if ``'AR'``, it denotes the
-        the third column value is the autoregressive parameter. The
+        autoregressive parameter. Variables can be projected using an
-        variables can be projected with an AR(p) model, if the third
+        AR(p) model if the third column contains a 1×p vector of
-        column contains a 1×p vector of doubles. The stationarity of
+        doubles. Note that the stationarity of the AR(p) model is not
-        the AR(p) model is not tested. The case of the constant
+        tested. For constant projections, one can use either `'Trend'`
-        projection, using the last value of the variable, is covered
+        with a growth factor of 1 or `'AR'` with an autoregressive
-        with 'Trend' and a growth factor equal to 1, or 'AR' with an
+        parameter of one (indicating a random walk). This projection
-        autoregressive parameter equal to one (random walk).  This
+        routine solely addresses exponential trends.
-        projection routine only deals with exponential trends.
        *Example*
@@ -2568,8 +2866,8 @@ The dseries class
    .. dseriesmethod:: B = qdiff (A)
                       B = qgrowth (A)
-                       qdiff_ (A)
+                       qdiff_ ()
-                       qgrowth_ (A)
+                       qgrowth_ ()
        |br| Computes quarterly differences or growth rates.
@@ -2603,12 +2901,14 @@ The dseries class
    .. dseriesmethod:: C = remove (A, B)
-                       remove_ (A, B)
+                       remove_ (B)
-        |br| If ``B`` is a row char array, the name of a variable, these methods
+        |br| If ``B`` is a row character array representing the name
-        are aliases for ``pop`` and ``pop_`` methods with two arguments. They
+        of a variable, these methods serve as aliases for the ``pop``
-        remove variable ``B`` from ``dseries`` object ``A``. To remove more than
+        and ``pop_`` methods that accept two arguments. They remove
-        one variable, one can pass a cell of row char arrays for ``B``.
+        the variable ``B`` from the ``dseries`` object ``A``. To
+        remove multiple variables, you can pass a cell array of row
+        character arrays for ``B``.
        *Example*
@@ -2624,11 +2924,11 @@ The dseries class
                2Y | 1          | 1
                3Y | 1          | 1
-            A shorter syntax is available: ``remove(ts,'Variable_2')``
+            A more concise syntax is available: ``remove(ts,
-            is equivalent to ``ts{'Variable_2'} = []`` (``[]`` can be
+            'Variable_2')``, which is equivalent to ``ts{'Variable_2'}
-            replaced by any empty object). This alternative syntax is
+            = []`` (where ``[]`` can be substituted with any empty
-            useful if more than one variable has to be removed. For
+            object). This alternative syntax proves useful when
-            instance::
+            removing multiple variables. For example::
                ts{'Variable_@2,3,4@'} = [];
@@ -2639,12 +2939,13 @@ The dseries class
    .. dseriesmethod:: B = rename (A, oldname, newname)
-                       rename_ (A, oldname, newname)
+                       rename_ (oldname, newname)
-        |br| Rename variable ``oldname`` to ``newname`` in ``dseries``
+        |br| Renames the variable ``oldname`` to ``newname`` in the
-        object ``A``. Returns a ``dseries`` object. If more than one
+        ``dseries`` object ``A``. This function returns a ``dseries``
-        variable needs to be renamed, it is possible to pass cells of
+        object. If multiple variables need to be renamed, you can
-        char arrays as second and third arguments.
+        provide cell arrays of row character arrays as the second and
+        third arguments.
        *Example*
@@ -2661,12 +2962,12 @@ The dseries class
    .. dseriesmethod:: C = rename (A, newname)
-                       rename_ (A, newname)
+                       rename_ (newname)
-        |br| Replace the names in ``A`` with those passed in the cell
+        |br| Replace the names in ``A`` with those specified in the
-        string array ``newname``. ``newname`` must have the same
+        cell of row character arrays ``newname``. The cell ``newname`` must contain
-        number of elements as ``dseries`` object ``A`` has
+        the same number of elements as there are variables in the
-        variables. Returns a ``dseries`` object.
+        ``dseries`` object ``A``.
        *Example*
@@ -2693,12 +2994,12 @@ The dseries class
    .. dseriesmethod:: B = round (A[, n])
-                       round_ (A[, n])
+                       round_ ([n])
-        |br| Rounds to the nearest decimal or integer. ``n`` is the
+        |br| Rounds each value to the nearest decimal or integer. The
-        precision parameter (number of decimals), default value is 0
+        parameter ``n`` specifies the precision (number of decimal
-        meaning that that by default the method rounds to the nearest
+        places), with a default value of 0, indicating that the method
-        integer.
+        will round to the nearest integer by default.
        *Example*
@@ -2722,11 +3023,12 @@ The dseries class
    .. dseriesmethod:: save (A, basename[, format])
-        |br| Overloads the MATLAB/Octave ``save`` function and saves
+        |br| Overloads the MATLAB/Octave ``save`` function to save the
-        ``dseries`` object ``A`` to disk. Possible formats are ``mat``
+        ``dseries`` object ``A`` to disk. The available formats
-        (this is the default), ``m`` (MATLAB/Octave script), and
+        include ``mat`` (default, MATLAB binary data file), ``m``
-        ``csv`` (MATLAB binary data file). The name of the file
+        (MATLAB/Octave script), and ``csv`` (comma-separated values
-        without extension is specified by ``basename``.
+        file). The base name of the file, excluding the extension, is
+        specified by ``basename``.
        *Example*
@@ -2773,13 +3075,14 @@ The dseries class
    .. dseriesmethod:: B = set_names(A, s1, s2, ...)
-        |br| Renames variables in ``dseries`` object ``A`` and returns
+        |br| Renames the variables in the ``dseries`` object ``A`` and
-        a ``dseries`` object ``B`` with new names ``s1``, ``s2``,
+        returns a new ``dseries`` object ``B`` with the updated names
-        ... The number of input arguments after the first one
+        ``s1``, ``s2``, and so forth. The number of input arguments
-        (``dseries`` object ``A``) must be equal to ``A.vobs`` (the
+        following the first one (the ``dseries`` object ``A``) must be
-        number of variables in ``A``). ``s1`` will be the name of the
+        equal to ``A.vobs`` (the total number of variables in
-        first variable in ``B``, ``s2`` the name of the second
+        ``A``). The name ``s1`` will correspond to the first variable
-        variable in ``B``, and so on.
+        in ``B``, ``s2`` to the second variable in ``B``, and this
+        pattern continues for the remaining variables.
        *Example*
@@ -2796,13 +3099,13 @@ The dseries class
    .. dseriesmethod:: [T, N ] = size(A[, dim])
-        Overloads the MATLAB/Octave’s ``size`` function. Returns the
+        Overloads the MATLAB/Octave ``size`` function to return the
-        number of observations in ``dseries`` object ``A``
+        number of observations in the ``dseries`` object ``A`` (i.e.,
-        (i.e. ``A.nobs``) and the number of variables
+        ``A.nobs``) as well as the number of variables (i.e.,
-        (i.e. ``A.vobs``). If a second input argument is passed, the
+        ``A.vobs``). If a second input argument is provided, the
-        ``size`` function returns the number of observations if
+        ``size`` function will return the number of observations when
-        ``dim=1`` or the number of variables if ``dim=2`` (for all
+        ``dim=1`` or the number of variables when ``dim=2``. An error
-        other values of ``dim`` an error is issued).
+        will be issued for any other values of ``dim``.
        *Example*
@@ -2819,18 +3122,20 @@ The dseries class
    .. dseriesmethod:: B = std (A[, geometric])
        |br| Overloads the MATLAB/Octave ``std`` function for
-        ``dseries`` objects. Returns the standard deviation of each
+        ``dseries`` objects. This function returns the standard
-        variable in ``dseries`` object ``A``. If the second argument
+        deviation of each variable within the ``dseries`` object
-        is ``true`` the geometric standard deviation is computed
+        ``A``. If the second argument is set to ``true``, the
-        (default value of the second argument is ``false``).
+        geometric standard deviation is calculated (the default value
+        for the second argument is ``false``).
    .. dseriesmethod:: B = subsample (A, d1, d2)
-        |br| Returns a subsample, for periods between ``dates`` ``d1``
+        |br| Returns a subsample for the period between ``d1`` and
-        and ``d2``. The same can be achieved by indexing a
+        ``d2``. While you can achieve the same result by indexing a
-        ``dseries`` object with a ``dates`` object, but the
+        ``dseries`` object with a ``dates`` object, the ``subsample``
-        ``subsample`` method is easier to use programmatically.
+        method offers a more straightforward approach for programmatic
+        use.
        *Example*
@@ -2849,7 +3154,7 @@ The dseries class
    .. dseriesmethod:: A = tag (A, a[, b, c])
-        |br| Add a tag to a variable in ``dseries`` object ``A``.
+        |br| Adds a tag to a variable in ``dseries`` object ``A``.
        *Example*
@@ -2864,23 +3169,24 @@ The dseries class
    .. dseriesmethod:: B = tex_rename (A, name, newtexname)
                       B = tex_rename (A, newtexname)
-                       tex_rename_ (A, name, newtexname)
+                       tex_rename_ (name, newtexname)
-                       tex_rename_ (A, newtexname)
+                       tex_rename_ (newtexname)
-        |br| Redefines the tex name of variable ``name`` to
+        |br| Updates the TeX name of the variable ``name`` to
-        ``newtexname`` in ``dseries`` object ``A``. Returns a
+        ``newtexname`` in the ``dseries`` object ``A``. Returns an
-        ``dseries`` object.
+        updated ``dseries`` object.
-        With only two arguments ``A`` and ``newtexname``, it redefines
+        With just two arguments, ``A`` and ``newtexname``, this
-        the tex names of the ``A`` to those contained in
+        function redefines the TeX names of the entries in ``A`` to
-        ``newtexname``. Here, ``newtexname`` is a cell string array
+        those specified in ``newtexname``. The ``newtexname`` argument
-        with the same number of entries as variables in ``A``.
+        must be a cell row character arrays containing the same number of
+        entries as there are variables in ``A``.
    .. dseriesmethod:: B = uminus(A)
-        |br| Overloads ``uminus`` (``-``, unary minus) for ``dseries``
+        |br| Overloads the ``uminus`` operator (``-``, unary minus)
-        object.
+        for the ``dseries`` object.
        *Example*
@@ -2903,13 +3209,13 @@ The dseries class
    .. dseriesmethod:: D = vertcat (A, B[, ...])
-        |br| Overloads the ``vertcat`` MATLAB/Octave method for
+        |br| Overloads the ``vertcat`` method in MATLAB/Octave for
-        ``dseries`` objects. This method is used to append more
+        ``dseries`` objects. This method facilitates the appending of
-        observations to a ``dseries`` object. Returns a ``dseries``
+        additional observations to a ``dseries`` object. It returns a
-        object ``D`` containing the variables in ``dseries`` objects
+        new ``dseries`` object, ``D``, which contains the variables
-        passed as inputs. All the input arguments must be ``dseries``
+        from the input ``dseries`` objects. All input arguments must
-        objects with the same variables defined on different time
+        be ``dseries`` objects that share the same variables but are
-        ranges.
+        defined over different time ranges.
        *Example*
@@ -2930,8 +3236,7 @@ The dseries class
    .. dseriesmethod:: B = vobs (A)
-        |br| Returns the number of variables in ``dseries`` object
+        |br| Returns the count of variables in the ``dseries`` object ``A``.
-        ``A``.
        *Example*
@@ -2947,10 +3252,10 @@ The dseries class
    .. dseriesmethod:: B = ydiff (A)
                       B = ygrowth (A)
-                       ydiff_ (A)
+                       ydiff_ ()
-                       ygrowth_ (A)
+                       ygrowth_ ()
-        |br| Computes yearly differences or growth rates.
+        |br| Calculates annual differences or growth rates.
 .. _x13-members:
@@ -3212,10 +3517,10 @@ X-13 ARIMA-SEATS interface
        The above example shows how to remove a seasonal pattern from a time series.
-        ``o.transform('function','auto','savelog','atr')`` instructs the subsequent 
+        ``o.transform('function','auto','savelog','atr')`` instructs the subsequent
-        ``o.automdl()`` command to check whether an additional or a multiplicative 
+        ``o.automdl()`` command to check whether an additional or a multiplicative
-        pattern fits the data better and to save the result. The result is saved in 
+        pattern fits the data better and to save the result. The result is saved in
-        `o.results.autotransform`, which in the present example indicates that a 
+        `o.results.autotransform`, which in the present example indicates that a
        log transformation, i.e. a multiplicative model was preferred. The ``o.automdl('savelog','all')`` automatically selects a fitting
        ARIMA model and saves all relevant output to the .log-file. The ``o.x11('save','(d11, d10)')`` instructs
        ``x11`` to save both the final seasonally adjusted series ``d11`` and the final seasonal factor ``d10``

--- a/doc/manual/utils/dynare_dom.py
+++ b/doc/manual/utils/dynare_dom.py
@@ -97,7 +97,7 @@ class DynObject(ObjectDescription):
                self.state_machine.reporter.warning(
                    'duplicate object description of %s, ' % fullname +
                    'other instance in ' +
-                    self.env.doc2path(objects[fullname][0]),line=self.lineno)
+                    str(self.env.doc2path(objects[fullname][0])),line=self.lineno)
            objects[fullname] = (self.env.docname, self.objtype)
        indextext = self.get_index_text(fullname,name_obj)
@@ -242,7 +242,7 @@ class DynSimpleObject(ObjectDescription):
                self.state_machine.reporter.warning(
                    'duplicate object description of %s, ' % fullname +
                    'other instance in ' +
-                    self.env.doc2path(objects[fullname][0]), line=self.lineno)
+                    str(self.env.doc2path(objects[fullname][0])), line=self.lineno)
            objects[fullname] = self.env.docname, self.objtype
        indextext = self.get_index_text(fullname,name_obj)

--- a/license.txt
+++ b/license.txt
 Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
 Upstream-Name: Dynare
-Upstream-Contact: Dynare Team, whose members in 2024 are:
+Upstream-Contact: Dynare Team, whose members in 2025 are:
                  - Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
                  - Michel Juillard <michel.juillard@mjui.fr>
                  - Frédéric Karamé <frederic.karame@univ-lemans.fr>
@@ -23,7 +23,7 @@ Upstream-Contact: Dynare Team, whose members in 2024 are:
 Source: https://www.dynare.org
 Files: *
-Copyright: 1996-2024 Dynare Team
+Copyright: 1996-2025 Dynare Team
 License: GPL-3+
 Files: matlab/+occbin/IVF_core.m
@@ -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/*

--- a/macOS/build.sh
+++ b/macOS/build.sh
@@ -40,11 +40,6 @@ else
 fi
 MATLAB_PATH=/Applications/"$PKG_ARCH"/MATLAB_R2024b.app
-# Workaround for bug in Xcode 15.3 which does not include m4
-# See https://github.com/Homebrew/homebrew-core/issues/165388
-# and https://trac.macports.org/ticket/69639
-path_prepend PATH "$BREWDIR"/opt/m4/bin
 # Append texbin to PATH to access latexmk and friends
 path_prepend PATH /Library/TeX/texbin

--- a/matlab/+gsa/run.m
+++ b/matlab/+gsa/run.m
@@ -360,7 +360,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
        if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
            fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
            fprintf('To obtain it, go to:\n\n')
-            fprintf('https://ec.europa.eu/jrc/en/macro-econometric-statistical-software/ss-anova-r-downloads \n\n')
+            fprintf('https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip \n\n')
            fprintf('and follow the instructions there.\n')
            fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
            error('SS-ANOVA-R Toolbox missing!')

--- a/matlab/+gsa/stability_mapping.m
+++ b/matlab/+gsa/stability_mapping.m
@@ -92,6 +92,10 @@ xparam1=[];
 [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
 if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
    % Set prior bounds
+    if options_.prior_trunc==0
+        fprintf('\nstability_mapping: GSA with priors requires bounded support. Setting options_.prior_trunc=1e-10.\n')
+        options_.prior_trunc=1e-10;
+    end
    bounds = prior_bounds(bayestopt_, options_.prior_trunc);
    bounds.lb = max(bounds.lb,lb);
    bounds.ub = min(bounds.ub,ub);

--- a/matlab/+hank/check_steady_state_input.m
+++ b/matlab/+hank/check_steady_state_input.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/>.
+%
+% Compute the stochastic simulations of heterogeneus-agent models
+%
+% INPUTS
+% - M_       [structure] Matlab's structure describing the model
+% - options_ [structure] Matlab's structure describing the current options
+% - ss       [structure] Matlab's structure with the model steady-state
+%                       information. It contains the following fields:
+%    - pol [structure] Matlab's structure containing the policy functions
+%                      discretization
+%       - pol.grids [structure]: Matlab's structure containing the nodes of the
+%                                state grids as column vectors.
+%       - pol.values [structure]: Matlab's structure containing the policy
+%                                 function as matrices. Row indices and column
+%                                 indices follow the lexicographic order
+%                                 specified in pol.shocks and pol.states
+%       - pol.shocks [array]: a list containing the order of shock/innovation
+%                             variables used for the rows of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement (resp. varexo(heterogeneity=)
+%                             statement) for discretizes AR(1) processes (resp.
+%                             discretized gaussian i.i.d innovations).
+%       - pol.states [array]: a list containing the order of state variables
+%                             used for the columns of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement.
+%    - shocks [structure]: Matlab's stucture describing the discretization of
+%                          individual shocks or innovations:
+%       - shocks.grids [structure]: Matlab's structure containing the nodes of
+%                                   the shock grids
+%       - shocks.Pi [structure]: Matlab's structure containing the Markov
+%                                matrices if the shock processes are discretized
+%                                AR(1) processes. The field should be absent
+%                                otherwise.
+%       - shocks.w [structure]: Matlab's structure containing the Gauss-Hermite
+%                               weights if the i.i.d gaussian innovation
+%                               processes are discretized
+%    - d [structure]: Matlab's structure describing the steady-state
+%                     distribution
+%       - d.grids [structure]: structure containing the states grids as column
+%                              vectors. If one of the states grid is not
+%                              specified, it is assumed to be identical to
+%                              its policy-function value stored in pol.grids
+%       - d.hist [matrix]: the histogram of the distribution as a matrix with
+%                          shocks/innovation indices as rows and state indices
+%                          as columns. Row and column indices follow the
+%                          lexicographic order specified in d.shocks and
+%                          d.states
+%       - d.shocks [array]: a list containing the order of shock/innovation
+%                           variables used for the rows of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.shocks.
+%       - d.states [array]: a list containing the order of state variables used
+%                           for the columns of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.states.
+%    - agg [structure]: Matlab's tructure containing the steady-state values of
+%                       aggregate variables
+% OUTPUTS
+% - out_ss [structure]: validated and normalized steady-state input structure.
+%                       It has a similar structure as the ss structure in inputs.
+% - sizes [structure]: structure containing sizes information
+%    - n_e [integer]: number of shocks
+%    - N_e [integer]: total number of nodes in the shocks grids
+%    - shocks [structure]: number of nodes for each shock grid
+%    - n_a [integer]: number of states
+%    - pol [structure]: 
+%       - pol.N_a [integer]: total number of nodes in the policy state grids
+%       - pol.states [structure]: number of nodes for each state grid
+%    - d [structure]:
+%       - d.states [structure]: number of nodes for each distribution state grid
+%       - d.N_a [integer]: total number of nodes in the distribution grids
+%    - n_pol [integer]: number of policy variables
+%    - agg [integer]: number of aggregate variables
+function [out_ss, sizes] = check_steady_state_input(M_, options_, ss)
+   %% Checks
+   % Retrieve variables information from M_
+   state_symbs = M_.heterogeneity(1).endo_names(M_.heterogeneity(1).state_var);
+   pol_symbs= M_.heterogeneity(1).endo_names;
+   inn_symbs = M_.heterogeneity(1).exo_names;
+   agg_symbs = M_.endo_names(1:M_.orig_endo_nbr);
+   % Initialize output variables
+   sizes = struct;
+   out_ss = struct;
+   if ~isstruct(ss)
+      error('Misspecified steady-state input `ss`: the steady-state input `ss` is not a structure.')
+   end
+   %% Shocks
+   % Check that the field `ss.shocks` exists
+   check_isfield('shocks', ss, 'ss.shocks', ' Models without idiosyncratic shocks but with ex-post heterogeneity over individual state variables cannot be solved yet.')
+   shocks = ss.shocks;
+   % Check that `ss.shocks` is a structure
+   check_isstruct(shocks, 'ss.shocks');
+   % Check that the field `ss.shocks.grids` exists
+   check_isfield('grids', ss.shocks, 'ss.shocks.grids');
+   % Check that `ss.shocks.grids` is a structure
+   check_isstruct(shocks.grids, 'ss.shocks.grids');
+   shock_symbs = fieldnames(shocks.grids);
+   % Check the types of `ss.shocks.grids` values
+   check_fun(shocks.grids, 'ss.shocks.grids', shock_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(shock_symbs)
+      s = shock_symbs{i};
+      if isrow(ss.shocks.grids.(s))
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s)';
+      else
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s);
+      end
+   end
+   % Check shock discretization
+   flag_ar1 = isfield(shocks, 'Pi');
+   flag_gh = isfield(shocks, 'w');
+   if ~flag_ar1 && ~flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both non-specified, while exactly one of them should be set.');
+   end
+   % Check that either individual AR(1) processes or i.i.d gaussion innovation
+   % processes are discretized
+   if flag_ar1 && flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both specified, while only one of them should be.');
+   end
+   % Case of discretized AR(1) exogenous processes
+   if flag_ar1
+      % Check that the grids specification for AR(1) shocks and the
+      % var(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, pol_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).endo_names');
+      % Check that shocks.Pi is a structure
+      check_isstruct(shocks.Pi, 'ss.shocks.Pi');
+      % Check that the grids specification for individual AR(1) shocks, the
+      % information in the M_ structure and the Markov transition matrices
+      % specification are compatible and display a warning if redundancies are
+      % detected
+      check_missingredundant(shocks.Pi, 'ss.shocks.Pi', shock_symbs, options_.hank.nowarningredundant);
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the type of shocks.Pi field values
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) ismatrix(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real matrices');
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      % Check that the matrices provided in shocks.Pi are Markov matrices:
+      %     - all elements are non-negative
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(x >= 0, 'all'), 'contain negative elements')
+      %     - the sum of each row equals 1.
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(abs(sum(x,2)-1) < options_.hank.tol_check_sum), 'have row sums different from 1.')
+      % Remove AR(1) shock processes from state and policy variables
+      state_symbs = setdiff(state_symbs, shock_symbs);
+      pol_symbs = setdiff(pol_symbs, shock_symbs);
+      % Copy relevant shock-related elements in out_ss 
+      out_ss.shocks.Pi = ss.shocks.Pi;
+   end
+   % Case of discretized i.i.d gaussian innovations
+   if flag_gh
+      % Check that shocks.w is a structure
+      check_isstruct(shocks.w, 'shocks.w');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the varexo(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, inn_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).exo_names');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the Gauss-Hermite weights specification are compatible and display a
+      % warning if redundancies are detected
+      check_missingredundant(shocks.w, 'ss.shocks.w', shock_symbs, options_.hank.nowarningredundant);
+      % Check the type of `ss.shocks.w` elements
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) isvector(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real vectors');
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids 
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) numel(x), 'have incompatible sizes with those of `ss.shocks.grids`', grid_nb);
+      % Check that the arrays provided in shocks.w have the properties of
+      % Gauss-Hermite weights:
+      %     - all elements are non-negative
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(x >= 0.), 'contain negative elements');
+      %     - the sum of Gauss-Hermite weights is 1
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(abs(sum(x)-1) < options_.hank.tol_check_sum), 'have sums different from 1.');
+      % Make `ss.shocks.w` values column vectors in the output steady-state
+      % structure
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         if isrow(ss.shocks.w.(s))
+            out_ss.shocks.w.(s) = ss.shocks.w.(s)';
+         else
+            out_ss.shocks.w.(s) = ss.shocks.w.(s);
+         end
+      end
+   end
+   %% Policy functions
+   % Check that the `ss.pol` field exists
+   check_isfield('pol', ss, 'ss.pol');
+   pol = ss.pol;
+   % Check that `ss.pol` is a structure
+   check_isstruct(ss.pol, 'ss.pol');
+   % Check that the `ss.pol`.grids field exists
+   check_isfield('grids', ss.pol, 'ss.pol.grids');
+   % Check that `ss.pol.grids` is a structure
+   check_isstruct(ss.pol.grids, 'ss.pol.grids');
+   % Check that the state grids specification and the var(heterogeneity=)
+   % statement are compatible
+   check_missingredundant(pol.grids, 'ss.pol.grids', state_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.grids` values are dense real vectors
+   check_fun(pol.grids, 'ss.pol.grids', state_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Store:
+   %    - the number of states
+   sizes.n_a = numel(state_symbs);
+   %    - the size of the tensor product of the states grids 
+   grid_nb = cellfun(@(s) numel(pol.grids.(s)), state_symbs);
+   sizes.pol.N_a = prod(grid_nb);
+   %    - the size of each state grid
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      sizes.pol.states.(s) = numel(pol.grids.(s));
+   end
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      if isrow(ss.pol.grids.(s))
+         out_ss.pol.grids.(s) = ss.pol.grids.(s)';
+      else
+         out_ss.pol.grids.(s) = ss.pol.grids.(s);
+      end
+   end
+   % Check that the field `ss.pol.values` exists 
+   check_isfield('values', pol, 'ss.pol.values');
+   % Check that `ss.pol.values` is a struct
+   check_isstruct(pol.values, 'ss.pol.values');
+   % Check the missing and redundant variables in `ss.pol.values`
+   check_missingredundant(pol.values, 'ss.pol.values', pol_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.values` values are dense real matrices
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) isnumeric(x) && isreal(x) && ismatrix(x) && ~issparse(x), 'are not dense real matrices'); 
+   % Check the internal size compatibility of `ss.pol.values`
+   sizes.n_pol = numel(pol_symbs);
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the sizes of `ss.shocks.grids` elements', sizes.N_e);
+   check_fun(pol.values, 'ss.pol.values', pol_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the sizes of `ss.pol.grids` elements', sizes.pol.N_a);
+   % Copy `ss.pol.values` in `out_ss`
+   out_ss.pol.values = ss.pol.values;
+   % Check the permutation of state variables for policy functions
+   [out_ss.pol.states] = check_permutation(pol, 'states', 'ss.pol', state_symbs);
+   % Check the permutation of shock variables for policy functions
+   [out_ss.pol.shocks] = check_permutation(pol, 'shocks', 'ss.pol', shock_symbs);
+   %% Distribution
+   % Check that the field `ss.d` exists 
+   check_isfield('d', ss, 'ss.d');
+   d = ss.d;
+   % Check that the field `ss.d.hist` exists 
+   check_isfield('hist', ss.d, 'ss.d.hist');
+   % Check the type of `ss.d.hist`
+   if ~(ismatrix(d.hist) && isnumeric(d.hist) && isreal(d.hist) && ~issparse(d.hist))
+      error('Misspecified steady-state input `ss`: `ss.d.hist` is not a dense real matrix.');
+   end
+   % Check the consistency of `ss.d.grids` 
+   if ~isfield(d, 'grids')
+      if ~options_.hank.nowarningdgrids
+         warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+      end
+      % Copy the relevant sizes from the pol field
+      sizes.d = sizes.pol;
+      out_ss.d.grids = out_ss.pol.grids;
+   else
+      % Check that `ss.d.grids` is a struct
+      check_isstruct(d.grids, 'ss.d.grids');
+      % Check redundant variables in `ss.d.grids`
+      d_grids_symbs = fieldnames(d.grids);
+      if isempty(d_grids_symbs)
+         if ~options_.hank.nowarningredundant
+            warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+         end
+         % Copy the relevant sizes from the pol field
+         sizes.d = sizes.pol;
+         out_ss.d.grids = out_ss.pol.grids;
+      else
+         d_grids_symbs_in_states = ismember(d_grids_symbs, state_symbs);
+         if ~all(d_grids_symbs_in_states)
+            if ~options_.hank.nowarningredundant
+               warning('In the steady-state input `ss.d.states`, the following specification for the states grids in the distribution structure are not useful: %s', strjoin(d_grids_symbs(~d_grids_symbs_in_states)));
+            end
+            d_grids_symbs = d_grids_symbs(d_grids_symbs_in_states);
+         end
+         % Check the types of `ss.d.grids` elements
+         check_fun(pol.grids, 'ss.d.grids', d_grids_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+         % Store the size of the states grids for the distribution
+         for i=1:numel(d_grids_symbs)
+            s = d_grids_symbs{i};
+            sizes.d.states.(s) = numel(d.grids.(s));
+            out_ss.d.grids.(s) = d.grids.(s);
+         end
+         states_out_of_d = setdiff(state_symbs, d_grids_symbs);
+         if ~isempty(states_out_of_d)
+            if ~options_.hank.nowarningdgrids
+               warning('hank.bbeg.het_stoch_simul: in the steady-state structure, no distribution-specific grid set for states %s. The policy grids specified in pol.grids shall be used.', strjoin(states_out_of_d));
+            end
+            % Store the sizes of the unspecified states grids from the pol field
+            for i=1:numel(states_out_of_d)
+               s = states_out_of_d{i};
+               sizes.d.states.(s) = sizes.pol.states.(s);
+               out_ss.d.grids.(s) = pol.grids.(s);
+            end
+         end
+      end
+   end
+   % Check the internal size compatibility of the distribution histogram
+   if size(d.hist,1) ~= sizes.N_e
+      error('Misspecified steady-state input `ss`: the number of rows of the histogram matrix `ss.d.hist` is not consistent with the sizes of shocks grids `ss.shocks.grids`');
+   end
+   sizes.d.N_a = prod(structfun(@(x) x, sizes.d.states));
+   if size(d.hist,2) ~= sizes.d.N_a
+      error('Misspecified steady-state input `ss`: the number of columns of the histogram matrix `ss.d.hist` is not consistent with the sizes of states grids `ss.d.grids`/`ss.pol.grids`');
+   end
+   % Copy `ss.d.hist` in `out_ss` 
+   out_ss.d.hist = ss.d.hist;
+   % Check the permutation of state variables in the distribution
+   [out_ss.d.states] = check_permutation(d, 'states', 'ss.d', state_symbs);
+   % Check the permutation of shock variables in the distribution
+   [out_ss.d.shocks] = check_permutation(d, 'shocks', 'ss.d', shock_symbs);
+   %% Aggregate variables
+   % Check that the field `ss.agg` exists
+   check_isfield('agg', ss, 'ss.agg');
+   agg = ss.agg;
+   % Check that `ss.agg` is a structure
+   check_isstruct(agg, 'ss.agg');
+   % Check that the aggregate variables specification and the var statement are
+   % compatible
+   check_missingredundant(agg, 'ss.agg', agg_symbs, options_.hank.nowarningredundant);
+   % Store the number of aggregate variables
+   sizes.agg = numel(fieldnames(agg));
+   % Check the types of `ss.agg` values
+   check_fun(agg, 'ss.agg', agg_symbs, @(x) isreal(x) && isscalar(x), 'are not real scalars');
+   % Copy `ss.agg` into `out_ss`
+   out_ss.agg = agg;
+end
--- a/matlab/+hank/private/check_consistency.m
+++ b/matlab/+hank/private/check_consistency.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 consistency between two sets of variable names.
+%
+% INPUTS
+% - f1 [cell array of strings]: list of variable names to be checked for consistency
+% - f2 [cell array of strings]: reference list of valid variable names
+% - f1_name [char]: name of the first set (for error messages)
+% - f2_name [char]: name of the second set (for error messages)
+%
+% OUTPUTS
+% - (none) This function throws an error if there are elements in `f1` that are not present in `f2`.
+%
+% DESCRIPTION
+% Checks that all elements of `f1` are included in `f2`. If not, throws a descriptive error 
+% mentioning the missing variables and the corresponding structure names for easier debugging.
+function [] = check_consistency(f1, f2, f1_name, f2_name)
+   f1_in_f2 = ismember(f1,f2);
+   if ~all(f1_in_f2)
+      error('Misspecified steady-state input `ss`: the following variables are mentioned in `%s`, but are missing in `%s`: %s', f1_name, f2_name, strjoin(f1(~f1_in_f2)));
+   end
+end
\ No newline at end of file
--- a/matlab/+hank/private/check_fun.m
+++ b/matlab/+hank/private/check_fun.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/>.
+%
+% Apply a validation function to a set of fields within a structure.
+%
+% INPUTS
+% - f [structure]: structure containing fields to be validated
+% - f_name [char]: name of the structure (for error messages)
+% - symbs [cell array of strings]: list of field names to validate
+% - fun [function handle]: validation function applied to each field
+% - text [char]: descriptive text for the error message in case of validation failure
+% - eq [optional, scalar]: if provided, the output of `fun` is compared to `eq`
+%
+% OUTPUTS
+% - (none) This function throws an error if any field fails the validation.
+%
+% DESCRIPTION
+% For each symbol in `symbs`, applies the validation function `fun` to the
+% corresponding field in `f`. If an optional `eq` argument is provided, the
+% output of `fun` is compared to `eq`. Throws a descriptive error listing all
+% fields that do not satisfy the validation criteria.
+function [] = check_fun(f, f_name, symbs, fun, text, eq)
+   elt_check = cellfun(@(s) fun(f.(s)), symbs);
+   if nargin == 6
+      elt_check = elt_check == eq;
+   end
+   if ~all(elt_check)
+      error('Misspecified steady-state input `ss`: the following fields in `%s` %s: %s', f_name, text, strjoin(symbs(~elt_check)));
+   end
+end
\ No newline at end of file
--- a/matlab/+hank/private/check_isfield.m
+++ b/matlab/+hank/private/check_isfield.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 the existence of a field in a structure.
+%
+% INPUTS
+% - f [char]: name of the field to check
+% - s [structure]: structure in which the field should be present
+% - f_name [optional, char]: name to display in the error message (defaults to `f`)
+% - details [optional, char]: additional details to append to the error message
+%
+% OUTPUTS
+% - (none) This function throws an error if the specified field is missing.
+%
+% DESCRIPTION
+% Checks whether the field `f` exists in the structure `s`.
+% If not, throws an error indicating the missing field.
+% If provided, `details` is appended to the error message for additional context.
+function [] = check_isfield(f, s, f_name, details)
+   if nargin < 4
+      details = '';
+   end
+   if nargin < 3
+      f_name = f;
+   end
+   if ~isfield(s, f)
+      error('Misspecified steady-state input `ss`: the `%s` field is missing.%s', f_name, details);
+   end
+end
\ No newline at end of file
No results found