src/HeterogeneousModel.hh

+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef HETEROGENEOUS_MODEL_HH
+#define HETEROGENEOUS_MODEL_HH
+
+#include <string>
+
+#include "ModelTree.hh"
+
+using namespace std;
+
+class HeterogeneousModel : public ModelTree
+{
+public:
+  const int heterogeneity_dimension;
+
+  HeterogeneousModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
+                     ExternalFunctionsTable& external_functions_table_arg,
+                     HeterogeneityTable& heterogeneity_table_arg, int heterogeneity_dimension_arg);
+
+  HeterogeneousModel(const HeterogeneousModel& m);
+  HeterogeneousModel& operator=(const HeterogeneousModel& m);
+
+  void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);
+
+  void writeModelFiles(const string& basename, bool julia) const;
+  void writeDriverOutput(ostream& output) const;
+
+  [[nodiscard]] int getJacobianCol(int deriv_id, bool sparse) const override;
+  [[nodiscard]] int getJacobianColsNbr(bool sparse) const override;
+
+#if 0
+  void substituteEndoLeadGreaterThanTwo();
+
+  //! Transforms the model by removing all lags greater or equal than 2 on endos
+  void substituteEndoLagGreaterThanTwo();
+
+  //! Transforms the model by removing all leads on exos
+  /*! Note that this can create new lags on endos and exos */
+  void substituteExoLead();
+
+  //! Transforms the model by removing all lags on exos
+  void substituteExoLag();
+
+  //! Transforms the model by removing all UnaryOpcode::expectation
+  void substituteExpectation(bool partial_information_model);
+
+  //! Transforms the model by decreasing the lead/lag of predetermined variables in model equations
+  //! by one
+  void transformPredeterminedVariables();
+
+  //! Substitutes out all model-local variables
+  void substituteModelLocalVariables();
+#endif
+
+  // FIXME: the following 5 functions are identical to those in DynamicModel. Factorization?
+  [[nodiscard]] int getDerivID(int symb_id, int lead_lag) const noexcept(false) override;
+  [[nodiscard]] SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getLagByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getTypeSpecificIDByDerivID(int deriv_id) const override;
+
+protected:
+  void computeChainRuleJacobian() override;
+  int getBlockJacobianEndoCol(int blk, int var, int lead_lag) const override;
+  string
+  modelClassName() const override
+  {
+    return "dynamic model for heterogeneity dimension '"
+           + heterogeneity_table.getName(heterogeneity_dimension) + "'";
+  }
+  int getMFS() const override;
+
+private:
+  // Maps a pair (symbol ID, lead/lag) to a deriv ID
+  map<pair<int, int>, int> deriv_id_table;
+  // Maps a deriv ID to a pair (symbol ID, lead/lag)
+  vector<pair<int, int>> inv_deriv_id_table;
+
+  // Allocates the derivation IDs for all endogenous variables for this heterogeneity dimension
+  void computeDerivIDs();
+};
+
+#endif

src/MacroExpandModFile.cc

@@ -41,7 +41,7 @@ macroExpandModFile(const filesystem::path& filename, const istream& modfile, boo
   if (save_macro)
     {
       if (save_macro_file.empty())
-        save_macro_file = filename.stem().string() + "-macroexp.mod";
+        save_macro_file = filename.stem().string() + "_macroexp.mod";
       ofstream macro_output_file {save_macro_file};
       if (macro_output_file.fail())
         {

src/ModFile.cc

@@ -31,25 +31,48 @@
 #include "Shocks.hh"
 
 ModFile::ModFile(WarningConsolidation& warnings_arg) :
+    symbol_table {heterogeneity_table},
     var_model_table {symbol_table},
     trend_component_model_table {symbol_table},
     var_expectation_model_table {symbol_table},
     pac_model_table {symbol_table},
-    expressions_tree {symbol_table, num_constants, external_functions_table},
-    original_model {symbol_table, num_constants, external_functions_table,
-                    trend_component_model_table, var_model_table},
-    dynamic_model {symbol_table, num_constants, external_functions_table,
-                   trend_component_model_table, var_model_table},
-    trend_dynamic_model {symbol_table, num_constants, external_functions_table,
-                         trend_component_model_table, var_model_table},
-    orig_ramsey_dynamic_model {symbol_table, num_constants, external_functions_table,
-                               trend_component_model_table, var_model_table},
-    epilogue {symbol_table, num_constants, external_functions_table, trend_component_model_table,
+    expressions_tree {symbol_table, num_constants, external_functions_table, heterogeneity_table},
+    original_model {symbol_table,
+                    num_constants,
+                    external_functions_table,
+                    heterogeneity_table,
+                    trend_component_model_table,
+                    var_model_table},
+    dynamic_model {symbol_table,
+                   num_constants,
+                   external_functions_table,
+                   heterogeneity_table,
+                   trend_component_model_table,
+                   var_model_table},
+    trend_dynamic_model {symbol_table,
+                         num_constants,
+                         external_functions_table,
+                         heterogeneity_table,
+                         trend_component_model_table,
+                         var_model_table},
+    orig_ramsey_dynamic_model {symbol_table,
+                               num_constants,
+                               external_functions_table,
+                               heterogeneity_table,
+                               trend_component_model_table,
+                               var_model_table},
+    epilogue {symbol_table,
+              num_constants,
+              external_functions_table,
+              heterogeneity_table,
+              trend_component_model_table,
               var_model_table},
-    static_model {symbol_table, num_constants, external_functions_table},
-    steady_state_model {symbol_table, num_constants, external_functions_table, static_model},
+    static_model {symbol_table, num_constants, external_functions_table, heterogeneity_table},
+    steady_state_model {symbol_table, num_constants, external_functions_table, heterogeneity_table,
+                        static_model},
     warnings {warnings_arg}
 {
+  heterogeneity_table.setSymbolTable(&symbol_table);
 }
 
 void
@@ -355,11 +378,9 @@ ModFile::checkPass(bool nostrict, bool stochastic)
   // Test if some estimated parameters are used within the values of shocks
   // statements (see issue #469)
   set<int> parameters_intersect;
-  set_intersection(mod_file_struct.parameters_within_shocks_values.begin(),
-                   mod_file_struct.parameters_within_shocks_values.end(),
-                   mod_file_struct.estimated_parameters.begin(),
-                   mod_file_struct.estimated_parameters.end(),
-                   inserter(parameters_intersect, parameters_intersect.begin()));
+  ranges::set_intersection(mod_file_struct.parameters_within_shocks_values,
+                           mod_file_struct.estimated_parameters,
+                           inserter(parameters_intersect, parameters_intersect.begin()));
   if (parameters_intersect.size() > 0)
     {
       cerr << "ERROR: some estimated parameters (";
@@ -378,8 +399,8 @@ ModFile::checkPass(bool nostrict, bool stochastic)
   // Check if some exogenous is not used in the model block, Issue #841
   set<int> unusedExo0 = dynamic_model.findUnusedExogenous();
   set<int> unusedExo;
-  set_difference(unusedExo0.begin(), unusedExo0.end(), mod_file_struct.pac_params.begin(),
-                 mod_file_struct.pac_params.end(), inserter(unusedExo, unusedExo.begin()));
+  ranges::set_difference(unusedExo0, mod_file_struct.pac_params,
+                         inserter(unusedExo, unusedExo.begin()));
   if (unusedExo.size() > 0)
     {
       ostringstream unused_exos;
@@ -398,6 +419,120 @@ ModFile::checkPass(bool nostrict, bool stochastic)
           exit(EXIT_FAILURE);
         }
     }
+
+  if (!heterogeneity_table.empty())
+    {
+      if (block)
+        {
+          cerr << "ERROR: the 'block' option of the 'model' block is not supported for "
+                  "heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.check_present)
+        {
+          cerr << "ERROR: The 'check' command is not supported for heterogeneous models" << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.steady_present)
+        {
+          cerr << "ERROR: The 'steady' command is not supported for heterogeneous models" << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.perfect_foresight_solver_present)
+        {
+          cerr << "ERROR: The 'perfect_foresight_solver' command is not supported for "
+                  "heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.perfect_foresight_with_expectation_errors_solver_present)
+        {
+          cerr << "ERROR: The 'perfect_foresight_with_expectation_errors_solver' command is not "
+                  "supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.stoch_simul_present)
+        {
+          cerr << "ERROR: The 'stoch_simul' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.estimation_present)
+        {
+          cerr << "ERROR: The 'estimation' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.osr_present)
+        {
+          cerr << "ERROR: The 'osr' command is not supported for heterogeneous models" << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.osr_params_present)
+        {
+          cerr << "ERROR: The 'osr_params' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.optim_weights_present)
+        {
+          cerr << "ERROR: The 'optim_weights' block is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.ramsey_model_present)
+        {
+          cerr << "ERROR: The 'ramsey_model' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.discretionary_policy_present)
+        {
+          cerr << "ERROR: The 'discretionary_policy' command is not supported for heterogeneous "
+                  "models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.planner_objective_present)
+        {
+          cerr << "ERROR: The 'planner_objective' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.extended_path_present)
+        {
+          cerr << "ERROR: The 'extended_path' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.identification_present)
+        {
+          cerr << "ERROR: The 'identification' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.sensitivity_present)
+        {
+          cerr << "ERROR: The 'sensitivity' command is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.mom_estimation_present)
+        {
+          cerr
+              << "ERROR: The 'methods_of_moments' command is not supported for heterogeneous models"
+              << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (mod_file_struct.occbin_constraints_present)
+        {
+          cerr << "ERROR: The 'occbin_constraints' block is not supported for heterogeneous models"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
 }
 
 void
@@ -531,9 +666,12 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
       */
       if (linear)
         orig_ramsey_dynamic_model = dynamic_model;
-      DynamicModel ramsey_FOC_equations_dynamic_model {
-          symbol_table, num_constants, external_functions_table, trend_component_model_table,
-          var_model_table};
+      DynamicModel ramsey_FOC_equations_dynamic_model {symbol_table,
+                                                       num_constants,
+                                                       external_functions_table,
+                                                       heterogeneity_table,
+                                                       trend_component_model_table,
+                                                       var_model_table};
       ramsey_FOC_equations_dynamic_model = dynamic_model;
       auto clone_if_not_null
           = [&](expr_t e) { return e ? e->clone(ramsey_FOC_equations_dynamic_model) : nullptr; };
@@ -552,6 +690,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
   // Must come after detrending of variables and Ramsey policy transformation
   dynamic_model.substituteLogTransform();
 
+  if (!heterogeneity_table.empty())
+    dynamic_model.substituteAggregationOperators();
+
   /* Create auxiliary vars for leads and lags greater than 2, on both endos and
      exos. The transformation is not exactly the same on stochastic and
      deterministic models, because there is no need to take into account the
@@ -704,6 +845,16 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
           exit(EXIT_FAILURE);
         }
     }
+
+  for (int dim {0}; dim < heterogeneity_table.size(); dim++)
+    if (heterogeneous_models.at(dim).equation_number() != symbol_table.het_endo_nbr(dim))
+      {
+        cerr << "ERROR: There are " << heterogeneous_models.at(dim).equation_number()
+             << " equations but " << symbol_table.het_endo_nbr(dim)
+             << " endogenous variables in the model for heterogeneity dimension '"
+             << heterogeneity_table.getName(dim) << "'!" << endl;
+        exit(EXIT_FAILURE);
+      }
 }
 
 void
@@ -824,6 +975,9 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
   // Those matrices can only be filled here, because we use derivatives
   dynamic_model.fillVarModelTableMatrices();
 
+  for (auto& hm : heterogeneous_models)
+    hm.computingPass(mod_file_struct.order_option, no_tmp_terms, use_dll);
+
   for (auto& statement : statements)
     statement->computingPass(mod_file_struct);
 
@@ -866,6 +1020,11 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
 
   auto plusfolder {DataTree::packageDir(basename)};
 
+  if (check_model_changes && !heterogeneity_table.empty())
+    {
+      cerr << "ERROR: the 'fast' option is not supported for heterogeneous models" << endl;
+      exit(EXIT_FAILURE);
+    }
   bool hasModelChanged = !dynamic_model.isChecksumMatching(basename) || !check_model_changes;
   if (hasModelChanged)
     {
@@ -906,9 +1065,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
     mOutputFile << "clearvars -global" << endl
                 << "clear_persistent_variables(fileparts(which('dynare')), false)" << endl;
   else if (clear_global)
-    mOutputFile
-        << "clear M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info estimation_info;"
-        << endl;
+    mOutputFile << "clearvars -global M_ options_ oo_ estim_params_ bayestopt_ dataset_ "
+                   "dataset_info estimation_info;"
+                << endl;
 
   if (!notime)
     mOutputFile << "tic0 = tic;" << endl;
@@ -941,6 +1100,7 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
     mOutputFile << "M_.parameter_used_with_lead_lag = true;" << endl;
 
   symbol_table.writeOutput(mOutputFile);
+  heterogeneity_table.writeOutput(mOutputFile);
 
   var_model_table.writeOutput(basename, mOutputFile);
   trend_component_model_table.writeOutput(basename, mOutputFile);
@@ -952,6 +1112,10 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
               << ");" << endl
               << "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", "
               << symbol_table.exo_nbr() << ");" << endl;
+  for (int hd {0}; hd < heterogeneity_table.size(); hd++)
+    mOutputFile << "M_.heterogeneity(" << hd + 1 << ").Sigma_e = zeros("
+                << symbol_table.het_exo_nbr(hd) << ", " << symbol_table.het_exo_nbr(hd) << ");"
+                << endl;
 
   if (mod_file_struct.calibrated_measurement_errors)
     mOutputFile << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", "
@@ -1040,6 +1204,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
         }
     }
 
+  for (const auto& hm : heterogeneous_models)
+    hm.writeDriverOutput(mOutputFile);
+
   if (onlymodel || gui)
     for (const auto& statement : statements)
       {
@@ -1188,6 +1355,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
       epilogue.writeEpilogueFile(basename);
 
       pac_model_table.writeTargetCoefficientsFile(basename);
+
+      for (const auto& hm : heterogeneous_models)
+        hm.writeModelFiles(basename, false);
     }
 }
 
@@ -1265,6 +1435,11 @@ ModFile::writeJsonOutputParsingCheck(const string& basename, JsonFileOutputType
 
   symbol_table.writeJsonOutput(output);
   output << ", ";
+  if (!heterogeneity_table.empty())
+    {
+      heterogeneity_table.writeJsonOutput(output);
+      output << ", ";
+    }
   dynamic_model.writeJsonOutput(output);
   output << ", ";
   static_model.writeJsonOutput(output);

src/ModFile.hh

@@ -30,6 +30,8 @@
 #include "DynamicModel.hh"
 #include "ExtendedPreprocessorTypes.hh"
 #include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
+#include "HeterogeneousModel.hh"
 #include "ModelEquationBlock.hh"
 #include "NumericalConstants.hh"
 #include "NumericalInitialization.hh"
@@ -46,6 +48,8 @@ class ModFile
 {
 public:
   explicit ModFile(WarningConsolidation& warnings_arg);
+  // For heterogeneity dimensions
+  HeterogeneityTable heterogeneity_table;
   //! Symbol table
   SymbolTable symbol_table;
   //! External Functions table
@@ -76,6 +80,8 @@ public:
   StaticModel static_model;
   //! Static model, as declared in the "steady_state_model" block if present
   SteadyStateModel steady_state_model;
+  // Heterogeneous model blocks, ordered per heterogeneity dimension ID
+  vector<HeterogeneousModel> heterogeneous_models;
 
   //! Option linear
   bool linear {false};

src/ModelEquationBlock.cc

@@ -26,8 +26,10 @@
 
 PlannerObjective::PlannerObjective(SymbolTable& symbol_table_arg,
                                    NumericalConstants& num_constants_arg,
-                                   ExternalFunctionsTable& external_functions_table_arg) :
-    StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg}
+                                   ExternalFunctionsTable& external_functions_table_arg,
+                                   HeterogeneityTable& heterogeneity_table_arg) :
+    StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+                 heterogeneity_table_arg}
 {
 }
 
@@ -38,7 +40,7 @@ PlannerObjective::writeDriverOutput(ostream& output) const
   for (const auto& it : temporary_terms_derivatives)
     output << it.size() << "; ";
   output << "];" << endl;
-  writeDriverSparseIndicesHelper<false, true>(output);
+  writeDriverSparseIndicesHelper("objective", output);
 }
 
 void
@@ -51,9 +53,14 @@ PlannerObjective::computingPassBlock([[maybe_unused]] const eval_context_t& eval
 OrigRamseyDynamicModel::OrigRamseyDynamicModel(
     SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
     ExternalFunctionsTable& external_functions_table_arg,
+    HeterogeneityTable& heterogeneity_table_arg,
     TrendComponentModelTable& trend_component_model_table_arg, VarModelTable& var_model_table_arg) :
-    DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
-                  trend_component_model_table_arg, var_model_table_arg}
+    DynamicModel {symbol_table_arg,
+                  num_constants_arg,
+                  external_functions_table_arg,
+                  heterogeneity_table_arg,
+                  trend_component_model_table_arg,
+                  var_model_table_arg}
 {
 }
 
@@ -67,8 +74,10 @@ OrigRamseyDynamicModel::operator=(const DynamicModel& m)
 SteadyStateModel::SteadyStateModel(SymbolTable& symbol_table_arg,
                                    NumericalConstants& num_constants_arg,
                                    ExternalFunctionsTable& external_functions_table_arg,
+                                   HeterogeneityTable& heterogeneity_table_arg,
                                    const StaticModel& static_model_arg) :
-    DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg},
+    DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+              heterogeneity_table_arg},
     static_model {static_model_arg}
 {
 }
@@ -327,10 +336,15 @@ SteadyStateModel::writeJsonSteadyStateFile(ostream& output, bool transformComput
 
 Epilogue::Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
                    ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg,
                    TrendComponentModelTable& trend_component_model_table_arg,
                    VarModelTable& var_model_table_arg) :
-    DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
-                  trend_component_model_table_arg, var_model_table_arg}
+    DynamicModel {symbol_table_arg,
+                  num_constants_arg,
+                  external_functions_table_arg,
+                  heterogeneity_table_arg,
+                  trend_component_model_table_arg,
+                  var_model_table_arg}
 {
 }
 

src/ModelEquationBlock.hh

@@ -30,7 +30,8 @@ class PlannerObjective : public StaticModel
 {
 public:
   PlannerObjective(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                   ExternalFunctionsTable& external_functions_table_arg);
+                   ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg);
   // NB: masks the method with the same name in StaticModel (not in a virtual fashion)
   void writeDriverOutput(ostream& output) const;
 
@@ -50,6 +51,7 @@ class OrigRamseyDynamicModel : public DynamicModel
 public:
   OrigRamseyDynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
                          ExternalFunctionsTable& external_functions_table_arg,
+                         HeterogeneityTable& heterogeneity_table_arg,
                          TrendComponentModelTable& trend_component_model_table_arg,
                          VarModelTable& var_model_table_arg);
   OrigRamseyDynamicModel& operator=(const DynamicModel& m);
@@ -75,6 +77,7 @@ private:
 public:
   SteadyStateModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
                    ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg,
                    const StaticModel& static_model_arg);
 
   SteadyStateModel(const SteadyStateModel& m);
@@ -108,6 +111,7 @@ private:
 public:
   Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
            ExternalFunctionsTable& external_functions_table_arg,
+           HeterogeneityTable& heterogeneity_table_arg,
            TrendComponentModelTable& trend_component_model_table_arg,
            VarModelTable& var_model_table_arg);
 

src/ModelTree.cc

@@ -144,8 +144,10 @@ ModelTree::copyHelper(const ModelTree& m)
 }
 
 ModelTree::ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                     ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
-    DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
+                     ExternalFunctionsTable& external_functions_table_arg,
+                     HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
+    DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+              heterogeneity_table_arg, is_dynamic_arg},
     derivatives(4),
     NNZDerivatives(4, 0),
     temporary_terms_derivatives(4)
@@ -275,8 +277,8 @@ ModelTree::computeNormalization(const jacob_map_t& contemporaneous_jacobian)
   // Create the resulting map, by copying the n first elements of mate_map, and substracting n to
   // them
   endo2eq.resize(equations.size());
-  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(),
-            [=](vertex_descriptor_t i) { return i - n; });
+  ranges::transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(),
+                    [=](vertex_descriptor_t i) { return i - n; });
 }
 
 bool
@@ -983,9 +985,9 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
   temporary_terms_derivatives.clear();
   temporary_terms_derivatives.resize(derivatives.size());
   for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
-    copy(temp_terms_map[{order, 0}].begin(), temp_terms_map[{order, 0}].end(),
-         inserter(temporary_terms_derivatives.at(order),
-                  temporary_terms_derivatives.at(order).begin()));
+    ranges::copy(temp_terms_map[{order, 0}],
+                 inserter(temporary_terms_derivatives.at(order),
+                          temporary_terms_derivatives.at(order).begin()));
 
   // Compute indices in MATLAB/Julia vector
   for (int order {0}, idx {0}; order < static_cast<int>(derivatives.size()); order++)
@@ -1038,9 +1040,9 @@ ModelTree::computeBlockTemporaryTerms(bool no_tmp_terms)
     {
       blocks_temporary_terms.at(blk).resize(temp_terms.at(blk).size());
       for (size_t i {0}; i < temp_terms.at(blk).size(); i++)
-        copy(temp_terms.at(blk).at(i).begin(), temp_terms.at(blk).at(i).end(),
-             inserter(blocks_temporary_terms.at(blk).at(i),
-                      blocks_temporary_terms.at(blk).at(i).begin()));
+        ranges::copy(temp_terms.at(blk).at(i),
+                     inserter(blocks_temporary_terms.at(blk).at(i),
+                              blocks_temporary_terms.at(blk).at(i).begin()));
     }
 
   // Compute indices in the temporary terms vector
@@ -1544,9 +1546,8 @@ ModelTree::computeParamsDerivativesTemporaryTerms()
       d->computeTemporaryTerms(order, temp_terms_map, reference_count, true);
 
   for (const auto& [order, tts] : temp_terms_map)
-    copy(temp_terms_map[order].begin(), temp_terms_map[order].end(),
-         inserter(params_derivs_temporary_terms[order],
-                  params_derivs_temporary_terms[order].begin()));
+    ranges::copy(temp_terms_map[order], inserter(params_derivs_temporary_terms[order],
+                                                 params_derivs_temporary_terms[order].begin()));
 
   for (int idx {0}; const auto& [order, tts] : params_derivs_temporary_terms)
     for (const auto& tt : tts)
@@ -1674,7 +1675,7 @@ ModelTree::findCompilerOnMacos(const string& mexext)
      Apple’s clang is located both in /usr/bin/gcc and /usr/bin/clang, it
      automatically selects x86_64 or arm64 depending on the compile-time
      environment. */
-  const string macos_gcc_version {"13"};
+  const string macos_gcc_version {"14"};
 
   if (filesystem::path global_gcc_path {"/usr/local/bin/gcc-" + macos_gcc_version};
       exists(global_gcc_path) && mexext == "mexmaci64")
@@ -1686,14 +1687,14 @@ ModelTree::findCompilerOnMacos(const string& mexext)
     return {global_clang_path, true};
   else
     {
-      cerr << "ERROR: You must install gcc-" << macos_gcc_version
+      cerr << "ERROR: You must install gcc@" << macos_gcc_version
            << " on your system before using the `use_dll` option of Dynare. "
            << "You should install Homebrew";
       if (mexext == "mexmaca64")
         cerr << " for arm64";
       else if (mexext == "mexmaci64")
         cerr << " for x86_64";
-      cerr << " and run `brew install gcc-" << macos_gcc_version << "` in a terminal." << endl;
+      cerr << " and run `brew install gcc@" << macos_gcc_version << "` in a terminal." << endl;
       exit(EXIT_FAILURE);
     }
 }
@@ -1820,8 +1821,8 @@ ModelTree::compileMEX(const filesystem::path& output_dir, const string& output_b
 
   // The prerequisites are the object files among the input files
   set<filesystem::path> prerequisites;
-  copy_if(input_files.begin(), input_files.end(), inserter(prerequisites, prerequisites.end()),
-          [](const auto& p) { return p.extension() == ".o"; });
+  ranges::copy_if(input_files, inserter(prerequisites, prerequisites.end()),
+                  [](const auto& p) { return p.extension() == ".o"; });
 
   unique_lock<mutex> lk {mex_compilation_mut};
   mex_compilation_queue.emplace_back(output_filename, prerequisites, cmd.str());
@@ -1981,8 +1982,7 @@ ModelTree::initializeMEXCompilationWorkers(int numworkers, const filesystem::pat
       auto pick_job = [&cmd, &output] {
         for (auto it {mex_compilation_queue.begin()}; it != mex_compilation_queue.end(); ++it)
           if (const auto& prerequisites {get<1>(*it)}; // Will become dangling after erase
-              includes(mex_compilation_done.begin(), mex_compilation_done.end(),
-                       prerequisites.begin(), prerequisites.end()))
+              ranges::includes(mex_compilation_done, prerequisites))
             {
               output = get<0>(*it);
               cmd = get<2>(*it);
@@ -2074,6 +2074,8 @@ ModelTree::waitForMEXCompilationWorkers()
 void
 ModelTree::computingPassBlock(const eval_context_t& eval_context, bool no_tmp_terms)
 {
+  if (!heterogeneity_table.empty())
+    return;
   if (!computeNonSingularNormalization(eval_context))
     return;
   auto [prologue, epilogue] = computePrologueAndEpilogue();
@@ -2150,3 +2152,33 @@ ModelTree::computeMCPEquationsReordering()
         swap(mcp_equations_reordering[endo_id], *it);
       }
 }
+
+void
+ModelTree::writeDriverSparseIndicesHelper(const string& prefix, ostream& output) const
+{
+  // Write indices for the sparse Jacobian (both naive and CSC storage)
+  output << "M_." << prefix << "_g1_sparse_rowval = int32([";
+  for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
+    output << indices.first + 1 << ' ';
+  output << "]);" << endl << "M_." << prefix << "_g1_sparse_colval = int32([";
+  for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
+    output << indices.second + 1 << ' ';
+  output << "]);" << endl << "M_." << prefix << "_g1_sparse_colptr = int32([";
+  for (int it : jacobian_sparse_colptr)
+    output << it + 1 << ' ';
+  output << "]);" << endl;
+
+  // Write indices for the sparse higher-order derivatives
+  for (int i {2}; i <= computed_derivs_order; i++)
+    {
+      output << "M_." << prefix << "_g" << i << "_sparse_indices = int32([";
+      for (const auto& [vidx, d] : derivatives[i])
+        {
+          for (bool row_number {true}; // First element of vidx is row number
+               int it : vidx)
+            output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 1 << ' ';
+          output << ';' << endl;
+        }
+      output << "]);" << endl;
+    }
+}

src/ParsingDriver.cc

@@ -22,6 +22,7 @@
 #include <fstream>
 #include <iostream>
 #include <numeric>
+#include <ranges>
 #include <sstream>
 
 #include "ExprNode.hh"
@@ -71,6 +72,7 @@ ParsingDriver::set_current_data_tree(DataTree* data_tree_arg)
   data_tree = data_tree_arg;
   model_tree = dynamic_cast<ModelTree*>(data_tree_arg);
   dynamic_model = dynamic_cast<DynamicModel*>(data_tree_arg);
+  heterogeneous_model = dynamic_cast<HeterogeneousModel*>(data_tree_arg);
 }
 
 void
@@ -183,12 +185,14 @@ ParsingDriver::warning(const string& m)
 
 int
 ParsingDriver::declare_symbol(const string& name, SymbolType type, const string& tex_name,
-                              const vector<pair<string, string>>& partition_value)
+                              const vector<pair<string, string>>& partition_value,
+                              const optional<int>& heterogeneity_dimension)
 {
   int symb_id;
   try
     {
-      symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value);
+      symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value,
+                                                 heterogeneity_dimension);
     }
   catch (SymbolTable::AlreadyDeclaredException& e)
     {
@@ -207,16 +211,30 @@ int
 ParsingDriver::declare_endogenous(const string& name, const string& tex_name,
                                   const vector<pair<string, string>>& partition_value)
 {
-  return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value);
+  return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value, {});
 }
 
 void
 ParsingDriver::var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
-                   bool log_option)
+                   const optional<string>& heterogeneity_dimension, bool log_option)
 {
   for (auto& [name, tex_name, partition] : symbol_list)
     {
-      int symb_id = declare_endogenous(name, tex_name, partition);
+      int symb_id {[&] {
+        if (heterogeneity_dimension)
+          try
+            {
+              return declare_symbol(name, SymbolType::heterogeneousEndogenous, tex_name, partition,
+                                    mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+            }
+          catch (HeterogeneityTable::UnknownDimensionNameException&)
+            {
+              error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+            }
+        else
+          return declare_endogenous(name, tex_name, partition);
+      }()};
+
       if (log_option)
         mod_file->symbol_table.markWithLogTransform(symb_id);
     }
@@ -226,15 +244,27 @@ int
 ParsingDriver::declare_exogenous(const string& name, const string& tex_name,
                                  const vector<pair<string, string>>& partition_value)
 {
-  return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value);
+  return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value, {});
 }
 
 void
 ParsingDriver::varexo(
-    const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
+    const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+    const optional<string>& heterogeneity_dimension)
 {
   for (auto& [name, tex_name, partition] : symbol_list)
-    declare_exogenous(name, tex_name, partition);
+    if (heterogeneity_dimension)
+      try
+        {
+          declare_symbol(name, SymbolType::heterogeneousExogenous, tex_name, partition,
+                         mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+        }
+      catch (HeterogeneityTable::UnknownDimensionNameException&)
+        {
+          error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+        }
+    else
+      declare_exogenous(name, tex_name, partition);
 }
 
 void
@@ -242,22 +272,34 @@ ParsingDriver::varexo_det(
     const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
 {
   for (auto& [name, tex_name, partition] : symbol_list)
-    declare_symbol(name, SymbolType::exogenousDet, tex_name, partition);
+    declare_symbol(name, SymbolType::exogenousDet, tex_name, partition, {});
 }
 
 int
 ParsingDriver::declare_parameter(const string& name, const string& tex_name,
                                  const vector<pair<string, string>>& partition_value)
 {
-  return declare_symbol(name, SymbolType::parameter, tex_name, partition_value);
+  return declare_symbol(name, SymbolType::parameter, tex_name, partition_value, {});
 }
 
 void
 ParsingDriver::parameters(
-    const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
+    const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+    const optional<string>& heterogeneity_dimension)
 {
   for (auto& [name, tex_name, partition] : symbol_list)
-    declare_parameter(name, tex_name, partition);
+    if (heterogeneity_dimension)
+      try
+        {
+          declare_symbol(name, SymbolType::heterogeneousParameter, tex_name, partition,
+                         mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+        }
+      catch (HeterogeneityTable::UnknownDimensionNameException&)
+        {
+          error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+        }
+    else
+      declare_parameter(name, tex_name, partition);
 }
 
 void
@@ -266,7 +308,7 @@ ParsingDriver::declare_statement_local_variable(const string& name)
   if (mod_file->symbol_table.exists(name))
     error("Symbol " + name + " cannot be assigned within a statement "
           + "while being assigned elsewhere in the modfile");
-  declare_symbol(name, SymbolType::statementDeclaredVariable, "", {});
+  declare_symbol(name, SymbolType::statementDeclaredVariable, "", {}, {});
 }
 
 void
@@ -293,7 +335,7 @@ ParsingDriver::end_trend_var(bool log_trend, expr_t growth_factor,
   for (auto& [name, tex_name] : symbol_list)
     {
       int symb_id = declare_symbol(name, log_trend ? SymbolType::logTrend : SymbolType::trend,
-                                   tex_name, {});
+                                   tex_name, {}, {});
       declared_trend_vars.push_back(symb_id);
     }
 
@@ -418,10 +460,6 @@ ParsingDriver::add_model_variable(int symb_id, int lag)
     error("Exogenous deterministic variable " + mod_file->symbol_table.getName(symb_id)
           + " cannot be given a lead or a lag.");
 
-  if (type == SymbolType::modelLocalVariable && lag != 0)
-    error("Model local variable " + mod_file->symbol_table.getName(symb_id)
-          + " cannot be given a lead or a lag.");
-
   if (data_tree == planner_objective.get())
     {
       if (lag != 0)
@@ -840,7 +878,7 @@ ParsingDriver::end_epilogue()
 void
 ParsingDriver::add_epilogue_variable(const string& name)
 {
-  declare_symbol(name, SymbolType::epilogue, "", {});
+  declare_symbol(name, SymbolType::epilogue, "", {}, {});
 }
 
 void
@@ -855,6 +893,22 @@ ParsingDriver::begin_model()
   set_current_data_tree(&mod_file->dynamic_model);
 }
 
+void
+ParsingDriver::begin_heterogeneous_model(const string& heterogeneity_dimension)
+{
+  int het_dim_id {[&] {
+    try
+      {
+        return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
+      }
+    catch (HeterogeneityTable::UnknownDimensionNameException&)
+      {
+        error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
+      }
+  }()};
+  set_current_data_tree(&mod_file->heterogeneous_models.at(het_dim_id));
+}
+
 void
 ParsingDriver::end_model()
 {
@@ -904,6 +958,28 @@ ParsingDriver::end_shocks(bool overwrite)
   corr_shocks.clear();
 }
 
+void
+ParsingDriver::end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite)
+{
+  int het_dim_id {[&] {
+    try
+      {
+        return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
+      }
+    catch (HeterogeneityTable::UnknownDimensionNameException&)
+      {
+        error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
+      }
+  }()};
+  mod_file->addStatement(make_unique<HeterogeneousShocksStatement>(
+      het_dim_id, overwrite, move(var_shocks), move(std_shocks), move(covar_shocks),
+      move(corr_shocks), mod_file->symbol_table, mod_file->heterogeneity_table));
+  var_shocks.clear();
+  std_shocks.clear();
+  covar_shocks.clear();
+  corr_shocks.clear();
+}
+
 void
 ParsingDriver::end_mshocks(bool overwrite, bool relative_to_initval)
 {
@@ -1292,7 +1368,7 @@ void
 ParsingDriver::add_positive_restriction_element(expr_t value, const string& variable,
                                                 const string& lag)
 {
-  // if the expression is not on the left handside, change its sign
+  // if the expression is not on the left-hand side, change its sign
   if (!svar_left_handside)
     value = add_uminus(value);
 
@@ -1304,7 +1380,7 @@ ParsingDriver::add_positive_restriction_element(const string& variable, const st
 {
   expr_t value(data_tree->One);
 
-  // if the expression is not on the left handside, change its sign
+  // if the expression is not on the left-hand side, change its sign
   if (!svar_left_handside)
     value = add_uminus(value);
 
@@ -1315,7 +1391,7 @@ void
 ParsingDriver::add_negative_restriction_element(expr_t value, const string& variable,
                                                 const string& lag)
 {
-  // if the expression is on the left handside, change its sign
+  // if the expression is on the left-hand side, change its sign
   if (svar_left_handside)
     value = add_uminus(value);
 
@@ -1327,7 +1403,7 @@ ParsingDriver::add_negative_restriction_element(const string& variable, const st
 {
   expr_t value(data_tree->One);
 
-  // if the expression is on the left handside, change its sign
+  // if the expression is on the left-hand side, change its sign
   if (svar_left_handside)
     value = add_uminus(value);
 
@@ -2224,7 +2300,8 @@ void
 ParsingDriver::begin_planner_objective()
 {
   planner_objective = make_unique<PlannerObjective>(mod_file->symbol_table, mod_file->num_constants,
-                                                    mod_file->external_functions_table);
+                                                    mod_file->external_functions_table,
+                                                    mod_file->heterogeneity_table);
   set_current_data_tree(planner_objective.get());
 }
 
@@ -2657,6 +2734,12 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
           error("'mcp' tag does not contain an inequality");
       }()};
 
+      // Trim whitespace
+      var_name.erase(var_name.find_last_not_of(" \n\t") + 1);
+      var_name.erase(0, var_name.find_first_not_of(" \n\t"));
+      constant.erase(constant.find_last_not_of(" \n\t") + 1);
+      constant.erase(0, constant.find_first_not_of(" \n\t"));
+
       int symb_id {[&] {
         try
           {
@@ -2664,18 +2747,18 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
           }
         catch (SymbolTable::UnknownSymbolNameException&)
           {
-            error("Left hand-side of expression in 'mcp' tag is not a variable");
+            error("Left-hand side of expression in 'mcp' tag is not a variable");
           }
       }()};
 
       if (mod_file->symbol_table.getType(symb_id) != SymbolType::endogenous)
-        error("Left hand-side of expression in 'mcp' tag is not an endogenous variable");
+        error("Left-hand side of expression in 'mcp' tag is not an endogenous variable");
 
       expr_t matched_constant {[&] {
         char* str_end;
         double d = strtod(constant.c_str(), &str_end);
         if (str_end == constant.c_str())
-          error("Right hand-side of expression in 'mcp' tag should be a constant");
+          error("Right-hand side of expression in 'mcp' tag should be a constant");
         return data_tree->AddPossiblyNegativeConstant(d);
       }()};
 
@@ -2795,7 +2878,7 @@ void
 ParsingDriver::model_local_variable(const vector<pair<string, string>>& symbol_list)
 {
   for (auto& [name, tex_name] : symbol_list)
-    declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {});
+    declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {}, {});
 }
 
 void
@@ -3032,9 +3115,8 @@ ParsingDriver::add_diff(expr_t arg1)
 expr_t
 ParsingDriver::add_adl(expr_t arg1, const string& name, const string& lag)
 {
-  vector<int> lags(stoi(lag));
-  iota(lags.begin(), lags.end(), 1);
-  return add_adl(arg1, name, lags);
+  auto lags = views::iota(1, stoi(lag) + 1);
+  return add_adl(arg1, name, vector<int> {lags.begin(), lags.end()});
 }
 
 expr_t
@@ -3215,6 +3297,23 @@ ParsingDriver::add_steady_state(expr_t arg1)
   return data_tree->AddSteadyState(arg1);
 }
 
+expr_t
+ParsingDriver::add_sum(expr_t arg)
+{
+  if (heterogeneous_model)
+    error("The SUM() operator cannot be used inside a model(heterogeneity=...) block");
+
+  VariableNode* varg {dynamic_cast<VariableNode*>(arg)};
+  if (!varg)
+    error("The argument to the SUM() operator must be a single variable");
+  if (varg->lag != 0)
+    error("The argument to the SUM() operator must not have a lead or lag");
+  if (mod_file->symbol_table.getType(varg->symb_id) != SymbolType::heterogeneousEndogenous)
+    error("The argument to the SUM() operator must be a heterogeneous endogenous variable");
+
+  return data_tree->AddSum(arg);
+}
+
 void
 ParsingDriver::external_function_option(const string& name_option, const string& opt)
 {
@@ -3223,7 +3322,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
       if (opt.empty())
         error("An argument must be passed to the 'name' option of the external_function() "
               "statement.");
-      declare_symbol(opt, SymbolType::externalFunction, "", {});
+      declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
       current_external_function_id = mod_file->symbol_table.getID(opt);
     }
   else if (name_option == "first_deriv_provided")
@@ -3233,7 +3332,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
             = ExternalFunctionsTable::IDSetButNoNameProvided;
       else
         {
-          int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
+          int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
           current_external_function_options.firstDerivSymbID = symb_id;
         }
     }
@@ -3244,7 +3343,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
             = ExternalFunctionsTable::IDSetButNoNameProvided;
       else
         {
-          int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
+          int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
           current_external_function_options.secondDerivSymbID = symb_id;
         }
     }
@@ -3410,7 +3509,7 @@ ParsingDriver::add_model_var_or_external_function(const string& function_name, b
                   + ") within the model block, you must first declare it via the "
                     "external_function() statement.");
         }
-      int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {});
+      int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {}, {});
       current_external_function_options.nargs = stack_external_function_args.top().size();
       mod_file->external_functions_table.addExternalFunction(
           symb_id, current_external_function_options, in_model_block);
@@ -3836,7 +3935,8 @@ ParsingDriver::begin_occbin_constraints()
      if added to the main DynamicModel tree. It also simplifies the
      enforcement of various constraints at parsing time. */
   occbin_constraints_tree = make_unique<DataTree>(mod_file->symbol_table, mod_file->num_constants,
-                                                  mod_file->external_functions_table, false);
+                                                  mod_file->external_functions_table,
+                                                  mod_file->heterogeneity_table, false);
   set_current_data_tree(occbin_constraints_tree.get());
 }
 
@@ -3971,3 +4071,26 @@ ParsingDriver::matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_we
 {
   mod_file->addStatement(make_unique<MatchedIrfsWeightsStatement>(move(weights), overwrite));
 }
+
+void
+ParsingDriver::heterogeneity_dimension(const vector<string>& dims)
+{
+  for (const auto& dim : dims)
+    {
+      int het_dim_id {[&] {
+        try
+          {
+            return mod_file->heterogeneity_table.addDimension(dim);
+          }
+        catch (HeterogeneityTable::AlreadyDeclaredDimensionException&)
+          {
+            error("Heterogeneity dimension '" + dim + "' already declared");
+          }
+      }()};
+
+      assert(static_cast<int>(mod_file->heterogeneous_models.size()) == het_dim_id);
+      mod_file->heterogeneous_models.emplace_back(mod_file->symbol_table, mod_file->num_constants,
+                                                  mod_file->external_functions_table,
+                                                  mod_file->heterogeneity_table, het_dim_id);
+    }
+}

src/ParsingDriver.hh

@@ -114,7 +114,8 @@ private:
 
   //! Helper to add a symbol declaration (returns its symbol ID)
   int declare_symbol(const string& name, SymbolType type, const string& tex_name,
-                     const vector<pair<string, string>>& partition_value);
+                     const vector<pair<string, string>>& partition_value,
+                     const optional<int>& heterogeneity_dimension);
 
   //! Temporary store for the planner objective
   unique_ptr<PlannerObjective> planner_objective;
@@ -137,6 +138,11 @@ private:
    * DynamicModel instance */
   DynamicModel* dynamic_model;
 
+  //! The heterogeneous model tree in which to add expressions currently parsed
+  /*! It is only a dynamic cast of data_tree pointer, and is therefore null if data_tree is not a
+   * HeterogeneousModel instance */
+  HeterogeneousModel* heterogeneous_model;
+
   //! Sets data_tree and model_tree pointers
   void set_current_data_tree(DataTree* data_tree_arg);
 
@@ -314,9 +320,9 @@ public:
   expr_t var_expectation_model_discount {nullptr};
 
   //! Error handler with explicit location
-  void error(const Dynare::parser::location_type& l, const string& m) __attribute__((noreturn));
+  [[noreturn]] void error(const Dynare::parser::location_type& l, const string& m);
   //! Error handler using saved location
-  void error(const string& m) __attribute__((noreturn));
+  [[noreturn]] void error(const string& m);
   //! Warning handler using saved location
   void warning(const string& m);
 
@@ -371,19 +377,21 @@ public:
                          const vector<pair<string, string>>& partition_value = {});
   // Handles a “var” or “var(log)” statement (without “deflator” or “log_deflator” options)
   void var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
-           bool log_option);
+           const optional<string>& heterogeneity_dimension, bool log_option);
   //! Declares an exogenous variable (and returns its symbol ID)
   int declare_exogenous(const string& name, const string& tex_name = "",
                         const vector<pair<string, string>>& partition_value = {});
   // Handles a “varexo” statement
-  void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
+  void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+              const optional<string>& heterogeneity_dimension);
   // Handles a “varexo_det” statement
   void varexo_det(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
   //! Declares a parameter (and returns its symbol ID)
   int declare_parameter(const string& name, const string& tex_name = "",
                         const vector<pair<string, string>>& partition_value = {});
   // Handles a “parameters” statement
-  void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
+  void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+                  const optional<string>& heterogeneity_dimension);
   // Handles a “model_local_variable” statement
   void model_local_variable(const vector<pair<string, string>>& symbol_list);
   //! Declares a statement local variable
@@ -454,10 +462,11 @@ public:
   void add_epilogue_variable(const string& varname);
   //! Add equation in epilogue block
   void add_epilogue_equal(const string& varname, expr_t expr);
-  /* Begin a model or model_replace block, or an expression as an option value
-     of some statement.
-     Must be followed by a call to reset_data_tree(). */
+  /* Begin a model block (without heterogeneity option), or an expression as an option value of some
+     statement. Must be followed by a call to reset_data_tree(). */
   void begin_model();
+  // Begin a model(heterogeneity=…) block
+  void begin_heterogeneous_model(const string& heterogeneity_dimension);
   //! End a model or model_replace block, printing errors that were encountered in parsing
   void end_model();
   //! Writes a shocks statement
@@ -468,6 +477,8 @@ public:
   void end_shocks_surprise(bool overwrite);
   //! Writes a shocks(learnt_in=…) block
   void end_shocks_learnt_in(const string& learnt_in_period, bool overwrite);
+  // For a shocks(heterogeneity=…) block
+  void end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite);
   //! Writes a mshocks(learnt_in=…) block
   void end_mshocks_learnt_in(const string& learnt_in_period, bool overwrite,
                              bool relative_to_initval);
@@ -844,6 +855,8 @@ public:
   expr_t add_erfc(expr_t arg);
   //! Writes token "steadyState(arg1)" to model tree
   expr_t add_steady_state(expr_t arg1);
+  // Add a “sum(arg)” node to model tree
+  expr_t add_sum(expr_t arg);
   //! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)
   void push_external_function_arg_vector_onto_stack();
   //! Adds an external function argument
@@ -952,6 +965,8 @@ public:
   // Add a matched_irfs_weights block
   void matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_weights_t weights,
                             bool overwrite);
+  void heterogeneity_dimension(const vector<string>& dims);
+
   // Returns true iff the string is a legal symbol identifier (see NAME token in lexer)
   static bool isSymbolIdentifier(const string& str);
   // Given an Occbin regime name, returns the corresponding auxiliary parameter

src/Shocks.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -25,6 +25,7 @@
 #include <vector>
 
 #include "ExprNode.hh"
+#include "HeterogeneityTable.hh"
 #include "Statement.hh"
 #include "SymbolTable.hh"
 
@@ -155,6 +156,47 @@ public:
   void writeJsonOutput(ostream& output) const override;
 };
 
+class HeterogeneousShocksStatement : public Statement
+{
+public:
+  const int heterogeneity_dimension;
+  const bool overwrite;
+
+  using var_and_std_shocks_t = map<int, expr_t>;
+  using covar_and_corr_shocks_t = map<pair<int, int>, expr_t>;
+
+  const var_and_std_shocks_t var_shocks, std_shocks;
+  const covar_and_corr_shocks_t covar_shocks, corr_shocks;
+
+private:
+  const SymbolTable& symbol_table;
+  const HeterogeneityTable& heterogeneity_table;
+
+  void writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it, bool stddev) const;
+  void writeVarAndStdShocks(ostream& output) const;
+  void writeCovarOrCorrShock(ostream& output, const pair<pair<int, int>, expr_t>& it,
+                             bool corr) const;
+  void writeCovarAndCorrShocks(ostream& output) const;
+
+  string
+  sigmaeName() const
+  {
+    return "M_.heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").Sigma_e"s;
+  }
+
+public:
+  HeterogeneousShocksStatement(int heterogeneity_dimension_arg, bool overwrite_arg,
+                               var_and_std_shocks_t var_shocks_arg,
+                               var_and_std_shocks_t std_shocks_arg,
+                               covar_and_corr_shocks_t covar_shocks_arg,
+                               covar_and_corr_shocks_t corr_shocks_arg,
+                               const SymbolTable& symbol_table_arg,
+                               const HeterogeneityTable& heterogeneity_table_arg);
+  void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
+  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream& output) const override;
+};
+
 class ConditionalForecastPathsStatement : public Statement
 {
 private:

src/Statement.cc

@@ -50,6 +50,11 @@ NativeStatement::writeOutput(ostream& output, [[maybe_unused]] const string& bas
      boost::xpressive does not look for the longest match in an alternation, but stops at the first
      match from left to right. */
   string date_regex = R"((-?\d+([YyAa]|[Mm](1[0-2]|[1-9])|[Qq][1-4]|[SsHh][1-2])))";
+  /* NB: the following dance around the dollar sign (exclude it from lookbehind, then use it in a
+     temporary string after the first replace, then remove it in the second replace) has a purpose:
+     it allows the user to disable the substitution mechanism. For example, if the user writes
+     “$2024Q4” in a native statement, it will be transformed into “2024Q4” and not
+     “$dates('2024Q4')”. */
   sregex regex_lookbehind = sregex::compile(R"((?<!\$|\d|[a-zA-Z_]|-|'))" + date_regex);
   sregex regex_dollar = sregex::compile(R"((\$))" + date_regex);
 

src/StaticModel.hh

@@ -121,7 +121,8 @@ protected:
 
 public:
   StaticModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants,
-              ExternalFunctionsTable& external_functions_table_arg);
+              ExternalFunctionsTable& external_functions_table_arg,
+              HeterogeneityTable& heterogeneity_table_arg);
 
   StaticModel(const StaticModel& m);
   StaticModel& operator=(const StaticModel& m);

src/SymbolList.cc

 /*
- * Copyright © 2003-2022 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -57,8 +57,8 @@ SymbolList::checkPass(WarningConsolidation& warnings, const vector<SymbolType>&
             throw SymbolListException {"Variable " + symbol + " was not declared."};
         }
 
-      if (none_of(types.begin(), types.end(),
-                  [&](SymbolType type) { return symbol_table.getType(symbol) == type; }))
+      if (ranges::none_of(types,
+                          [&](SymbolType type) { return symbol_table.getType(symbol) == type; }))
         {
           string valid_types;
           for (auto type : types)
@@ -103,6 +103,15 @@ SymbolList::checkPass(WarningConsolidation& warnings, const vector<SymbolType>&
               case SymbolType::excludedVariable:
                 valid_types += "excludedVariable, ";
                 break;
+              case SymbolType::heterogeneousEndogenous:
+                valid_types += "heterogeneousEndogenous, ";
+                break;
+              case SymbolType::heterogeneousExogenous:
+                valid_types += "heterogeneousExogenous, ";
+                break;
+              case SymbolType::heterogeneousParameter:
+                valid_types += "heterogeneousParameter, ";
+                break;
               }
           valid_types = valid_types.erase(valid_types.size() - 2, 2);
           throw SymbolListException {"Variable " + symbol + " is not one of {" + valid_types + "}"};

src/macro/Directives.cc

@@ -189,10 +189,7 @@ void
 EchoMacroVars::interpret(ostream& output, Environment& env,
                          [[maybe_unused]] vector<filesystem::path>& paths)
 {
-  if (save)
-    env.print(output, vars, location.begin.line, true);
-  else
-    env.print(cout, vars);
+  env.print(save ? output : cout, vars, location.begin.line, save);
   printEndLineInfo(output);
 }