src/DynareMain.cc

@@ -56,7 +56,7 @@ usage()
           "[conffile=path_to_config_file] [parallel_follower_open_mode] "
           "[parallel_test] [parallel_use_psexec=true|false]"
        << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] "
-          "[compute_xrefs] [output=second|third] [language=matlab|julia]"
+          "[compute_xrefs] [output=first|second|third] [language=matlab|julia]"
        << " [params_derivs_order=0|1|2] [transform_unary_ops] "
           "[exclude_eqs=<equation_tag_list_or_file>] [include_eqs=<equation_tag_list_or_file>]"
        << " [json=parse|check|transform|compute] [jsonstdout] [onlyjson] [jsonderivsimple] "
@@ -320,7 +320,9 @@ main(int argc, char** argv)
 
           s.erase(0, 7);
 
-          if (s == "second")
+          if (s == "first")
+            output_mode = OutputType::first;
+          else if (s == "second")
             output_mode = OutputType::second;
           else if (s == "third")
             output_mode = OutputType::third;

src/EquationTags.cc

 /*
- * Copyright © 2020-2023 Dynare Team
+ * Copyright © 2020-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -74,10 +74,9 @@ EquationTags::erase(const set<int>& eqns, const map<int, int>& old_eqn_num_2_new
     eqn_tags.erase(eqn);
 
   for (const auto& [oldeqn, neweqn] : old_eqn_num_2_new)
-    for (auto& [eqn, tags] : eqn_tags)
-      if (eqn == oldeqn)
+    if (eqn_tags.contains(oldeqn))
       {
-          auto tmp = eqn_tags.extract(eqn);
+        auto tmp = eqn_tags.extract(oldeqn);
         tmp.key() = neweqn;
         eqn_tags.insert(move(tmp));
       }

src/ExprNode.hh

 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,6 +20,7 @@
 #ifndef EXPR_NODE_HH
 #define EXPR_NODE_HH
 
+#include <concepts>
 #include <functional>
 #include <map>
 #include <optional>
@@ -282,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;
@@ -571,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;
@@ -941,6 +939,19 @@ public:
 
   // Substitutes orig_symb_id(±l) with exp(aux_symb_id(±l)) (used for “var(log)”)
   [[nodiscard]] virtual expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const = 0;
+
+  /* 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 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)
@@ -1052,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
@@ -1084,6 +1097,10 @@ protected:
 
 public:
   VariableNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg, int lag_arg);
+  [[nodiscard]] SymbolType get_type() const;
+  [[nodiscard]] string getName() const;
+  [[nodiscard]] int getDerivID() const;
+  [[nodiscard]] int getTypeSpecificID() const;
   void writeOutput(ostream& output, ExprNodeOutputType output_type,
                    const temporary_terms_t& temporary_terms,
                    const temporary_terms_idxs_t& temporary_terms_idxs,
@@ -1101,7 +1118,6 @@ public:
                            const deriv_node_temp_terms_t& tef_terms) const override;
   expr_t toStatic(DataTree& static_datatree) const override;
   void computeXrefs(EquationInfo& ei) const override;
-  [[nodiscard]] SymbolType get_type() const;
   BinaryOpNode* normalizeEquationHelper(const set<expr_t>& contain_var, expr_t rhs) const override;
   [[nodiscard]] int maxEndoLead() const override;
   [[nodiscard]] int maxExoLead() const override;
@@ -1156,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
@@ -1173,6 +1191,7 @@ protected:
   // Returns the node obtained by applying a transformation recursively on the argument (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1302,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
@@ -1339,6 +1360,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1451,7 +1473,7 @@ public:
   [[nodiscard]] expr_t unpackPowerDeriv() const;
   //! Returns MULT_i*(lhs-rhs) = 0, creating multiplier MULT_i
   expr_t addMultipliersToConstraints(int i);
-  //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
+  //! Returns the non-zero hand side of an equation (that must have a hand side equal to zero)
   [[nodiscard]] expr_t getNonZeroPartofEquation() const;
   [[nodiscard]] bool isInStaticForm() const override;
   void fillAutoregressiveRow(int eqn, const vector<int>& lhs,
@@ -1496,6 +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]] 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
@@ -1532,6 +1559,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1638,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;
 };
 
@@ -1660,6 +1690,7 @@ private:
   // Returns the node obtained by applying a transformation recursively on the arguments (in same
   // datatree)
   template<typename Callable, typename... Args>
+    requires invocable<Callable, expr_t, Args...>
   expr_t
   recurseTransform(Callable&& op, Args&&... args) const
   {
@@ -1709,6 +1740,7 @@ protected:
 public:
   AbstractExternalFunctionNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg,
                                vector<expr_t> arguments_arg);
+  [[nodiscard]] string getName() const;
   void computeTemporaryTerms(const pair<int, int>& derivOrder,
                              map<pair<int, int>, unordered_set<expr_t>>& temp_terms_map,
                              unordered_map<expr_t, pair<int, pair<int, int>>>& reference_count,
@@ -1811,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
@@ -2013,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;

src/ExtendedPreprocessorTypes.hh

@@ -24,6 +24,7 @@
 enum class OutputType
 {
   standard, // Default value, infer the derivation order from .mod file only
+  first,    // Output only 1st dynamic derivatives with no other computations
   second,   // Output at least 2nd dynamic derivatives
   third,    // Output at least 3rd dynamic derivatives
 };

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;
+}

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

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);
+    }
+}
+
+/*
+ * Unfold complementarity conditions: (i) declare the multipliers associated
+ * with each bound constraint μ_l and μ_u ; (ii) add or substract the
+ * multiplier into the associated condition; (iii) add the the complementarity
+ * slackness conditions into the set of equations. For example,
+ * households choose {cₜ, aₜ₊₁} to maximize expected lifetime utility:
+ *     max 𝐸ₜ [∑ₛ₌₀^∞ βˢ · u(cₜ₊ₛ)]
+ *
+ * Subject to:
+ *   1. Budget constraint:      cₜ + aₜ₊₁ = yₜ + (1 + rₜ) · aₜ
+ *   2. Borrowing constraint:   aₜ₊₁ ≥ aₘᵢₙ
+ *
+ * Let u'(cₜ) denote the marginal utility of consumption.
+ * Let μₜ ≥ 0 be the Lagrange multiplier on the borrowing constraint.
+ *
+ * Then, the Euler equation becomes:
+ *     u′(cₜ) = β · (1 + rₜ₊₁) · u′(cₜ₊₁) − μₜ
+ *
+ * Together with:
+ *     aₜ₊₁ ≥ aₘᵢₙ                 [primal feasibility]
+ *     μₜ ≥ 0                      [dual feasibility]
+ *     μₜ · (aₜ₊₁ − aₘᵢₙ) = 0      [complementarity slackness]
+ * Note that the primal feasibility and dual feasibility constraints are not
+ * introduced here, but Bhandari et al. (2023) show in Appendix B.1 that they
+ * are redundant.
+ */
+void
+HeterogeneousModel::transformPass()
+{
+  for (int i = 0; i < static_cast<int>(equations.size()); ++i)
+    {
+      if (!complementarity_conditions[i])
+        continue;
+
+      /*
+       * `const auto& [symb_id, lb, ub] = *complementarity_conditions[i];` was not used here because
+       * the call to `addEquation` may eventually lead to a resize of the
+       * `complementarity_conditions` vector, which may invalidate the reference to its element. We
+       * take a copy instead for safety.
+       */
+      auto [symb_id, lb, ub] = *complementarity_conditions[i];
+
+      VariableNode* var = getVariable(symb_id);
+      if (lb)
+        {
+          int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
+              heterogeneity_dimension, i, "MULT_L_" + symbol_table.getName(symb_id));
+          expr_t mu_L = AddVariable(mu_id);
+          auto substeq = AddEqual(AddPlus(equations[i]->arg1, mu_L), equations[i]->arg2);
+          assert(substeq);
+          equations[i] = substeq;
+          addEquation(AddEqual(AddTimes(mu_L, AddMinus(var, lb)), Zero), nullopt);
+        }
+      if (ub)
+        {
+          int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
+              heterogeneity_dimension, i, "MULT_U_" + symbol_table.getName(symb_id));
+          auto mu_U = AddVariable(mu_id);
+          auto substeq = AddEqual(AddMinus(equations[i]->arg1, mu_U), equations[i]->arg2);
+          assert(substeq);
+          equations[i] = substeq;
+          addEquation(AddEqual(AddTimes(mu_U, AddMinus(ub, var)), Zero), nullopt);
+        }
+    }
+}
+
+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;
+}

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 transformPass();
+
+  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

 /*
- * Copyright © 2006-2023 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},
-    static_model {symbol_table, num_constants, external_functions_table},
-    steady_state_model {symbol_table, num_constants, external_functions_table, static_model},
+    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, 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
@@ -167,7 +190,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
       exit(EXIT_FAILURE);
     }
 
-  if (mod_file_struct.ramsey_constraints_present && !mod_file_struct.ramsey_model_present)
+  if (!ramsey_constraints.empty() && !mod_file_struct.ramsey_model_present)
     {
       cerr << "ERROR: A ramsey_constraints block requires the presence of a ramsey_model or "
               "ramsey_policy statement"
@@ -279,6 +302,10 @@ ModFile::checkPass(bool nostrict, bool stochastic)
         }
     }
 
+  /* This check must come before checking that there are as many static-only as dynamic-only
+     equations, see #103 */
+  dynamic_model.checkOccbinRegimes();
+
   if (dynamic_model.staticOnlyEquationsNbr() != dynamic_model.dynamicOnlyEquationsNbr())
     {
       cerr << "ERROR: the number of equations marked [static] must be equal to the number of "
@@ -351,10 +378,8 @@ 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(),
+  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)
     {
@@ -374,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;
@@ -394,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
@@ -527,12 +666,22 @@ 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,
+      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; };
+      map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints;
+      for (const auto& [symb_id, bounds] : ramsey_constraints)
+        cloned_ramsey_constraints.try_emplace(symb_id, clone_if_not_null(bounds.first),
+                                              clone_if_not_null(bounds.second));
       mod_file_struct.ramsey_orig_endo_nbr
-          = ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective);
+          = ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective,
+                                                                       cloned_ramsey_constraints);
       ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
     }
 
@@ -541,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
@@ -579,6 +731,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
 
   dynamic_model.reorderAuxiliaryEquations();
 
+  symbol_table.resizeHetAuxVars();
+  for (auto& hm : heterogeneous_models)
+    hm.transformPass();
+
   // Freeze the symbol table
   symbol_table.freeze();
 
@@ -693,6 +849,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
@@ -798,8 +964,22 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
                                                         no_tmp_terms, block, use_dll);
             }
         }
-      else // No computing task requested, compute derivatives up to 2nd order by default
+      else // No computing task requested, compute derivatives up to 2nd order by default unless
+           // output=first (preprocessor#100) or third (preprocessor#121) is requested
+        {
+          switch (output)
+            {
+            case OutputType::first:
+              dynamic_model.computingPass(1, 0, global_eval_context, no_tmp_terms, block, use_dll);
+              break;
+            case OutputType::third:
+              dynamic_model.computingPass(3, 0, global_eval_context, no_tmp_terms, block, use_dll);
+              break;
+            default:
               dynamic_model.computingPass(2, 0, global_eval_context, no_tmp_terms, block, use_dll);
+              break;
+            }
+        }
 
       if (linear)
         {
@@ -813,6 +993,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);
 
@@ -855,6 +1038,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)
     {
@@ -895,8 +1083,8 @@ 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;"
+    mOutputFile << "clearvars -global M_ options_ oo_ estim_params_ bayestopt_ dataset_ "
+                   "dataset_info estimation_info;"
                 << endl;
 
   if (!notime)
@@ -930,6 +1118,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);
@@ -941,6 +1130,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() << ", "
@@ -962,7 +1155,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
@@ -1029,6 +1223,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)
       {
@@ -1133,8 +1330,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
@@ -1177,6 +1374,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);
     }
 }
 
@@ -1254,6 +1454,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

 /*
- * Copyright © 2006-2023 Dynare Team
+ * Copyright © 2006-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -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};
@@ -117,6 +123,11 @@ public:
   /*! (i.e. option parallel_local_files of model block) */
   vector<string> parallel_local_files;
 
+  /* Contents of the ramsey_constraints block.
+     Maps symb_id → (lower_bound, upper_bound).
+     NB: The two expr_t live in dynamic_model */
+  map<int, pair<expr_t, expr_t>> ramsey_constraints;
+
 private:
   //! List of statements
   vector<unique_ptr<Statement>> statements;

src/ModelEquationBlock.cc

 /*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -26,9 +26,21 @@
 
 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}
+{
+}
+
+void
+PlannerObjective::writeDriverOutput(ostream& output) const
 {
+  output << "M_.objective_tmp_nbr = [";
+  for (const auto& it : temporary_terms_derivatives)
+    output << it.size() << "; ";
+  output << "];" << endl;
+  writeDriverSparseIndicesHelper("objective", output);
 }
 
 void
@@ -41,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}
 {
 }
 
@@ -57,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}
 {
 }
@@ -317,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}
 {
 }
 
@@ -367,7 +391,7 @@ Epilogue::checkPass(ModFileStructure& mod_file_struct) const
   for (const auto& [symb_id, expr] : dynamic_def_table)
     if (so_far_defined.contains(symb_id))
       {
-        cerr << "WARNING: in the 'epilogue' block, variable '" << symbol_table.getName(symb_id)
+        cerr << "ERROR: in the 'epilogue' block, variable '" << symbol_table.getName(symb_id)
              << "' is declared twice" << endl;
         exit(EXIT_FAILURE);
       }
@@ -386,10 +410,10 @@ void
 Epilogue::detrend(const map<int, expr_t>& trend_symbols_map,
                   const nonstationary_symbols_map_t& nonstationary_symbols_map)
 {
-  for (const auto& it : ranges::reverse_view(nonstationary_symbols_map))
+  for (const auto& [symb_id, deflator] : ranges::reverse_view(nonstationary_symbols_map))
     for (auto& [symb_id, expr] : dynamic_def_table)
       {
-        expr = expr->detrend(it.first, it.second.first, it.second.second);
+        expr = expr->detrend(symb_id, deflator.first, deflator.second);
         assert(expr);
       }
 

src/ModelEquationBlock.hh

 /*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -30,7 +30,10 @@ 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;
 
 protected:
   string
@@ -48,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);
@@ -58,6 +62,17 @@ protected:
   {
     return "original Ramsey model";
   }
+  int
+  getJacobianCol(int deriv_id, [[maybe_unused]] bool sparse) const override
+  {
+    /* Override the DynamicModel method by returning a dummy Jacobian column number.
+       The override is necessary because the method from DynamicModel fails with
+       endos with lag/lead greater than 1 or exos with a lag/lead, while substitutions
+       are by definition not done for an original model.
+       In particular, this fixes dynare#1960 (equation derivatives are computed for models declared
+       as linear, to check whether they are truly linear). */
+    return deriv_id;
+  }
 };
 
 class SteadyStateModel : public DataTree
@@ -73,6 +88,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);
@@ -106,6 +122,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);