Compare revisions

Sébastien Villemot · Stéphane Adjemian · Sébastien Villemot · Sébastien Villemot · Sébastien Villemot · Sébastien Villemot
--- a/.clang-format
+++ b/.clang-format
@@ -11,9 +11,12 @@ BreakInheritanceList: AfterColon
 Cpp11BracedListStyle: true
 DeriveLineEnding: false
 IndentPPDirectives: AfterHash
+InsertNewlineAtEOF: true
 PackConstructorInitializers: NextLine
 PPIndentWidth: 1
 PointerAlignment: Left
+# RemoveParentheses: ReturnStatement
+# RemoveSemicolon: true
 SpaceAfterTemplateKeyword: false
 SpaceBeforeParens: ControlStatements
 SpaceBeforeCpp11BracedList: true
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
 variables:
  TERM: linux
-  MINGW64_BOOST_VERSION: 1.85.0-2
+  MINGW64_BOOST_VERSION: 1.86.0-7
  WGET_OPTIONS: '--no-verbose --no-use-server-timestamps --retry-connrefused --retry-on-host-error'
  # To ensure that "false && true" fails, see https://gitlab.com/gitlab-org/gitlab-runner/-/issues/25394#note_412609647
  FF_ENABLE_BASH_EXIT_CODE_CHECK: 'true'
@@ -77,3 +77,16 @@ test_clang_format:
    - meson setup build-clang-format
    - ninja -C build-clang-format clang-format-check
  needs: []
+
+test_clang_tidy:
+  stage: test
+  script:
+    # Hack needed for meson < 1.6.0 which only looks for unversioned clang-tidy
+    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-19 ~/.local/bin/clang-tidy
+    - export PATH="$HOME/.local/bin:$PATH"
+    - meson setup build-clang-tidy
+    # Generate Flex and Bison files
+    - meson compile -C build-clang-tidy
+    - ninja -C build-clang-tidy clang-tidy
+  needs: []
+  when: manual
--- a/scripts/wasm.ini
+++ b/scripts/wasm.ini
+# Meson cross file for creating a WebAssembly version of the preprocessor.
+#
+# Requires emscripten to be installed.
+# Was successfully tested with emscripten 3.1.69 installed through emsdk
+# tool, as described on: https://emscripten.org/docs/getting_started/downloads.html
+# Don’t forget to source script snippet in current shell before running meson.
+#
+# Compilation creates a .wasm and .js wrapper under <builddir>/src/
+#
+# Can be run locally with node.js using:
+#  node dynare-preprocessor.js file.mod
+# NB: a version of node.js is shipped with emscripten (under the node/
+# subdirectory), but another version should also work.
+
+[binaries]
+cpp = 'em++'
+
+[host_machine]
+system = 'emscripten'
+# Could be changed to wasm64 if 4GB memory constraint is hit
+# Some background: https://v8.dev/blog/4gb-wasm-memory
+cpu_family = 'wasm32'
+cpu = 'wasm32'
+endian = 'little'
+
+[built-in options]
+# Never do a debug build, because otherwise the lack of optimisations can
+# overflow the memory capacities.
+buildtype = 'release'
+# The -fexceptions flag (for both compilation and linking) is needed for an
+# unknown reason (C++ compilers are supposed to always add exception support).
+# The -Wno-unqualified-std-cast-call flag removes many warnings about “move”
+# not being qualified with “std::” namespace.
+# The -fexperimental-library flag is needed to get std::jthread support (it was
+# supposed to no longer be necessary for LLVM 20, but for some reason we still
+# need it).
+cpp_args = [ '-fexceptions', '-Wno-unqualified-std-cast-call', '-fexperimental-library' ]
+# NODERAWFS=1 is needed for accessing the local filesystem
+cpp_link_args = [ '-s', 'NODERAWFS=1', '-fexceptions' ]
+
+[properties]
+# It’s necessary to use a different copy of Boost than the one under
+# /usr/include, because otherwise GCC headers confuse Clang
+boost_root = '/tmp/boost_1_86_0'
--- a/src/CommonEnums.hh
+++ b/src/CommonEnums.hh
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -25,14 +25,17 @@
 * command */
 enum class SymbolType
 {
-  endogenous = 0,            //!< Endogenous
-  exogenous = 1,             //!< Exogenous
-  exogenousDet = 2,          //!< Exogenous deterministic
-  parameter = 4,             //!< Parameter
-  modelLocalVariable = 10,   //!< Local variable whose scope is model (pound expression)
-  modFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
-  externalFunction = 12,     //!< External (user-defined) function
-  trend = 13,                //!< Trend variable
+  endogenous = 0,              // Endogenous (non-heterogeneous)
+  exogenous = 1,               // Exogenous (non-heterogeneous)
+  exogenousDet = 2,            // Exogenous deterministic (non-heterogeneous)
+  parameter = 4,               // Parameter (non-heterogeneous)
+  heterogeneousEndogenous = 5, // Endogenous that is heterogeneous across some dimension
+  heterogeneousExogenous = 6,  // Exogenous that is heterogeneous across some dimension
+  heterogeneousParameter = 7,  // Parameter that is heterogeneous across some dimension
+  modelLocalVariable = 10,     // Local variable whose scope is model (pound expression)
+  modFileLocalVariable = 11,   // Local variable whose scope is mod file (model excluded)
+  externalFunction = 12,       // External (user-defined) function
+  trend = 13,                  // Trend variable
  statementDeclaredVariable
  = 14, //!< Local variable assigned within a Statement (see subsample statement for example)
  logTrend = 15, //!< Log-trend variable
@@ -45,6 +48,13 @@ enum class SymbolType
  excludedVariable = 19 //!< Variable excluded via model_remove/var_remove/include_eqs/exclude_eqs
 };

+constexpr bool
+isHeterogeneous(SymbolType type)
+{
+  return type == SymbolType::heterogeneousEndogenous || type == SymbolType::heterogeneousExogenous
+         || type == SymbolType::heterogeneousParameter;
+}
+
 enum class UnaryOpcode
 {
  uminus,
@@ -75,7 +85,8 @@ enum class UnaryOpcode
  erf,
  erfc,
  diff,
-  adl
+  adl,
+  sum
 };

 enum class BinaryOpcode

--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -283,7 +283,9 @@ PerfectForesightWithExpectationErrorsSolverStatement::writeOutput(
    [[maybe_unused]] bool minimal_workspace) const
 {
  options_list.writeOutput(output);
-  output << "oo_ = perfect_foresight_with_expectation_errors_solver(M_, options_, oo_);" << endl;
+  output << "[oo_, Simulated_time_series] = perfect_foresight_with_expectation_errors_solver(M_, "
+            "options_, oo_);"
+         << endl;
 }

 void
@@ -3713,8 +3715,7 @@ SvarGlobalIdentificationCheckStatement::writeJsonOutput(ostream& output) const
  output << R"({"statementName": "svar_global_identification"})";
 }

-SetTimeStatement::SetTimeStatement(OptionsList options_list_arg) :
-    options_list {move(options_list_arg)}
+SetTimeStatement::SetTimeStatement(string period_arg) : period {move(period_arg)}
 {
 }

@@ -3722,19 +3723,13 @@ void
 SetTimeStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                              [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output);
+  output << "options_.initial_period = " << period << endl;
 }

 void
 SetTimeStatement::writeJsonOutput(ostream& output) const
 {
-  output << R"({"statementName": "set_time")";
-  if (!options_list.empty())
-    {
-      output << ", ";
-      options_list.writeJsonOutput(output);
-    }
-  output << "}";
+  output << R"({"statementName": "set_time", "period": ")" << period << R"("})";
 }

 EstimationDataStatement::EstimationDataStatement(OptionsList options_list_arg) :
@@ -5347,8 +5342,8 @@ void
 ResidStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                            [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output, "options_resid_");
-  output << "display_static_residuals(M_, options_, oo_, options_resid_);" << endl;
+  options_list.writeOutput(output);
+  output << "display_static_residuals(M_, options_, oo_);" << endl;
 }

 void

--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -524,7 +524,7 @@ public:
  blockName() const override
  {
    return "estimated_params";
-  };
+  }
  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;
@@ -542,7 +542,7 @@ public:
  blockName() const override
  {
    return "estimated_params_init";
-  };
+  }
  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;
@@ -557,7 +557,7 @@ public:
  blockName() const override
  {
    return "estimated_params_bounds";
-  };
+  }
  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;
@@ -1007,10 +1007,10 @@ public:
 class SetTimeStatement : public Statement
 {
 private:
-  const OptionsList options_list;
+  const string period;

 public:
-  explicit SetTimeStatement(OptionsList options_list_arg);
+  explicit SetTimeStatement(string period_arg);
  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
  void writeJsonOutput(ostream& output) const override;
 };

--- a/src/Configuration.hh
+++ b/src/Configuration.hh
@@ -47,7 +47,7 @@ private:
    get_paths() const
    {
      return paths;
-    };
+    }

  private:
    map<string, vector<string>> paths;

--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -24,6 +24,7 @@
 #include <fstream>
 #include <iostream>
 #include <iterator>
+#include <ranges>

 #include "DataTree.hh"

@@ -47,10 +48,12 @@ DataTree::initConstants()
 }

 DataTree::DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                   ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
+                   ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
    symbol_table {symbol_table_arg},
    num_constants {num_constants_arg},
    external_functions_table {external_functions_table_arg},
+    heterogeneity_table {heterogeneity_table_arg},
    is_dynamic {is_dynamic_arg}
 {
  initConstants();
@@ -60,19 +63,21 @@ DataTree::DataTree(const DataTree& d) :
    symbol_table {d.symbol_table},
    num_constants {d.num_constants},
    external_functions_table {d.external_functions_table},
+    heterogeneity_table {d.heterogeneity_table},
    is_dynamic {d.is_dynamic},
    local_variables_vector {d.local_variables_vector}
 {
  // Constants must be initialized first because they are used in some Add* methods
  initConstants();

+  // See commment in DataTree::operator=() for the rationale
+  for (int symb_id : d.local_variables_vector)
+    local_variables_table[symb_id] = d.local_variables_table.at(symb_id)->clone(*this);
+
  for (const auto& it : d.node_list)
    it->clone(*this);

  assert(node_list.size() == d.node_list.size());
-
-  for (const auto& [symb_id, value] : d.local_variables_table)
-    local_variables_table[symb_id] = value->clone(*this);
 }

 DataTree&
@@ -81,6 +86,7 @@ DataTree::operator=(const DataTree& d)
  assert(&symbol_table == &d.symbol_table);
  assert(&num_constants == &d.num_constants);
  assert(&external_functions_table == &d.external_functions_table);
+  assert(&heterogeneity_table == &d.heterogeneity_table);
  assert(is_dynamic == d.is_dynamic);

  num_const_node_map.clear();
@@ -791,12 +797,18 @@ DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<exp
  return p;
 }

+expr_t
+DataTree::AddSum(expr_t arg)
+{
+  return AddUnaryOp(UnaryOpcode::sum, arg);
+}
+
 bool
 DataTree::isSymbolUsed(int symb_id) const
 {
-  for (const auto& [symb_lag, expr] : variable_node_map)
-    if (symb_lag.first == symb_id)
-      return true;
+  if (ranges::any_of(views::keys(variable_node_map),
+                     [=](const auto& symb_lag) { return symb_lag.first == symb_id; }))
+    return true;

  if (local_variables_table.contains(symb_id))
    return true;
@@ -842,18 +854,18 @@ DataTree::addAllParamDerivId([[maybe_unused]] set<int>& deriv_id_set)
 bool
 DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 {
-  return ranges::any_of(unary_op_node_map,
-                        [=](const auto& it) { return get<1>(it.first) == opcode; });
+  return ranges::any_of(views::keys(unary_op_node_map),
+                        [=](const auto& key) { return get<1>(key) == opcode; });
 }

 bool
 DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 {
  set<int> var;
-  for (const auto& it : unary_op_node_map)
-    if (get<1>(it.first) == opcode)
+  for (const auto& [key, value] : unary_op_node_map)
+    if (get<1>(key) == opcode)
      {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
        if (!var.empty())
          return true;
      }
@@ -863,18 +875,18 @@ DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 bool
 DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
 {
-  return ranges::any_of(binary_op_node_map,
-                        [=](const auto& it) { return get<2>(it.first) == opcode; });
+  return ranges::any_of(views::keys(binary_op_node_map),
+                        [=](const auto& key) { return get<2>(key) == opcode; });
 }

 bool
 DataTree::isBinaryOpUsedOnType(SymbolType type, BinaryOpcode opcode) const
 {
  set<int> var;
-  for (const auto& it : binary_op_node_map)
-    if (get<2>(it.first) == opcode)
+  for (const auto& [key, value] : binary_op_node_map)
+    if (get<2>(key) == opcode)
      {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
        if (!var.empty())
          return true;
      }

--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -33,6 +33,7 @@

 #include "ExprNode.hh"
 #include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
 #include "NumericalConstants.hh"
 #include "SubModel.hh"
 #include "SymbolTable.hh"
@@ -48,6 +49,8 @@ public:
  NumericalConstants& num_constants;
  //! A reference to the external functions table
  ExternalFunctionsTable& external_functions_table;
+  // A reference to the heterogeneity table
+  HeterogeneityTable& heterogeneity_table;
  //! Is it possible to use leads/lags on variable nodes?
  /* NB: This data member cannot be replaced by a virtual method, because this information is needed
     in AddVariable(), which itself can be called from the copy constructor. */
@@ -139,7 +142,8 @@ private:

 public:
  DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-           ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
+           ExternalFunctionsTable& external_functions_table_arg,
+           HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);

  virtual ~DataTree() = default;

@@ -275,6 +279,9 @@ public:
  //! Adds an external function node for the second derivative of an external function
  expr_t AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t>& arguments,
                                        int input_index1, int input_index2);
+  // Adds "SUM(arg)" to model tree
+  expr_t AddSum(expr_t arg);
+
  //! Checks if a given symbol is used somewhere in the data tree
  [[nodiscard]] bool isSymbolUsed(int symb_id) const;
  //! Checks if a given unary op is used somewhere in the data tree
@@ -360,7 +367,13 @@ public:
    if (it == local_variables_table.end())
      throw UnknownLocalVariableException {symb_id};

-    return it->second->decreaseLeadsLags(-lead_lag);
+    /* In the following, the case without lead/lag is optimized. It makes a difference on models
+       with many nested model-local variables, see e.g.
+       https://forum.dynare.org/t/pre-processing-takes-very-long/26865 */
+    if (lead_lag == 0)
+      return it->second;
+    else
+      return it->second->decreaseLeadsLags(-lead_lag);
  }

  static void
@@ -389,7 +402,7 @@ DataTree::AddPossiblyNegativeConstant(double v)
  if (isnan(v))
    return NaN;
  if (isinf(v))
-    return (v < 0 ? MinusInfinity : Infinity);
+    return v < 0 ? MinusInfinity : Infinity;

  bool neg = false;
  if (v < 0)

--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -50,9 +50,11 @@ DynamicModel::copyHelper(const DynamicModel& m)

 DynamicModel::DynamicModel(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) :
-    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg, true},
    trend_component_model_table {trend_component_model_table_arg},
    var_model_table {var_model_table_arg}
 {
@@ -1041,14 +1043,8 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
         << (static_only_equations.size() > 0) << ";" << endl;

  // Say if model contains an external function call
-  bool has_external_function = false;
-  for (auto equation : equations)
-    if (equation->containsExternalFunction())
-      {
-        has_external_function = true;
-        break;
-      }
-  output << "M_.has_external_function = " << boolalpha << has_external_function << ';' << endl;
+  output << "M_.has_external_function = " << boolalpha
+         << ranges::any_of(equations, &ExprNode::containsExternalFunction) << ';' << endl;

  // Compute list of state variables, ordered in block-order
  vector<int> state_var;
@@ -1133,7 +1129,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
    output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
  output << "];" << endl;

-  writeDriverSparseIndicesHelper<true, false>(output);
+  writeDriverSparseIndicesHelper("dynamic", output);

  // Write LHS of each equation in text form
  output << "M_.lhs = {" << endl;
@@ -1357,8 +1353,7 @@ DynamicModel::fillVarModelTableFromOrigModel() const
          // save lhs variables
          equations[eqn]->arg1->collectVARLHSVariable(lhs);

-          equations[eqn]->arg1->countDiffs() > 0 ? diff_vec.push_back(true)
-                                                 : diff_vec.push_back(false);
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
          if (diff_vec.back())
            {
              set<pair<int, int>> diff_set;
@@ -1702,10 +1697,7 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
          // save lhs variables
          equations[eqn]->arg1->collectVARLHSVariable(lhs);

-          if (equations[eqn]->arg1->countDiffs() > 0)
-            diff_vec.push_back(true);
-          else
-            diff_vec.push_back(false);
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
          if (diff_vec.back())
            {
              set<pair<int, int>> diff_set;
@@ -1778,7 +1770,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
          exit(EXIT_FAILURE);
        }

-      if (diff.at(i) != true)
+      if (!diff.at(i))
        {
          cerr << "ERROR: the variable on the LHS of equation #" << eqn
               << " does not have the diff operator applied to it yet you are trying to undiff it."
@@ -1786,7 +1778,6 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
          exit(EXIT_FAILURE);
        }

-      bool printerr = false;
      expr_t node = nullptr;
      expr_t aux_var = lhs_expr_t.at(i);
      for (const auto& it : diff_subst_table)
@@ -1803,11 +1794,8 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
        }

      node = node->undiff();
-      auto it1 = diff_subst_table.find(node);
-      if (it1 == diff_subst_table.end())
-        printerr = true;

-      if (printerr)
+      if (auto it1 = diff_subst_table.find(node); it1 == diff_subst_table.end())
        { // we have undiffed something like diff(x), hence x is not in diff_subst_table
          lhs_expr_t.at(i) = node;
          lhs.at(i) = dynamic_cast<VariableNode*>(node)->symb_id;
@@ -1828,7 +1816,7 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
  for (auto& equation : equations)
    if (equation->containsPacExpectation(name))
      {
-        if (!pac_eq_name[name].empty())
+        if (pac_eq_name.contains(name))
          {
            cerr << "It is not possible to use 'pac_expectation(" << name
                 << ")' in several equations." << endl;
@@ -1927,6 +1915,13 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
                                   move(additive_vars_params_and_constants),
                                   move(optim_additive_vars_params_and_constants)};
      }
+
+  if (!pac_eq_name.contains(name))
+    {
+      cerr << "ERROR: the model does not contain the 'pac_expectation(" << name << ")' operator."
+           << endl;
+      exit(EXIT_FAILURE);
+    }
 }

 int
@@ -3697,6 +3692,36 @@ DynamicModel::substituteLogTransform()
    }
 }

+void
+DynamicModel::substituteAggregationOperators()
+{
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode*> neweqs;
+
+  for (auto& [symb_id, expr] : local_variables_table)
+    expr = expr->substituteAggregationOperators(subst_table, neweqs);
+
+  for (auto& equation : equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+
+  for (auto& equation : static_only_equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+
+  for (auto neweq : neweqs)
+    {
+      addEquation(neweq, nullopt);
+      addAuxEquation(neweq);
+    }
+}
+
 void
 DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
 {
@@ -3723,18 +3748,18 @@ void
 DynamicModel::detrendEquations()
 {
  // We go backwards in the list of trend_vars, to deal correctly with I(2) processes
-  for (const auto& it : std::ranges::reverse_view(nonstationary_symbols_map))
+  for (const auto& [symb_id, deflator] : std::ranges::reverse_view(nonstationary_symbols_map))
    {
      for (auto& equation : equations)
        {
          equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
          assert(equation);
        }
      for (auto& equation : static_only_equations)
        {
          equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
          assert(equation);
        }
    }

--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -73,7 +73,8 @@ public:
    void checkAllRegimesPresent() const noexcept(false);

  private:
-    pair<vector<string>, vector<string>> convertBitVectorToRegimes(const vector<bool>& r) const;
+    [[nodiscard]] pair<vector<string>, vector<string>>
+    convertBitVectorToRegimes(const vector<bool>& r) const;
  };

 private:
@@ -325,6 +326,7 @@ protected:
 public:
  DynamicModel(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);

@@ -566,6 +568,10 @@ public:
  // Performs the transformations associated to variables declared with “var(log)”
  void substituteLogTransform();

+  /* Performs the transformations associated to aggregation operators in heterogeneous models, such
+     as SUM(…) */
+  void substituteAggregationOperators();
+
  // Check that no variable was declared with “var(log)” in the given equations
  void checkNoWithLogTransform(const set<int>& eqnumbers);

@@ -708,7 +714,7 @@ public:
  {
    return tuple {static_only_equations, static_only_equations_lineno,
                  static_only_complementarity_conditions, static_only_equations_equation_tags};
-  };
+  }

  //! Returns true if a parameter was used in the model block with a lead or lag
  bool ParamUsedWithLeadLag() const;

--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
 /* -*- C++ -*- */
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -20,7 +20,6 @@


 %{
-#include <cstring>
 #include "ParsingDriver.hh"

 using namespace std;
@@ -60,7 +59,6 @@ string eofbuff;
 %x VERBATIM_BLOCK
 %x NATIVE
 %x NATIVE_COMMENT
-%x DATES_STATEMENT
 %x LINE1
 %x LINE2
 %x LINE3
@@ -70,7 +68,9 @@ string eofbuff;
 #define YY_USER_ACTION location_increment(yylloc, yytext);
 %}

-DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
+NAME [a-z_][a-z0-9_]*
+FLOAT_NUMBER ((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+)
+DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|[sh][12])

 %%
 /* Code put at the beginning of yylex() */
@@ -92,12 +92,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
                }

 /* spaces, tabs and carriage returns are ignored */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }

 /* Comments */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>%.*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"//".*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>%.*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"//".*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}

 <COMMENT>"*/"        {BEGIN comment_caller;}
 <COMMENT>.
@@ -111,6 +111,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>predetermined_variables {BEGIN DYNARE_STATEMENT; return token::PREDETERMINED_VARIABLES;}
 <INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;}
 <INITIAL>model_local_variable {BEGIN DYNARE_STATEMENT; return token::MODEL_LOCAL_VARIABLE;}
+<INITIAL>heterogeneity_dimension {BEGIN DYNARE_STATEMENT; return token::HETEROGENEITY_DIMENSION;}
 <INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
 <INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
 <INITIAL>set_time {BEGIN DYNARE_STATEMENT; return token::SET_TIME;}
@@ -235,6 +236,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>pac_target_info {BEGIN DYNARE_BLOCK; return token::PAC_TARGET_INFO;}
 <INITIAL>matched_irfs {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS;}
 <INITIAL>matched_irfs_weights {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS_WEIGHTS;}
+<INITIAL>perfect_foresight_controlled_paths {BEGIN DYNARE_BLOCK; return token::PERFECT_FORESIGHT_CONTROLLED_PATHS;}

 /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -245,7 +247,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>subsamples {return token::SUBSAMPLES;}
 <DYNARE_STATEMENT>options {return token::OPTIONS;}
 <DYNARE_STATEMENT>prior {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::PRIOR;
 }
 <INITIAL>std {BEGIN DYNARE_STATEMENT; return token::STD;}
@@ -255,31 +257,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>prior_function {BEGIN DYNARE_STATEMENT; return token::PRIOR_FUNCTION;}
 <INITIAL>posterior_function {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_FUNCTION;}

- /* Inside  of a Dynare statement */
-<DYNARE_STATEMENT>{DATE} {
-                           char* yycopy = strdup(yytext);
-                           char* uput = yycopy + yyleng;
-                           unput(')');
-                           unput('\'');
-                           while (uput > yycopy)
-                             unput(*--uput);
-                           unput('\'');
-                           unput('(');
-                           unput('s');
-                           unput('e');
-                           unput('t');
-                           unput('a');
-                           unput('d');
-                           free( yycopy );
-                         }
-<DYNARE_STATEMENT>${DATE} { yylloc->step();
-#if (YY_FLEX_MAJOR_VERSION > 2) || (YY_FLEX_MAJOR_VERSION == 2 && YY_FLEX_MINOR_VERSION >= 6)
-                            yyout << yytext + 1;
-#else
-                            *yyout << yytext + 1;
-#endif
-                          }
-<DYNARE_STATEMENT>dates  {dates_parens_nb=0; BEGIN DATES_STATEMENT; yylval->build<string>("dates");}
+<DYNARE_STATEMENT,DYNARE_BLOCK>{DATE} {
+  yylval->emplace<string>(yytext);
+  return token::DATE;
+}
+
+ /* Inside a Dynare statement */
 <DYNARE_STATEMENT>file                  {return token::FILE;}
 <DYNARE_STATEMENT>datafile 		{return token::DATAFILE;}
 <DYNARE_STATEMENT>dirname       {return token::DIRNAME;}
@@ -312,10 +295,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>posterior_nograph   		{return token::POSTERIOR_NOGRAPH;}
 <DYNARE_STATEMENT>nodisplay     {return token::NODISPLAY;}
 <DYNARE_STATEMENT>graph_format  {return token::GRAPH_FORMAT;}
-<DYNARE_STATEMENT>eps  {yylval->build<string>(yytext); return token::EPS;}
-<DYNARE_STATEMENT>pdf  {yylval->build<string>(yytext); return token::PDF;}
-<DYNARE_STATEMENT>fig  {yylval->build<string>(yytext); return token::FIG;}
-<DYNARE_STATEMENT>none  {yylval->build<string>(yytext); return token::NONE;}
+<DYNARE_STATEMENT>eps  {yylval->emplace<string>(yytext); return token::EPS;}
+<DYNARE_STATEMENT>pdf  {yylval->emplace<string>(yytext); return token::PDF;}
+<DYNARE_STATEMENT>fig  {yylval->emplace<string>(yytext); return token::FIG;}
+<DYNARE_STATEMENT>none  {yylval->emplace<string>(yytext); return token::NONE;}
 <DYNARE_STATEMENT>print   		{return token::PRINT;}
 <DYNARE_STATEMENT>noprint   		{return token::NOPRINT;}
 <DYNARE_STATEMENT>conf_sig  		{return token::CONF_SIG;}
@@ -423,6 +406,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>logarithmic_reduction {return token::LOGARITHMIC_REDUCTION;}
 <DYNARE_STATEMENT>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
 <DYNARE_STATEMENT>hybrid {return token::HYBRID;}
+<DYNARE_STATEMENT>use_first_order_solution {return token::USE_FIRST_ORDER_SOLUTION;}
 <DYNARE_STATEMENT>default {return token::DEFAULT;}
 <DYNARE_STATEMENT>init2shocks {return token::INIT2SHOCKS;}

@@ -435,6 +419,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>kitagawa {return token::KITAGAWA;}
 <DYNARE_STATEMENT>smooth {return token::SMOOTH;}
 <DYNARE_STATEMENT>stratified {return token::STRATIFIED;}
+<DYNARE_STATEMENT>residual {
+  yylval->emplace<string>(yytext);
+  return token::RESIDUAL;
+}
 <DYNARE_STATEMENT>cpf_weights {return token::CPF_WEIGHTS;}
 <DYNARE_STATEMENT>amisanotristani {return token::AMISANOTRISTANI;}
 <DYNARE_STATEMENT>murrayjonesparslow {return token::MURRAYJONESPARSLOW;}
@@ -456,43 +444,40 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>fsolve_options {return token::FSOLVE_OPTIONS;}

 <DYNARE_STATEMENT>alpha {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ALPHA;
 }
 <DYNARE_STATEMENT>beta {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::BETA;
 }
 <DYNARE_STATEMENT>gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma1 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA1;
 }
 <DYNARE_STATEMENT>inv_gamma2 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA2;
 }
 <DYNARE_STATEMENT>dirichlet {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DIRICHLET;
 }
-<DYNARE_STATEMENT>weibull {
-  yylval->build<string>(yytext);
-  return token::WEIBULL;
-}
+<DYNARE_STATEMENT>weibull {return token::WEIBULL;}
 <DYNARE_STATEMENT>normal {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NORMAL;
 }
 <DYNARE_STATEMENT>uniform {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::UNIFORM;
 }
 <DYNARE_STATEMENT>gsig2_lmdm {return token::GSIG2_LMDM;}
@@ -503,13 +488,13 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>ncsk {return token::NCSK;}
 <DYNARE_STATEMENT>nstd {return token::NSTD;}
 <DYNARE_STATEMENT>ninv {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NINV;
 }
 <DYNARE_STATEMENT>indxparr {return token::INDXPARR;}
 <DYNARE_STATEMENT>indxovr {return token::INDXOVR;}
 <DYNARE_STATEMENT>aband {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ABAND;
 }
 <DYNARE_STATEMENT>write_equation_tags {return token::WRITE_EQUATION_TAGS;}
@@ -526,18 +511,18 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>indxgdls {return token::INDXGDLS;}
 <DYNARE_STATEMENT>eq_ms {return token::EQ_MS;}
 <DYNARE_STATEMENT>cms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::CMS;
 }
 <DYNARE_STATEMENT>ncms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NCMS;
 }
 <DYNARE_STATEMENT>eq_cms {return token::EQ_CMS;}
 <DYNARE_STATEMENT>tlindx {return token::TLINDX;}
 <DYNARE_STATEMENT>tlnumber {return token::TLNUMBER;}
 <DYNARE_STATEMENT>cnum {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::CNUM;
 }
 <DYNARE_STATEMENT>nodecomposition {return token::NODECOMPOSITION;};
@@ -619,18 +604,9 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>substitute_libs {return token::SUBSTITUTE_LIBS;}
 <DYNARE_STATEMENT>compiler {return token::COMPILER;}
 <DYNARE_STATEMENT>instruments {return token::INSTRUMENTS;}
-<DYNARE_STATEMENT>hessian  {
-  yylval->build<string>(yytext);
-  return token::HESSIAN;
-}
-<DYNARE_STATEMENT>prior_variance  {
-  yylval->build<string>(yytext);
-  return token::PRIOR_VARIANCE;
-}
-<DYNARE_STATEMENT>identity_matrix  {
-  yylval->build<string>(yytext);
-  return token::IDENTITY_MATRIX;
-}
+<DYNARE_STATEMENT>hessian  {return token::HESSIAN;}
+<DYNARE_STATEMENT>prior_variance  {return token::PRIOR_VARIANCE;}
+<DYNARE_STATEMENT>identity_matrix  {return token::IDENTITY_MATRIX;}
 <DYNARE_STATEMENT>mcmc_jumping_covariance {return token::MCMC_JUMPING_COVARIANCE;}

 /* These four (var, varexo, varexo_det, parameters) are for change_type */
@@ -730,14 +706,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
 <DYNARE_STATEMENT>additional_optimizer_steps	{return token::ADDITIONAL_OPTIMIZER_STEPS;}
 <DYNARE_STATEMENT>bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; }
-<DYNARE_STATEMENT>optimal {
-  yylval->build<string>(yytext);
-  return token::OPTIMAL;
-}
-<DYNARE_STATEMENT>diagonal  {
-  yylval->build<string>(yytext);
-  return token::DIAGONAL;
-}
 <DYNARE_STATEMENT>gmm {return token::GMM;}
 <DYNARE_STATEMENT>smm {return token::SMM;}
 <DYNARE_STATEMENT>irf_matching {return token::IRF_MATCHING;}
@@ -783,8 +751,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>non_zero {return token::NON_ZERO;}

 <DYNARE_STATEMENT>\$[^$]*\$ {
-  strtok(yytext + 1, "$");
-  yylval->build<string>(yytext + 1);
+  yylval->emplace<string>(yytext + 1).pop_back();
  return token::TEX_NAME;
 }

@@ -797,20 +764,21 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>periods {return token::PERIODS;}
 <DYNARE_BLOCK>scales {return token::SCALES;}
 <DYNARE_BLOCK>add {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ADD;
 }
 <DYNARE_BLOCK>multiply {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::MULTIPLY;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>cutoff {return token::CUTOFF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>mfs {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::MFS;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>static_mfs {return token::STATIC_MFS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>balanced_growth_test_tol {return token::BALANCED_GROWTH_TEST_TOL;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>heterogeneity {return token::HETEROGENEITY;}
 <DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
 <DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}
 <DYNARE_BLOCK>normal_pdf {return token::NORMAL_PDF;}
@@ -822,19 +790,19 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>dsge_prior_weight {return token::DSGE_PRIOR_WEIGHT;}
 <DYNARE_BLOCK>surprise {return token::SURPRISE;}
 <DYNARE_BLOCK>bind {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
    return token::BIND;
 }
 <DYNARE_BLOCK>relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::RELAX;
 }
 <DYNARE_BLOCK>error_bind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ERROR_BIND;
 }
 <DYNARE_BLOCK>error_relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ERROR_RELAX;
 }
 <DYNARE_BLOCK>relative_to_initval {return token::RELATIVE_TO_INITVAL;}
@@ -848,22 +816,24 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK,DYNARE_STATEMENT>auxname {return token::AUXNAME;}
 <DYNARE_BLOCK>auxname_target_nonstationary {return token::AUXNAME_TARGET_NONSTATIONARY;}
 <DYNARE_BLOCK,DYNARE_STATEMENT>kind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::KIND;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>ll {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::LL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dl {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dd {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DD;
 }
 <DYNARE_BLOCK>weights {return token::WEIGHTS;}
+<DYNARE_BLOCK>exogenize {return token::EXOGENIZE;}
+<DYNARE_BLOCK>endogenize {return token::ENDOGENIZE;}

 /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}
@@ -875,7 +845,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>order {return token::ORDER;}
 <DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
 <DYNARE_STATEMENT>dr {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DR;
 }
 <DYNARE_STATEMENT>lyapunov_complex_threshold {return token::LYAPUNOV_COMPLEX_THRESHOLD;}
@@ -939,11 +909,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>time_shift {return token::TIME_SHIFT;}
 <DYNARE_STATEMENT>structural {return token::STRUCTURAL;}
 <DYNARE_STATEMENT>true {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::TRUE;
 }
 <DYNARE_STATEMENT>false {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
    return token::FALSE;
 }

@@ -1011,6 +981,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
 <DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
 <DYNARE_BLOCK>pac_target_nonstationary {return token::PAC_TARGET_NONSTATIONARY;}
+<DYNARE_BLOCK>sum {return token::SUM;}
 <DYNARE_STATEMENT>discount {return token::DISCOUNT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
@@ -1066,38 +1037,27 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>use_shock_groups {return token::USE_SHOCK_GROUPS;}
 <DYNARE_STATEMENT>colormap {return token::COLORMAP;}

-<DYNARE_STATEMENT,DYNARE_BLOCK>[a-z_][a-z0-9_]* {
-  yylval->build<string>(yytext);
+<DYNARE_STATEMENT,DYNARE_BLOCK>{NAME} {
+  yylval->emplace<string>(yytext);
  return token::NAME;
 }

-<DYNARE_STATEMENT,DYNARE_BLOCK>((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+) {
-  yylval->build<string>(yytext);
+<DYNARE_STATEMENT,DYNARE_BLOCK>{FLOAT_NUMBER} {
+  yylval->emplace<string>(yytext);
  return token::FLOAT_NUMBER;
 }

 <DYNARE_STATEMENT,DYNARE_BLOCK>[0-9]+ {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INT_NUMBER;
 }

-<DATES_STATEMENT>\( { yylval->as<string>().append(yytext); dates_parens_nb++; }
-<DATES_STATEMENT>\) {
-                      yylval->as<string>().append(yytext);
-                      if (--dates_parens_nb == 0)
-                      {
-                        BEGIN DYNARE_STATEMENT;
-                        return token::DATES;
-                      }
-                    }
-<DATES_STATEMENT>.  { yylval->as<string>().append(yytext); }
-
 <DYNARE_BLOCK>\|e { return token::PIPE_E; }
 <DYNARE_BLOCK>\|x { return token::PIPE_X; }
 <DYNARE_BLOCK>\|p { return token::PIPE_P; }

 <DYNARE_STATEMENT,DYNARE_BLOCK>\'[^\']*\' {
-  yylval->build<string>(yytext + 1).pop_back();
+  yylval->emplace<string>(yytext + 1).pop_back();
  return token::QUOTED_STRING;
 }

@@ -1129,11 +1089,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
    element in initval (in which case Dynare recognizes the matrix name as an external
    function symbol), and may want to modify the matrix later with Matlab statements.
 */
-<INITIAL>[a-z_][a-z0-9_]* {
+<INITIAL>{NAME} {
  if (driver.symbol_exists_and_is_not_modfile_local_or_external_function(yytext))
    {
      BEGIN DYNARE_STATEMENT;
-      yylval->build<string>(yytext);
+      yylval->emplace<string>(yytext);
      return token::NAME;
    }
  else
@@ -1154,7 +1114,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
    be able to back out of the statement if we realize it's a native statement
    and move to the NATIVE context
 */
-<INITIAL>\[([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*){1}\] {
+<INITIAL>\[([[:space:]]*{NAME}[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*{NAME}[[:space:]]*){1}\] {
  string yytextcpy{yytext};
  yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), '['), yytextcpy.end());
  yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), ']'), yytextcpy.end());
@@ -1178,7 +1138,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
  if (dynare_statement)
    {
      BEGIN DYNARE_STATEMENT;
-      yylval->build<vector<string>>(val);
+      yylval->emplace<vector<string>>(val);
      return token::SYMBOL_VEC;
    }
 }
@@ -1223,7 +1183,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <NATIVE_COMMENT>"*/"[[:space:]]*\n   { BEGIN NATIVE; }
 <NATIVE_COMMENT>.

-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }

 <*>.      { driver.error(*yylloc, "character unrecognized by lexer"); }
 %%

--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -283,7 +283,7 @@ protected:
  min_cost(bool is_matlab)
  {
    return is_matlab ? min_cost_matlab : min_cost_c;
-  };
+  }

  //! Initializes data member non_null_derivatives
  virtual void prepareForDerivation() = 0;
@@ -572,38 +572,35 @@ public:
  {
  };

-  //! Returns the maximum lead of endogenous in this expression
+  //! Returns the maximum lead of endogenous in this expression (not incl. heterogeneous endo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxEndoLead() const = 0;

-  //! Returns the maximum lead of exogenous in this expression
+  //! Returns the maximum lead of exogenous in this expression (not incl. heterogeneous exo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxExoLead() const = 0;

-  //! Returns the maximum lag of endogenous in this expression
+  //! Returns the maximum lag of endogenous in this expression (not incl. heterogeneous endo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxEndoLag() const = 0;

-  //! Returns the maximum lag of exogenous in this expression
+  //! Returns the maximum lag of exogenous in this expression (not incl. heterogeneous exo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxExoLag() const = 0;

-  //! Returns the maximum lead of endo/exo/exodet in this expression
-  /*! A negative value means that the expression contains only lagged
-    variables. A value of numeric_limits<int>::min() means that there is
-    no variable. */
+  /* Returns the maximum lead of endo/exo/exodet in this expression (including heterogeneous
+     endo/exo). A negative value means that the expression contains only lagged variables. A value
+     of numeric_limits<int>::min() means that there is no variable. */
  [[nodiscard]] virtual int maxLead() const = 0;

-  //! Returns the maximum lag of endo/exo/exodet in this expression
-  /*! A negative value means that the expression contains only leaded
-    variables. A value of numeric_limits<int>::min() means that there is
-    no variable. */
+  /* Returns the maximum lag of endo/exo/exodet in this expression (including heterogeneous
+     endo/exo). A negative value means that the expression contains only leaded variables. A value
+     of numeric_limits<int>::min() means that there is no variable. */
  [[nodiscard]] virtual int maxLag() const = 0;

-  //! Returns the maximum lag of endo/exo/exodet, as if diffs were expanded
-  /*! This function behaves as maxLag(), except that it treats diff()
-    differently. For e.g., on diff(diff(x(-1))), maxLag() returns 1 while
-    maxLagWithDiffsExpanded() returns 3. */
+  /* Returns the maximum lag of endo/exo/exodet (including heterogeneous endo/exo), as if diffs were
+     expanded. This function behaves as maxLag(), except that it treats diff() differently. For
+     e.g., on diff(diff(x(-1))), maxLag() returns 1 while maxLagWithDiffsExpanded() returns 3. */
  [[nodiscard]] virtual int maxLagWithDiffsExpanded() const = 0;

  [[nodiscard]] virtual expr_t undiff() const = 0;
@@ -946,7 +943,15 @@ public:
  /* Matches an expression that constitutes a complementarity condition.
     If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
     Otherwise, throws a MatchFailureException. */
-  [[nodiscard]] virtual tuple<int, expr_t, expr_t> matchComplementarityCondition() const;
+  [[nodiscard]] virtual tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension = nullopt) const;
+
+  /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
+     Also declares those aggregation operators in the HeterogeneityTable, so as to
+     compute their index in the dedicated vector in argument of the dynamic/static files. */
+  [[nodiscard]] virtual expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                              vector<BinaryOpNode*>& neweqs) const
+      = 0;
 };

 //! Object used to compare two nodes (using their indexes)
@@ -1058,6 +1063,8 @@ public:
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
 };

 //! Symbol or variable node
@@ -1165,6 +1172,8 @@ public:
                          vector<int>& powers) const override;
  [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
 };

 //! Unary operator node
@@ -1312,6 +1321,8 @@ public:
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  void decomposeAdditiveTerms(vector<pair<expr_t, int>>& terms, int current_sign) const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
 };

 //! Binary operator node
@@ -1507,7 +1518,11 @@ public:
                          vector<int>& powers) const override;
  [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
-  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
+  [[nodiscard]] tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension
+                                = nullopt) const override;
 };

 //! Trinary operator node
@@ -1651,6 +1666,8 @@ public:
  [[nodiscard]] bool containsPacTargetNonstationary(const string& pac_model_name
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
 };

@@ -1826,6 +1843,8 @@ public:
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
 };

 class ExternalFunctionNode : public AbstractExternalFunctionNode
@@ -2028,6 +2047,8 @@ public:
  expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
  [[nodiscard]] expr_t removeTrendLeadLag(const map<int, expr_t>& trend_symbols_map) const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;

 protected:
  void prepareForDerivation() override;

--- a/src/HeterogeneityTable.cc
+++ b/src/HeterogeneityTable.cc
+/*
+ * 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/>.
+ */
+
+#include <cassert>
+#include <utility>
+
+#include "HeterogeneityTable.hh"
+#include "SymbolTable.hh"
+
+void
+HeterogeneityTable::setSymbolTable(SymbolTable* symbol_table_arg)
+{
+  symbol_table = symbol_table_arg;
+}
+
+int
+HeterogeneityTable::addDimension(string name)
+{
+  if (name_to_id.contains(name))
+    throw AlreadyDeclaredDimensionException {move(name)};
+
+  int id {static_cast<int>(id_to_name.size())};
+  name_to_id.emplace(name, id);
+  id_to_name.push_back(move(name));
+  return id;
+}
+
+bool
+HeterogeneityTable::exists(const string& name) const
+{
+  return name_to_id.contains(name);
+}
+
+int
+HeterogeneityTable::getID(const string& name) const
+{
+  if (auto it = name_to_id.find(name); it != name_to_id.end())
+    return it->second;
+  else
+    throw UnknownDimensionNameException {name};
+}
+
+string
+HeterogeneityTable::getName(int id) const
+{
+  if (id < 0 || id >= static_cast<int>(id_to_name.size()))
+    throw UnknownDimensionIDException {id};
+  else
+    return id_to_name[id];
+}
+
+bool
+HeterogeneityTable::empty() const
+{
+  return name_to_id.empty();
+}
+
+vector<string>
+HeterogeneityTable::getDimensions() const
+{
+  return id_to_name;
+}
+
+int
+HeterogeneityTable::size() const
+{
+  return static_cast<int>(name_to_id.size());
+}
+
+void
+HeterogeneityTable::addSummedHeterogeneousEndogenous(int symb_id)
+{
+  assert(symbol_table->getType(symb_id) == SymbolType::heterogeneousEndogenous);
+  if (summed_het_endo_to_index.contains(symb_id))
+    throw AlreadyDeclaredSummedHeterogeneousEndogenousException {symb_id};
+
+  int index {static_cast<int>(index_to_summed_het_endo.size())};
+  summed_het_endo_to_index.emplace(symb_id, index);
+  index_to_summed_het_endo.push_back(symb_id);
+}
+
+int
+HeterogeneityTable::getSummedHeterogenousEndogenousIndex(int symb_id) const
+{
+  if (auto it = summed_het_endo_to_index.find(symb_id); it != summed_het_endo_to_index.end())
+    return it->second;
+  else
+    throw UnknownSummedHeterogeneousEndogenousException {symb_id};
+}
+
+int
+HeterogeneityTable::aggregateEndoSize() const
+{
+  return index_to_summed_het_endo.size();
+}
+
+void
+HeterogeneityTable::writeOutput(ostream& output) const
+{
+  for (size_t id {0}; id < id_to_name.size(); id++)
+    output << "M_.heterogeneity(" << id + 1 << ").dimension_name = '" << id_to_name[id] << "';"
+           << endl;
+
+  output << "M_.heterogeneity_aggregates = {" << endl;
+  for (int symb_id : index_to_summed_het_endo)
+    output << "'sum', " << symbol_table->getHeterogeneityDimension(symb_id) + 1 << ", "
+           << symbol_table->getTypeSpecificID(symb_id) + 1 << ";" << endl;
+  output << "};" << endl;
+}
+
+void
+HeterogeneityTable::writeJsonOutput(ostream& output) const
+{
+  assert(!empty());
+  output << R"("heterogeneity_dimension": [)";
+  for (bool first_written {false}; const auto& dim : id_to_name)
+    {
+      if (exchange(first_written, true))
+        output << ", ";
+      output << '"' << dim << '"';
+    }
+  output << "]" << endl;
+}
--- a/src/HeterogeneityTable.hh
+++ b/src/HeterogeneityTable.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 HETEROGENEITY_TABLE_HH
+#define HETEROGENEITY_TABLE_HH
+
+#include <map>
+#include <ostream>
+#include <string>
+#include <vector>
+
+using namespace std;
+
+class SymbolTable; // Forward declaration, to avoid circularity
+
+/*
+  There is a guarantee that heterogeneity IDs are increasing, i.e. if dimension A is added after
+  dimension B, then the ID of A is greater than the ID of B.
+  Moreover, the IDs form a contiguous interval starting at 0.
+*/
+class HeterogeneityTable
+{
+private:
+  // Maps dimension names to IDs
+  map<string, int> name_to_id;
+  // Maps dimension IDs to names
+  vector<string> id_to_name;
+
+  SymbolTable* symbol_table {nullptr}; // Cannot be a reference, because of circularity
+
+  /* Keeps track of the SUM() operator instances.
+     Maps a symbol ID that appears inside a SUM() operator into an index in
+     M_.heterogeneity_aggregates */
+  map<int, int> summed_het_endo_to_index;
+
+  // Maps an index in M_.heterogeneity_aggregates into a symbol ID
+  vector<int> index_to_summed_het_endo;
+
+public:
+  void setSymbolTable(SymbolTable* symbol_table_arg);
+
+  struct AlreadyDeclaredDimensionException
+  {
+    // Dimension name
+    const string name;
+  };
+  struct UnknownDimensionNameException
+  {
+    // Dimension name
+    const string name;
+  };
+  struct UnknownDimensionIDException
+  {
+    // Dimension ID
+    const int id;
+  };
+
+  // Returns the dimension ID
+  int addDimension(string name);
+  [[nodiscard]] bool exists(const string& name) const;
+  [[nodiscard]] int getID(const string& name) const;
+  [[nodiscard]] string getName(int id) const;
+  [[nodiscard]] bool empty() const;
+  [[nodiscard]] vector<string> getDimensions() const;
+  [[nodiscard]] int size() const;
+
+  struct AlreadyDeclaredSummedHeterogeneousEndogenousException
+  {
+    const int symb_id;
+  };
+  struct UnknownSummedHeterogeneousEndogenousException
+  {
+    const int symb_id;
+  };
+
+  void addSummedHeterogeneousEndogenous(int symb_id);
+  [[nodiscard]] int getSummedHeterogenousEndogenousIndex(int symb_id) const;
+  [[nodiscard]] int aggregateEndoSize() const;
+
+  void writeOutput(ostream& output) const;
+  void writeJsonOutput(ostream& output) const;
+};
+
+#endif
--- a/src/HeterogeneousModel.cc
+++ b/src/HeterogeneousModel.cc
+/*
+ * Copyright © 2024-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/>.
+ */
+
+#include <cstdlib>
+#include <iostream>
+
+#include "HeterogeneousModel.hh"
+
+HeterogeneousModel::HeterogeneousModel(SymbolTable& symbol_table_arg,
+                                       NumericalConstants& num_constants_arg,
+                                       ExternalFunctionsTable& external_functions_table_arg,
+                                       HeterogeneityTable& heterogeneity_table_arg,
+                                       int heterogeneity_dimension_arg) :
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg, true},
+    heterogeneity_dimension {heterogeneity_dimension_arg}
+{
+}
+
+HeterogeneousModel&
+HeterogeneousModel::operator=(const HeterogeneousModel& m)
+{
+  ModelTree::operator=(m);
+
+  assert(heterogeneity_dimension == m.heterogeneity_dimension);
+
+  deriv_id_table = m.deriv_id_table;
+  inv_deriv_id_table = m.inv_deriv_id_table;
+
+  return *this;
+}
+
+void
+HeterogeneousModel::computeChainRuleJacobian()
+{
+  cerr << "Heterogeneous::computeChainRuleJacobian(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getBlockJacobianEndoCol([[maybe_unused]] int blk, [[maybe_unused]] int var,
+                                            [[maybe_unused]] int lead_lag) const
+{
+  cerr << "Heterogeneous::getBlockJacobianEndoCol(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getMFS() const
+{
+  cerr << "Heterogeneous::getMFS(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+void
+HeterogeneousModel::computeDerivIDs()
+{
+  set<pair<int, int>> dynvars;
+
+  for (auto& equation : equations)
+    {
+      equation->collectDynamicVariables(SymbolType::heterogeneousEndogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::heterogeneousExogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::parameter, dynvars);
+    }
+
+  for (const auto& [symb_id, lead_lag] : dynvars)
+    {
+      auto type {symbol_table.getType(symb_id)};
+      if (isHeterogeneous(type))
+        assert(symbol_table.getHeterogeneityDimension(symb_id) == heterogeneity_dimension);
+      if (type == SymbolType::heterogeneousEndogenous || type == SymbolType::endogenous)
+        assert(abs(lead_lag) <= 1);
+      if (type == SymbolType::heterogeneousExogenous || type == SymbolType::exogenous)
+        assert(lead_lag == 0);
+      int deriv_id {static_cast<int>(deriv_id_table.size())};
+      deriv_id_table.emplace(pair {symb_id, lead_lag}, deriv_id);
+      inv_deriv_id_table.emplace_back(symb_id, lead_lag);
+    }
+}
+
+void
+HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll)
+{
+  assert(!use_dll); // Not yet implemented
+
+  computeDerivIDs();
+
+  set<int> vars;
+  for (auto& [symb_lag, deriv_id] : deriv_id_table)
+    if (symbol_table.getType(symb_lag.first) != SymbolType::parameter)
+      vars.insert(deriv_id);
+
+  cout << "Computing " << modelClassName() << " derivatives (order " << derivsOrder << ")." << endl;
+
+  computeDerivatives(derivsOrder, vars);
+
+  computeTemporaryTerms(!use_dll, no_tmp_terms);
+
+  computeMCPEquationsReordering(heterogeneity_dimension);
+}
+
+void
+HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
+{
+  assert(!julia); // Not yet implemented
+  writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+  writeComplementarityConditionsFile<true>(basename, heterogeneity_dimension);
+}
+
+int
+HeterogeneousModel::getJacobianCol(int deriv_id, bool sparse) const
+{
+  assert(sparse);
+
+  SymbolType type {getTypeByDerivID(deriv_id)};
+  int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+  int lag {getLagByDerivID(deriv_id)};
+
+  if (type == SymbolType::heterogeneousEndogenous)
+    return tsid + (lag + 1) * symbol_table.het_endo_nbr(heterogeneity_dimension);
+
+  int shift {3 * symbol_table.het_endo_nbr(heterogeneity_dimension)};
+
+  if (type == SymbolType::heterogeneousExogenous)
+    return shift + tsid;
+
+  shift += symbol_table.het_exo_nbr(heterogeneity_dimension);
+
+  if (type == SymbolType::endogenous)
+    return shift + tsid + (lag + 1) * symbol_table.endo_nbr();
+
+  shift += symbol_table.endo_nbr();
+
+  if (type == SymbolType::exogenous)
+    return shift + tsid;
+
+  throw UnknownDerivIDException();
+}
+
+int
+HeterogeneousModel::getJacobianColsNbr(bool sparse) const
+{
+  assert(sparse);
+  return 3 * (symbol_table.het_endo_nbr(heterogeneity_dimension) + symbol_table.endo_nbr())
+         + symbol_table.het_exo_nbr(heterogeneity_dimension) + symbol_table.exo_nbr();
+}
+
+SymbolType
+HeterogeneousModel::getTypeByDerivID(int deriv_id) const noexcept(false)
+{
+  return symbol_table.getType(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getLagByDerivID(int deriv_id) const noexcept(false)
+{
+  if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].second;
+}
+
+int
+HeterogeneousModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
+{
+  if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].first;
+}
+
+int
+HeterogeneousModel::getTypeSpecificIDByDerivID(int deriv_id) const
+{
+  return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getDerivID(int symb_id, int lead_lag) const noexcept(false)
+{
+  if (auto it = deriv_id_table.find({symb_id, lead_lag}); it == deriv_id_table.end())
+    throw UnknownDerivIDException();
+  else
+    return it->second;
+}
+
+void
+HeterogeneousModel::writeDriverOutput(ostream& output) const
+{
+  std::vector<int> state_var;
+  for (int endoID = 0; endoID < symbol_table.het_endo_nbr(heterogeneity_dimension); endoID++)
+    try
+      {
+        getDerivID(symbol_table.getID(SymbolType::heterogeneousEndogenous, endoID,
+                                      heterogeneity_dimension),
+                   -1);
+        if (ranges::find(state_var, endoID) == state_var.end())
+          state_var.push_back(endoID);
+      }
+    catch (UnknownDerivIDException& e)
+      {
+      }
+
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").state_var = [";
+  for (int it : state_var)
+    output << it + 1 << " ";
+  output << "];" << endl;
+
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").dynamic_tmp_nbr = [";
+  for (const auto& it : temporary_terms_derivatives)
+    output << it.size() << "; ";
+  output << "];" << endl;
+  writeDriverSparseIndicesHelper(
+      "heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1
+         << ").dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
+}
--- a/src/HeterogeneousModel.hh
+++ b/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) = default;
+  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
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
 /*
- * Copyright © 2006-2024 Dynare Team
+ * Copyright © 2006-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -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
@@ -396,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
@@ -529,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; };
@@ -550,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
@@ -702,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
@@ -822,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);

@@ -864,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)
    {
@@ -904,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;
@@ -939,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);
@@ -950,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() << ", "
@@ -971,7 +1137,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
              << "M_.heteroskedastic_shocks.Qvalue_orig = [];" << endl
              << "M_.heteroskedastic_shocks.Qscale_orig = [];" << endl
              << "M_.matched_irfs = {};" << endl
-              << "M_.matched_irfs_weights = {};" << endl;
+              << "M_.matched_irfs_weights = {};" << endl
+              << "M_.perfect_foresight_controlled_paths = [];" << endl;

  // NB: options_.{ramsey,discretionary}_policy should rather be fields of M_
  mOutputFile << boolalpha << "options_.linear = " << linear << ";" << endl
@@ -1038,6 +1205,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)
      {
@@ -1142,8 +1312,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,

      if (!no_warn)
        {
-          if (warnings.countWarnings() > 0)
-            mOutputFile << "disp('Note: " << warnings.countWarnings()
+          if (int num_warnings {warnings.numWarnings()}; num_warnings > 0)
+            mOutputFile << "disp('Note: " << num_warnings
                        << " warning(s) encountered in the preprocessor')" << endl;

          mOutputFile << "if ~isempty(lastwarn)" << endl
@@ -1186,6 +1356,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);
    }
 }

@@ -1263,6 +1436,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);
No results found