SMC: gracefully exit with unsupported options

(cherry picked from commit 5d47ac2aa9)

.gitlab-ci.yml

@@ -11,12 +11,10 @@ variables:
 # - if VERSION was already set (when manually running a pipeline), use it
 # - if we are in the official Dynare repository:
 #   + if on a tag: use the tag
-#   + if on master: use 7-unstable-$TIMESTAMP-$COMMIT
 #   + on another branch: use $BRANCH-$TIMESTAMP-$COMMIT
 # - if in a personal repository: use $USER-$TIMESTAMP-$COMMIT
 before_script:
   - 'if [[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]] && [[ -n $CI_COMMIT_TAG ]]; then export VERSION=$CI_COMMIT_TAG; fi'
-  - 'if [[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]] && [[ $CI_COMMIT_REF_NAME == master ]]; then export VERSION=7-unstable-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA; fi'
   - 'if [[ -z $VERSION ]] && [[ $CI_PROJECT_NAMESPACE == Dynare ]]; then export VERSION=$CI_COMMIT_REF_NAME-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA; fi'
   - 'if [[ -z $VERSION ]]; then export VERSION=$CI_PROJECT_NAMESPACE-$(date +%F-%H%M)-$CI_COMMIT_SHORT_SHA; fi'
 
@@ -162,7 +160,6 @@ test_old_matlab:
     paths:
       - build-old-matlab/meson-logs/testlog.txt
     when: always
-  when: manual
 
 test_octave:
   stage: test
@@ -175,7 +172,6 @@ test_octave:
       - build-octave/meson-logs/testlog.txt
     when: always
   needs: [ "build_octave" ]
-  when: manual
 
 test_clang_format:
   stage: test
@@ -191,7 +187,7 @@ test_clang_format:
 sign_windows:
   stage: sign
   rules:
-    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_REF_NAME == "master"'
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^6/'
       when: on_success
     - when: never
   tags:
@@ -205,10 +201,10 @@ sign_windows:
       - windows/exe-signed/*
     expire_in: 3 days
 
-deploy_manual_unstable:
+deploy_manual_stable:
   stage: deploy
   rules:
-    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_REF_NAME == "master"'
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^6\.[0-9]+$/'
       when: on_success
     - when: never
   tags:
@@ -216,12 +212,13 @@ deploy_manual_unstable:
   dependencies:
     - build_doc
   script:
-    - rsync --recursive --links --delete build-doc/dynare-manual.html/ /srv/www.dynare.org/manual-unstable/
+    - rsync --recursive --links --delete build-doc/dynare-manual.html/ /srv/www.dynare.org/manual/
+    - cp build-doc/dynare-manual.pdf /srv/www.dynare.org/manual.pdf
 
-deploy_snapshot_unstable:
+deploy_release_stable:
   stage: deploy
   rules:
-    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_REF_NAME == "master"'
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^6\.[0-9]+$/'
       when: on_success
     - when: never
   tags:
@@ -233,11 +230,33 @@ deploy_snapshot_unstable:
     - pkg_macOS_arm64
     - pkg_macOS_x86_64
   script:
-    - cp build-src/meson-dist/*.tar.xz /srv/www.dynare.org/snapshot/source/ && ln -sf *.tar.xz /srv/www.dynare.org/snapshot/source/dynare-latest-src.tar.xz
+    - cp build-src/meson-dist/*.tar.xz /srv/www.dynare.org/release/source/
-    - f=(windows/exe-signed/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows/dynare-latest-win.exe
+    - cp windows/exe-signed/* /srv/www.dynare.org/release/windows/
-    - f=(windows/7z/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows-7z/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows-7z/dynare-latest-win.7z
+    - cp windows/7z/* /srv/www.dynare.org/release/windows-7z/
-    - f=(windows/zip/*) && cp ${f[0]} /srv/www.dynare.org/snapshot/windows-zip/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/windows-zip/dynare-latest-win.zip
+    - cp windows/zip/* /srv/www.dynare.org/release/windows-zip/
-    - f=(macOS/pkg/*-arm64.pkg) && cp ${f[0]} /srv/www.dynare.org/snapshot/macos-arm64/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/macos-arm64/dynare-latest-macos-arm64.pkg
+    - cp macOS/pkg/*-arm64.pkg /srv/www.dynare.org/release/macos-arm64/
-    - f=(macOS/pkg/*-x86_64.pkg) && cp ${f[0]} /srv/www.dynare.org/snapshot/macos-x86_64/ && ln -sf ${f[0]##*/} /srv/www.dynare.org/snapshot/macos-x86_64/dynare-latest-macos-x86_64.pkg
+    - cp macOS/pkg/*-x86_64.pkg /srv/www.dynare.org/release/macos-x86_64/
-    - ~/update-snapshot-list.sh
+    - ~/update-release-list.sh
     - curl -X POST -F token="$WEBSITE_PIPELINE_TRIGGER_TOKEN" -F ref=master https://git.dynare.org/api/v4/projects/40/trigger/pipeline
+
+deploy_beta_stable:
+  stage: deploy
+  rules:
+    - if: '$CI_PROJECT_NAMESPACE == "Dynare" && $CI_COMMIT_TAG =~ /^6(\.[0-9]+)?-(beta|rc)[0-9]+$/'
+      when: on_success
+    - when: never
+  tags:
+    - deploy
+  dependencies:
+    - pkg_source
+    - pkg_windows
+    - sign_windows
+    - pkg_macOS_arm64
+    - pkg_macOS_x86_64
+  script:
+    - cp build-src/meson-dist/*.tar.xz /srv/www.dynare.org/beta/source/
+    - cp windows/exe-signed/* /srv/www.dynare.org/beta/windows/
+    - cp windows/7z/* /srv/www.dynare.org/beta/windows-7z/
+    - cp windows/zip/* /srv/www.dynare.org/beta/windows-zip/
+    - cp macOS/pkg/*-arm64.pkg /srv/www.dynare.org/beta/macos-arm64/
+    - cp macOS/pkg/*-x86_64.pkg /srv/www.dynare.org/beta/macos-x86_64/

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

@@ -7493,11 +7493,6 @@ observed variables.
                standard Random-Walk Metropolis-Hastings. Does not yet support
                ``moments_varendo``, ``bayesian_irf``, and ``smoother``.
 
-           ``'dsmh'``
-
-               Instructs Dynare to use the Dynamic Striated Metropolis Hastings
-               sampler proposed by *Waggoner, Wu and Zha (2016)* instead of the
-               standard Random-Walk Metropolis-Hastings.
 
     .. option:: posterior_sampler_options = (NAME, VALUE, ...)
 
@@ -16020,7 +16015,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``, ``amssymb``, 
+    ``booktabs``, ``caption``, ``float,`` ``amsmath``, ``amsfonts``,
     and ``morefloats``.
 
 

macOS/deps/Makefile

@@ -1,4 +1,4 @@
-# Copyright © 2019-2024 Dynare Team
+# Copyright © 2019-2023 Dynare Team
 #
 # This file is part of Dynare.
 #
@@ -53,7 +53,7 @@ clean-all: clean-lib clean-src clean-tar
 #
 tarballs/slicot-$(SLICOT_VERSION).tar.gz:
 	mkdir -p tarballs
-	wget $(WGET_OPTIONS) -O $@ https://deb.debian.org/debian/pool/main/s/slicot/slicot_$(SLICOT_VERSION).orig.tar.xz
+	wget $(WGET_OPTIONS) -O $@ https://deb.debian.org/debian/pool/main/s/slicot/slicot_$(SLICOT_VERSION).orig.tar.gz
 
 $(DEPS_ARCH)/sources64/slicot-$(SLICOT_VERSION): tarballs/slicot-$(SLICOT_VERSION).tar.gz
 	rm -rf $(DEPS_ARCH)/sources64/slicot-*
@@ -62,7 +62,7 @@ $(DEPS_ARCH)/sources64/slicot-$(SLICOT_VERSION): tarballs/slicot-$(SLICOT_VERSIO
 	touch $@
 
 $(DEPS_ARCH)/lib64/slicot/libslicot64_pic.a: $(DEPS_ARCH)/sources64/slicot-$(SLICOT_VERSION)
-	make -C $< -f makefile_Unix FORTRAN=$(BREWDIR)/bin/gfortran LOADER=$(BREWDIR)/bin/gfortran SLICOTLIB=../libslicot64_pic.a OPTS="-O3 -fdefault-integer-8" lib -j$(NTHREADS)
+	make -C $< FORTRAN=$(BREWDIR)/bin/gfortran LOADER=$(BREWDIR)/bin/gfortran SLICOTLIB=../libslicot64_pic.a OPTS="-O3 -fdefault-integer-8" lib -j$(NTHREADS)
 	strip -S $</libslicot64_pic.a
 	mkdir -p $(dir $@)
 	cp $</libslicot64_pic.a $@

macOS/deps/versions.mk

@@ -1,2 +1,2 @@
-SLICOT_VERSION = 5.9~20240205.gita037f7e
+SLICOT_VERSION = 5.0+20101122
 X13AS_VERSION = 1-1-b60

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=plot_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
+    dyn_graph=dynare_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,31 +183,3 @@ 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

matlab/+gsa/map_identification.m

@@ -272,8 +272,80 @@ elseif opt_gsa.morris==3
     return
 elseif opt_gsa.morris==2   % ISKREV stuff
     return
-else  
+else  % main effects analysis
-    error('gsa/map_identification: unsupported option morris=%u',opt_gsa.morris)
+    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)
 end
 
 function []=create_TeX_loader(options_,figpath,ifig_number,caption,label_name,scale_factor)

matlab/+identification/get_perturbation_params_derivs.m

@@ -93,8 +93,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr
 % This function calls
 %   * [fname,'.dynamic']
 %   * [fname,'.dynamic_params_derivs']
-%   * [fname,'.sparse.static_g1']
+%   * [fname,'.static']
-%   * [fname,'.sparse.static_g2']
 %   * [fname,'.static_params_derivs']
 %   * commutation
 %   * dyn_vech
@@ -110,7 +109,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr
 %   * sylvester3a
 %   * get_perturbation_params_derivs_numerical_objective
 % =========================================================================
-% Copyright © 2019-2024 Dynare Team
+% Copyright © 2019-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -455,13 +454,12 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
         error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
     end
     %% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss
-    [g1_static, T_order_static, T_static] = feval([fname,'.sparse.static_g1'], ys, exo_steady_state', params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order
+    [~, g1_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order
     rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order
     dys = -g1_static\rp_static; %use implicit function theorem (equation 5 of Ratto and Iskrev (2012) to compute [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of steady state for all endogenous variables analytically, note that dys is in declaration order
     d2ys = zeros(endo_nbr, param_nbr, param_nbr); %initialize in tensor notation, note that d2ys is only needed for d2flag, i.e. for g1pp
     if d2flag
-        g2_static_v = feval([fname,'.sparse.static_g2'], ys, exo_steady_state', params, T_order_static, T_static);
+        [~, ~, g2_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order
-        g2_static = build_two_dim_hessian(M_.static_g2_sparse_indices, g2_static_v, endo_nbr, endo_nbr); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order
         if order < 3
             [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements
             g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3
@@ -507,7 +505,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
     end
     %handling of steady state for nonstationary variables
     if any(any(isnan(dys)))
-        [U,T] = schur(full(g1_static));
+        [U,T] = schur(g1_static);
         e1 = abs(ordeig(T)) < qz_criterium-1;
         k = sum(e1);       % Number of non stationary variables.
                            % Number of stationary variables: n = length(e1)-k

matlab/+identification/run.m

@@ -736,7 +736,7 @@ if iload <=0
                 end
                 run_index = 0; % reset index
             end
-            if SampleSize > 1 && mod(iteration,3)
+            if SampleSize > 1
                 dyn_waitbar(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]);
             end
         end

matlab/backward/simul_static_model.m

@@ -15,7 +15,7 @@ function simulation = simul_static_model(samplesize, innovations)
 % [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided
 %     through the shocks block.
 
-% Copyright © 2019-2024 Dynare Team
+% Copyright © 2019-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -83,17 +83,12 @@ else
     oo_.exo_simul = Innovations;
 end
 
-static_resid = str2func(sprintf('%s.sparse.static_resid', M_.fname));
+staticmodel = str2func(sprintf('%s.static', M_.fname));
-static_g1 = str2func(sprintf('%s.sparse.static_g1', M_.fname));
-function [resid, g1] = staticmodel(y, x, params)
-    [resid, T_order, T] = static_resid(y, x, params);
-    g1 = static_g1(y, x, params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
-end
 
 % Simulations (call a Newton-like algorithm for each period).
 for t=1:samplesize
     y = zeros(M_.endo_nbr, 1);
-    [oo_.endo_simul(:,t), errorflag, ~, ~, errorcode] = dynare_solve(@staticmodel, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, oo_.exo_simul(t,:), M_.params);
+    [oo_.endo_simul(:,t), errorflag, ~, ~, errorcode] = dynare_solve(staticmodel, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, oo_.exo_simul(t,:), M_.params);
     if errorflag
         dprintf('simul_static_mode: Nonlinear solver failed with errorcode=%i in period %i.', errorcode, t)
         oo_.endo_simul(:,t) = nan;
@@ -108,6 +103,4 @@ if isdseries(innovations)
     initperiod = innovations.dates(1);
 end
 
-simulation = [dseries(ysim', initperiod, M_.endo_names(1:M_.orig_endo_nbr)), dseries(xsim, initperiod, M_.exo_names)];
+simulation = [dseries(ysim', initperiod, M_.endo_names(1:M_.orig_endo_nbr)), dseries(xsim, initperiod, M_.exo_names)];
-
-end

matlab/build_two_dim_hessian.m

@@ -1,58 +0,0 @@
-function H = build_two_dim_hessian(sparse_indices, g2_v, neq, nvar)
-% Creates a 2D Hessian (equations in rows, pairs of variables in columns),
-% given the output from the sparse {static,dynamic}_g2.m
-%
-% – sparse_indices is typically equal to M_.{static,dynamic}_g2_sparse_indices
-% – g2_v is the vector of non zero values returned by {static,dynamic}_g2.m
-% – neq is the number of equations (equal to number of rows of the output matrix)
-% – nvar is the number of variables (the output matrix will have nvar² columns)
-
-% 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/>.
-
-nnz = size(sparse_indices, 1);
-
-%% The g2_* arrays may be expanded if there are symmetric elements added
-g2_i = int32(zeros(nnz, 1));
-g2_j = int32(zeros(nnz, 1));
-
-next_sym_idx = nnz + 1; % Index of next symmetric element to be added
-
-for k = 1:length(g2_v)
-    eq = sparse_indices(k, 1);
-    var1 = sparse_indices(k, 2)-1;
-    var2 = sparse_indices(k, 3)-1;
-
-    g2_i(k) = eq;
-    g2_j(k) = var1 * nvar + var2 + 1;
-
-    %% Add symmetric elements, which are not included by sparse {static,dynamic}_g2.m
-    if var1 ~= var2
-        g2_i(next_sym_idx) = eq;
-        g2_j(next_sym_idx) = var2 * nvar + var1 + 1;
-        g2_v(next_sym_idx) = g2_v(k);
-        next_sym_idx = next_sym_idx + 1;
-    end
-end
-
-%% On MATLAB < R2020a, sparse() does not accept int32 indices
-if ~isoctave && matlab_ver_less_than('9.8')
-    g2_i = double(g2_i);
-    g2_j = double(g2_j);
-end
-
-H = sparse(g2_i, g2_j, g2_v, neq, nvar*nvar);

matlab/discretionary_policy/discretionary_policy_1.m

@@ -13,7 +13,7 @@ function [dr, info, params]=discretionary_policy_1(M_, options_, dr, endo_steady
 % - info          [integer]       scalar or vector, error code.
 % - params        [double]        vector of potentially updated parameters
 
-% Copyright © 2007-2024 Dynare Team
+% Copyright © 2007-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -49,9 +49,7 @@ else
 end
 params=M_.params;
 
-y = zeros(M_.endo_nbr,1);
+[~,Uy,W] = feval([M_.fname,'.objective.static'],zeros(M_.endo_nbr,1),[], M_.params);
-[Uy, T_order, T] = feval([M_.fname,'.objective.sparse.static_g1'], y, [], params, M_.objective_g1_sparse_rowval, M_.objective_g1_sparse_colval, M_.objective_g1_sparse_colptr);
-
 if any(any(isnan(Uy)))
     info = 64 ; %the derivatives of the objective function contain NaN
     return;
@@ -72,8 +70,6 @@ if any(any(Uy~=0))
     return;
 end
 
-g2_v = feval([M_.fname,'.objective.sparse.static_g2'], y, [], params, T_order, T);
-W = build_two_dim_hessian(M_.objective_g2_sparse_indices, g2_v, 1, M_.endo_nbr);
 W=reshape(W,M_.endo_nbr,M_.endo_nbr);
 
 klen = M_.maximum_lag + M_.maximum_lead + 1;
@@ -126,4 +122,4 @@ dr.ghu=G(dr.order_var,:);
 if M_.maximum_endo_lag
     Selection=M_.lead_lag_incidence(1,dr.order_var)>0;%select state variables
 end
-dr.ghx=T(:,Selection);
+dr.ghx=T(:,Selection);

matlab/dynare.m

@@ -122,7 +122,7 @@ if isempty(dot_location)
     fnamelength = length(fname);
     fname1 = [fname '.dyn'];
     d = dir(fname1);
-    if isempty(d)
+    if length(d) == 0
         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 contains(fname,filesep)
+if ~isempty(strfind(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,13 +297,17 @@ 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(1:end-4) '/driver'])
+clear(['+' fname '/driver'])
 
 try
     cTic = tic;
-    evalin('base',[fname(1:end-4) '.driver']);
+    evalin('base',[fname '.driver']);
     cToc = toc(cTic);
     if nargout
         DynareInfo.time.compute = cToc;
@@ -311,7 +315,7 @@ try
 catch ME
     W = evalin('caller','whos');
     diary off
-    if ismember(fname(1:end-4),{W(:).name})
+    if ismember(fname,{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)

matlab/dynare_graph.m

@@ -0,0 +1,54 @@
+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

matlab/dynare_graph_close.m

@@ -0,0 +1,29 @@
+function dynare_graph_close()
+% function dynare_graph_close()
+% close a figure
+%
+% INPUT
+%   none
+%
+% OUTPUT
+%   none
+%
+% 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/>.

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 ~isempty(title)
+if length(title) > 0
     fprintf('\n\n%s\n',title);
 end
 
@@ -72,7 +72,7 @@ if nargin==9
     disp(optional_header)
 end
 
-if ~isempty(headers)
+if length(headers) > 0
     hh_tbl = sprintf(label_format_leftbound , headers{1});
     for i=2:length(headers)
         hh_tbl  = [hh_tbl sprintf(header_string_format, headers{i})];

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 negative of the posterior kernel of a DSGE model using the specified
+% Evaluates 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 minus the likelihood or posterior kernel.
+% - fval                    [double]        scalar, value of 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-2024 Dynare Team
+% Copyright © 2004-2023 Dynare Team
 %
 % This file is part of Dynare.
 %

matlab/estimation/smc/dsmh.m

@@ -1,299 +0,0 @@
-function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
-
-% Dynamic Striated Metropolis-Hastings algorithm.
-%
-% INPUTS
-%   o TargetFun  [char]     string specifying the name of the objective
-%                           function (posterior kernel).
-%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
-%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
-%   o dataset_              data structure
-%   o dataset_info          dataset info structure
-%   o options_              options structure
-%   o M_                    model structure
-%   o estim_params_         estimated parameters structure
-%   o bayestopt_            estimation options structure
-%   o oo_                   outputs structure
-%
-% SPECIAL REQUIREMENTS
-%   None.
-%
-% PARALLEL CONTEXT
-% The most computationally intensive part of this function may be executed
-% in parallel. The code suitable to be executed in
-% parallel on multi core or cluster machine (in general a 'for' cycle)
-% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
-%
-% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
-% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
-% to manage the parallel computations.
-%
-% This function was the first function to be parallelized. Later, other
-% functions have been parallelized using the same methodology.
-% Then the comments write here can be used for all the other pairs of
-% parallel functions and also for management functions.
-
-% 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_.posterior_sampler_options.dsmh;
-
-lambda = exp(bsxfun(@minus,options_.posterior_sampler_options.dsmh.H,1:1:options_.posterior_sampler_options.dsmh.H)/(options_.posterior_sampler_options.dsmh.H-1)*log(options_.posterior_sampler_options.dsmh.lambda1));
-c = 0.055 ;
-MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/10) ;
-
-% Step 0: Initialization of the sampler
-[param, tlogpost_iminus1, loglik, bayestopt_] = ...
-    smc_samplers_initialization(TargetFun, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
-
-ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
-zhat = 1 ;
-
-% The DSMH starts here
-for i=2:options_.posterior_sampler_options.dsmh.H
-    disp('');
-    disp('Tempered iteration');
-    disp(i) ;
-    % Step 1: sort the densities and compute IS weigths
-    [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
-    [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS) ;
-    % Step 2: tune c_i
-    c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-    % Step 3: Metropolis step
-    [param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-end
-
-weights = exp(loglik*(lambda(end)-lambda(end-1)));
-weights = weights/sum(weights);
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.particles);
-distrib_param = param(:,indx_resmpl);
-
-mean_xparam = mean(distrib_param,2);
-npar  = length(xparam1);
-lb95_xparam = zeros(npar,1) ;
-ub95_xparam = zeros(npar,1) ;
-for i=1:npar
-    temp = sortrows(distrib_param(i,:)') ;
-    lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.particles) ;
-    ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.particles) ;
-end
-
-TeX = options_.TeX;
-
-str = sprintf(' Param. \t Lower Bound (95%%) \t Mean \t Upper Bound (95%%)');
-for l=1:npar
-    name = get_the_name(l,TeX,M_,estim_params_,options_.varobs);
-    str = sprintf('%s\n %s \t\t %5.4f \t\t %7.5f \t\t %5.4f', str, name, lb95_xparam(l), mean_xparam(l), ub95_xparam(l));
-end
-disp(str)
-disp('')
-
-%% Plot parameters densities
-
-if TeX
-    fidTeX = fopen([M_.fname '_param_density.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by DSMH.m (Dynare).\n');
-    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
-    fprintf(fidTeX,' \n');
-end
-
-number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
-bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
-kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourier Transform approximation.
-
-plt = 1 ;
-%for plt = 1:nbplt,
-if TeX
-    NAMES = [];
-    TeXNAMES = [];
-end
-hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
-for k=1:npar %min(nstar,npar-(plt-1)*nstar)
-    subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
-    %kk = (plt-1)*nstar+k;
-    [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
-    optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.particles,bandwidth,kernel_function);
-    [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
-        options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
-    plot(density(:,1),density(:,2));
-    hold on
-    if TeX
-        title(texname,'interpreter','latex')
-    else
-        title(name,'interpreter','none')
-    end
-    hold off
-    axis tight
-    drawnow
-end
-dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
-if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    % TeX eps loader file
-    fprintf(fidTeX,'\\begin{figure}[H]\n');
-    fprintf(fidTeX,'\\centering \n');
-    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/floor(sqrt(npar)),1),M_.fname,int2str(plt));
-    fprintf(fidTeX,'\\caption{Parameter densities based on the Dynamic Striated Metropolis-Hastings algorithm.}');
-    fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
-    fprintf(fidTeX,'\\end{figure}\n');
-    fprintf(fidTeX,' \n');
-end
-%end
-
-function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param)
-[~,indx_ord] = sortrows(tlogpost_iminus1);
-tlogpost_iminus1 = tlogpost_iminus1(indx_ord);
-param = param(:,indx_ord);
-loglik = loglik(indx_ord);
-
-function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS)
-if i==1
-    tlogpost_i = tlogpost_iminus1 + loglik*lambda(i);
-else
-    tlogpost_i = tlogpost_iminus1 + loglik*(lambda(i)-lambda(i-1));
-end
-weights = exp(tlogpost_i-tlogpost_iminus1);
-zhat = (mean(weights))*zhat ;
-weights = weights/sum(weights);
-ESS(i) = 1/sum(weights.^2);
-% estimates of mean and variance
-mu = param*weights;
-z = bsxfun(@minus,param,mu);
-Omega = z*diag(weights)*z';
-Omegachol = chol(Omega)';
-
-function c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-disp('tuning c_i...');
-disp('Initial value =');
-disp(c) ;
-npar = size(param,1);
-lower_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^5;
-upper_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^(1/5);
-stop=0 ;
-while stop==0
-    acpt = 0.0;
-    indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
-    param0 = param(:,indx_resmpl);
-    tlogpost0 = tlogpost_i(indx_resmpl);
-    for j=1:options_.posterior_sampler_options.dsmh.G
-        for l=1:options_.posterior_sampler_options.dsmh.K
-            validate = 0;
-            while validate == 0
-                candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
-                if all(candidate >= mh_bounds.lb) && all(candidate <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-                    if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
-                        validate = 1;
-                        if rand(1,1)<exp(tlogpostx-tlogpost0(j)) % accept
-                            acpt = acpt + 1/(options_.posterior_sampler_options.dsmh.G*options_.posterior_sampler_options.dsmh.K);
-                            param0(:,j)= candidate;
-                            tlogpost0(j) = tlogpostx;
-                        end
-                    end
-                end
-            end
-        end
-    end
-    disp('Acceptation rate =') ;
-    disp(acpt) ;
-    if options_.posterior_sampler_options.dsmh.alpha0<=acpt && acpt<=options_.posterior_sampler_options.dsmh.alpha1
-        disp('done!');
-        stop=1;
-    else
-        if acpt<lower_prob
-            c = c/5;
-        elseif lower_prob<=acpt && acpt<=upper_prob
-            c = c*log(.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))/log(acpt);
-        else
-            c = 5*c;
-        end
-        disp('Trying with c= ') ;
-        disp(c)
-    end
-end
-
-function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-indx_levels = (1:1:MM-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM;
-npar = size(param,1) ;
-p = 1/(10*options_.posterior_sampler_options.dsmh.tau);
-disp('Metropolis step...');
-% build the dynamic grid of levels
-levels = [0.0;tlogpost_iminus1(indx_levels)];
-% initialize the outputs
-out_param = param;
-out_tlogpost_iminus1 = tlogpost_i;
-out_loglik = loglik;
-% resample and initialize the starting groups
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
-param0 = param(:,indx_resmpl);
-tlogpost_iminus10 = tlogpost_iminus1(indx_resmpl);
-tlogpost_i0 = tlogpost_i(indx_resmpl);
-loglik0 = loglik(indx_resmpl);
-% Start the Metropolis
-for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.tau
-    for j=1:options_.posterior_sampler_options.dsmh.G
-        u1 = rand(1,1);
-        u2 = rand(1,1);
-        if u1<p
-            k=1 ;
-            for m=1:MM-1
-                if levels(m)<=tlogpost_iminus10(j) && tlogpost_iminus10(j)<levels(m+1)
-                    k = m+1;
-                    break
-                end
-            end
-            indx = floor( (k-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM+1 + u2*(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM-1) );
-            if i==1
-                alp = (loglik(indx)-loglik0(j))*lambda(i);
-            else
-                alp = (loglik(indx)-loglik0(j))*(lambda(i)-lambda(i-1));
-            end
-            if u2<exp(alp)
-                param0(:,j) = param(:,indx);
-                tlogpost_i0(j) = tlogpost_i(indx);
-                loglik0(j) = loglik(indx);
-                tlogpost_iminus10(j) = tlogpost_iminus1(indx);
-            end
-        else
-            validate= 0;
-            while validate==0
-                candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
-                if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-                    if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
-                        validate = 1;
-                        if u2<exp(tlogpostx-tlogpost_i0(j)) % accept
-                            param0(:,j) = candidate;
-                            tlogpost_i0(j) = tlogpostx;
-                            loglik0(j) = loglikx;
-                            if i==1
-                                tlogpost_iminus10(j) = tlogpostx-loglikx*lambda(i);
-                            else
-                                tlogpost_iminus10(j) = tlogpostx-loglikx*(lambda(i)-lambda(i-1));
-                            end
-                        end
-                    end
-                end
-            end
-        end
-    end
-    if mod(l,options_.posterior_sampler_options.dsmh.tau)==0
-        rang = (l/options_.posterior_sampler_options.dsmh.tau-1)*options_.posterior_sampler_options.dsmh.G+1:l*options_.posterior_sampler_options.dsmh.G/options_.posterior_sampler_options.dsmh.tau;
-        out_param(:,rang) = param0;
-        out_tlogpost_iminus1(rang) = tlogpost_i0;
-        out_loglik(rang) = loglik0;
-    end
-end

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,amssymb}');
+fprintf(fid,'%s \n','\usepackage{amsmath,amsfonts}');
 fprintf(fid,'%s \n','\usepackage{breqn}');
 fprintf(fid,'%s \n','\usepackage{float,morefloats,caption}');
 fprintf(fid,'%s \n','\begin{document}');

matlab/optimal_policy/evaluate_planner_objective.m

@@ -59,7 +59,7 @@ function planner_objective_value = evaluate_planner_objective(M_,options_,oo_)
 
 % In the deterministic case, resorting to approximations for welfare is no longer required as it is possible to simulate the model given initial conditions for pre-determined variables and terminal conditions for forward-looking variables, whether these initial and terminal conditions are explicitly or implicitly specified. Assuming that the number of simulated periods is high enough for the new steady-state to be reached, the new unconditional welfare is thus the last period's welfare. As for the conditional welfare, it can be derived using backward recursions on the equation W = U + beta*W(+1) starting from the final unconditional steady-state welfare.
 
-% Copyright © 2007-2024 Dynare Team
+% Copyright © 2007-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -92,11 +92,11 @@ end
 
 if options_.ramsey_policy && oo_.gui.ran_perfect_foresight
     T = size(oo_.endo_simul,2);
-    U_term = feval([M_.fname '.objective.sparse.static_resid'], oo_.endo_simul(:,T-M_.maximum_lead), oo_.exo_simul(T-M_.maximum_lead,:), M_.params);
+    [U_term] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,T-M_.maximum_lead),oo_.exo_simul(T-M_.maximum_lead,:), M_.params);
     EW = U_term/(1-beta);
     W = EW;
     for t=T-M_.maximum_lead:-1:1+M_.maximum_lag
-        U = feval([M_.fname '.objective.sparse.static_resid'], oo_.endo_simul(:,t), oo_.exo_simul(t,:), M_.params);
+        [U] = feval([M_.fname '.objective.static'],oo_.endo_simul(:,t),oo_.exo_simul(t,:), M_.params);
         W = U + beta*W;
     end
     planner_objective_value = struct('conditional', W, 'unconditional', EW);
@@ -108,8 +108,7 @@ else
         ys = oo_.dr.ys;
     end
     if options_.order == 1 && ~options_.discretionary_policy
-        [U, T_order, T] = feval([M_.fname '.objective.sparse.static_resid'], ys, zeros(1,exo_nbr), M_.params);
+        [U,Uy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-        Uy = feval([M_.fname '.objective.sparse.static_g1'], ys, zeros(1,exo_nbr), M_.params, M_.objective_g1_sparse_rowval, M_.objective_g1_sparse_colval, M_.objective_g1_sparse_colptr, T_order, T);
 
         Gy = dr.ghx(nstatic+(1:nspred),:);
         Gu = dr.ghu(nstatic+(1:nspred),:);
@@ -138,9 +137,7 @@ else
         planner_objective_value.conditional.zero_initial_multiplier = W_L_0;
 
     elseif options_.order == 2 && ~M_.hessian_eq_zero %full second order approximation
-        [U, T_order, T] = feval([M_.fname '.objective.sparse.static_resid'], ys, zeros(1,exo_nbr), M_.params);
+        [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-        [Uy, T_order, T] = feval([M_.fname '.objective.sparse.static_g1'], ys, zeros(1,exo_nbr), M_.params, M_.objective_g1_sparse_rowval, M_.objective_g1_sparse_colval, M_.objective_g1_sparse_colptr, T_order, T);
-        Uyy_v = feval([M_.fname '.objective.sparse.static_g2'], ys, zeros(1,exo_nbr), M_.params, T_order, T);
 
         Gy = dr.ghx(nstatic+(1:nspred),:);
         Gu = dr.ghu(nstatic+(1:nspred),:);
@@ -156,7 +153,6 @@ else
         guu(dr.order_var,:) = dr.ghuu;
         gss(dr.order_var,:) = dr.ghs2;
 
-        Uyy = build_two_dim_hessian(M_.objective_g2_sparse_indices, Uyy_v, 1, M_.endo_nbr);
         Uyy = full(Uyy);
 
         Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
@@ -232,16 +228,13 @@ else
         planner_objective_value.conditional.steady_initial_multiplier = W_L_SS;
         planner_objective_value.conditional.zero_initial_multiplier = W_L_0;
     elseif (options_.order == 2 && M_.hessian_eq_zero) || options_.discretionary_policy %linear quadratic problem
-        [U, T_order, T] = feval([M_.fname '.objective.sparse.static_resid'], ys, zeros(1,exo_nbr), M_.params);
+        [U,Uy,Uyy] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
-        [Uy, T_order, T] = feval([M_.fname '.objective.sparse.static_g1'], ys, zeros(1,exo_nbr), M_.params, M_.objective_g1_sparse_rowval, M_.objective_g1_sparse_colval, M_.objective_g1_sparse_colptr, T_order, T);
-        Uyy_v = feval([M_.fname '.objective.sparse.static_g2'], ys, zeros(1,exo_nbr), M_.params, T_order, T);
 
         Gy = dr.ghx(nstatic+(1:nspred),:);
         Gu = dr.ghu(nstatic+(1:nspred),:);
         gy(dr.order_var,:) = dr.ghx;
         gu(dr.order_var,:) = dr.ghu;
 
-        Uyy = build_two_dim_hessian(M_.objective_g2_sparse_indices, Uyy_v, 1, M_.endo_nbr);
         Uyy = full(Uyy);
 
         Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
@@ -312,7 +305,7 @@ else
             dr.(['g_' num2str(i)]) = [dr.(['g_' num2str(i)]); W.(['W_' num2str(i)])];
         end
         % Amends the steady-state vector accordingly
-        U = feval([M_.fname '.objective.sparse.static_resid'], ys, zeros(1,exo_nbr), M_.params);
+        [U] = feval([M_.fname '.objective.static'],ys,zeros(1,exo_nbr), M_.params);
         ysteady = [ys(oo_.dr.order_var); U/(1-beta)];
 
         % Generates the sequence of shocks to compute unconditional welfare

matlab/optimization/solve1.m

@@ -123,11 +123,7 @@ for its = 1:maxit
         fjac2=fjac'*fjac;
         temp=max(sum(abs(fjac2)));
         if temp>0
-            if issparse(fjac)
+            p=-(fjac2+sqrt(nn*eps)*temp*eye(nn))\(fjac'*fvec);
-                p=-(fjac2+sqrt(nn*eps)*temp*speye(nn))\(fjac'*fvec);
-            else
-                p=-(fjac2+sqrt(nn*eps)*temp*eye(nn))\(fjac'*fvec);
-            end
         else
             errorflag = true;
             errorcode = 5;

matlab/partial_information/PCL_resol.m

@@ -25,7 +25,7 @@ function [dr,info]=PCL_resol(ys,check_flag)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2001-2024 Dynare Team
+% Copyright © 2001-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -65,13 +65,7 @@ end
 dr.ys = ys;
 check1 = 0;
 % testing for steadystate file
-static_resid = str2func(sprintf('%s.sparse.static_resid', M_.fname));
+fh = str2func([M_.fname '.static']);
-static_g1 = str2func(sprintf('%s.sparse.static_g1', M_.fname));
-function [resid, g1] = static_resid_g1(y, x, params)
-    [resid, T_order, T] = static_resid(y, x, params);
-    g1 = static_g1(y, x, params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
-end
-
 if options_.steadystate_flag
     [dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,...
                            [oo_.exo_steady_state; ...
@@ -93,17 +87,17 @@ else
     if ~options_.ramsey_policy
         if options_.linear == 0
             % nonlinear models
-            if max(abs(static_resid_g1(dr.ys, [oo_.exo_steady_state; ...
+            if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; ...
-                                               oo_.exo_det_steady_state], M_.params))) > options_.solve_tolf
+                                    oo_.exo_det_steady_state], M_.params))) > options_.solve_tolf
                 opt = options_;
                 opt.jacobian_flag = false;
-                [dr.ys, check1] = dynare_solve(static_resid_g1, dr.ys, options.steady_.maxit, options_.solve_tolf, options_.solve_tolx, ...
+                [dr.ys, check1] = dynare_solve(fh,dr.ys, options.steady_.maxit, options_.solve_tolf, options_.solve_tolx, ...
                                                opt, [oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params);
             end
         else
             % linear models
-            [fvec,jacob] = static_resid_g1(dr.ys, [oo_.exo_steady_state;...
+            [fvec,jacob] = feval(fh,dr.ys,[oo_.exo_steady_state;...
-                                                   oo_.exo_det_steady_state], M_.params);
+                                oo_.exo_det_steady_state], M_.params);
             if max(abs(fvec)) > 1e-12
                 dr.ys = dr.ys-jacob\fvec;
             end
@@ -117,7 +111,7 @@ if check1
         resid = check1 ;
     else
         info(1)= 20;
-        resid = static_resid(ys, oo_.exo_steady_state, M_.params);
+        resid = feval(fh,ys,oo_.exo_steady_state, M_.params);
     end
     info(2) = resid'*resid ;
     return
@@ -146,8 +140,6 @@ end
 oo_.exo_simul = tempex;
 tempex = [];
 
-end
-
 % 01/01/2003 MJ added dr_algo == 1
 % 08/24/2001 MJ uses Schmitt-Grohe and Uribe (2001) constant correction
 %               in dr.ghs2

matlab/partial_information/add_auxiliary_variables_to_steadystate.m

@@ -2,7 +2,7 @@ function ys1 = add_auxiliary_variables_to_steadystate(ys,aux_vars,fname, ...
                                                   exo_steady_state, exo_det_steady_state,params, byte_code)
 % Add auxiliary variables to the steady state vector
 
-% Copyright © 2009-2024 Dynare Team
+% Copyright © 2009-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,7 +30,7 @@ for i=1:n+1
                        [exo_steady_state; ...
                         exo_det_steady_state],params);
     else
-        res = feval([fname '.sparse.static_resid'], ys1,...
+        res = feval([fname '.static'],ys1,...
                     [exo_steady_state; ...
                      exo_det_steady_state],params);
     end

matlab/perfect-foresight-models/perfect_foresight_solver.m

@@ -23,7 +23,7 @@ function [oo_, ts]=perfect_foresight_solver(M_, options_, oo_, no_error_if_learn
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 1996-2024 Dynare Team
+% Copyright © 1996-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -55,7 +55,7 @@ end
 periods = options_.periods;
 
 if options_.debug
-    model_static = str2func([M_.fname,'.sparse.static_resid']);
+    model_static = str2func([M_.fname,'.static']);
     for ii=1:size(oo_.exo_simul,1)
         [residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params);
     end

meson.build

@@ -1,6 +1,6 @@
 project('dynare',
         'cpp', 'fortran', 'c',
-        version : '7-unstable',
+        version : '6.1',
         # NB: update C++ standard in .clang-format whenever the following is modified
         default_options : [ 'cpp_std=gnu++20', 'fortran_std=f2018',
                             'c_std=gnu17', 'warning_level=2' ],
@@ -365,10 +365,6 @@ shared_module('logarithmic_reduction', [ 'mex/sources/logarithmic_reduction/mexF
 shared_module('disclyap_fast', [ 'mex/sources/disclyap_fast/disclyap_fast.f08' ] + mex_blas_fortran_iface,
               kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ])
 
-# TODO: Same remark as A_times_B_kronecker_C
-shared_module('riccati_update', [ 'mex/sources/riccati_update/mexFunction.f08' ] + mex_blas_fortran_iface,
-              kwargs : mex_kwargs, dependencies : [ blas_dep, lapack_dep ])
-
 qmc_sequence_src = [ 'mex/sources/sobol/qmc_sequence.cc',
                      'mex/sources/sobol/sobol.f08' ]
 # Hack for statically linking libgfortran
@@ -1830,7 +1826,6 @@ mod_and_m_tests = [
                 'solver-test-functions/wood.m',] },
   { 'test' : [ 'cyclereduction.m' ] },
   { 'test' : [ 'logarithmicreduction.m' ] },
-  { 'test' : [ 'riccatiupdate.m' ] },
   { 'test' : [ 'kalman/likelihood/test_kalman_mex.m' ] },
   { 'test' : [ 'contribs.m' ],
     'extra' : [ 'sandbox.mod',

mex/sources/k_order_welfare/k_ord_objective.cc

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021-2024 Dynare Team
+ * Copyright © 2021-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -18,23 +18,26 @@
  */
 
 #include "k_ord_objective.hh"
+#include "objective_abstract_class.hh"
 
 #include <cassert>
 #include <utility>
 
-KordwDynare::KordwDynare(KordpDynare& m, Journal& jr, Vector& inParams,
+KordwDynare::KordwDynare(KordpDynare& m, ConstVector& NNZD_arg, Journal& jr, Vector& inParams,
-                         std::unique_ptr<ObjectiveMFile> objectiveFile_arg,
+                         std::unique_ptr<ObjectiveAC> objectiveFile_arg,
                          const std::vector<int>& dr_order) :
     model {m},
+    NNZD {NNZD_arg},
     journal {jr},
     params {inParams},
     resid(1),
     ud {1},
     objectiveFile {std::move(objectiveFile_arg)}
 {
+  dynppToDyn = dr_order;
   dynToDynpp.resize(model.ny());
   for (int i = 0; i < model.ny(); i++)
-    dynToDynpp[dr_order[i]] = i;
+    dynToDynpp[dynppToDyn[i]] = i;
 }
 
 void
@@ -43,10 +46,71 @@ KordwDynare::calcDerivativesAtSteady()
 
   assert(ud.begin() == ud.end());
 
+  std::vector<TwoDMatrix> dyn_ud;     // Planner's objective derivatives, in Dynare form
+  dyn_ud.emplace_back(1, model.ny()); // Allocate Jacobian
+  dyn_ud.back().zeros();
+
+  for (int i = 2; i <= model.order(); i++)
+    {
+      // Higher order derivatives, as sparse (3-column) matrices
+      dyn_ud.emplace_back(static_cast<int>(NNZD[i - 1]), 3);
+      dyn_ud.back().zeros();
+    }
+
   Vector xx(model.nexog());
   xx.zeros();
   resid.zeros();
-  objectiveFile->eval(model.getSteady(), xx, params, resid, dynToDynpp, ud);
+  objectiveFile->eval(model.getSteady(), xx, params, resid, dyn_ud);
+
+  for (int i = 1; i <= model.order(); i++)
+    populateDerivativesContainer(dyn_ud, i);
+}
+
+void
+KordwDynare::populateDerivativesContainer(const std::vector<TwoDMatrix>& dyn_ud, int ord)
+{
+  const TwoDMatrix& u = dyn_ud[ord - 1];
+
+  // utility derivatives FSSparseTensor instance
+  auto udTi = std::make_unique<FSSparseTensor>(ord, model.ny(), 1);
+
+  IntSequence s(ord, 0);
+
+  if (ord == 1)
+    for (int i = 0; i < u.ncols(); i++)
+      {
+        for (int j = 0; j < u.nrows(); j++)
+          {
+            double x = u.get(j, dynppToDyn[s[0]]);
+            if (x != 0.0)
+              udTi->insert(s, j, x);
+          }
+        s[0]++;
+      }
+  else // ord ≥ 2
+    for (int i = 0; i < u.nrows(); i++)
+      {
+        int j = static_cast<int>(u.get(i, 0)) - 1;
+        int i1 = static_cast<int>(u.get(i, 1)) - 1;
+        if (j < 0 || i1 < 0)
+          continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
+
+        for (int k = 0; k < ord; k++)
+          {
+            s[k] = dynToDynpp[i1 % model.ny()];
+            i1 /= model.ny();
+          }
+
+        if (ord == 2 && !s.isSorted())
+          continue; // Skip symmetric elements (only needed at order 2)
+        else if (ord > 2)
+          s.sort(); // For higher order, canonicalize the multi-index
+
+        double x = u.get(i, 2);
+        udTi->insert(s, j, x);
+      }
+
+  ud.insert(std::move(udTi));
 }
 
 template<>

mex/sources/k_order_welfare/k_ord_objective.hh

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021-2024 Dynare Team
+ * Copyright © 2021-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,7 +21,7 @@
 #define K_ORD_OBJECTIVE_HH
 
 #include "k_ord_dynare.hh"
-#include "objective_m.hh"
+#include "objective_abstract_class.hh"
 
 class KordwDynare;
 
@@ -29,19 +29,22 @@ class KordwDynare
 {
 public:
   KordpDynare& model;
+  const ConstVector& NNZD;
 
 private:
   Journal& journal;
   Vector& params;
   Vector resid;
   TensorContainer<FSSparseTensor> ud; // planner's objective derivatives, in Dynare++ form
+  std::vector<int> dynppToDyn;        // Maps Dynare++ jacobian variable indices to Dynare ones
   std::vector<int> dynToDynpp;        // Maps Dynare jacobian variable indices to Dynare++ ones
-  std::unique_ptr<ObjectiveMFile> objectiveFile;
+  std::unique_ptr<ObjectiveAC> objectiveFile;
 
 public:
-  KordwDynare(KordpDynare& m, Journal& jr, Vector& inParams,
+  KordwDynare(KordpDynare& m, ConstVector& NNZD_arg, Journal& jr, Vector& inParams,
-              std::unique_ptr<ObjectiveMFile> objectiveFile_arg, const std::vector<int>& varOrder);
+              std::unique_ptr<ObjectiveAC> objectiveFile_arg, const std::vector<int>& varOrder);
   void calcDerivativesAtSteady();
+  void populateDerivativesContainer(const std::vector<TwoDMatrix>& dyn_ud, int ord);
   [[nodiscard]] const TensorContainer<FSSparseTensor>&
   getPlannerObjDerivatives() const
   {

mex/sources/k_order_welfare/k_order_welfare.cc

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021-2024 Dynare Team
+ * Copyright © 2021-2022 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -31,7 +31,6 @@
 
 #include <algorithm>
 #include <cassert>
-#include <string>
 
 #include "dynmex.h"
 
@@ -88,8 +87,8 @@ The routine proceeds in several steps:
    to the approxAtSteady method in the ApproximationWelfare class carries out the necessary
    operation
    - Importing the derivatives of the felicity function with the calcDerivativesAtSteady() method of
-     the KordwDynare class. It relies on the MATLAB-generated files, which are handled by the
+     the KordwDynare class. It relies on the Matlab-generated files, which are handled by the
-     ObjectiveMFile class
+     ObjectiveAC and ObjectiveMFile classes
    - Pinpointing the derivatives of the felicity and welfare functions. The performStep method of
      the KOrderWelfare class carries out the calculations,resorting to the FaaDiBruno class and its
      methods to get the needed intermediary results.
@@ -215,6 +214,15 @@ extern "C"
       mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you "
                    "used the right order option inside the `stoch_simul' command");
 
+    const mxArray* nnzderivatives_obj_mx = mxGetField(M_mx, 0, "NNZDerivatives_objective");
+    if (!(nnzderivatives_obj_mx && mxIsDouble(nnzderivatives_obj_mx)
+          && !mxIsComplex(nnzderivatives_obj_mx) && !mxIsSparse(nnzderivatives_obj_mx)))
+      mexErrMsgTxt("M_.NNZDerivatives should be a real dense array");
+    ConstVector NNZD_obj {nnzderivatives_obj_mx};
+    if (NNZD.length() < kOrder || NNZD_obj[kOrder - 1] == -1)
+      mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you "
+                   "used the right order option inside the `stoch_simul' command");
+
     const mxArray* endo_names_mx = mxGetField(M_mx, 0, "endo_names");
     if (!(endo_names_mx && mxIsCell(endo_names_mx)
           && mxGetNumberOfElements(endo_names_mx) == static_cast<size_t>(nEndo)))
@@ -253,32 +261,6 @@ extern "C"
     std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx) + nEndo, dr_order.begin(),
                    [](double x) { return static_cast<int>(x) - 1; });
 
-    const mxArray* objective_g1_sparse_rowval_mx
-        = mxGetField(M_mx, 0, "objective_g1_sparse_rowval");
-    if (!(objective_g1_sparse_rowval_mx && mxIsInt32(objective_g1_sparse_rowval_mx)))
-      mexErrMsgTxt("M_.objective_g1_sparse_rowval should be an int32 array");
-
-    const mxArray* objective_g1_sparse_colval_mx
-        = mxGetField(M_mx, 0, "objective_g1_sparse_colval");
-    if (!(objective_g1_sparse_colval_mx && mxIsInt32(objective_g1_sparse_colval_mx)))
-      mexErrMsgTxt("M_.objective_g1_sparse_colval should be an int32 array");
-
-    const mxArray* objective_g1_sparse_colptr_mx
-        = mxGetField(M_mx, 0, "objective_g1_sparse_colptr");
-    if (!(objective_g1_sparse_colptr_mx && mxIsInt32(objective_g1_sparse_colptr_mx)))
-      mexErrMsgTxt("M_.objective_g1_sparse_colptr should be an int32 array");
-
-    std::vector<const mxArray*> objective_gN_sparse_indices;
-    for (int o {2}; o <= kOrder; o++)
-      {
-        using namespace std::string_literals;
-        auto fieldname {"objective_g"s + std::to_string(o) + "_sparse_indices"};
-        const mxArray* indices = mxGetField(M_mx, 0, fieldname.c_str());
-        if (!(indices && mxIsInt32(indices)))
-          mexErrMsgTxt(("M_."s + fieldname + " should be an int32 array").c_str());
-        objective_gN_sparse_indices.push_back(indices);
-      }
-
     const int nSteps
         = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
 
@@ -309,14 +291,21 @@ extern "C"
     // run stochastic steady
     app.walkStochSteady();
 
+    const mxArray* objective_tmp_nbr_mx = mxGetField(M_mx, 0, "objective_tmp_nbr");
+    if (!(objective_tmp_nbr_mx && mxIsDouble(objective_tmp_nbr_mx)
+          && !mxIsComplex(objective_tmp_nbr_mx) && !mxIsSparse(objective_tmp_nbr_mx)
+          && mxGetNumberOfElements(objective_tmp_nbr_mx) >= static_cast<size_t>(kOrder + 1)))
+      mexErrMsgTxt("M_.objective_tmp_nbr should be a real dense array with strictly more elements "
+                   "than the order of derivation");
+    int ntt_objective = std::accumulate(mxGetPr(objective_tmp_nbr_mx),
+                                        mxGetPr(objective_tmp_nbr_mx) + kOrder + 1, 0);
+
     // Getting derivatives of the planner's objective function
-    std::unique_ptr<ObjectiveMFile> objectiveFile;
+    std::unique_ptr<ObjectiveAC> objectiveFile;
-    objectiveFile = std::make_unique<ObjectiveMFile>(
+    objectiveFile = std::make_unique<ObjectiveMFile>(fname, ntt_objective);
-        fname, kOrder, objective_g1_sparse_rowval_mx, objective_g1_sparse_colval_mx,
-        objective_g1_sparse_colptr_mx, objective_gN_sparse_indices);
 
     // make KordwDynare object
-    KordwDynare welfare(dynare, journal, modParams, std::move(objectiveFile), dr_order);
+    KordwDynare welfare(dynare, NNZD_obj, journal, modParams, std::move(objectiveFile), dr_order);
 
     // construct main K-order approximation class of welfare
     ApproximationWelfare appwel(welfare, discount_factor, app.get_rule_ders(),

mex/sources/k_order_welfare/objective_abstract_class.hh

@@ -0,0 +1,38 @@
+/*
+ * 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/>.
+ */
+
+#ifndef OBJECTIVE_ABSTRACT_CLASS_HH
+#define OBJECTIVE_ABSTRACT_CLASS_HH
+
+#include <vector>
+
+#include "twod_matrix.hh"
+
+class ObjectiveAC
+{
+protected:
+  int ntt; // Size of vector of temporary terms
+public:
+  ObjectiveAC(int ntt_arg) : ntt {ntt_arg} {};
+  virtual ~ObjectiveAC() = default;
+  virtual void eval(const Vector& y, const Vector& x, const Vector& params, Vector& residual,
+                    std::vector<TwoDMatrix>& md)
+      = 0;
+};
+#endif

mex/sources/k_order_welfare/objective_m.cc

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021-2024 Dynare Team
+ * Copyright © 2021-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -26,41 +26,88 @@
 
 #include "objective_m.hh"
 
-ObjectiveMFile::ObjectiveMFile(const std::string& modName, int kOrder_arg,
+ObjectiveMFile::ObjectiveMFile(const std::string& modName, int ntt_arg) :
-                               const mxArray* objective_g1_sparse_rowval_mx_arg,
+    ObjectiveAC(ntt_arg), ObjectiveMFilename {modName + ".objective.static"}
-                               const mxArray* objective_g1_sparse_colval_mx_arg,
-                               const mxArray* objective_g1_sparse_colptr_mx_arg,
-                               const std::vector<const mxArray*> objective_gN_sparse_indices_arg) :
-    ObjectiveMFilename {modName + ".objective.sparse.static"},
-    kOrder {kOrder_arg},
-    objective_g1_sparse_rowval_mx {objective_g1_sparse_rowval_mx_arg},
-    objective_g1_sparse_colval_mx {objective_g1_sparse_colval_mx_arg},
-    objective_g1_sparse_colptr_mx {objective_g1_sparse_colptr_mx_arg},
-    objective_gN_sparse_indices {objective_gN_sparse_indices_arg}
 {
 }
 
 void
-ObjectiveMFile::eval(const Vector& y, const Vector& x, const Vector& modParams, Vector& residual,
+ObjectiveMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm)
-                     const std::vector<int>& dynToDynpp,
-                     TensorContainer<FSSparseTensor>& derivatives) const
 {
-  mxArray* y_mx = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
+  int totalCols = mxGetN(sparseMat);
-  std::copy_n(y.base(), y.length(), mxGetPr(y_mx));
+  mwIndex* rowIdxVector = mxGetIr(sparseMat);
+  mwIndex* colIdxVector = mxGetJc(sparseMat);
 
-  mxArray* x_mx = mxCreateDoubleMatrix(1, x.length(), mxREAL);
+  assert(tdm.ncols() == 3);
-  std::copy_n(x.base(), x.length(), mxGetPr(x_mx));
+  /* Under MATLAB, the following check always holds at equality; under Octave,
+     there may be an inequality, because Octave diminishes nzmax if one gives
+     zeros in the values vector when calling sparse(). */
+  assert(tdm.nrows() >= static_cast<int>(mxGetNzmax(sparseMat)));
 
-  mxArray* params_mx = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL);
+  double* ptr = mxGetPr(sparseMat);
-  std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_mx));
 
-  mxArray *T_order_mx, *T_mx;
+  int rind = 0;
+  int output_row = 0;
+
+  for (int i = 0; i < totalCols; i++)
+    for (int j = 0; j < static_cast<int>((colIdxVector[i + 1] - colIdxVector[i])); j++, rind++)
+      {
+        tdm.get(output_row, 0) = rowIdxVector[rind] + 1;
+        tdm.get(output_row, 1) = i + 1;
+        tdm.get(output_row, 2) = ptr[rind];
+        output_row++;
+      }
+
+  /* If there are less elements than expected (that might happen if some
+     derivative is symbolically not zero but numerically zero at the evaluation
+     point), then fill in the matrix with empty entries, that will be
+     recognized as such by KordpDynare::populateDerivativesContainer() */
+  while (output_row < tdm.nrows())
+    {
+      tdm.get(output_row, 0) = 0;
+      tdm.get(output_row, 1) = 0;
+      tdm.get(output_row, 2) = 0;
+      output_row++;
+    }
+}
+
+void
+ObjectiveMFile::eval(const Vector& y, const Vector& x, const Vector& modParams, Vector& residual,
+                     std::vector<TwoDMatrix>& md) noexcept(false)
+{
+  mxArray* T_m = mxCreateDoubleMatrix(ntt, 1, mxREAL);
+
+  mxArray* y_m = mxCreateDoubleMatrix(y.length(), 1, mxREAL);
+  std::copy_n(y.base(), y.length(), mxGetPr(y_m));
+
+  mxArray* x_m = mxCreateDoubleMatrix(1, x.length(), mxREAL);
+  std::copy_n(x.base(), x.length(), mxGetPr(x_m));
+
+  mxArray* params_m = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL);
+  std::copy_n(modParams.base(), modParams.length(), mxGetPr(params_m));
+
+  mxArray* T_flag_m = mxCreateLogicalScalar(false);
+
+  {
+    // Compute temporary terms (for all orders)
+    std::string funcname = ObjectiveMFilename + "_g" + std::to_string(md.size()) + "_tt";
+    std::array<mxArray*, 1> plhs;
+    std::array prhs {T_m, y_m, x_m, params_m};
+
+    int retVal
+        = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
+    if (retVal != 0)
+      throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
+
+    mxDestroyArray(T_m);
+    T_m = plhs[0];
+  }
 
   {
     // Compute residuals
     std::string funcname = ObjectiveMFilename + "_resid";
-    std::array<mxArray*, 3> plhs;
+    std::array<mxArray*, 1> plhs;
-    std::array prhs {y_mx, x_mx, params_mx};
+    std::array prhs {T_m, y_m, x_m, params_m, T_flag_m};
 
     int retVal
         = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
@@ -69,103 +116,35 @@ ObjectiveMFile::eval(const Vector& y, const Vector& x, const Vector& modParams,
 
     residual = Vector {plhs[0]};
     mxDestroyArray(plhs[0]);
-
-    T_order_mx = plhs[1];
-    T_mx = plhs[2];
   }
 
-  {
+  for (size_t i = 1; i <= md.size(); i++)
-    // Compute Jacobian
-    std::string funcname = ObjectiveMFilename + "_g1";
-    std::array<mxArray*, 3> plhs;
-    std::array prhs {y_mx,
-                     x_mx,
-                     params_mx,
-                     const_cast<mxArray*>(objective_g1_sparse_rowval_mx),
-                     const_cast<mxArray*>(objective_g1_sparse_colval_mx),
-                     const_cast<mxArray*>(objective_g1_sparse_colptr_mx),
-                     T_order_mx,
-                     T_mx};
-
-    int retVal
-        = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
-    if (retVal != 0)
-      throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
-
-    assert(static_cast<int>(mxGetN(plhs[0])) == y.length());
-    double* g1_v {mxGetPr(plhs[0])};
-    mwIndex* g1_ir {mxGetIr(plhs[0])};
-    mwIndex* g1_jc {mxGetJc(plhs[0])};
-
-    IntSequence s(1, 0);
-    auto tensor = std::make_unique<FSSparseTensor>(1, y.length(), 1);
-
-    for (int j {0}; j < y.length(); j++)
-      for (mwIndex k {g1_jc[j]}; k < g1_jc[j + 1]; k++)
-        {
-          s[0] = dynToDynpp[j];
-          tensor->insert(s, g1_ir[k], g1_v[k]);
-        }
-
-    mxDestroyArray(plhs[0]);
-
-    mxDestroyArray(T_order_mx);
-    T_order_mx = plhs[1];
-
-    mxDestroyArray(T_mx);
-    T_mx = plhs[2];
-
-    derivatives.insert(std::move(tensor));
-  }
-
-  for (int o {2}; o <= kOrder; o++)
     {
-      // Compute higher derivatives
+      // Compute model derivatives
-      std::string funcname = ObjectiveMFilename + "_g" + std::to_string(o);
+      std::string funcname = ObjectiveMFilename + "_g" + std::to_string(i);
-      std::array<mxArray*, 3> plhs;
+      std::array<mxArray*, 1> plhs;
-      std::array prhs {y_mx, x_mx, params_mx, T_order_mx, T_mx};
+      std::array prhs {T_m, y_m, x_m, params_m, T_flag_m};
 
       int retVal
           = mexCallMATLAB(plhs.size(), plhs.data(), prhs.size(), prhs.data(), funcname.c_str());
       if (retVal != 0)
         throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
-      const mxArray* sparse_indices_mx {objective_gN_sparse_indices[o - 2]};
+      if (i == 1)
-      size_t nnz {mxGetM(sparse_indices_mx)};
-#if MX_HAS_INTERLEAVED_COMPLEX
-      const int32_T* sparse_indices {mxGetInt32s(sparse_indices_mx)};
-#else
-      const int32_T* sparse_indices {static_cast<const int32_T*>(mxGetData(sparse_indices_mx))};
-#endif
-
-      assert(mxGetNumberOfElements(plhs[0]) == nnz);
-      double* gN_v {mxGetPr(plhs[0])};
-
-      IntSequence s(o, 0);
-      auto tensor = std::make_unique<FSSparseTensor>(o, y.length(), 1);
-
-      for (size_t k {0}; k < nnz; k++)
         {
-          for (int i {0}; i < o; i++)
+          assert(static_cast<int>(mxGetM(plhs[0])) == md[i - 1].nrows());
-            s[i] = dynToDynpp[sparse_indices[k + (i + 1) * nnz] - 1];
+          assert(static_cast<int>(mxGetN(plhs[0])) == md[i - 1].ncols());
-          assert(s.isSorted());
+          std::copy_n(mxGetPr(plhs[0]), mxGetM(plhs[0]) * mxGetN(plhs[0]), md[i - 1].base());
-          tensor->insert(s, sparse_indices[k] - 1, gN_v[k]);
         }
+      else
+        unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[0], md[i - 1]);
 
       mxDestroyArray(plhs[0]);
-
-      mxDestroyArray(T_order_mx);
-      T_order_mx = plhs[1];
-
-      mxDestroyArray(T_mx);
-      T_mx = plhs[2];
-
-      derivatives.insert(std::move(tensor));
     }
 
-  mxDestroyArray(y_mx);
+  mxDestroyArray(T_m);
-  mxDestroyArray(x_mx);
+  mxDestroyArray(y_m);
-  mxDestroyArray(params_mx);
+  mxDestroyArray(x_m);
-  mxDestroyArray(T_order_mx);
+  mxDestroyArray(params_m);
-  mxDestroyArray(T_mx);
+  mxDestroyArray(T_flag_m);
 }

mex/sources/k_order_welfare/objective_m.hh

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021-2024 Dynare Team
+ * Copyright © 2021-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,34 +20,24 @@
 #ifndef OBJECTIVE_M_HH
 #define OBJECTIVE_M_HH
 
-#include <vector>
+#include "objective_abstract_class.hh"
 
-#include "dynmex.h"
+#include "mex.h"
+#include <dynmex.h>
 
-#include "Vector.hh"
+/**
-#include "sparse_tensor.hh"
+ * handles calls to <model>/+objective/static.m
-#include "t_container.hh"
+ *
-
+ **/
-// Handles calls to <model>/+objective/+sparse/static*.m
+class ObjectiveMFile : public ObjectiveAC
-class ObjectiveMFile
 {
 private:
   const std::string ObjectiveMFilename;
-  const int kOrder;
+  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray* sparseMat, TwoDMatrix& tdm);
-  const mxArray *const objective_g1_sparse_rowval_mx, *const objective_g1_sparse_colval_mx,
-                                                          *const objective_g1_sparse_colptr_mx;
-  // Stores M_.objective_gN_sparse_indices, starting from N=2
-  const std::vector<const mxArray*> objective_gN_sparse_indices;
 
 public:
-  ObjectiveMFile(const std::string& modName, int kOrder_arg,
+  explicit ObjectiveMFile(const std::string& modName, int ntt_arg);
-                 const mxArray* objective_g1_sparse_rowval_mx_arg,
-                 const mxArray* objective_g1_sparse_colval_mx_arg,
-                 const mxArray* objective_g1_sparse_colptr_mx_arg,
-                 const std::vector<const mxArray*> objective_gN_sparse_indices_arg);
   void eval(const Vector& y, const Vector& x, const Vector& params, Vector& residual,
-            const std::vector<int>& dynToDynpp, TensorContainer<FSSparseTensor>& derivatives) const
+            std::vector<TwoDMatrix>& md) override;
-      noexcept(false);
 };
-
 #endif

mex/sources/libkorder/sylv/QuasiTriangular.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -240,11 +240,16 @@ public:
     real = r;
   }
   virtual ~_matrix_iter() = default;
-  [[nodiscard]] bool
+  bool
   operator==(const _Self& it) const
   {
     return ptr == it.ptr;
   }
+  bool
+  operator!=(const _Self& it) const
+  {
+    return ptr != it.ptr;
+  }
   _TRef
   operator*() const
   {

mex/sources/libkorder/sylv/Vector.cc

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -114,6 +114,42 @@ 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()
 {
@@ -287,10 +323,17 @@ ConstVector::operator==(const ConstVector& y) const
   return i == len;
 }
 
-std::partial_ordering
+bool
-ConstVector::operator<=>(const ConstVector& y) const
+ConstVector::operator<(const ConstVector& y) const
 {
-  return std::lexicographical_compare_three_way(data, data + len, y.data, y.data + y.len);
+  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;
 }
 
 double

mex/sources/libkorder/sylv/Vector.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -25,7 +25,6 @@
    to avoid running virtual method invokation mechanism. Some
    members, and methods are thus duplicated */
 
-#include <compare>
 #include <complex>
 #include <utility>
 
@@ -109,6 +108,15 @@ 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)
@@ -211,9 +219,29 @@ public:
     return s;
   }
   // Exact equality
-  [[nodiscard]] bool operator==(const ConstVector& y) const;
+  bool operator==(const ConstVector& y) const;
+  bool
+  operator!=(const ConstVector& y) const
+  {
+    return !operator==(y);
+  }
   // Lexicographic ordering
-  [[nodiscard]] std::partial_ordering operator<=>(const ConstVector& y) const;
+  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]] double getNorm() const;
   [[nodiscard]] double getMax() const;

mex/sources/libkorder/tl/equivalence.cc

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -35,14 +35,12 @@ 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. */
 
-std::partial_ordering
+bool
-OrdSequence::operator<=>(const OrdSequence& s) const
+OrdSequence::operator<(const OrdSequence& s) const
 {
   double ta = average();
   double sa = s.average();
-  if (auto cmp1 = ta <=> sa; cmp1 != 0)
+  return (ta < sa || ((ta == sa) && (operator[](0) > s[0])));
-    return cmp1;
-  return operator[](0) <=> s[0];
 }
 
 bool

mex/sources/libkorder/tl/equivalence.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -55,7 +55,6 @@
 
 #include "int_sequence.hh"
 
-#include <compare>
 #include <list>
 #include <string>
 #include <vector>
@@ -63,7 +62,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. */
 
@@ -75,9 +74,9 @@ public:
   OrdSequence() : data()
   {
   }
-  [[nodiscard]] bool operator==(const OrdSequence& s) const;
+  bool operator==(const OrdSequence& s) const;
   int operator[](int i) const;
-  [[nodiscard]] std::partial_ordering operator<=>(const OrdSequence& s) const;
+  bool operator<(const OrdSequence& s) const;
   [[nodiscard]] const std::vector<int>&
   getData() const
   {
@@ -121,7 +120,12 @@ public:
   // Copy constructor plus gluing i1 and i2 in one class
   Equivalence(const Equivalence& e, int i1, int i2);
 
-  [[nodiscard]] bool operator==(const Equivalence& e) const;
+  bool operator==(const Equivalence& e) const;
+  bool
+  operator!=(const Equivalence& e) const
+  {
+    return !operator==(e);
+  }
   [[nodiscard]] int
   getN() const
   {

mex/sources/libkorder/tl/gs_tensor.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -76,11 +76,16 @@ public:
   // Constructs the tensor dimensions for slicing (see the implementation for details)
   TensorDimens(const IntSequence& ss, const IntSequence& coor);
 
-  [[nodiscard]] bool
+  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

mex/sources/libkorder/tl/int_sequence.cc

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 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);
 }
 
-std::strong_ordering
+bool
-IntSequence::operator<=>(const IntSequence& s) const
+IntSequence::operator<(const IntSequence& s) const
 {
-  return std::lexicographical_compare_three_way(data, data + length, s.data, s.data + s.length);
+  return std::lexicographical_compare(data, data + length, s.data, s.data + s.length);
 }
 
 bool

mex/sources/libkorder/tl/int_sequence.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -43,7 +43,6 @@
 #define INT_SEQUENCE_HH
 
 #include <algorithm>
-#include <compare>
 #include <initializer_list>
 #include <utility>
 #include <vector>
@@ -120,7 +119,12 @@ public:
     if (destroy)
       delete[] data;
   }
-  [[nodiscard]] bool operator==(const IntSequence& s) const;
+  bool operator==(const IntSequence& s) const;
+  bool
+  operator!=(const IntSequence& s) const
+  {
+    return !operator==(s);
+  }
   int&
   operator[](int i)
   {
@@ -137,10 +141,15 @@ 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. */
-  [[nodiscard]] std::strong_ordering operator<=>(const IntSequence& s) const;
+  bool operator<(const IntSequence& s) const;
+  bool
+  operator<=(const IntSequence& s) const
+  {
+    return (operator==(s) || operator<(s));
+  }
   [[nodiscard]] bool lessEq(const IntSequence& s) const;
   [[nodiscard]] bool less(const IntSequence& s) const;
 

mex/sources/libkorder/tl/kron_prod.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -86,7 +86,7 @@ public:
 
   KronProdDimens& operator=(const KronProdDimens& kd) = default;
   KronProdDimens& operator=(KronProdDimens&& kd) = default;
-  [[nodiscard]] bool
+  bool
   operator==(const KronProdDimens& kd) const
   {
     return rows == kd.rows && cols == kd.cols;

mex/sources/libkorder/tl/permutation.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 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))
   {
   }
-  [[nodiscard]] bool
+  bool
   operator==(const Permutation& p) const
   {
     return permap == p.permap;

mex/sources/libkorder/tl/ps_tensor.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -115,7 +115,7 @@ public:
   {
     per.apply(nvmax);
   }
-  [[nodiscard]] bool
+  bool
   operator==(const PerTensorDimens& td) const
   {
     return TensorDimens::operator==(td) && per == td.per;

mex/sources/libkorder/tl/symmetry.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -126,11 +126,16 @@ public:
   {
     return run;
   }
-  [[nodiscard]] bool
+  bool
-  operator==(const symiterator& it) const
+  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

mex/sources/libkorder/tl/tensor.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2024 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -150,11 +150,16 @@ public:
     {
       return offset;
     }
-    [[nodiscard]] bool
+    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
     {

mex/sources/riccati_update/mexFunction.f08

@@ -1,98 +0,0 @@
-! 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/>.
-
-! Implements Ptmp = T*(P-K*Z*P)*transpose(T)+Q where 
-! P is the (r x r) variance-covariance matrix of the state vector
-! T is the (r x r) transition matrix of the state vector
-! K is the (r x n) gain matrix
-! Z is the (n x r) matrix linking observable variables to state variables
-! Q is the (r x r) variance-covariance matrix of innovations in the state equation
-! and accounting for different properties:
-! P is a (symmetric) positive semi-definite matrix
-! T can be triangular
-
-subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
-   use matlab_mex
-   use blas
-   implicit none (type, external)
-
-   type(c_ptr), dimension(*), intent(in), target :: prhs
-   type(c_ptr), dimension(*), intent(out) :: plhs
-   integer(c_int), intent(in), value :: nlhs, nrhs
-
-   real(real64), dimension(:,:), pointer, contiguous :: P, T, K, Z, Q, Pnew
-   real(real64), dimension(:,:), allocatable :: tmp1, tmp2
-   integer :: i, n, r
-   character(kind=c_char, len=2) :: num2str 
-
-   ! 0. Checking the consistency and validity of input arguments
-   if (nrhs /= 5_c_int) then
-      call mexErrMsgTxt("Must have 5 input arguments")
-   end if
-   if (nlhs > 1_c_int) then
-      call mexErrMsgTxt("Too many output arguments")
-   end if
-
-   do i=1,5
-      if (.not. (c_associated(prhs(i)) .and. mxIsDouble(prhs(i)) .and. & 
-          (.not. mxIsComplex(prhs(i))) .and. (.not. mxIsSparse(prhs(i))))) then
-            write (num2str,"(i2)") i
-            call mexErrMsgTxt("Argument " // trim(num2str) // " should be a real dense matrix")
-      end if
-   end do
-
-   r = int(mxGetM(prhs(1)))            ! Number of states
-   n = int(mxGetN(prhs(3)))            ! Number of observables
-
-   if ((r /= mxGetN(prhs(1))) &        ! Number of columns of P
-      &.or. (r /= mxGetM(prhs(2))) &   ! Number of lines of T
-      &.or. (r /= mxGetN(prhs(2))) &   ! Number of columns of T
-      &.or. (r /= mxGetM(prhs(3))) &   ! Number of lines of K
-      &.or. (n /= mxGetM(prhs(4))) &   ! Number of lines of Z
-      &.or. (r /= mxGetN(prhs(4))) &   ! Number of columns of Z
-      &.or. (r /= mxGetM(prhs(5))) &   ! Number of lines of Q
-      &.or. (r /= mxGetN(prhs(5))) &   ! Number of columns of Q
-      ) then  
-      call mexErrMsgTxt("Input dimension mismatch")
-   end if
-
-   ! 1. Storing the relevant information in Fortran format
-   P(1:r,1:r) => mxGetPr(prhs(1))
-   T(1:r,1:r) => mxGetPr(prhs(2))
-   K(1:r,1:n) => mxGetPr(prhs(3))
-   Z(1:n,1:r) => mxGetPr(prhs(4))
-   Q(1:r,1:r) => mxGetPr(prhs(5))
-
-   plhs(1) = mxCreateDoubleMatrix(int(r, mwSize), int(r, mwSize), mxREAL)
-   Pnew(1:r, 1:r) => mxGetPr(plhs(1))
-
-   ! 2. Computing the Riccati update of the P matrix 
-   allocate(tmp1(r,r), tmp2(r,r))
-   ! Pnew <- Q
-   Pnew = Q
-   ! tmp1 <- K*Z
-   call matmul_add("N", "N", 1._real64, K, Z, 0._real64, tmp1)
-   ! tmp2 <- P
-   tmp2 = P
-   ! tmp2 <- tmp2 - tmp1*P
-   call matmul_add("N", "N", -1._real64, tmp1, P, 1._real64, tmp2)
-   ! tmp1 <- T*tmp2
-   call matmul_add("N", "N", 1._real64, T, tmp2, 0._real64, tmp1)
-   ! Pnew <- tmp1*T' + Pnew
-   call matmul_add("N", "T", 1._real64, tmp1, T, 1._real64, Pnew)
-
-end subroutine mexFunction

mex/sources/sobol/gaussian.hh

@@ -2,7 +2,7 @@
 **
 ** Pseudo code of the algorithm is given at http://home.online.no/~pjacklam/notes/invnorm
 **
-** Copyright © 2010-2024 Dynare Team
+** Copyright © 2010-2023 Dynare Team
 **
 ** This file is part of Dynare.
 **
@@ -104,8 +104,9 @@ icdf(const T uniform)
   return gaussian;
 }
 
+template<typename T>
 void
-icdfm(int n, auto* U)
+icdfm(int n, T* U)
 {
 #pragma omp parallel for
   for (int i = 0; i < n; i++)
@@ -113,8 +114,9 @@ icdfm(int n, auto* U)
   return;
 }
 
+template<typename T>
 void
-icdfmSigma(int d, int n, auto* U, const double* LowerCholSigma)
+icdfmSigma(int d, int n, T* U, const double* LowerCholSigma)
 {
   double one = 1.0;
   double zero = 0.0;
@@ -126,8 +128,9 @@ icdfmSigma(int d, int n, auto* U, const double* LowerCholSigma)
   copy_n(tmp.begin(), d * n, U);
 }
 
+template<typename T>
 void
-usphere(int d, int n, auto* U)
+usphere(int d, int n, T* U)
 {
   icdfm(n * d, U);
 #pragma omp parallel for
@@ -144,8 +147,9 @@ usphere(int d, int n, auto* U)
     }
 }
 
+template<typename T>
 void
-usphereRadius(int d, int n, double radius, auto* U)
+usphereRadius(int d, int n, double radius, T* U)
 {
   icdfm(n * d, U);
 #pragma omp parallel for

preprocessor

@@ -1 +1 @@
-Subproject commit d8866fbec8a09df0d8fa1d2e7e60b2aae8685ea8
+Subproject commit 20acdbc119648f8ae6438e7f3a18943b4c23e7d2

tests/riccatiupdate.m

@@ -1,87 +0,0 @@
-source_dir = getenv('source_root');
-addpath([source_dir filesep 'matlab']);
-
-dynare_config;
-
-testFailed = 0;
-
-skipline()
-disp('***  TESTING: riccatiupdate.m ***');
-
-t0 = clock;
-
-% Set the number of experiments for time measurement
-N = 5000;
-% Set the dimension of the problem to be solved.
-r = 50;
-n = 100;
-tol = 1e-15;
-% Set the input arguments
-% P, Q: use the fact that for any real matrix A, A'*A is positive semidefinite
-P = rand(n,r);
-P = P'*P;
-Q = rand(n,r);
-Q = Q'*Q;
-K = rand(r,n);
-Z = rand(n,r);
-T = rand(r,r);
-% Computing an upperbound for the norm the updated variance-covariance matrix
-ub = norm(T,1)^2*norm(P,1)*(1+norm(K*Z,1))+norm(Q,1);
-% Weighting the P and Q matrices to keep the norm of the variance-covariance matrix below 1
-P = 0.5*P/ub;
-Q = 0.5*Q/ub;
-
-% 1. Update the state vairance-covariance matrix with Matlab
-tElapsed1 = 0.;
-tic;
-for i=1:N
-   Ptmp_matlab = T*(P-K*Z*P)*transpose(T)+Q;
-end
-tElapsed1 = toc;
-disp(['Elapsed time for the Matlab Riccati update is: ' num2str(tElapsed1) ' (N=' int2str(N) ').'])
-
-% 2. Update the state varance-covariance matrix with the mex routine
-tElapsed2 = 0.;
-Ptmp_fortran = P;
-try
-   tic;
-   for i=1:N
-      Ptmp_fortran = riccati_update(P, T, K, Z, Q);
-   end
-   tElapsed2 = toc;
-   disp(['Elapsed time for the Fortran Riccati update is: ' num2str(tElapsed2) ' (N=' int2str(N) ').'])
-   R = norm(Ptmp_fortran-Ptmp_matlab,1);
-   if (R > tol)
-      testFailed = testFailed+1;
-      dprintf('The Fortran Riccati update is wrong')
-   end
-catch
-   testFailed = testFailed+1;
-   dprintf('Fortran Riccati update failed')
-end
-
-% Compare the Fortran and Matlab execution time
-if tElapsed1<tElapsed2
-    skipline()
-    dprintf('Matlab Riccati update is %5.2f times faster than its Fortran counterpart.', tElapsed2/tElapsed1)
-    skipline()
-else
-    skipline()
-    dprintf('Fortran Riccati update is %5.2f times faster than its Matlab counterpart.', tElapsed1/tElapsed2)
-    skipline()
-end
-
-% Compare results after multiple calls
-N = 50;
-disp(['After 1 update using the Riccati formula, the norm-1 discrepancy is ' num2str(norm(Ptmp_fortran-Ptmp_matlab,1)) '.']);
-for i=2:N
-   Ptmp_matlab = T*(Ptmp_matlab-K*Z*Ptmp_matlab)*transpose(T)+Q;
-   Ptmp_fortran = riccati_update(Ptmp_fortran, T, K, Z, Q);
-   disp(['After ' int2str(i) ' updates using the Riccati formula, the norm-1 discrepancy is ' num2str(norm(Ptmp_fortran-Ptmp_matlab,1)) '.'])
-end
-
-t1 = clock;
-
-fprintf('\n*** Elapsed time (in seconds): %.1f\n\n', etime(t1, t0));
-
-quit(testFailed > 0)

windows/deps/Makefile

@@ -1,4 +1,4 @@
-# Copyright © 2017-2024 Dynare Team
+# Copyright © 2017-2023 Dynare Team
 #
 # This file is part of Dynare.
 #
@@ -49,7 +49,7 @@ clean-all: clean-lib clean-src clean-tar
 
 tarballs/slicot-$(SLICOT_VERSION).tar.gz:
 	mkdir -p tarballs
-	wget $(WGET_OPTIONS) -O $@ https://deb.debian.org/debian/pool/main/s/slicot/slicot_$(SLICOT_VERSION).orig.tar.xz
+	wget $(WGET_OPTIONS) -O $@ https://deb.debian.org/debian/pool/main/s/slicot/slicot_$(SLICOT_VERSION).orig.tar.gz
 
 sources64/slicot-$(SLICOT_VERSION)-with-64bit-integer: tarballs/slicot-$(SLICOT_VERSION).tar.gz
 	rm -rf sources64/slicot-*-with-64bit-integer
@@ -64,13 +64,13 @@ sources64/slicot-$(SLICOT_VERSION)-with-32bit-integer-and-underscore: tarballs/s
 	touch $@
 
 lib64/Slicot/without-underscore/libslicot64_pic.a: sources64/slicot-$(SLICOT_VERSION)-with-64bit-integer
-	make -C $< -f makefile_Unix lib SLICOTLIB=../libslicot64_pic.a OPTS="$(FFLAGS) -fno-underscoring -fdefault-integer-8" FORTRAN=x86_64-w64-mingw32-gfortran LOADER=x86_64-w64-mingw32-gfortran ARCH=x86_64-w64-mingw32-ar
+	make -C $< lib SLICOTLIB=../libslicot64_pic.a OPTS="$(FFLAGS) -fno-underscoring -fdefault-integer-8" FORTRAN=x86_64-w64-mingw32-gfortran LOADER=x86_64-w64-mingw32-gfortran ARCH=x86_64-w64-mingw32-ar
 	x86_64-w64-mingw32-strip --strip-debug $</libslicot64_pic.a
 	mkdir -p $(dir $@)
 	cp $</libslicot64_pic.a $@
 
 lib64/Slicot/with-underscore/libslicot_pic.a: sources64/slicot-$(SLICOT_VERSION)-with-32bit-integer-and-underscore
-	make -C $< -f makefile_Unix lib SLICOTLIB=../libslicot_pic.a OPTS="$(FFLAGS)" FORTRAN=x86_64-w64-mingw32-gfortran LOADER=x86_64-w64-mingw32-gfortran ARCH=x86_64-w64-mingw32-ar
+	make -C $< lib SLICOTLIB=../libslicot_pic.a OPTS="$(FFLAGS)" FORTRAN=x86_64-w64-mingw32-gfortran LOADER=x86_64-w64-mingw32-gfortran ARCH=x86_64-w64-mingw32-ar
 	x86_64-w64-mingw32-strip --strip-debug $</libslicot_pic.a
 	mkdir -p $(dir $@)
 	cp $</libslicot_pic.a $@
@@ -158,7 +158,7 @@ mingw64: tarballs/mingw-w64-x86_64-gcc-$(MINGW64_GCC_VERSION)-any.pkg.tar.zst ta
 
 tarballs/mingw-w64-x86_64-%-any.pkg.tar.zst:
 	mkdir -p tarballs
-	wget $(WGET_OPTIONS) -O $@ http://repo.msys2.org/mingw/x86_64/$(notdir $@)
+	wget $(WGET_OPTIONS) -O $@ https://www.dynare.org/windows-pkg-build/msys2/6.x/$(notdir $@)
 
 clean-msys2:
 	rm -rf lib64-msys2

windows/deps/versions.mk

@@ -1,4 +1,4 @@
-SLICOT_VERSION = 5.9~20240205.gita037f7e
+SLICOT_VERSION = 5.0+20101122
 X13AS_VERSION = 1-1-b60
 
 OCTAVE_VERSION = 8.4.0