macOS/homebrew-native-arm64.ini

 # Meson native file for compiling under Homebrew for arm64 architecture
 
 [binaries]
-cpp = '/opt/homebrew/bin/g++-14'
-c = '/opt/homebrew/bin/gcc-14'
+cpp = '/opt/homebrew/bin/g++-15'
+c = '/opt/homebrew/bin/gcc-15'
 flex = '/opt/homebrew/opt/flex/bin/flex'
 bison = '/opt/homebrew/opt/bison/bin/bison'
 

macOS/homebrew-native-x86_64.ini

 # Meson native file for compiling under Homebrew for x86_64 architecture
 
 [binaries]
-cpp = '/usr/local/bin/g++-14'
-c = '/usr/local/bin/gcc-14'
+cpp = '/usr/local/bin/g++-15'
+c = '/usr/local/bin/gcc-15'
 flex = '/usr/local/opt/flex/bin/flex'
 bison = '/usr/local/opt/bison/bin/bison'
 

matlab/+bgp/write.m

@@ -7,7 +7,7 @@ function write(M_)
 % sharing the same trend(s) are relevant.
 %
 % INPUTS
-% - M_     [struct]   Dynare generated stucture describing the model (M_).
+% - M_     [struct]   Dynare generated structure describing the model
 %
 % OUTPUTS
 % None
@@ -15,7 +15,7 @@ function write(M_)
 % REMARKS
 % - The trends are assumed to be multiplicative.
 
-% Copyright © 2019-2023 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -75,7 +75,7 @@ fid = fopen(sprintf('+%s/bgpfun.m', M_.fname), 'w');
 
 % Write header. 
 fprintf(fid, 'function [F, JAC] = bgpfun(z)\n\n');
-fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
+fprintf(fid, '%% This file has been generated by Dynare (%s).\n\n', datestr(now));
 
 % The function admits a unique vector as input argument. The first
 % half of the elements are for the levels of the endogenous
@@ -84,37 +84,40 @@ fprintf(fid, 'y = z(1:%u);\n\n', M_.endo_nbr);
 fprintf(fid, 'g = z(%u:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr);
 
 % Define the point where the dynamic model is to be evaluated.
-fprintf(fid, 'Y = zeros(%u, 1);\n', 2*(n0+n1+n2));
+% In period t, then in period t+1
+fprintf(fid, 'Y0 = NaN(%u, 1);\n', 3*M_.endo_nbr);
+fprintf(fid, 'Y1 = NaN(%u, 1);\n', 3*M_.endo_nbr);
 for i=1:length(I0) % period t equations, lagged variables.
     if I0(i)
-        fprintf(fid, 'Y(%u) = y(%u);\n', I0(i), i);
+        fprintf(fid, 'Y0(%u) = y(%u);\n', i+(purely_forward_model*M_.endo_nbr), i);
     end
 end
 for i=1:length(I1) % period t equations, current variables.
     if I1(i)
-        fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', I1(i), i, i);
+        fprintf(fid, 'Y0(%u) = y(%u)*g(%u);\n', i+M_.endo_nbr+(purely_forward_model*M_.endo_nbr), i, i);
     end
 end
 for i=1:length(I2) % period t equations, leaded variables.
     if I2(i)
-        fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u);\n', I2(i), i, i, i);
+        fprintf(fid, 'Y0(%u) = y(%u)*g(%u)*g(%u);\n', i+2*M_.endo_nbr, i, i, i);
     end
 end
 for i=1:length(I0) % period t+1 equations lagged variables.
     if I0(i)
-        fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', n0+n1+n2+I0(i), i, i);
+        fprintf(fid, 'Y1(%u) = y(%u)*g(%u);\n', i+(purely_forward_model*M_.endo_nbr), i, i);
     end
 end
 for i=1:length(I1) % period t+1 equations current variables.
     if I1(i)
-        fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u);\n', n0+n1+n2+I1(i), i, i, i);
+        fprintf(fid, 'Y1(%u) = y(%u)*g(%u)*g(%u);\n', i+M_.endo_nbr+(purely_forward_model*M_.endo_nbr), i, i, i);
     end
 end
 for i=1:length(I2) % period t+1 equations leaded variables.
     if I2(i)
-        fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u)*g(%u);\n', n0+n1+n2+I2(i), i, i, i, i);
+        fprintf(fid, 'Y1(%u) = y(%u)*g(%u)*g(%u)*g(%u);\n', i+2*M_.endo_nbr, i, i, i, i);
     end
 end
+
 fprintf(fid, '\n');
 
 % Define the vector of parameters.
@@ -130,20 +133,37 @@ fprintf(fid, 'F = NaN(%u, 1);\n', 2*M_.endo_nbr);
 % Set vector of exogenous variables to 0.
 fprintf(fid, 'x = zeros(1, %u);\n\n', M_.exo_nbr);
 
-% Evaluate the residuals and jacobian of the dynamic model in periods t and t+1.
+% Evaluate the residuals and Jacobian of the dynamic model in periods t and t+1.
+fprintf(fid, '[F(1:%u), T0_order, T0] = %s.sparse.dynamic_resid(Y0, x, p, y);\n', M_.endo_nbr, M_.fname);
+fprintf(fid, '[F(%u:%u), T1_order, T1] = %s.sparse.dynamic_resid(Y1, x, p, y);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname);
 fprintf(fid, 'if nargout>1\n');
+fprintf(fid, '    sparse_rowval = [');
+fprintf(fid, '%u ', M_.dynamic_g1_sparse_rowval);
+fprintf(fid, '];\n');
+fprintf(fid, '    sparse_colval = [');
+fprintf(fid, '%u ', M_.dynamic_g1_sparse_colval);
+fprintf(fid, '];\n');
+fprintf(fid, '    sparse_colptr = [');
+fprintf(fid, '%u ', M_.dynamic_g1_sparse_colptr);
+fprintf(fid, '];\n');
+fprintf(fid, '    J0 = %s.sparse.dynamic_g1(Y0, x, p, y, sparse_rowval, sparse_colval, sparse_colptr, T0_order, T0);\n', M_.fname);
+fprintf(fid, '    J1 = %s.sparse.dynamic_g1(Y1, x, p, y, sparse_rowval, sparse_colval, sparse_colptr, T1_order, T1);\n', M_.fname);
+
+% Transform back the Jacobians J0 and J1 in the legacy format (non-sparse)
+% NB: it is probably possible to simplify the rest of this file, but maintaining
+% decent performance does not seem straightforward.
+lli = find(M_.lead_lag_incidence');
+if purely_forward_model
+    lli = lli + M_.endo_nbr;
+end
+fprintf(fid, '    lli = [');
+fprintf(fid, '%u ', lli);
+fprintf(fid, '];\n');
+
 fprintf(fid, '    J = zeros(%u, %u);\n', 2*M_.endo_nbr, n0+n1+n2+M_.endo_nbr);
-fprintf(fid, '    [F(1:%u), tmp] = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
-fprintf(fid, '    J(1:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr, n0+n1+n2, n0+n1+n2);
-fprintf(fid, '    [F(%u:%u), tmp] = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
-fprintf(fid, '    J(%u:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2, n0+n1+n2);
-fprintf(fid, 'else\n');
-fprintf(fid, '    F(1:%u) = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
-fprintf(fid, '    F(%u:%u) = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
-fprintf(fid, 'end\n\n');
+fprintf(fid, '    J(1:%u,1:%u) = full(J0(:,lli));\n', M_.endo_nbr, n0+n1+n2);
+fprintf(fid, '    J(%u:%u,1:%u) = full(J1(:,lli));\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2);
 
-% Compute the jacobian if required.
-fprintf(fid, 'if nargout>1\n');
 fprintf(fid, '    JAC = zeros(%u,%u);\n', 2*M_.endo_nbr, 2*M_.endo_nbr);
 
 % Compute the derivatives of the first block of equations (period t)
@@ -279,7 +299,7 @@ else
 end
 
 % Compute the derivatives of the second block of equations (period t+1)
-% with respect to the endogenous variables.
+% with respect to the growth factors.
 if purely_backward_model || purely_forward_model
     for i=M_.eq_nbr+1:2*M_.eq_nbr
         for j=1:M_.endo_nbr

matlab/+bvar/density.m

 function density(maxnlags)
 % function density(maxnlags)
-% computes the density of a bayesian var
+% computes the density of a Bayesian VAR
 %
 % INPUTS
-%    maxnlags:      maximum number of lags in the bvar
+%    maxnlags:      maximum number of lags in the BVAR
 %
 % OUTPUTS
 %    none

matlab/+bvar/forecast.m

 function forecast(nlags)
 % function forecast(nlags)
-% builds forecats for a bvar model
+% builds forecasts for a BVAR model
 %
 % INPUTS
-%    nlags:     number of lags for the bvar
+%    nlags:     number of lags for the BVAR
 %
 % OUTPUTS
 %    none
@@ -41,7 +41,7 @@ sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
 S_inv_upper_chol = chol(inv(posterior.S));
 
 % Option 'lower' of chol() not available in old versions of
-% Matlab, so using transpose
+% MATLAB, so using transpose
 XXi_lower_chol = chol(posterior.XXi)';
 
 k = ny*nlags+nx;
@@ -58,7 +58,7 @@ while d <= options_.bvar_replic
     Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
 
     % Option 'lower' of chol() not available in old versions of
-    % Matlab, so using transpose
+    % MATLAB, so using transpose
     Sigma_lower_chol = chol(Sigma)';
 
     Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);

matlab/+bvar/irf.m

@@ -3,7 +3,7 @@ function irf(nlags,identification)
 % builds IRFs for a BVAR model
 %
 % INPUTS
-%    nlags            [integer]     number of lags for the bvar
+%    nlags            [integer]     number of lags for the BVAR
 %    identification   [string]      identification scheme ('Cholesky' or 'SquareRoot')
 %
 % OUTPUTS
@@ -40,7 +40,7 @@ end
 S_inv_upper_chol = chol(inv(posterior.S));
 
 % Option 'lower' of chol() not available in old versions of
-% Matlab, so using transpose
+% MATLAB, so using transpose
 XXi_lower_chol = chol(posterior.XXi)';
 
 k = ny*nlags+nx;

matlab/+bvar/toolbox.m

@@ -143,7 +143,7 @@ if prior.df < ny
     error('Too few degrees of freedom in the inverse-Wishart part of prior distribution. You should increase training sample size.')
 end
 
-% Add forecast informations
+% Add forecast information
 if nargout >= 5
     forecast_data.xdata = ones(options_.forecast, nx);
     forecast_data.initval = ydata(end-nlags+1:end, :);

matlab/+estimate/nls.m

@@ -332,7 +332,7 @@ C = C/T;
 % Save results
 dlhs = eval(strrep(lhs, 'data', 'data(range(1)-real(islaggedvariables):range(end)).data')); % REMARK: If lagged variables are present in the estimated equation (ie rhs) then an observation
                                                                                             % has been appended to the dataset (period range(1)-1) and the left hand side has been replaced
-                                                                                            % by data(2:end,lhs_id) in the matlab routine generated to evaluate the residuals (it is also the value
+                                                                                            % by data(2:end,lhs_id) in the MATLAB routine generated to evaluate the residuals (it is also the value
                                                                                             % hold by string variable lhs). Hence to evaluate the left hand side variable used for estimation,
                                                                                             % between range(1) and range(end) we need here to append the same observation in period range(1)-1.
 oo_.nls.(eqname).lhs = dseries(dlhs, range(1), sprintf('%s_lhs', eqname));

matlab/+gsa/map_calibration.m

@@ -105,7 +105,7 @@ if init
     end
 
     irestrictions = 1:Nsam;
-    h = dyn_waitbar(0,'Please wait...');
+    h = waitbar.run(0,'Please wait...');
     for j=1:Nsam
         M_ = set_all_parameters(lpmat(j,:)',estim_params_,M_);
         if nbr_moment_restrictions
@@ -131,10 +131,10 @@ if init
             irestrictions(j)=0;
         end
         if mod(j,3)==0
-            dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
+            waitbar.run(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
         end
     end
-    dyn_waitbar_close(h);
+    waitbar.close(h);
 
     irestrictions=irestrictions(find(irestrictions));
     xmat=lpmat(irestrictions,:);

matlab/+gsa/map_identification.m

@@ -3,9 +3,9 @@ function map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_pa
 % Inputs
 %  - OutputDirectoryName [string]    name of the output directory
 %  - opt_gsa             [structure]     GSA options structure
-%  - M_                  [structure]     Matlab's structure describing the model
-%  - oo_                 [structure]     Matlab's structure describing the results
-%  - options_            [structure]     Matlab's structure describing the current options
+%  - M_                  [structure]     MATLAB's structure describing the model
+%  - oo_                 [structure]     MATLAB's structure describing the results
+%  - options_            [structure]     MATLAB's structure describing the current options
 %  - estim_params_       [structure]     characterizing parameters to be estimated
 %  - bayestopt_          [structure]     describing the priors
 

matlab/+gsa/monte_carlo_filtering.m

@@ -4,10 +4,10 @@ function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dat
 %  - OutputDirectoryName [string]       name of the output directory
 %  - options_gsa_        [structure]    GSA options
 %  - dataset_            [dseries]      object storing the dataset
-%  - dataset_info        [structure]    storing informations about the sample.
-%  - M_                  [structure]    Matlab's structure describing the model
+%  - dataset_info        [structure]    storing information about the sample.
+%  - M_                  [structure]    MATLAB's structure describing the model
 %  - oo_                 [structure]    storing the results
-%  - options_            [structure]    Matlab's structure describing the current options
+%  - options_            [structure]    MATLAB's structure describing the current options
 %  - bayestopt_          [structure]    describing the priors
 %  - estim_params_       [structure]    characterizing parameters to be estimated
 %
@@ -17,7 +17,7 @@ function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dat
 %                                       indices (descending order of RMSEs)
 %
 % Notes: the R^2 definition is 1-var(ymodel-ydata)/var(ydata). It ranges
-% bewteen (-inf, 1], with negative values indicating that themodel is a worse 
+% between (-inf, 1], with negative values indicating that the model is a worse 
 % predictor than the sample mean of the data
 
 % inputs (from opt_gsa structure)

matlab/+gsa/monte_carlo_moments.m

@@ -4,7 +4,7 @@ function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_pa
 % Inputs:
 %  - dr                 [structure]     decision rules
 %  - M_                 [structure]     model structure
-%  - options_           [structure]     Matlab's structure describing the current options
+%  - options_           [structure]     MATLAB's structure describing the current options
 %  - estim_params_      [structure]     characterizing parameters to be estimated
 %
 % Outputs:
@@ -33,7 +33,7 @@ function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_pa
 [~, nc1, nsam] = size(mm);
 nobs=length(options_.varobs);
 disp('monte_carlo_moments: Computing theoretical moments ...')
-h = dyn_waitbar(0,'Theoretical moments ...');
+h = waitbar.run(0,'Theoretical moments ...');
 vdec = zeros(nobs,M_.exo_nbr,nsam);
 cc = zeros(nobs,nobs,nsam);
 ac = zeros(nobs,nobs*options_.ar,nsam);
@@ -52,9 +52,9 @@ for j=1:nsam
     end
     ac(:,:,j)=dum;
     if mod(j,3)==0
-        dyn_waitbar(j/nsam,h)
+        waitbar.run(j/nsam,h)
     end
 end
-dyn_waitbar_close(h)
+waitbar.close(h)
 skipline()
 disp('... done !')

matlab/+gsa/run.m

@@ -2,12 +2,12 @@ function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
 % x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
 % Frontend to the Sensitivity Analysis Toolbox for DYNARE
 % Inputs:
-%  - M_                     [structure]     Matlab's structure describing the model
-%  - oo_                    [structure]     Matlab's structure describing the results
-%  - options_               [structure]     Matlab's structure describing the current options
+%  - M_                     [structure]     MATLAB's structure describing the model
+%  - oo_                    [structure]     MATLAB's structure describing the results
+%  - options_               [structure]     MATLAB's structure describing the current options
 %  - bayestopt_             [structure]     describing the priors
 %  - estim_params_          [structure]     characterizing parameters to be estimated
-%  - options_gsa            [structure]     Matlab's structure describing the GSA options
+%  - options_gsa            [structure]     MATLAB's structure describing the GSA options
 %
 % Reference:
 % M. Ratto (2008), Analysing DSGE Models with Global Sensitivity Analysis, 
@@ -360,9 +360,9 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
         if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
             fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
             fprintf('To obtain it, go to:\n\n')
-            fprintf('https://ec.europa.eu/jrc/en/macro-econometric-statistical-software/ss-anova-r-downloads \n\n')
+            fprintf('https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip \n\n')
             fprintf('and follow the instructions there.\n')
-            fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
+            fprintf('After obtaining the files, you need to unpack them and set a MATLAB Path to those files.\n')
             error('SS-ANOVA-R Toolbox missing!')
         end
         gsa.reduced_form_mapping(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);

matlab/+gsa/set_shocks_param.m

@@ -2,11 +2,11 @@ function M_=set_shocks_param(M_,estim_params_,xparam1)
 % function M_=set_shocks_param(M_,estim_params_,xparam1)
 % Set the structural and measurement error variances and covariances
 % Inputs
-%  - M_                     [structure]     Matlab's structure describing the model
+%  - M_                     [structure]     MATLAB's structure describing the model
 %  - estim_params_          [structure]     characterizing parameters to be estimated
 %  - xparam1                [double]        parameter vector
 % Outputs:
-%  - M_                     [structure]     Matlab's structure describing the model
+%  - M_                     [structure]     MATLAB's structure describing the model
 %
 % Notes: closely follows set_all_parameters.m
 

matlab/+gsa/stability_mapping.m

@@ -6,9 +6,9 @@ function x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,baye
 % Inputs
 %  - OutputDirectoryName    [string]        name of the output directory
 %  - opt_gsa                [structure]     GSA options structure
-%  - M_                     [structure]     Matlab's structure describing the model
-%  - oo_                    [structure]     Matlab's structure describing the results
-%  - options_               [structure]     Matlab's structure describing the current options
+%  - M_                     [structure]     MATLAB's structure describing the model
+%  - oo_                    [structure]     MATLAB's structure describing the results
+%  - options_               [structure]     MATLAB's structure describing the current options
 %  - bayestopt_             [structure]     describing the priors
 %  - estim_params_          [structure]     characterizing parameters to be estimated
 %
@@ -18,7 +18,7 @@ function x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,baye
 %
 % Inputs from opt_gsa structure
 % Nsam = MC sample size
-% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
+% fload = 0 to run new MC; 1 to load previously generated analysis
 % alpha2 =  significance level for bivariate sensitivity analysis
 % [abs(corrcoef) > alpha2]
 % prepSA = 1: save transition matrices for mapping reduced form
@@ -30,7 +30,7 @@ function x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,baye
 %
 % GRAPHS
 % 1) Pdf's of marginal distributions under the stability (dotted
-%     lines) and unstability (solid lines) regions
+%     lines) and instability (solid lines) regions
 % 2) Cumulative distributions of:
 %   - stable subset (dotted lines)
 %   - unacceptable subset (solid lines)
@@ -92,6 +92,10 @@ xparam1=[];
 [~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
 if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
     % Set prior bounds
+    if options_.prior_trunc==0
+        fprintf('\nstability_mapping: GSA with priors requires bounded support. Setting options_.prior_trunc=1e-10.\n')
+        options_.prior_trunc=1e-10;
+    end
     bounds = prior_bounds(bayestopt_, options_.prior_trunc);
     bounds.lb = max(bounds.lb,lb);
     bounds.ub = min(bounds.ub,ub);
@@ -232,7 +236,7 @@ if fload==0 %run new MC
         end
     end
     %
-    h = dyn_waitbar(0,'Please wait...');
+    h = waitbar.run(0,'Please wait...');
     istable=1:Nsam;
     jstab=0;
     iunstable=1:Nsam;
@@ -329,10 +333,10 @@ if fload==0 %run new MC
         ys_=real(dr_.ys);
         yys(:,j) = ys_;
         if mod(j,3)
-            dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
+            waitbar.run(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
         end
     end
-    dyn_waitbar_close(h);
+    waitbar.close(h);
     if prepSA && jstab
         T=T(:,:,1:jstab);
     else
@@ -343,7 +347,7 @@ if fload==0 %run new MC
     inorestriction=inorestriction(inorestriction~=0);  % stable params violating restrictions
     iunstable=iunstable(iunstable~=0);   % violation of BK & restrictions & solution could not be found (whatever goes wrong)
     iindeterm=iindeterm(iindeterm~=0);  % indeterminacy
-    iwrong=iwrong(iwrong~=0);  % dynare could not find solution
+    iwrong=iwrong(iwrong~=0);  % Dynare could not find solution
     ixun=iunstable(~ismember(iunstable,[iindeterm,iwrong,inorestriction])); % explosive roots
 
     bkpprior.pshape=bayestopt_.pshape;
@@ -385,7 +389,7 @@ else %load old run
     end
 
     if prepSA && isempty(strmatch('T',who('-file', filetoload),'exact'))
-        h = dyn_waitbar(0,'Please wait...');
+        h = waitbar.run(0,'Please wait...');
         options_.periods=0;
         options_.nomoments=1;
         options_.irf=0;
@@ -404,10 +408,10 @@ else %load old run
             ys_=real(dr_.ys);
             yys(:,j) = ys_;
             if mod(j,3)
-                dyn_waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
+                waitbar.run(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
             end
         end
-        dyn_waitbar_close(h);
+        waitbar.close(h);
         save(filetoload,'T','-append')
     end
 end
@@ -447,7 +451,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
         end
         if ~isempty(iwrong)
             skipline()
-            disp(['For ',num2str(length(iwrong)/Nsam*100,'%4.1f'),'% of the prior support dynare could not find a solution.'])
+            disp(['For ',num2str(length(iwrong)/Nsam*100,'%4.1f'),'% of the prior support Dynare could not find a solution.'])
             skipline()
         end
         if any(infox==1)
@@ -457,7 +461,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
             disp(['    For ',num2str(length(find(infox==2))/Nsam*100,'%4.1f'),'% MJDGGES returned an error code.'])
         end
         if any(infox==6)
-            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%4.1f'),'% The jacobian evaluated at the deterministic steady state is complex.'])
+            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%4.1f'),'% The Jacobian evaluated at the deterministic steady state is complex.'])
         end
         if any(infox==19)
             disp(['    For ',num2str(length(find(infox==19))/Nsam*100,'%4.1f'),'% The steadystate routine has thrown an exception (inconsistent deep parameters).'])
@@ -485,7 +489,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
     skipline()
     if length(iunstable)<Nsam || length(istable)>1
         itot = 1:Nsam;
-        isolve = itot(~ismember(itot,iwrong)); % dynare could find a solution
+        isolve = itot(~ismember(itot,iwrong)); % Dynare could find a solution
                                                      % Blanchard Kahn
         if neighborhood_width
             options_mcf.xparam1 = xparam1(nshock+1:end);

matlab/+gsa/stability_mapping_bivariate.m

@@ -4,8 +4,8 @@ function indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options
 %  - x
 %  - alpha2
 %  - pvalue_crit
-%  - M_                     [structure]     Matlab's structure describing the model
-%  - options_               [structure]     Matlab's structure describing the current options
+%  - M_                     [structure]     MATLAB's structure describing the model
+%  - options_               [structure]     MATLAB's structure describing the current options
 %  - bayestopt_             [structure]     describing the priors
 %  - estim_params_          [structure]     characterizing parameters to be estimated
 %  - fnam                   [string]        file name

matlab/+gsa/th_moments.m

@@ -4,8 +4,8 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,options_,M_)
 %
 % INPUTS
 % - dr            [structure]     model information structure
-% - options_      [structure]     Matlab's structure describing the current options
-% - M_            [structure]     Matlab's structure describing the model
+% - options_      [structure]     MATLAB's structure describing the current options
+% - M_            [structure]     MATLAB's structure describing the model
 %
 % OUTPUTS
 % - vdec          [double]        variance decomposition matrix

matlab/+gui/+perfect_foresight/run.m

@@ -13,7 +13,7 @@ function run(json)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2019-2023 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -103,7 +103,6 @@ if ~isempty(jm.anticipated_transitory_shocks)
             'periods', s.start_date:s.end_date, ...
             'value', s.value)];
     end
-    M_.exo_det_length = 0;
 end
 
 %% Make unanticipated shock map

matlab/+gui/+stochastic-simulation/read.m

@@ -13,7 +13,7 @@ function read(json)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2019-2020 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +36,6 @@ global M_ options_ oo_
 jm = loadjson_(json, 'SimplifyCell', 1);
 data2json=struct();
 
-M_.exo_det_length = 0;
 for nshocks = 1:length(jm.stochasticshocksdescription)
 	covartype=jm.stochasticshocksdescription{nshocks}.shockattributevalue;
 	thisshock=(jm.stochasticshocksdescription{nshocks}.shockindex)+1;

matlab/+hank/check_steady_state_input.m

+% Copyright © 2025 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+%
+% Compute the stochastic simulations of heterogeneus-agent models
+%
+% INPUTS
+% - M_       [structure] MATLAB's structure describing the model
+% - options_ [structure] MATLAB's structure describing the current options
+% - ss       [structure] MATLAB's structure with the model steady-state
+%                       information. It contains the following fields:
+%    - pol [structure] MATLAB's structure containing the policy functions
+%                      discretization
+%       - pol.grids [structure]: MATLAB's structure containing the nodes of the
+%                                state grids as column vectors.
+%       - pol.values [structure]: MATLAB's structure containing the policy
+%                                 function as matrices. Row indices and column
+%                                 indices follow the lexicographic order
+%                                 specified in pol.shocks and pol.states
+%       - pol.shocks [array]: a list containing the order of shock/innovation
+%                             variables used for the rows of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement (resp. varexo(heterogeneity=)
+%                             statement) for discretizes AR(1) processes (resp.
+%                             discretized gaussian i.i.d innovations).
+%       - pol.states [array]: a list containing the order of state variables
+%                             used for the columns of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement.
+%    - shocks [structure]: MATLAB's structure describing the discretization of
+%                          individual shocks or innovations:
+%       - shocks.grids [structure]: MATLAB's structure containing the nodes of
+%                                   the shock grids
+%       - shocks.Pi [structure]: MATLAB's structure containing the Markov
+%                                matrices if the shock processes are discretized
+%                                AR(1) processes. The field should be absent
+%                                otherwise.
+%       - shocks.w [structure]: MATLAB's structure containing the Gauss-Hermite
+%                               weights if the i.i.d gaussian innovation
+%                               processes are discretized
+%    - d [structure]: MATLAB's structure describing the steady-state
+%                     distribution
+%       - d.grids [structure]: structure containing the states grids as column
+%                              vectors. If one of the states grid is not
+%                              specified, it is assumed to be identical to
+%                              its policy-function value stored in pol.grids
+%       - d.hist [matrix]: the histogram of the distribution as a matrix with
+%                          shocks/innovation indices as rows and state indices
+%                          as columns. Row and column indices follow the
+%                          lexicographic order specified in d.shocks and
+%                          d.states
+%       - d.shocks [array]: a list containing the order of shock/innovation
+%                           variables used for the rows of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.shocks.
+%       - d.states [array]: a list containing the order of state variables used
+%                           for the columns of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.states.
+%    - agg [structure]: MATLAB's structure containing the steady-state values of
+%                       aggregate variables
+% OUTPUTS
+% - out_ss [structure]: validated and normalized steady-state input structure.
+%                       It has a similar structure as the ss structure in inputs.
+% - sizes [structure]: structure containing sizes information
+%    - n_e [integer]: number of shocks
+%    - N_e [integer]: total number of nodes in the shocks grids
+%    - shocks [structure]: number of nodes for each shock grid
+%    - n_a [integer]: number of states
+%    - pol [structure]: 
+%       - pol.N_a [integer]: total number of nodes in the policy state grids
+%       - pol.states [structure]: number of nodes for each state grid
+%    - d [structure]:
+%       - d.states [structure]: number of nodes for each distribution state grid
+%       - d.N_a [integer]: total number of nodes in the distribution grids
+%    - n_pol [integer]: number of policy variables
+%    - agg [integer]: number of aggregate variables
+function [out_ss, sizes] = check_steady_state_input(M_, options_, ss)
+   %% Checks
+   % Retrieve variables information from M_
+   H_ = M_.heterogeneity(1);
+   state_symbs = H_.endo_names(H_.state_var);
+   inn_symbs = H_.exo_names;
+   agg_symbs = M_.endo_names(1:M_.orig_endo_nbr);
+   required_pol_symbs = H_.endo_names(1:H_.orig_endo_nbr);
+   mult_symbs = {};
+   for i = 1:numel(H_.aux_vars)
+      if (H_.aux_vars(i).type == 15)
+         mult_symbs{end+1} = H_.endo_names{H_.aux_vars(i).endo_index};
+      end
+   end
+   % Initialize output variables
+   sizes = struct;
+   out_ss = struct;
+   if ~isstruct(ss)
+      error('Misspecified steady-state input `ss`: the steady-state input `ss` is not a structure.')
+   end
+   %% Shocks
+   % Check that the field `ss.shocks` exists
+   check_isfield('shocks', ss, 'ss.shocks', ' Models without idiosyncratic shocks but with ex-post heterogeneity over individual state variables cannot be solved yet.')
+   shocks = ss.shocks;
+   % Check that `ss.shocks` is a structure
+   check_isstruct(shocks, 'ss.shocks');
+   % Check that the field `ss.shocks.grids` exists
+   check_isfield('grids', ss.shocks, 'ss.shocks.grids');
+   % Check that `ss.shocks.grids` is a structure
+   check_isstruct(shocks.grids, 'ss.shocks.grids');
+   shock_symbs = fieldnames(shocks.grids);
+   % Check the types of `ss.shocks.grids` values
+   check_fun(shocks.grids, 'ss.shocks.grids', shock_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(shock_symbs)
+      s = shock_symbs{i};
+      if isrow(ss.shocks.grids.(s))
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s)';
+      else
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s);
+      end
+   end
+   % Check shock discretization
+   flag_ar1 = isfield(shocks, 'Pi');
+   flag_gh = isfield(shocks, 'w');
+   if ~flag_ar1 && ~flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both non-specified, while exactly one of them should be set.');
+   end
+   % Check that either individual AR(1) processes or i.i.d gaussion innovation
+   % processes are discretized
+   if flag_ar1 && flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both specified, while only one of them should be.');
+   end
+   % Case of discretized AR(1) exogenous processes
+   if flag_ar1
+      % Check that the grids specification for AR(1) shocks and the
+      % var(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, required_pol_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).endo_names');
+      % Check that shocks.Pi is a structure
+      check_isstruct(shocks.Pi, 'ss.shocks.Pi');
+      % Check that the grids specification for individual AR(1) shocks, the
+      % information in the M_ structure and the Markov transition matrices
+      % specification are compatible and display a warning if redundancies are
+      % detected
+      check_missingredundant(shocks.Pi, 'ss.shocks.Pi', shock_symbs, options_.hank.nowarningredundant);
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the type of shocks.Pi field values
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) ismatrix(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real matrices');
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      % Check that the matrices provided in shocks.Pi are Markov matrices:
+      %     - all elements are non-negative
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(x >= 0, 'all'), 'contain negative elements')
+      %     - the sum of each row equals 1.
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(abs(sum(x,2)-1) < options_.hank.tol_check_sum), 'have row sums different from 1.')
+      % Remove AR(1) shock processes from state and policy variables
+      state_symbs = setdiff(state_symbs, shock_symbs);
+      required_pol_symbs = setdiff(required_pol_symbs, shock_symbs);
+      % Copy relevant shock-related elements in out_ss 
+      out_ss.shocks.Pi = ss.shocks.Pi;
+   end
+   relevant_pol_symbs = [required_pol_symbs; mult_symbs];
+   % Case of discretized i.i.d gaussian innovations
+   if flag_gh
+      % Check that shocks.w is a structure
+      check_isstruct(shocks.w, 'shocks.w');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the varexo(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, inn_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).exo_names');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the Gauss-Hermite weights specification are compatible and display a
+      % warning if redundancies are detected
+      check_missingredundant(shocks.w, 'ss.shocks.w', shock_symbs, options_.hank.nowarningredundant);
+      % Check the type of `ss.shocks.w` elements
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) isvector(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real vectors');
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids 
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) numel(x), 'have incompatible sizes with those of `ss.shocks.grids`', grid_nb);
+      % Check that the arrays provided in shocks.w have the properties of
+      % Gauss-Hermite weights:
+      %     - all elements are non-negative
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(x >= 0.), 'contain negative elements');
+      %     - the sum of Gauss-Hermite weights is 1
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(abs(sum(x)-1) < options_.hank.tol_check_sum), 'have sums different from 1.');
+      % Make `ss.shocks.w` values column vectors in the output steady-state
+      % structure
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         if isrow(ss.shocks.w.(s))
+            out_ss.shocks.w.(s) = ss.shocks.w.(s)';
+         else
+            out_ss.shocks.w.(s) = ss.shocks.w.(s);
+         end
+      end
+   end
+   %% Policy functions
+   % Check that the `ss.pol` field exists
+   check_isfield('pol', ss, 'ss.pol');
+   pol = ss.pol;
+   % Check that `ss.pol` is a structure
+   check_isstruct(ss.pol, 'ss.pol');
+   % Check that the `ss.pol`.grids field exists
+   check_isfield('grids', ss.pol, 'ss.pol.grids');
+   % Check that `ss.pol.grids` is a structure
+   check_isstruct(ss.pol.grids, 'ss.pol.grids');
+   % Check that the state grids specification and the var(heterogeneity=)
+   % statement are compatible
+   check_missingredundant(pol.grids, 'ss.pol.grids', state_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.grids` values are dense real vectors
+   check_fun(pol.grids, 'ss.pol.grids', state_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Store:
+   %    - the number of states
+   sizes.n_a = numel(state_symbs);
+   %    - the size of the tensor product of the states grids 
+   grid_nb = cellfun(@(s) numel(pol.grids.(s)), state_symbs);
+   sizes.pol.N_a = prod(grid_nb);
+   %    - the size of each state grid
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      sizes.pol.states.(s) = numel(pol.grids.(s));
+   end
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      if isrow(ss.pol.grids.(s))
+         out_ss.pol.grids.(s) = ss.pol.grids.(s)';
+      else
+         out_ss.pol.grids.(s) = ss.pol.grids.(s);
+      end
+   end
+   % Check that the field `ss.pol.values` exists 
+   check_isfield('values', pol, 'ss.pol.values');
+   % Check that `ss.pol.values` is a struct
+   check_isstruct(pol.values, 'ss.pol.values');
+   % Check the missing and redundant variables in `ss.pol.values`
+   check_missingredundant(pol.values, 'ss.pol.values', required_pol_symbs, true);
+   pol_values = fieldnames(pol.values);
+   pol_values_in_relevant_pol_symbs = ismember(pol_values, relevant_pol_symbs);
+   if ~options_.hank.nowarningredundant
+      if ~all(pol_values_in_relevant_pol_symbs)
+         warning('Steady-state input `ss`. The following fields are redundant in `%s`: %s.', 'ss.pol.values', strjoin(pol_values(~pol_values_in_relevant_pol_symbs)));
+      end
+   end
+   provided_relevant_pol_symbs = pol_values(pol_values_in_relevant_pol_symbs);
+   % Check that `ss.pol.values` values are dense real matrices
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) isnumeric(x) && isreal(x) && ismatrix(x) && ~issparse(x), 'are not dense real matrices');
+   % Check the internal size compatibility of `ss.pol.values`
+   sizes.n_pol = H_.endo_nbr;
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the sizes of `ss.shocks.grids` elements', sizes.N_e);
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the sizes of `ss.pol.grids` elements', sizes.pol.N_a);
+   % Copy `ss.pol.values` in `out_ss`
+   out_ss.pol.values = ss.pol.values;
+   % Check the permutation of state variables for policy functions
+   out_ss.pol.states = check_permutation(pol, 'states', 'ss.pol', state_symbs);
+   % Check the permutation of shock variables for policy functions
+   out_ss.pol.shocks = check_permutation(pol, 'shocks', 'ss.pol', shock_symbs);
+   %% Distribution
+   % Check that the field `ss.d` exists 
+   check_isfield('d', ss, 'ss.d');
+   d = ss.d;
+   % Check that the field `ss.d.hist` exists 
+   check_isfield('hist', ss.d, 'ss.d.hist');
+   % Check the type of `ss.d.hist`
+   if ~(ismatrix(d.hist) && isnumeric(d.hist) && isreal(d.hist) && ~issparse(d.hist))
+      error('Misspecified steady-state input `ss`: `ss.d.hist` is not a dense real matrix.');
+   end
+   % Check the consistency of `ss.d.grids` 
+   if ~isfield(d, 'grids')
+      if ~options_.hank.nowarningdgrids
+         warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+      end
+      % Copy the relevant sizes from the pol field
+      sizes.d = sizes.pol;
+      out_ss.d.grids = out_ss.pol.grids;
+   else
+      % Check that `ss.d.grids` is a struct
+      check_isstruct(d.grids, 'ss.d.grids');
+      % Check redundant variables in `ss.d.grids`
+      d_grids_symbs = fieldnames(d.grids);
+      if isempty(d_grids_symbs)
+         if ~options_.hank.nowarningredundant
+            warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+         end
+         % Copy the relevant sizes from the pol field
+         sizes.d = sizes.pol;
+         out_ss.d.grids = out_ss.pol.grids;
+      else
+         d_grids_symbs_in_states = ismember(d_grids_symbs, state_symbs);
+         if ~all(d_grids_symbs_in_states)
+            if ~options_.hank.nowarningredundant
+               warning('In the steady-state input `ss.d.states`, the following specification for the states grids in the distribution structure are not useful: %s.', strjoin(d_grids_symbs(~d_grids_symbs_in_states)));
+            end
+            d_grids_symbs = d_grids_symbs(d_grids_symbs_in_states);
+         end
+         % Check the types of `ss.d.grids` elements
+         check_fun(pol.grids, 'ss.d.grids', d_grids_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+         % Store the size of the states grids for the distribution
+         for i=1:numel(d_grids_symbs)
+            s = d_grids_symbs{i};
+            sizes.d.states.(s) = numel(d.grids.(s));
+            out_ss.d.grids.(s) = d.grids.(s);
+         end
+         states_out_of_d = setdiff(state_symbs, d_grids_symbs);
+         if ~isempty(states_out_of_d)
+            if ~options_.hank.nowarningdgrids
+               warning('hank.bbeg.het_stoch_simul: in the steady-state structure, no distribution-specific grid set for states %s. The policy grids specified in pol.grids shall be used.', strjoin(states_out_of_d));
+            end
+            % Store the sizes of the unspecified states grids from the pol field
+            for i=1:numel(states_out_of_d)
+               s = states_out_of_d{i};
+               sizes.d.states.(s) = sizes.pol.states.(s);
+               out_ss.d.grids.(s) = pol.grids.(s);
+            end
+         end
+      end
+   end
+   % Check the internal size compatibility of the distribution histogram
+   if size(d.hist,1) ~= sizes.N_e
+      error('Misspecified steady-state input `ss`: the number of rows of the histogram matrix `ss.d.hist` is not consistent with the sizes of shocks grids `ss.shocks.grids`.');
+   end
+   sizes.d.N_a = prod(structfun(@(x) x, sizes.d.states));
+   if size(d.hist,2) ~= sizes.d.N_a
+      error('Misspecified steady-state input `ss`: the number of columns of the histogram matrix `ss.d.hist` is not consistent with the sizes of states grids `ss.d.grids`/`ss.pol.grids`.');
+   end
+   % Copy `ss.d.hist` in `out_ss` 
+   out_ss.d.hist = ss.d.hist;
+   % Check the permutation of state variables in the distribution
+   out_ss.d.states = check_permutation(d, 'states', 'ss.d', state_symbs);
+   % Check the permutation of shock variables in the distribution
+   out_ss.d.shocks = check_permutation(d, 'shocks', 'ss.d', shock_symbs);
+   %% Aggregate variables
+   % Check that the field `ss.agg` exists
+   check_isfield('agg', ss, 'ss.agg');
+   agg = ss.agg;
+   % Check that `ss.agg` is a structure
+   check_isstruct(agg, 'ss.agg');
+   % Check that the aggregate variables specification and the var statement are
+   % compatible
+   check_missingredundant(agg, 'ss.agg', agg_symbs, options_.hank.nowarningredundant);
+   % Store the number of aggregate variables
+   sizes.agg = numel(fieldnames(agg));
+   % Check the types of `ss.agg` values
+   check_fun(agg, 'ss.agg', agg_symbs, @(x) isreal(x) && isscalar(x), 'are not real scalars');
+   % Copy `ss.agg` into `out_ss`
+   out_ss.agg = agg;
+end
\ No newline at end of file