matlab/estimation/smc/ishssmc.m

+function bool = ishssmc(options_)
+
+% Copyright © 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/>.
+
+bool = false;
+if isfield(options_, 'posterior_sampler_options')
+    if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'hssmc')
+        bool = true;
+    end
+end

matlab/estimation/smc/logsumexp.m

+function [lse,sm] = logsumexp(x)
+%LOGSUMEXP  Log-sum-exp function.
+%    lse = LOGSUMEXP(x) returns the log-sum-exp function evaluated at 
+%    the vector x, defined by lse = log(sum(exp(x)).
+%    [lse,sm] = LOGSUMEXP(x) also returns the softmax function evaluated
+%    at x, defined by sm = exp(x)/sum(exp(x)).
+%    The functions are computed in a way that avoids overflow and 
+%    optimizes numerical stability.   
+%
+%    Reference:
+%    P. Blanchard, D. J. Higham, and N. J. Higham.  
+%    Accurately computing the log-sum-exp and softmax functions. 
+%    IMA J. Numer. Anal., Advance access, 2020.
+%
+%    Taken from https://de.mathworks.com/matlabcentral/fileexchange/84892-logsumexp-softmax
+%
+% Copyright (c) 2020, Nicholas J. Higham 
+% All rights reserved.
+%
+% 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.
+%
+% 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 HOLDER 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.
+
+if ~isvector(x), error('Input x must be a vector.'), end
+
+n = length(x);
+s = 0; e = zeros(n,1);
+[xmax,k] = max(x); a = xmax;
+s = 0;
+for i = 1:n
+    e(i) = exp(x(i)-xmax);
+    if i ~= k
+       s = s + e(i);
+    end   
+end
+lse = a + log1p(s);
+if nargout > 1
+   sm = e/(1+s);
+end  

matlab/estimation/smc/smc_samplers_initialization.m

+function [particles, tlogpostkernel, loglikelihood] = smc_samplers_initialization(objective_function, sampler, n, Prior, SimulationFolder, nsteps)
+% function [particles, tlogpostkernel, loglikelihood] = smc_samplers_initialization(objective_function, sampler, n, Prior, SimulationFolder, nsteps)
+% Initialize SMC samplers by drawing initial particles in the prior distribution.
+%
+% INPUTS
+% - objective_function [char]     string specifying the name of the objective function (posterior kernel).
+% - sampler          [char]     name of the sampler.
+% - n                [integer]  scalar, number of particles.
+% - Prior            [class]    prior information
+% - SimulationFolder [char]     name of output folder
+% - nsteps           [integer]  number of steps
+%
+% OUTPUTS
+% - particles             [double]    p×n matrix of particles
+% - tlogpostkernel        [double]    n×1 vector of posterior kernel values for the particles
+% - loglikelihood         [double]    n×1 vector of likelihood values for the particles
+%
+% SPECIAL REQUIREMENTS
+%   None.
+
+% Copyright © 2022-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/>.
+
+dprintf('Estimation:%s: Initialization...', sampler)
+
+% Delete old mat files storign particles if any...
+matfiles = sprintf('%s%sparticles*.mat', SimulationFolder, filesep());
+files = dir(matfiles);
+if ~isempty(files)
+    delete(matfiles);
+    dprintf('Estimation:%s: Old %s-files successfully erased.', sampler, sampler)
+end
+
+% Simulate a pool of particles characterizing the prior distribution (with the additional constraint that the likelihood is finite)
+set_dynare_seed('default');
+dprintf('Estimation:%s: Searching for initial values...', sampler);
+particles = zeros(Prior.length(), n);
+tlogpostkernel = zeros(n, 1);
+loglikelihood = zeros(n, 1);
+
+t0 = tic;
+parfor j=1:n
+    notvalid = true;
+    while notvalid
+        candidate = Prior.draw();
+        if Prior.admissible(candidate)
+            particles(:,j) = candidate;
+            [tlogpostkernel(j), loglikelihood(j)] = tempered_likelihood(objective_function, candidate, 0.0, Prior);
+            if isfinite(loglikelihood(j)) % if returned log-density is Inf or Nan (penalized value)
+                notvalid = false;
+            end
+        end
+    end
+end
+tt = toc(t0);
+
+save(sprintf('%s%sparticles-1-%u.mat', SimulationFolder, filesep(), nsteps), 'particles', 'tlogpostkernel', 'loglikelihood')
+dprintf('Estimation:%s: Initial values found (%.2fs)', sampler, tt)
+skipline()

matlab/estimation/smc/tempered_likelihood.m

+function [tlogpostkernel,loglikelihood] = tempered_likelihood(objective_function, xparam, lambda, Prior)
+
+% Evaluate tempered likelihood (posterior kernel)
+%
+% INPUTS
+% - objective_function  [handle]   Function handle for the opposite of the  posterior kernel.
+% - xparam              [double]   n×1 vector of parameters.
+% - lambda              [double]   scalar between 0 and 1, weight on the posterior kernel.
+% - Prior               [dprior]   Prior specification.
+%
+% OUTPUTS
+% - tlogpostkernel      [double]   scalar, value of the tempered posterior kernel.
+% - loglikelihood       [double]   scalar, value of the log likelihood.
+
+% 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/>.
+
+logpostkernel = -objective_function(xparam);
+logprior = Prior.density(xparam);
+loglikelihood = logpostkernel-logprior;
+tlogpostkernel = lambda*loglikelihood + logprior;

matlab/estimation/smc/trace_plot_dime.m

+function lprobs_sample = trace_plot_dime(options_, M_)
+
+% Plot the history of the densities of an ensemble to visually inspect convergence.
+%
+% INPUTS
+% - options_         [struct]   Dynare's options
+% - M_               [struct]   model description
+% - oo_              [struct]   outputs
+%
+% SPECIAL REQUIREMENTS
+% None.
+
+% Copyright © 2022-2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+lprobs = GetAllPosteriorDraws(options_, M_.dname, [], 0);
+tune = options_.posterior_sampler_options.current_options.tune;
+[niter, nchain] = size(lprobs);
+
+if max(lprobs(end,:)) > 0
+    ulim = max(lprobs(end,:))*1.2;
+else
+    ulim = max(lprobs(end,:))/1.2;
+end
+if min(lprobs(end,:)) > 0
+    llim = min(lprobs(end,:))/5;
+else
+    llim = min(lprobs(end,:))*5;
+end
+
+graphFolder = CheckPath('graphs',M_.dname);
+hh_fig = dyn_figure(options_.nodisplay,'Name','DIME Convergence Diagnostics');
+
+hold on
+lines_tune = plot(1:niter-tune, lprobs(1:end-tune,:), 'Color', '#D95319');
+lines_sample = plot(niter-tune:niter, lprobs(end-tune:end,:), 'Color', '#0072BD');
+if ~isoctave
+    % Set the transparency (alpha channel) of line objects (undocumented MATLAB feature)
+    % See https://fr.mathworks.com/matlabcentral/discussions/ideas/833472-add-alpha-capability-to-line-class
+    for i = 1:nchain
+        lines_tune(i).Color(4) = min(1,10/nchain);
+        lines_sample(i).Color(4) = min(1,10/nchain);
+    end
+end
+ylim([llim, ulim])
+hold off
+
+dyn_saveas(hh_fig,[graphFolder '/' M_.fname '_trace_lprob'],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=0.8\\textwidth]{%s_trace_lprob}\n',[graphFolder '/' M_.fname]);
+    fprintf(fidTeX,'\\caption{Ensemble traces of posterior densities for DIME.\n');
+    fprintf(fidTeX,'\\label{Fig:DIME_trace}\n');
+    fprintf(fidTeX,'\\end{figure}\n');
+    fprintf(fidTeX,'\n');
+end
+
+lprobs_sample = lprobs(end-tune:end,:);

matlab/trace_plot.m → matlab/estimation/trace_plot.m

@@ -52,22 +52,22 @@ if isempty(column)
     return
 end
 
-% Get informations about the posterior draws:
-MetropolisFolder = CheckPath('metropolis',M_.dname);
-record=load_last_mh_history_file(MetropolisFolder, M_.fname);
-
-FirstMhFile = 1;
-FirstLine = 1;
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-[mh_nblck] = size(record.LastParameters,2);
-clear record;
-
+if ~issmc(options_)
+    options_.mh_drop=0; % locally, do not drop as we want all draws
     n_nblocks_to_plot=length(blck);
+else
+    if ishssmc(options_)
+        n_nblocks_to_plot=1;
+    elseif isdime(options_)
+        error('trace_plot:: DIME does not support the trace_plot command')
+    end
+end
+
+[~, ~, TotalNumberOfMhDraws]=set_number_of_subdraws(M_,options_); 
 
 if n_nblocks_to_plot==1
 % Get all the posterior draws:
-PosteriorDraws = GetAllPosteriorDraws(M_.dname,M_.fname,column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
+    PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, 1, 1, blck);
 else
     PosteriorDraws=NaN(TotalNumberOfMhDraws,n_nblocks_to_plot);
     save_string='';
@@ -75,7 +75,7 @@ else
         title_string_tex='';
     end
     for block_iter=1:n_nblocks_to_plot
-        PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
+        PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, 1, 1, blck(block_iter));
         save_string=[save_string,'_',num2str(blck(block_iter))];
         if options_.TeX
             title_string_tex=[title_string_tex, ', ' num2str(blck(block_iter))];
@@ -153,7 +153,7 @@ end
 if strcmpi(type,'PosteriorDensity')
     plot_name='Posterior';
 else
-    plot_name=get_the_name(column,0,M_,estim_params_,options_);
+    plot_name=get_the_name(column,0,M_,estim_params_,options_.varobs);
 end
 if n_nblocks_to_plot==1
     plot_name=[plot_name,'_blck_',num2str(blck)];

matlab/update_last_mh_history_file.m → matlab/estimation/update_last_mh_history_file.m

@@ -29,7 +29,7 @@ BaseName = [MetropolisFolder filesep ModelName];
 % Get the list of all the mh_history files.
 mh_history_files = dir([BaseName '_mh_history_*.mat']);
 
-% Check the existence of mh-files (assuming version 2, ie dynare version greater than 4.3.x).
+% Check the existence of mh-files (assuming version 2, i.e., Dynare version greater than 4.3.x).
 if isequal(length(mh_history_files),0)
     error(['update_last_mh_history_file: I cannot find any mh-history file in ' MetropolisFolder '!'])
 end

matlab/var_sample_moments.m → matlab/estimation/var_sample_moments.m

-function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)%datafile,varobs,xls_sheet,xls_range)
+function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)
+% dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, dataset_info)
 % Computes the sample moments of a VAR model.
 %
 % The VAR(p) model is defined by:
@@ -51,7 +52,7 @@ function dataset_info=var_sample_moments(nlag, var_trend_order, dataset_, datase
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright © 2007-2021 Dynare Team
+% Copyright © 2007-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -80,14 +81,11 @@ if isequal(var_trend_order,-1)
 elseif isequal(var_trend_order,0)
     % Constant and no linear trend case.
     X = ones(NumberOfObservations,NumberOfVariables*nlag+1);
-    indx = NumberOfVariables*nlag+1;
 elseif isequal(var_trend_order,1)
     % Constant and linear trend case.
     X = ones(NumberOfObservations,NumberOfVariables*nlag+2);
-    indx = NumberOfVariables*nlag+1:NumberOfVariables*nlag+2;
 else
     error('Estimation::var_sample_moments: trend must be equal to -1,0 or 1!')
-    return
 end
 
 % I build matrices Y and X
@@ -109,9 +107,3 @@ dataset_info.mXY=X'*Y;
 dataset_info.mXX=X'*X;
 dataset_info.Ydata=Y;
 dataset_info.Xdata=X;
\ No newline at end of file
-% assignin('base', 'mYY', Y'*Y);
-% assignin('base', 'mYX', Y'*X);
-% assignin('base', 'mXY', X'*Y);
-% assignin('base', 'mXX', X'*X);
-% assignin('base', 'Ydata', Y);
-% assignin('base', 'Xdata', X);
\ No newline at end of file

matlab/variance_decomposition_ME_mc_analysis.m → matlab/estimation/variance_decomposition_ME_mc_analysis.m

@@ -25,7 +25,7 @@ function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dn
 
 
 
-% Copyright © 2017 Dynare Team
+% Copyright © 2017-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -43,8 +43,9 @@ function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dn
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
@@ -87,7 +88,7 @@ ListOfFiles = dir([ PATH  fname '_' TYPE 'VarianceDecompME*.mat']);
 i1 = 1; tmp = zeros(NumberOfSimulations,1);
 indice = (indx-1)*rows(exonames)+jndx;
 for file = 1:length(ListOfFiles)
-    load([ PATH ListOfFiles(file).name ]);
+    load([ PATH ListOfFiles(file).name ],'Decomposition_array_ME');
     i2 = i1 + rows(Decomposition_array_ME) - 1;
     tmp(i1:i2) = Decomposition_array_ME(:,indice);
     i1 = i2+1;
@@ -95,10 +96,10 @@ end
 
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 
 oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Mean.(var).(exo) = p_mean;

matlab/variance_decomposition_mc_analysis.m → matlab/estimation/variance_decomposition_mc_analysis.m

@@ -25,7 +25,7 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
 
 
 
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -43,8 +43,9 @@ function oo_ = variance_decomposition_mc_analysis(NumberOfSimulations,type,dname
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if strcmpi(type,'posterior')
+    folder_name=get_posterior_folder_name(options_);
     TYPE = 'Posterior';
-    PATH = [dname '/metropolis/'];
+    PATH = [dname filesep folder_name filesep];
 else
     TYPE = 'Prior';
     PATH = [dname '/prior/moments/'];
@@ -85,7 +86,7 @@ ListOfFiles = dir([ PATH  fname '_' TYPE 'VarianceDecomposition*.mat']);
 i1 = 1; tmp = zeros(NumberOfSimulations,1);
 indice = (indx-1)*rows(exonames)+jndx;
 for file = 1:length(ListOfFiles)
-    load([ PATH ListOfFiles(file).name ]);
+    load([ PATH ListOfFiles(file).name ],'Decomposition_array');
     i2 = i1 + rows(Decomposition_array) - 1;
     tmp(i1:i2) = Decomposition_array(:,indice);
     i1 = i2+1;
@@ -93,10 +94,10 @@ end
 
 if options_.estimation.moments_posterior_density.indicator
     [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
-        posterior_moments(tmp,1,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 else
     [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
-        posterior_moments(tmp,0,mh_conf_sig);
+        posterior_moments(tmp,mh_conf_sig);
 end
 
 oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecomposition.Mean.(var).(exo) = p_mean;

matlab/evaluate_static_model.m

-function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,options)
-
-% function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,options)
+function [residuals,jacob] = evaluate_static_model(ys,exo_ss,params,M_,options_)
+% function [residuals,jacob] = evaluate_static_model(ys,exo_ss,params,M_,options_)
 % Evaluates the static model
 %
 % INPUTS
@@ -8,19 +7,18 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
 %                                                 state
 %   exo_ss                    vector           exogenous steady state
 %   params                    vector           parameters
-%   M                         struct           model structure
-%   options                   struct           options
+%   M_                        struct           model structure
+%   options_                  struct           options
 %
 % OUTPUTS
 %   residuals                 vector           residuals when ys is not
 %                                              the steady state
-%   check1                    scalar           error flag
 %   jacob                     matrix           Jacobian of static model
 %
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,19 +35,18 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-check1 = 0;
-if options.bytecode
+if options_.bytecode
     if nargout<3
-        [residuals]= bytecode('evaluate','static',ys,...
+        [residuals]= bytecode('evaluate', 'static', M_, options_, ys, ...
                          exo_ss, params, ys, 1);
     else
-        [residuals, junk]= bytecode('evaluate','static',ys,...
+        [residuals, junk]= bytecode('evaluate', 'static', M_, options_, ys, ...
             exo_ss, params, ys, 1);
         jacob = junk.g1;
     end      
 else
-    [residuals, T_order, T] = feval([M.fname '.sparse.static_resid'], ys, exo_ss, params);
+    [residuals, T_order, T] = feval([M_.fname '.sparse.static_resid'], ys, exo_ss, params);
     if nargout >= 3
-        jacob = feval([M.fname '.sparse.static_g1'], ys, exo_ss, params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
+        jacob = feval([M_.fname '.sparse.static_g1'], ys, exo_ss, params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
     end
 end

matlab/evaluate_steady_state.m

-function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,steadystate_check_flag)
-% function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,steadystate_check_flag)
+function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M_,options_,steadystate_check_flag)
+% function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M_,options_,steadystate_check_flag)
 % Computes the steady state
 %
 % INPUTS
 %   ys_init                   vector           initial values used to compute the steady
 %                                                 state
 %   exo_ss                    vector           exogenous steady state (incl. deterministic exogenous)
-%   M                         struct           model structure
-%   options                   struct           options
+%   M_                        struct           model structure
+%   options_                  struct           options
 %   steadystate_check_flag    boolean          if true, check that the
 %                                              steadystate verifies the
 %                                              static model
@@ -22,7 +22,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,stead
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,21 +39,16 @@ function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M,options,stead
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if options.solve_algo < 0 || options.solve_algo > 14
+if options_.solve_algo < 0 || options_.solve_algo > 14
     error('STEADY: solve_algo must be between 0 and 14')
 end
 
-if ~options.bytecode && ~options.block && options.solve_algo > 4 && ...
-        options.solve_algo < 9
-    error('STEADY: you can''t use solve_algo = {5,6,7,8} without block nor bytecode options')
-end
-
-if ~options.bytecode && options.block && options.solve_algo == 5
+if ~options_.bytecode && options_.block && options_.solve_algo == 5
     error('STEADY: you can''t use solve_algo = 5 without bytecode option')
 end
 
-if isoctave && options.solve_algo == 11
-    error(['STEADY: you can''t use solve_algo = %u under Octave'],options.solve_algo)
+if isoctave && options_.solve_algo == 11
+    error('STEADY: you can''t use solve_algo = %u under Octave',options_.solve_algo)
 end
 
 % To ensure that the z and zx matrices constructed by repmat and passed to bytecode
@@ -68,36 +63,36 @@ end
 info = 0;
 check = 0;
 
-steadystate_flag = options.steadystate_flag;
-params = M.params;
+steadystate_flag = options_.steadystate_flag;
+params = M_.params;
 
-if length(M.aux_vars) > 0 && ~steadystate_flag && M.set_auxiliary_variables
-    h_set_auxiliary_variables = str2func([M.fname '.set_auxiliary_variables']);
+if length(M_.aux_vars) > 0 && ~steadystate_flag && M_.set_auxiliary_variables
+    h_set_auxiliary_variables = str2func([M_.fname '.set_auxiliary_variables']);
     ys_init = h_set_auxiliary_variables(ys_init,exo_ss,params);
 end
 
-if options.ramsey_policy
-    if ~isfinite(M.params(strmatch('optimal_policy_discount_factor',M.param_names,'exact')))
+if options_.ramsey_policy
+    if ~isfinite(M_.params(strmatch('optimal_policy_discount_factor',M_.param_names,'exact')))
         fprintf('\nevaluate_steady_state: the planner_discount is NaN/Inf. That will cause problems.\n')
     end
     if steadystate_flag
         % explicit steady state file
-        [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, ...
-                                                      options,steadystate_check_flag);
+        [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_, ...
+                                                      options_,steadystate_check_flag);
         %test whether it solves model conditional on the instruments
-        if ~options.debug
-            resids = evaluate_static_model(ys,exo_ss,params,M,options);
+        if ~options_.debug
+            resids = evaluate_static_model(ys,exo_ss,params,M_,options_);
         else
-            [resids, ~ , jacob]= evaluate_static_model(ys,exo_ss,params,M,options);
+            [resids, jacob]= evaluate_static_model(ys,exo_ss,params,M_,options_);
         end
-        nan_indices=find(isnan(resids(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr))));
+        nan_indices=find(isnan(resids(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr))));
 
         if ~isempty(nan_indices)
-            if options.debug
+            if options_.debug
                 fprintf('\nevaluate_steady_state: The steady state file computation for the Ramsey problem resulted in NaNs.\n')
                 fprintf('evaluate_steady_state: The steady state was computed conditional on the following initial instrument values: \n')
-                for ii = 1:size(options.instruments,1)
-                    fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
+                for ii = 1:size(options_.instruments,1)
+                    fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
                 end
                 fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
                 fprintf('\t Equation(s): ')
@@ -113,12 +108,12 @@ if options.ramsey_policy
             return
         end
 
-        if any(imag(ys(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr))))
-            if options.debug
+        if any(imag(ys(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr))))
+            if options_.debug
                 fprintf('\nevaluate_steady_state: The steady state file computation for the Ramsey problem resulted in complex numbers.\n')
                 fprintf('evaluate_steady_state: The steady state was computed conditional on the following initial instrument values: \n')
-                for ii = 1:size(options.instruments,1)
-                    fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
+                for ii = 1:size(options_.instruments,1)
+                    fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
                 end
                 fprintf('evaluate_steady_state: If those initial values are not admissable, change them using an initval-block.\n')
                 skipline(2)
@@ -128,17 +123,17 @@ if options.ramsey_policy
             return
         end
 
-        if max(abs(resids(M.ramsey_orig_endo_nbr+(1:M.ramsey_orig_eq_nbr)))) > options.solve_tolf %does it solve for all variables except for the Lagrange multipliers
-            if options.debug
+        if max(abs(resids(M_.ramsey_orig_endo_nbr+(1:M_.ramsey_orig_eq_nbr)))) > options_.solve_tolf %does it solve for all variables except for the Lagrange multipliers
+            if options_.debug
                 fprintf('\nevaluate_steady_state: The steady state file does not solve the steady state for the Ramsey problem.\n')
                 fprintf('evaluate_steady_state: Conditional on the following instrument values: \n')
-                for ii = 1:size(options.instruments,1)
-                    fprintf('\t %s \t %f \n',options.instruments{ii},ys_init(strmatch(options.instruments{ii},M.endo_names,'exact')))
+                for ii = 1:size(options_.instruments,1)
+                    fprintf('\t %s \t %f \n',options_.instruments{ii},ys_init(strmatch(options_.instruments{ii},M_.endo_names,'exact')))
                 end
                 fprintf('evaluate_steady_state: the following equations have non-zero residuals: \n')
-                for ii=M.ramsey_orig_endo_nbr+1:M.endo_nbr
-                    if abs(resids(ii)) > options.solve_tolf
-                        fprintf('\t Equation number %d: %f\n',ii-M.ramsey_orig_endo_nbr, resids(ii))
+                for ii=M_.ramsey_orig_endo_nbr+1:M_.endo_nbr
+                    if abs(resids(ii)) > options_.solve_tolf
+                        fprintf('\t Equation number %d: %f\n',ii-M_.ramsey_orig_endo_nbr, resids(ii))
                     end
                 end
                 skipline(2)
@@ -148,31 +143,31 @@ if options.ramsey_policy
             return
         end
     end
-    if options.debug
+    if options_.debug
         if steadystate_flag
-            infrow=find(isinf(ys_init(1:M.orig_endo_nbr)));
+            infrow=find(isinf(ys_init(1:M_.orig_endo_nbr)));
         else
             infrow=find(isinf(ys_init));
         end
         if ~isempty(infrow)
             fprintf('\nevaluate_steady_state: The initial values for the steady state of the following variables are Inf:\n');
             for iter=1:length(infrow)
-                fprintf('%s\n',M.endo_names{infrow(iter)});
+                fprintf('%s\n',M_.endo_names{infrow(iter)});
             end
         end
         if steadystate_flag
-            nanrow=find(isnan(ys_init(1:M.orig_endo_nbr)));
+            nanrow=find(isnan(ys_init(1:M_.orig_endo_nbr)));
         else
             nanrow=find(isnan(ys_init));
         end
         if ~isempty(nanrow)
             fprintf('\nevaluate_steady_state: The initial values for the steady state of the following variables are NaN:\n');
             for iter=1:length(nanrow)
-                fprintf('%s\n',M.endo_names{nanrow(iter)});
+                fprintf('%s\n',M_.endo_names{nanrow(iter)});
             end
         end
         if steadystate_flag
-            nan_indices_mult=find(isnan(resids(1:M.ramsey_orig_endo_nbr)));
+            nan_indices_mult=find(isnan(resids(1:M_.ramsey_orig_endo_nbr)));
             if any(nan_indices_mult)
                 fprintf('evaluate_steady_state: The steady state results NaN for auxiliary equation %u.\n',nan_indices_mult);
                 fprintf('evaluate_steady_state: This is often a sign of problems.\n');
@@ -181,41 +176,41 @@ if options.ramsey_policy
             
             if ~isempty(infrow)
                 fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains Inf. The problem is associated with:\n\n')
-                display_problematic_vars_Jacobian(infrow,infcol,M,ys,'static','evaluate_steady_state: ')
+                display_problematic_vars_Jacobian(infrow,infcol,M_,ys,'static','evaluate_steady_state: ')
             end
             
             if ~isreal(jacob)
                 [imagrow,imagcol]=find(abs(imag(jacob))>1e-15);
                 fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
-                display_problematic_vars_Jacobian(imagrow,imagcol,M,ys,'static','evaluate_steady_state: ')
+                display_problematic_vars_Jacobian(imagrow,imagcol,M_,ys,'static','evaluate_steady_state: ')
             end
             
             [nanrow,nancol]=find(isnan(jacob));
             if ~isempty(nanrow)
                 fprintf('\nevaluate_steady_state: The Jacobian of the dynamic model contains NaN. The problem is associated with:\n\n')
-                display_problematic_vars_Jacobian(nanrow,nancol,M,ys,'static','evaluate_steady_state: ')
+                display_problematic_vars_Jacobian(nanrow,nancol,M_,ys,'static','evaluate_steady_state: ')
             end
             
         end
 
     end
     %either if no steady state file or steady state file without problems
-    [ys,params,info] = dyn_ramsey_static(ys_init,exo_ss,M,options);
+    [ys,params,info] = dyn_ramsey_static(ys_init,exo_ss,M_,options_);
     if info
         return
     end
     %check whether steady state really solves the model
-    resids = evaluate_static_model(ys,exo_ss,params,M,options);
+    resids = evaluate_static_model(ys,exo_ss,params,M_,options_);
 
-    nan_indices_multiplier=find(isnan(resids(1:M.ramsey_orig_endo_nbr)));
-    nan_indices=find(isnan(resids(M.ramsey_orig_endo_nbr+1:end)));
+    nan_indices_multiplier=find(isnan(resids(1:M_.ramsey_orig_endo_nbr)));
+    nan_indices=find(isnan(resids(M_.ramsey_orig_endo_nbr+1:end)));
 
     if ~isempty(nan_indices)
-        if options.debug
+        if options_.debug
             fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs.\n')
             fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n')
-            for i = 1:size(options.instruments,1)
-                fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
+            for i = 1:size(options_.instruments,1)
+                fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
             end
             fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
             fprintf('\t Equation(s): ')
@@ -229,11 +224,11 @@ if options.ramsey_policy
     end
 
     if ~isempty(nan_indices_multiplier)
-        if options.debug
+        if options_.debug
             fprintf('\nevaluate_steady_state: The steady state computation for the Ramsey problem resulted in NaNs in the auxiliary equations.\n')
             fprintf('evaluate_steady_state: The steady state computation resulted in the following instrument values: \n')
-            for i = 1:size(options.instruments,1)
-                fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
+            for i = 1:size(options_.instruments,1)
+                fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
             end
             fprintf('evaluate_steady_state: The problem occured in the following equations: \n')
             fprintf('\t Auxiliary equation(s): ')
@@ -246,22 +241,22 @@ if options.ramsey_policy
         return
     end
 
-    if max(abs(resids)) > options.solve_tolf %does it solve for all variables including the auxiliary ones
-        if options.debug
+    if max(abs(resids)) > options_.solve_tolf %does it solve for all variables including the auxiliary ones
+        if options_.debug
             fprintf('\nevaluate_steady_state: The steady state for the Ramsey problem could not be computed.\n')
             fprintf('evaluate_steady_state: The steady state computation stopped with the following instrument values:: \n')
-            for i = 1:size(options.instruments,1)
-                fprintf('\t %s \t %f \n',options.instruments{i},ys(strmatch(options.instruments{i},M.endo_names,'exact')))
+            for i = 1:size(options_.instruments,1)
+                fprintf('\t %s \t %f \n',options_.instruments{i},ys(strmatch(options_.instruments{i},M_.endo_names,'exact')))
             end
             fprintf('evaluate_steady_state: The following equations have non-zero residuals: \n')
-            for ii=1:M.ramsey_orig_endo_nbr
-                if abs(resids(ii)) > options.solve_tolf/100
+            for ii=1:M_.ramsey_orig_endo_nbr
+                if abs(resids(ii)) > options_.solve_tolf/100
                     fprintf('\t Auxiliary Ramsey equation number %d: %f\n',ii, resids(ii))
                 end
             end
-            for ii=M.ramsey_orig_endo_nbr+1:M.endo_nbr
-                if abs(resids(ii)) > options.solve_tolf/100
-                    fprintf('\t Equation number %d: %f\n',ii-M.ramsey_orig_endo_nbr, resids(ii))
+            for ii=M_.ramsey_orig_endo_nbr+1:M_.endo_nbr
+                if abs(resids(ii)) > options_.solve_tolf/100
+                    fprintf('\t Equation number %d: %f\n',ii-M_.ramsey_orig_endo_nbr, resids(ii))
                 end
             end
             skipline(2)
@@ -272,46 +267,46 @@ if options.ramsey_policy
     end
 elseif steadystate_flag
     % explicit steady state file
-    [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, options,steadystate_check_flag);
+    [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_, options_,steadystate_check_flag);
     if size(ys,2)>size(ys,1)
         error('STEADY: steady_state-file must return a column vector, not a row vector.')
     end
     if info(1)
         return
     end
-elseif ~options.bytecode && ~options.block
-    static_resid = str2func(sprintf('%s.sparse.static_resid', M.fname));
-    static_g1 = str2func(sprintf('%s.sparse.static_g1', M.fname));
-    if ~options.linear
+elseif ~options_.bytecode && ~options_.block
+    static_resid = str2func(sprintf('%s.sparse.static_resid', M_.fname));
+    static_g1 = str2func(sprintf('%s.sparse.static_g1', M_.fname));
+    if ~options_.linear
         % non linear model
-        if  ismember(options.solve_algo,[10,11])
-            [lb,ub,eq_index] = get_complementarity_conditions(M,options.ramsey_policy);
-            if options.solve_algo == 10
-                options.lmmcp.lb = lb;
-                options.lmmcp.ub = ub;
-            elseif options.solve_algo == 11
-                options.mcppath.lb = lb;
-                options.mcppath.ub = ub;
+        if  ismember(options_.solve_algo,[10,11])
+            [lb, ub] = feval(sprintf('%s.static_complementarity_conditions', M_.fname), M_.params);
+            if options_.solve_algo == 10
+                options_.lmmcp.lb = lb;
+                options_.lmmcp.ub = ub;
+            elseif options_.solve_algo == 11
+                options_.mcppath.lb = lb;
+                options_.mcppath.ub = ub;
             end
             [ys,check,fvec, fjac, errorcode] = dynare_solve(@static_mcp_problem,...
                 ys_init,...
-                options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
-                options, exo_ss, params,...
-                M.endo_nbr, static_resid, static_g1, ...
-                M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, eq_index);
+                options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
+                options_, exo_ss, params,...
+                M_.endo_nbr, static_resid, static_g1, ...
+                M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, M_.static_mcp_equations_reordering);
         else
             [ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
-                options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
-                options, exo_ss, params, M.endo_nbr, static_resid, static_g1, ...
-                M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
+                options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
+                options_, exo_ss, params, M_.endo_nbr, static_resid, static_g1, ...
+                M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr);
         end
-        if check && options.debug
+        if check && options_.debug
             dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
             skipline()
             [infrow,infcol]=find(isinf(fjac) | isnan(fjac));
             if ~isempty(infrow)
                 fprintf('\nSTEADY:  The Jacobian at the initial values contains Inf or NaN. The problem arises from: \n')
-                display_problematic_vars_Jacobian(infrow,infcol,M,ys_init,'static','STEADY: ')
+                display_problematic_vars_Jacobian(infrow,infcol,M_,ys_init,'static','STEADY: ')
             end
             problematic_equation = find(~isfinite(fvec));
             if ~isempty(problematic_equation)
@@ -325,7 +320,7 @@ elseif ~options.bytecode && ~options.block
     else
         % linear model
         [fvec, T_order, T] = static_resid(ys_init, exo_ss, params);
-        jacob = static_g1(ys_init, exo_ss, params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T);
+        jacob = static_g1(ys_init, exo_ss, params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
 
         ii = find(~isfinite(fvec));
         if ~isempty(ii)
@@ -337,7 +332,7 @@ elseif ~options.bytecode && ~options.block
             disp('Check whether your model is truly linear. Put "resid(1);" before "steady;" to see the problematic equations.')
         elseif isempty(ii) && max(abs(fvec)) > 1e-12
             ys = ys_init-jacob\fvec;
-            resid = evaluate_static_model(ys,exo_ss,params,M,options);
+            resid = evaluate_static_model(ys,exo_ss,params,M_,options_);
             if max(abs(resid)) > 1e-6
                 check=1;
                 fprintf('STEADY: No steady state for your model could be found\n')
@@ -346,39 +341,39 @@ elseif ~options.bytecode && ~options.block
         else
             ys = ys_init;
         end
-        if options.debug
+        if options_.debug
             if any(any(isinf(jacob) | isnan(jacob)))
                 [infrow,infcol]=find(isinf(jacob) | isnan(jacob));
                 fprintf('\nSTEADY:  The Jacobian contains Inf or NaN. The problem arises from: \n\n')
                 for ii=1:length(infrow)
-                    fprintf('STEADY:  Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',infrow(ii),M.endo_names{infcol(ii),:},M.endo_names{infcol(ii),:},ys_init(infcol(ii)))
+                    fprintf('STEADY:  Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',infrow(ii),M_.endo_names{infcol(ii),:},M_.endo_names{infcol(ii),:},ys_init(infcol(ii)))
                 end
                 fprintf('Check whether your model is truly linear. Put "resid(1);" before "steady;" to see the problematic equations.\n')
             end
         end
     end
-elseif ~options.bytecode && options.block
+elseif ~options_.bytecode && options_.block
     ys = ys_init;
-    T = NaN(M.block_structure_stat.tmp_nbr, 1);
-    for b = 1:length(M.block_structure_stat.block)
-        fh_static = str2func(sprintf('%s.sparse.block.static_%d', M.fname, b));
-        if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
-                M.block_structure_stat.block(b).Simulation_Type ~= 2
-            mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
-            if options.solve_algo <= 4 || options.solve_algo >= 9
+    T = NaN(M_.block_structure_stat.tmp_nbr, 1);
+    for b = 1:length(M_.block_structure_stat.block)
+        fh_static = str2func(sprintf('%s.sparse.block.static_%d', M_.fname, b));
+        if M_.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
+                M_.block_structure_stat.block(b).Simulation_Type ~= 2
+            mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
+            if options_.solve_algo <= 4 || options_.solve_algo >= 9
                 [ys(mfs_idx), errorflag] = dynare_solve(@block_mfs_steadystate, ys(mfs_idx), ...
-                                                        options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
-                                                        options, fh_static, b, ys, exo_ss, params, T, M);
+                                                        options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
+                                                        options_, fh_static, b, ys, exo_ss, params, T, M_);
                 if errorflag
                     check = 1;
                     break
                 end
             else
-                nze = length(M.block_structure_stat.block(b).g1_sparse_rowval);
+                nze = length(M_.block_structure_stat.block(b).g1_sparse_rowval);
                 [ys, T, success] = solve_one_boundary(fh_static, ys, exo_ss, ...
-                                                      params, [], T, mfs_idx, nze, 1, false, b, 0, options.steady.maxit, ...
-                                                      options.solve_tolf, ...
-                                                      0, options.solve_algo, true, false, false, M, options);
+                                                      params, [], T, mfs_idx, nze, 1, false, b, 0, options_.steady.maxit, ...
+                                                      options_.solve_tolf, ...
+                                                      0, options_.solve_algo, true, false, false, M_, options_);
                 if ~success
                     check = 1;
                     break
@@ -387,36 +382,36 @@ elseif ~options.bytecode && options.block
         end
         % Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block.
         % Also update the temporary terms vector (needed for the dynare_solve case)
-        [ys, T] = fh_static(ys, exo_ss, params, M.block_structure_stat.block(b).g1_sparse_rowval, ...
-                            M.block_structure_stat.block(b).g1_sparse_colval, ...
-                            M.block_structure_stat.block(b).g1_sparse_colptr, T);
+        [ys, T] = fh_static(ys, exo_ss, params, M_.block_structure_stat.block(b).g1_sparse_rowval, ...
+                            M_.block_structure_stat.block(b).g1_sparse_colval, ...
+                            M_.block_structure_stat.block(b).g1_sparse_colptr, T);
     end
-elseif options.bytecode
-    if options.solve_algo >= 5 && options.solve_algo <= 8
+elseif options_.bytecode
+    if options_.solve_algo >= 5 && options_.solve_algo <= 8
         try
-            if options.block
-                ys = bytecode('static', 'block_decomposed', ys_init, exo_ss, params);
+            if options_.block
+                ys = bytecode('static', 'block_decomposed', M_, options_, ys_init, exo_ss, params);
             else
-                ys = bytecode('static', ys_init, exo_ss, params);
+                ys = bytecode('static', M_, options_, ys_init, exo_ss, params);
             end
         catch ME
-            if options.verbosity >= 1
+            if options_.verbosity >= 1
                 disp(ME.message);
             end
             ys = ys_init;
             check = 1;
         end
-    elseif options.block
+    elseif options_.block
         ys = ys_init;
-        T = NaN(M.block_structure_stat.tmp_nbr, 1);
-        for b = 1:length(M.block_structure_stat.block)
-            if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
-                    M.block_structure_stat.block(b).Simulation_Type ~= 2
-                mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
+        T = NaN(M_.block_structure_stat.tmp_nbr, 1);
+        for b = 1:length(M_.block_structure_stat.block)
+            if M_.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
+                    M_.block_structure_stat.block(b).Simulation_Type ~= 2
+                mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
                 [ys(mfs_idx), errorflag] = dynare_solve(@block_bytecode_mfs_steadystate, ...
-                                                        ys(mfs_idx), options.steady.maxit, ...
-                                                        options.solve_tolf, options.solve_tolx, ...
-                                                        options, b, ys, exo_ss, params, T, M);
+                                                        ys(mfs_idx), options_.steady.maxit, ...
+                                                        options_.solve_tolf, options_.solve_tolx, ...
+                                                        options_, b, ys, exo_ss, params, T, M_, options_);
                 if errorflag
                     check = 1;
                     break
@@ -425,10 +420,10 @@ elseif options.bytecode
             % Compute endogenous if the block is of type evaluate forward/backward or if there are recursive variables in a solve block.
             % Also update the temporary terms vector (needed for the dynare_solve case)
             try
-                [~, ~, ys, T] = bytecode(ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ...
+                [~, ~, ys, T] = bytecode(M_, options_, ys, exo_ss, params, ys, 1, ys, T, 'evaluate', 'static', ...
                                          'block_decomposed', ['block=' int2str(b)]);
             catch ME
-                if options.verbosity >= 1
+                if options_.verbosity >= 1
                     disp(ME.message);
                 end
                 check = 1;
@@ -437,20 +432,20 @@ elseif options.bytecode
         end
     else
         [ys, check] = dynare_solve(@bytecode_steadystate, ys_init, ...
-                                   options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
-                                   options, exo_ss, params);
+                                   options_.steady.maxit, options_.solve_tolf, options_.solve_tolx, ...
+                                   options_, exo_ss, params, M_, options_);
     end
 end
 
 if check
     info(1)= 20;
     %make sure ys contains auxiliary variables in case of problem with dynare_solve
-    if length(M.aux_vars) > 0 && ~steadystate_flag
-        if M.set_auxiliary_variables
+    if length(M_.aux_vars) > 0 && ~steadystate_flag
+        if M_.set_auxiliary_variables
             ys = h_set_auxiliary_variables(ys,exo_ss,params);
         end
     end
-    resid = evaluate_static_model(ys,exo_ss,params,M,options);
+    resid = evaluate_static_model(ys,exo_ss,params,M_,options_);
     info(2) = resid'*resid ;
     if isnan(info(2))
         info(1)=22;
@@ -459,18 +454,18 @@ if check
 end
 
 % If some equations are tagged [static] or [dynamic], verify consistency
-if M.static_and_dynamic_models_differ
+if M_.static_and_dynamic_models_differ
     % Evaluate residual of *dynamic* model using the steady state
     % computed on the *static* one
-    if options.bytecode
-        z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
-        zx = repmat(exo_ss', M.maximum_lead + M.maximum_lag + 1, 1);
-        r = bytecode('dynamic','evaluate', z, zx, params, ys, 1);
+    if options_.bytecode
+        z = repmat(ys,1,M_.maximum_lead + M_.maximum_lag + 1);
+        zx = repmat(exo_ss', M_.maximum_lead + M_.maximum_lag + 1, 1);
+        r = bytecode('dynamic','evaluate', M_, options_, z, zx, params, ys, 1);
     else
-        r = feval([M.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
+        r = feval([M_.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
     end
     % Fail if residual greater than tolerance
-    if max(abs(r)) > options.solve_tolf
+    if max(abs(r)) > options_.solve_tolf
         info(1) = 25;
         return
     end
@@ -505,21 +500,21 @@ j = fh_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr, T_or
 resids = r(eq_index);
 jac = j(eq_index,1:nvar);
 
-function [r, g1] = block_mfs_steadystate(y, fh_static, b, y_all, exo, params, T, M)
+function [r, g1] = block_mfs_steadystate(y, fh_static, b, y_all, exo, params, T, M_)
 % Wrapper around the static files, for block without bytecode
-mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
+mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
 y_all(mfs_idx) = y;
-[~,~,r,g1] = fh_static(y_all, exo, params, M.block_structure_stat.block(b).g1_sparse_rowval, ...
-                       M.block_structure_stat.block(b).g1_sparse_colval, ...
-                       M.block_structure_stat.block(b).g1_sparse_colptr, T);
+[~,~,r,g1] = fh_static(y_all, exo, params, M_.block_structure_stat.block(b).g1_sparse_rowval, ...
+                       M_.block_structure_stat.block(b).g1_sparse_colval, ...
+                       M_.block_structure_stat.block(b).g1_sparse_colptr, T);
 
-function [r, g1] = bytecode_steadystate(y, exo, params)
+function [r, g1] = bytecode_steadystate(y, exo, params, M_, options_)
 % Wrapper around the static file, for bytecode (without block)
-[r, g1] = bytecode(y, exo, params, y, 1, exo, 'evaluate', 'static');
+[r, g1] = bytecode(M_, options_, y, exo, params, y, 1, exo, 'evaluate', 'static');
 
-function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M)
-% Wrapper around the static files, for block without bytecode
-mfs_idx = M.block_structure_stat.block(b).variable(end-M.block_structure_stat.block(b).mfs+1:end);
+function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M_, options_)
+% Wrapper around the static files, for bytecode with block
+mfs_idx = M_.block_structure_stat.block(b).variable(end-M_.block_structure_stat.block(b).mfs+1:end);
 y_all(mfs_idx) = y;
-[r, g1] = bytecode(y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block = ' int2str(b) ]);
-g1 = g1(:,end-M.block_structure_stat.block(b).mfs+1:end); % Make Jacobian square if mfs>0
+[r, g1] = bytecode(M_, options_, y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', 'block_decomposed', ['block=' int2str(b) ]);
+g1 = g1(:,end-M_.block_structure_stat.block(b).mfs+1:end); % Make Jacobian square if mfs>0

matlab/evaluate_steady_state_file.m

-function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,steady_state_checkflag)
-% function [ys,params1,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,steady_state_checkflag)
+function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M_,options_,steady_state_checkflag)
+% function [ys,params1,info] = evaluate_steady_state_file(ys_init,exo_ss,M_,options_,steady_state_checkflag)
 % Evaluates steady state files
 %
 % INPUTS
 %   ys_init                   vector           initial values used to compute the steady
 %                                                 state
 %   exo_ss                    vector           exogenous steady state (incl. deterministic exogenous)
-%   M                         struct           model parameters
-%   options                   struct           options
+%   M_                        struct           model parameters
+%   options_                  struct           options
 %   steady_state_checkflag    boolean          indicator whether to check steady state returned
 % OUTPUTS
 %   ys                        vector           steady state
@@ -36,17 +36,15 @@ function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-params = M.params;
+params = M_.params;
 info = 0;
 
-fname = M.fname;
+fname = M_.fname;
 
-if options.steadystate_flag == 1
+if options_.steadystate_flag == 1
     % old format
-    assignin('base','tmp_00_',params);
-    evalin('base','M_.params=tmp_00_; clear(''tmp_00_'')');
     h_steadystate = str2func([fname '_steadystate']);
-    [ys,params1,check] = h_steadystate(ys_init, exo_ss,M,options);
+    [ys,params1,check] = h_steadystate(ys_init, exo_ss,M_,options_);
 else % steadystate_flag == 2
      % new format
     h_steadystate = str2func([fname '.steadystate']);
@@ -59,19 +57,19 @@ if check
     return
 end
 
-if M.param_nbr > 0
+if M_.param_nbr > 0
     updated_params_flag = max(abs(params1-params)) > 1e-12 ...
         || ~isequal(isnan(params1),isnan(params)); %checks whether numbers or NaN changed
 else
     updated_params_flag = 0;
 end
 
-h_set_auxiliary_variables = str2func([M.fname '.set_auxiliary_variables']);
+h_set_auxiliary_variables = str2func([M_.fname '.set_auxiliary_variables']);
 
 if  isnan(updated_params_flag) || (updated_params_flag  && any(isnan(params(~isnan(params))-params1(~isnan(params))))) %checks if new NaNs were added
     info(1) = 24;
     info(2) = NaN;
-    if M.set_auxiliary_variables
+    if M_.set_auxiliary_variables
         ys = h_set_auxiliary_variables(ys,exo_ss,params);
     end
     return
@@ -79,8 +77,12 @@ end
 
 if updated_params_flag && ~isreal(params1)
     info(1) = 23;
+    if ~isoctave
         info(2) = sum(imag(params).^2,'omitnan');
-    if M.set_auxiliary_variables
+    else
+        info(2) = nansum(imag(params).^2);
+    end
+    if M_.set_auxiliary_variables
         ys = h_set_auxiliary_variables(ys,exo_ss,params);
     end
     return
@@ -91,21 +93,16 @@ if updated_params_flag
 end
 
 % adding values for auxiliary variables
-if ~isempty(M.aux_vars) && ~options.ramsey_policy
-    if M.set_auxiliary_variables
+if ~isempty(M_.aux_vars) && ~options_.ramsey_policy
+    if M_.set_auxiliary_variables
         ys = h_set_auxiliary_variables(ys,exo_ss,params);
     end
 end
 
 if steady_state_checkflag
     % Check whether the steady state obtained from the _steadystate file is a steady state.
-    [residuals, check] = evaluate_static_model(ys, exo_ss, params, M, options);
-    if check
-        info(1) = 19;
-        info(2) = check; % to be improved
-        return
-    end
-    if max(abs(residuals)) > options.solve_tolf
+    residuals = evaluate_static_model(ys, exo_ss, params, M_, options_);
+    if max(abs(residuals)) > options_.solve_tolf
         info(1) = 19;
         info(2) = residuals'*residuals;
         return

matlab/forcst.m

-function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,options_)
-% function [yf,int_width,int_width_ME]=forecst(dr,y0,horizon,var_list,M_,oo_,options_)
+function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,options_)
+% function [yf,int_width,int_width_ME]=forecst(dr,y0,horizon,var_list,M_,options_)
 %   computes mean forecast for a given value of the parameters
 %   computes also confidence band for the forecast
 %
@@ -10,7 +10,6 @@ function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,option
 %   var_list:    list of variables (character matrix)
 %   M_:          Dynare model structure
 %   options_:    Dynare options structure
-%   oo_:         Dynare results structure
 
 % OUTPUTS:
 %   yf:          mean forecast
@@ -22,7 +21,7 @@ function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,option
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2003-2019 Dynare Team
+% Copyright © 2003-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,6 +38,10 @@ function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,option
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+if options_.order>1
+    error('forcst.m: Only order=1 is supported. Skipping computuations.')
+end
+
 yf = simult_(M_,options_,y0,dr,zeros(horizon,M_.exo_nbr),1);
 nstatic = M_.nstatic;
 nspred = M_.nspred;

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);

matlab/ftest.m

-function ftest (s1,s2)
-
-% Copyright © 2001-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/>.
-
-global nvx nvy x y lag1
-
-if size(s1,1) ~= 2
-    error ('Spécifiez deux fichiers pour la comparaison.') ;
-end
-
-for i = 1:2
-    if ~ isempty(find(abs(s1(i,:)) == 46))
-        error ('Entrez les noms de fichiers sans extensions.') ;
-    end
-end
-
-s1 = [s1 [' ';' ']] ;
-file1 = [s1(1,1:min(find(abs(s1(1,:)) == 32))-1) '.BIN'] ;
-file2 = [s1(2,1:min(find(abs(s1(2,:)) == 32))-1) '.BIN'] ;
-
-fid=fopen(file1,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag1 = fread(fid,4,'int') ;
-nvx = fread(fid,[n1,n3],'int') ;
-x = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvx = char(nvx) ;
-
-fid=fopen(file2,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag2 = fread(fid,4,'int') ;
-nvy = fread(fid,[n1,n3],'int') ;
-y = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvy = char(nvy) ;
-
-if size(x,1) ~= size(y,1)
-    error ('FTEST: The two files don''t have the same number of variables.');
-end
-
-for i = 1:size(x,1)
-    if ~ strcmp(nvx(i,:),nvy(i,:))
-        error ('FTEST: The two files don''t have the same  variables.') ;
-    end
-end
-
-if nnz(lag1 - lag2) > 0
-    error ('FTEST: Leads and lags aren''t the same in both files.') ;
-end
-
-j = zeros(size(s2,1),1);
-for i=1:size(s2,1)
-    k = strmatch(s2(i,:),nvx,'exact') ;
-    if isempty(k)
-        t = ['FTEST: Variable ' s2(i) 'doesn''t exist'] ;
-        error (t) ;
-    else
-        j(i) =k;
-    end
-end
-
-y = y(j,:) ;
-x = x(j,:) ;
-
-%06/18/01 MJ replaced beastr by strmatch

matlab/functionSignatures.json

+{
+  "_schemaVersion": "1.0.0",
+  "dynare":
+  {
+     "inputs":
+     [
+        {"name":"fname", "kind":"required", "type":["file=*.mod,*.dyn"], "purpose":"Model file"},
+        {"name":"options", "kind":"flag", "type":["char", "choices={'debug', 'noclearall', 'onlyclearglobals', 'savemacro[=macro_file]', 'onlymacro', 'linemacro', 'notmpterms', 'nolog', 'warn_uninit', 'console', 'nograph', 'nointeractive', 'parallel[=cluster_name]', 'conffile=parallel_config_path_and_filename', 'parallel_slave_open_mode', 'parallel_test', 'parallel_use_psexec=true|false', '-D<variable>[=<value>]', '-I/path', 'nostrict', 'stochastic', 'fast', 'minimal_workspace', 'compute_xrefs', 'output=second|third', 'params_derivs_order=0|1|2', 'transform_unary_ops', 'exclude_eqs=<equation_tag_list_or_file>', 'include_eqs=<equation_tag_list_or_file>', 'json=parse|check|transform|compute', 'jsonstdout', 'onlyjson', 'jsonderivsimple', 'nopathchange', 'nopreprocessoroutput', 'mexext=<extension>', 'matlabroot=<path>', 'onlymodel', 'notime', 'use_dll', 'nocommutativity'}"],  "repeating":true}
+     ]
+  },
+  "internals":
+  {
+     "inputs":
+     [
+        {"name":"flag", "kind":"required", "type":["char", "choices={'--test'}"], "purpose":"Test routine"},
+        {"name":"fname", "kind":"required", "type":["file=*.m"], "purpose":"Routine name"}
+     ]
+  }, 
+  "internals":
+  {
+     "inputs":
+     [
+        {"name":"flag", "kind":"required", "type":["char", "choices={'--display-mh-history','--load-mh-history'}"], "purpose":"Model file"},
+        {"name":"fname", "kind":"required", "type":["file=*.mod,*.dyn"], "purpose":"Model name"}
+     ]
+  }
+}
\ No newline at end of file

matlab/gcompare.m

-function gcompare(s1,s2)
-% GCOMPARE :    GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
-%               This optional command plots the trajectories of a list of
-%               variables in two different simulations. One plot is drawn
-%               for each variable. The trajectories must have been previously
-%               saved by the instruction DYNASAVE. The simulation in file1
-%               is refered to as the base simulation.
-
-% Copyright © 2001-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/>.
-
-global options_ M_
-global nvx nvy x y lag1
-
-ftest(s1,s2) ;
-
-ix = [1-lag1(1):size(x,2)-lag1(1)]' ;
-i = [lag1(1):size(ix,1)-lag1(2)+1]' ;
-
-if options_.smpl == 0
-    i = [M_.maximum_lag:size(y,2)]' ;
-else
-    i = [options_.smpl(1):options_.smpl(2)]' ;
-end
-
-for k = 1:size(x,1)
-    figure ;
-    plot (ix(i),x(k,i),ix(i),y(k,i)) ;
-    xlabel (['Periods']) ;
-    title (['Variable ' s2(k,:)]) ;
-    l = min(i) + 1;
-    ll = max(i) - 1 ;
-    text (l,x(k,l),s1(1,:)) ;
-    text (ll,y(k,ll),s1(2,:)) ;
-end
-
-% 06/18/01 MJ corrected treatment of options_.smpl
-% 06/24/01 MJ removed color specification