[WIP] Add Sequential Monte Carlo samplers.

Name of the parameter.

.clang-tidy

@@ -0,0 +1,9 @@
+# NB: to use clang-tidy on the MEX source code, make sure that you have
+# libomp-dev installed (the LLVM implementation of OpenMP)
+
+# TODO: add the following check families:
+# - performance-*
+# - bugprone-*
+# - cppcoreguidelines-
+
+Checks: 'modernize-*,-modernize-use-trailing-return-type,-clang-diagnostic-unqualified-std-cast-call'

README.md

@@ -479,7 +479,7 @@ If you want a certain version (e.g. 5.x) , then add `--single-branch --branch 5.
 ```sh
 export BUILDDIR=build-matlab
 export MATLABPATH=/Applications/MATLAB_R2023b.app
-arch -$ARCH meson setup --native-file scripts/homebrew-native-$ARCH.ini -Dmatlab_path=$MATLABPATH -Dbuildtype=debugoptimized -Dfortran_args="['-B','$DYNAREDIR/slicot/lib']" $BUILDDIR
+arch -$ARCH meson setup --native-file macOS/homebrew-native-$ARCH.ini -Dmatlab_path=$MATLABPATH -Dbuildtype=debugoptimized -Dfortran_args="['-B','$DYNAREDIR/slicot/lib']" $BUILDDIR
 ```
 where you need to adapt the path to MATLAB.
 Similarly, if you want to compile for Octave, replace the `-Dmatlab_path` option by `-Dbuild_for=octave`, and change the build directory to `build-octave`.
@@ -514,4 +514,4 @@ e.g. by adding this to your mod file. Alternatively, you can create a `startup.m
 ## Docker
 We offer a variety of pre-configured Docker containers for Dynare, pre-configured with Octave and MATLAB including all recommended toolboxes.
 These are readily available for your convenience on [Docker Hub](https://hub.docker.com/r/dynare/dynare).
-The docker folder contains [information and instructions](docker/README.md) to interact, built and customize the containers.
+The `scripts/docker` folder contains [information and instructions](scripts/docker/README.md) to interact, built and customize the containers.

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

@@ -15,11 +15,16 @@ related to the model (and hence not placed in the model file). At the
 moment, it is only used when using Dynare to run parallel
 computations.
 
-On Linux and macOS, the default location of the configuration file is
-``$HOME/.dynare``, while on Windows it is ``%APPDATA%\dynare.ini``
-(typically ``c:\Users\USERNAME\AppData\dynare.ini``). You
-can specify a non standard location using the ``conffile`` option of
-the ``dynare`` command (see :ref:`dyn-invoc`).
+On Linux and macOS, the configuration file is searched by default under
+``dynare/dynare.ini`` in the configuration directories defined by the XDG
+specification (typically ``$HOME/.config/dynare/dynare.ini`` for the
+user-specific configuration and ``/etc/xdg/dynare/dynare.ini`` for the
+system-wide configuration, the former having precedence over the latter). Under
+Windows, the configuration file is searched by default in
+``%APPDATA%\dynare\dynare.ini`` (typically
+``c:\Users\USERNAME\AppData\Roaming\dynare\dynare.ini``). You can specify a non
+standard location using the ``conffile`` option of the ``dynare`` command (see
+:ref:`dyn-invoc`).
 
 The parsing of the configuration file is case-sensitive and it should
 take the following form, with each option/choice pair placed on a
@@ -76,8 +81,15 @@ processing. Currently, there is only one option available.
 
     .. option:: GlobalInitFile = PATH_AND_FILE
 
-        The location of the global initialization file to be run at
-        the end of ``global_initialization.m``.
+        The location of a global initialization file that can be used to
+        customize some Dynare internals (typically default option values). This
+        is a MATLAB/Octave script.
+
+        If this option is not specified, Dynare will look for a
+        ``global_init.m`` file in its configuration directory (typically
+        ``$HOME/.config/dynare/global_init.m`` under Linux and macOS, and
+        ``c:\Users\USERNAME\AppData\Roaming\dynare\global_init.m`` under
+        Windows).
 
     *Example*
 

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

@@ -4787,31 +4787,6 @@ Computing the stochastic solution
        and this option has no effect. More references can be found
        `here <https://archives.dynare.org/DynareWiki/PartialInformation>`__ .
 
-    .. option:: sylvester = OPTION
-
-       Determines the algorithm used to solve the Sylvester equation
-       for block decomposed model. Possible values for OPTION are:
-
-           ``default``
-
-                Uses the default solver for Sylvester equations
-                (``gensylv``) based on Ondra Kamenik’s algorithm (see
-                `here
-                <https://www.dynare.org/assets/team-presentations/sylvester.pdf>`__
-                for more information).
-
-           ``fixed_point``
-
-                Uses a fixed point algorithm to solve the Sylvester
-                equation (``gensylv_fp``). This method is faster than
-                the default one for large scale models.
-
-       |br| Default value is ``default``.
-
-    .. option:: sylvester_fixed_point_tol = DOUBLE
-
-       The convergence criterion used in the fixed point
-       Sylvester solver. Its default value is ``1e-12``.
 
     .. option:: dr = OPTION
 
@@ -7935,15 +7910,6 @@ observed variables.
 
        See :opt:`aim_solver`.
 
-    .. option:: sylvester = OPTION
-
-       See :opt:`sylvester <sylvester = OPTION>`.
-
-    .. option:: sylvester_fixed_point_tol = DOUBLE
-
-       See :opt:`sylvester_fixed_point_tol <sylvester_fixed_point_tol
-       = DOUBLE>` .
-
     .. option:: lyapunov = OPTION
 
        Determines the algorithm used to solve the Lyapunov equation to
@@ -9487,16 +9453,6 @@ adding prior information comes at the cost of a loss in efficiency of the estima
         See :opt:`lyapunov_doubling_tol <lyapunov_doubling_tol = DOUBLE>`.
         Default: ``1e-16``.
 
-    .. option:: sylvester = OPTION
-
-        See :opt:`sylvester <sylvester = OPTION>`.
-        Default: ``default``, i.e. uses ``gensylv``.
-
-    .. option:: sylvester_fixed_point_tol = DOUBLE
-
-        See :opt:`sylvester_fixed_point_tol <sylvester_fixed_point_tol = DOUBLE>`.
-        Default: ``1e-12``.
-
     .. option:: qz_criterium = DOUBLE
 
         See :opt:`qz_criterium <qz_criterium = DOUBLE>`.

examples/dynare++/README.txt

@@ -1,6 +0,0 @@
-Assuming that the dynare++ binary is in your PATH, you can run the example by using the following command
-in a Command Prompt Window:
-
-... > dynare++ example1.mod
-
-Please, read the manual (doc\dynare++\dynare++-tutorial.pdf) for a description of the generated output.

examples/dynare++/example1.mod

@@ -1,45 +0,0 @@
-/*
- * This Dynare++ mod-file implements the RBC model with time-to-build
- * described in Kamenik (2011): "DSGE Models with Dynare++. A Tutorial."
- * Note that Dynare++ uses the same stock-at-the-end-of-period timing convention
- * as the regular Dynare
-*/
-
-var Y, C, K, A, H, B;
-
-varexo EPS, NU;
-
-parameters beta, rho, alpha, delta, theta, psi, tau;
-
-alpha = 0.36;
-rho = 0.95;
-tau = 0.025;
-beta = 1/(1.03^0.25);
-delta = 0.025;
-psi = 0;
-theta = 2.95;
-
-model;
-  C*theta*H^(1+psi) = (1-alpha)*Y;
-  beta*exp(B)*C/exp(B(1))/C(1)*
-  (exp(B(1))*alpha*Y(1)/K(1)+1-delta) = 1;
-  Y = exp(A)*K^alpha*H^(1-alpha);
-  K = exp(B(-1))*(Y(-1)-C(-1)) + (1-delta)*K(-1);
-  A = rho*A(-1) + tau*B(-1) + EPS;
-  B = tau*A(-1) + rho*B(-1) + NU;
-end;
-
-initval;
-  A = 0;
-  B = 0;
-  H = ((1-alpha)/(theta*(1-(delta*alpha)/(1/beta-1+delta))))^(1/(1+psi));
-  Y = (alpha/(1/beta-1+delta))^(alpha/(1-alpha))*H;
-  K = alpha/(1/beta-1+delta)*Y;
-  C = Y - delta*K;
-end;
-
-vcov = [0.0002 0.00005;
-0.00005 0.0001
-];
-
-order = 7;

macOS/build.sh

@@ -54,7 +54,7 @@ LIB64="$ROOTDIR"/macOS/deps/"$PKG_ARCH"/lib64
 ## - the macOS linker is different from GNU ld and does not have the equivalent of -Bstatic/-Bdynamic
 ## - libgfortran.spec does not include --as-needed on macOS, hence it will link the library anyways
 ## Also, it does not seem possible to override libgfortran.spec with the --specs option.
-GCC_VERSION=$(sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" "$ROOTDIR"/scripts/homebrew-native-"$PKG_ARCH".ini)
+GCC_VERSION=$(sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" "$ROOTDIR"/macOS/homebrew-native-"$PKG_ARCH".ini)
 
 QUADMATH_DIR=$(mktemp -d)
 ln -s "$BREWDIR"/opt/gcc/lib/gcc/"$GCC_VERSION"/libquadmath.a "$QUADMATH_DIR"
@@ -67,7 +67,7 @@ cd "$ROOTDIR"
 # NB: the addition of -Wl,-ld_classic is a workaround for https://github.com/mesonbuild/meson/issues/12282 (see also the native file)
 common_meson_opts=(-Dbuild_for=matlab -Dbuildtype=release -Dprefer_static=true -Dfortran_args="[ '-B', '$LIB64/Slicot/' ]" \
                    -Dc_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dcpp_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" -Dfortran_link_args="[ '-Wl,-ld_classic', '-L$QUADMATH_DIR' ]" \
-                   --native-file scripts/homebrew-native-$PKG_ARCH.ini)
+                   --native-file macOS/homebrew-native-$PKG_ARCH.ini)
 
 # Build for MATLAB ⩾ R2018b (x86_64) and MATLAB ⩾ R2023b (arm64)
 arch -"$PKG_ARCH" meson setup "${common_meson_opts[@]}" -Dmatlab_path="$MATLAB_PATH" build-matlab --wipe

macOS/deps/Makefile

@@ -22,7 +22,7 @@ DEPS_ARCH ?= x86_64 # use x86_64 by default
 
 BREWDIR := $(if $(filter arm64,$(DEPS_ARCH)),/opt/homebrew,/usr/local)
 
-GCC_VERSION = $(shell sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" ../../scripts/homebrew-native-$(DEPS_ARCH).ini)
+GCC_VERSION = $(shell sed -En "/^c[[:space:]]*=/s/c[[:space:]]*=[[:space:]]*'.*gcc-([0-9]+)'/\1/p" ../homebrew-native-$(DEPS_ARCH).ini)
 
 ROOT_PATH = $(realpath .)
 

matlab/+mom/default_option_mom_values.m

@@ -160,8 +160,6 @@ options_mom_ = set_default_option(options_mom_,'lyapunov_srs',false);
 options_mom_ = set_default_option(options_mom_,'lyapunov_complex_threshold',1e-15);     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
 options_mom_ = set_default_option(options_mom_,'lyapunov_fixed_point_tol',1e-10);       % convergence criterion used in the fixed point Lyapunov solver
 options_mom_ = set_default_option(options_mom_,'lyapunov_doubling_tol',1e-16);          % convergence criterion used in the doubling algorithm
-options_mom_ = set_default_option(options_mom_,'sylvester_fp',false);                   % determines whether to use fixed point algorihtm to solve Sylvester equation (gensylv_fp), faster for large scale models
-options_mom_ = set_default_option(options_mom_,'sylvester_fixed_point_tol',1e-12);      % convergence criterion used in the fixed point Sylvester solver
 
 % mode check plot
 options_mom_.mode_check.nolik = false;                                                          % we don't do likelihood (also this initializes mode_check substructure)

matlab/AHessian.m

@@ -1,152 +0,0 @@
-function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-%
-% computes the asymptotic hessian matrix of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% Thanks to  Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-
-% Copyright © 2011-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/>.
-
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-pp   = size(Y,1);                               % Maximum number of observed variables.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-lik  = zeros(smpl,1);                           % Initialization of the vector gathering the densities.
-LIK  = Inf;                                     % Default value of the log likelihood.
-if nargout > 1
-    DLIK  = zeros(k,1);                             % Initialization of the score.
-end
-AHess  = zeros(k,k);                             % Initialization of the Hessian
-Da    = zeros(mm,k);                             % State vector.
-Dv = zeros(length(mf),k);
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady && t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-        lik(t) = log(det(F))+transpose(v)*iF*v;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start && nargout > 1
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        vecDPmf = reshape(DP(mf,mf,:),[],k);
-        %                     iPmf = inv(P(mf,mf));
-        if t>=start
-            AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start && nargout >1
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        if t>=start
-            AHess = AHess + Dv'*iF*Dv;
-        end
-        a = T*(a+K*v);
-        lik(t) = transpose(v)*iF*v;
-    end
-    AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
-    if nargout > 1
-        for ii = 1:k
-            %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
-        end
-    end
-    lik(t0:smpl) = lik(t0:smpl) + log(det(F));
-    %         for ii = 1:k;
-    %             for jj = 1:ii
-    %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
-    %             end
-    %         end
-end
-
-AHess = -AHess;
-if nargout > 1
-    DLIK = DLIK/2;
-end
-% adding log-likelihhod constants
-lik = (lik + pp*log(2*pi))/2;
-
-LIK = sum(lik(start:end)); % Minus the log-likelihood.
-                           % end of main function
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman

matlab/DsgeSmoother.m

@@ -275,7 +275,7 @@ if kalman_algo == 1 || kalman_algo == 3
             kalman_algo = 2;
         elseif kalman_algo == 3
             fprintf('\nDsgeSmoother: Switching to univariate filter. This is usually due to co-integration in diffuse filter,\n')
-            fprintf(' otherwise it may be a sign of stochastic singularity.\n')
+            fprintf('otherwise it may be a sign of stochastic singularity.\n')
             kalman_algo = 4;
         else
             error('This case shouldn''t happen')

matlab/GetPosteriorParametersStatistics.m

@@ -79,7 +79,7 @@ else
 end
 
 if ishssmc(options_)
-    num_draws = options_.posterior_sampler_options.hssmc.nparticles;
+    num_draws = options_.posterior_sampler_options.hssmc.particles;
     hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
 else
     num_draws=NumberOfDraws*mh_nblck;

matlab/bseastr.m

@@ -1,47 +0,0 @@
-function x = bseastr(s1,s2)
-
-% Copyright © 2001-2017 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-m = size(s1,1) ;
-x = zeros(m,1) ;
-s1=upper(deblank(s1));
-s2=upper(deblank(s2));
-
-for im = 1:m
-    key = s1(im,:) ;
-    h = size(s2,1) ;
-    l = 1 ;
-    while l <= h
-        mid = round((h+l)/2) ;
-        temp = s2(mid,:) ;
-        if ~ strcmp(key,temp)
-            for i = 1:min(length(key),length(temp))
-                if temp(i) > key(i)
-                    h = mid - 1 ;
-                    break
-                else
-                    l = mid + 1 ;
-                    break
-                end
-            end
-        else
-            x(im) = mid ;
-            break
-        end
-    end
-end

matlab/check_posterior_sampler_options.m

@@ -45,7 +45,9 @@ end
 if init
     % set default options and user defined options
     posterior_sampler_options.posterior_sampling_method = options_.posterior_sampler_options.posterior_sampling_method;
-    posterior_sampler_options.bounds = bounds;
+    if ~issmc(options_)
+        posterior_sampler_options.bounds = bounds;
+    end
 
     switch posterior_sampler_options.posterior_sampling_method
 
@@ -338,87 +340,73 @@ if init
         end
 
         % moreover slice must be associated to:
-        %     options_.mh_posterior_mode_estimation = false;
-        % this is done below, but perhaps preprocessing should do this?
+            %     options_.mh_posterior_mode_estimation = false;
+            % this is done below, but perhaps preprocessing should do this?
 
-        if ~isempty(posterior_sampler_options.mode)
-            % multimodal case
-            posterior_sampler_options.rotated = 1;
-            posterior_sampler_options.WR=[];
-        end
-        %     posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
-
-
-        posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
-        if options_.load_mh_file
-            posterior_sampler_options.slice_initialize_with_mode = 0;
-        else
-            if ~posterior_sampler_options.slice_initialize_with_mode
-                posterior_sampler_options.invhess=[];
+            if ~isempty(posterior_sampler_options.mode)
+                % multimodal case
+                posterior_sampler_options.rotated = 1;
+                posterior_sampler_options.WR=[];
             end
-        end
+            %     posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
 
-        if ~isempty(posterior_sampler_options.mode_files) % multimodal case
-            modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
-            load(modes, 'xparams')
-            if size(xparams,2)<2
-                error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
-            end
-            for j=1:size(xparams,2)
-                mode(j).m=xparams(:,j);
-            end
-            posterior_sampler_options.mode = mode;
-            posterior_sampler_options.rotated = 1;
-            posterior_sampler_options.WR=[];
-        end
 
-      case 'hssmc'
-
-        posterior_sampler_options.parallel_bar_refresh_rate=1;
-        posterior_sampler_options.serial_bar_title='HS Sequential Monte-Carlo';
-
-        % default options
-        posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
-
-        % user defined options
-        if ~isempty(options_.posterior_sampler_options.sampling_opt)
-            options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
-            for i=1:rows(options_list)
-                switch options_list{i,1}
-                  case 'proposal_distribution'
-                    if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
-                         strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
-                        error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
-                               'rand_multivariate_student or rand_multivariate_normal as options']);
-                    else
-                        posterior_sampler_options.proposal_distribution=options_list{i,2};
-                    end
-                  case 'student_degrees_of_freedom'
-                    if options_list{i,2} <= 0
-                        error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
-                    else
-                        posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
-                    end
-                  case 'save_tmp_file'
-                    posterior_sampler_options.save_tmp_file = options_list{i,2};
-                  case 'number_of_particles'
-                    posterior_sampler_options.nparticles = options_list{i,2};
-                  otherwise
-                    warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
+            posterior_sampler_options.W1=posterior_sampler_options.initial_step_size*(bounds.ub-bounds.lb);
+            if options_.load_mh_file
+                posterior_sampler_options.slice_initialize_with_mode = 0;
+            else
+                if ~posterior_sampler_options.slice_initialize_with_mode
+                    posterior_sampler_options.invhess=[];
                 end
             end
-        end
 
-        options_.mode_compute = 0;
-        options_.cova_compute = 0;
-        options_.mh_replic = 0;
-        options_.mh_posterior_mode_estimation = true;
+            if ~isempty(posterior_sampler_options.mode_files) % multimodal case
+                modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
+                load(modes, 'xparams')
+                if size(xparams,2)<2
+                    error(['check_posterior_sampler_options:: Variable xparams loaded in file <' modes '> has size [' int2str(size(xparams,1)) 'x' int2str(size(xparams,2)) ']: it must contain at least two columns, to allow multi-modal sampling.'])
+                end
+                for j=1:size(xparams,2)
+                    mode(j).m=xparams(:,j);
+                end
+                posterior_sampler_options.mode = mode;
+                posterior_sampler_options.rotated = 1;
+                posterior_sampler_options.WR=[];
+            end
+
+          case 'hssmc'
+
+            % default options
+            posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
+
+            % user defined options
+            if ~isempty(options_.posterior_sampler_options.sampling_opt)
+                options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
+                for i=1:rows(options_list)
+                    switch options_list{i,1}
+                      case 'target'
+                        posterior_sampler_options.target = options_list{i,2};
+                      case 'steps'
+                        posterior_sampler_options.steps = options_list{i,2};
+                      case 'scale'
+                        posterior_sampler_options.scale = options_list{i,2};
+                      case 'particles'
+                        posterior_sampler_options.particles = options_list{i,2};
+                      case 'lambda'
+                        posterior_sampler_options.lambda = options_list{i,2};
+                      otherwise
+                        warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
+                    end
+                end
+            end
+
+            options_.mode_compute = 0;
+            options_.cova_compute = 0;
+            options_.mh_replic = 0;
+            options_.mh_posterior_mode_estimation = false;
 
       case 'dsmh'
 
-        posterior_sampler_options.parallel_bar_refresh_rate=1;
-        posterior_sampler_options.serial_bar_title='Dynamic Striated Metropolis Hastings';
-
         % default options
         posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
 
@@ -444,7 +432,7 @@ if init
                   case 'save_tmp_file'
                     posterior_sampler_options.save_tmp_file = options_list{i,2};
                   case 'number_of_particles'
-                    posterior_sampler_options.nparticles = options_list{i,2};
+                    posterior_sampler_options.particles = options_list{i,2};
                   otherwise
                     warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
                 end
@@ -456,11 +444,11 @@ if init
         options_.mh_replic = 0;
         options_.mh_posterior_mode_estimation = true;
 
-      otherwise
-        error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
-    end
+          otherwise
+            error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
+        end
 
-    return
+        return
 end
 
 % here are all samplers requiring a proposal distribution

matlab/default_option_values.m

@@ -492,21 +492,19 @@ options_.posterior_sampler_options.imh.proposal_distribution = 'rand_multivariat
 options_.posterior_sampler_options.imh.use_mh_covariance_matrix=0;
 options_.posterior_sampler_options.imh.save_tmp_file=0;
 % Herbst and Schorfeide SMC Sampler
-options_.posterior_sampler_options.hssmc.nphi= 25 ;
-options_.posterior_sampler_options.hssmc.lambda = 2 ;
-options_.posterior_sampler_options.hssmc.nparticles = 20000 ;
-options_.posterior_sampler_options.hssmc.c = 0.5 ;
-options_.posterior_sampler_options.hssmc.acpt = 0.25 ;
-options_.posterior_sampler_options.hssmc.trgt = 0.25 ;
-options_.posterior_sampler_options.hssmc.option_mutation = true ;
-options_.posterior_sampler_options.hssmc.alp = .9 ;
+options_.posterior_sampler_options.hssmc.steps= 25;
+options_.posterior_sampler_options.hssmc.lambda = 2;
+options_.posterior_sampler_options.hssmc.particles = 20000;
+options_.posterior_sampler_options.hssmc.scale = 0.5;
+options_.posterior_sampler_options.hssmc.acpt = 1.00;
+options_.posterior_sampler_options.hssmc.target = 0.25;
 % DSMH: Dynamic Striated Metropolis-Hastings algorithm
 options_.posterior_sampler_options.dsmh.H = 25 ;
 options_.posterior_sampler_options.dsmh.N = 20 ;
 options_.posterior_sampler_options.dsmh.G = 10 ;
 options_.posterior_sampler_options.dsmh.K = 50 ;
 options_.posterior_sampler_options.dsmh.lambda1 = 0.1 ;
-options_.posterior_sampler_options.dsmh.nparticles = 20000 ;
+options_.posterior_sampler_options.dsmh.particles = 20000 ;
 options_.posterior_sampler_options.dsmh.alpha0 = 0.2 ;
 options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
 options_.posterior_sampler_options.dsmh.tau = 10 ;
@@ -663,12 +661,6 @@ options_.use_dll = false;
 % model evaluated using bytecode.dll
 options_.bytecode = false;
 
-% if true, use a fixed point method to solve Sylvester equation (for large scale models)
-options_.sylvester_fp = false;
-
-% convergence criteria to solve iteratively a sylvester equations
-options_.sylvester_fixed_point_tol = 1e-12;
-
 % if true, use a fixed point method to solve Lyapunov equation (for large scale models)
 options_.lyapunov_fp = false;
 % if true, use a doubling algorithm to solve Lyapunov equation (for large scale models)

matlab/dyn_forecast.m

@@ -11,7 +11,7 @@ function forecast = dyn_forecast(var_list,M_,options_,oo_,task,dataset_info)
 %   task:         indicates how to initialize the forecast
 %                 either 'simul' or 'smoother'
 %   dataset_info:   Various informations about the dataset (descriptive statistics and missing observations).
-
+%
 % OUTPUTS
 %   forecast:   structure containing fields
 %                   Mean:       point estimate

matlab/dynare_estimation_1.m

@@ -393,7 +393,8 @@ end
 %
 
 if ishssmc(options_)
-    options_.posterior_sampler_options.hssmc.nphi=10;
+    [posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
+    options_.posterior_sampler_options.current_options = posterior_sampler_options;
     oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
 elseif isdsmh(options_)
     dsmh(objective_function, xparam1, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)

matlab/estimation/smc/dsmh.m

@@ -58,7 +58,7 @@ MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_
 
 % Step 0: Initialization of the sampler
 [param, tlogpost_iminus1, loglik, bayestopt_] = ...
-    smc_samplers_initialization(TargetFun, 'dsmh', opts.nparticles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
+    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 ;
@@ -79,20 +79,17 @@ 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.nparticles);
+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);
-%mat_var_cov = bsxfun(@minus,distrib_param,mean_xparam) ;
-%mat_var_cov = (mat_var_cov*mat_var_cov')/(options_.HSsmc.nparticles-1) ;
-%std_xparam = sqrt(diag(mat_var_cov)) ;
 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.nparticles) ;
-    ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.nparticles) ;
+    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;
@@ -131,9 +128,9 @@ 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.nparticles,bandwidth,kernel_function);
+    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.nparticles,optimal_bandwidth,kernel_function);
+        options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
     plot(density(:,1),density(:,2));
     hold on
     if TeX

matlab/estimation/smc/hssmc.m

@@ -34,7 +34,7 @@ function mdd = hssmc(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_,
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-    smcopt = options_.posterior_sampler_options.hssmc;
+    smcopt = options_.posterior_sampler_options.current_options;
 
     % Set location for the simulated particles.
     SimulationFolder = CheckPath('hssmc', M_.dname);
@@ -48,73 +48,72 @@ function mdd = hssmc(TargetFun, mh_bounds, dataset_, dataset_info, options_, M_,
     mlogit = @(x) .95 + .1/(1 + exp(-16*x)); % Update of the scale parameter
 
     % Create the tempering schedule
-    phi = ((0:smcopt.nphi-1)/(smcopt.nphi-1)).^smcopt.lambda;
+    phi = ((0:smcopt.steps-1)/(smcopt.steps-1)).^smcopt.lambda;
 
     % Initialise the estimate of the marginal density of the data
     mdd = .0;
 
     % tuning for MH algorithms matrices
-    scl = zeros(smcopt.nphi, 1);        % scale parameter
-    ESS = zeros(smcopt.nphi, 1);        % ESS
-    acpt = zeros(smcopt.nphi, 1);       % average acceptance rate
+    scl = zeros(smcopt.steps, 1);        % scale parameter
+    ESS = zeros(smcopt.steps, 1);        % ESS
+    acpt = zeros(smcopt.steps, 1);       % average acceptance rate
 
     % Initialization of the sampler (draws from the prior distribution with finite logged likelihood)
     t0 = tic;
     [particles, tlogpostkernel, loglikelihood] = ...
-        smc_samplers_initialization(funobj, 'hssmc', smcopt.nparticles, Prior, SimulationFolder, smcopt.nphi);
+        smc_samplers_initialization(funobj, 'hssmc', smcopt.particles, Prior, SimulationFolder, smcopt.steps);
     tt = toc(t0);
 
     dprintf('#Iter.       lambda        ESS          Acceptance rate   scale   resample   seconds')
     dprintf('%3u          %5.4f     %9.5E         %5.4f        %4.3f     %3s       %5.2f', 1, 0, 0, 0, 0, 'no', tt)
 
-    weights = ones(smcopt.nparticles, 1)/smcopt.nparticles;
+    weights = ones(smcopt.particles, 1)/smcopt.particles;
 
     resampled_particle_swarm = false;
 
-    for i=2:smcopt.nphi % Loop over the weight on the liklihood (phi)
+    for i=2:smcopt.steps % Loop over the weight on the liklihood (phi)
         weights = weights.*exp((phi(i)-phi(i-1))*loglikelihood);
         sweight = sum(weights);
         weights = weights/sweight;
         mdd = mdd + log(sweight);
         ESS(i) = 1.0/sum(weights.^2);
-        if (2*ESS(i) < smcopt.nparticles) % Resampling
+        if (2*ESS(i) < smcopt.particles) % Resampling
             resampled_particle_swarm = true;
             iresample = kitagawa(weights);
             particles = particles(:,iresample);
             loglikelihood = loglikelihood(iresample);
             tlogpostkernel = tlogpostkernel(iresample);
-            weights = ones(smcopt.nparticles, 1)/smcopt.nparticles;
+            weights = ones(smcopt.particles, 1)/smcopt.particles;
         end
-        smcopt.c = smcopt.c*mlogit(smcopt.acpt-smcopt.trgt); % Adjust the scale parameter
-        scl(i) = smcopt.c; % Scale parameter (for the jumping distribution in MH mutation step).
+        smcopt.scale = smcopt.scale*mlogit(smcopt.acpt-smcopt.target); % Adjust the scale parameter
+        scl(i) = smcopt.scale; % Scale parameter (for the jumping distribution in MH mutation step).
         mu = particles*weights; % Weighted average of the particles.
         z = particles-mu;
         R = z*(z'.*weights); % Weighted covariance matrix of the particles.
         t0 = tic;
-        acpt_ = zeros(smcopt.nparticles, 1);
+        acpt_ = false(smcopt.particles, 1);
         tlogpostkernel = tlogpostkernel + (phi(i)-phi(i-1))*loglikelihood;
         [acpt_, particles, loglikelihood, tlogpostkernel] = ...
-            randomwalk(funobj, chol(R, 'lower'), mu, smcopt.c, phi(i), acpt_, Prior, particles, loglikelihood, tlogpostkernel);
-        smcopt.acpt = sum(acpt_)/smcopt.nparticles; % Acceptance rate.
+            randomwalk(funobj, chol(R, 'lower'), mu, scl(i), phi(i), acpt_, Prior, particles, loglikelihood, tlogpostkernel);
+        smcopt.acpt = sum(acpt_)/smcopt.particles; % Acceptance rate.
         tt = toc(t0);
         acpt(i) = smcopt.acpt;
-        scl(i) = smcopt.c;
         if resampled_particle_swarm
             dprintf('%3u          %5.4f     %9.5E         %5.4f        %4.3f     %3s       %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'yes', tt)
         else
             dprintf('%3u          %5.4f     %9.5E         %5.4f        %4.3f     %3s       %5.2f', i, phi(i), ESS(i), acpt(i), scl(i), 'no', tt)
         end
-        if i==smcopt.nphi
+        if i==smcopt.steps
             iresample = kitagawa(weights);
             particles = particles(:,iresample);
         end
-        save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, smcopt.nphi), 'particles', 'tlogpostkernel', 'loglikelihood')
+        save(sprintf('%s%sparticles-%u-%u.mat', SimulationFolder, filesep(), i, smcopt.steps), 'particles', 'tlogpostkernel', 'loglikelihood')
         resampled_particle_swarm = false;
     end
 
 end
 
-function [acpt, particles, loglikelihood, tlogpostkernel] = randomwalk(funobj, RR, mu, scale, phi, acpt, Prior, particles, loglikelihood, tlogpostkernel)
+function [accept, particles, loglikelihood, tlogpostkernel] = randomwalk(funobj, RR, mu, scale, phi, accept, Prior, particles, loglikelihood, tlogpostkernel)
 
     parfor j=1:size(particles, 2)
         notvalid= true;
@@ -125,7 +124,7 @@ function [acpt, particles, loglikelihood, tlogpostkernel] = randomwalk(funobj, R
                 if isfinite(loglik)
                     notvalid = false;
                     if rand<exp(tlogpost-tlogpostkernel(j))
-                        acpt(j) = 1 ;
+                        accept(j) = true ;
                         particles(:,j) = candidate;
                         loglikelihood(j) = loglik;
                         tlogpostkernel(j) = tlogpost;

matlab/ftest.m

@@ -1,84 +0,0 @@
-function ftest (s1,s2)
-
-% Copyright © 2001-2017 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-global nvx nvy x y lag1
-
-if size(s1,1) ~= 2
-    error ('Spécifiez deux fichiers pour la comparaison.') ;
-end
-
-for i = 1:2
-    if ~ isempty(find(abs(s1(i,:)) == 46))
-        error ('Entrez les noms de fichiers sans extensions.') ;
-    end
-end
-
-s1 = [s1 [' ';' ']] ;
-file1 = [s1(1,1:min(find(abs(s1(1,:)) == 32))-1) '.BIN'] ;
-file2 = [s1(2,1:min(find(abs(s1(2,:)) == 32))-1) '.BIN'] ;
-
-fid=fopen(file1,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag1 = fread(fid,4,'int') ;
-nvx = fread(fid,[n1,n3],'int') ;
-x = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvx = char(nvx) ;
-
-fid=fopen(file2,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag2 = fread(fid,4,'int') ;
-nvy = fread(fid,[n1,n3],'int') ;
-y = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvy = char(nvy) ;
-
-if size(x,1) ~= size(y,1)
-    error ('FTEST: The two files don''t have the same number of variables.');
-end
-
-for i = 1:size(x,1)
-    if ~ strcmp(nvx(i,:),nvy(i,:))
-        error ('FTEST: The two files don''t have the same  variables.') ;
-    end
-end
-
-if nnz(lag1 - lag2) > 0
-    error ('FTEST: Leads and lags aren''t the same in both files.') ;
-end
-
-j = zeros(size(s2,1),1);
-for i=1:size(s2,1)
-    k = strmatch(s2(i,:),nvx,'exact') ;
-    if isempty(k)
-        t = ['FTEST: Variable ' s2(i) 'doesn''t exist'] ;
-        error (t) ;
-    else
-        j(i) =k;
-    end
-end
-
-y = y(j,:) ;
-x = x(j,:) ;
-
-%06/18/01 MJ replaced beastr by strmatch

matlab/gcompare.m

@@ -1,52 +0,0 @@
-function gcompare(s1,s2)
-% GCOMPARE :    GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
-%               This optional command plots the trajectories of a list of
-%               variables in two different simulations. One plot is drawn
-%               for each variable. The trajectories must have been previously
-%               saved by the instruction DYNASAVE. The simulation in file1
-%               is refered to as the base simulation.
-
-% Copyright © 2001-2017 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
-
-global options_ M_
-global nvx nvy x y lag1
-
-ftest(s1,s2) ;
-
-ix = [1-lag1(1):size(x,2)-lag1(1)]' ;
-i = [lag1(1):size(ix,1)-lag1(2)+1]' ;
-
-if options_.smpl == 0
-    i = [M_.maximum_lag:size(y,2)]' ;
-else
-    i = [options_.smpl(1):options_.smpl(2)]' ;
-end
-
-for k = 1:size(x,1)
-    figure ;
-    plot (ix(i),x(k,i),ix(i),y(k,i)) ;
-    xlabel (['Periods']) ;
-    title (['Variable ' s2(k,:)]) ;
-    l = min(i) + 1;
-    ll = max(i) - 1 ;
-    text (l,x(k,l),s1(1,:)) ;
-    text (ll,y(k,ll),s1(2,:)) ;
-end
-
-% 06/18/01 MJ corrected treatment of options_.smpl
-% 06/24/01 MJ removed color specification

matlab/get_Hessian.m

@@ -1,255 +0,0 @@
-function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
-% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
-%
-% computes the hessian matrix of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% Thanks to  Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
-% two dimensions equal to the number of structural parameters
-
-% Copyright © 2011-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/>.
-
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-pp   = size(Y,1);                               % Maximum number of observed variables.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-Hess  = zeros(k,k);                             % Initialization of the Hessian
-Da    = zeros(mm,k);                             % State vector.
-Dv = zeros(length(mf),k);
-D2a    = zeros(mm,k,k);                             % State vector.
-D2v = zeros(length(mf),k,k);
-
-C = zeros(length(mf),mm);
-for ii=1:length(mf); C(ii,mf(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
-dC = zeros(length(mf),mm,k);
-d2C = zeros(length(mf),mm,k,k);
-
-s   = zeros(pp,1);                      % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
-ds  = zeros(pp,1,k);
-d2s = zeros(pp,1,k,k);
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady & t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-        [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
-        tmp = (a+K*v);
-
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
-            %  dai = da(:,:,ii);
-            dKi  = DK(:,:,ii);
-            diFi = -iF*DF(:,:,ii)*iF;
-            dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
-
-
-            for jj = 1:ii
-                dFj    = DF(:,:,jj);
-                diFj   = -iF*DF(:,:,jj)*iF;
-                dKj  = DK(:,:,jj);
-                d2Kij  = D2K(:,:,jj,ii);
-                d2Fij  = D2F(:,:,jj,ii);
-                d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-                dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
-
-                d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
-                d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
-                D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
-
-                Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
-            end
-            Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
-        end
-        %                     vecDPmf = reshape(DP(mf,mf,:),[],k);
-        %                     iPmf = inv(P(mf,mf));
-        if t>=start
-            Hess = Hess + Hesst;
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-        D2P     = D2P1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        tmp = (a+K*v);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            dKi  = DK(:,:,ii);
-            diFi = -iF*DF(:,:,ii)*iF;
-            dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
-
-            for jj = 1:ii
-                dFj    = DF(:,:,jj);
-                diFj   = -iF*DF(:,:,jj)*iF;
-                dKj  = DK(:,:,jj);
-                d2Kij  = D2K(:,:,jj,ii);
-                d2Fij  = D2F(:,:,jj,ii);
-                d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-                dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
-
-                d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
-                d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
-                D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
-
-                Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
-            end
-            Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
-        end
-        if t>=start
-            Hess = Hess + Hesst;
-        end
-        a = T*(a+K*v);
-    end
-    %         Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
-    %         for ii = 1:k;
-    %             for jj = 1:ii
-    %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
-    %             end
-    %         end
-end
-
-Hess = Hess + tril(Hess,-1)';
-
-Hess = -Hess/2;
-% end of main function
-
-function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
-% computes (i,j) term in the Hessian
-
-Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
-
-% end of getHesst_ij
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman
-
-function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0)
-% computes the second derivatives of the Kalman matrices
-% note: A=T in main func.
-
-k      = size(dA,3);
-tmp    = P0-K0*P0(mf,:);
-[ns,no] = size(K0);
-
-% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
-% APC = A*P0*C';
-% APA = A*P0*A';
-
-
-d2K  = zeros(ns,no,k,k);
-d2S  = zeros(no,no,k,k);
-d2P1 = zeros(ns,ns,k,k);
-
-for ii = 1:k
-    dAi = dA(:,:,ii);
-    dFi = dP0(mf,mf,ii);
-    d2Omi = d2Om(:,:,ii);
-    diFi = -iF*dFi*iF;
-    dKi = dK0(:,:,ii);
-    for jj = 1:k
-        dAj = dA(:,:,jj);
-        dFj = dP0(mf,mf,jj);
-        d2Omj = d2Om(:,:,jj);
-        dFj = dP0(mf,mf,jj);
-        diFj = -iF*dFj*iF;
-        dKj = dK0(:,:,jj);
-
-        d2Aij = d2A(:,:,jj,ii);
-        d2Pij = d2P0(:,:,jj,ii);
-        d2Omij = d2Om(:,:,jj,ii);
-
-        % second order
-
-        d2Fij = d2Pij(mf,mf) ;
-
-        %     d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi';
-        d2APC = d2Pij(:,mf);
-
-        d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-
-        d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj;
-
-        d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ;
-
-        dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii);
-        dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj);
-        d2tmp = d2Pij - d2KCP;
-
-        d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi';
-
-        d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
-        d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
-        d2S(:,:,ii,jj)  = d2Fij;
-        %     d2iS(:,:,ii,jj) = d2iF;
-    end
-end
-
-% end of computeD2Kalman

matlab/gsa/smirnov.m

@@ -34,11 +34,16 @@ end
 % empirical cdfs.
 xmix= [x1;x2];
 bin = [-inf ; sort(xmix) ; inf];
-
-ncount1 = histc (x1 , bin);
-ncount1 = ncount1(:);
-ncount2 = histc (x2 , bin);
-ncount2 = ncount2(:);
+if isoctave
+    ncount1 = histc(x1 , bin);
+else
+    ncount1 = histcounts(x1 , bin);
+end
+if isoctave
+    ncount2 = histc(x2 , bin);
+else
+    ncount2 = histcounts(x2 , bin);
+end
 
 cum1  =  cumsum(ncount1)./sum(ncount1);
 cum1  =  cum1(1:end-1);

matlab/identification_checks.m

@@ -170,7 +170,7 @@ end
 jweak      = zeros(1,param_nbr);
 jweak_pair = zeros(param_nbr,param_nbr);
 
-if test_flag ~= 0 || test_flag ~= 0
+if test_flag ~= 0
     % these tests only apply to Jacobians, not to Gram matrices, i.e. Hessian-type or 'covariance' matrices
     Pco = NaN(param_nbr,param_nbr);
     for ii = 1:size(Xparnonzero,2)

matlab/load_last_mh_history_file.m

@@ -7,10 +7,8 @@ function record = load_last_mh_history_file(MetropolisFolder, ModelName)
 % Outputs:
 %   record              [struct]    structure storing the MH history
 %
-% Notes: The record structure is written to the caller workspace via an
-% assignin statement.
 
-% Copyright © 2013-2017 Dynare Team
+% Copyright © 2013-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +34,7 @@ mh_history_files = dir([BaseName '_mh_history_*.mat']);
 
 % Consistency with older versions of Dynare.
 if isequal(length(mh_history_files),0)
-    if exist([BaseName '_mh_history.mat'])
+    if exist([BaseName '_mh_history.mat'],'file')
         format_mh_history_file = 1; % old Dynare format
     else
         error(['Estimation::load_mh_file: I cannot find any mh-history file in ' MetropolisFolder '!'])
@@ -46,7 +44,7 @@ else
 end
 
 if format_mh_history_file %needed to preserve backward compatibility
-    load([BaseName '_mh_history.mat']);
+    load([BaseName '_mh_history.mat'],'record');
     record.LastLogPost = record.LastLogLiK;
     record.InitialLogPost = record.InitialLogLiK;
     record.LastSeeds = record.Seeds;
@@ -60,7 +58,7 @@ if format_mh_history_file %needed to preserve backward compatibility
     record = rmfield(record,'AcceptationRates');
     save([BaseName '_mh_history_0.mat'],'record');
 else
-    load([BaseName '_mh_history_' num2str(length(mh_history_files)-1) '.mat']);
+    load([BaseName '_mh_history_' num2str(length(mh_history_files)-1) '.mat'],'record');
     % add fields that have later been introduced
     if ~isfield(record,'MCMCConcludedSuccessfully')
         record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...

matlab/mcpath_function.m

@@ -1,8 +0,0 @@
-function [res,jac,domerr] = mcpath_function(func,z,jacflag,varargin)
-domerr = 0;
-
-if jacflag
-    [res,jac] = func(z,varargin{:});
-else
-    res = func(z,varargin{:});
-end

matlab/metropolis_run_analysis.m

@@ -1,55 +0,0 @@
-function metropolis_run_analysis(M_,basetopt,j)
-%function metropolis_run_analysis(M_,basetopt,j
-% analizes Metropolis runs
-%
-% INPUTS
-%   M_:        (struct)  Model structure
-%   basetopt:  (struct)  Estimated parameters structure
-%   j:         (int)     Index of estimated paramter
-%
-% OUTPUTS
-%   none
-%
-% SPECIAL REQUIREMENTS
-%   none
-
-% Copyright © 2003-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/>.
-
-load([M_.fname '/metropolis/' M_.fname '_mh_history'])
-nblck = record.Nblck;
-ndraws = sum(record.MhDraws(:,1));
-
-logPost = [];
-params = [];
-blck = 1;
-for i=1:record.LastFileNumber
-    fname = [M_.fname '/metropolis/' M_.fname '_mh' int2str(i) '_blck' ...
-             int2str(blck) '.mat'];
-    if exist(fname,'file')
-        o=load(fname);
-        logPost = [logPost; o.logpo2];
-        params  = [params; o.x2];
-    end
-end
-
-figure;
-subplot(2,1,1)
-plot(logPost)
-subplot(2,1,2)
-plot(params(:,j))
-title(['parameter ' int2str(j)])

matlab/modules/dseries

@@ -1 +1 @@
-Subproject commit 80446e7cdcf392167e91b428b61d77e028e47674
+Subproject commit a3874403fe5b403f74d727006cb4ccf198cf63c8

matlab/plot_shock_decomposition.m

@@ -536,6 +536,11 @@ b = size(z,3);
 if ~isempty(options_.plot_shock_decomp.plot_init_date)
     my_initial_date = max(initial_date,options_.plot_shock_decomp.plot_init_date);
     a = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_init_date);
+    if isempty(a)
+        warning(['You set plot_init_date larger than the last observation. The last observation is %s,\n' ...
+            'while you requested as the first observation %s.Exiting because there is nothing to do.'],initial_date+b,options_.plot_shock_decomp.plot_init_date);
+        return;
+    end
 end
 if ~isempty(options_.plot_shock_decomp.plot_end_date)
     if options_.plot_shock_decomp.plot_end_date<=(max(initial_date:initial_date+b-1))

matlab/score.m

@@ -1,123 +0,0 @@
-function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-% function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-%
-% computes the derivative of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% thanks to Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-
-% Copyright © 2009-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 <http://www.gnu.org/licen
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-DLIK  = zeros(k,1);                             % Initialization of the score.
-Da    = zeros(mm,k);                            % State vector.
-Dv    = zeros(length(mf),k);                    % observation vector.
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady & t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-for ii = 1:k
-    tmp0(:,:,ii) = iF*DF(:,:,ii)*iF;
-end
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        a = T*(a+K*v);
-    end
-    for ii = 1:k
-        %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
-    end
-
-end
-
-DLIK = DLIK/2;
-
-% end of main function
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman

matlab/shiftS.m

@@ -1,28 +0,0 @@
-function S = shiftS(S,n)
-%function S = shiftS(S,n)
-%
-% Removes the first n elements of a one dimensional cell array.
-
-% Copyright © 2013-2014 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 length(S) >= n+1
-    S = S(n+1:end);
-else
-    S = {};
-end
-end

matlab/var_sample_moments.m

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

meson.build

@@ -457,7 +457,8 @@ shared_module('k_order_welfare', k_order_welfare_src, kwargs : korder_mex_kwargs
 
 # Unit tests
 
-k_order_test_kwargs = { 'include_directories' : [ mex_incdir, korder_incdir ],
+k_order_test_kwargs = { 'include_directories' : mex_kwargs['include_directories'] + korder_incdir,
+                        'cpp_args' : mex_kwargs['cpp_args'],
                         'link_args' : exe_link_args,
                         'build_rpath' : exe_rpath,
                         'link_with' : korder_lib }

mex/sources/libkorder/k_ord_dynare.cc

@@ -1,5 +1,5 @@
 /*
- * Copyright © 2008-2022 Dynare Team
+ * Copyright © 2008-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -248,7 +248,7 @@ DynareStateNameList::DynareStateNameList(const KordpDynare& dynare, const Dynare
                                          const DynareNameList& denl)
 {
   for (int i = 0; i < dynare.nYs; i++)
-    names.emplace_back(dnl.getName(i + dynare.nstat()));
-  for (int i = 0; i < dynare.nexog(); i++)
+    names.emplace_back(dnl.getName(i + dynare.nStat));
+  for (int i = 0; i < dynare.nExog; i++)
     names.emplace_back(denl.getName(i));
 }

mex/sources/libkorder/sylv/GeneralMatrix.hh

@@ -90,12 +90,10 @@ public:
   ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
   // Create submatrix (with data sharing)
   ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols);
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit ConstGeneralMatrix(const mxArray* p) :
       data(p), rows {static_cast<int>(mxGetM(p))}, cols {static_cast<int>(mxGetN(p))}, ld {rows}
   {
   }
-#endif
   virtual ~ConstGeneralMatrix() = default;
 
   ConstGeneralMatrix& operator=(const ConstGeneralMatrix& v) = delete;
@@ -234,12 +232,10 @@ public:
   // Create submatrix (with data sharing)
   GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols);
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit GeneralMatrix(mxArray* p) :
       data(p), rows {static_cast<int>(mxGetM(p))}, cols {static_cast<int>(mxGetN(p))}, ld {rows}
   {
   }
-#endif
 
   virtual ~GeneralMatrix() = default;
   GeneralMatrix& operator=(const GeneralMatrix& m) = default;

mex/sources/libkorder/sylv/SylvParams.cc

@@ -125,7 +125,6 @@ SylvParams::getArrayNames() const
   return names;
 }
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 mxArray*
 SylvParams::DoubleParamItem::createMatlabArray() const
 {
@@ -136,11 +135,11 @@ mxArray*
 SylvParams::IntParamItem::createMatlabArray() const
 {
   mxArray* res = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
-# if MX_HAS_INTERLEAVED_COMPLEX
+#if MX_HAS_INTERLEAVED_COMPLEX
   *mxGetInt32s(res) = value;
-# else
+#else
   *static_cast<int*>(mxGetData(res)) = value;
-# endif
+#endif
   return res;
 }
 
@@ -225,4 +224,3 @@ SylvParams::createStructArray() const
 
   return res;
 }
-#endif

mex/sources/libkorder/sylv/SylvParams.hh

@@ -25,9 +25,7 @@
 #include <string>
 #include <vector>
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
-# include <dynmex.h>
-#endif
+#include <dynmex.h>
 
 enum class status
 {
@@ -111,9 +109,7 @@ protected:
       ParamItem<double>::operator=(val);
       return *this;
     }
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
     [[nodiscard]] mxArray* createMatlabArray() const;
-#endif
   };
 
   class IntParamItem : public ParamItem<int>
@@ -132,9 +128,7 @@ protected:
       ParamItem<int>::operator=(val);
       return *this;
     }
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
     [[nodiscard]] mxArray* createMatlabArray() const;
-#endif
   };
 
   class BoolParamItem : public ParamItem<bool>
@@ -153,9 +147,7 @@ protected:
       ParamItem<bool>::operator=(val);
       return *this;
     }
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
     [[nodiscard]] mxArray* createMatlabArray() const;
-#endif
   };
 
   class MethodParamItem : public ParamItem<solve_method>
@@ -174,9 +166,7 @@ protected:
       ParamItem<solve_method>::operator=(val);
       return *this;
     }
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
     [[nodiscard]] mxArray* createMatlabArray() const;
-#endif
   };
 
 public:
@@ -224,9 +214,8 @@ public:
   void print(const std::string& prefix) const;
   void print(std::ostream& fdesc, const std::string& prefix) const;
   [[nodiscard]] std::vector<const char*> getArrayNames() const;
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   [[nodiscard]] mxArray* createStructArray() const;
-#endif
+
 private:
   void copy(const SylvParams& p);
 };

mex/sources/libkorder/sylv/Vector.cc

@@ -107,14 +107,12 @@ Vector::Vector(const Vector& v, int off_arg, int skip, int l) : len(l), data {ne
   copy(v.data + off_arg * v.s, v.s * skip);
 }
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 Vector::Vector(mxArray* p) :
     len {static_cast<int>(mxGetNumberOfElements(p))}, data {mxGetPr(p)}, destroy {false}
 {
   if (!mxIsDouble(p) || mxIsComplex(p) || mxIsSparse(p))
     throw SYLV_MES_EXCEPTION("This is not a dense array of real doubles.");
 }
-#endif
 
 bool
 Vector::operator==(const Vector& y) const
@@ -305,14 +303,12 @@ ConstVector::ConstVector(const double* d, int skip, int l) : len {l}, s {skip},
 {
 }
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 ConstVector::ConstVector(const mxArray* p) :
     len {static_cast<int>(mxGetNumberOfElements(p))}, data {mxGetPr(p)}
 {
   if (!mxIsDouble(p))
     throw SYLV_MES_EXCEPTION("This is not a MATLAB array of doubles.");
 }
-#endif
 
 bool
 ConstVector::operator==(const ConstVector& y) const

mex/sources/libkorder/sylv/Vector.hh

@@ -28,9 +28,7 @@
 #include <complex>
 #include <utility>
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
-# include <dynmex.h>
-#endif
+#include <dynmex.h>
 
 class GeneralMatrix;
 class ConstVector;
@@ -72,9 +70,7 @@ public:
   Vector(const Vector& v, int off_arg, int l);
   Vector(Vector& v, int off_arg, int skip, int l);
   Vector(const Vector& v, int off_arg, int skip, int l);
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit Vector(mxArray* p);
-#endif
   Vector& operator=(const Vector& v);
   /* The move-assignment operator is not implemented, because moving pointers
      across class instances would break the “reference semantics” that the
@@ -198,9 +194,7 @@ public:
   ConstVector(const ConstVector& v, int off_arg, int l);
   ConstVector(const ConstVector& v, int off_arg, int skip, int l);
   ConstVector(const double* d, int skip, int l);
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit ConstVector(const mxArray* p);
-#endif
   virtual ~ConstVector() = default;
   ConstVector& operator=(const ConstVector& v) = delete;
   ConstVector& operator=(ConstVector&& v) = delete;

mex/sources/libkorder/tl/t_polynomial.hh

@@ -216,7 +216,7 @@ public:
                   auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
                   ten_smart->zeros();
                   ten = ten_smart.get();
-                  insert(std::move(ten_smart));
+                  _Tparent::insert(std::move(ten_smart));
                 }
 
               Symmetry sym {i, j};
@@ -245,7 +245,7 @@ public:
               auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
               ten_smart->zeros();
               ten = ten_smart.get();
-              insert(std::move(ten_smart));
+              _Tparent::insert(std::move(ten_smart));
             }
 
           Symmetry sym {0, j};

mex/sources/libkorder/tl/twod_matrix.hh

@@ -1,6 +1,6 @@
 /*
  * Copyright © 2004-2011 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -63,11 +63,9 @@ public:
       ConstGeneralMatrix(m, first_row, first_col, rows, cols)
   {
   }
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit ConstTwoDMatrix(const mxArray* p) : ConstGeneralMatrix(p)
   {
   }
-#endif
   ~ConstTwoDMatrix() override = default;
 
   ConstTwoDMatrix& operator=(const ConstTwoDMatrix& v) = delete;
@@ -125,11 +123,9 @@ public:
   {
   }
 
-#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   explicit TwoDMatrix(mxArray* p) : GeneralMatrix(p)
   {
   }
-#endif
   ~TwoDMatrix() override = default;
 
   TwoDMatrix& operator=(const TwoDMatrix& m) = default;

preprocessor

@@ -1 +1 @@
-Subproject commit 378d00fc3a50a11c3dce615ea337ad55989c938e
+Subproject commit cd86b1895d97d1b195849072a8d26fd372ca5d8c

scripts/dynare.el

@@ -55,39 +55,40 @@
 ;; Also include "end" in this list
 (defvar dynare-statements
   '("var" "varexo" "varexo_det" "trend_var" "log_trend_var"
-    "predetermined_variables" "parameters" "model_local_variable"
-    "model_info" "estimation" "set_time" "data" "varobs"
-    "varexobs" "unit_root_vars" "rplot" "osr_params" "osr" "dynatype"
-    "dynasave" "model_comparison" "change_type" "load_params_and_steady_state"
-    "save_params_and_steady_state" "write_latex_dynamic_model"
-    "write_latex_static_model" "write_latex_original_model"
-    "write_latex_steady_state_model" "steady" "check" "simul" "stoch_simul"
-    "var_model" "trend_component_model" "var_expectation_model" "pac_model"
-    "dsample" "planner_objective" "ramsey_model" "ramsey_policy"
-    "evaluate_planner_objective" "occbin_setup" "occbin_solver"
-    "occbin_write_regimes" "occbin_graph"
+    "predetermined_variables" "parameters" "model_local_variable" "model_info"
+    "estimation" "set_time" "data" "varobs" "varexobs" "unit_root_vars" "rplot"
+    "osr_params" "osr" "dynatype" "dynasave" "model_comparison" "change_type"
+    "load_params_and_steady_state" "save_params_and_steady_state"
+    "write_latex_dynamic_model" "write_latex_static_model"
+    "write_latex_original_model" "write_latex_steady_state_model" "steady"
+    "check" "simul" "stoch_simul" "var_model" "trend_component_model"
+    "var_expectation_model" "pac_model" "dsample" "planner_objective"
+    "ramsey_model" "ramsey_policy" "evaluate_planner_objective" "occbin_setup"
+    "occbin_solver" "occbin_write_regimes" "occbin_graph"
     "discretionary_policy" "identification" "bvar_density" "bvar_forecast"
-    "dynare_sensitivity" "initval_file" "histval_file" "forecast"
+    "bvar_irf" "dynare_sensitivity" "initval_file" "histval_file" "forecast"
     "shock_decomposition" "realtime_shock_decomposition"
     "plot_shock_decomposition" "initial_condition_decomposition"
-    "squeeze_shock_decomposition" "sbvar"
-    "ms_estimation" "ms_simulation" "ms_compute_mdd" "ms_compute_probabilities"
-    "ms_forecast" "ms_irf" "ms_variance_decomposition" "conditional_forecast"
-    "plot_conditional_forecast" "method_of_moments"
-    "markov_switching" "svar" "svar_global_identification_check"
-    "external_function" "calib_smoother" "model_diagnostics" "extended_path"
-    "smoother2histval" "perfect_foresight_setup" "perfect_foresight_solver"
+    "squeeze_shock_decomposition" "sbvar" "ms_estimation" "ms_simulation"
+    "ms_compute_mdd" "ms_compute_probabilities" "ms_forecast" "ms_irf"
+    "ms_variance_decomposition" "conditional_forecast"
+    "plot_conditional_forecast" "method_of_moments" "markov_switching" "svar"
+    "svar_global_identification_check" "external_function" "calib_smoother"
+    "model_diagnostics" "extended_path" "smoother2histval"
+    "perfect_foresight_setup" "perfect_foresight_solver"
     "perfect_foresight_with_expectation_errors_setup"
-    "perfect_foresight_with_expectation_errors_solver"
-    "compilation_setup" "model_remove" "model_options" "var_remove" "resid"
-    "std" "corr" "prior_function" "posterior_function" "end")
+    "perfect_foresight_with_expectation_errors_solver" "compilation_setup"
+    "model_remove" "model_options" "var_remove" "resid" "std" "corr"
+    "prior_function" "posterior_function" "end")
   "Dynare statement keywords.")
 
 ;; Keywords that may appear in blocks, and that begin a statement which will be
 ;; closed by a semicolon
 (defvar dynare-statements-like
-  '("stderr" "values" "periods" "scales" "restriction" "exclusion" "upper_cholesky" "lower_cholesky"
-    "bind" "relax" "error_bind" "error_relax" "add" "multiply")
+  '("stderr" "values" "periods" "scales" "restriction" "exclusion"
+    "upper_cholesky" "lower_cholesky" "equation" "bind" "relax" "error_bind"
+    "error_relax" "add" "multiply" "target" "auxname_target_nonstationary"
+    "component" "growth" "auxname" "kind")
   "Dynare statements-like keywords.")
 
 ;; Those keywords that makes the lexer enter the DYNARE_BLOCK start condition
@@ -96,13 +97,15 @@
 ;; referenced in another eval-when-compile statement in dynare-calculate-indentation
 (eval-when-compile
   (defvar dynare-blocks
-    '("model" "steady_state_model" "initval" "endval" "histval" "shocks" "heteroskedastic_shocks"
-      "shock_groups" "init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
-      "estimated_params_init" "estimated_params_bounds" "estimated_params_remove"
-      "osr_params_bounds" "observation_trends" "deterministic_trends" "optim_weights"
-      "homotopy_setup" "conditional_forecast_paths" "svar_identification"
-      "moment_calibration" "irf_calibration" "ramsey_constraints" "generate_irfs"
-      "matched_moments" "occbin_constraints" "model_replace" "verbatim" "pac_target_info")
+    '("model" "steady_state_model" "initval" "endval" "histval"
+      "filter_initial_state" "shocks" "heteroskedastic_shocks" "shock_groups"
+      "init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
+      "estimated_params_init" "estimated_params_bounds"
+      "estimated_params_remove" "osr_params_bounds" "observation_trends"
+      "deterministic_trends" "optim_weights" "homotopy_setup"
+      "conditional_forecast_paths" "svar_identification" "moment_calibration"
+      "irf_calibration" "ramsey_constraints" "generate_irfs" "matched_moments"
+      "occbin_constraints" "model_replace" "pac_target_info" "verbatim")
     "Dynare block keywords."))
 
 ;; Mathematical functions and operators used in model equations (see "hand_side" in Bison file)

tests/estimation/hssmc/fs2000.mod

@@ -84,4 +84,8 @@ options_.solve_tolf = 1e-12;
 
 estimation(order=1, datafile='../fsdat_simul.m', nobs=192, loglinear,
            posterior_sampling_method='hssmc',
-           posterior_sampler_options=('nphi',10));
+           posterior_sampler_options=('steps',10,
+                                      'lambda',2,
+                                      'particles', 20000,
+                                      'scale',.5,
+                                      'target', .25));

tests/estimation/method_of_moments/AnScho/AnScho_MoM_common.inc

@@ -269,8 +269,6 @@ method_of_moments(
     % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
     % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
     % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
-    % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
-    % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
     % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
     % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
     % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

tests/estimation/method_of_moments/RBC/RBC_MoM_Andreasen.mod

@@ -218,8 +218,6 @@ method_of_moments(
         % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
         % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
         % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
-        % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
-        % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
         % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
         % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
         % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

tests/estimation/method_of_moments/RBC/RBC_MoM_SMM_ME.mod

@@ -205,8 +205,6 @@ end
         % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
         % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
         % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
-        % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
-        % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
         % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
         % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
         % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

tests/estimation/method_of_moments/RBC/RBC_MoM_prefilter.mod

@@ -177,8 +177,6 @@ save('test_matrix.mat','weighting_matrix')
         % , lyapunov_complex_threshold = 1e-15     % complex block threshold for the upper triangular matrix in symmetric Lyapunov equation solver
         % , lyapunov_fixed_point_tol = 1e-10       % convergence criterion used in the fixed point Lyapunov solver
         % , lyapunov_doubling_tol = 1e-16          % convergence criterion used in the doubling algorithm
-        % , sylvester = default                    % algorithm to solve Sylvester equation; possible values are DEFAULT, FIXED_POINT
-        % , sylvester_fixed_point_tol = 1e-12      % convergence criterion used in the fixed point Sylvester solver
         % , qz_criterium = 0.999999                % value used to split stable from unstable eigenvalues in reordering the Generalized Schur decomposition used for solving first order problems
         % , qz_zero_threshold = 1e-6               % value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition
         % , schur_vec_tol=1e-11                    % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix

windows/build.sh

@@ -53,15 +53,15 @@ ln -s "$ROOT_DIRECTORY"/deps/mkoctfile64 /tmp/windeps/
 # Go to source root directory
 cd ..
 
-common_meson_opts=(-Dbuildtype=release --cross-file scripts/windows-cross.ini)
+common_meson_opts=(-Dbuildtype=release --cross-file windows/mingw-cross.ini)
 
 # Create Windows 64-bit DLL binaries for MATLAB ≥ R2018b
-meson setup --cross-file scripts/windows-cross-matlab.ini -Dmatlab_path=/tmp/windeps/matlab64/R2018b \
+meson setup --cross-file windows/mingw-cross-matlab.ini -Dmatlab_path=/tmp/windeps/matlab64/R2018b \
       "${common_meson_opts[@]}" build-win-matlab
 meson compile -v -C build-win-matlab
 
 # Create Windows DLL binaries for Octave/MinGW (64bit)
-meson setup --cross-file scripts/windows-cross-octave.ini \
+meson setup --cross-file windows/mingw-cross-octave.ini \
       "${common_meson_opts[@]}" build-win-octave
 meson compile -v -C build-win-octave
 

windows/deps/versions.mk

@@ -18,7 +18,7 @@ MATLAB64_VERSION = 20231122
 ## Build dependencies
 
 # pacman -Ss mingw-w64-x86_64-boost
-MINGW64_BOOST_VERSION = 1.83.0-1
+MINGW64_BOOST_VERSION = 1.83.0-2
 
 # pacman -Ss mingw-w64-x86_64-gsl
 MINGW64_GSL_VERSION = 2.7.1-2
@@ -40,11 +40,11 @@ MINGW64_LIBAEC_VERSION = 1.0.6-2
 
 # Dependency of HDF5
 # pacman -Ss mingw-w64-x86_64-openssl
-MINGW64_OPENSSL_VERSION = 3.1.4-1
+MINGW64_OPENSSL_VERSION = 3.2.0-1
 
 # Dependency of HDF5
 # pacman -Ss mingw-w64-x86_64-curl
-MINGW64_CURL_VERSION = 8.4.0-1
+MINGW64_CURL_VERSION = 8.5.0-1
 
 # Dependency of curl (and of the MinGW compiler)
 # pacman -Ss mingw-w64-x86_64-zstd
@@ -81,22 +81,22 @@ MINGW64_LIBUNISTRING_VERSION = 1.1-1
 ## MinGW packages for the embedded compiler
 
 # pacman -Ss mingw-w64-x86_64-gcc$
-MINGW64_GCC_VERSION = 13.2.0-2
+MINGW64_GCC_VERSION = 13.2.0-3
 
 # pacman -Ss mingw-w64-x86_64-gmp
 MINGW64_GMP_VERSION = 6.3.0-2
 
 # pacman -Ss mingw-w64-x86_64-binutils
-MINGW64_BINUTILS_VERSION = 2.41-2
+MINGW64_BINUTILS_VERSION = 2.41-3
 
 # pacman -Ss mingw-w64-x86_64-headers-git
-MINGW64_HEADERS_VERSION = 11.0.0.r404.g3a137bd87-1
+MINGW64_HEADERS_VERSION = 11.0.0.r442.ga27e7b27e-1
 
 # pacman -Ss mingw-w64-x86_64-crt-git
-MINGW64_CRT_VERSION = 11.0.0.r404.g3a137bd87-1
+MINGW64_CRT_VERSION = 11.0.0.r442.ga27e7b27e-1
 
 # pacman -Ss mingw-w64-x86_64-winpthreads-git
-MINGW64_WINPTHREADS_VERSION = 11.0.0.r404.g3a137bd87-1
+MINGW64_WINPTHREADS_VERSION = 11.0.0.r442.ga27e7b27e-1
 
 # pacman -Ss mingw-w64-x86_64-isl
 MINGW64_ISL_VERSION = 0.26-1

windows/mingw-cross.ini

@@ -1,5 +1,5 @@
 # Meson common cross file for targeting Windows from Linux
-# To be included after windows-cross-{matlab,octave}.ini
+# To be included after mingw-cross-{matlab,octave}.ini
 
 [constants]
 # For the architectural baseline, we follow MSYS2: