Bug fix (wrong condition)

We should test the numbe of output arguments, not the number of input arguments. This was bug was probably not affecting the outcomes since the number of input arguments is always greater than 1.
Cosmetic changes.
2024-01-22 16:54:02 +01:00 · 2024-01-22 16:54:01 +01:00 · 2024-01-22 16:47:32 +01:00 · 2024-01-22 16:47:23 +01:00 · 2024-01-22 16:00:39 +01:00 · 2024-01-19 20:40:56 +01:00
55 changed files with 829 additions and 948 deletions
--- a/doc/manual/source/conf.py
+++ b/doc/manual/source/conf.py
@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-

-# Copyright © 2018-2023 Dynare Team
+# Copyright © 2018-2024 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–2023 Dynare Team'
+copyright = u'1996–2024 Dynare Team'
 author = u'Dynare Team'

 add_function_parentheses = False
@ -71,7 +71,6 @@ latex_elements = {
                    warningBorderColor={RGB}{255,50,50},OuterLinkColor={RGB}{34,139,34}, \
                    InnerLinkColor={RGB}{51,51,255},TitleColor={RGB}{51,51,255}',
    'papersize': 'a4paper',
-    'preamble': r'\DeclareUnicodeCharacter{200B}{}', # Part of the workaround for #1707
 }

 latex_documents = [
--- a/doc/manual/source/dynare-misc-commands.rst
+++ b/doc/manual/source/dynare-misc-commands.rst
@ -8,17 +8,17 @@
 Dynare misc commands
 ####################

-.. matcomm:: send_endogenous_variables_to_workspace
+.. matcomm:: send_endogenous_variables_to_workspace ;

    Puts the simulation results for the endogenous variables stored in ``oo_.endo_simul`` 
    into vectors with the same name as the respective variables into the base workspace.

-.. matcomm:: send_exogenous_variables_to_workspace
+.. matcomm:: send_exogenous_variables_to_workspace ;

    Puts the simulation results for the exogenous variables stored in ``oo_.exo_simul`` 
    into vectors with the same name as the respective variables into the base workspace.

-.. matcomm:: send_irfs_to_workspace
+.. matcomm:: send_irfs_to_workspace ;

    Puts the IRFs stored in ``oo_.irfs`` into vectors with the same name into the base workspace.

--- a/doc/manual/source/index.rst
+++ b/doc/manual/source/index.rst
@ -11,7 +11,7 @@ Currently the development team of Dynare is composed of:
 * Willi Mutschler (University of Tübingen)
 * Johannes Pfeifer (University of the Bundeswehr Munich)
 * Marco Ratto (European Commission, Joint Research Centre - JRC)
-* Normann Rion (CY Cergy Paris Université and CEPREMAP)
+* Normann Rion (CEPREMAP)
 * Sébastien Villemot (CEPREMAP)

 The following people used to be members of the team:
@ -26,7 +26,7 @@ The following people used to be members of the team:
 * Ferhat Mihoubi
 * George Perendia

-Copyright © 1996-2023, Dynare Team.
+Copyright © 1996-2024, 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/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@ -1259,7 +1259,8 @@ command, the list of transformed model equations using the
 ``write_latex_dynamic_model command``, and the list of static model
 equations using the ``write_latex_static_model`` command.

-.. command:: write_latex_original_model (OPTIONS);
+.. command:: write_latex_original_model ;
+             write_latex_original_model (OPTIONS);

    |br| This command creates two LaTeX files: one
    containing the model as defined in the model block and one
@ -1334,7 +1335,8 @@ equations using the ``write_latex_static_model`` command.
        See :opt:`write_equation_tags`


-.. command:: write_latex_static_model (OPTIONS);
+.. command:: write_latex_static_model ;
+             write_latex_static_model (OPTIONS);

    |br| This command creates two LaTeX files: one
    containing the static model and one containing the LaTeX
@ -1369,7 +1371,7 @@ equations using the ``write_latex_static_model`` command.

        See :opt:`write_equation_tags`.

-.. command:: write_latex_steady_state_model
+.. command:: write_latex_steady_state_model ;

    |br| This command creates two LaTeX files: one containing the steady
    state model and one containing the LaTeX document header
@ -2945,8 +2947,8 @@ Finding the steady state with Dynare nonlinear solver

           ``5``

-                Newton algorithm with a sparse Gaussian elimination
-                (SPE) (requires ``bytecode`` option, see
+                Newton algorithm with a sparse Gaussian elimination (SPE)
+                solver at each iteration (requires ``bytecode`` option, see
                :ref:`model-decl`).

           ``6``
@ -2964,7 +2966,7 @@ Finding the steady state with Dynare nonlinear solver
           ``8``

                Newton algorithm with a Stabilized Bi-Conjugate
-                Gradient (BICGSTAB) solver at each iteration (requires
+                Gradient (BiCGStab) solver at each iteration (requires
                bytecode and/or block option, see :ref:`model-decl`).

           ``9``
@ -3680,60 +3682,64 @@ speed-up on large models.
           ``0``

               Use a Newton algorithm with a direct sparse LU solver at each
-               iteration, applied on the stacked system of all the equations at
-               every period (Default).
+               iteration, applied to the stacked system of all equations in all
+               periods (Default).

           ``1``

               Use the Laffargue-Boucekkine-Juillard (LBJ) algorithm proposed
-               in *Juillard (1996)*. It is slower than ``stack_solve_algo=0``,
-               but may be less memory consuming on big models. Note that if the
-               ``block`` option is used (see :ref:`model-decl`), a simple
-               Newton algorithm with sparse matrices is used for blocks which
-               are purely backward or forward (of type ``SOLVE BACKWARD`` or
-               ``SOLVE FORWARD``, see :comm:`model_info`), since LBJ only makes
-               sense on blocks with both leads and lags (of type ``SOLVE TWO
-               BOUNDARIES``).
+               in *Juillard (1996)* on top of a LU solver. It is slower
+               than ``stack_solve_algo=0``, but may be less memory consuming on
+               big models. Note that if the ``block`` option is used (see
+               :ref:`model-decl`), a simple Newton algorithm with sparse
+               matrices, applied to the stacked system of all block equations
+               in all periods, is used for blocks which are purely backward or
+               forward (of type ``SOLVE BACKWARD`` or ``SOLVE FORWARD``, see
+               :comm:`model_info`), since LBJ only makes sense on blocks with
+               both leads and lags (of type ``SOLVE TWO BOUNDARIES``).

           ``2``

-               Use a Newton algorithm with a Generalized Minimal
-               Residual (GMRES) solver at each iteration (requires
-               ``bytecode`` and/or ``block`` option, see
-               :ref:`model-decl`)
+               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 (requires ``bytecode`` and/or
+               ``block`` option, see :ref:`model-decl`)

           ``3``

-               Use a Newton algorithm with a Stabilized Bi-Conjugate
-               Gradient (BICGSTAB) solver at each iteration (requires
-               ``bytecode`` and/or ``block`` option, see
-               :ref:`model-decl`).
+               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 (requires ``bytecode``
+               and/or ``block`` option, see :ref:`model-decl`).

           ``4``

-               Use a Newton algorithm with an optimal path length at
-               each iteration (requires ``bytecode`` and/or ``block``
-               option, see :ref:`model-decl`).
+               Use a Newton algorithm with a direct sparse LU solver and an
+               optimal path length at each iteration, applied on the stacked
+               system of all equations in all periods (requires ``bytecode``
+               and/or ``block`` option, see :ref:`model-decl`).

           ``5``

-               Use a Newton algorithm with a sparse Gaussian
-               elimination (SPE) solver at each iteration (requires
-               ``bytecode`` option, see :ref:`model-decl`).
+               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
+               (requires ``bytecode`` option, see :ref:`model-decl`).

           ``6``

-               Synonymous for ``stack_solve_algo=1``. Kept for historical
-               reasons.
+               Synonymous for ``stack_solve_algo=1``. Kept for backward
+               compatibility.

           ``7``

-               Allows the user to solve the perfect foresight model
-               with the solvers available through option
-               ``solve_algo`` (See :ref:`solve_algo <solvalg>` for a
-               list of possible values, note that values 5, 6, 7 and
-               8, which require ``bytecode`` and/or ``block`` options,
-               are not allowed). For instance, the following
+               Allows the user to solve the perfect foresight model with the
+               solvers available through option ``solve_algo``, applied on the
+               stacked system of all equations in all periods (See
+               :ref:`solve_algo <solvalg>` for a list of possible values, note
+               that values 5, 6, 7 and 8, which require ``bytecode`` and/or
+               ``block`` options, are not allowed). For instance, the following
               commands::

                    perfect_foresight_setup(periods=400);
@ -11739,7 +11745,7 @@ with ``discretionary_policy`` or for optimal simple rules with ``osr``

    With ``discretionary_policy``, the objective function must be quadratic.

-.. command:: evaluate_planner_objective;
+.. command:: evaluate_planner_objective ;
             evaluate_planner_objective (OPTIONS...);

    This command computes, displays, and stores the value of the
@ -14293,7 +14299,7 @@ a trend target to which the endogenous variables may be attracted in the long-ru
   :math:`n\times 1` vector of parameters, :math:`A_i` (:math:`i=0,\ldots,p`)
   are :math:`n\times n` matrices of parameters, and :math:`A_0` is non
   singular square matrix. Vector :math:`\mathbf{c}` and matrices :math:`A_i`
-   (:math:`i=0,\ldots,p`) are set by Dynare by parsing the equations in the
+   (:math:`i=0,\ldots,p`) are set by parsing the equations in the
   ``model`` block. Then, Dynare builds a VAR(1)-companion form model for
   :math:`\mathcal{Y}_t = (1, Y_t, \ldots, Y_{t-p+1})'` as:

@ -14504,7 +14510,7 @@ up to time :math:`t-\tau`, :math:`\mathcal{Y}_{\underline{t-\tau}}`) is:

 In a semi-structural model, variables appearing in :math:`t+h` (*e.g.*
 the expected output gap in a dynamic IS curve or expected inflation in a
-(New Keynesian) Phillips curve) will be replaced by the expectation implied by an auxiliary VAR
+New Keynesian Phillips curve) will be replaced by the expectation implied by an auxiliary VAR
 model. Another use case is for the computation of permanent
 incomes. Typically, consumption will depend on something like:

@ -14512,13 +14518,13 @@ incomes. Typically, consumption will depend on something like:

      \sum_{h=0}^{\infty} \beta^h y_{t+h|t-\tau}

-Assuming that $0<\beta<1$ and knowing the limit of geometric series, the conditional expectation of this variable can be evaluated based on the same auxiliary model:
+Assuming that :math:`0<\beta<1` and knowing the limit of geometric series, the conditional expectation of this variable can be evaluated based on the same auxiliary model:

   .. math ::

      \mathbb E \left[\sum_{h=0}^{\infty} \beta^h y_{t+h}\Biggl| \mathcal{Y}_{\underline{t-\tau}}\right] = \alpha \mathcal{C}^\tau(I-\beta\mathcal{C})^{-1}\mathcal{Y}_{t-\tau}

-More generally, it is possible to consider finite discounted sums.
+Finite discounted sums can also be considered.

 .. command:: var_expectation_model (OPTIONS...);

@ -15861,7 +15867,7 @@ Misc commands
    ``pdflatex`` and automatically tries to load all required
    packages. Requires the following LaTeX packages:
    ``breqn``, ``psfrag``, ``graphicx``, ``epstopdf``, ``longtable``,
-    ``booktabs``, ``caption``, ``float,`` ``amsmath``, ``amsfonts``,
+    ``booktabs``, ``caption``, ``float,`` ``amsmath``, ``amsfonts``, ``amssymb``, 
    and ``morefloats``.


--- a/doc/manual/utils/dynare_dom.py
+++ b/doc/manual/utils/dynare_dom.py
@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-

-# Copyright © 2018-2019 Dynare Team
+# Copyright © 2018-2024 Dynare Team
 #
 # This file is part of Dynare.
 #
@ -80,9 +80,7 @@ class DynObject(ObjectDescription):
        signode += addnodes.desc_name(name, name)

        if self.has_arguments:
-            if not arglist:
-                signode += addnodes.desc_parameterlist()
-            else:
+            if arglist:
                signode += addnodes.desc_addname(arglist,arglist)
        return fullname, prefix

--- a/license.txt
+++ b/license.txt
@ -1,6 +1,6 @@
 Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
 Upstream-Name: Dynare
-Upstream-Contact: Dynare Team, whose members in 2023 are:
+Upstream-Contact: Dynare Team, whose members in 2024 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 2023 are:
 Source: https://www.dynare.org

 Files: *
-Copyright: 1996-2023 Dynare Team
+Copyright: 1996-2024 Dynare Team
 License: GPL-3+

 Files: matlab/+occbin/IVF_core.m
@ -86,7 +86,7 @@ License: public-domain-aim
 Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
 pages 472-489

-Files: matlab/optimization/bfgsi1.m matlab/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
+Files: matlab/optimization/bfgsi1.m matlab/optimization/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m
 matlab/optimization/numgrad3.m matlab/optimization/numgrad3_.m matlab/optimization/numgrad5.m
 matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/+bvar/density.m
 matlab/+bvar/toolbox.m matlab/partial_information/PI_gensys.m matlab/partial_information/qzswitch.m
@ -123,6 +123,11 @@ Copyright: 2010-2015 Alexander Meyer-Gohde
           2015-2017 Dynare Team
 License: GPL-3+

+Files: matlab/collapse_figures_in_tabgroup.m
+Copyright: 2023 Eduard Benet Cerda
+           2024 Dynare Team
+License: GPL-3+
+
 Files: matlab/convergence_diagnostics/raftery_lewis.m
 Copyright: 2016 Benjamin Born and Johannes Pfeifer
           2016-2017 Dynare Team
@ -172,7 +177,7 @@ Comment: Written by Jessica Cariboni and Francesca Campolongo
 Files: matlab/+gsa/cumplot.m
       matlab/+gsa/monte_carlo_filtering.m
       matlab/+gsa/skewness.m
-       matlab/+gsa/log_trans_.m
+       matlab/+gsa/log_transform.m
       matlab/+gsa/map_calibration.m
       matlab/+gsa/map_identification.m
       matlab/+gsa/monte_carlo_filtering_analysis.m
@ -247,7 +252,7 @@ License: BSD-2-clause

 Files: examples/fs2000_data.m
 Copyright: 2000-2022 Frank Schorfheide
-Copyright: 2023 Dynare Team
+           2023 Dynare Team
 License: CC-BY-SA-4.0

 Files: doc/*.rst doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
@ -292,28 +297,6 @@ Files: preprocessor/doc/preprocessor/*
 Copyright: 2007-2019 Dynare Team
 License: CC-BY-SA-4.0

-Files: contrib/jsonlab/*
-Copyright: 2011-2020 Qianqian Fang <q.fang at neu.edu>
-           2016 Bastian Bechtold
-License: GPL-3+ or BSD-3-clause
-
-Files: contrib/jsonlab/base64decode.m
-       contrib/jsonlab/base64encode.m
-       contrib/jsonlab/gzipdecode.m
-       contrib/jsonlab/gzipencode.m
-       contrib/jsonlab/zlibdecode.m
-       contrib/jsonlab/zlibencode.m
-Copyright: 2012 Kota Yamaguchi
-           2011-2020 Qianqian Fang <q.fang at neu.edu>
-License: GPL-3+ or BSD-2-clause
-
-Files: contrib/jsonlab/loadjson.m
-Copyright: 2011-2020 Qianqian Fang
-           2009 Nedialko Krouchev
-           2009 François Glineur
-           2008 Joel Feenstra
-License: GPL-3+ or BSD-2-clause or BSD-3-clause
-
 Files: contrib/ms-sbvar/utilities_dw/*
 Copyright: 1996-2011 Daniel Waggoner
 License: GPL-3+
@ -420,32 +403,6 @@ License: BSD-2-clause
 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.

-License: BSD-3-clause
- Redistribution and use in source and binary forms, with or without modification,
- are permitted provided that the following conditions are met:
- .
- * Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
- .
- * Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
- .
- * Neither the name of the copyright holder nor the names of its contributors
-   may be used to endorse or promote products derived from this software without
-   specific prior written permission.
- .
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
 License: GFDL-NIV-1.3+
 Permission is granted to copy, distribute and/or modify this document
 under the terms of the GNU Free Documentation License, Version 1.3 or
--- a/matlab/+bvar/forecast.m
+++ b/matlab/+bvar/forecast.m
@ -119,7 +119,7 @@ OutputDirectoryName = CheckPath('graphs',M_.dname);
 dyn_graph=bvar.graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});

 for i = 1:ny
-    dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
+    dyn_graph=plot_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
                        sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
                        sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
                           options_.varobs{i});
@ -183,3 +183,31 @@ for i = 1:length(options_.varobs)
        oo_.bvar.forecast.rmse.(name) = rmse(i);
    end
 end
+
+
+function dyn_graph=plot_graph(dyn_graph,y,tit,x)
+% function plot_graph(dyn_graph, y,tit,x)
+
+if nargin < 4
+    x = (1:size(y,1))';
+end
+nplot = dyn_graph.plot_nbr + 1;
+if nplot > dyn_graph.max_nplot
+    figure('Name',dyn_graph.figure_name);
+    nplot = 1;
+end
+dyn_graph.plot_nbr = nplot;
+subplot(dyn_graph.nr,dyn_graph.nc,nplot);
+
+line_types = dyn_graph.line_types;
+line_type = line_types{1};
+for i=1:size(y,2)
+    if length(line_types) > 1
+        line_type = line_types{i};
+    end
+
+    plot(x,y(:,i),line_type);
+    hold on
+end
+title(tit);
+hold off
--- a/matlab/+gsa/map_identification.m
+++ b/matlab/+gsa/map_identification.m
@ -272,80 +272,8 @@ elseif opt_gsa.morris==3
    return
 elseif opt_gsa.morris==2   % ISKREV stuff
    return
-else  % main effects analysis
-    if itrans==0
-        fsuffix = '';
-    elseif itrans==1
-        fsuffix = '_log';
-    else
-        fsuffix = '_rank';
-    end
-
-    imap=1:npT;
-
-    if isempty(lpmat0)
-        x0=lpmat(istable,:);
-    else
-        x0=[lpmat0(istable,:), lpmat(istable,:)];
-    end
-    nrun=length(istable);
-    nest=min(250,nrun);
-
-    if opt_gsa.load_ident_files==0
-        try
-            EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAcc','ir_cc','ic_cc');
-        catch
-            EET=[];
-        end
-        ccac = gsa.standardize_columns([mss cc ac]);
-        [pcc, dd] = eig(cov(ccac(istable,:)));
-        [latent, isort] = sort(-diag(dd));
-        latent = -latent;
-        figure, bar(latent)
-        title('Eigenvalues in PCA')
-        pcc=pcc(:,isort);
-        ccac = ccac*pcc;
-        npca = max(find(cumsum(latent)./length(latent)<0.99))+1;
-        siPCA = (EET.SAcc'*abs(pcc'))';
-        siPCA = siPCA./(max(siPCA,[],2)*ones(1,npT));
-        SAcc=zeros(size(ccac,2),npT);
-        gsa_=NaN(npca);
-        for j=1:npca %size(ccac,2),
-            if itrans==0
-                y0 = ccac(istable,j);
-            elseif itrans==1
-                y0 = gsa.log_transform(ccac(istable,j));
-            else
-                y0 = trank(ccac(istable,j));
-            end
-            if ~isempty(EET)
-                imap=find(siPCA(j,:) >= (0.1.*max(siPCA(j,:))) );
-            end
-            gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
-                              2, [],[],[],0,[OutputDirectoryName,'/map_cc',fsuffix,int2str(j)], pnames);
-            SAcc(j,imap)=gsa_(j).si;
-            imap_cc{j}=imap;
-        end
-        save([OutputDirectoryName,'/map_cc',fsuffix,'.mat'],'gsa_')
-        save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'imap_cc','SAcc','ccac','-append')
-    else
-        load([OutputDirectoryName,'/',fname_,'_main_eff'],'SAcc')
-    end
-
-    hh_fig=dyn_figure(options_.nodisplay,'Name',['Identifiability indices in the ',fsuffix,' moments.']);
-    bar(sum(SAcc))
-    set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
-    set(gca,'xlim',[0.5 npT+0.5])
-    ydum = get(gca,'ylim');
-    set(gca,'ylim',[0 ydum(2)])
-    set(gca,'position',[0.13 0.2 0.775 0.7])
-    for ip=1:npT
-        text(ip,-0.02*(ydum(2)),bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
-    end
-    xlabel(' ')
-    title(['Identifiability indices in the ',fsuffix,' moments.'],'interpreter','none')
-    dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_ident_ALL',fsuffix],options_.nodisplay,options_.graph_format);
-    create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_ident_ALL',fsuffix],1,['Identifiability indices in the ',fsuffix,' moments.'],['ident_ALL',fsuffix]',1)
+else  
+    error('gsa/map_identification: unsupported option morris=%u',opt_gsa.morris)
 end

 function []=create_TeX_loader(options_,figpath,ifig_number,caption,label_name,scale_factor)
--- a/matlab/+gsa/run.m
+++ b/matlab/+gsa/run.m
@ -13,7 +13,7 @@ function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
 % M. Ratto (2008), Analysing DSGE Models with Global Sensitivity Analysis, 
 % Computational Economics (2008), 31, pp. 115–139

-% Copyright © 2008-2023 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -101,6 +101,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
        disp('must be specified for RMSE analysis!');
        error('Sensitivity anaysis error!')
    end
+    if isfield(options_gsa,'nobs')
+        options_.nobs=options_gsa.nobs;
+    end
    if ~isempty(options_.nobs) && length(options_.nobs)~=1
        error('dynare_sensitivity does not support recursive estimation. Please specify nobs as a scalar, not a vector.')
    end
@ -108,9 +111,6 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
    if isfield(options_gsa,'first_obs')
        options_.first_obs=options_gsa.first_obs;
    end
-    if isfield(options_gsa,'nobs')
-        options_.nobs=options_gsa.nobs;
-    end
    if isfield(options_gsa,'presample')
        options_.presample=options_gsa.presample;
    end
@ -150,6 +150,11 @@ end

 [~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);

+if isfield(oo_.dr,'eigval') && any(abs(oo_.dr.eigval-1)<abs(1-options_.qz_criterium)) && options_.qz_criterium<1
+    fprintf('\ngsa: The model features a unit root, but qz_criterium<1. Check whether that is intended.')
+    fprintf('\ngsa: If not, use the diffuse_filter option.\n')
+end
+
 options_gsa = set_default_option(options_gsa,'identification',0);
 if options_gsa.identification
    options_gsa.redform=0;
@ -522,4 +527,4 @@ if options_gsa.glue
            save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
        end
    end
-end
+end
--- a/matlab/+identification/run.m
+++ b/matlab/+identification/run.m
@ -736,7 +736,7 @@ if iload <=0
                end
                run_index = 0; % reset index
            end
-            if SampleSize > 1
+            if SampleSize > 1 && mod(iteration,3)
                dyn_waitbar(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]);
            end
        end
--- a/matlab/+identification/simulated_moment_uncertainty.m
+++ b/matlab/+identification/simulated_moment_uncertainty.m
@ -12,7 +12,7 @@ function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic,options_
 %   - cmm:      [n_moments by n_moments] covariance matrix of simulated moments
 %   - mm:       [n_moments by replic] matrix of moments

-% Copyright © 2009-2018 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -60,15 +60,10 @@ end

 oo_.dr=set_state_space(oo_.dr,M_);

-
 if options_.logged_steady_state %if steady state was previously logged, undo this
    oo_.dr.ys=exp(oo_.dr.ys);
    oo_.steady_state=exp(oo_.steady_state);
    options_.logged_steady_state=0;
-    logged_steady_state_indicator=1;
-    evalin('base','options_.logged_steady_state=0;')
-else
-    logged_steady_state_indicator=0;
 end

 [oo_.dr,info,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
@ -92,7 +87,6 @@ else
    end
 end

-
 for j=1:replic
    [ys, oo_.exo_simul] = simult(y0,oo_.dr,M_,options_);%do simulation
    oo_=disp_moments(ys, options_.varobs,M_,options_,oo_); %get moments
@ -106,8 +100,5 @@ for j=1:replic
 end
 dyn_waitbar_close(h);

-if logged_steady_state_indicator
-    evalin('base','options_.logged_steady_state=1;') %reset base workspace option to conform to base oo_
-end
 cmm = cov(mm');
 disp('Simulated moment uncertainty ... done!')
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@ -479,7 +479,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
                        fprintf(fidTeX,'\\begin{figure}[H]\n');
                        fprintf(fidTeX,'\\centering \n');
                        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
-                        fprintf(fidTeX,'\\caption{Check plots.}');
+                        fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
                        fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
                        fprintf(fidTeX,'\\end{figure}\n');
                        fprintf(fidTeX,' \n');
@ -488,7 +488,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
            end
        end

-            if mod(j1,9)~=0 && j==M_.exo_nbr
+        if mod(j1,9)~=0 && j==M_.exo_nbr
                annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
                annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
                dyn_saveas(hh_fig,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
@ -497,7 +497,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
                fprintf(fidTeX,'\\begin{figure}[H]\n');
                fprintf(fidTeX,'\\centering \n');
                fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_smoothedshocks_occbin%s}\n',options_.figures.textwidth*min(j1/3,1),[GraphDirectoryName '/' M_.fname],int2str(ifig)); % don't use filesep as it will create issues with LaTeX on Windows
-                fprintf(fidTeX,'\\caption{Check plots.}');
+                fprintf(fidTeX,'\\caption{OccBin smoothed shocks.}');
                fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
                fprintf(fidTeX,'\\end{figure}\n');
                fprintf(fidTeX,' \n');
@ -505,6 +505,6 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
        end
        if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
            fclose(fidTeX);
+        end
    end
 end
-end
--- a/matlab/+occbin/forecast.m
+++ b/matlab/+occbin/forecast.m
@ -1,143 +1,169 @@
-function [oo_, error_flag] = forecast(options_,M_,oo_,forecast) %,hist_period)
-%function oo_ = forecast(options_,M_,oo_,forecast)
-% forecast
+function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
+% [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
+% Occbin forecasts
+%
+% INPUTS
+% - options_                [structure]     Matlab's structure describing the current options
+% - M_                      [structure]     Matlab's structure describing the model
+% - dr_in                   [structure]     model information structure
+% - endo_steady_state       [double]        steady state value for endogenous variables
+% - exo_steady_state        [double]        steady state value for exogenous variables
+% - exo_det_steady_state    [double]        steady state value for exogenous deterministic variables                                    
+% - forecast_horizon        [integer]       forecast horizon
+%
+% OUTPUTS
+% - forecast                [structure]     forecast results
+% - error_flag              [integer]       error code
+%
+% SPECIAL REQUIREMENTS
+%   none.
+
+% Copyright © 2022-2023 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/>.

 opts = options_.occbin.forecast;

 options_.occbin.simul.maxit = opts.maxit;
 options_.occbin.simul.check_ahead_periods = opts.check_ahead_periods;
-options_.occbin.simul.periods = forecast;
-SHOCKS0 = opts.SHOCKS0;
-if ~isempty(SHOCKS0)
-    if iscell(SHOCKS0)
-        for j=1:length(SHOCKS0)
-            sname = SHOCKS0{j}{1};
-            inds(j)=strmatch(sname,M_.exo_names);
-            SHOCKS1(j,:)=SHOCKS0{j}{2};
+options_.occbin.simul.periods = forecast_horizon;
+shocks_input = opts.SHOCKS0;
+
+if ~isempty(shocks_input)
+    n_shocks=size(shocks_input,1);
+    if iscell(shocks_input)
+        inds=NaN(n_shocks,1);
+        periods=length(shocks_input{1}{2});
+        shock_mat=NaN(n_shocks,periods);
+        for j=1:n_shocks
+            exo_pos=strmatch(shocks_input{j}{1},M_.exo_names,'exact');
+            if isempty(exo_pos)
+        	    error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
+            else
+                inds(j)=exo_pos;
+            end
+            if length(shocks_input{j}{2})~=periods
+        	    error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
+            else
+                shock_mat(j,:)=shocks_input{j}{2};
+            end
        end
-    elseif isreal(SHOCKS0)
-        SHOCKS1=SHOCKS0;
-        inds = 1:M_.exo_nbr;
+    elseif isreal(shocks_input)
+        shock_mat=shocks_input;
+        inds = (1:M_.exo_nbr)';
    end
 end

+options_.occbin.simul.endo_init = M_.endo_histval(:,1)-endo_steady_state; %initial condition
+options_.occbin.simul.init_regime = opts.frcst_regimes;
+options_.occbin.simul.init_binding_indicator = [];
+
+shocks_base = zeros(forecast_horizon,M_.exo_nbr);
+if ~isempty(shocks_input)
+    for j=1:n_shocks
+        shocks_base(:,inds(j))=shock_mat(j,:);
+    end
+end

 if opts.replic
    h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
    ishock = find(sqrt(diag((M_.Sigma_e))));
+    options_.occbin.simul.exo_pos=ishock;
    effective_exo_nbr=  length(ishock);
-    effective_exo_names = M_.exo_names(ishock);
-    effective_Sigma_e = M_.Sigma_e(ishock,ishock);
+    effective_Sigma_e = M_.Sigma_e(ishock,ishock);  % does not take heteroskedastic shocks into account
    [U,S] = svd(effective_Sigma_e);
+    % draw random shocks
    if opts.qmc
        opts.replic =2^(round(log2(opts.replic+1)))-1;
-        SHOCKS_ant =   qmc_sequence(forecast*effective_exo_nbr, int64(1), 1, opts.replic)';
+        SHOCKS_add = qmc_sequence(forecast_horizon*effective_exo_nbr, int64(1), 1, opts.replic);
    else
-        SHOCKS_ant = randn(forecast*effective_exo_nbr,opts.replic)';
+        SHOCKS_add = randn(forecast_horizon*effective_exo_nbr,opts.replic);
    end
-    zlin0=zeros(forecast,M_.endo_nbr,opts.replic);
-    zpiece0=zeros(forecast,M_.endo_nbr,opts.replic);
+    SHOCKS_add=reshape(SHOCKS_add,effective_exo_nbr,forecast_horizon,opts.replic);
+    z.linear=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
+    z.piecewise=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
+    error_flag=true(opts.replic,1);
+    simul_SHOCKS=NaN(forecast_horizon,M_.exo_nbr,opts.replic);
    for iter=1:opts.replic
-
-        if ~isempty(SHOCKS0)
-            for j=1:length(SHOCKS0)
-                SHOCKS(:,inds(j))=SHOCKS1(j,:);
-            end
-        end
-
-        error_flagx=1;
-        % while error_flagx,
-        % SHOCKS=transpose(sqrt(diag(diag(effective_Sigma_e)))*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
-        SHOCKS=transpose(U*sqrt(S)*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
-        %             SHOCKS=transpose(U*sqrt(S)*randn(forecast,M_.exo_nbr)');   %realized shocks
-        options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
-        options_.occbin.simul.init_regime = opts.frcst_regimes;
-        options_.occbin.simul.init_binding_indicator = [];
-        options_.occbin.simul.exo_pos=ishock;
-        options_.occbin.simul.SHOCKS = SHOCKS;
+        options_.occbin.simul.SHOCKS = shocks_base+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
        options_.occbin.simul.waitbar=0;
-        [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-        zlin0(:,:,iter)=out.linear;
-        zpiece0(:,:,iter)=out.piecewise;
-        ys=out.ys;
-        frcst_regime_history(iter,:)=out.regime_history;
+        [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
        error_flag(iter)=out.error_flag;
-        error_flagx = error_flag(iter);
-        % end
-        simul_SHOCKS(:,:,iter) = SHOCKS;
-
-        if error_flag(iter) && debug_flag
-            %             display('no solution')
-
-%             keyboard;
-            save no_solution SHOCKS zlin0 zpiece0 iter frcst_regime_history
+        if ~error_flag(iter)
+            z.linear(:,:,iter)=out.linear;
+            z.piecewise(:,:,iter)=out.piecewise;
+            frcst_regime_history(iter,:)=out.regime_history;
+            error_flag(iter)=out.error_flag;
+            simul_SHOCKS(:,:,iter) = shocks_base;
+        else
+            if options_.debug
+                save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag','out','shocks_base')
+            end
        end
        dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
    end
    dyn_waitbar_close(h);
-    save temp zlin0 zpiece0 simul_SHOCKS error_flag
+    if options_.debug
+         save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag')
+    end
    inx=find(error_flag==0);
-    zlin0=zlin0(:,:,inx);
-    zpiece0=zpiece0(:,:,inx);
-    zlin   = mean(zlin0,3);
-    zpiece = mean(zpiece0,3);
-    zpiecemin = min(zpiece0,[],3);
-    zpiecemax = max(zpiece0,[],3);
-    zlinmin = min(zlin0,[],3);
-    zlinmax = max(zlin0,[],3);
-
-    for i=1:M_.endo_nbr
-        for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),options_.forecasts.conf_sig);
+    z.linear=z.linear(:,:,inx);
+    z.piecewise=z.piecewise(:,:,inx);
+    z.min.piecewise = min(z.piecewise,[],3);
+    z.max.piecewise = max(z.piecewise,[],3);
+    z.min.linear = min(z.linear,[],3);
+    z.max.linear = max(z.linear,[],3);
+    
+    field_names={'linear','piecewise'};
+    post_mean=NaN(forecast_horizon,1);
+    post_median=NaN(forecast_horizon,1);
+    post_var=NaN(forecast_horizon,1);
+    hpd_interval=NaN(forecast_horizon,2);
+    post_deciles=NaN(forecast_horizon,9);
+    for field_iter=1:2
+        for i=1:M_.endo_nbr
+            for j=1:forecast_horizon
+                [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(z.(field_names{field_iter})(j,i,:)),options_.forecasts.conf_sig);
+            end
+            forecast.(field_names{field_iter}).Mean.(M_.endo_names{i})=post_mean;
+            forecast.(field_names{field_iter}).Median.(M_.endo_names{i})=post_median;
+            forecast.(field_names{field_iter}).Var.(M_.endo_names{i})=post_var;
+            forecast.(field_names{field_iter}).HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
+            forecast.(field_names{field_iter}).HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
+            forecast.(field_names{field_iter}).Deciles.(M_.endo_names{i})=post_deciles;
+            forecast.(field_names{field_iter}).Min.(M_.endo_names{i})=z.min.(field_names{field_iter})(:,i);
+            forecast.(field_names{field_iter}).Max.(M_.endo_names{i})=z.max.(field_names{field_iter})(:,i);
        end
-        oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
-        oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
-        oo_.occbin.linear_forecast.Var.(M_.endo_names{i})=post_var;
-        oo_.occbin.linear_forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
-        oo_.occbin.linear_forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
-        oo_.occbin.linear_forecast.Deciles.(M_.endo_names{i})=post_deciles;
-        oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
-        oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
-        for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),options_.forecasts.conf_sig);
-        end
-        oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
-        oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;
-        oo_.occbin.forecast.Var.(M_.endo_names{i})=post_var;
-        oo_.occbin.forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
-        oo_.occbin.forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
-        oo_.occbin.forecast.Deciles.(M_.endo_names{i})=post_deciles;
-        oo_.occbin.forecast.Min.(M_.endo_names{i})=zpiecemin(:,i);
-        oo_.occbin.forecast.Max.(M_.endo_names{i})=zpiecemax(:,i);
-        %   eval([M_.endo_names{i},'_ss=zdatass(i);']);
    end
-
 else
-    SHOCKS = zeros(forecast,M_.exo_nbr);
-    if ~isempty(SHOCKS0)
-        for j=1:length(SHOCKS0)
-            SHOCKS(:,inds(j))=SHOCKS1(j,:);
-        end
-    end
-    options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
-    options_.occbin.simul.init_regime = opts.frcst_regimes;
-    options_.occbin.simul.init_violvecbool = [];
    options_.occbin.simul.irfshock = M_.exo_names;
-    options_.occbin.simul.SHOCKS = SHOCKS;
-    [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+    options_.occbin.simul.SHOCKS = shocks_base;
+    [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
+    error_flag=out.error_flag;
+    if out.error_flag
+        fprintf('occbin.forecast: forecast simulation failed.')
+        return;
+    end

-    zlin=out.linear;
-    zpiece=out.piecewise;
    frcst_regime_history=out.regime_history;
    error_flag=out.error_flag;
    for i=1:M_.endo_nbr
-        oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})= zlin(:,i);
-        oo_.occbin.forecast.Mean.(M_.endo_names{i})= zpiece(:,i);
-        oo_.occbin.forecast.HPDinf.(M_.endo_names{i})= nan;
-        oo_.occbin.forecast.HPDsup.(M_.endo_names{i})= nan;
+        forecast.linear.Mean.(M_.endo_names{i})= out.linear(:,i);
+        forecast.piecewise.Mean.(M_.endo_names{i})= out.piecewise(:,i);
    end
-
 end

-oo_.occbin.forecast.regimes=frcst_regime_history;
-
+forecast.regimes=frcst_regime_history;
--- a/matlab/+occbin/irf.m
+++ b/matlab/+occbin/irf.m
@ -1,140 +1,124 @@
-function [oo_] = irf(M_,oo_,options_)
+function irfs = irf(M_,oo_,options_)
+% irfs = irf(M_,oo_,options_)
+% Calls a minimizer
+%
+% INPUTS
+% - M_                  [structure]     Matlab's structure describing the model
+% - oo_                 [structure]     Matlab's structure containing the results
+% - options_            [structure]     Matlab's structure describing the current options
+%
+% OUTPUTS
+% - irfs                [structure]     IRF results
+%
+% SPECIAL REQUIREMENTS
+%   none.
+%
+%
+% Copyright © 2022-2023 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/>.

 shocknames = options_.occbin.irf.exo_names;
-shocksigns = options_.occbin.irf.shocksigns;
+shocksigns = options_.occbin.irf.shocksigns; %'pos','neg'
 shocksize = options_.occbin.irf.shocksize;
-t0 = options_.occbin.irf.t0;
+t_0 = options_.occbin.irf.t0;
+
+%% set simulation options based on IRF options
 options_.occbin.simul.init_regime = options_.occbin.irf.init_regime;
 options_.occbin.simul.check_ahead_periods = options_.occbin.irf.check_ahead_periods;
 options_.occbin.simul.maxit = options_.occbin.irf.maxit;
 options_.occbin.simul.periods = options_.irf;

-% Run inital conditions + other shocks
-if t0 == 0
-    shocks0= zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
+%% Run initial conditions + other shocks
+if t_0 == 0
+    shocks_base = zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
    options_.occbin.simul.endo_init = [];
 else
-    % girf conditional to smoothed states in t0 and shocks in t0+1
-    shocks0= [oo_.occbin.smoother.etahat(:,t0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
-    options_.occbin.simul.SHOCKS=shocks0;
-    options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
+    if ~isfield(oo_.occbin,'smoother')
+        error('occbin.irfs: smoother must be run before requesting GIRFs based on smoothed results')
+    end
+    % GIRF conditional on smoothed states in t_0 and shocks in t_0+1
+    shocks_base= [oo_.occbin.smoother.etahat(:,t_0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
+    options_.occbin.simul.SHOCKS=shocks_base;
+    options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
 end
-options_.occbin.simul.SHOCKS=shocks0;
-[~, out0] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-zlin0 = out0.linear;
-zpiece0 = out0.piecewise;
+options_.occbin.simul.SHOCKS=shocks_base;

-% Select shocks of interest
-jexo_all =zeros(size(shocknames,1),1);
+[~, out_base] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+if out_base.error_flag
+    error('occbin.irfs: could not compute the solution')
+end

+irfs.linear = struct();
+irfs.piecewise = struct();
+
+% Get indices of shocks of interest
+exo_index =zeros(size(shocknames,1),1);
 for i=1:length(shocknames)
-    jexo_all(i) = strmatch(shocknames{i},M_.exo_names,'exact');
+    exo_index(i) = strmatch(shocknames{i},M_.exo_names,'exact');
 end

-oo_.occbin.linear_irfs = struct();
-oo_.occbin.irfs = struct();
+% cs=get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
+% irf_shocks_indx = getIrfShocksIndx(M_, options_);
+
 % Set shock size
 if isempty(shocksize)
-    %         if isfield(oo_.posterior_mode.shocks_std,M_.exo_names{jexo})
-    shocksize = sqrt(diag(M_.Sigma_e(jexo_all,jexo_all))); %oo_.posterior_mode.shocks_std.(M_.exo_names{jexo});
+    shocksize = sqrt(diag(M_.Sigma_e(exo_index,exo_index)));
    if any(shocksize < 1.e-9)
        shocksize(shocksize < 1.e-9) = 0.01;
    end
 end
+
 if numel(shocksize)==1
    shocksize=repmat(shocksize,[length(shocknames),1]);
 end

 % Run IRFs
-for counter = 1:length(jexo_all)
-  
-    jexo = jexo_all(counter);
-    
-    if ~options_.noprint
-        fprintf('Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},counter,size(jexo_all,1));
-    end 
-    
-    if ismember('pos',shocksigns)
-        % (+) shock
-        shocks1=shocks0;
-        shocks1(1,jexo)=shocks0(1,jexo)+shocksize(counter);
-        if t0 == 0
-            options_.occbin.simul.SHOCKS=shocks1;
+for sign_iter=1:length(shocksigns)
+    for IRF_counter = 1:length(exo_index)
+        jexo = exo_index(IRF_counter);
+        if ~options_.noprint && options_.debug
+            fprintf('occbin.irf: Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},IRF_counter,size(exo_index,1));
+        end
+        shocks1=shocks_base;
+        if ismember('pos',shocksigns{sign_iter})
+            shocks1(1,jexo)=shocks_base(1,jexo)+shocksize(IRF_counter);
+        elseif ismember('neg',shocksigns{sign_iter})
+            shocks1(1,jexo)=shocks_base(1,jexo)-shocksize(IRF_counter);
+        end
+        options_.occbin.simul.SHOCKS=shocks1;
+        if t_0 == 0
            options_.occbin.simul.endo_init = [];
-            [~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
        else
-            options_.occbin.simul.SHOCKS=shocks1;
-            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
-            [~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
        end
-        if out_pos.error_flag
-            warning('Occbin error.')
-            return
+        [~, out_sim] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+        if out_sim.error_flag
+            warning('occbin.irfs: simulation failed')
+            skip
        end
-        zlin_pos = out_pos.linear;
-        zpiece_pos = out_pos.piecewise;
        % Substract inital conditions + other shocks
-        zlin_pos_diff       = zlin_pos-zlin0;
-        zpiece_pos_diff      = zpiece_pos-zpiece0;
-    end
-    
-    if ismember('neg',shocksigns)
-        % (-) shock
-        shocks_1=shocks0;
-        shocks_1(1,jexo)=shocks0(1,jexo)-shocksize(counter);
-        if t0 == 0
-            options_.occbin.simul.SHOCKS=shocks_1;
-            options_.occbin.simul.endo_init = [];
-            [~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-        else
-            options_.occbin.simul.SHOCKS=shocks_1;
-            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
-            [~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+        zdiff.linear.(shocksigns{sign_iter}) = out_sim.linear-out_base.linear;
+        zdiff.piecewise.(shocksigns{sign_iter}) = out_sim.piecewise-out_base.piecewise;
+
+        for j_endo=1:M_.endo_nbr
+            if ismember('pos',shocksigns)
+                irfs.piecewise.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}])   = zdiff.piecewise.(shocksigns{sign_iter})(:,j_endo);
+                irfs.linear.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}])   = zdiff.linear.(shocksigns{sign_iter})(:,j_endo);
+            end
        end
-        if out_neg.error_flag
-            warning('Occbin error.')
-            return
-        end
-        zlin_neg    = out_neg.linear;
-        zpiece_neg  = out_neg.piecewise;
-        zlin_neg_diff    = zlin_neg-zlin0;
-        zpiece_neg_diff  = zpiece_neg-zpiece0;
    end
-
-    % Save
-    if ~isfield(oo_.occbin,'linear_irfs')
-        oo_.occbin.linear_irfs = struct();
-    end
-    if ~isfield(oo_.occbin,'irfs')
-        oo_.occbin.irfs = struct();
-    end
-
-    for jendo=1:M_.endo_nbr
-%         oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zpiece_pos (:,jendo);
-%         oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zpiece_neg (:,jendo);
-%         oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zlin_pos (:,jendo);
-%         oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zlin_neg(:,jendo);
-   
-        if ismember('pos',shocksigns)
-        oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos'])   = zpiece_pos_diff (:,jendo);
-        oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos'])   = zlin_pos_diff (:,jendo);
-        end
-        
-        if ismember('neg',shocksigns)
-            oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg'])  = zpiece_neg_diff (:,jendo);
-            oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg'])  = zlin_neg_diff (:,jendo);
-        end
-        
-% %         
-%         oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zpiece0(:,jendo);
-%         oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])  = zlin0(:,jendo);
-%         oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])   = zpiece0(:,jendo);
-%         oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zlin0(:,jendo);
-
-    end
-end
-
-
-
-end
-
+end
--- a/matlab/+occbin/plot_irfs.m
+++ b/matlab/+occbin/plot_irfs.m
@ -1,34 +1,74 @@
-function  plot_irfs(M_,oo_,options_,irf3,irf4)
+function plot_irfs(M_,irfs,options_,var_list)
+% plot_irfs(M_,irfs,options_,var_list)
+%
+% INPUTS
+% - M_                      [structure]     Matlab's structure describing the model
+% - irfs                    [structure]     IRF results
+% - options_                [structure]     Matlab's structure describing the current options
+% - var_list                [character array]  list of endogenous variables specified
+%
+% OUTPUTS
+%   none
+%
+% SPECIAL REQUIREMENTS
+%   none.

-shocknames = options_.occbin.plot_irf.exo_names;
-simulname = options_.occbin.plot_irf.simulname;
-if isempty(simulname)
-   simulname_ = simulname;
-else 
-   simulname_ = [ simulname '_' ];
+% Copyright © 2022-2023 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/>.
+
+if (isfield(options_,'irf_shocks')==0)
+    shock_names  = M_.exo_names;
+else
+    shock_names  = options_.irf_shocks;
 end
-vars_irf = options_.occbin.plot_irf.endo_names;
-endo_names_long = options_.occbin.plot_irf.endo_names_long;
+
+simul_name = options_.occbin.plot_irf.simulname;
+if isempty(simul_name)
+    save_name = simul_name;
+else
+    save_name = [ simul_name '_' ];
+end
+
+if isempty(var_list)
+    var_list = M_.endo_names(1:M_.orig_endo_nbr);
+end
+
+[i_var, ~, index_uniques] = varlist_indices(var_list, M_.endo_names);
+vars_irf=var_list(index_uniques);
+
 endo_scaling_factor = options_.occbin.plot_irf.endo_scaling_factor;
-length_irf = options_.occbin.plot_irf.tplot;
-if isempty(length_irf)
-    length_irf = options_.irf;
-end
+length_irf = options_.irf;

-irflocation_lin   = oo_.occbin.linear_irfs;
-irflocation_piece = oo_.occbin.irfs;
-
-
-steps_irf = 1; 
-warning('off','all')
+steps_irf = 1;

 DirectoryName = CheckPath('graphs',M_.dname);
+latexFolder = CheckPath('latex',M_.dname);
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fidTeX = fopen([latexFolder '/' M_.fname '_occbin_irfs.tex'],'w');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by occbin.plot_irfs.m (Dynare).\n');
+    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+    fprintf(fidTeX,' \n');
+end

 iexo=[];
-for i=1:size(shocknames,1)
-    itemp = strmatch(shocknames{i},M_.exo_names,'exact');
+for var_iter=1:size(shock_names,1)
+    itemp = strmatch(shock_names{var_iter},M_.exo_names,'exact');
    if isempty(itemp)
-        error(['Shock ',shocknames{i},' is not defined!'])
+        error(['Shock ',shock_names{var_iter},' is not defined!'])
    else
        iexo=[iexo, itemp];
    end
@ -38,104 +78,102 @@ ncols     = options_.occbin.plot_irf.ncols;
 nrows     = options_.occbin.plot_irf.nrows;
 npan      = ncols*nrows;

-plot_grid = options_.occbin.plot_irf.grid;
-shocksigns = options_.occbin.plot_irf.shocksigns; 
-threshold = options_.occbin.plot_irf.threshold; 
-
-% Add steady_state
-if options_.occbin.plot_irf.add_steadystate
-    add_stst = options_.occbin.plot_irf.add_steadystate;
-else
-    add_stst = 0;
-end 
-for sss = 1:numel(shocksigns)
-    
-    shocksign = shocksigns{sss};
-    
-    for j=1:size(shocknames,1)
-        %shocknames = M_.exo_names{j};
+shocksigns = options_.occbin.plot_irf.shocksigns;
+threshold = options_.impulse_responses.plot_threshold;

+for shock_sign_iter = 1:numel(shocksigns)
+    shocksign = shocksigns{shock_sign_iter};
+    if strcmp(shocksign,'pos')
+        plot_title_sign='positive';
+    elseif strcmp(shocksign,'neg')
+        plot_title_sign='negative';
+    else
+       error('Unknown shock sign %s',shocksign);
+    end
+        
+    for shock_iter=1:size(shock_names,1)
        j1   = 0;
        isub = 0;
        ifig = 0;
-
        % Variables
-        % ----------------------
-        for i = 1:length(vars_irf)
-
+        for var_iter = 1:length(vars_irf)
            j1=j1+1;
            if mod(j1,npan)==1
-               % vector corresponds to [left bottom width height]. 680 and 678 for the left and bottom elements correspond to the default values used by MATLAB while creating a figure and width, .
-                hfig = dyn_figure(options_.nodisplay,'name',['OccbinIRFs ' shocknames{j} ' ' simulname ' ' shocksign],'PaperPositionMode', 'auto','PaperType','A4','PaperOrientation','portrait','renderermode','auto','position',[10 10 950 650]);
+                % vector corresponds to [left bottom width height]. 680 and 678 for the left and bottom elements correspond to the default values used by MATLAB while creating a figure and width, .
+                screensize = get( groot, 'Screensize' );
+                hfig = dyn_figure(options_.nodisplay,'name',['OccBin IRFs to ' plot_title_sign ' ' shock_names{shock_iter} ' shock ' simul_name],'OuterPosition' ,[50 50 min(1000,screensize(3)-50) min(750,screensize(4)-50)]);
                ifig=ifig+1;
                isub=0;
            end
-            isub=isub+1;      
-            
+            isub=isub+1;
            if isempty(endo_scaling_factor)
                exofactor = 1;
-            else               
-                exofactor = endo_scaling_factor{i};
-            end 
-
-            subplot(nrows,ncols,isub)            
-            irf_field   = strcat(vars_irf{i,1},'_',shocknames{j},'_',shocksign);
-            
-            irfvalues   = irflocation_lin.(irf_field);
-            if add_stst
-                irfvalues = irfvalues + get_mean(vars_irf{i,1});
-            end          
-            
-            irfvalues(abs(irfvalues) <threshold) = 0;
-            
-            plot(irfvalues(1:steps_irf:length_irf)*exofactor,'linewidth',2);
-            hold on
-
-            irfvalues   = irflocation_piece.(irf_field);  
-            if add_stst
-                irfvalues = irfvalues + get_mean(vars_irf{i,1});
-            end
-            irfvalues(abs(irfvalues) <threshold) = 0;
-            plot(irfvalues(1:steps_irf:length_irf)*exofactor,'r--','linewidth',2);
-            
-            hold on
-            plot(irfvalues(1:steps_irf:length_irf)*0,'k-','linewidth',1.5);
-            % Optional additional IRFs
-            if nargin > 10
-                  irfvalues   = irf3.(irf_field) ;
-                  irfvalues(abs(irfvalues) <threshold) = 0;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'k:','linewidth',2);
-            end
-            if nargin > 11
-                  irfvalues   = irf4.(irf_field)   ;   
-                  irfvalues(abs(irfvalues) <threshold) = 0;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'g-.','linewidth',2);
-            end
-      
-
-            if plot_grid
-                grid on
-            end
-
-            xlim([1 (length_irf/steps_irf)]);
-            
-            % title
-            if isempty(endo_names_long)
-                title(regexprep(vars_irf{i},'_',' '))
            else
-                title(endo_names_long{i})
+                exofactor = endo_scaling_factor{var_iter};
            end

+            subplot(nrows,ncols,isub)
+            irf_field   = strcat(vars_irf{var_iter},'_',shock_names{shock_iter},'_',shocksign);
+
+            %%linear IRFs
+            if ~isfield(irfs.linear,irf_field)
+                warning('occbin.plot_irfs: no linear IRF for %s to %s detected',vars_irf{var_iter,1},shock_names{shock_iter})
+            else
+                irfvalues   = irfs.linear.(irf_field);
+                irfvalues(abs(irfvalues) <threshold) = 0;
+                if options_.occbin.plot_irf.add_steadystate
+                    irfvalues = irfvalues + get_mean(vars_irf{var_iter,1});
+                end
+                max_irf_length_1=min(length_irf,length(irfvalues));
+                plot(irfvalues(1:steps_irf:max_irf_length_1)*exofactor,'linewidth',2);
+            end
+            hold on
+
+            if ~isfield(irfs.piecewise,irf_field)
+                warning('occbin.plot_irfs: no piecewise linear IRF for %s to %s detected',vars_irf{var_iter,1},shock_names{shock_iter})
+            else
+                irfvalues   = irfs.piecewise.(irf_field);
+                irfvalues(abs(irfvalues) <threshold) = 0;
+                if options_.occbin.plot_irf.add_steadystate
+                    irfvalues = irfvalues + get_mean(vars_irf{var_iter,1});
+                end
+                max_irf_length_2=min(length_irf,length(irfvalues));
+                plot(irfvalues(1:steps_irf:max_irf_length_2)*exofactor,'r--','linewidth',2);
+            end
+
+            plot(irfvalues(1:steps_irf:max(max_irf_length_1,max_irf_length_2))*0,'k-','linewidth',1.5);
+
+            if options_.occbin.plot_irf.grid
+                grid on
+            end     
+            xlim([1 max(max_irf_length_1,max_irf_length_2)]);
+            if options_.TeX
+                title(['$' M_.endo_names_tex{i_var(var_iter)}, '$'],'Interpreter','latex')
+            else
+                title(M_.endo_names{i_var(var_iter)},'Interpreter','none')
+            end
            % Annotation Box + save figure
            % ----------------------
-            if mod(j1,npan)==0 || (mod(j1,npan)~=0 && i==length(vars_irf))
+            if mod(j1,npan)==0 || (mod(j1,npan)~=0 && var_iter==length(vars_irf))
                annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
                annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
-                dyn_saveas(hfig,[DirectoryName,filesep,M_.fname,'_irf_occbin_',simulname_,shocknames{j},'_',shocksign,'_',int2str(ifig)],options_.nodisplay,options_.graph_format);
+                dyn_saveas(hfig,[DirectoryName,filesep,M_.fname,'_irf_occbin_',save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)],options_.nodisplay,options_.graph_format);
+                if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+                    fprintf(fidTeX,'\\begin{figure}[H]\n');
+                    fprintf(fidTeX,'\\centering \n');
+                    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_irf_occbin_%s}\n',options_.figures.textwidth*min((j1-1)/ncols,1),...
+                        [DirectoryName '/' ,M_.fname],[save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)]);
+                    fprintf(fidTeX,'\\caption{OccBin IRFs to  %s shock to %s}\n',plot_title_sign,shock_names{shock_iter});
+                    fprintf(fidTeX,'\\label{Fig:occbin_irfs:%s}\n',[save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)]);
+                    fprintf(fidTeX,'\\end{figure}\n');
+                    fprintf(fidTeX,'\n');
+                end
            end
-        end 
+        end
    end
 end

-warning('on','all')
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fprintf(fidTeX,'%% End Of TeX file.');
+    fclose(fidTeX);
 end
--- a/matlab/+occbin/plot_regimes.m
+++ b/matlab/+occbin/plot_regimes.m
@ -1,4 +1,26 @@
 function plot_regimes(regimes,M_,options_)
+% plot_regimes(regimes,M_,options_)
+% Inputs: 
+% - regimes     	[structure]     OccBin regime information
+% - M_              [structure]     Matlab's structure describing the model
+% - options_        [structure]     Matlab's structure containing the options
+
+% Copyright © 2021-2023 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/>.

 nperiods = size(regimes,2);
 nconstr = length(fieldnames(regimes(1)))/2;
@ -13,9 +35,15 @@ else
 end    

 GraphDirectoryName = CheckPath('graphs',M_.dname);
+fhandle = dyn_figure(options_.nodisplay,'Name',[M_.fname ': OccBin regimes']);

-fhandle = dyn_figure(options_.nodisplay,'Name',[M_.fname ' occbin regimes']);
-
+latexFolder = CheckPath('latex',M_.dname);
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fidTeX = fopen([latexFolder '/' M_.fname '_occbin_regimes.tex'],'w');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by occbin.plot_regimes.m (Dynare).\n');
+    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+    fprintf(fidTeX,' \n');
+end

 for k=1:nconstr
    subplot(nconstr,1,k)
@ -36,12 +64,23 @@ for k=1:nconstr
            end
        end
    end
-    title(['regime ' int2str(k)])
-    xlabel('historic period')
-    ylabel('regime expected start')
+    xlim([1 nperiods])
+    title(['Regime ' int2str(k)])
+    xlabel('Historic period')
+    ylabel('Regime: expected start')
 end
-annotation('textbox',[.25,0,.15,.05],'String','Unbinding','Color','blue');
-annotation('textbox',[.65,0,.15,.05],'String','Binding','Color','red');
-
+annotation('textbox',[.25,0,.15,.05],'String','Slack','Color','blue');
+annotation('textbox',[.65,0,.2,.05],'String','Binding','Color','red');

 dyn_saveas(fhandle,[GraphDirectoryName, filesep, M_.fname '_occbin_regimes'],options_.nodisplay,options_.graph_format);
+
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fprintf(fidTeX,'\\begin{figure}[H]\n');
+    fprintf(fidTeX,'\\centering \n');
+    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_occbin_regimes}\n',options_.figures.textwidth,[GraphDirectoryName '/' M_.fname]);
+    fprintf(fidTeX,'\\caption{OccBin: regime evolution over time.}\n');
+    fprintf(fidTeX,'\\label{Fig:occbin_regimes}\n');
+    fprintf(fidTeX,'\\end{figure}\n');
+    fprintf(fidTeX,'\n');
+    fclose(fidTeX);
+end
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@ -55,9 +55,7 @@ if ismember(flag,{'forecast','all'})
    options_occbin_.forecast.maxit=30;
    options_occbin_.forecast.qmc=0;
    options_occbin_.forecast.replic=0;
-    options_occbin_.forecast.sepath=0;    
    options_occbin_.forecast.SHOCKS0=[];
-    options_occbin_.forecast.treepath=1; % number of branches
 end

 if ismember(flag,{'irf','all'})    
@ -98,17 +96,12 @@ end

 if ismember(flag,{'plot_irf','all'})
    options_occbin_.plot_irf.add_steadystate = 0;
-    options_occbin_.plot_irf.exo_names = [];
-    options_occbin_.plot_irf.endo_names = M_.endo_names;
-    options_occbin_.plot_irf.endo_names_long = [];
    options_occbin_.plot_irf.endo_scaling_factor = [];
    options_occbin_.plot_irf.grid            = true;
    options_occbin_.plot_irf.ncols            = 3;
    options_occbin_.plot_irf.nrows            = 3;
    options_occbin_.plot_irf.shocksigns = {'pos','neg'}; 
    options_occbin_.plot_irf.simulname='';
-    options_occbin_.plot_irf.threshold = 10^-6;
-    options_occbin_.plot_irf.tplot = [];
 end

 if ismember(flag,{'plot_shock_decomp','all'})
--- a/matlab/collapse_figures_in_tabgroup.m
+++ b/matlab/collapse_figures_in_tabgroup.m
@ -0,0 +1,42 @@
+function collapse_figures_in_tabgroup
+
+% Copyright © 2023 Eduard Benet Cerda
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+% Create a new figure with results
+fig = uifigure(Name = 'Dynare Results');
+
+% Add a grid layout to make sure it spans the entire width
+g = uigridlayout(fig, [1,1], Padding = 0);
+
+% Add a tabgroup
+tg = uitabgroup(g);
+
+% Find all figures with Dynare Tag 
+f = findobj('-regexp','tag','dynare-figure');
+
+% Loop over all figures and reparent them to a tab. Avoid legends, they are
+% automatically tied.
+for j = 1 : numel(f)
+    t = uitab(tg);
+    types = arrayfun(@class, f(j).Children, 'UniformOutput', false);
+    idx = ismember(types, 'matlab.graphics.illustration.Legend'); % no need to reparent legends
+    set(f(j).Children(~idx),'Parent',t)
+    t.Title = f(j).Name;
+    delete(f(j))
+end
--- a/matlab/dyn_figure.m
+++ b/matlab/dyn_figure.m
@ -12,7 +12,7 @@ function h = dyn_figure(nodisplay, varargin)
 % SPECIAL REQUIREMENTS
 %    none

-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -30,7 +30,7 @@ function h = dyn_figure(nodisplay, varargin)
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

 if nodisplay
-    h = figure(varargin{:},'visible','off');
+    h = figure(varargin{:},'visible','off','tag','dynare-figure');
 else
-    h = figure(varargin{:});
+    h = figure(varargin{:},'tag','dynare-figure');
 end
--- a/matlab/dynare.m
+++ b/matlab/dynare.m
@ -122,7 +122,7 @@ if isempty(dot_location)
    fnamelength = length(fname);
    fname1 = [fname '.dyn'];
    d = dir(fname1);
-    if length(d) == 0
+    if isempty(d)
        fname1 = [fname '.mod'];
    end
    fname = fname1;
@ -138,7 +138,7 @@ if fnamelength + length('.set_auxiliary_variables') > namelengthmax()
    error('Dynare: the name of your .mod file is too long, please shorten it')
 end

-if ~isempty(strfind(fname,filesep))
+if contains(fname,filesep)
    fprintf('\nIt seems you are trying to call a .mod file not located in the "Current Folder". This is not possible (the %s symbol is not allowed in the name of the .mod file).\n', filesep)
    [pathtomodfile,basename] = fileparts(fname);
    if exist(pathtomodfile,'dir')
@ -297,17 +297,13 @@ if status
    error('Dynare: preprocessing failed')
 end

-if ~ isempty(find(abs(fname) == 46))
-    fname = fname(:,1:find(abs(fname) == 46)-1) ;
-end
-
 % We need to clear the driver (and only the driver, because the "clear all"
 % within the driver will clean the rest)
-clear(['+' fname '/driver'])
+clear(['+' fname(1:end-4) '/driver'])

 try
    cTic = tic;
-    evalin('base',[fname '.driver']);
+    evalin('base',[fname(1:end-4) '.driver']);
    cToc = toc(cTic);
    if nargout
        DynareInfo.time.compute = cToc;
@ -315,7 +311,7 @@ try
 catch ME
    W = evalin('caller','whos');
    diary off
-    if ismember(fname,{W(:).name})
+    if ismember(fname(1:end-4),{W(:).name})
        error('Your workspace already contains a variable with the same name as the mod-file. You need to delete it or rename the mod-file.')
    else
        rethrow(ME)
--- a/matlab/dynare_graph.m
+++ b/matlab/dynare_graph.m
@ -1,54 +0,0 @@
-function dyn_graph=dynare_graph(dyn_graph,y,tit,x)
-% function dynare_graph(y,tit,x)
-% graphs
-%
-% INPUT
-%   figure_name: name of the figures
-%   colors: line colors
-%
-% OUTPUT
-%   dyn_graph: structure with figure information
-%
-% SPECIAL REQUIREMENT
-%   none
-
-% Copyright © 2006-2017 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/>.
-
-if nargin < 4
-    x = (1:size(y,1))';
-end
-nplot = dyn_graph.plot_nbr + 1;
-if nplot > dyn_graph.max_nplot
-    figure('Name',dyn_graph.figure_name);
-    nplot = 1;
-end
-dyn_graph.plot_nbr = nplot;
-subplot(dyn_graph.nr,dyn_graph.nc,nplot);
-
-line_types = dyn_graph.line_types;
-line_type = line_types{1};
-for i=1:size(y,2)
-    if length(line_types) > 1
-        line_type = line_types{i};
-    end
-
-    plot(x,y(:,i),line_type);
-    hold on
-end
-title(tit);
-hold off
--- a/matlab/dyntable.m
+++ b/matlab/dyntable.m
@ -63,7 +63,7 @@ end
 value_format  = sprintf('%%%d.%df', val_width, val_precis);
 header_string_format  = sprintf('%%%ds', val_width);

-if length(title) > 0
+if ~isempty(title)
    fprintf('\n\n%s\n',title);
 end

@ -72,7 +72,7 @@ if nargin==9
    disp(optional_header)
 end

-if length(headers) > 0
+if ~isempty(headers)
    hh_tbl = sprintf(label_format_leftbound , headers{1});
    for i=2:length(headers)
        hh_tbl  = [hh_tbl sprintf(header_string_format, headers{i})];
--- a/matlab/estimation/dsge_likelihood.m
+++ b/matlab/estimation/dsge_likelihood.m
@ -1,6 +1,6 @@
 function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,derivatives_info)
 % [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,oo_] = dsge_likelihood(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,oo_,derivatives_info)
-% Evaluates the posterior kernel of a DSGE model using the specified
+% Evaluates the negative of the posterior kernel of a DSGE model using the specified
 % kalman_algo; the resulting posterior includes the 2*pi constant of the
 % likelihood function
 %
@ -21,7 +21,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
 % - derivatives_info    [structure]     derivative info for identification
 %
 % OUTPUTS
-% - fval                    [double]        scalar, value of the likelihood or posterior kernel.
+% - fval                    [double]        scalar, value of minus the likelihood or posterior kernel.
 % - info                    [integer]       4×1 vector, informations resolution of the model and evaluation of the likelihood.
 % - exit_flag               [integer]       scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
 % - DLIK                    [double]        Vector with score of the likelihood
@ -37,7 +37,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
 % This function calls: dynare_resolve, lyapunov_symm, lyapunov_solver, compute_Pinf_Pstar, kalman_filter_d, missing_observations_kalman_filter_d,
 % univariate_kalman_filter_d, kalman_steady_state, get_perturbation_params_deriv, kalman_filter, missing_observations_kalman_filter, univariate_kalman_filter, priordens

-% Copyright © 2004-2023 Dynare Team
+% Copyright © 2004-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
--- a/matlab/forecast_graphs.m
+++ b/matlab/forecast_graphs.m
@ -76,7 +76,7 @@ for j= 1:nvar
            fprintf(fidTeX,' \n');
        end
        n_fig =n_fig+1;
-        eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+        hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts (' int2str(n_fig) ')']);
        m = 1;
    end
    subplot(nr, nc, m);
@ -138,7 +138,7 @@ if isfield(oo_.forecast,'HPDinf_ME')
                fprintf(fidTeX,' \n');
            end
            n_fig =n_fig+1;
-            eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+            hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts including ME (' int2str(n_fig) ')']);
            m = 1;
        end
        subplot(nr,nc,m);
--- a/matlab/latex/collect_latex_files.m
+++ b/matlab/latex/collect_latex_files.m
@ -35,7 +35,7 @@ fprintf(fid,'%s \n','\usepackage{psfrag}');
 fprintf(fid,'%s \n','\usepackage{graphicx}');
 fprintf(fid,'%s \n','\usepackage{epstopdf}');
 fprintf(fid,'%s \n','\usepackage{longtable,booktabs}');
-fprintf(fid,'%s \n','\usepackage{amsmath,amsfonts}');
+fprintf(fid,'%s \n','\usepackage{amsmath,amsfonts,amssymb}');
 fprintf(fid,'%s \n','\usepackage{breqn}');
 fprintf(fid,'%s \n','\usepackage{float,morefloats,caption}');
 fprintf(fid,'%s \n','\begin{document}');
--- a/matlab/model_info.m
+++ b/matlab/model_info.m
@ -280,20 +280,10 @@ function print_line(names,var_index,lead_lag,M_)
    else
        aux_index=find([M_.aux_vars(:).endo_index]==var_index);
        aux_type=M_.aux_vars(aux_index).type;
-        if ~isfield(M_.aux_vars(aux_index),'orig_lead_lag') || isempty(M_.aux_vars(aux_index).orig_lead_lag)
-            if ismember(aux_type,[1,3])
-                str = subst_auxvar(var_index, -1, M_);
-            elseif ismember(aux_type,[0,2])
-                str = subst_auxvar(var_index, 1, M_);
-            else
-                if lead_lag==0
-                    str = subst_auxvar(var_index, [], M_);
-                else
-                    str = subst_auxvar(var_index, lead_lag, M_);
-                end
-            end
+        if lead_lag==0
+            str = subst_auxvar(var_index, [], M_);
        else
-            str = subst_auxvar(var_index, M_.aux_vars(aux_index).orig_lead_lag, M_);
+            str = subst_auxvar(var_index, lead_lag, M_);
        end
        aux_orig_expression=M_.aux_vars(aux_index).orig_expr;
        if isempty(aux_orig_expression)
--- a/matlab/moments/compute_variance_decomposition.m
+++ b/matlab/moments/compute_variance_decomposition.m
@ -56,7 +56,7 @@ else
        var_decomp(stationary_vars,i) = vx2;
        variance_sum_loop = variance_sum_loop +vx2; %track overall variance over shocks
    end
-    if ~options_.pruning && max(abs(variance_sum_loop-var_stationary)./var_stationary) > 1e-4
+    if ~options_.pruning && max(abs(variance_sum_loop-var_stationary)./var_stationary) > 1e-4 && max(abs(variance_sum_loop-var_stationary))>1e-7
        warning(['Aggregate variance and sum of variances by shocks ' ...
            'differ by more than 0.01 %'])
    end
--- a/matlab/pac-tools/hVectors.m
+++ b/matlab/pac-tools/hVectors.m
@ -18,7 +18,7 @@ function [h, lrcp] = hVectors(params, H, auxmodel, kind, id)
 %    params(2:end-1) ⟶ Autoregressive parameters.
 %    params(end)     ⟶ Discount factor.

-% Copyright © 2018-2021 Dynare Team
+% Copyright © 2018-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -52,21 +52,21 @@ n = length(H);
 tmp = eye(n*m)-kron(G, transpose(H)); % inv(W2)

 switch kind
-  case 'll'
+  case 'll' % (A.84), page 28 in Brayton, Davis and Tulip (2000)
    h = A_1*A_b*((kron(iota(m, m), H))'*(tmp\kron(iota(m, m), iota(n, id))));
-  case 'dd'
+  case 'dd' % (A.79), page 26 in Brayton, Davis and Tulip (2000)
    h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), H')*(tmp\kron(iota(m, m), iota(n, id))));
-  case 'dl'
+  case 'dl' % (A.74), page 24 in Brayton, Davis and Tulip (2000)
    h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), (H'-eye(length(H))))*(tmp\kron(iota(m, m), iota(n, id))));
  otherwise
    error('Unknown kind value in PAC model.')
 end

-if nargin>1
+if nargout>1
    if isequal(kind, 'll')
        lrcp = NaN;
    else
        d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
        lrcp = (1-sum(params(2:end-1))-d);
    end
-end
+end
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
@ -14,7 +14,7 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl
 %   maxerror         [double]    ∞-norm of the residual
 %   per_block_status [struct]    vector structure with per-block information about convergence

-% Copyright © 2020-2023 Dynare Team
+% Copyright © 2020-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -33,18 +33,21 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl

 cutoff = 1e-15;

-if options_.stack_solve_algo==0
-    mthd='Sparse LU';
-elseif options_.stack_solve_algo==1 || options_.stack_solve_algo==6
-    mthd='LBJ';
-elseif options_.stack_solve_algo==2
-    mthd='GMRES';
-elseif options_.stack_solve_algo==3
-    mthd='BICGSTAB';
-elseif options_.stack_solve_algo==4
-    mthd='OPTIMPATH';
-else
-    mthd='UNKNOWN';
+switch options_.stack_solve_algo
+    case 0
+        mthd='Sparse LU on stacked system';
+    case {1,6}
+        mthd='LBJ with LU solver';
+    case 2
+        mthd='GMRES on stacked system';
+    case 3
+        mthd='BiCGStab on stacked system';
+    case 4
+        mthd='Sparse LU solver with optimal path length on stacked system';
+    case 7
+        mthd='Solver from solve_algo option on stacked system'
+    otherwise
+        error('Unsupported stack_solve_algo value')
 end
 if options_.verbosity
    printline(41)
--- a/matlab/send_endogenous_variables_to_workspace.m
+++ b/matlab/send_endogenous_variables_to_workspace.m
@ -20,5 +20,5 @@ function send_endogenous_variables_to_workspace()
 global M_ oo_

 for idx = 1:M_.endo_nbr
-    assignin('base', M_.endo_names{idx}, oo_.endo_simul(idx,:))
+    assignin('base', M_.endo_names{idx}, oo_.endo_simul(idx,:)')
 end
--- a/matlab/send_exogenous_variables_to_workspace.m
+++ b/matlab/send_exogenous_variables_to_workspace.m
@ -21,5 +21,5 @@ function send_exogenous_variables_to_workspace()
 global M_ oo_

 for idx = 1:M_.exo_nbr
-    assignin('base', M_.exo_names{idx}, oo_.exo_simul(:,idx))
+    assignin('base', M_.exo_names{idx}, oo_.exo_simul(:,idx)')
 end
--- a/matlab/stochastic_solver/k_order_pert.m
+++ b/matlab/stochastic_solver/k_order_pert.m
@ -1,7 +1,7 @@
 function [dr,info] = k_order_pert(dr,M_,options_)
 % Compute decision rules using the k-order DLL from Dynare++

-% Copyright © 2009-2023 Dynare Team
+% Copyright © 2009-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -30,6 +30,16 @@ if M_.maximum_endo_lead == 0 && order>1
           'backward models'])
 end

+if options_.aim_solver
+    error('Option aim_solver is not compatible with k_order_solver')
+end
+if options_.dr_cycle_reduction
+    error('Option dr=cycle_reduction is not compatible with k_order_solver')
+end
+if options_.dr_logarithmic_reduction
+    error('Option dr=logarithmic_reduction is not compatible with k_order_solver')
+end
+
 try
    [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M_,options_);
 catch ME
--- a/matlab/utilities/dataset/makedataset.m
+++ b/matlab/utilities/dataset/makedataset.m
@ -242,8 +242,14 @@ if newdatainterface
    end
 else
    % ...  or check that nobs is smaller than the number of observations in dataset_.
-    if nobs>dataset_.nobs
-        error('makedataset: nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(dataset_.nobs))
+    if FIRSTOBS>dataset_.dates(1)
+        if FIRSTOBS+nobs-1>dataset_.dates(end)
+            error('makedataset: Given first_obs=%u and %u total observations in the dataset, the current nobs of %s must not be greater than %s!', options_.first_obs, dataset_.nobs, num2str(nobs), num2str(dataset_.nobs-find(dataset_.dates==FIRSTOBS)+1))
+        end
+    else
+        if nobs>dataset_.nobs
+            error('makedataset: nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(dataset_.nobs))
+        end
    end
 end

--- a/matlab/utilities/general/display_parameter_values.m
+++ b/matlab/utilities/general/display_parameter_values.m
@ -1,17 +1,6 @@
-function dynare_graph_close()
-% function dynare_graph_close()
-% close a figure
-%
-% INPUT
-%   none
-%
-% OUTPUT
-%   none
-%
-% SPECIAL REQUIREMENT
-%   none
+function display_parameter_values

-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@ -27,3 +16,10 @@ function dynare_graph_close()
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+global M_
+my_title='Current parameter values:';
+labels = M_.param_names;
+headers = {'Parameter'; 'Value'};
+lh = cellofchararraymaxlength(labels)+2;
+options_.noprint=false;
+dyntable(options_, my_title, headers, labels, M_.params, lh, 10, 6);
--- a/meson.build
+++ b/meson.build
@ -21,9 +21,14 @@ c_compiler = meson.get_compiler('c')

 subdir('preprocessor/src')

-# Compatibility symlink
-# NB: the following my be improved if this wishlist item is done:
+# NB: the following two symlinks my be improved if this wishlist item is done:
 #  https://github.com/mesonbuild/meson/issues/11519
+
+install_symlink('dynare-preprocessor', install_dir : 'lib/dynare/preprocessor/',
+                pointing_to : '../../../bin/dynare-preprocessor'
+                              + (host_machine.system() == 'windows' ? '.exe' : ''))
+
+# Compatibility symlink
 install_symlink('dynare_m' + (host_machine.system() == 'windows' ? '.exe' : ''),
                install_dir : 'lib/dynare/matlab/preprocessor64',
                pointing_to : '../../../../bin/dynare-preprocessor'
@ -34,17 +39,12 @@ install_symlink('dynare_m' + (host_machine.system() == 'windows' ? '.exe' : ''),
 install_subdir('matlab', install_dir : 'lib/dynare',
               exclude_files : [ 'utilities/tests/.git' ,
                                 'utilities/tests/.gitignore',
-                                 'modules/reporting/.git',
-                                 'modules/reporting/.gitignore',
-                                 'modules/reporting/.gitlab-ci.yml',
-                                 'modules/dseries/.git',
-                                 'modules/dseries/.gitignore',
-                                 'modules/dseries/src/modules/matlab-fame-io/.git',
-                                 'modules/dseries/src/modules/matlab-fame-io/.gitignore',
-                                 'modules/dseries/src/modules/matlab-fame-io/doc/.gitignore',
-                                 'modules/dseries/src/modules/matlab-fame-io/tests/FameDatabases/.gitignore',
-                                 'modules/dseries/.gitmodules',
-                                 'modules/dseries/.gitlab-ci.yml' ])
+                                 'dseries/.git',
+                                 'dseries/.gitignore',
+                                 'dseries/src/modules/matlab-fame-io/.git',
+                                 'dseries/src/modules/matlab-fame-io/.gitignore',
+                                 'dseries/.gitmodules',
+                                 'dseries/.gitlab-ci.yml' ])

 sed_exe = find_program('sed')
 custom_target(output : 'dynare_version.m', input : 'matlab/dynare_version.m.in',
@ -53,6 +53,9 @@ custom_target(output : 'dynare_version.m', input : 'matlab/dynare_version.m.in',
              install : true,
              install_dir : 'lib/dynare/matlab')

+install_subdir('contrib/ms-sbvar/TZcode/MatlabFiles',
+               install_dir : 'lib/dynare/contrib/ms-sbvar/TZcode')
+
 ### MEX files

 mex_incdir = include_directories('mex/sources')
@ -157,7 +160,20 @@ if get_option('build_for') == 'matlab'
  umfpack_dep = declare_dependency(link_args : '-lmwumfpack', dependencies : blas_dep)
  ut_dep = declare_dependency(link_args : '-lut')

-  slicot_dep = declare_dependency(dependencies : [ fortran_compiler.find_library('slicot64_pic'), blas_dep, lapack_dep ])
+  # Workaround for Meson bug https://github.com/mesonbuild/meson/issues/12757
+  # Use the C compiler as a fallback for detecting SLICOT under Linux with
+  # prefer_static=true (but still try the Fortran compiler to honour the -B
+  # option in fortran_args, as documented). Needed for building the MATLAB
+  # Online package.
+  if get_option('prefer_static') and host_machine.system() == 'linux'
+    slicot_dep_tmp = fortran_compiler.find_library('slicot64_pic', required : false)
+    if not slicot_dep_tmp.found()
+      slicot_dep_tmp = c_compiler.find_library('slicot64_pic')
+    endif
+    slicot_dep = declare_dependency(dependencies : [ slicot_dep_tmp, blas_dep, lapack_dep ])
+  else
+    slicot_dep = declare_dependency(dependencies : [ fortran_compiler.find_library('slicot64_pic'), blas_dep, lapack_dep ])
+  endif
 else # Octave build
  octave_exe = find_program('octave', required : not meson.is_cross_build(), disabler : true)
  mkoctfile_exe = find_program('mkoctfile')
@ -291,8 +307,11 @@ if get_option('prefer_static')

  # NB: constructing a dependency object with link_args : ['-Wl,-Bstatic', '-lgomp', '-Wl,-Bdynamic'] does not work,
  # because it reorders the three arguments and puts -lgomp at the end
-  openmp_dep_tmp = cpp_compiler.find_library('gomp', static : true)
-  openmp_dep = declare_dependency(dependencies : [ openmp_dep, openmp_dep_tmp ])
+  if host_machine.system() != 'linux'
+    # Under Debian 12, trying to link (static) libgomp.a in a MEX fails.
+    openmp_dep_tmp = cpp_compiler.find_library('gomp', static : true)
+    openmp_dep = declare_dependency(dependencies : [ openmp_dep, openmp_dep_tmp ])
+  endif
 endif

 # For use when creating intermediate static libraries to be incorporated in MEX files
@ -1892,7 +1911,9 @@ endforeach
 git_exe = find_program('git', required : false)
 etags_exe = find_program('etags', required : false)

-if git_exe.found() and etags_exe.found()
+fs = import('fs')
+
+if fs.is_dir('.git') and git_exe.found() and etags_exe.found()
  all_files = run_command(git_exe,
                          [ '--git-dir=@0@/.git'.format(meson.project_source_root()),
                            'ls-files', '--recurse-submodules',
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -4710,23 +4710,26 @@ Interpreter::Simulate_Newton_Two_Boundaries(
          switch (stack_solve_algo)
            {
            case 0:
-              mexPrintf("MODEL SIMULATION: (method=Sparse LU)\n");
+              mexPrintf("MODEL SIMULATION: (method=Sparse LU solver on stacked system)\n");
              break;
            case 2:
-              mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GMRES)\n",
-                                                 preconditioner, false)
-                            .c_str());
+              mexPrintf(
+                  preconditioner_print_out("MODEL SIMULATION: (method=GMRES on stacked system)\n",
+                                           preconditioner, false)
+                      .c_str());
              break;
            case 3:
-              mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=BiCGStab)\n",
-                                                 preconditioner, false)
+              mexPrintf(preconditioner_print_out(
+                            "MODEL SIMULATION: (method=BiCGStab on stacked system)\n",
+                            preconditioner, false)
                            .c_str());
              break;
            case 4:
-              mexPrintf("MODEL SIMULATION: (method=Sparse LU & optimal path length)\n");
+              mexPrintf("MODEL SIMULATION: (method=Sparse LU solver with optimal path length on "
+                        "stacked system)\n");
              break;
            case 5:
-              mexPrintf("MODEL SIMULATION: (method=Sparse Gaussian Elimination)\n");
+              mexPrintf("MODEL SIMULATION: (method=LBJ with Sparse Gaussian Elimination)\n");
              break;
            }
        }
--- a/mex/sources/libkorder/sylv/QuasiTriangular.hh
+++ b/mex/sources/libkorder/sylv/QuasiTriangular.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -240,16 +240,11 @@ public:
    real = r;
  }
  virtual ~_matrix_iter() = default;
-  bool
+  [[nodiscard]] bool
  operator==(const _Self& it) const
  {
    return ptr == it.ptr;
  }
-  bool
-  operator!=(const _Self& it) const
-  {
-    return ptr != it.ptr;
-  }
  _TRef
  operator*() const
  {
--- a/mex/sources/libkorder/sylv/Vector.cc
+++ b/mex/sources/libkorder/sylv/Vector.cc
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -114,42 +114,6 @@ Vector::Vector(mxArray* p) :
    throw SYLV_MES_EXCEPTION("This is not a dense array of real doubles.");
 }

-bool
-Vector::operator==(const Vector& y) const
-{
-  return ConstVector(*this) == y;
-}
-
-bool
-Vector::operator!=(const Vector& y) const
-{
-  return ConstVector(*this) != y;
-}
-
-bool
-Vector::operator<(const Vector& y) const
-{
-  return ConstVector(*this) < y;
-}
-
-bool
-Vector::operator<=(const Vector& y) const
-{
-  return ConstVector(*this) <= y;
-}
-
-bool
-Vector::operator>(const Vector& y) const
-{
-  return ConstVector(*this) > y;
-}
-
-bool
-Vector::operator>=(const Vector& y) const
-{
-  return ConstVector(*this) >= y;
-}
-
 void
 Vector::zeros()
 {
@ -323,17 +287,10 @@ ConstVector::operator==(const ConstVector& y) const
  return i == len;
 }

-bool
-ConstVector::operator<(const ConstVector& y) const
+std::partial_ordering
+ConstVector::operator<=>(const ConstVector& y) const
 {
-  int i = std::min(len, y.len);
-  int ii = 0;
-  while (ii < i && operator[](ii) == y[ii])
-    ii++;
-  if (ii < i)
-    return operator[](ii) < y[ii];
-  else
-    return len < y.len;
+  return std::lexicographical_compare_three_way(data, data + len, y.data, y.data + y.len);
 }

 double
--- a/mex/sources/libkorder/sylv/Vector.hh
+++ b/mex/sources/libkorder/sylv/Vector.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -25,6 +25,7 @@
   to avoid running virtual method invokation mechanism. Some
   members, and methods are thus duplicated */

+#include <compare>
 #include <complex>
 #include <utility>

@ -108,15 +109,6 @@ public:
    return s;
  }

-  // Exact equality.
-  bool operator==(const Vector& y) const;
-  bool operator!=(const Vector& y) const;
-  // Lexicographic ordering.
-  bool operator<(const Vector& y) const;
-  bool operator<=(const Vector& y) const;
-  bool operator>(const Vector& y) const;
-  bool operator>=(const Vector& y) const;
-
  virtual ~Vector()
  {
    if (destroy)
@ -219,29 +211,9 @@ public:
    return s;
  }
  // Exact equality
-  bool operator==(const ConstVector& y) const;
-  bool
-  operator!=(const ConstVector& y) const
-  {
-    return !operator==(y);
-  }
+  [[nodiscard]] bool operator==(const ConstVector& y) const;
  // Lexicographic ordering
-  bool operator<(const ConstVector& y) const;
-  bool
-  operator<=(const ConstVector& y) const
-  {
-    return operator<(y) || operator==(y);
-  }
-  bool
-  operator>(const ConstVector& y) const
-  {
-    return !operator<=(y);
-  }
-  bool
-  operator>=(const ConstVector& y) const
-  {
-    return !operator<(y);
-  }
+  [[nodiscard]] std::partial_ordering operator<=>(const ConstVector& y) const;

  [[nodiscard]] double getNorm() const;
  [[nodiscard]] double getMax() const;
--- a/mex/sources/libkorder/tl/equivalence.cc
+++ b/mex/sources/libkorder/tl/equivalence.cc
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -35,12 +35,14 @@ OrdSequence::operator[](int i) const
   orderings can be used for different problem sizes. We order them
   according to the average, and then according to the first item. */

-bool
-OrdSequence::operator<(const OrdSequence& s) const
+std::partial_ordering
+OrdSequence::operator<=>(const OrdSequence& s) const
 {
  double ta = average();
  double sa = s.average();
-  return (ta < sa || ((ta == sa) && (operator[](0) > s[0])));
+  if (auto cmp1 = ta <=> sa; cmp1 != 0)
+    return cmp1;
+  return operator[](0) <=> s[0];
 }

 bool
--- a/mex/sources/libkorder/tl/equivalence.hh
+++ b/mex/sources/libkorder/tl/equivalence.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -55,6 +55,7 @@

 #include "int_sequence.hh"

+#include <compare>
 #include <list>
 #include <string>
 #include <vector>
@ -62,7 +63,7 @@
 /* Here is the abstraction for an equivalence class. We implement it as
   vector<int>. We have a constructor for empty class, copy
   constructor. What is important here is the ordering operator
-   operator<() and methods for addition of an integer, and addition of
+   operator<=>() and methods for addition of an integer, and addition of
   another sequence. Also we provide method has() which returns true if a
   given integer is contained. */

@ -74,9 +75,9 @@ public:
  OrdSequence() : data()
  {
  }
-  bool operator==(const OrdSequence& s) const;
+  [[nodiscard]] bool operator==(const OrdSequence& s) const;
  int operator[](int i) const;
-  bool operator<(const OrdSequence& s) const;
+  [[nodiscard]] std::partial_ordering operator<=>(const OrdSequence& s) const;
  [[nodiscard]] const std::vector<int>&
  getData() const
  {
@ -120,12 +121,7 @@ public:
  // Copy constructor plus gluing i1 and i2 in one class
  Equivalence(const Equivalence& e, int i1, int i2);

-  bool operator==(const Equivalence& e) const;
-  bool
-  operator!=(const Equivalence& e) const
-  {
-    return !operator==(e);
-  }
+  [[nodiscard]] bool operator==(const Equivalence& e) const;
  [[nodiscard]] int
  getN() const
  {
--- a/mex/sources/libkorder/tl/gs_tensor.hh
+++ b/mex/sources/libkorder/tl/gs_tensor.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -76,16 +76,11 @@ public:
  // Constructs the tensor dimensions for slicing (see the implementation for details)
  TensorDimens(const IntSequence& ss, const IntSequence& coor);

-  bool
+  [[nodiscard]] bool
  operator==(const TensorDimens& td) const
  {
    return nvs == td.nvs && sym == td.sym;
  }
-  bool
-  operator!=(const TensorDimens& td) const
-  {
-    return !operator==(td);
-  }

  [[nodiscard]] int
  dimen() const
--- a/mex/sources/libkorder/tl/int_sequence.cc
+++ b/mex/sources/libkorder/tl/int_sequence.cc
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -103,10 +103,10 @@ IntSequence::operator==(const IntSequence& s) const
  return std::equal(data, data + length, s.data, s.data + s.length);
 }

-bool
-IntSequence::operator<(const IntSequence& s) const
+std::strong_ordering
+IntSequence::operator<=>(const IntSequence& s) const
 {
-  return std::lexicographical_compare(data, data + length, s.data, s.data + s.length);
+  return std::lexicographical_compare_three_way(data, data + length, s.data, s.data + s.length);
 }

 bool
--- a/mex/sources/libkorder/tl/int_sequence.hh
+++ b/mex/sources/libkorder/tl/int_sequence.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -43,6 +43,7 @@
 #define INT_SEQUENCE_HH

 #include <algorithm>
+#include <compare>
 #include <initializer_list>
 #include <utility>
 #include <vector>
@ -119,12 +120,7 @@ public:
    if (destroy)
      delete[] data;
  }
-  bool operator==(const IntSequence& s) const;
-  bool
-  operator!=(const IntSequence& s) const
-  {
-    return !operator==(s);
-  }
+  [[nodiscard]] bool operator==(const IntSequence& s) const;
  int&
  operator[](int i)
  {
@ -141,15 +137,10 @@ public:
    return length;
  }

-  /* We provide two orderings. The first operator<() is the linear
+  /* We provide two orderings. The first operator<=>() is the linear
     lexicographic ordering, the second less() is the non-linear Cartesian
     ordering. */
-  bool operator<(const IntSequence& s) const;
-  bool
-  operator<=(const IntSequence& s) const
-  {
-    return (operator==(s) || operator<(s));
-  }
+  [[nodiscard]] std::strong_ordering operator<=>(const IntSequence& s) const;
  [[nodiscard]] bool lessEq(const IntSequence& s) const;
  [[nodiscard]] bool less(const IntSequence& s) const;

--- a/mex/sources/libkorder/tl/kron_prod.hh
+++ b/mex/sources/libkorder/tl/kron_prod.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -86,7 +86,7 @@ public:

  KronProdDimens& operator=(const KronProdDimens& kd) = default;
  KronProdDimens& operator=(KronProdDimens&& kd) = default;
-  bool
+  [[nodiscard]] bool
  operator==(const KronProdDimens& kd) const
  {
    return rows == kd.rows && cols == kd.cols;
--- a/mex/sources/libkorder/tl/permutation.hh
+++ b/mex/sources/libkorder/tl/permutation.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -95,7 +95,7 @@ public:
  Permutation(const Permutation& p, int i) : permap(p.permap.insert(p.size(), i))
  {
  }
-  bool
+  [[nodiscard]] bool
  operator==(const Permutation& p) const
  {
    return permap == p.permap;
--- a/mex/sources/libkorder/tl/ps_tensor.hh
+++ b/mex/sources/libkorder/tl/ps_tensor.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -115,7 +115,7 @@ public:
  {
    per.apply(nvmax);
  }
-  bool
+  [[nodiscard]] bool
  operator==(const PerTensorDimens& td) const
  {
    return TensorDimens::operator==(td) && per == td.per;
--- a/mex/sources/libkorder/tl/symmetry.hh
+++ b/mex/sources/libkorder/tl/symmetry.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -126,16 +126,11 @@ public:
  {
    return run;
  }
-  bool
-  operator=(const symiterator& it)
+  [[nodiscard]] bool
+  operator==(const symiterator& it)
  {
    return dim == it.dim && run == it.run;
  }
-  bool
-  operator!=(const symiterator& it)
-  {
-    return !operator=(it);
-  }
 };

 /* The class SymmetrySet defines a set of symmetries of the given length
--- a/mex/sources/libkorder/tl/tensor.hh
+++ b/mex/sources/libkorder/tl/tensor.hh
@ -1,6 +1,6 @@
 /*
 * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@ -150,16 +150,11 @@ public:
    {
      return offset;
    }
-    bool
+    [[nodiscard]] bool
    operator==(const index& n) const
    {
      return offset == n.offset;
    }
-    bool
-    operator!=(const index& n) const
-    {
-      return offset != n.offset;
-    }
    [[nodiscard]] const IntSequence&
    getCoor() const
    {
--- a/scripts/matlab-online/packageDynare.sh
+++ b/scripts/matlab-online/packageDynare.sh
@ -4,12 +4,8 @@ set -exo pipefail
 # Creates a dynare-X.Y.mltbx in the current repository, using the settings below.
 # Needs to be run from Ubuntu 22.04 LTS, with the needed packages installed.

-DYNAREVER=5.5
 X13ASVER=1-1-b60
-LIBLOCATION="/MATLAB Add-Ons/Toolboxes/Dynare/mex/matlab/libs"
 MATLABPATH=/opt/MATLAB/R2023b
-# MatIO has been recompiled by hand to avoid the dependency on HDF5, which is a nightmare
-MATIO_PREFIX=/home/sebastien/usr
 # TODO: change size and put white background for better rendering in MATLAB Add-Ons browser
 DYNARE_PNG_LOGO=../../preprocessor/doc/logos/dlogo.png

@ -21,48 +17,31 @@ cleanup ()
 }
 trap cleanup EXIT

-pushd "$tmpdir"
+pushd ../..
+meson setup -Dmatlab_path="$MATLABPATH" -Dbuildtype=release -Dprefer_static=true "$tmpdir"/build-matlab-online

-# Get Dynare
-wget -q https://www.dynare.org/release/source/dynare-$DYNAREVER.tar.xz
-tar xf dynare-$DYNAREVER.tar.xz
+cd "$tmpdir"/build-matlab-online
+meson compile
+meson install --destdir "$tmpdir"
+DYNAREVER=$(meson introspect --projectinfo | jq -r '.version')

-# Build Dynare
-cd dynare-$DYNAREVER
-
-# Use static libstdc++ otherwise the one shipped with MATLAB creates problem (since it overrides the systemd-wide one from Ubuntu)
-./configure --with-matlab="$MATLABPATH" --with-matio="$MATIO_PREFIX" --disable-octave --disable-doc --disable-dynare++ LDFLAGS=-static-libstdc++
-make -j$(nproc)
-
-strip preprocessor/dynare-preprocessor
-strip mex/matlab/*.mexa64
-
-# Patch mex files to look into ./lib folder
-for f in mex/matlab/*.mexa64
-do
-    patchelf --set-rpath "$LIBLOCATION" $f
-done
-
-# Grab the shared libraries needed
-mkdir mex/matlab/libs
-for l in libgsl.so.27 libgslcblas.so.0
-do
-    cp /usr/lib/x86_64-linux-gnu/$l mex/matlab/libs/
-    # Patch rpath to find dependencies at runtime
-    patchelf --set-rpath "$LIBLOCATION" mex/matlab/libs/$l
-done
-cp "$MATIO_PREFIX"/lib/libmatio.so.11 mex/matlab/libs/
-patchelf --set-rpath "$LIBLOCATION" mex/matlab/libs/libmatio.so.11 # Probably not needed
+cd ..
+strip usr/local/bin/dynare-preprocessor
+strip usr/local/lib/dynare/mex/matlab/*.mexa64

 # Get X13 binary from the Census Bureau website
 # The binary from Ubuntu has some shared library dependencies, so it is safer to use a static binary
 wget -q https://www2.census.gov/software/x-13arima-seats/x13as/unix-linux/program-archives/x13as_ascii-v${X13ASVER}.tar.gz
 tar xf x13as_ascii-v${X13ASVER}.tar.gz
-mkdir -p matlab/modules/dseries/externals/x13/linux/64
-cp x13as/x13as_ascii matlab/modules/dseries/externals/x13/linux/64/x13as
+
+# Populate staging area for the zip
+cp -pRL usr/local/lib/dynare dynare # -L is needed to dereference the preprocessor symlink
+mkdir -p dynare/matlab/dseries/externals/x13/linux/64
+cp -p x13as/x13as_ascii dynare/matlab/dseries/externals/x13/linux/64/x13as

 # zip dynare
-zip -q -r "$tmpdir"/dynare.zip contrib/jsonlab contrib/ms-sbvar/TZcode/MatlabFiles examples matlab mex/matlab preprocessor/dynare-preprocessor license.txt NEWS.md README.md VERSION
+cd dynare
+zip -q -r "$tmpdir"/dynare.zip *

 # make toolbox
 popd
--- a/tests/deterministic_simulations/purely_backward/ar1.mod
+++ b/tests/deterministic_simulations/purely_backward/ar1.mod
@ -33,6 +33,6 @@ if ~oo_.deterministic_simulation.status
 end
 send_endogenous_variables_to_workspace;

-if max(abs(y'-[1; exp(cumprod([1; rho*ones(9, 1)]))]))>options_.dynatol.x
+if max(abs(y-[1; exp(cumprod([1; rho*ones(9, 1)]))]))>options_.dynatol.x
    error('Wrong solution!')
 end
--- a/tests/moments/example1_hp_test.mod
+++ b/tests/moments/example1_hp_test.mod
@ -77,12 +77,12 @@ send_endogenous_variables_to_workspace;
 options_.nomoments=0;
 oo_unfiltered_all_shocks=oo_;

-[junk, y_filtered]=sample_hp_filter(y',1600);
-[junk, c_filtered]=sample_hp_filter(c',1600);
-[junk, k_filtered]=sample_hp_filter(k',1600);
-[junk, a_filtered]=sample_hp_filter(a',1600);
-[junk, h_filtered]=sample_hp_filter(h',1600);
-[junk, b_filtered]=sample_hp_filter(b',1600);
+[junk, y_filtered]=sample_hp_filter(y,1600);
+[junk, c_filtered]=sample_hp_filter(c,1600);
+[junk, k_filtered]=sample_hp_filter(k,1600);
+[junk, a_filtered]=sample_hp_filter(a,1600);
+[junk, h_filtered]=sample_hp_filter(h,1600);
+[junk, b_filtered]=sample_hp_filter(b,1600);

 verbatim;
 total_std_all_shocks_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]);
@ -112,12 +112,12 @@ stoch_simul(order=1,nofunctions,hp_filter=0,periods=2500000,nomoments);
 send_endogenous_variables_to_workspace;
 oo_unfiltered_one_shock=oo_;

-[junk, y_filtered]=sample_hp_filter(y',1600);
-[junk, c_filtered]=sample_hp_filter(c',1600);
-[junk, k_filtered]=sample_hp_filter(k',1600);
-[junk, a_filtered]=sample_hp_filter(a',1600);
-[junk, h_filtered]=sample_hp_filter(h',1600);
-[junk, b_filtered]=sample_hp_filter(b',1600);
+[junk, y_filtered]=sample_hp_filter(y,1600);
+[junk, c_filtered]=sample_hp_filter(c,1600);
+[junk, k_filtered]=sample_hp_filter(k,1600);
+[junk, a_filtered]=sample_hp_filter(a,1600);
+[junk, h_filtered]=sample_hp_filter(h,1600);
+[junk, b_filtered]=sample_hp_filter(b,1600);

 verbatim;
 total_std_one_shock_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]);
--- a/tests/occbin/filter/NKM.mod
+++ b/tests/occbin/filter/NKM.mod
@ -13,66 +13,66 @@

 // ----------------- Defintions -----------------------------------------//
 var        
-    c          //1  Consumption 
-    n          //2  Labor
-    y          //5  Output
-    yf         //6  Final goods       
-    yg         //11 Output growth gap
-    w          //12 Real wage rate
-    wf         //13 Flexible real wage
-    pigap      //15 Inflation rate -> pi(t)/pibar = pigap
-    inom    ${i^{nom}}$       //16 Nominal interest rate
-    inomnot    //17 Notional interest rate
-    mc         //19 Real marginal cost
-    lam  ${\lambda}$      //20 Inverse marginal utility of wealth  
-    g          //21 Growth shock       
-    s          //22 Risk premium shock
-    mp         //23 Monetary policy shock    
-    pi   ${\pi}$   //24 Observed inflation    
+    c       $c$         (long_name='Consumption')
+    n                   (long_name='Labor')
+    y       $y$         (long_name='Output')
+    yf                  (long_name='Final goods')
+    yg                  (long_name='Output growth gap')
+    w                   (long_name='Real wage rate')
+    wf                  (long_name='Flexible real wage')
+    pigap               (long_name='Inflation rate -> pi(t)/pibar = pigap')
+    inom    ${i^{nom}}$ (long_name='Nominal interest rate')
+    inomnot             (long_name='Notional interest rate')
+    mc                  (long_name='Real marginal cost')
+    lam  ${\lambda}$    (long_name='Inverse marginal utility of wealth')
+    g                   (long_name='Growth shock')
+    s                   (long_name='Risk premium shock')
+    mp                  (long_name='Monetary policy shock')
+    pi   ${\pi}$        (long_name='Observed inflation')
    @#if !(small_model)
-        x          //3  Investment
-        k          //4  Capital    
-        u          //7  Utilization cost
-        ups        //8  Utilization choice
-        wg         //9  Real wage growth gap    
-        xg         //10 Investment growth
-        rk         //14 Real rental rate
-        q          //18 Tobins q   
+        x               (long_name='Investment')
+        k               (long_name='Capital')
+        u               (long_name='Utilization cost')
+        ups             (long_name='Utilization choice')
+        wg              (long_name='Real wage growth gap')    
+        xg              (long_name='Investment growth')
+        rk              (long_name='Real rental rate')
+        q               (long_name='Tobins q')
    @#endif
 ;
 varexo          
-    epsg    ${\varepsilon_g}$   // Productivity growth shock
-    epsi       // Notional interest rate shock
-    epss       // Risk premium shock
+    epsg    ${\varepsilon_g}$   (long_name='Productivity growth shock')
+    epsi                        (long_name='Notional interest rate shock')
+    epss                        (long_name='Risk premium shock')
 ;        
 parameters
    // Calibrated Parameters    
-    beta    $\beta$    // Discount factor
-    chi         // Labor disutility scale
-    thetap      // Elasticity of subs. between intermediate goods
-    thetaw      // Elasticity of subs. between labor types
-    nbar        // Steady state labor
-    eta         // Inverse frish elasticity of labor supply
-    delta       // Depreciation
-    alpha       // Capital share
-    gbar        // Mean growth rate
-    pibar       // Inflation target
-    inombar     // Steady gross nom interest rate    
-    inomlb      // Effective lower bound on gross nominal interest rate    
-    sbar        // Average risk premium
+    beta    $\beta$    (long_name='Discount factor')
+    chi                 (long_name='Labor disutility scale')
+    thetap              (long_name='Elasticity of subs. between intermediate goods')
+    thetaw              (long_name='Elasticity of subs. between labor types')
+    nbar                (long_name='Steady state labor')
+    eta                 (long_name='Inverse frish elasticity of labor supply')
+    delta               (long_name='Depreciation')
+    alpha               (long_name='Capital share')
+    gbar                (long_name='Mean growth rate')
+    pibar               (long_name='Inflation target')
+    inombar             (long_name='Steady gross nom interest rate')
+    inomlb              (long_name='Effective lower bound on gross nominal interest rate')
+    sbar                (long_name='Average risk premium')
    // Parameters for DGP and Estimated parameters
-    varphip     // Rotemberg price adjustment cost
-    varphiw     // Rotemberg wage adjustment cost
-    h          	// Habit persistence
-    rhos        // Persistence
-    rhoi    	// Persistence
-    sigz        // Standard deviation technology
-    sigs        // Standard deviation risk premia
-    sigi        // Standard deviation mon pol
-    phipi       // Inflation responsiveness
-    phiy        // Output responsiveness
-    nu          // Investment adjustment cost
-    sigups      // Utilization    
+    varphip             (long_name='Rotemberg price adjustment cost')
+    varphiw             (long_name='Rotemberg wage adjustment cost')
+    h          	        (long_name='Habit persistence')
+    rhos                (long_name='Persistence')
+    rhoi    	        (long_name='Persistence')
+    sigz                (long_name='Standard deviation technology')
+    sigs                (long_name='Standard deviation risk premia')
+    sigi                (long_name='Standard deviation mon pol')
+    phipi               (long_name='Inflation responsiveness')
+    phiy                (long_name='Output responsiveness')
+    nu                  (long_name='Investment adjustment cost')
+    sigups              (long_name='Utilization')
 ;


@ -316,36 +316,34 @@ varobs yg inom pi;

    // forecast starting from period 42, zero shocks (default)
    smoother2histval(period=42);
-    [oo, error_flag] = occbin.forecast(options_,M_,oo_,8);
+    [forecast, error_flag] = occbin.forecast(options_,M_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,8);
    // forecast with stochastic shocks
    options_.occbin.forecast.qmc=true;
    options_.occbin.forecast.replic=127;
-    [oo1, error_flag] = occbin.forecast(options_,M_,oo_,8);
+    [forecast1, error_flag] = occbin.forecast(options_,M_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,8);
+
+    figure('Name','OccBin: Forecasts')
+    subplot(2,1,1)
+    plot(1:8,forecast.piecewise.Mean.inom,'b-',1:8,forecast1.piecewise.Mean.inom,'r--')
+    subplot(2,1,2)
+    plot(1:8,forecast.piecewise.Mean.y,'b-',1:8,forecast1.piecewise.Mean.y,'r--')

    // GIRF given states in 42 and shocks in 43
    t0=42;
    options_.occbin.irf.exo_names=M_.exo_names;
    options_.occbin.irf.t0=t0;
-    oo_ = occbin.irf(M_,oo_,options_);
+    oo_.occbin.irfs = occbin.irf(M_,oo_,options_);

-    vars_irf = {
-    'c', 'consumption'    
-    'n', 'labor'    
-    'y', 'output'    
-    'pigap', 'inflation rate'
-    'inom', 'interest rate'  
-    'inomnot', 'shadow rate'
-    };
+    var_list_ = {'c','n','y','pigap','inom','inomnot'};

-    options_.occbin.plot_irf.exo_names = M_.exo_names;
-    options_.occbin.plot_irf.endo_names = vars_irf(:,1);
-    options_.occbin.plot_irf.endo_names_long = vars_irf(:,2);
    // if you want to scale ...
    // options_occbin_.plot_irf.endo_scaling_factor = vars_irf(:,3);
    options_.occbin.plot_irf.simulname = ['t0_' int2str(t0)];
-    options_.occbin.plot_irf.tplot = min(40,options_.irf);
-    occbin.plot_irfs(M_,oo_,options_);
-
+    options_.irf=40;
+    occbin.plot_irfs(M_,oo_.occbin.irfs,options_,var_list_);
+    var_list_={};
+    options_.occbin.plot_irf.simulname = ['t0_' int2str(t0) '_full'];
+    occbin.plot_irfs(M_,oo_.occbin.irfs,options_,var_list_);
    oo0=oo_;
    // use smoother_redux
    estimation(
@ -374,7 +372,7 @@ varobs yg inom pi;
            consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty);
            
    // show initial condition effect of IF
-    figure,
+    figure('Name','OccBin: Smoothed shocks')
    subplot(221)
    plot([oo0.SmoothedShocks.epsg oo_.SmoothedShocks.epsg]), title('epsg')
    subplot(222)
@ -382,7 +380,7 @@ varobs yg inom pi;
    subplot(223)
    plot([oo0.SmoothedShocks.epss oo_.SmoothedShocks.epss]), title('epss')
    legend('PKF','IF')
-    figure,
+    figure('Name','OccBin: Smoothed Variables')
    subplot(221)
    plot([oo0.SmoothedVariables.inom oo_.SmoothedVariables.inom]), title('inom')
    subplot(222)
--- a/windows/deps/Makefile
+++ b/windows/deps/Makefile
@ -157,6 +157,7 @@ mingw64: tarballs/mingw-w64-x86_64-gcc-$(MINGW64_GCC_VERSION)-any.pkg.tar.zst ta
 	touch $@

 tarballs/mingw-w64-x86_64-%-any.pkg.tar.zst:
+	mkdir -p tarballs
 	wget $(WGET_OPTIONS) -O $@ http://repo.msys2.org/mingw/x86_64/$(notdir $@)

 clean-msys2:
--- a/windows/dynare.nsi
+++ b/windows/dynare.nsi
@ -52,7 +52,7 @@ Section "Dynare core (preprocessor and M-files)"
 File README.txt ..\NEWS.md ..\license.txt

 SetOutPath $INSTDIR\matlab
- File /r ..\matlab\*.m ..\build-win-matlab\dynare_version.m
+ File /r ..\matlab\*.m ..\matlab\*.json ..\build-win-matlab\dynare_version.m

 SetOutPath $INSTDIR\preprocessor
 File ..\build-win-matlab\preprocessor\src\dynare-preprocessor.exe
Author	SHA1	Message	Date
Stéphane Adjemian (Argos)	c33b53b045	Bug fix (wrong condition) We should test the numbe of output arguments, not the number of input arguments. This was bug was probably not affecting the outcomes since the number of input arguments is always greater than 1.	2024-01-22 16:54:02 +01:00
Stéphane Adjemian (Argos)	50beb8000d	Cosmetic changes.	2024-01-22 16:54:01 +01:00
Sébastien Villemot	881f5f2e62	Merge branch 'dynare-gsa_qz' Ref. !2270	2024-01-22 16:47:32 +01:00
Johannes Pfeifer	b1cb309a73	gsa: add proper check for correctness of qz_criterium with unit roots Critical for stability mapping	2024-01-22 16:47:23 +01:00
Sébastien Villemot	330542a79f	Build system: workaround for Meson bug, needed for building MATLAB Online package	2024-01-22 16:00:39 +01:00
Sébastien Villemot	b2f603091a	MATLAB Online package: update build script for meson Also, use the local git checkout instead of downloading a source tarball from the website.	2024-01-19 20:40:56 +01:00
Sébastien Villemot	88236b1cc0	Build system: under Linux, do not try to statically link libgomp even with -Dprefer_static=true Under Debian 12, it fails with: /usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/12/libgomp.a(team.o): réadressage R_X86_64_TPOFF32 vers symbole caché « gomp_tls_data » ne peut pas être utilisé en créant un objet partagé	2024-01-19 16:22:25 +01:00
Sébastien Villemot	c14d410699	Merge branch 'gsa_recursive' into 'master' GSA_recursive: make sure nobs is correctly set before checking for recursive estimation Closes #1611 See merge request Dynare/dynare!2269	2024-01-17 20:05:34 +00:00
Johannes Pfeifer	619de017d6	GSA_recursive: make sure nobs is correctly set before checking for recursive estimation Closes #1611	2024-01-17 17:35:50 +01:00
Sébastien Villemot	85c637d3d1	Manual: do not add parentheses to synopsis of functions without arguments Closes: #1707	2024-01-17 16:31:41 +01:00
Sébastien Villemot	6fe43545d8	Remove remnants of workaround for incorrect display of macro-directives without arguments This workaround was implemented in `cd195279e9`. The zero-width spaces were inadvertently removed in `248ad18846`. Ref. #1707	2024-01-17 16:31:40 +01:00
Sébastien Villemot	6b6b648255	Manual: minor fixes to function synopses	2024-01-17 16:31:37 +01:00
Sébastien Villemot	aa76e06099	Merge branch 'param_table' of git.dynare.org:JohannesPfeifer/dynare Ref. !2268	2024-01-17 12:09:01 +01:00
Sébastien Villemot	e68be9ec55	Merge branch 'latex_doc' of git.dynare.org:JohannesPfeifer/dynare Ref. !2267	2024-01-17 12:08:22 +01:00
Sébastien Villemot	2c6df2b668	Merge branch 'dsge_likelihood' of git.dynare.org:JohannesPfeifer/dynare Ref. !2266	2024-01-17 12:08:00 +01:00
Sébastien Villemot	c6c179d8dd	Merge branch 'uifigures' of git.dynare.org:JohannesPfeifer/dynare Ref. !2263	2024-01-17 12:06:19 +01:00
Sébastien Villemot	aceebe9769	Merge branch 'variance_warning' of git.dynare.org:JohannesPfeifer/dynare Ref. !2261	2024-01-17 12:05:22 +01:00
Johannes Pfeifer	28df34df06	Add display_parameter_values.m utility	2024-01-16 17:25:09 +01:00
Johannes Pfeifer	d7242c056d	dyntable.m: cosmetic fix	2024-01-16 17:24:44 +01:00
Johannes Pfeifer	6ba9cb06e7	manual: add missing package from `e5951650`	2024-01-16 16:13:14 +01:00
Johannes Pfeifer	bcd57142c6	dsge_likelihood.m: clarify that output is minus the posterior	2024-01-16 16:03:29 +01:00
Sébastien Villemot	6868e2f27e	Merge branch 'occbin_tools' into 'master' OccBin tools: rework codes See merge request Dynare/dynare!2251	2024-01-12 19:09:39 +00:00
Johannes Pfeifer	248d8ae84f	occbin.DSGE_smoother.m: correct figure caption	2024-01-12 18:16:30 +01:00
Sébastien Villemot	6780b9a009	Merge branch 'collect_latex' of git.dynare.org:JohannesPfeifer/dynare Ref. !2265	2024-01-12 17:50:35 +01:00
Sébastien Villemot	6235de3a82	Merge branch 'model_info' of git.dynare.org:JohannesPfeifer/dynare Ref. !2264	2024-01-12 17:50:07 +01:00
Sébastien Villemot	e1833c475d	Merge branch 'gsa' of git.dynare.org:JohannesPfeifer/dynare Ref. !2262	2024-01-12 17:48:54 +01:00
Johannes Pfeifer	e5951650d9	collect_latex_files.m: add amssymb package required for some symbols users regularly use	2024-01-12 11:22:46 +01:00
Johannes Pfeifer	1b181fca57	OccBin tools: rework codes	2024-01-10 14:42:21 +01:00
Johannes Pfeifer	45e8ab14dc	🐛 model_info.m: fix display of lagged states preprocessor increments lags/leads always only by lead_lag	2024-01-09 17:12:04 +01:00
Johannes Pfeifer	6bb973a0fa	map_identification.m: remove unsupported morris=0 case	2024-01-09 16:03:33 +01:00
Johannes Pfeifer	23225aca1b	simulated_moment_uncertainty.m: remove evalin following removal of globals transformed steady state is not passed back anymore	2024-01-09 16:03:10 +01:00
Sébastien Villemot	2d90ef9890	Merge branch 'dynare_m' of git.dynare.org:JohannesPfeifer/dynare Ref. !2256	2024-01-09 15:33:37 +01:00
Sébastien Villemot	5d169d658e	Manual: update Normann’s affiliation	2024-01-09 13:08:17 +01:00
Sébastien Villemot	4683b9aeb6	C++17 modernization: add [[nodiscard]] to operator==()	2024-01-09 13:08:17 +01:00
Sébastien Villemot	cb7bd0778b	C++20 modernization: drop comparison operators made obsolete by rewriting rules The implicit conversion from Vector to ConstVector, in addition to the rewriting rules and ConstVector::operator<=>, ensure that those are no longer needed.	2024-01-09 13:08:17 +01:00
Sébastien Villemot	2fd4a6fac4	C++20 modernization: use new spaceship operator to simplify comparisons	2024-01-09 13:08:15 +01:00
Sébastien Villemot	2bfa311636	C++20 modernization: remove trivial != operator definitions Those are no longer needed thanks to expression rewriting.	2024-01-09 12:12:38 +01:00
Sébastien Villemot	3ab9be71c8	🐛 k-order MEX: equality operator on symiterator had wrong name	2024-01-09 12:12:25 +01:00
Johannes Pfeifer	0187ebe0a2	Tag Dynare figures and add utility for moving figures to uitabgroup	2024-01-08 14:54:15 +01:00
Johannes Pfeifer	6ff924550c	compute_variance_decomposition.m: only print warning if absolute difference is meaningful Prevents warnings if relative difference involves division by almost 0	2024-01-07 10:24:13 +01:00
Sébastien Villemot	346ee01107	Merge branch 'beta_fixes' of git.dynare.org:JohannesPfeifer/dynare Ref. !2260	2024-01-05 16:12:25 +01:00
Sébastien Villemot	cc02690acf	Forbid alternative 1st order solvers with k_order_solver option	2024-01-05 16:08:09 +01:00
Sébastien Villemot	ffb578276e	Improve naming and description of various stack_solve_algo values Also minor improvements to solve_algo description.	2024-01-05 16:03:19 +01:00
Johannes Pfeifer	4a7851b069	windows installer: add json-file	2024-01-05 11:50:50 +01:00
Johannes Pfeifer	46d7e155d9	send_endogenous_variables_to_workspace.m and friends: output column instead of row vectors	2024-01-05 11:47:35 +01:00
Sébastien Villemot	f9cd465fea	Build system: don’t try to create TAGS file when not in a git working directory	2024-01-03 18:08:51 +01:00
Sébastien Villemot	53d8278d8a	license.txt: fix various issues detected by lintian	2024-01-03 17:22:49 +01:00
Sébastien Villemot	049006a1bf	Build system: install preprocessor symlink under libdir	2024-01-03 16:14:29 +01:00
Sébastien Villemot	e7cd6eb408	Build system: install .m files for MS-SBVAR	2024-01-03 15:58:28 +01:00
Sébastien Villemot	0679da4cba	Build system: update list of ignored files under matlab/	2024-01-03 15:53:38 +01:00
Sébastien Villemot	a99beac083	Windows package: add missing rule for creating tarballs directory [skip ci]	2024-01-03 10:23:53 +01:00
Sébastien Villemot	44044904c1	Merge branch 'forecast_graph' of git.dynare.org:JohannesPfeifer/dynare Ref. !2259	2024-01-03 10:12:57 +01:00
Sébastien Villemot	b89e79862a	Merge branch 'bvar' of git.dynare.org:JohannesPfeifer/dynare Ref. !2258	2024-01-03 10:11:33 +01:00
Sébastien Villemot	3249d74220	Merge branch 'first_obs' of git.dynare.org:JohannesPfeifer/dynare Ref. !2257	2024-01-03 10:10:47 +01:00
Sébastien Villemot	8a7440c6ac	Update copyright years	2024-01-03 10:02:21 +01:00
Johannes Pfeifer	90266c0bfc	gsa: decrease waitbar refresh frequency	2023-12-28 14:52:08 +01:00
Johannes Pfeifer	02d1e8d3ed	forecast_graphs.m: fix wrong naming Also removes eval	2023-12-27 20:54:43 +01:00
Johannes Pfeifer	a6ad26788e	BVAR: remove dedicated graph functions one was empty, the other was moved to an inline function	2023-12-27 20:43:28 +01:00
Johannes Pfeifer	8f07f37138	🐛 makedataset.m: correct error message with first_obs specified	2023-12-27 17:51:57 +01:00
Johannes Pfeifer	e9d79796cf	dynare.m: cosmetic changes	2023-12-27 17:12:27 +01:00
Johannes Pfeifer	b7e9d69874	dynare.m: remove redundant legacy check for file ending	2023-12-27 17:12:15 +01:00