src/Shocks.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -22,9 +22,11 @@
 
 #include <map>
 #include <string>
+#include <variant>
 #include <vector>
 
 #include "ExprNode.hh"
+#include "HeterogeneityTable.hh"
 #include "Statement.hh"
 #include "SymbolTable.hh"
 
@@ -33,8 +35,10 @@ using namespace std;
 class AbstractShocksStatement : public Statement
 {
 public:
-  // The tuple is (period1, period2, value)
-  using det_shocks_t = map<int, vector<tuple<int, int, expr_t>>>;
+  // A period range is either two indices (1-based), or two dates (from dseries)
+  using period_range_t = variant<pair<int, int>, pair<string, string>>;
+  // The pair is (period range, value)
+  using det_shocks_t = map<int, vector<pair<period_range_t, expr_t>>>;
   enum class ShockType
   {
     level,               // The value is the level of the exogenous (“values” statement in “shocks”)
@@ -122,7 +126,7 @@ public:
 class ShocksLearntInStatement : public Statement
 {
 public:
-  const int learnt_in_period;
+  const variant<int, string> learnt_in_period;
   //! Does this “shocks(learnt_in=…)” or “mshocks(learnt_in=…)” block replace the previous ones?
   const bool overwrite;
   enum class LearntShockType
@@ -139,8 +143,9 @@ public:
                                // state as anticipated in the same informational period (“values”
                                // statement in “mshocks(learnt_in=…, relative_to_initval)”)
   };
-  // The tuple is (type, period1, period2, value)
-  using learnt_shocks_t = map<int, vector<tuple<LearntShockType, int, int, expr_t>>>;
+  // The tuple is (type, period range, value)
+  using learnt_shocks_t
+      = map<int, vector<tuple<LearntShockType, AbstractShocksStatement::period_range_t, expr_t>>>;
   const learnt_shocks_t learnt_shocks;
 
 private:
@@ -148,13 +153,54 @@ private:
   static string typeToString(LearntShockType type);
 
 public:
-  ShocksLearntInStatement(int learnt_in_period_arg, bool overwrite_arg,
+  ShocksLearntInStatement(variant<int, string> learnt_in_period_arg, bool overwrite_arg,
                           learnt_shocks_t learnt_shocks_arg, const SymbolTable& symbol_table_arg);
   void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
   void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream& output) const override;
 };
 
+class HeterogeneousShocksStatement : public Statement
+{
+public:
+  const int heterogeneity_dimension;
+  const bool overwrite;
+
+  using var_and_std_shocks_t = map<int, expr_t>;
+  using covar_and_corr_shocks_t = map<pair<int, int>, expr_t>;
+
+  const var_and_std_shocks_t var_shocks, std_shocks;
+  const covar_and_corr_shocks_t covar_shocks, corr_shocks;
+
+private:
+  const SymbolTable& symbol_table;
+  const HeterogeneityTable& heterogeneity_table;
+
+  void writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it, bool stddev) const;
+  void writeVarAndStdShocks(ostream& output) const;
+  void writeCovarOrCorrShock(ostream& output, const pair<pair<int, int>, expr_t>& it,
+                             bool corr) const;
+  void writeCovarAndCorrShocks(ostream& output) const;
+
+  string
+  sigmaeName() const
+  {
+    return "M_.heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").Sigma_e"s;
+  }
+
+public:
+  HeterogeneousShocksStatement(int heterogeneity_dimension_arg, bool overwrite_arg,
+                               var_and_std_shocks_t var_shocks_arg,
+                               var_and_std_shocks_t std_shocks_arg,
+                               covar_and_corr_shocks_t covar_shocks_arg,
+                               covar_and_corr_shocks_t corr_shocks_arg,
+                               const SymbolTable& symbol_table_arg,
+                               const HeterogeneityTable& heterogeneity_table_arg);
+  void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
+  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream& output) const override;
+};
+
 class ConditionalForecastPathsStatement : public Statement
 {
 private:
@@ -170,6 +216,26 @@ public:
   static int computePathLength(const AbstractShocksStatement::det_shocks_t& paths);
 };
 
+class PerfectForesightControlledPathsStatement : public Statement
+{
+public:
+  // (exogenize_id, vector of (period range, value), endogenize_id)
+  using paths_t
+      = vector<tuple<int, vector<pair<AbstractShocksStatement::period_range_t, expr_t>>, int>>;
+
+private:
+  const paths_t paths;
+  const variant<int, string> learnt_in_period;
+  const SymbolTable& symbol_table;
+
+public:
+  PerfectForesightControlledPathsStatement(paths_t paths_arg,
+                                           variant<int, string> learnt_in_period_arg,
+                                           const SymbolTable& symbol_table_arg);
+  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream& output) const override;
+};
+
 class MomentCalibration : public Statement
 {
 public:

src/Statement.cc

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

src/Statement.hh

 /*
- * Copyright © 2006-2023 Dynare Team
+ * Copyright © 2006-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -30,6 +30,10 @@
 #include "SymbolList.hh"
 #include "WarningConsolidation.hh"
 
+// Helper constant for visitors
+template<class>
+inline constexpr bool always_false_v {false};
+
 struct ModFileStructure
 {
   //! Whether check is present
@@ -169,8 +173,6 @@ struct ModFileStructure
   bool endval_learnt_in_present {false};
   // Whether an occbin_constraints block appears
   bool occbin_constraints_present {false};
-  // Whether a ramsey_constraints block appears
-  bool ramsey_constraints_present {false};
 };
 
 class Statement
@@ -254,6 +256,12 @@ public:
   {
   };
 
+  // pair<string, string> corresponds to a pair of numerical values
+  // vector<vector<string>> corresponds to a vector of vectors of numerical values
+  using OptionValue
+      = variant<NumVal, pair<string, string>, StringVal, DateVal, SymbolListVal, vector<int>,
+                VecStrVal, VecCellStrVal, VecValueVal, vector<vector<string>>>;
+
   [[nodiscard]] bool empty() const;
   void clear();
   // Whether there is an option with that name that has been given a value
@@ -281,7 +289,8 @@ public:
   /* Retrieves the value of the option with that name.
      Throws UnknownOptionException if there is no option with that name.
      Throws bad_variant_access if the option has a value of a different type. */
-  template<class T>
+  template<typename T>
+    requires requires(T p) { std::get<T>(OptionValue {}); }
   T
   get(const string& name) const
   {
@@ -295,7 +304,8 @@ public:
   /* Retrieves the value of the option with that name.
      Returns nullopt if there is no option with that name.
      Throws bad_variant_access if the option has a value of a different type. */
-  template<class T>
+  template<typename T>
+    requires requires(T p) { std::get<T>(OptionValue {}); }
   optional<T>
   get_if(const string& name) const
   {
@@ -309,6 +319,7 @@ public:
   /* Applies a variant visitor to the value of the option with that name.
      Throws UnknownOptionException if there is no option with that name. */
   template<class Visitor>
+    requires invocable<Visitor, OptionValue>
   decltype(auto)
   visit(const string& name, Visitor&& vis) const
   {
@@ -324,15 +335,8 @@ public:
   void writeJsonOutput(ostream& output) const;
 
 private:
-  // pair<string, string> corresponds to a pair of numerical values
-  // vector<vector<string>> corresponds to a vector of vectors of numerical values
-  map<string, variant<NumVal, pair<string, string>, StringVal, DateVal, SymbolListVal, vector<int>,
-                      VecStrVal, VecCellStrVal, VecValueVal, vector<vector<string>>>>
-      options;
+  map<string, OptionValue> options;
   void writeOutputCommon(ostream& output, const string& option_group) const;
-  // Helper constant for visitors
-  template<class>
-  static constexpr bool always_false_v {false};
 };
 
 #endif

src/StaticModel.cc

 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -23,6 +23,7 @@
 #include <cstdlib>
 #include <iostream>
 #include <numeric>
+#include <ranges>
 #include <sstream>
 #include <unordered_map>
 
@@ -30,8 +31,10 @@
 #include "StaticModel.hh"
 
 StaticModel::StaticModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                         ExternalFunctionsTable& external_functions_table_arg) :
-    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg}
+                         ExternalFunctionsTable& external_functions_table_arg,
+                         HeterogeneityTable& heterogeneity_table_arg) :
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg}
 {
 }
 
@@ -71,7 +74,7 @@ StaticModel::operator=(const StaticModel& m)
 }
 
 StaticModel::StaticModel(const DynamicModel& m) :
-    ModelTree {m.symbol_table, m.num_constants, m.external_functions_table}
+    ModelTree {m.symbol_table, m.num_constants, m.external_functions_table, m.heterogeneity_table}
 {
   // Convert model local variables (need to be done first)
   for (int it : m.local_variables_vector)
@@ -87,17 +90,18 @@ StaticModel::StaticModel(const DynamicModel& m) :
         if (dynamic_equations.contains(i))
           {
             auto [static_only_equations, static_only_equations_lineno,
-                  static_only_equations_equation_tags]
+                  static_only_complementarity_conditions, static_only_equations_equation_tags]
                 = m.getStaticOnlyEquationsInfo();
 
             addEquation(static_only_equations[static_only_index]->toStatic(*this),
                         static_only_equations_lineno[static_only_index],
+                        static_only_complementarity_conditions[static_only_index],
                         static_only_equations_equation_tags.getTagsByEqn(static_only_index));
             static_only_index++;
           }
         else
           addEquation(m.equations[i]->toStatic(*this), m.equations_lineno[i],
-                      m.equation_tags.getTagsByEqn(i));
+                      m.complementarity_conditions[i], m.equation_tags.getTagsByEqn(i));
       }
     catch (DataTree::DivisionByZeroException)
       {
@@ -132,10 +136,8 @@ StaticModel::writeStaticBytecode(const string& basename) const
   int u_count_int {writeBytecodeBinFile(basename + "/model/bytecode/static.bin", false)};
 
   Bytecode::Writer code_file {basename + "/model/bytecode/static.cod"};
-  vector<int> eq_idx(equations.size());
-  iota(eq_idx.begin(), eq_idx.end(), 0);
-  vector<int> endo_idx(symbol_table.endo_nbr());
-  iota(endo_idx.begin(), endo_idx.end(), 0);
+  auto eq_idx = views::iota(0, static_cast<int>(equations.size()));
+  auto endo_idx = views::iota(0, symbol_table.endo_nbr());
 
   // Declare temporary terms and the (single) block
   code_file << Bytecode::FDIMST {static_cast<int>(temporary_terms_derivatives[0].size()
@@ -144,8 +146,8 @@ StaticModel::writeStaticBytecode(const string& basename) const
                                       BlockSimulationType::solveForwardComplete,
                                       0,
                                       symbol_table.endo_nbr(),
-                                      endo_idx,
-                                      eq_idx,
+                                      {endo_idx.begin(), endo_idx.end()},
+                                      {eq_idx.begin(), eq_idx.end()},
                                       false,
                                       u_count_int,
                                       symbol_table.endo_nbr()};
@@ -250,6 +252,8 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder,
            << endl;
       exit(EXIT_FAILURE);
     }
+
+  computeMCPEquationsReordering();
 }
 
 void
@@ -283,10 +287,6 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string&
     }
   create_directories(model_dir / "bytecode" / "block");
 
-  // Legacy representation
-  if (use_dll)
-    writeModelCFile<false>(basename, mexext, matlabroot);
-
   /* PlannerObjective subclass or discretionary optimal policy models don’t
      have as many variables as equations; bytecode does not support that
      case */
@@ -305,10 +305,13 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string&
 
   writeSetAuxiliaryVariablesFile<false>(basename, julia);
 
-  // Support for model debugging
   if (!julia)
+    {
+      writeComplementarityConditionsFile<false>(basename);
+      // Support for model debugging
       writeDebugModelMFiles<false>(basename);
     }
+}
 
 bool
 StaticModel::exoPresentInEqs() const
@@ -330,7 +333,12 @@ StaticModel::writeDriverOutput(ostream& output) const
   if (block_decomposed)
     writeBlockDriverOutput(output);
 
-  writeDriverSparseIndicesHelper<false, false>(output);
+  writeDriverSparseIndicesHelper("static", output);
+
+  output << "M_.static_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
 
 void
@@ -634,7 +642,7 @@ StaticModel::computeRamseyMultipliersDerivatives(int ramsey_orig_endo_nbr, bool
          Lagrange multiplier. We use the guarantee given by SymbolTable that
          symbol IDs are increasing. */
       if (varexpr->symb_id > *mult_symb_ids.crbegin())
-        recursive_variables.emplace(getDerivID(varexpr->symb_id, 0), aux_eq);
+        recursive_variables.emplace(varexpr->getDerivID(), aux_eq);
     }
 
   // Compute the chain rule derivatives w.r.t. multipliers
@@ -658,7 +666,7 @@ StaticModel::computeRamseyMultipliersDerivatives(int ramsey_orig_endo_nbr, bool
   if (no_tmp_terms)
     for (auto& it : temp_terms_map)
       erase_if(it.second, [](expr_t e) { return !dynamic_cast<AbstractExternalFunctionNode*>(e); });
-  copy(temp_terms_map[{1, 0}].begin(), temp_terms_map[{1, 0}].end(),
+  ranges::copy(temp_terms_map[{1, 0}],
                inserter(ramsey_multipliers_derivatives_temporary_terms,
                         ramsey_multipliers_derivatives_temporary_terms.begin()));
   for (int idx {0}; auto it : ramsey_multipliers_derivatives_temporary_terms)
@@ -704,12 +712,7 @@ StaticModel::writeRamseyMultipliersDerivativesMFile(const string& basename,
 
   writeRamseyMultipliersDerivativesHelper<output_type>(output_file);
 
-  // On MATLAB < R2020a, sparse() does not accept int32 indices
-  output_file << "if ~isoctave && matlab_ver_less_than('9.8')" << endl
-              << "    sparse_rowval = double(sparse_rowval);" << endl
-              << "    sparse_colval = double(sparse_colval);" << endl
-              << "end" << endl
-              << "g1m = sparse(sparse_rowval, sparse_colval, g1m_v, " << ramsey_orig_endo_nbr
+  output_file << "g1m = sparse(sparse_rowval, sparse_colval, g1m_v, " << ramsey_orig_endo_nbr
               << ", " << symbol_table.getLagrangeMultipliers().size() << ");" << endl
               << "end" << endl;
   output_file.close();

src/StaticModel.hh

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

src/SubModel.cc

 /*
- * Copyright © 2018-2023 Dynare Team
+ * Copyright © 2018-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -58,8 +58,7 @@ TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg,
     {
       vector<int> nontrend_vec;
       for (auto eq : it.second)
-        if (find(target_eqnums[it.first].begin(), target_eqnums[it.first].end(), eq)
-            == target_eqnums[it.first].end())
+        if (ranges::find(target_eqnums[it.first], eq) == target_eqnums[it.first].end())
           nontrend_vec.push_back(eq);
       nontarget_eqnums[it.first] = nontrend_vec;
     }
@@ -71,12 +70,12 @@ TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg,
       vector<int> eqnumsv = getEqNums(name);
       for (int nontrend_it : getNonTargetEqNums(name))
         nontarget_lhs_vec.push_back(
-            lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), nontrend_it))));
+            lhsv.at(ranges::distance(eqnumsv.begin(), ranges::find(eqnumsv, nontrend_it))));
       nontarget_lhs[name] = nontarget_lhs_vec;
 
       for (int trend_it : getTargetEqNums(name))
         target_lhs_vec.push_back(
-            lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), trend_it))));
+            lhsv.at(ranges::distance(eqnumsv.begin(), ranges::find(eqnumsv, trend_it))));
       target_lhs[name] = target_lhs_vec;
     }
 }
@@ -274,8 +273,7 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
         output << it + 1 << " ";
       output << "];" << endl << "M_.trend_component." << name << ".targets = [";
       for (auto it : eqnums.at(name))
-        if (find(target_eqnums.at(name).begin(), target_eqnums.at(name).end(), it)
-            == target_eqnums.at(name).end())
+        if (ranges::find(target_eqnums.at(name), it) == target_eqnums.at(name).end())
           output << "false ";
         else
           output << "true ";
@@ -322,8 +320,7 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
       vector<string> eqtags_vec = eqtags.at(name);
       output << "M_.trend_component." << name << ".target_eqn = [";
       for (const auto& it : target_eqtags_vec)
-        output << distance(eqtags_vec.begin(), find(eqtags_vec.begin(), eqtags_vec.end(), it)) + 1
-               << " ";
+        output << ranges::distance(eqtags_vec.begin(), ranges::find(eqtags_vec, it)) + 1 << " ";
       output << "];" << endl;
 
       vector<int> target_lhs_vec = getTargetLhs(name);
@@ -336,9 +333,8 @@ TrendComponentModelTable::writeOutput(const string& basename, ostream& output) c
       for (const auto& [key, expr] : AR.at(name))
         {
           auto [eqn, lag, lhs_symb_id] = key;
-          int colidx = static_cast<int>(
-              distance(nontarget_lhs_vec.begin(),
-                       find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), lhs_symb_id)));
+          int colidx = static_cast<int>(ranges::distance(
+              nontarget_lhs_vec.begin(), ranges::find(nontarget_lhs_vec, lhs_symb_id)));
           ar_ec_output << "    AR(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
           expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
           ar_ec_output << ";" << endl;
@@ -488,7 +484,7 @@ VarModelTable::writeOutput(const string& basename, ostream& output) const
         {
           auto [eqn, lag, lhs_symb_id] = key;
           int colidx
-              = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
+              = static_cast<int>(ranges::distance(lhs.begin(), ranges::find(lhs, lhs_symb_id)));
           ar_output << "    ar(" << eqn + 1 << "," << colidx + 1 << "," << lag << ") = ";
           expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
           ar_output << ";" << endl;
@@ -498,7 +494,7 @@ VarModelTable::writeOutput(const string& basename, ostream& output) const
         {
           auto [eqn, lhs_symb_id] = key;
           int colidx
-              = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
+              = static_cast<int>(ranges::distance(lhs.begin(), ranges::find(lhs, lhs_symb_id)));
           if (eqn != colidx)
             {
               ar_output << "        a0(" << eqn + 1 << "," << colidx + 1 << ") = ";
@@ -803,9 +799,8 @@ VarExpectationModelTable::writeOutput(ostream& output) const
              << mstruct << ".expr.constants = [ " << constants_list.str() << " ];" << endl;
 
       if (auto disc_var = dynamic_cast<const VariableNode*>(discount.at(name)); disc_var)
-        output << mstruct
-               << ".discount_index = " << symbol_table.getTypeSpecificID(disc_var->symb_id) + 1
-               << ';' << endl;
+        output << mstruct << ".discount_index = " << disc_var->getTypeSpecificID() + 1 << ';'
+               << endl;
       else
         {
           output << mstruct << ".discount_value = ";
@@ -1285,8 +1280,8 @@ PacModelTable::transformPass(const lag_equivalence_table_t& unary_ops_nodes,
 
           // Associate the coefficients of the linear combination with the right components
           for (auto [var, coeff] : terms)
-            if (auto it = find_if(
-                    components.begin(), components.end(),
+            if (auto it = ranges::find_if(
+                    components,
                     [&](const auto& v) { return get<0>(v) == dynamic_model.AddVariable(var); });
                 it != components.end())
               get<4>(*it) = coeff;
@@ -1384,11 +1379,11 @@ PacModelTable::transformPass(const lag_equivalence_table_t& unary_ops_nodes,
         {
           if (target_info.contains(name))
             {
-              cerr << "ERROR: the block 'pac_target_info(" << name
-                   << ")' is not supported in the context of a PAC model with model-consistent "
-                      "expectations (MCE)."
-                   << endl;
-              exit(EXIT_FAILURE);
+              assert(growth_correction_term == dynamic_model.Zero);
+              dynamic_model.computePacModelConsistentExpectationSubstitutionWithComponents(
+                  name, symbol_table.getID(discount[name]), pacEquationMaxLag(name),
+                  diff_subst_table, aux_param_symb_ids, get<2>(target_info[name]),
+                  pac_expectation_substitution);
             }
           else
             dynamic_model.computePacModelConsistentExpectationSubstitution(
@@ -1507,7 +1502,7 @@ PacModelTable::writeOutput(ostream& output) const
 
   // Write the auxiliary variable IDs created for the pac_expectation operator
   for (auto& [name, id] : aux_var_symb_ids)
-    output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z1" : "aux_id")
+    output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z" : "aux_id")
            << " = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
 
   // Write PAC equation name info
@@ -1687,9 +1682,8 @@ PacModelTable::writeOutput(ostream& output) const
         string fieldname = "M_.pac." + name + ".components(" + to_string(component_idx) + ")";
         output << fieldname << ".aux_id = " << symbol_table.getTypeSpecificID(auxname) + 1 << ";"
                << endl
-               << fieldname << ".endo_var = "
-               << symbol_table.getTypeSpecificID(dynamic_cast<VariableNode*>(component)->symb_id)
-                      + 1
+               << fieldname
+               << ".endo_var = " << dynamic_cast<VariableNode*>(component)->getTypeSpecificID() + 1
                << ";" << endl
                << fieldname << ".kind = '" << kindToString(kind) << "';" << endl
                << fieldname << ".h_param_indices = [";
@@ -1709,6 +1703,21 @@ PacModelTable::writeOutput(ostream& output) const
           }
         component_idx++;
       }
+
+  for (auto& [name, val] : target_info)
+    if (aux_model_name.at(name).empty())
+      {
+        string pac_model_name = "M_.pac." + name + ".mce.";
+        output << pac_model_name << "z = NaN(" << get<2>(val).size() << ",1);" << endl;
+        for (int component_idx {1};
+             auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                    h_indices, original_growth_component, growth_component_info] : get<2>(val))
+          {
+            output << pac_model_name << "z(" << component_idx
+                   << ") = " << symbol_table.getTypeSpecificID(auxname) + 1 << ";" << endl;
+            component_idx++;
+          }
+      }
 }
 
 void

src/SymbolList.cc

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

src/SymbolList.hh

@@ -63,7 +63,7 @@ public:
   empty() const
   {
     return symbols.empty();
-  };
+  }
   //! Return the list of symbols
   [[nodiscard]] vector<string> getSymbols() const;
 };

src/SymbolTable.cc

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -31,7 +31,8 @@
 
 int
 SymbolTable::addSymbol(const string& name, SymbolType type, const string& tex_name,
-                       const vector<pair<string, string>>& partition_value) noexcept(false)
+                       const vector<pair<string, string>>& partition_value,
+                       const optional<int>& heterogeneity_dimension) noexcept(false)
 {
   if (frozen)
     throw FrozenException();
@@ -78,13 +79,20 @@ SymbolTable::addSymbol(const string& name, SymbolType type, const string& tex_na
         pmv[it.first] = it.second;
       partition_value_map[id] = pmv;
     }
+
+  assert(!isHeterogeneous(type)
+         || (heterogeneity_dimension.has_value() && *heterogeneity_dimension >= 0
+             && *heterogeneity_dimension < heterogeneity_table.size()));
+  if (isHeterogeneous(type))
+    heterogeneity_dimensions.emplace(id, *heterogeneity_dimension);
+
   return id;
 }
 
 int
 SymbolTable::addSymbol(const string& name, SymbolType type) noexcept(false)
 {
-  return addSymbol(name, type, "", {});
+  return addSymbol(name, type, "", {}, {});
 }
 
 void
@@ -95,6 +103,10 @@ SymbolTable::freeze() noexcept(false)
 
   frozen = true;
 
+  het_endo_ids.resize(heterogeneity_table.size());
+  het_exo_ids.resize(heterogeneity_table.size());
+  het_param_ids.resize(heterogeneity_table.size());
+
   for (int i = 0; i < static_cast<int>(symbol_table.size()); i++)
     {
       int tsi;
@@ -116,6 +128,18 @@ SymbolTable::freeze() noexcept(false)
           tsi = param_ids.size();
           param_ids.push_back(i);
           break;
+        case SymbolType::heterogeneousEndogenous:
+          tsi = het_endo_ids.at(heterogeneity_dimensions.at(i)).size();
+          het_endo_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+          break;
+        case SymbolType::heterogeneousExogenous:
+          tsi = het_exo_ids.at(heterogeneity_dimensions.at(i)).size();
+          het_exo_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+          break;
+        case SymbolType::heterogeneousParameter:
+          tsi = het_param_ids.at(heterogeneity_dimensions.at(i)).size();
+          het_param_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+          break;
         default:
           continue;
         }
@@ -131,12 +155,18 @@ SymbolTable::unfreeze()
   exo_ids.clear();
   exo_det_ids.clear();
   param_ids.clear();
+  het_endo_ids.clear();
+  het_exo_ids.clear();
+  het_param_ids.clear();
   type_specific_ids.clear();
 }
 
 void
 SymbolTable::changeType(int id, SymbolType newtype) noexcept(false)
 {
+  // FIXME: implement switch to heterogeneous variable; dimension will have to be provided
+  assert(!isHeterogeneous(newtype));
+
   if (frozen)
     throw FrozenException();
 
@@ -146,7 +176,8 @@ SymbolTable::changeType(int id, SymbolType newtype) noexcept(false)
 }
 
 int
-SymbolTable::getID(SymbolType type, int tsid) const noexcept(false)
+SymbolTable::getID(SymbolType type, int tsid, const optional<int>& heterogeneity_dimension) const
+    noexcept(false)
 {
   if (!frozen)
     throw NotYetFrozenException();
@@ -155,26 +186,44 @@ SymbolTable::getID(SymbolType type, int tsid) const noexcept(false)
     {
     case SymbolType::endogenous:
       if (tsid < 0 || tsid >= static_cast<int>(endo_ids.size()))
-        throw UnknownTypeSpecificIDException {tsid, type};
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return endo_ids[tsid];
     case SymbolType::exogenous:
       if (tsid < 0 || tsid >= static_cast<int>(exo_ids.size()))
-        throw UnknownTypeSpecificIDException {tsid, type};
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return exo_ids[tsid];
     case SymbolType::exogenousDet:
       if (tsid < 0 || tsid >= static_cast<int>(exo_det_ids.size()))
-        throw UnknownTypeSpecificIDException {tsid, type};
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return exo_det_ids[tsid];
     case SymbolType::parameter:
       if (tsid < 0 || tsid >= static_cast<int>(param_ids.size()))
-        throw UnknownTypeSpecificIDException {tsid, type};
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return param_ids[tsid];
+    case SymbolType::heterogeneousEndogenous:
+      assert(heterogeneity_dimension.has_value());
+      if (tsid < 0 || tsid >= static_cast<int>(het_endo_ids.at(*heterogeneity_dimension).size()))
+        throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+      else
+        return het_endo_ids.at(*heterogeneity_dimension).at(tsid);
+    case SymbolType::heterogeneousExogenous:
+      assert(heterogeneity_dimension.has_value());
+      if (tsid < 0 || tsid >= static_cast<int>(het_exo_ids.at(*heterogeneity_dimension).size()))
+        throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+      else
+        return het_exo_ids.at(*heterogeneity_dimension).at(tsid);
+    case SymbolType::heterogeneousParameter:
+      assert(heterogeneity_dimension.has_value());
+      if (tsid < 0 || tsid >= static_cast<int>(het_param_ids.at(*heterogeneity_dimension).size()))
+        throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+      else
+        return het_param_ids.at(*heterogeneity_dimension).at(tsid);
     default:
-      throw UnknownTypeSpecificIDException {tsid, type};
+      throw UnknownTypeSpecificIDException {tsid, type, {}};
     }
 }
 
@@ -348,6 +397,8 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
           case AuxVarType::expectation:
           case AuxVarType::pacExpectation:
           case AuxVarType::pacTargetNonstationary:
+          case AuxVarType::aggregationOp:
+          case AuxVarType::heterogeneousMultiplier:
             break;
           case AuxVarType::endoLag:
           case AuxVarType::exoLag:
@@ -419,6 +470,54 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
         output << getTypeSpecificID(varexob) + 1 << " ";
       output << " ];" << endl;
     }
+
+  // Heterogeneous symbols
+  // FIXME: the following helper could be used to simplify non-heterogenous variables
+  auto print_symb_names = [this, &output](const string& field, const auto& symb_ids) {
+    auto helper = [this, &output, &symb_ids](auto nameMethod) {
+      for (bool first_printed {false}; int symb_id : symb_ids)
+        {
+          if (exchange(first_printed, true))
+            output << "; ";
+          output << "'" << (this->*nameMethod)(symb_id) << "'";
+        }
+    };
+    output << field << " = {";
+    helper(&SymbolTable::getName);
+    output << "};" << endl << field << "_tex = {";
+    helper(&SymbolTable::getTeXName);
+    output << "};" << endl << field << "_long = {";
+    helper(&SymbolTable::getLongName);
+    output << "};" << endl;
+  };
+  for (int het_dim {0}; het_dim < heterogeneity_table.size(); het_dim++)
+    {
+      const string basefield {"M_.heterogeneity(" + to_string(het_dim + 1) + ")."};
+
+      output << basefield << "endo_nbr = " << het_endo_nbr(het_dim) << ";" << endl;
+      print_symb_names(basefield + "endo_names", het_endo_ids.at(het_dim));
+
+      output << basefield << "orig_endo_nbr = " << het_orig_endo_nbr(het_dim) << ";" << endl;
+
+      output << basefield << "exo_nbr = " << het_exo_nbr(het_dim) << ";" << endl;
+      print_symb_names(basefield + "exo_names", het_exo_ids.at(het_dim));
+
+      output << basefield << "param_nbr = " << het_param_nbr(het_dim) << ";" << endl;
+      print_symb_names(basefield + "param_names", het_param_ids.at(het_dim));
+
+      const vector<AuxVarInfo>& hav = het_aux_vars[het_dim];
+      if (hav.size() == 0)
+        output << basefield << "aux_vars = [];";
+      else
+        for (int i = 0; i < static_cast<int>(hav.size()); i++)
+          {
+            output << basefield << "aux_vars(" << i + 1
+                   << ").endo_index = " << getTypeSpecificID(hav[i].symb_id) + 1 << ";" << endl
+                   << basefield << "aux_vars(" << i + 1 << ").type = " << hav[i].get_type_id()
+                   << ";" << basefield << "aux_vars(" << i + 1
+                   << ").eq_nbr = " << hav[i].equation_number_for_multiplier + 1 << ";" << endl;
+          }
+    }
 }
 
 int
@@ -630,6 +729,27 @@ SymbolTable::addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op,
   return symb_id;
 }
 
+int
+SymbolTable::addHeterogeneousMultiplierAuxiliaryVar(int het_dim, int index,
+                                                    const string& varname) noexcept(false)
+{
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::heterogeneousEndogenous, "", {}, het_dim);
+    }
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  het_aux_vars[het_dim].emplace_back(symb_id, AuxVarType::heterogeneousMultiplier, 0, 0, index, 0,
+                                     nullptr, "");
+  return symb_id;
+}
+
 int
 SymbolTable::addMultiplierAuxiliaryVar(int index) noexcept(false)
 {
@@ -714,6 +834,28 @@ SymbolTable::addPacTargetNonstationaryAuxiliaryVar(const string& name, expr_t ex
   return symb_id;
 }
 
+int
+SymbolTable::addAggregationOpAuxiliaryVar(const string& name, expr_t expr_arg)
+{
+  int symb_id {[&] {
+    try
+      {
+        return addSymbol(name, SymbolType::endogenous);
+      }
+    catch (AlreadyDeclaredException&)
+      {
+        cerr << "ERROR: the variable/parameter '" << name
+             << "' conflicts with a variable that will be generated for an aggregation operator. "
+                "Please rename it."
+             << endl;
+        exit(EXIT_FAILURE);
+      }
+  }()};
+
+  aux_vars.emplace_back(symb_id, AuxVarType::aggregationOp, 0, 0, 0, 0, expr_arg, "");
+  return symb_id;
+}
+
 int
 SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const noexcept(false)
 {
@@ -812,13 +954,13 @@ SymbolTable::observedVariablesNbr() const
 bool
 SymbolTable::isObservedVariable(int symb_id) const
 {
-  return find(varobs.begin(), varobs.end(), symb_id) != varobs.end();
+  return ranges::find(varobs, symb_id) != varobs.end();
 }
 
 int
 SymbolTable::getObservedVariableIndex(int symb_id) const
 {
-  auto it = find(varobs.begin(), varobs.end(), symb_id);
+  auto it = ranges::find(varobs, symb_id);
   assert(it != varobs.end());
   return static_cast<int>(it - varobs.begin());
 }
@@ -827,7 +969,7 @@ void
 SymbolTable::addObservedExogenousVariable(int symb_id) noexcept(false)
 {
   validateSymbID(symb_id);
-  assert(getType(symb_id) != SymbolType::endogenous);
+  assert(getType(symb_id) == SymbolType::exogenous);
   varexobs.push_back(symb_id);
 }
 
@@ -840,13 +982,13 @@ SymbolTable::observedExogenousVariablesNbr() const
 bool
 SymbolTable::isObservedExogenousVariable(int symb_id) const
 {
-  return find(varexobs.begin(), varexobs.end(), symb_id) != varexobs.end();
+  return ranges::find(varexobs, symb_id) != varexobs.end();
 }
 
 int
 SymbolTable::getObservedExogenousVariableIndex(int symb_id) const
 {
-  auto it = find(varexobs.begin(), varexobs.end(), symb_id);
+  auto it = ranges::find(varexobs, symb_id);
   assert(it != varexobs.end());
   return static_cast<int>(it - varexobs.begin());
 }
@@ -895,14 +1037,13 @@ SymbolTable::getEndogenous() const
 bool
 SymbolTable::isAuxiliaryVariable(int symb_id) const
 {
-  return any_of(aux_vars.begin(), aux_vars.end(),
-                [=](const auto& av) { return av.symb_id == symb_id; });
+  return ranges::any_of(aux_vars, [=](const auto& av) { return av.symb_id == symb_id; });
 }
 
 bool
 SymbolTable::isDiffAuxiliaryVariable(int symb_id) const
 {
-  return any_of(aux_vars.begin(), aux_vars.end(), [=](const auto& av) {
+  return ranges::any_of(aux_vars, [=](const auto& av) {
     return av.symb_id == symb_id
            && (av.type == AuxVarType::diff || av.type == AuxVarType::diffLag
                || av.type == AuxVarType::diffLead);
@@ -995,6 +1136,8 @@ SymbolTable::writeJsonOutput(ostream& output) const
             case AuxVarType::expectation:
             case AuxVarType::pacExpectation:
             case AuxVarType::pacTargetNonstationary:
+            case AuxVarType::aggregationOp:
+            case AuxVarType::heterogeneousMultiplier:
               break;
             case AuxVarType::endoLag:
             case AuxVarType::exoLag:
@@ -1029,6 +1172,25 @@ SymbolTable::writeJsonOutput(ostream& output) const
         }
       output << "]" << endl;
     }
+
+  if (!heterogeneity_table.empty())
+    {
+      output << R"(, "heterogeneous_symbols": [)";
+      for (int i {0}; i < heterogeneity_table.size(); i++)
+        {
+          if (i != 0)
+            output << ", ";
+          output << R"({ "dimension": ")" << heterogeneity_table.getName(i)
+                 << R"(", "endogenous": )";
+          writeJsonVarVector(output, het_endo_ids.at(i));
+          output << R"(, "exogenous": )";
+          writeJsonVarVector(output, het_exo_ids.at(i));
+          output << R"(, "parameters": )";
+          writeJsonVarVector(output, het_param_ids.at(i));
+          output << "}";
+        }
+      output << "]" << endl;
+    }
 }
 
 void
@@ -1088,3 +1250,14 @@ SymbolTable::getLagrangeMultipliers() const
       r.insert(aux_var.symb_id);
   return r;
 }
+
+int
+SymbolTable::getHeterogeneityDimension(int symb_id) const
+{
+  validateSymbID(symb_id);
+  auto it = heterogeneity_dimensions.find(symb_id);
+  if (it != heterogeneity_dimensions.end())
+    return it->second;
+  else
+    throw NonHeteregeneousSymbolException {symb_id};
+}

src/SymbolTable.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -30,6 +30,7 @@
 
 #include "CommonEnums.hh"
 #include "ExprNode.hh"
+#include "HeterogeneityTable.hh"
 
 using namespace std;
 
@@ -53,7 +54,11 @@ enum class AuxVarType
   diffLead = 11, //!< Variable for timing between Diff operators (lead)
   pacExpectation = 12, //!< Variable created for the substitution of the pac_expectation operator
   pacTargetNonstationary
-  = 13 //!< Variable created for the substitution of the pac_target_nonstationary operator
+  = 13, //!< Variable created for the substitution of the pac_target_nonstationary operator
+  aggregationOp
+  = 14, // Substitute for an aggregation operator in a heterogeneous setup, such as SUM()
+  heterogeneousMultiplier = 15 /* Multiplier for bound conditions of complementarity conditions in
+                                  the heterogeneous equations */
 };
 
 //! Information on some auxiliary variables
@@ -106,6 +111,8 @@ struct AuxVarInfo
 class SymbolTable
 {
 private:
+  HeterogeneityTable& heterogeneity_table;
+
   //! Has method freeze() been called?
   bool frozen {false};
 
@@ -123,6 +130,8 @@ private:
   map<int, map<string, string>> partition_value_map;
   //! Maps IDs to types
   vector<SymbolType> type_table;
+  // Maps IDs of heterogenous symbols to heterogeneity dimension IDs
+  map<int, int> heterogeneity_dimensions;
 
   //! Maps symbol IDs to type specific IDs
   map<int, int> type_specific_ids;
@@ -135,9 +144,22 @@ private:
   vector<int> exo_det_ids;
   //! Maps type specific IDs of parameters to symbol IDs
   vector<int> param_ids;
+  /* Maps type specific IDs of heterogeneous endogenous to symbol IDs (outer vector is for
+     heterogeneity dimensions) */
+  vector<vector<int>> het_endo_ids;
+  /* Maps type specific IDs of heterogeneous exogenous to symbol IDs (outer vector is for
+     heterogeneity dimensions) */
+  vector<vector<int>> het_exo_ids;
+  /* Maps type specific IDs of heterogeneous parameters to symbol IDs (outer vector is for
+     heterogeneity dimensions) */
+  vector<vector<int>> het_param_ids;
+
   //! Information about auxiliary variables
   vector<AuxVarInfo> aux_vars;
 
+  //! Information about heterogeneous auxiliary variables
+  vector<vector<AuxVarInfo>> het_aux_vars;
+
   //! Stores the predetermined variables (by symbol IDs)
   set<int> predetermined_variables;
 
@@ -168,6 +190,7 @@ public:
   {
     const int tsid;
     const SymbolType type;
+    const optional<int> heterogeneity_dimension;
   };
   /* Thrown when requesting the type specific ID of a symbol which doesn’t
      have one */
@@ -202,6 +225,12 @@ public:
     }
   };
 
+  // Thrown by getHeterogeneityDimension() on non-heterogeneous symbols
+  struct NonHeteregeneousSymbolException
+  {
+    const int id;
+  };
+
 private:
   //! Factorized code for adding aux lag variables
   int addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag,
@@ -214,11 +243,19 @@ private:
   inline void validateSymbID(int symb_id) const noexcept(false);
 
 public:
+  SymbolTable(HeterogeneityTable& heterogeneity_table_arg) :
+      heterogeneity_table {heterogeneity_table_arg}
+  {
+  }
+
   //! Add a symbol
-  /*! Returns the symbol ID */
+  /* Returns the symbol ID.
+     heterogeneity_dimension must be defined if this is a heterogeneous symbol (otherwise it is
+     ignored) */
   int addSymbol(const string& name, SymbolType type, const string& tex_name,
-                const vector<pair<string, string>>& partition_value) noexcept(false);
-  //! Add a symbol without its TeX name (will be equal to its name)
+                const vector<pair<string, string>>& partition_value,
+                const optional<int>& heterogeneity_dimension) noexcept(false);
+  //! Add a (non-heterogenous) symbol without its TeX name (will be equal to its name)
   /*! Returns the symbol ID */
   int addSymbol(const string& name, SymbolType type) noexcept(false);
   //! Adds an auxiliary variable for endogenous with lead >= 2
@@ -250,6 +287,14 @@ public:
     \return the symbol ID of the new symbol
   */
   int addExpectationAuxiliaryVar(int information_set, int index, expr_t arg) noexcept(false);
+  /* Adds an auxiliary variable for the bound-specific multiplier of the MCPs and returns its
+     symbolID.
+     – het_dim is the heterogeneity dimension of the model
+     – index is the equation's index
+     – varname is the multiplier name
+  */
+  int addHeterogeneousMultiplierAuxiliaryVar(int het_dim, int index,
+                                             const string& varname) noexcept(false);
   //! Adds an auxiliary variable for the multiplier for the FOCs of the Ramsey Problem
   /*!
     \param[in] index Used to construct the variable name
@@ -313,12 +358,16 @@ public:
   int addPacExpectationAuxiliaryVar(const string& name, expr_t expr_arg);
   //! An auxiliary variable for a pac_target_nonstationary operator
   int addPacTargetNonstationaryAuxiliaryVar(const string& name, expr_t expr_arg);
+  // An auxiliary variable for an aggregation operator (e.g. SUM(yh) where yh is heterogeneous)
+  int addAggregationOpAuxiliaryVar(const string& name, expr_t expr_arg);
+  //! Set the size of het_aux_vars
+  inline void resizeHetAuxVars();
   //! Returns the number of auxiliary variables
   [[nodiscard]] int
   AuxVarsSize() const
   {
     return aux_vars.size();
-  };
+  }
   //! Tests if symbol already exists
   [[nodiscard]] inline bool exists(const string& name) const;
   //! Get symbol name (by ID)
@@ -340,7 +389,9 @@ public:
   //! Get ID (by name)
   [[nodiscard]] inline int getID(const string& name) const noexcept(false);
   //! Get ID (by type specific ID)
-  [[nodiscard]] int getID(SymbolType type, int tsid) const noexcept(false);
+  [[nodiscard]] int getID(SymbolType type, int tsid,
+                          const optional<int>& heterogeneity_dimension = nullopt) const
+      noexcept(false);
   //! Freeze symbol table
   void freeze() noexcept(false);
   //! unreeze symbol table
@@ -360,6 +411,15 @@ public:
   [[nodiscard]] inline int exo_det_nbr() const noexcept(false);
   //! Get number of parameters
   [[nodiscard]] inline int param_nbr() const noexcept(false);
+  //! Get number of heterogeneous endogenous variables along a given dimension
+  [[nodiscard]] inline int het_endo_nbr(int het_dim) const noexcept(false);
+  //! Get number of user-declared heterogeneous endogenous variables (without
+  //! the auxiliary variables)
+  [[nodiscard]] inline int het_orig_endo_nbr(int het_dim) const noexcept(false);
+  //! Get number of heterogeneous exogenous variables along a given dimension
+  [[nodiscard]] inline int het_exo_nbr(int het_dim) const noexcept(false);
+  //! Get number of heterogeneous parameters along a given dimension
+  [[nodiscard]] inline int het_param_nbr(int het_dim) const noexcept(false);
   //! Returns the greatest symbol ID (the smallest is zero)
   [[nodiscard]] inline int maxID() const;
   //! Get number of user-declared endogenous variables (without the auxiliary variables)
@@ -423,6 +483,9 @@ public:
   [[nodiscard]] const set<int>& getVariablesWithLogTransform() const;
   // Returns all Lagrange multipliers
   [[nodiscard]] set<int> getLagrangeMultipliers() const;
+  /* Get heterogeneity dimension of a given symbol. Throws NonHeterogeneousSymbolException
+     if there is no such dimension. */
+  [[nodiscard]] int getHeterogeneityDimension(int symb_id) const;
 };
 
 inline void
@@ -432,6 +495,12 @@ SymbolTable::validateSymbID(int symb_id) const noexcept(false)
     throw UnknownSymbolIDException {symb_id};
 }
 
+inline void
+SymbolTable::resizeHetAuxVars()
+{
+  het_aux_vars.resize(heterogeneity_table.size());
+}
+
 inline bool
 SymbolTable::exists(const string& name) const
 {
@@ -537,6 +606,42 @@ SymbolTable::param_nbr() const noexcept(false)
   return param_ids.size();
 }
 
+inline int
+SymbolTable::het_endo_nbr(int het_dim) const noexcept(false)
+{
+  if (!frozen)
+    throw NotYetFrozenException();
+
+  return het_endo_ids.at(het_dim).size();
+}
+
+inline int
+SymbolTable::het_orig_endo_nbr(int het_dim) const noexcept(false)
+{
+  if (!frozen)
+    throw NotYetFrozenException();
+
+  return het_endo_ids.at(het_dim).size() - het_aux_vars.at(het_dim).size();
+}
+
+inline int
+SymbolTable::het_exo_nbr(int het_dim) const noexcept(false)
+{
+  if (!frozen)
+    throw NotYetFrozenException();
+
+  return het_exo_ids.at(het_dim).size();
+}
+
+inline int
+SymbolTable::het_param_nbr(int het_dim) const noexcept(false)
+{
+  if (!frozen)
+    throw NotYetFrozenException();
+
+  return het_param_ids.at(het_dim).size();
+}
+
 inline int
 SymbolTable::maxID() const
 {

src/WarningConsolidation.cc

 /*
- * Copyright © 2012-2022 Dynare Team
+ * Copyright © 2012-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -18,39 +18,33 @@
  */
 
 #include "WarningConsolidation.hh"
-#include <ostream>
 
-WarningConsolidation&
-operator<<(WarningConsolidation& wcc, const string& warning)
+ostream&
+operator<<(ostream& stream, const Dynare::location& l)
 {
-  if (wcc.no_warn)
-    return wcc;
+  stream << *l.begin.filename << ": line " << l.begin.line;
+  if (l.begin.line == l.end.line)
+    if (l.begin.column == l.end.column - 1)
+      stream << ", col " << l.begin.column;
+    else
+      stream << ", cols " << l.begin.column << "-" << l.end.column - 1;
+  else
+    stream << ", col " << l.begin.column << " -"
+           << " line " << l.end.line << ", col " << l.end.column - 1;
 
-  cerr << warning;
-  wcc.addWarning(warning);
-  return wcc;
-};
+  return stream;
+}
 
-WarningConsolidation&
-operator<<(WarningConsolidation& wcc, const Dynare::location& loc)
+void
+WarningConsolidation::incrementWarnings(const string& msg)
 {
-  if (wcc.no_warn)
-    return wcc;
-
-  stringstream ostr;
-  Dynare::position last = loc.end - 1;
-  ostr << loc.begin;
-  if (last.filename && (!loc.begin.filename || *loc.begin.filename != *last.filename))
-    ostr << '-' << last;
-  else if (loc.begin.line != last.line)
-    ostr << '-' << last.line << '.' << last.column;
-  else if (loc.begin.column != last.column)
-    ostr << '-' << last.column;
-
-  cerr << ostr.str();
-  wcc.addWarning(ostr.str());
-  return wcc;
-};
+  size_t p {0};
+  while ((p = msg.find('\n', p)) != string::npos)
+    {
+      p++;
+      num_warnings++;
+    }
+}
 
 WarningConsolidation&
 operator<<(WarningConsolidation& wcc, ostream& (*pf)(ostream&))
@@ -58,49 +52,12 @@ operator<<(WarningConsolidation& wcc, ostream& (*pf)(ostream&))
   if (wcc.no_warn)
     return wcc;
 
-  cerr << pf;
-  wcc.addWarning(pf);
-  return wcc;
-}
+  ostringstream ostr;
+  ostr << pf;
 
-void
-WarningConsolidation::writeOutput(ostream& output) const
-{
-  if (warnings.str().empty())
-    return;
-
-  output << "disp([char(10) 'Dynare Preprocessor Warning(s) Encountered:']);" << endl;
+  cerr << ostr.str();
 
-  bool writedisp = true;
-  string warningsstr = warnings.str();
-  for (size_t i = 0; i < warningsstr.length(); i++)
-    {
-      if (writedisp)
-        {
-          output << "disp('     ";
-          writedisp = false;
-        }
+  wcc.incrementWarnings(ostr.str());
 
-      if (warningsstr[i] != '\n')
-        output << warningsstr[i];
-      else
-        {
-          output << "');" << endl;
-          if (i + 1 < warningsstr.length())
-            writedisp = true;
-        }
-    }
-}
-
-int
-WarningConsolidation::countWarnings() const
-{
-  size_t p = 0;
-  int n = 0;
-  while ((p = warnings.str().find('\n', p)) != string::npos)
-    {
-      p++;
-      n++;
-    }
-  return n;
+  return wcc;
 }

src/WarningConsolidation.hh

 /*
- * Copyright © 2012-2023 Dynare Team
+ * Copyright © 2012-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,45 +21,64 @@
 #define WARNING_CONSOLIDATION_HH
 
 #include "DynareBisonLocation.hh"
+
+#include <iostream>
+#include <ostream>
 #include <sstream>
 #include <string>
 
 using namespace std;
 
-//! Stores Warnings issued by the Preprocessor
+/* Provide our implementation of operator<< with locations in DynareBison.hh. Note that the
+   following is a template specialization of the version provided in DynareBisonLocation.hh.
+
+   Ideally it should go into DynareBisonLocation.hh, but there does not seem to be a way to achieve
+   that. */
+ostream& operator<<(ostream& stream, const Dynare::location& l);
+
 class WarningConsolidation
 {
 private:
-  stringstream warnings;
-  bool no_warn;
+  const bool no_warn;
+  int num_warnings {0};
+
+  // Increases the warning counter by as many newlines as there are in the message
+  void incrementWarnings(const string& msg);
 
 public:
   explicit WarningConsolidation(bool no_warn_arg) : no_warn {no_warn_arg}
   {
   }
 
-  //! Add A Warning to the StringStream
-  friend WarningConsolidation& operator<<(WarningConsolidation& wcc, const string& warning);
-  friend WarningConsolidation& operator<<(WarningConsolidation& wcc, const Dynare::location& loc);
+  // Generic function to print something to the warning stream
+  friend WarningConsolidation& operator<<(WarningConsolidation& wcc, auto&& warning);
+
+  /* Print std::endl to the warning stream. Unfortunately, since std::endl is a template of
+     functions, it cannot be bound to the universal reference of the generic function, hence the
+     need for this specialization. */
   friend WarningConsolidation& operator<<(WarningConsolidation& wcc, ostream& (*pf)(ostream&));
 
-  void
-  addWarning(const string& w)
+  int
+  numWarnings() const
   {
-    warnings << w;
+    return num_warnings;
+  }
 };
-  void
-  addWarning(ostream& (*pf)(ostream&))
+
+WarningConsolidation&
+operator<<(WarningConsolidation& wcc, auto&& warning)
 {
-    warnings << pf;
-  };
+  if (wcc.no_warn)
+    return wcc;
 
-  //! Write Warnings to m file
-  void writeOutput(ostream& output) const;
-  //! Count warnings
-  /*! This is done in a very lousy way, by counting newlines in the
-    stringstream... */
-  int countWarnings() const;
-};
+  ostringstream ostr;
+  ostr << warning;
+
+  cerr << ostr.str();
+
+  wcc.incrementWarnings(ostr.str());
+
+  return wcc;
+}
 
 #endif

src/macro/Directives.cc

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

src/macro/Environment.cc

 /*
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -17,7 +17,9 @@
  * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
+#include <algorithm>
 #include <cassert>
+#include <ranges>
 
 #include "Environment.hh"
 #include "Expressions.hh"
@@ -54,7 +56,7 @@ Environment::getVariable(const string& name) const
   return getGlobalEnv()->getVariable(name);
 }
 
-tuple<FunctionPtr, ExpressionPtr>
+pair<FunctionPtr, ExpressionPtr>
 Environment::getFunction(const string& name) const
 {
   if (auto it = functions.find(name); it != functions.end())
@@ -101,41 +103,54 @@ Environment::isFunctionDefined(const string& name) const noexcept
 }
 
 void
-Environment::print(ostream& output, const vector<string>& vars, const optional<int>& line,
-                   bool save) const
+Environment::print(ostream& output, const vector<string>& vars, int line, bool save) const
 {
   if (!save && !variables.empty())
-    output << "Macro Variables:" << endl;
+    output << "Macro Variables (at line " << line << "):" << endl;
+
+  // For sorting the symbols in a case-insensitive way, see #128
+  auto case_insensitive_string_less = [](const string& a, const string& b) {
+    return ranges::lexicographical_compare(
+        a, b, [](char c1, char c2) { return tolower(c1) < tolower(c2); });
+  };
 
   if (vars.empty())
-    for (const auto& it : variables)
-      printVariable(output, it.first, line, save);
+    {
+      vector<string> variables_sorted;
+      ranges::copy(views::keys(variables), back_inserter(variables_sorted));
+      ranges::sort(variables_sorted, case_insensitive_string_less);
+      for (const auto& it : variables_sorted)
+        printVariable(output, it, line, save);
+    }
   else
     for (const auto& it : vars)
       if (isVariableDefined(it))
         printVariable(output, it, line, save);
 
   if (!save && !functions.empty())
-    output << "Macro Functions:" << endl;
+    output << "Macro Functions (at line " << line << "):" << endl;
 
   if (vars.empty())
-    for (const auto& it : functions)
-      printFunction(output, it.second, line, save);
+    {
+      vector<string> functions_sorted;
+      ranges::copy(views::keys(functions), back_inserter(functions_sorted));
+      ranges::sort(functions_sorted, case_insensitive_string_less);
+      for (const auto& it : functions_sorted)
+        printFunction(output, it, line, save);
+    }
   else
     for (const auto& it : vars)
       if (isFunctionDefined(it))
-        printFunction(output, functions.at(it), line, save);
+        printFunction(output, it, line, save);
 
   if (parent)
     parent->print(output, vars, line, save);
 }
 
 void
-Environment::printVariable(ostream& output, const string& name, const optional<int>& line,
-                           bool save) const
+Environment::printVariable(ostream& output, const string& name, int line, bool save) const
 {
-  assert(!save || line);
-  output << (save ? "options_.macrovars_line_" + to_string(*line) + "." : "  ") << name << " = ";
+  output << (save ? "options_.macrovars_line_" + to_string(line) + "." : "  ") << name << " = ";
   getVariable(name)->eval(const_cast<Environment&>(*this))->print(output, save);
   if (save)
     output << ";";
@@ -143,20 +158,19 @@ Environment::printVariable(ostream& output, const string& name, const optional<i
 }
 
 void
-Environment::printFunction(ostream& output, const tuple<FunctionPtr, ExpressionPtr>& function,
-                           const optional<int>& line, bool save) const
+Environment::printFunction(ostream& output, const string& name, int line, bool save) const
 {
-  assert(!save || line);
-  output << (save ? "options_.macrovars_line_" + to_string(*line) + ".function." : "  ");
+  auto [func_signature, func_body] = getFunction(name);
+  output << (save ? "options_.macrovars_line_" + to_string(line) + ".function." : "  ");
   if (save)
     {
-      get<0>(function)->printName(output);
+      func_signature->printName(output);
       output << " = '";
     }
 
-  get<0>(function)->print(output);
+  func_signature->print(output);
   output << " = ";
-  get<1>(function)->print(output);
+  func_body->print(output);
 
   if (save)
     output << "';";

src/macro/Environment.hh

 /*
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -33,7 +33,7 @@ class Environment
 private:
   const Environment* parent {nullptr};
   map<string, ExpressionPtr> variables;
-  map<string, tuple<FunctionPtr, ExpressionPtr>> functions;
+  map<string, pair<FunctionPtr, ExpressionPtr>> functions;
 
 public:
   Environment() = default;
@@ -45,7 +45,7 @@ public:
   /* The following two functions are not marked [[nodiscard]], because they are used without output
      to check whether they return an exception or not. */
   ExpressionPtr getVariable(const string& name) const; // NOLINT(modernize-use-nodiscard)
-  tuple<FunctionPtr, ExpressionPtr>                    // NOLINT(modernize-use-nodiscard)
+  pair<FunctionPtr, ExpressionPtr>                     // NOLINT(modernize-use-nodiscard)
   getFunction(const string& name) const;
   [[nodiscard]] codes::BaseType getType(const string& name) const;
   [[nodiscard]] bool isVariableDefined(const string& name) const noexcept;
@@ -55,12 +55,9 @@ public:
   {
     return isVariableDefined(name) || isFunctionDefined(name);
   }
-  void print(ostream& output, const vector<string>& vars, const optional<int>& line = nullopt,
-             bool save = false) const;
-  void printVariable(ostream& output, const string& name, const optional<int>& line,
-                     bool save) const;
-  void printFunction(ostream& output, const tuple<FunctionPtr, ExpressionPtr>& function,
-                     const optional<int>& line, bool save) const;
+  void print(ostream& output, const vector<string>& vars, int line, bool save) const;
+  void printVariable(ostream& output, const string& name, int line, bool save) const;
+  void printFunction(ostream& output, const string& name, int line, bool save) const;
   [[nodiscard]] size_t
   size() const noexcept
   {

src/macro/Expressions.cc

 /*
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -241,9 +241,9 @@ Real::normpdf(const BaseTypePtr& btp1, const BaseTypePtr& btp2) const
   auto btp22 = dynamic_pointer_cast<Real>(btp2);
   if (!btp12 || !btp22)
     throw StackTrace("Type mismatch for operands of `normpdf` operator");
-  return make_shared<Real>((1
+  return make_shared<Real>(1
                            / (btp22->value * std::sqrt(2 * numbers::pi)
-                               * std::exp(pow((value - btp12->value) / btp22->value, 2) / 2))));
+                              * std::exp(pow((value - btp12->value) / btp22->value, 2) / 2)));
 }
 
 RealPtr
@@ -254,7 +254,7 @@ Real::normcdf(const BaseTypePtr& btp1, const BaseTypePtr& btp2) const
   if (!btp12 || !btp22)
     throw StackTrace("Type mismatch for operands of `normpdf` operator");
   return make_shared<Real>(
-      (0.5 * (1 + std::erf((value - btp12->value) / btp22->value / numbers::sqrt2))));
+      0.5 * (1 + std::erf((value - btp12->value) / btp22->value / numbers::sqrt2)));
 }
 
 BaseTypePtr
@@ -314,12 +314,12 @@ String::is_equal(const BaseTypePtr& btp) const
 BoolPtr
 String::cast_bool([[maybe_unused]] Environment& env) const
 {
-  auto f = [](const char& a, const char& b) { return (tolower(a) == tolower(b)); };
+  auto f = [](const char& a, const char& b) { return tolower(a) == tolower(b); };
 
-  if (string tf = "true"; equal(value.begin(), value.end(), tf.begin(), tf.end(), f))
+  if (ranges::equal(value, "true"s, f))
     return make_shared<Bool>(true);
 
-  if (string tf = "false"; equal(value.begin(), value.end(), tf.begin(), tf.end(), f))
+  if (ranges::equal(value, "false"s, f))
     return make_shared<Bool>(false);
 
   try

src/macro/Parser.yy

 // -*- C++ -*-
 /*
- * Copyright © 2019-2023 Dynare Team
+ * Copyright © 2019-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -28,10 +28,6 @@
 %define parse.error verbose
 %define parse.trace
 
-%code requires {
-namespace macro { class Driver; }
-}
-
 %param { macro::Driver& driver }
 
 %locations
@@ -43,6 +39,9 @@ namespace macro { class Driver; }
 
 %code requires {
 #include "Directives.hh"
+
+namespace macro { class Driver; }
+
 using namespace macro;
 }
 

src/meson.build

@@ -53,6 +53,8 @@ preprocessor_src = [ 'ComputingTasks.cc',
 	             'ModelEquationBlock.cc',
 	             'WarningConsolidation.cc',
 	             'SubModel.cc',
+                     'HeterogeneityTable.cc',
+                     'HeterogeneousModel.cc',
                      'macro/Driver.cc',
 	             'macro/Environment.cc',
 	             'macro/Expressions.cc',