src/SubModel.cc

+/*
+ * Copyright © 2018-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 <algorithm>
+#include <cassert>
+#include <numeric>
+
+#include "DynamicModel.hh"
+#include "SubModel.hh"
+
+TrendComponentModelTable::TrendComponentModelTable(SymbolTable& symbol_table_arg) :
+    symbol_table {symbol_table_arg}
+{
+}
+
+void
+TrendComponentModelTable::addTrendComponentModel(string name_arg, vector<string> eqtags_arg,
+                                                 vector<string> target_eqtags_arg)
+{
+  if (isExistingTrendComponentModelName(name_arg))
+    {
+      cerr << "Error: a trend component model already exists with the name " << name_arg << endl;
+      exit(EXIT_FAILURE);
+    }
+  eqtags[name_arg] = move(eqtags_arg);
+  target_eqtags[name_arg] = move(target_eqtags_arg);
+  names.insert(move(name_arg));
+}
+
+void
+TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg,
+                                  map<string, vector<int>> target_eqnums_arg,
+                                  map<string, vector<int>> lhs_arg,
+                                  map<string, vector<expr_t>> lhs_expr_t_arg)
+{
+  eqnums = move(eqnums_arg);
+  target_eqnums = move(target_eqnums_arg);
+  lhs = move(lhs_arg);
+  lhs_expr_t = move(lhs_expr_t_arg);
+
+  for (const auto& it : eqnums)
+    {
+      vector<int> nontrend_vec;
+      for (auto eq : it.second)
+        if (ranges::find(target_eqnums[it.first], eq) == target_eqnums[it.first].end())
+          nontrend_vec.push_back(eq);
+      nontarget_eqnums[it.first] = nontrend_vec;
+    }
+
+  for (const auto& name : names)
+    {
+      vector<int> nontarget_lhs_vec, target_lhs_vec;
+      vector<int> lhsv = getLhs(name);
+      vector<int> eqnumsv = getEqNums(name);
+      for (int nontrend_it : getNonTargetEqNums(name))
+        nontarget_lhs_vec.push_back(
+            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(ranges::distance(eqnumsv.begin(), ranges::find(eqnumsv, trend_it))));
+      target_lhs[name] = target_lhs_vec;
+    }
+}
+
+void
+TrendComponentModelTable::setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg)
+{
+  rhs = move(rhs_arg);
+}
+
+void
+TrendComponentModelTable::setTargetVar(map<string, vector<optional<int>>> target_vars_arg)
+{
+  target_vars = move(target_vars_arg);
+}
+
+void
+TrendComponentModelTable::setMaxLags(map<string, vector<int>> max_lags_arg)
+{
+  max_lags = move(max_lags_arg);
+}
+
+void
+TrendComponentModelTable::setDiff(map<string, vector<bool>> diff_arg)
+{
+  diff = move(diff_arg);
+}
+
+void
+TrendComponentModelTable::setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg)
+{
+  orig_diff_var = move(orig_diff_var_arg);
+}
+
+void
+TrendComponentModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
+{
+  AR = move(AR_arg);
+}
+
+void
+TrendComponentModelTable::setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg,
+                                map<string, map<tuple<int, int>, expr_t>> A0star_arg)
+{
+  A0 = move(A0_arg);
+  A0star = move(A0star_arg);
+}
+
+const map<string, vector<string>>&
+TrendComponentModelTable::getEqTags() const
+{
+  return eqtags;
+}
+
+const vector<string>&
+TrendComponentModelTable::getEqTags(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return eqtags.at(name_arg);
+}
+
+void
+TrendComponentModelTable::checkModelName(const string& name_arg) const
+{
+  if (!isExistingTrendComponentModelName(name_arg))
+    {
+      cerr << name_arg << " is not a recognized equation tag of a trend component model equation"
+           << endl;
+      exit(EXIT_FAILURE);
+    }
+}
+
+const vector<int>&
+TrendComponentModelTable::getNonTargetLhs(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return nontarget_lhs.at(name_arg);
+}
+
+const vector<int>&
+TrendComponentModelTable::getTargetLhs(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return target_lhs.at(name_arg);
+}
+
+const vector<int>&
+TrendComponentModelTable::getLhs(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return lhs.at(name_arg);
+}
+
+const vector<expr_t>&
+TrendComponentModelTable::getLhsExprT(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return lhs_expr_t.at(name_arg);
+}
+
+const map<string, vector<string>>&
+TrendComponentModelTable::getTargetEqTags() const
+{
+  return target_eqtags;
+}
+
+const map<string, vector<int>>&
+TrendComponentModelTable::getEqNums() const
+{
+  return eqnums;
+}
+
+const map<string, vector<int>>&
+TrendComponentModelTable::getTargetEqNums() const
+{
+  return target_eqnums;
+}
+
+const vector<int>&
+TrendComponentModelTable::getTargetEqNums(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return target_eqnums.at(name_arg);
+}
+
+const map<string, vector<int>>&
+TrendComponentModelTable::getNonTargetEqNums() const
+{
+  return nontarget_eqnums;
+}
+
+const vector<int>&
+TrendComponentModelTable::getNonTargetEqNums(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return nontarget_eqnums.at(name_arg);
+}
+
+const vector<int>&
+TrendComponentModelTable::getEqNums(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return eqnums.at(name_arg);
+}
+
+const vector<int>&
+TrendComponentModelTable::getMaxLags(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return max_lags.at(name_arg);
+}
+
+int
+TrendComponentModelTable::getMaxLag(const string& name_arg) const
+{
+  int max_lag_int = 0;
+  for (auto it : getMaxLags(name_arg))
+    max_lag_int = max(max_lag_int, it);
+  return max_lag_int;
+}
+
+const vector<bool>&
+TrendComponentModelTable::getDiff(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return diff.at(name_arg);
+}
+
+void
+TrendComponentModelTable::writeOutput(const string& basename, ostream& output) const
+{
+  if (names.empty())
+    return;
+
+  const filesystem::path filename {DataTree::packageDir(basename) / "trend_component_ar_a0.m"};
+  ofstream ar_ec_output {filename, ios::out | ios::binary};
+  if (!ar_ec_output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  ar_ec_output << "function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
+               << "%function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
+               << "% File automatically generated by the Dynare preprocessor" << endl
+               << endl;
+
+  for (const auto& name : names)
+    {
+      output << "M_.trend_component." << name << ".model_name = '" << name << "';" << endl
+             << "M_.trend_component." << name << ".eqtags = {";
+      for (const auto& it : eqtags.at(name))
+        output << "'" << it << "'; ";
+      output << "};" << endl << "M_.trend_component." << name << ".eqn = [";
+      for (auto it : eqnums.at(name))
+        output << it + 1 << " ";
+      output << "];" << endl << "M_.trend_component." << name << ".targets = [";
+      for (auto it : eqnums.at(name))
+        if (ranges::find(target_eqnums.at(name), it) == target_eqnums.at(name).end())
+          output << "false ";
+        else
+          output << "true ";
+      output << "];" << endl << "M_.trend_component." << name << ".lhs = [";
+      for (auto it : lhs.at(name))
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl << "M_.trend_component." << name << ".max_lag = [";
+      for (auto it : max_lags.at(name))
+        output << it << " ";
+      output << "];" << endl << "M_.trend_component." << name << ".diff = [";
+      for (bool it : diff.at(name))
+        output << boolalpha << it << " ";
+      output << "];" << endl << "M_.trend_component." << name << ".orig_diff_var = [";
+      for (const auto& it : orig_diff_var.at(name))
+        output << (it ? symbol_table.getTypeSpecificID(*it) + 1 : -1) << " ";
+      output << "];" << endl << "M_.trend_component." << name << ".nonstationary = [";
+      for (size_t i = 0; i < diff.at(name).size(); i++)
+        output << "true ";
+      output << "];" << endl;
+      for (int i {1}; const auto& it : rhs.at(name))
+        {
+          output << "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
+          for (auto [var, lag] : it)
+            output << symbol_table.getTypeSpecificID(var) + 1 << " ";
+          output << "];" << endl
+                 << "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
+          for (auto [var, lag] : it)
+            output << lag << " ";
+          output << "];" << endl;
+
+          i++;
+        }
+      output << "M_.trend_component." << name << ".target_vars = [";
+      for (const optional<int>& it : target_vars.at(name))
+        output << (it ? symbol_table.getTypeSpecificID(*it) + 1 : -1) << " ";
+      output << "];" << endl;
+
+      vector<string> target_eqtags_vec = target_eqtags.at(name);
+      output << "M_.trend_component." << name << ".target_eqtags = {";
+      for (const auto& it : target_eqtags_vec)
+        output << "'" << it << "';";
+      output << "};" << endl;
+
+      vector<string> eqtags_vec = eqtags.at(name);
+      output << "M_.trend_component." << name << ".target_eqn = [";
+      for (const auto& it : target_eqtags_vec)
+        output << ranges::distance(eqtags_vec.begin(), ranges::find(eqtags_vec, it)) + 1 << " ";
+      output << "];" << endl;
+
+      vector<int> target_lhs_vec = getTargetLhs(name);
+      vector<int> nontarget_lhs_vec = getNonTargetLhs(name);
+
+      ar_ec_output << "if strcmp(model_name, '" << name << "')" << endl
+                   << "    % AR" << endl
+                   << "    AR = zeros(" << nontarget_lhs_vec.size() << ", "
+                   << nontarget_lhs_vec.size() << ", " << getMaxLag(name) << ");" << endl;
+      for (const auto& [key, expr] : AR.at(name))
+        {
+          auto [eqn, lag, lhs_symb_id] = key;
+          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;
+        }
+
+      ar_ec_output << endl
+                   << "    % A0" << endl
+                   << "    A0 = zeros(" << nontarget_lhs_vec.size() << ", "
+                   << nontarget_lhs_vec.size() << ");" << endl;
+      for (const auto& [key, expr] : A0.at(name))
+        {
+          auto [eqn, colidx] = key;
+          ar_ec_output << "    A0(" << eqn + 1 << ", " << colidx + 1 << ") = ";
+          expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_ec_output << ";" << endl;
+        }
+
+      ar_ec_output << endl
+                   << "    % A0star" << endl
+                   << "    A0star = zeros(" << nontarget_lhs_vec.size() << ", "
+                   << target_lhs_vec.size() << ");" << endl;
+      for (const auto& [key, expr] : A0star.at(name))
+        {
+          auto [eqn, colidx] = key;
+          ar_ec_output << "    A0star(" << eqn + 1 << ", " << colidx + 1 << ") = ";
+          expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_ec_output << ";" << endl;
+        }
+
+      ar_ec_output << "    return" << endl << "end" << endl << endl;
+    }
+  ar_ec_output << "error([model_name ' is not a valid trend_component_model name'])" << endl
+               << "end" << endl;
+  ar_ec_output.close();
+}
+
+void
+TrendComponentModelTable::writeJsonOutput(ostream& output) const
+{
+  for (bool printed_something {false}; const auto& name : names)
+    {
+      if (exchange(printed_something, true))
+        output << ", ";
+      output << R"({"statementName": "trend_component_model",)"
+             << R"("model_name": ")" << name << R"(",)"
+             << R"("eqtags": [)";
+      for (bool printed_something2 {false}; const auto& it : eqtags.at(name))
+        {
+          if (exchange(printed_something2, true))
+            output << ", ";
+          output << R"(")" << it << R"(")";
+        }
+      output << R"(], "target_eqtags": [)";
+      for (bool printed_something2 {false}; const auto& it : target_eqtags.at(name))
+        {
+          if (exchange(printed_something2, true))
+            output << ", ";
+          output << R"(")" << it << R"(")";
+        }
+      output << "]}";
+    }
+}
+
+VarModelTable::VarModelTable(SymbolTable& symbol_table_arg) : symbol_table {symbol_table_arg}
+{
+}
+
+void
+VarModelTable::addVarModel(string name_arg, bool structural_arg, vector<string> eqtags_arg)
+{
+  if (isExistingVarModelName(name_arg))
+    {
+      cerr << "Error: a VAR model already exists with the name " << name_arg << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  structural[name_arg] = structural_arg;
+  eqtags[name_arg] = move(eqtags_arg);
+  names.insert(move(name_arg));
+}
+
+void
+VarModelTable::writeOutput(const string& basename, ostream& output) const
+{
+  if (names.empty())
+    return;
+
+  const filesystem::path filename {DataTree::packageDir(basename) / "varmatrices.m"};
+  ofstream ar_output {filename, ios::out | ios::binary};
+  if (!ar_output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  ar_output << "function [ar, a0, constants] = varmatrices(model_name, params, reducedform)" << endl
+            << "% File automatically generated by the Dynare preprocessor" << endl
+            << endl
+            << "if nargin<3" << endl
+            << "    reducedform = false;" << endl
+            << "end" << endl
+            << endl;
+
+  for (const auto& name : names)
+    {
+      output << "M_.var." << name << ".model_name = '" << name << "';" << endl
+             << "M_.var." << name << ".structural = " << boolalpha << structural.at(name) << ";"
+             << endl
+             << "M_.var." << name << ".eqtags = {";
+      for (const auto& it : eqtags.at(name))
+        output << "'" << it << "'; ";
+      output << "};" << endl << "M_.var." << name << ".eqn = [";
+      for (auto it : eqnums.at(name))
+        output << it + 1 << " ";
+      output << "];" << endl << "M_.var." << name << ".lhs = [";
+      for (auto it : lhs.at(name))
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl << "M_.var." << name << ".max_lag = [";
+      for (auto it : max_lags.at(name))
+        output << it << " ";
+      output << "];" << endl << "M_.var." << name << ".diff = [";
+      for (bool it : diff.at(name))
+        output << boolalpha << it << " ";
+      output << "];" << endl
+             << "M_.var." << name << ".nonstationary = M_.var." << name << ".diff;" << endl
+             << "M_.var." << name << ".orig_diff_var = [";
+      for (const auto& it : orig_diff_var.at(name))
+        output << (it ? symbol_table.getTypeSpecificID(*it) + 1 : -1) << " ";
+      output << "];" << endl;
+      for (int i {1}; const auto& it : rhs.at(name))
+        {
+          output << "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
+          for (auto [var, lag] : it)
+            output << symbol_table.getTypeSpecificID(var) + 1 << " ";
+          output << "];" << endl << "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
+          for (auto [var, lag] : it)
+            output << lag << " ";
+          output << "];" << endl;
+
+          i++;
+        }
+
+      vector<int> lhs = getLhsOrigIds(name);
+      ar_output << "if strcmp(model_name, '" << name << "')" << endl
+                << "    ar = zeros(" << lhs.size() << ", " << lhs.size() << ", " << getMaxLag(name)
+                << ");" << endl;
+      for (const auto& [key, expr] : AR.at(name))
+        {
+          auto [eqn, lag, lhs_symb_id] = key;
+          int colidx
+              = 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;
+        }
+      ar_output << "    if nargout>1" << endl << "        a0 = eye(" << lhs.size() << ");" << endl;
+      for (const auto& [key, expr] : A0.at(name))
+        {
+          auto [eqn, lhs_symb_id] = key;
+          int colidx
+              = static_cast<int>(ranges::distance(lhs.begin(), ranges::find(lhs, lhs_symb_id)));
+          if (eqn != colidx)
+            {
+              ar_output << "        a0(" << eqn + 1 << "," << colidx + 1 << ") = ";
+              expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
+              ar_output << ";" << endl;
+            }
+        }
+      ar_output << "        if reducedform" << endl
+                << "            for i=1:" << getMaxLag(name) << endl
+                << "                ar(:,:,i) = a0\\ar(:,:,i);" << endl
+                << "            end" << endl
+                << "            if nargout<3" << endl
+                << "                a0 = eye(" << lhs.size() << ");" << endl
+                << "            end" << endl
+                << "        end" << endl
+                << "        if nargout>2" << endl
+                << "            constants = zeros(" << lhs.size() << ",1);" << endl;
+      for (auto [eqn, expr] : constants.at(name))
+        {
+          ar_output << "            constants(" << eqn + 1 << ") = ";
+          expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_output << ";" << endl;
+        }
+      ar_output << "        end" << endl
+                << "        if reducedform" << endl
+                << "            constants = a0\\constants;" << endl
+                << "            a0 = eye(" << lhs.size() << ");" << endl
+                << "        end" << endl
+                << "    end" << endl
+                << "    return" << endl
+                << "end" << endl
+                << endl;
+    }
+  ar_output << "error('%s is not a valid var_model name', model_name)" << endl;
+  ar_output.close();
+}
+
+void
+VarModelTable::writeJsonOutput(ostream& output) const
+{
+  for (bool printed_something {false}; const auto& name : names)
+    {
+      if (exchange(printed_something, true))
+        output << ", ";
+      output << R"({"statementName": "var_model",)"
+             << R"("model_name": ")" << name << R"(",)"
+             << R"("eqtags": [)";
+      for (bool printed_something2 {false}; const auto& it : eqtags.at(name))
+        {
+          if (exchange(printed_something2, true))
+            output << ", ";
+          output << R"(")" << it << R"(")";
+        }
+      output << "]}";
+    }
+}
+
+const map<string, bool>&
+VarModelTable::getStructural() const
+{
+  return structural;
+}
+
+const map<string, vector<string>>&
+VarModelTable::getEqTags() const
+{
+  return eqtags;
+}
+
+const vector<string>&
+VarModelTable::getEqTags(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return eqtags.at(name_arg);
+}
+
+void
+VarModelTable::checkModelName(const string& name_arg) const
+{
+  if (!isExistingVarModelName(name_arg))
+    {
+      cerr << name_arg << " is not a recognized equation tag of a VAR model equation" << endl;
+      exit(EXIT_FAILURE);
+    }
+}
+
+void
+VarModelTable::setEqNums(map<string, vector<int>> eqnums_arg)
+{
+  eqnums = move(eqnums_arg);
+}
+
+void
+VarModelTable::setLhs(map<string, vector<int>> lhs_arg)
+{
+  lhs = move(lhs_arg);
+  for (const auto& it : lhs)
+    {
+      vector<int> lhsvec;
+      for (auto ids : it.second)
+        {
+          int lhs_last_orig_symb_id = ids;
+          int lhs_orig_symb_id = ids;
+          if (symbol_table.isDiffAuxiliaryVariable(lhs_orig_symb_id))
+            try
+              {
+                lhs_last_orig_symb_id = lhs_orig_symb_id;
+                lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_orig_symb_id);
+              }
+            catch (...)
+              {
+              }
+          lhsvec.emplace_back(lhs_last_orig_symb_id);
+        }
+      lhs_orig_symb_ids[it.first] = lhsvec;
+    }
+}
+
+void
+VarModelTable::setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg)
+{
+  rhs = move(rhs_arg);
+}
+
+void
+VarModelTable::setLhsExprT(map<string, vector<expr_t>> lhs_expr_t_arg)
+{
+  lhs_expr_t = move(lhs_expr_t_arg);
+}
+
+const map<string, vector<int>>&
+VarModelTable::getEqNums() const
+{
+  return eqnums;
+}
+
+const vector<bool>&
+VarModelTable::getDiff(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return diff.at(name_arg);
+}
+
+const vector<int>&
+VarModelTable::getEqNums(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return eqnums.at(name_arg);
+}
+
+void
+VarModelTable::setMaxLags(map<string, vector<int>> max_lags_arg)
+{
+  max_lags = move(max_lags_arg);
+}
+
+void
+VarModelTable::setDiff(map<string, vector<bool>> diff_arg)
+{
+  diff = move(diff_arg);
+}
+
+void
+VarModelTable::setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg)
+{
+  orig_diff_var = move(orig_diff_var_arg);
+}
+
+void
+VarModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
+{
+  AR = move(AR_arg);
+}
+
+void
+VarModelTable::setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg)
+{
+  A0 = move(A0_arg);
+}
+
+void
+VarModelTable::setConstants(map<string, map<int, expr_t>> constants_arg)
+{
+  constants = move(constants_arg);
+}
+
+const vector<int>&
+VarModelTable::getMaxLags(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return max_lags.at(name_arg);
+}
+
+int
+VarModelTable::getMaxLag(const string& name_arg) const
+{
+  vector<int> maxlags {getMaxLags(name_arg)};
+  return reduce(maxlags.begin(), maxlags.end(), 0, [](int a, int b) { return max(a, b); });
+}
+
+const vector<int>&
+VarModelTable::getLhs(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return lhs.at(name_arg);
+}
+
+const vector<int>&
+VarModelTable::getLhsOrigIds(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return lhs_orig_symb_ids.at(name_arg);
+}
+
+const vector<set<pair<int, int>>>&
+VarModelTable::getRhs(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return rhs.at(name_arg);
+}
+
+const vector<expr_t>&
+VarModelTable::getLhsExprT(const string& name_arg) const
+{
+  checkModelName(name_arg);
+  return lhs_expr_t.at(name_arg);
+}
+
+VarExpectationModelTable::VarExpectationModelTable(SymbolTable& symbol_table_arg) :
+    symbol_table {symbol_table_arg}
+{
+}
+
+void
+VarExpectationModelTable::addVarExpectationModel(string name_arg, expr_t expression_arg,
+                                                 string aux_model_name_arg, string horizon_arg,
+                                                 expr_t discount_arg, int time_shift_arg)
+{
+  if (isExistingVarExpectationModelName(name_arg))
+    {
+      cerr << "Error: a var_expectation_model already exists with the name " << name_arg << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  expression[name_arg] = expression_arg;
+  aux_model_name[name_arg] = move(aux_model_name_arg);
+  horizon[name_arg] = move(horizon_arg);
+  discount[name_arg] = discount_arg;
+  time_shift[name_arg] = time_shift_arg;
+  names.insert(move(name_arg));
+}
+
+bool
+VarExpectationModelTable::isExistingVarExpectationModelName(const string& name_arg) const
+{
+  return names.contains(name_arg);
+}
+
+bool
+VarExpectationModelTable::empty() const
+{
+  return names.empty();
+}
+
+void
+VarExpectationModelTable::writeOutput(ostream& output) const
+{
+  for (const auto& name : names)
+    {
+      string mstruct = "M_.var_expectation." + name;
+      output << mstruct << ".auxiliary_model_name = '" << aux_model_name.at(name) << "';" << endl
+             << mstruct << ".horizon = " << horizon.at(name) << ';' << endl
+             << mstruct << ".time_shift = " << time_shift.at(name) << ';' << endl;
+
+      auto& vpc = vars_params_constants.at(name);
+      if (!vpc.size())
+        {
+          cerr << "ERROR: VarExpectationModelStatement::writeOutput: matchExpression() has not "
+                  "been called"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      ostringstream vars_list, params_list, constants_list;
+      for (bool printed_something {false}; const auto& [variable_id, param_id, constant] : vpc)
+        {
+          if (exchange(printed_something, true))
+            {
+              vars_list << ", ";
+              params_list << ", ";
+              constants_list << ", ";
+            }
+          vars_list << symbol_table.getTypeSpecificID(variable_id) + 1;
+          if (param_id)
+            params_list << symbol_table.getTypeSpecificID(*param_id) + 1;
+          else
+            params_list << "NaN";
+          constants_list << constant;
+        }
+      output << mstruct << ".expr.vars = [ " << vars_list.str() << " ];" << endl
+             << mstruct << ".expr.params = [ " << params_list.str() << " ];" << endl
+             << mstruct << ".expr.constants = [ " << constants_list.str() << " ];" << endl;
+
+      if (auto disc_var = dynamic_cast<const VariableNode*>(discount.at(name)); disc_var)
+        output << mstruct << ".discount_index = " << disc_var->getTypeSpecificID() + 1 << ';'
+               << endl;
+      else
+        {
+          output << mstruct << ".discount_value = ";
+          discount.at(name)->writeOutput(output);
+          output << ';' << endl;
+        }
+      output << mstruct << ".param_indices = [ ";
+      for (int param_id : aux_param_symb_ids.at(name))
+        output << symbol_table.getTypeSpecificID(param_id) + 1 << ' ';
+      output << "];" << endl;
+    }
+}
+
+void
+VarExpectationModelTable::substituteUnaryOpsInExpression(const lag_equivalence_table_t& nodes,
+                                                         ExprNode::subst_table_t& subst_table,
+                                                         vector<BinaryOpNode*>& neweqs)
+{
+  for (const auto& name : names)
+    expression[name] = expression[name]->substituteUnaryOpNodes(nodes, subst_table, neweqs);
+}
+
+void
+VarExpectationModelTable::substituteDiffNodesInExpression(const lag_equivalence_table_t& nodes,
+                                                          ExprNode::subst_table_t& subst_table,
+                                                          vector<BinaryOpNode*>& neweqs)
+{
+  for (const auto& name : names)
+    expression[name] = expression[name]->substituteDiff(nodes, subst_table, neweqs);
+}
+
+void
+VarExpectationModelTable::transformPass(ExprNode::subst_table_t& diff_subst_table,
+                                        DynamicModel& dynamic_model,
+                                        const VarModelTable& var_model_table,
+                                        const TrendComponentModelTable& trend_component_model_table)
+{
+  map<string, expr_t> var_expectation_subst_table;
+
+  for (const auto& name : names)
+    {
+      // Collect information about the auxiliary model
+
+      int max_lag;
+      vector<int> lhs;
+      if (var_model_table.isExistingVarModelName(aux_model_name[name]))
+        {
+          max_lag = var_model_table.getMaxLag(aux_model_name[name]);
+          lhs = var_model_table.getLhs(aux_model_name[name]);
+        }
+      else if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name[name]))
+        {
+          max_lag = trend_component_model_table.getMaxLag(aux_model_name[name]) + 1;
+          lhs = dynamic_model.getUndiffLHSForPac(aux_model_name[name], diff_subst_table);
+        }
+      else
+        {
+          cerr << "ERROR: var_expectation_model " << name
+               << " refers to nonexistent auxiliary model " << aux_model_name[name] << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      // Match the linear combination in the expression option
+      try
+        {
+          auto vpc = expression[name]->matchLinearCombinationOfVariables();
+          for (const auto& [variable_id, lag, param_id, constant] : vpc)
+            {
+              if (lag != 0)
+                throw ExprNode::MatchFailureException {"lead/lags are not allowed"};
+              if (symbol_table.getType(variable_id) != SymbolType::endogenous)
+                throw ExprNode::MatchFailureException {"Variable is not an endogenous"};
+              vars_params_constants[name].emplace_back(variable_id, param_id, constant);
+            }
+        }
+      catch (ExprNode::MatchFailureException& e)
+        {
+          cerr << "ERROR: expression in var_expectation_model " << name
+               << " is not of the expected form: " << e.message << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      /* Create auxiliary parameters and the expression to be substituted into
+         the var_expectations statement */
+      expr_t subst_expr = dynamic_model.Zero;
+      if (var_model_table.isExistingVarModelName(aux_model_name[name]))
+        {
+          /* If the auxiliary model is a VAR, add a parameter corresponding to
+             the constant. */
+          string constant_param_name = "var_expectation_model_" + name + "_constant";
+          int constant_param_id
+              = symbol_table.addSymbol(constant_param_name, SymbolType::parameter);
+          aux_param_symb_ids[name].push_back(constant_param_id);
+          subst_expr
+              = dynamic_model.AddPlus(subst_expr, dynamic_model.AddVariable(constant_param_id));
+        }
+      for (int lag = 0; lag < max_lag; lag++)
+        for (auto variable : lhs)
+          {
+            string param_name = "var_expectation_model_" + name + '_'
+                                + symbol_table.getName(variable) + '_' + to_string(lag);
+            int new_param_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+            aux_param_symb_ids[name].push_back(new_param_id);
+
+            subst_expr = dynamic_model.AddPlus(
+                subst_expr, dynamic_model.AddTimes(
+                                dynamic_model.AddVariable(new_param_id),
+                                dynamic_model.AddVariable(variable, -lag + time_shift[name])));
+          }
+
+      if (var_expectation_subst_table.contains(name))
+        {
+          cerr << "ERROR: model name '" << name
+               << "' is used by several var_expectation_model statements" << endl;
+          exit(EXIT_FAILURE);
+        }
+      var_expectation_subst_table[name] = subst_expr;
+    }
+
+  // Actually substitute var_expectation statements
+  dynamic_model.substituteVarExpectation(var_expectation_subst_table);
+  /* At this point, we know that all var_expectation operators have been
+     substituted, because of the error check performed in
+     VarExpectationNode::substituteVarExpectation(). */
+}
+
+void
+VarExpectationModelTable::writeJsonOutput(ostream& output) const
+{
+  for (bool printed_something {false}; const auto& name : names)
+    {
+      if (exchange(printed_something, true))
+        output << ", ";
+      output << R"({"statementName": "var_expectation_model",)"
+             << R"("model_name": ")" << name << R"(", )"
+             << R"("expression": ")";
+      expression.at(name)->writeOutput(output);
+      output << R"(", )"
+             << R"("auxiliary_model_name": ")" << aux_model_name.at(name) << R"(", )"
+             << R"("horizon": ")" << horizon.at(name) << R"(", )"
+             << R"("discount": ")";
+      discount.at(name)->writeOutput(output);
+      output << R"(", )"
+             << R"("time_shift": )" << time_shift.at(name) << R"(})";
+    }
+}
+
+PacModelTable::PacModelTable(SymbolTable& symbol_table_arg) : symbol_table {symbol_table_arg}
+{
+}
+
+void
+PacModelTable::addPacModel(string name_arg, string aux_model_name_arg, string discount_arg,
+                           expr_t growth_arg, string auxname_arg, PacTargetKind kind_arg)
+{
+  if (isExistingPacModelName(name_arg))
+    {
+      cerr << "Error: a PAC model already exists with the name " << name_arg << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_model_name[name_arg] = move(aux_model_name_arg);
+  discount[name_arg] = move(discount_arg);
+  growth[name_arg] = growth_arg;
+  original_growth[name_arg] = growth_arg;
+  auxname[name_arg] = move(auxname_arg);
+  kind[name_arg] = kind_arg;
+  names.insert(move(name_arg));
+}
+
+bool
+PacModelTable::isExistingPacModelName(const string& name_arg) const
+{
+  return names.contains(name_arg);
+}
+
+bool
+PacModelTable::empty() const
+{
+  return names.empty();
+}
+
+void
+PacModelTable::checkPass(ModFileStructure& mod_file_struct)
+{
+  for (auto& [name, gv] : growth)
+    if (gv)
+      {
+        if (target_info.contains(name))
+          {
+            cerr << "ERROR: for PAC model '" << name
+                 << "', it is not possible to declare a 'growth' option in the 'pac_model' command "
+                    "when there is also a 'pac_target_info' block"
+                 << endl;
+            exit(EXIT_FAILURE);
+          }
+        gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+      }
+
+  for (auto& [name, auxn] : auxname)
+    if (!auxn.empty() && target_info.contains(name))
+      {
+        cerr << "ERROR: for PAC model '" << name
+             << "', it is not possible to declare an 'auxname' option in the 'pac_model' command "
+                "when there is also a 'pac_target_info' block"
+             << endl;
+        exit(EXIT_FAILURE);
+      }
+
+  for (auto& [name, k] : kind)
+    if (k != PacTargetKind::unspecified)
+      {
+        if (target_info.contains(name))
+          {
+            cerr << "ERROR: for PAC model '" << name
+                 << "', it is not possible to declare a 'kind' option in the 'pac_model' command "
+                    "when there is also a 'pac_target_info' block"
+                 << endl;
+            exit(EXIT_FAILURE);
+          }
+        if (aux_model_name[name].empty())
+          {
+            cerr << "ERROR: for PAC model '" << name
+                 << "', it is not possible to declare a 'kind' option in the 'pac_model' command "
+                    "since this is a MCE model"
+                 << endl;
+            exit(EXIT_FAILURE);
+          }
+      }
+
+  for (const auto& [name, ti] : target_info)
+    for (auto& [expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv,
+                gv_info] : get<2>(ti))
+      if (gv)
+        gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+
+  for (const auto& [name, ti] : target_info)
+    {
+      auto& [target, auxname_target_nonstationary, components] = ti;
+      if (!target)
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name
+               << ")' is missing the 'target' statement" << endl;
+          exit(EXIT_FAILURE);
+        }
+      if (auxname_target_nonstationary.empty())
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name
+               << ")' is missing the 'auxname_target_nonstationary' statement" << endl;
+          exit(EXIT_FAILURE);
+        }
+      int nonstationary_nb = 0;
+      for (auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                  h_indices, original_growth_component, growth_component_info] : components)
+        {
+          if (auxname.empty())
+            {
+              cerr << "ERROR: the block 'pac_target_info(" << name
+                   << ")' is missing the 'auxname' statement in some 'component'" << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::unspecified)
+            {
+              cerr << "ERROR: the block 'pac_target_info(" << name
+                   << ")' is missing the 'kind' statement in some 'component'" << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::ll && growth_component)
+            {
+              cerr << "ERROR: in the block 'pac_target_info(" << name
+                   << ")', a component of 'kind ll' (i.e. stationary) has a 'growth' option. This "
+                      "is not permitted."
+                   << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (kind == PacTargetKind::dd || kind == PacTargetKind::dl)
+            nonstationary_nb++;
+        }
+      if (!nonstationary_nb)
+        {
+          cerr << "ERROR: the block 'pac_target_info(" << name
+               << ")' must contain at least one nonstationary component (i.e. of 'kind' equal to "
+                  "either 'dd' or 'dl')."
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
+}
+
+void
+PacModelTable::substituteUnaryOpsInGrowth(const lag_equivalence_table_t& nodes,
+                                          ExprNode::subst_table_t& subst_table,
+                                          vector<BinaryOpNode*>& neweqs)
+{
+  for (auto& [name, gv] : growth)
+    if (gv)
+      gv = gv->substituteUnaryOpNodes(nodes, subst_table, neweqs);
+
+  for (auto& [name, ti] : target_info)
+    for (auto& [expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv,
+                gv_info] : get<2>(ti))
+      if (gv)
+        gv = gv->substituteUnaryOpNodes(nodes, subst_table, neweqs);
+}
+
+void
+PacModelTable::findDiffNodesInGrowth(lag_equivalence_table_t& diff_nodes) const
+{
+  for (auto& [name, gv] : growth)
+    if (gv)
+      gv->findDiffNodes(diff_nodes);
+
+  for (const auto& [name, ti] : target_info)
+    for (auto& [expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv,
+                gv_info] : get<2>(ti))
+      if (gv)
+        gv->findDiffNodes(diff_nodes);
+}
+
+void
+PacModelTable::substituteDiffNodesInGrowth(const lag_equivalence_table_t& diff_nodes,
+                                           ExprNode::subst_table_t& diff_subst_table,
+                                           vector<BinaryOpNode*>& neweqs)
+{
+  for (auto& [name, gv] : growth)
+    if (gv)
+      gv = gv->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+
+  for (auto& [name, ti] : target_info)
+    for (auto& [expr, gv, auxname, kind, coeff, growth_neutrality_param, h_indices, original_gv,
+                gv_info] : get<2>(ti))
+      if (gv)
+        gv = gv->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+}
+
+void
+PacModelTable::transformPass(const lag_equivalence_table_t& unary_ops_nodes,
+                             ExprNode::subst_table_t& unary_ops_subst_table,
+                             const lag_equivalence_table_t& diff_nodes,
+                             ExprNode::subst_table_t& diff_subst_table, DynamicModel& dynamic_model,
+                             const VarModelTable& var_model_table,
+                             const TrendComponentModelTable& trend_component_model_table)
+{
+  // model name → expression for pac_expectation
+  map<string, expr_t> pac_expectation_substitution;
+
+  for (const auto& name : names)
+    {
+      /* Fill the growth_info structure.
+         Cannot be done in an earlier pass since growth terms can be
+         transformed by DynamicModel::substituteDiff(). */
+      if (growth[name])
+        try
+          {
+            growth_info[name] = growth[name]->matchLinearCombinationOfVariablesPlusConstant();
+          }
+        catch (ExprNode::MatchFailureException& e)
+          {
+            cerr << "ERROR: PAC growth must be a linear combination of variables" << endl;
+            exit(EXIT_FAILURE);
+          }
+
+      // Perform transformations for the pac_target_info block (if any)
+      if (target_info.contains(name))
+        {
+          // Substitute unary ops and diffs in the target…
+          expr_t& target = get<0>(target_info[name]);
+          vector<BinaryOpNode*> neweqs;
+          target = target->substituteUnaryOpNodes(unary_ops_nodes, unary_ops_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr
+                  << "ERROR: the 'target' expression of 'pac_target_info(" << name
+                  << ")' contains a variable with a unary operator that is not present in the model"
+                  << endl;
+              exit(EXIT_FAILURE);
+            }
+          target = target->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr << "ERROR: the 'target' expression of 'pac_target_info(" << name
+                   << ")' contains a diff'd variable that is not present in the model" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // …and in component expressions
+          auto& components = get<2>(target_info[name]);
+          for (auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                      h_indices, original_growth_component, growth_component_info] : components)
+            {
+              component = component->substituteUnaryOpNodes(unary_ops_nodes, unary_ops_subst_table,
+                                                            neweqs);
+              if (neweqs.size() > 0)
+                {
+                  cerr << "ERROR: a 'component' expression of 'pac_target_info(" << name
+                       << ")' contains a variable with a unary operator that is not present in the "
+                          "model"
+                       << endl;
+                  exit(EXIT_FAILURE);
+                }
+              component = component->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+              if (neweqs.size() > 0)
+                {
+                  cerr << "ERROR: a 'component' expression of 'pac_target_info(" << name
+                       << ")' contains a diff'd variable that is not present in the model" << endl;
+                  exit(EXIT_FAILURE);
+                }
+            }
+
+          /* Fill the growth_info structure.
+             Cannot be done in an earlier pass since growth terms can be
+             transformed by DynamicModel::substituteDiff(). */
+          for (auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                      h_indices, original_growth_component, growth_component_info] : components)
+            {
+              if (growth_component)
+                try
+                  {
+                    growth_component_info
+                        = growth_component->matchLinearCombinationOfVariablesPlusConstant();
+                  }
+                catch (ExprNode::MatchFailureException& e)
+                  {
+                    cerr << "ERROR: PAC growth must be a linear combination of variables" << endl;
+                    exit(EXIT_FAILURE);
+                  }
+            }
+
+          // Identify the model equation defining the target
+          expr_t target_expr;
+          try
+            {
+              target_expr = dynamic_model.getRHSFromLHS(target);
+            }
+          catch (ExprNode::MatchFailureException)
+            {
+              cerr << "ERROR: there is no equation whose LHS is equal to the 'target' of "
+                      "'pac_target_info("
+                   << name << ")'" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // Substitute unary ops and diffs in that equation, before parsing (see dynare#1837)
+          target_expr
+              = target_expr->substituteUnaryOpNodes(unary_ops_nodes, unary_ops_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr
+                  << "ERROR: the equation defining the target of 'pac_target_info(" << name
+                  << ")' contains a variable with a unary operator that is not present in the model"
+                  << endl;
+              exit(EXIT_FAILURE);
+            }
+          target_expr = target_expr->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+          if (neweqs.size() > 0)
+            {
+              cerr << "ERROR: the equation defining the target of 'pac_target_info(" << name
+                   << ")' contains a diff'd variable that is not present in the model" << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // Parse that model equation
+          vector<pair<int, expr_t>> terms;
+          expr_t constant;
+          try
+            {
+              tie(terms, constant) = target_expr->matchLinearCombinationOfEndogenousWithConstant();
+            }
+          catch (ExprNode::MatchFailureException)
+            {
+              cerr << "ERROR: the model equation defining the 'target' of 'pac_target_info(" << name
+                   << ")' is not of the right form (should be a linear combination of endogenous "
+                      "variables)"
+                   << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          // Associate the coefficients of the linear combination with the right components
+          for (auto [var, coeff] : terms)
+            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;
+            else
+              {
+                cerr << "ERROR: the model equation defining the 'target' of 'pac_target_info("
+                     << name << ")' contains a variable (" << symbol_table.getName(var)
+                     << ") that is not declared as a 'component'" << endl;
+                exit(EXIT_FAILURE);
+              }
+
+          // Verify that all declared components appear in that equation
+          for (const auto& [component, growth_component, auxname, kind, coeff,
+                            growth_neutrality_param, h_indices, original_growth_component,
+                            growth_component_info] : components)
+            if (!coeff)
+              {
+                cerr << "ERROR: a 'component' of 'pac_target_info(" << name
+                     << ")' does not appear in the model equation defining the 'target'" << endl;
+                exit(EXIT_FAILURE);
+              }
+
+          /* Add the variable and equation defining the stationary part of the
+             target. Note that it includes the constant. */
+          expr_t yns = constant;
+          for (const auto& [component, growth_component, auxname, kind, coeff,
+                            growth_neutrality_param, h_indices, original_growth_component,
+                            growth_component_info] : components)
+            if (kind != PacTargetKind::ll)
+              yns = dynamic_model.AddPlus(yns, dynamic_model.AddTimes(coeff, component));
+          int target_nonstationary_id
+              = symbol_table.addPacTargetNonstationaryAuxiliaryVar(get<1>(target_info[name]), yns);
+          expr_t neweq
+              = dynamic_model.AddEqual(dynamic_model.AddVariable(target_nonstationary_id), yns);
+          dynamic_model.addEquation(neweq, nullopt);
+          dynamic_model.addAuxEquation(neweq);
+
+          /* Perform the substitution of the pac_target_nonstationary operator.
+             This needs to be done here, otherwise
+             DynamicModel::analyzePacEquationStructure() will not be able to
+             identify the error-correction part */
+          dynamic_model.substitutePacTargetNonstationary(
+              name, dynamic_model.AddVariable(target_nonstationary_id, -1));
+        }
+
+      // Collect some information about PAC models
+      int max_lag;
+      if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name[name]))
+        {
+          aux_model_type[name] = "trend_component";
+          max_lag = trend_component_model_table.getMaxLag(aux_model_name[name]) + 1;
+          lhs[name] = dynamic_model.getUndiffLHSForPac(aux_model_name[name], diff_subst_table);
+        }
+      else if (var_model_table.isExistingVarModelName(aux_model_name[name]))
+        {
+          aux_model_type[name] = "var";
+          max_lag = var_model_table.getMaxLag(aux_model_name[name]);
+          lhs[name] = var_model_table.getLhs(aux_model_name[name]);
+        }
+      else if (aux_model_name[name].empty())
+        max_lag = 0;
+      else
+        {
+          cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model"
+               << endl;
+          exit(EXIT_FAILURE);
+        }
+      dynamic_model.analyzePacEquationStructure(name, eq_name, equation_info);
+
+      // Declare parameter for growth neutrality correction
+      if (growth[name])
+        {
+          string param_name = name + "_pac_growth_neutrality_correction";
+          try
+            {
+              int param_idx = symbol_table.addSymbol(param_name, SymbolType::parameter);
+              growth_neutrality_params[name] = param_idx;
+            }
+          catch (SymbolTable::AlreadyDeclaredException)
+            {
+              cerr << "The variable/parameter '" << param_name
+                   << "' conflicts with the auxiliary parameter that will be generated for the "
+                      "growth neutrality correction of the '"
+                   << name << "' PAC model. Please rename that parameter." << endl;
+              exit(EXIT_FAILURE);
+            }
+        }
+
+      // Compute the expressions that will be substituted for the pac_expectation operators
+      expr_t growth_correction_term = dynamic_model.Zero;
+      if (growth[name])
+        growth_correction_term = dynamic_model.AddTimes(
+            growth[name], dynamic_model.AddVariable(growth_neutrality_params[name]));
+      if (aux_model_name[name].empty())
+        {
+          if (target_info.contains(name))
+            {
+              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(
+                name, symbol_table.getID(discount[name]), pacEquationMaxLag(name),
+                growth_correction_term, auxname[name], diff_subst_table, aux_var_symb_ids,
+                aux_param_symb_ids, pac_expectation_substitution);
+        }
+      else
+        {
+          if (target_info.contains(name))
+            {
+              assert(growth_correction_term == dynamic_model.Zero);
+              dynamic_model.computePacBackwardExpectationSubstitutionWithComponents(
+                  name, lhs[name], max_lag, aux_model_type[name], get<2>(target_info[name]),
+                  pac_expectation_substitution);
+            }
+          else
+            dynamic_model.computePacBackwardExpectationSubstitution(
+                name, lhs[name], max_lag, aux_model_type[name], growth_correction_term,
+                auxname[name], aux_var_symb_ids, aux_param_symb_ids, pac_expectation_substitution);
+        }
+    }
+
+  // Actually perform the substitution of pac_expectation
+  dynamic_model.substitutePacExpectation(pac_expectation_substitution, eq_name);
+  dynamic_model.checkNoRemainingPacExpectation();
+
+  // Check that there is no remaining pac_target_nonstationary operator
+  dynamic_model.checkNoRemainingPacTargetNonstationary();
+}
+
+void
+PacModelTable::writeOutput(ostream& output) const
+{
+  // Helper to print the “growth_info” structure (linear decomposition of growth)
+  auto growth_info_helper = [&](const string& fieldname, const growth_info_t& gi) {
+    for (int i {1}; auto [growth_symb_id, growth_lag, param_id, constant] : gi)
+      {
+        string structname = fieldname + "(" + to_string(i++) + ").";
+        if (growth_symb_id)
+          {
+            string var_field = "endo_id";
+            if (symbol_table.getType(*growth_symb_id) == SymbolType::exogenous)
+              {
+                var_field = "exo_id";
+                output << structname << "endo_id = 0;" << endl;
+              }
+            else
+              output << structname << "exo_id = 0;" << endl;
+            try
+              {
+                // case when this is not the highest lag of the growth variable
+                int aux_symb_id = symbol_table.searchAuxiliaryVars(*growth_symb_id, growth_lag);
+                output << structname << var_field << " = "
+                       << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
+                       << structname << "lag = 0;" << endl;
+              }
+            catch (...)
+              {
+                try
+                  {
+                    // case when this is the highest lag of the growth variable
+                    int tmp_growth_lag = growth_lag + 1;
+                    int aux_symb_id
+                        = symbol_table.searchAuxiliaryVars(*growth_symb_id, tmp_growth_lag);
+                    output << structname << var_field << " = "
+                           << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
+                           << structname << "lag = -1;" << endl;
+                  }
+                catch (...)
+                  {
+                    // case when there is no aux var for the variable
+                    output << structname << var_field << " = "
+                           << symbol_table.getTypeSpecificID(*growth_symb_id) + 1 << ";" << endl
+                           << structname << "lag = " << growth_lag << ";" << endl;
+                  }
+              }
+          }
+        else
+          output << structname << "endo_id = 0;" << endl
+                 << structname << "exo_id = 0;" << endl
+                 << structname << "lag = 0;" << endl;
+        output << structname
+               << "param_id = " << (param_id ? symbol_table.getTypeSpecificID(*param_id) + 1 : 0)
+               << ";" << endl
+               << structname << "constant = " << constant << ";" << endl;
+      }
+  };
+
+  for (const auto& name : names)
+    {
+      output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name.at(name) << "';"
+             << endl
+             << "M_.pac." << name
+             << ".discount_index = " << symbol_table.getTypeSpecificID(discount.at(name)) + 1 << ";"
+             << endl;
+
+      if (growth.at(name))
+        {
+          output << "M_.pac." << name << ".growth_str = '";
+          original_growth.at(name)->writeJsonOutput(output, {}, {}, true);
+          output << "';" << endl;
+          growth_info_helper("M_.pac." + name + ".growth_linear_comb", growth_info.at(name));
+        }
+    }
+
+  // Write the auxiliary parameter IDs created for the pac_expectation operator
+  for (auto& [name, ids] : aux_param_symb_ids)
+    {
+      output << "M_.pac." << name << "."
+             << (aux_model_name.at(name).empty() ? "mce.alpha" : "h_param_indices") << " = [";
+      for (auto id : ids)
+        output << symbol_table.getTypeSpecificID(id) + 1 << " ";
+      output << "];" << endl;
+    }
+
+  // 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.z" : "aux_id")
+           << " = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
+
+  // Write PAC equation name info
+  for (auto& [name, eq] : eq_name)
+    output << "M_.pac." << name << ".eq_name = '" << eq << "';" << endl;
+
+  for (auto& [model, growth_neutrality_param_index] : growth_neutrality_params)
+    output << "M_.pac." << model << ".growth_neutrality_param_index = "
+           << symbol_table.getTypeSpecificID(growth_neutrality_param_index) + 1 << ";" << endl;
+
+  for (auto& [model, type] : aux_model_type)
+    output << "M_.pac." << model << ".auxiliary_model_type = '" << type << "';" << endl;
+
+  for (auto& [name, k] : kind)
+    if (!aux_model_name.empty())
+      output << "M_.pac." << name << ".kind = '"
+             << (k == PacTargetKind::unspecified ? "" : kindToString(k)) << "';" << endl;
+
+  for (auto& [name, val] : equation_info)
+    {
+      auto& [lhs_pac_var, optim_share_index, ar_params_and_vars, ec_params_and_vars,
+             non_optim_vars_params_and_constants, additive_vars_params_and_constants,
+             optim_additive_vars_params_and_constants]
+          = val;
+      output << "M_.pac." << name
+             << ".lhs_var = " << symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";"
+             << endl;
+
+      if (optim_share_index)
+        output << "M_.pac." << name << ".share_of_optimizing_agents_index = "
+               << symbol_table.getTypeSpecificID(*optim_share_index) + 1 << ";" << endl;
+
+      output << "M_.pac." << name
+             << ".ec.params = " << symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1
+             << ";" << endl
+             << "M_.pac." << name << ".ec.vars = [";
+      for (auto& it : ec_params_and_vars.second)
+        output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+      output << "];" << endl << "M_.pac." << name << ".ec.istarget = [";
+      for (auto& it : ec_params_and_vars.second)
+        output << boolalpha << get<1>(it) << " ";
+      output << "];" << endl << "M_.pac." << name << ".ec.scale = [";
+      for (auto& it : ec_params_and_vars.second)
+        output << get<2>(it) << " ";
+      output << "];" << endl << "M_.pac." << name << ".ec.isendo = [";
+      for (auto& it : ec_params_and_vars.second)
+        switch (symbol_table.getType(get<0>(it)))
+          {
+          case SymbolType::endogenous:
+            output << "true ";
+            break;
+          case SymbolType::exogenous:
+            output << "false ";
+            break;
+          default:
+            cerr << "expecting endogenous or exogenous" << endl;
+            exit(EXIT_FAILURE);
+          }
+      output << "];" << endl << "M_.pac." << name << ".ar.params = [";
+      for (auto& [pid, vid, vlag] : ar_params_and_vars)
+        output << (pid ? symbol_table.getTypeSpecificID(*pid) + 1 : -1) << " ";
+      output << "];" << endl << "M_.pac." << name << ".ar.vars = [";
+      for (auto& [pid, vid, vlag] : ar_params_and_vars)
+        output << (vid ? symbol_table.getTypeSpecificID(*vid) + 1 : -1) << " ";
+      output << "];" << endl << "M_.pac." << name << ".ar.lags = [";
+      for (auto& [pid, vid, vlag] : ar_params_and_vars)
+        output << vlag << " ";
+      output << "];" << endl
+             << "M_.pac." << name << ".max_lag = " << pacEquationMaxLag(name) << ";" << endl;
+      if (!non_optim_vars_params_and_constants.empty())
+        {
+          output << "M_.pac." << name << ".non_optimizing_behaviour.params = [";
+          for (auto& it : non_optim_vars_params_and_constants)
+            if (get<2>(it))
+              output << symbol_table.getTypeSpecificID(*get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl << "M_.pac." << name << ".non_optimizing_behaviour.vars = [";
+          for (auto& it : non_optim_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl << "M_.pac." << name << ".non_optimizing_behaviour.isendo = [";
+          for (auto& it : non_optim_vars_params_and_constants)
+            switch (symbol_table.getType(get<0>(it)))
+              {
+              case SymbolType::endogenous:
+                output << "true ";
+                break;
+              case SymbolType::exogenous:
+                output << "false ";
+                break;
+              default:
+                cerr << "expecting endogenous or exogenous" << endl;
+                exit(EXIT_FAILURE);
+              }
+          output << "];" << endl << "M_.pac." << name << ".non_optimizing_behaviour.lags = [";
+          for (auto& it : non_optim_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl
+                 << "M_.pac." << name << ".non_optimizing_behaviour.scaling_factor = [";
+          for (auto& it : non_optim_vars_params_and_constants)
+            output << get<3>(it) << " ";
+          output << "];" << endl;
+        }
+      if (!additive_vars_params_and_constants.empty())
+        {
+          output << "M_.pac." << name << ".additive.params = [";
+          for (auto& it : additive_vars_params_and_constants)
+            if (get<2>(it))
+              output << symbol_table.getTypeSpecificID(*get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl << "M_.pac." << name << ".additive.vars = [";
+          for (auto& it : additive_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl << "M_.pac." << name << ".additive.isendo = [";
+          for (auto& it : additive_vars_params_and_constants)
+            switch (symbol_table.getType(get<0>(it)))
+              {
+              case SymbolType::endogenous:
+                output << "true ";
+                break;
+              case SymbolType::exogenous:
+                output << "false ";
+                break;
+              default:
+                cerr << "expecting endogenous or exogenous" << endl;
+                exit(EXIT_FAILURE);
+              }
+          output << "];" << endl << "M_.pac." << name << ".additive.lags = [";
+          for (auto& it : additive_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl << "M_.pac." << name << ".additive.scaling_factor = [";
+          for (auto& it : additive_vars_params_and_constants)
+            output << get<3>(it) << " ";
+          output << "];" << endl;
+        }
+      if (!optim_additive_vars_params_and_constants.empty())
+        {
+          output << "M_.pac." << name << ".optim_additive.params = [";
+          for (auto& it : optim_additive_vars_params_and_constants)
+            if (get<2>(it))
+              output << symbol_table.getTypeSpecificID(*get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl << "M_.pac." << name << ".optim_additive.vars = [";
+          for (auto& it : optim_additive_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl << "M_.pac." << name << ".optim_additive.isendo = [";
+          for (auto& it : optim_additive_vars_params_and_constants)
+            switch (symbol_table.getType(get<0>(it)))
+              {
+              case SymbolType::endogenous:
+                output << "true ";
+                break;
+              case SymbolType::exogenous:
+                output << "false ";
+                break;
+              default:
+                cerr << "expecting endogenous or exogenous" << endl;
+                exit(EXIT_FAILURE);
+              }
+          output << "];" << endl << "M_.pac." << name << ".optim_additive.lags = [";
+          for (auto& it : optim_additive_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl << "M_.pac." << name << ".optim_additive.scaling_factor = [";
+          for (auto& it : optim_additive_vars_params_and_constants)
+            output << get<3>(it) << " ";
+          output << "];" << endl;
+        }
+    }
+
+  for (auto& [name, val] : target_info)
+    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))
+      {
+        string fieldname = "M_.pac." + name + ".components(" + to_string(component_idx) + ")";
+        output << fieldname << ".aux_id = " << symbol_table.getTypeSpecificID(auxname) + 1 << ";"
+               << endl
+               << fieldname
+               << ".endo_var = " << dynamic_cast<VariableNode*>(component)->getTypeSpecificID() + 1
+               << ";" << endl
+               << fieldname << ".kind = '" << kindToString(kind) << "';" << endl
+               << fieldname << ".h_param_indices = [";
+        for (int id : h_indices)
+          output << symbol_table.getTypeSpecificID(id) + 1 << " ";
+        output << "];" << endl << fieldname << ".coeff_str = '";
+        coeff->writeJsonOutput(output, {}, {}, true);
+        output << "';" << endl;
+        if (growth_component)
+          {
+            output << fieldname << ".growth_neutrality_param_index = "
+                   << symbol_table.getTypeSpecificID(growth_neutrality_param) + 1 << ";" << endl
+                   << fieldname << ".growth_str = '";
+            original_growth_component->writeJsonOutput(output, {}, {}, true);
+            output << "';" << endl;
+            growth_info_helper(fieldname + ".growth_linear_comb", growth_component_info);
+          }
+        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
+PacModelTable::writeJsonOutput(ostream& output) const
+{
+  for (bool printed_something {false}; const auto& name : names)
+    {
+      /* The calling method has already added a comma, so don’t output one for
+         the first statement */
+      if (exchange(printed_something, true))
+        output << ", ";
+      output << R"({"statementName": "pac_model",)"
+             << R"("model_name": ")" << name << R"(",)"
+             << R"("auxiliary_model_name": ")" << aux_model_name.at(name) << R"(",)"
+             << R"("discount_index": )" << symbol_table.getTypeSpecificID(discount.at(name)) + 1;
+      if (growth.at(name))
+        {
+          output << R"(,"growth_str": ")";
+          original_growth.at(name)->writeJsonOutput(output, {}, {}, true);
+          output << R"(")";
+        }
+      output << "}" << endl;
+    }
+
+  for (auto& [name, val] : target_info)
+    {
+      output << R"(, {"statementName": "pac_target_info", "model_name": ")" << name
+             << R"(", "target": ")";
+      get<0>(val)->writeJsonOutput(output, {}, {}, true);
+      output << R"(", "auxname_target_nonstationary": ")" << get<1>(val) << R"(", "components": [)";
+      for (auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                  h_indices, original_growth_component, growth_component_info] : get<2>(val))
+        {
+          if (component != get<0>(get<2>(val).front()))
+            output << ", ";
+          output << R"({"component": ")";
+          component->writeJsonOutput(output, {}, {}, true);
+          output << R"(", "auxname": ")" << auxname << R"(", "kind": ")" << kindToString(kind);
+          if (growth_component)
+            {
+              output << R"(", "growth_str": ")";
+              original_growth_component->writeJsonOutput(output, {}, {}, true);
+            }
+          output << R"("})";
+        }
+      output << "]}" << endl;
+    }
+}
+
+int
+PacModelTable::pacEquationMaxLag(const string& name_arg) const
+{
+  return get<2>(equation_info.at(name_arg)).size();
+}
+
+string
+PacModelTable::kindToString(PacTargetKind kind)
+{
+  switch (kind)
+    {
+    case PacTargetKind::unspecified:
+      cerr << "Internal error: kind should not be unspecified" << endl;
+      exit(EXIT_FAILURE);
+    case PacTargetKind::ll:
+      return "ll";
+    case PacTargetKind::dl:
+      return "dl";
+    case PacTargetKind::dd:
+      return "dd";
+    }
+  __builtin_unreachable(); // Silence GCC warning
+}
+
+void
+PacModelTable::setTargetExpr(const string& name_arg, expr_t target)
+{
+  get<0>(target_info[name_arg]) = target;
+}
+
+void
+PacModelTable::setTargetAuxnameNonstationary(const string& name_arg, string auxname)
+{
+  get<1>(target_info[name_arg]) = move(auxname);
+}
+
+void
+PacModelTable::addTargetComponent(const string& name_arg, target_component_t component)
+{
+  get<7>(component) = get<1>(component); // original_growth = growth
+  get<2>(target_info[name_arg]).emplace_back(move(component));
+}
+
+void
+PacModelTable::writeTargetCoefficientsFile(const string& basename) const
+{
+  if (target_info.empty())
+    return;
+
+  filesystem::path filename {DataTree::packageDir(basename) / "pac_target_coefficients.m"};
+  ofstream output {filename, ios::out | ios::binary};
+  if (!output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  output << "function coeffs = pac_target_coefficients(model_name, params)" << endl;
+  for (auto& [model_name, val] : target_info)
+    {
+      output << "  if strcmp(model_name, '" << model_name << "')" << endl
+             << "    coeffs = NaN(" << get<2>(val).size() << ",1);" << endl;
+      for (int i {1};
+           auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+                  h_indices, original_growth_component, growth_component_info] : get<2>(val))
+        {
+          output << "    coeffs(" << i++ << ") = ";
+          coeff->writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
+          output << ";" << endl;
+        }
+      output << "    return" << endl << "  end" << endl;
+    }
+  output << "  error([ 'Unknown PAC model: ' model_name ])" << endl << "end" << endl;
+  output.close();
+}

src/SubModel.hh

+/*
+ * Copyright © 2018-2023 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.SS
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef SUB_MODEL_HH
+#define SUB_MODEL_HH
+
+#include <iostream>
+#include <map>
+#include <optional>
+#include <set>
+#include <vector>
+
+#include "ExprNode.hh"
+#include "Statement.hh"
+#include "SymbolList.hh"
+#include "SymbolTable.hh"
+
+// DynamicModel.hh can’t be included here, otherwise it would be a circular dependency
+class DynamicModel;
+
+using namespace std;
+
+//! A table with all Trend Component Models in the .mod file
+/*!
+  The unique name of the trend component model is the identifier
+*/
+class TrendComponentModelTable
+{
+private:
+  SymbolTable& symbol_table;
+  set<string> names;
+  map<string, vector<string>> eqtags, target_eqtags;
+  map<string, vector<int>> eqnums, target_eqnums, nontarget_eqnums, max_lags, lhs, target_lhs,
+      nontarget_lhs;
+  map<string, vector<optional<int>>> orig_diff_var;
+  map<string, vector<set<pair<int, int>>>> rhs;
+  map<string, vector<bool>> diff;
+  map<string, vector<expr_t>> lhs_expr_t;
+  map<string, vector<optional<int>>> target_vars;
+  map<string, map<tuple<int, int, int>, expr_t>> AR; // name -> (eqn, lag, lhs_symb_id) -> expr_t
+  /* Note that A0 in the trend-component model context is not the same thing as
+     in the structural VAR context. */
+  map<string, map<tuple<int, int>, expr_t>> A0, A0star; // name -> (eqn, col) -> expr_t
+public:
+  explicit TrendComponentModelTable(SymbolTable& symbol_table_arg);
+
+  //! Add a trend component model
+  void addTrendComponentModel(string name_arg, vector<string> eqtags_arg,
+                              vector<string> target_eqtags_arg);
+
+  [[nodiscard]] inline bool isExistingTrendComponentModelName(const string& name_arg) const;
+  [[nodiscard]] inline bool empty() const;
+
+  [[nodiscard]] const map<string, vector<string>>& getEqTags() const;
+  [[nodiscard]] const vector<string>& getEqTags(const string& name_arg) const;
+  [[nodiscard]] const map<string, vector<string>>& getTargetEqTags() const;
+  [[nodiscard]] const map<string, vector<int>>& getEqNums() const;
+  [[nodiscard]] const map<string, vector<int>>& getTargetEqNums() const;
+  [[nodiscard]] const vector<int>& getTargetEqNums(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getEqNums(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getMaxLags(const string& name_arg) const;
+  [[nodiscard]] int getMaxLag(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getLhs(const string& name_arg) const;
+  [[nodiscard]] const vector<expr_t>& getLhsExprT(const string& name_arg) const;
+  [[nodiscard]] const vector<bool>& getDiff(const string& name_arg) const;
+  [[nodiscard]] const map<string, vector<int>>& getNonTargetEqNums() const;
+  [[nodiscard]] const vector<int>& getNonTargetEqNums(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getNonTargetLhs(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getTargetLhs(const string& name_arg) const;
+
+  void setVals(map<string, vector<int>> eqnums_arg, map<string, vector<int>> target_eqnums_arg,
+               map<string, vector<int>> lhs_arg, map<string, vector<expr_t>> lhs_expr_t_arg);
+  void setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg);
+  void setMaxLags(map<string, vector<int>> max_lags_arg);
+  void setDiff(map<string, vector<bool>> diff_arg);
+  void setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg);
+  void setTargetVar(map<string, vector<optional<int>>> target_vars_arg);
+  void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
+  void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg,
+             map<string, map<tuple<int, int>, expr_t>> A0star_arg);
+
+  //! Write output of this class
+  void writeOutput(const string& basename, ostream& output) const;
+
+  //! Write JSON Output
+  void writeJsonOutput(ostream& output) const;
+
+private:
+  void checkModelName(const string& name_arg) const;
+  void setNonTargetEqnums();
+};
+
+inline bool
+TrendComponentModelTable::isExistingTrendComponentModelName(const string& name_arg) const
+{
+  return names.contains(name_arg);
+}
+
+inline bool
+TrendComponentModelTable::empty() const
+{
+  return names.empty();
+}
+
+class VarModelTable
+{
+private:
+  SymbolTable& symbol_table;
+  set<string> names;
+  map<string, bool> structural; // Whether VARs are structural or reduced-form
+  map<string, vector<string>> eqtags;
+  map<string, vector<int>> eqnums, max_lags, lhs, lhs_orig_symb_ids;
+  map<string, vector<optional<int>>> orig_diff_var;
+  map<string, vector<set<pair<int, int>>>>
+      rhs; // name -> for each equation: set of pairs (var, lag)
+  map<string, vector<bool>> diff;
+  map<string, vector<expr_t>> lhs_expr_t;
+  map<string, map<tuple<int, int, int>, expr_t>>
+      AR; // name -> (eqn, lag, lhs_symb_id) -> param_expr_t
+  /* The A0 matrix is mainly for structural VARs. For reduced-form VARs, it
+     will be equal to the identity matrix. Also note that A0 in the structural
+     VAR context is not the same thing as in the trend-component model
+     context. */
+  map<string, map<tuple<int, int>, expr_t>> A0; // name -> (eqn, lhs_symb_id) -> param_expr_t
+  map<string, map<int, expr_t>> constants;      // name -> eqn -> constant
+public:
+  explicit VarModelTable(SymbolTable& symbol_table_arg);
+
+  //! Add a VAR model
+  void addVarModel(string name, bool structural_arg, vector<string> eqtags);
+
+  [[nodiscard]] inline bool isExistingVarModelName(const string& name_arg) const;
+  [[nodiscard]] inline bool empty() const;
+
+  [[nodiscard]] const map<string, bool>& getStructural() const;
+  [[nodiscard]] const map<string, vector<string>>& getEqTags() const;
+  [[nodiscard]] const vector<string>& getEqTags(const string& name_arg) const;
+  [[nodiscard]] const map<string, vector<int>>& getEqNums() const;
+  [[nodiscard]] const vector<bool>& getDiff(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getEqNums(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getMaxLags(const string& name_arg) const;
+  [[nodiscard]] int getMaxLag(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getLhs(const string& name_arg) const;
+  [[nodiscard]] const vector<int>& getLhsOrigIds(const string& name_arg) const;
+  [[nodiscard]] const vector<set<pair<int, int>>>& getRhs(const string& name_arg) const;
+  [[nodiscard]] const vector<expr_t>& getLhsExprT(const string& name_arg) const;
+
+  void setEqNums(map<string, vector<int>> eqnums_arg);
+  void setLhs(map<string, vector<int>> lhs_arg);
+  void setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg);
+  void setLhsExprT(map<string, vector<expr_t>> lhs_expr_t_arg);
+  void setDiff(map<string, vector<bool>> diff_arg);
+  void setMaxLags(map<string, vector<int>> max_lags_arg);
+  void setOrigDiffVar(map<string, vector<optional<int>>> orig_diff_var_arg);
+  void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
+  void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg);
+  void setConstants(map<string, map<int, expr_t>> constants_arg);
+
+  //! Write output of this class
+  void writeOutput(const string& basename, ostream& output) const;
+
+  //! Write JSON Output
+  void writeJsonOutput(ostream& output) const;
+
+private:
+  void checkModelName(const string& name_arg) const;
+};
+
+inline bool
+VarModelTable::isExistingVarModelName(const string& name_arg) const
+{
+  return names.contains(name_arg);
+}
+
+inline bool
+VarModelTable::empty() const
+{
+  return names.empty();
+}
+
+class VarExpectationModelTable
+{
+private:
+  SymbolTable& symbol_table;
+  set<string> names;
+  map<string, expr_t> expression;
+  map<string, string> aux_model_name;
+  map<string, string> horizon;
+  map<string, expr_t> discount;
+  map<string, int> time_shift;
+  // For each model, list of generated auxiliary param ids, in variable-major order
+  map<string, vector<int>> aux_param_symb_ids;
+  // Decomposition of the expression
+  map<string, vector<tuple<int, optional<int>, double>>> vars_params_constants;
+
+public:
+  explicit VarExpectationModelTable(SymbolTable& symbol_table_arg);
+  void addVarExpectationModel(string name_arg, expr_t expression_arg, string aux_model_name_arg,
+                              string horizon_arg, expr_t discount_arg, int time_shift_arg);
+  [[nodiscard]] bool isExistingVarExpectationModelName(const string& name_arg) const;
+  [[nodiscard]] bool empty() const;
+  void substituteUnaryOpsInExpression(const lag_equivalence_table_t& nodes,
+                                      ExprNode::subst_table_t& subst_table,
+                                      vector<BinaryOpNode*>& neweqs);
+  // Called by DynamicModel::substituteDiff()
+  void substituteDiffNodesInExpression(const lag_equivalence_table_t& diff_nodes,
+                                       ExprNode::subst_table_t& diff_subst_table,
+                                       vector<BinaryOpNode*>& neweqs);
+  void transformPass(ExprNode::subst_table_t& diff_subst_table, DynamicModel& dynamic_model,
+                     const VarModelTable& var_model_table,
+                     const TrendComponentModelTable& trend_component_model_table);
+  void writeOutput(ostream& output) const;
+  void writeJsonOutput(ostream& output) const;
+};
+
+class PacModelTable
+{
+private:
+  SymbolTable& symbol_table;
+  set<string> names;
+  map<string, string> aux_model_name;
+  map<string, string> discount;
+  /* The growth expressions belong to the main dynamic_model from the ModFile
+     instance. The growth expression is necessarily nullptr for a model with a
+     pac_target_info block. */
+  map<string, expr_t> growth, original_growth;
+  /* Information about the structure of growth expressions (which must be a
+     linear combination of variables, possibly with additional constants).
+     Each tuple represents a term: (endo_id, lag, param_id, constant) */
+  using growth_info_t = vector<tuple<optional<int>, int, optional<int>, double>>;
+  map<string, growth_info_t> growth_info;
+  // The “auxname” option of pac_model (empty if not passed)
+  map<string, string> auxname;
+  // The “kind” option of pac_model (“undefined” if not passed)
+  map<string, PacTargetKind> kind;
+
+  /* Stores the name of the PAC equation associated to the model.
+     pac_model_name → eq_name */
+  map<string, string> eq_name;
+
+  /* Stores symb_ids for auxiliary endogenous created for the expression
+     substituted to the pac_expectation operator:
+     - in the backward case, this auxiliary contains exactly the
+     pac_expectation value
+     - in the MCE case, this auxiliary represents Z₁ (i.e. without the growth
+     correction term)
+     Note that this structure is not used in the presence of the
+     pac_target_info block.
+     pac_model_name → symb_id */
+  map<string, int> aux_var_symb_ids;
+  /* Stores symb_ids for auxiliary parameters created for the expression
+     substituted to the pac_expectation operator (excluding the growth
+     neutrality correction):
+     - in the backward case, contains the “h” parameters
+     - in the MCE case, contains the “α” parameters
+     Note that this structure is not used in the presence of the
+     pac_target_info block.
+     pac_model_name → symb_ids */
+  map<string, vector<int>> aux_param_symb_ids;
+  /* Stores indices for growth neutrality parameters
+     pac_model_name → growth_neutrality_param_index.
+     This map is not used for PAC models with a pac_target_info block. */
+  map<string, int> growth_neutrality_params;
+
+  // Stores LHS vars (only for backward PAC models)
+  map<string, vector<int>> lhs;
+
+  // Stores auxiliary model type (only for backward PAC models)
+  map<string, string> aux_model_type;
+
+public:
+  /* Stores info about PAC equations
+     pac_model_name →
+         (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars,
+     non_optim_vars_params_and_constants, additive_vars_params_and_constants,
+     optim_additive_vars_params_and_constants)
+  */
+  using equation_info_t
+      = map<string,
+            tuple<pair<int, int>, optional<int>, vector<tuple<optional<int>, optional<int>, int>>,
+                  pair<int, vector<tuple<int, bool, int>>>,
+                  vector<tuple<int, int, optional<int>, double>>,
+                  vector<tuple<int, int, optional<int>, double>>,
+                  vector<tuple<int, int, optional<int>, double>>>>;
+
+private:
+  equation_info_t equation_info;
+
+public:
+  /* (component variable/expr, growth, auxname, kind, coeff. in the linear
+     combination, growth_param ID (unused if growth is nullptr), vector of h parameters,
+     original_growth, growth_info) */
+  using target_component_t = tuple<expr_t, expr_t, string, PacTargetKind, expr_t, int, vector<int>,
+                                   expr_t, growth_info_t>;
+
+private:
+  // pac_model_name → (target variable/expr, auxname_target_nonstationary, target components)
+  map<string, tuple<expr_t, string, vector<target_component_t>>> target_info;
+
+  [[nodiscard]] int pacEquationMaxLag(const string& name_arg) const;
+
+  // Return a text representation of a kind (but fails on “unspecified” kind value)
+  static string kindToString(PacTargetKind kind);
+
+public:
+  explicit PacModelTable(SymbolTable& symbol_table_arg);
+  void addPacModel(string name_arg, string aux_model_name_arg, string discount_arg,
+                   expr_t growth_arg, string auxname_arg, PacTargetKind kind_arg);
+  [[nodiscard]] bool isExistingPacModelName(const string& name_arg) const;
+  [[nodiscard]] bool empty() const;
+  void checkPass(ModFileStructure& mod_file_struct);
+  // Called by DynamicModel::substituteUnaryOps()
+  void substituteUnaryOpsInGrowth(const lag_equivalence_table_t& nodes,
+                                  ExprNode::subst_table_t& subst_table,
+                                  vector<BinaryOpNode*>& neweqs);
+  void findDiffNodesInGrowth(lag_equivalence_table_t& diff_nodes) const;
+  // Called by DynamicModel::substituteDiff()
+  void substituteDiffNodesInGrowth(const lag_equivalence_table_t& diff_nodes,
+                                   ExprNode::subst_table_t& diff_subst_table,
+                                   vector<BinaryOpNode*>& neweqs);
+  // Must be called after substituteDiffNodesInGrowth() and substituteUnaryOpsInGrowth()
+  void transformPass(const lag_equivalence_table_t& unary_ops_nodes,
+                     ExprNode::subst_table_t& unary_ops_subst_table,
+                     const lag_equivalence_table_t& diff_nodes,
+                     ExprNode::subst_table_t& diff_subst_table, DynamicModel& dynamic_model,
+                     const VarModelTable& var_model_table,
+                     const TrendComponentModelTable& trend_component_model_table);
+  void writeOutput(ostream& output) const;
+  void writeJsonOutput(ostream& output) const;
+  void setTargetExpr(const string& name_arg, expr_t target);
+  void setTargetAuxnameNonstationary(const string& name_arg, string auxname);
+  /* Only the first four elements of the tuple are expected to be set by the
+     caller. The other ones will be filled by this class. */
+  void addTargetComponent(const string& name_arg, target_component_t component);
+  void writeTargetCoefficientsFile(const string& basename) const;
+};
+
+#endif

src/SymbolList.cc

 /*
- * Copyright (C) 2003-2018 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -14,27 +14,120 @@
  * 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 <http://www.gnu.org/licenses/>.
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
+#include <regex>
+
 #include "SymbolList.hh"
 
+SymbolList::SymbolList(vector<string> symbols_arg) : symbols {move(symbols_arg)}
+{
+}
+
 void
-SymbolList::addSymbol(const string &symbol)
+SymbolList::checkPass(WarningConsolidation& warnings, const vector<SymbolType>& types,
+                      const SymbolTable& symbol_table) const noexcept(false)
+{
+  if (types.empty())
+    return;
+
+  smatch m;
+  string regex_str = "AUX_EXPECT_|MULT_";
+  for (auto type : types)
+    if (type == SymbolType::endogenous)
       {
-  symbols.push_back(symbol);
+        regex_str += "|AUX_ENDO_|LOG_";
+        break;
+      }
+  regex re("^(" + regex_str + ")");
+  for (const auto& symbol : symbols)
+    {
+      if (!symbol_table.exists(symbol))
+        {
+          if (regex_search(symbol, m, re))
+            {
+              warnings
+                  << "WARNING: symbol_list variable " << symbol << " has not yet been declared. "
+                  << "This is being ignored because the variable name corresponds to a possible "
+                  << "auxiliary variable name." << endl;
+              return;
+            }
+          else
+            throw SymbolListException {"Variable " + symbol + " was not declared."};
+        }
+
+      if (ranges::none_of(types,
+                          [&](SymbolType type) { return symbol_table.getType(symbol) == type; }))
+        {
+          string valid_types;
+          for (auto type : types)
+            switch (type)
+              {
+              case SymbolType::endogenous:
+                valid_types += "endogenous, ";
+                break;
+              case SymbolType::exogenous:
+                valid_types += "exogenous, ";
+                break;
+              case SymbolType::epilogue:
+                valid_types += "epilogue, ";
+                break;
+              case SymbolType::parameter:
+                valid_types += "parameter, ";
+                break;
+              case SymbolType::exogenousDet:
+                valid_types += "exogenousDet, ";
+                break;
+              case SymbolType::trend:
+                valid_types += "trend, ";
+                break;
+              case SymbolType::logTrend:
+                valid_types += "logTrend, ";
+                break;
+              case SymbolType::modFileLocalVariable:
+                valid_types += "modFileLocalVariable, ";
+                break;
+              case SymbolType::modelLocalVariable:
+                valid_types += "modelLocalVariable, ";
+                break;
+              case SymbolType::externalFunction:
+                valid_types += "externalFunction, ";
+                break;
+              case SymbolType::statementDeclaredVariable:
+                valid_types += "statementDeclaredVariable, ";
+                break;
+              case SymbolType::unusedEndogenous:
+                valid_types += "unusedEndogenous, ";
+                break;
+              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 + "}"};
+        }
+    }
 }
 
 void
 SymbolList::writeOutput(const string& varname, ostream& output) const
 {
   output << varname << " = {";
-  for (vector<string>::const_iterator it = symbols.begin();
-       it != symbols.end(); ++it)
+  for (bool printed_something {false}; const auto& name : symbols)
     {
-      if (it != symbols.begin())
+      if (exchange(printed_something, true))
         output << ";";
-      output << "'" << *it << "'";
+      output << "'" << name << "'";
     }
   output << "};" << endl;
 }
@@ -42,19 +135,32 @@ SymbolList::writeOutput(const string &varname, ostream &output) const
 void
 SymbolList::writeJsonOutput(ostream& output) const
 {
-  output << "\"symbol_list\": [";
-  for (vector<string>::const_iterator it = symbols.begin();
-       it != symbols.end(); ++it)
+  output << R"("symbol_list": [)";
+  for (bool printed_something {false}; const auto& name : symbols)
     {
-      if (it != symbols.begin())
+      if (exchange(printed_something, true))
         output << ",";
-      output << "\"" << *it << "\"";
+      output << R"(")" << name << R"(")";
     }
   output << "]";
 }
 
+vector<string>
+SymbolList::getSymbols() const
+{
+  return symbols;
+}
+
 void
-SymbolList::clear()
+SymbolList::removeDuplicates(const string& dynare_command, WarningConsolidation& warnings)
 {
-  symbols.clear();
+  vector<string> unique_symbols;
+  for (const auto& it : symbols)
+    if (find(unique_symbols.begin(), unique_symbols.end(), it) == unique_symbols.end())
+      unique_symbols.push_back(it);
+    else
+      warnings << "WARNING: In " << dynare_command << ": " << it
+               << " found more than once in symbol list. Removing all but first occurrence."
+               << endl;
+  symbols = unique_symbols;
 }

src/SymbolList.hh

 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright © 2003-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -14,15 +14,19 @@
  * 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 <http://www.gnu.org/licenses/>.
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
-#ifndef _SYMBOL_LIST_HH
-#define _SYMBOL_LIST_HH
+#ifndef SYMBOL_LIST_HH
+#define SYMBOL_LIST_HH
 
+#include <algorithm>
+#include <ostream>
 #include <string>
 #include <vector>
-#include <ostream>
+
+#include "SymbolTable.hh"
+#include "WarningConsolidation.hh"
 
 using namespace std;
 
@@ -31,30 +35,37 @@ using namespace std;
 class SymbolList
 {
 private:
-  //! Internal container for symbol list
   vector<string> symbols;
+
 public:
-  //! Adds a symbol to the list
-  void addSymbol(const string &symbol);
+  SymbolList() = default;
+  // This constructor is deliberately not marked explicit, to allow implicit conversion
+  SymbolList(vector<string> symbols_arg);
+
+  struct SymbolListException
+  {
+    const string message;
+  };
+  //! Remove duplicate symbols
+  void removeDuplicates(const string& dynare_command, WarningConsolidation& warnings);
+  //! Check symbols to ensure variables have been declared and are endogenous
+  void checkPass(WarningConsolidation& warnings, const vector<SymbolType>& types,
+                 const SymbolTable& symbol_table) const noexcept(false);
   //! Output content in Matlab format
   /*! Creates a string array for Matlab, stored in variable "varname" */
   void writeOutput(const string& varname, ostream& output) const;
+  //! Output content in Matlab format without preceding varname of writeOutput
+  void write(ostream& output) const;
   //! Write JSON output
   void writeJsonOutput(ostream& output) const;
-  //! Clears all content
-  void clear();
-  //! Get a copy of the string vector
-  vector<string>
-  get_symbols() const
-  {
-    return symbols;
-  };
   //! Is Empty
-  int
+  [[nodiscard]] bool
   empty() const
   {
     return symbols.empty();
-  };
+  }
+  //! Return the list of symbols
+  [[nodiscard]] vector<string> getSymbols() const;
 };
 
 #endif

src/SymbolTable.cc

 /*
- * Copyright (C) 2003-2018 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -14,36 +14,25 @@
  * 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 <http://www.gnu.org/licenses/>.
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
 #include <algorithm>
-#include <sstream>
-#include <iostream>
 #include <cassert>
+#include <iostream>
+#include <sstream>
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wold-style-cast"
 #include <boost/algorithm/string/replace.hpp>
+#pragma GCC diagnostic pop
+#include <utility>
 
 #include "SymbolTable.hh"
 
-AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg,
-                       int equation_number_for_multiplier_arg, int information_set_arg,
-                       expr_t expr_node_arg) :
-  symb_id(symb_id_arg),
-  type(type_arg),
-  orig_symb_id(orig_symb_id_arg),
-  orig_lead_lag(orig_lead_lag_arg),
-  equation_number_for_multiplier(equation_number_for_multiplier_arg),
-  information_set(information_set_arg),
-  expr_node(expr_node_arg)
-{
-}
-
-SymbolTable::SymbolTable() : frozen(false)
-{
-}
-
 int
-SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_name, const vector<pair<string *, string *> *> *partition_value) throw (AlreadyDeclaredException, FrozenException)
+SymbolTable::addSymbol(const string& name, SymbolType type, const string& tex_name,
+                       const vector<pair<string, string>>& partition_value,
+                       const optional<int>& heterogeneity_dimension) noexcept(false)
 {
   if (frozen)
     throw FrozenException();
@@ -51,9 +40,9 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
   if (exists(name))
     {
       if (type_table[getID(name)] == type)
-        throw AlreadyDeclaredException(name, true);
+        throw AlreadyDeclaredException {name, true};
       else
-        throw AlreadyDeclaredException(name, false);
+        throw AlreadyDeclaredException {name, false};
     }
 
   string final_tex_name = tex_name;
@@ -63,18 +52,16 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
       size_t pos = 0;
       while ((pos = final_tex_name.find('_', pos)) != string::npos)
         {
-          final_tex_name.insert(pos, "\\");
+          final_tex_name.insert(pos, R"(\)");
           pos += 2;
         }
     }
 
   string final_long_name = name;
   bool non_long_name_partition_exists = false;
-  if (partition_value)
-    for (vector<pair<string *, string *> *>::const_iterator it = partition_value->begin();
-         it != partition_value->end(); it++)
-      if (*((*it)->first) == "long_name")
-        final_long_name = *((*it)->second);
+  for (const auto& it : partition_value)
+    if (it.first == "long_name")
+      final_long_name = it.second;
     else
       non_long_name_partition_exists = true;
 
@@ -88,54 +75,75 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
   if (non_long_name_partition_exists)
     {
       map<string, string> pmv;
-      for (vector<pair<string *, string *> *>::const_iterator it = partition_value->begin();
-           it != partition_value->end(); it++)
-        pmv[*((*it)->first)] = *((*it)->second);
+      for (const auto& it : partition_value)
+        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) throw (AlreadyDeclaredException, FrozenException)
+SymbolTable::addSymbol(const string& name, SymbolType type) noexcept(false)
 {
-  return addSymbol(name, type, "", NULL);
+  return addSymbol(name, type, "", {}, {});
 }
 
 void
-SymbolTable::freeze() throw (FrozenException)
+SymbolTable::freeze() noexcept(false)
 {
   if (frozen)
     throw FrozenException();
 
   frozen = true;
 
-  for (int i = 0; i < (int) symbol_table.size(); i++)
+  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;
       switch (getType(i))
         {
-        case eEndogenous:
+        case SymbolType::endogenous:
           tsi = endo_ids.size();
           endo_ids.push_back(i);
           break;
-        case eExogenous:
+        case SymbolType::exogenous:
           tsi = exo_ids.size();
           exo_ids.push_back(i);
           break;
-        case eExogenousDet:
+        case SymbolType::exogenousDet:
           tsi = exo_det_ids.size();
           exo_det_ids.push_back(i);
           break;
-        case eParameter:
+        case SymbolType::parameter:
           tsi = param_ids.size();
           param_ids.push_back(i);
           break;
-        default:
-          tsi = -1;
+        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;
         }
-      type_specific_ids.push_back(tsi);
+      type_specific_ids[i] = tsi;
     }
 }
 
@@ -147,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) throw (UnknownSymbolIDException, FrozenException)
+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();
 
@@ -162,57 +176,70 @@ SymbolTable::changeType(int id, SymbolType newtype) throw (UnknownSymbolIDExcept
 }
 
 int
-SymbolTable::getID(SymbolType type, int tsid) const throw (UnknownTypeSpecificIDException, NotYetFrozenException)
+SymbolTable::getID(SymbolType type, int tsid, const optional<int>& heterogeneity_dimension) const
+    noexcept(false)
 {
   if (!frozen)
     throw NotYetFrozenException();
 
   switch (type)
     {
-    case eEndogenous:
-      if (tsid < 0 || tsid >= (int) endo_ids.size())
-        throw UnknownTypeSpecificIDException(tsid, type);
+    case SymbolType::endogenous:
+      if (tsid < 0 || tsid >= static_cast<int>(endo_ids.size()))
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return endo_ids[tsid];
-    case eExogenous:
-      if (tsid < 0 || tsid >= (int) exo_ids.size())
-        throw UnknownTypeSpecificIDException(tsid, type);
+    case SymbolType::exogenous:
+      if (tsid < 0 || tsid >= static_cast<int>(exo_ids.size()))
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return exo_ids[tsid];
-    case eExogenousDet:
-      if (tsid < 0 || tsid >= (int) exo_det_ids.size())
-        throw UnknownTypeSpecificIDException(tsid, type);
+    case SymbolType::exogenousDet:
+      if (tsid < 0 || tsid >= static_cast<int>(exo_det_ids.size()))
+        throw UnknownTypeSpecificIDException {tsid, type, {}};
       else
         return exo_det_ids[tsid];
-    case eParameter:
-      if (tsid < 0 || tsid >= (int) param_ids.size())
-        throw UnknownTypeSpecificIDException(tsid, type);
+    case SymbolType::parameter:
+      if (tsid < 0 || tsid >= static_cast<int>(param_ids.size()))
+        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, {}};
     }
 }
 
 map<string, map<int, string>>
-SymbolTable::getPartitionsForType(enum SymbolType st) const throw (UnknownSymbolIDException)
+SymbolTable::getPartitionsForType(SymbolType st) const noexcept(false)
 {
   map<string, map<int, string>> partitions;
-  for (map<int, map<string, string> >::const_iterator it = partition_value_map.begin();
-       it != partition_value_map.end(); it++)
-    if (getType(it->first) == st)
-      for (map<string, string>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        {
-          if (partitions.find(it1->first) == partitions.end())
-            partitions[it1->first] = map<int, string> ();
-          partitions[it1->first][it->first] = it1->second;
-        }
+  for (const auto& it : partition_value_map)
+    if (getType(it.first) == st)
+      for (const auto& it1 : it.second)
+        partitions[it1.first][it.first] = it1.second;
   return partitions;
 }
 
 void
-SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
+SymbolTable::writeOutput(ostream& output) const noexcept(false)
 {
   if (!frozen)
     throw NotYetFrozenException();
@@ -224,26 +251,34 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       output << "M_.exo_names_long = cell(" << exo_nbr() << ",1);" << endl;
       for (int id = 0; id < exo_nbr(); id++)
         output << "M_.exo_names(" << id + 1 << ") = {'" << getName(exo_ids[id]) << "'};" << endl
-               << "M_.exo_names_tex(" << id+1 << ") = {'" << getTeXName(exo_ids[id]) << "'};" << endl
-               << "M_.exo_names_long(" << id+1 << ") = {'" << getLongName(exo_ids[id]) << "'};" << endl;
-      map<string, map<int, string> > partitions = getPartitionsForType(eExogenous);
-      for (map<string, map<int, string> >::const_iterator it = partitions.begin();
-           it != partitions.end(); it++)
-        if (it->first != "long_name")
-          {
-            map<int, string>::const_iterator it1;
-            output << "M_.exo_partitions." << it->first << " = { ";
+               << "M_.exo_names_tex(" << id + 1 << ") = {'" << getTeXName(exo_ids[id]) << "'};"
+               << endl
+               << "M_.exo_names_long(" << id + 1 << ") = {'" << getLongName(exo_ids[id]) << "'};"
+               << endl;
+      for (auto& partition : getPartitionsForType(SymbolType::exogenous))
+        if (partition.first != "long_name")
+          {
+            output << "M_.exo_partitions." << partition.first << " = { ";
             for (int id = 0; id < exo_nbr(); id++)
               {
                 output << "'";
-                it1 = it->second.find(exo_ids[id]);
-                if (it1 != it->second.end())
+                if (auto it1 = partition.second.find(exo_ids[id]); it1 != partition.second.end())
                   output << it1->second;
                 output << "' ";
               }
             output << "};" << endl;
+            if (partition.first == "status")
+              output << "M_ = set_observed_exogenous_variables(M_);" << endl;
+            if (partition.first == "used")
+              output << "M_ = set_exogenous_variables_for_simulation(M_);" << endl;
           }
     }
+  else
+    {
+      output << "M_.exo_names = {};" << endl;
+      output << "M_.exo_names_tex = {};" << endl;
+      output << "M_.exo_names_long = {};" << endl;
+    }
 
   if (exo_det_nbr() > 0)
     {
@@ -251,22 +286,22 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       output << "M_.exo_det_names_tex = cell(" << exo_det_nbr() << ",1);" << endl;
       output << "M_.exo_det_names_long = cell(" << exo_det_nbr() << ",1);" << endl;
       for (int id = 0; id < exo_det_nbr(); id++)
-        output << "M_.exo_det_names(" << id+1 << ") = {'" << getName(exo_det_ids[id]) << "'};" << endl
-               << "M_.exo_det_names_tex(" << id+1 << ") = {'" << getTeXName(exo_det_ids[id]) << "'};" << endl
-               << "M_.exo_det_names_long(" << id+1 << ") = {'" << getLongName(exo_det_ids[id]) << "'};" << endl;
+        output << "M_.exo_det_names(" << id + 1 << ") = {'" << getName(exo_det_ids[id]) << "'};"
+               << endl
+               << "M_.exo_det_names_tex(" << id + 1 << ") = {'" << getTeXName(exo_det_ids[id])
+               << "'};" << endl
+               << "M_.exo_det_names_long(" << id + 1 << ") = {'" << getLongName(exo_det_ids[id])
+               << "'};" << endl;
       output << "M_.exo_det_partitions = struct();" << endl;
-      map<string, map<int, string> > partitions = getPartitionsForType(eExogenousDet);
-      for (map<string, map<int, string> >::const_iterator it = partitions.begin();
-           it != partitions.end(); it++)
-        if (it->first != "long_name")
+      for (auto& partition : getPartitionsForType(SymbolType::exogenousDet))
+        if (partition.first != "long_name")
           {
-            map<int, string>::const_iterator it1;
-            output << "M_.exo_det_partitions." << it->first << " = { ";
+            output << "M_.exo_det_partitions." << partition.first << " = { ";
             for (int id = 0; id < exo_det_nbr(); id++)
               {
                 output << "'";
-                it1 = it->second.find(exo_det_ids[id]);
-                if (it1 != it->second.end())
+                if (auto it1 = partition.second.find(exo_det_ids[id]);
+                    it1 != partition.second.end())
                   output << it1->second;
                 output << "' ";
               }
@@ -281,21 +316,19 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       output << "M_.endo_names_long = cell(" << endo_nbr() << ",1);" << endl;
       for (int id = 0; id < endo_nbr(); id++)
         output << "M_.endo_names(" << id + 1 << ") = {'" << getName(endo_ids[id]) << "'};" << endl
-               << "M_.endo_names_tex(" << id+1 << ") = {'" << getTeXName(endo_ids[id]) << "'};" << endl
-               << "M_.endo_names_long(" << id+1 << ") = {'" << getLongName(endo_ids[id]) << "'};" << endl;
+               << "M_.endo_names_tex(" << id + 1 << ") = {'" << getTeXName(endo_ids[id]) << "'};"
+               << endl
+               << "M_.endo_names_long(" << id + 1 << ") = {'" << getLongName(endo_ids[id]) << "'};"
+               << endl;
       output << "M_.endo_partitions = struct();" << endl;
-      map<string, map<int, string> > partitions = getPartitionsForType(eEndogenous);
-      for (map<string, map<int, string> >::const_iterator it = partitions.begin();
-           it != partitions.end(); it++)
-        if (it->first != "long_name")
+      for (auto& partition : getPartitionsForType(SymbolType::endogenous))
+        if (partition.first != "long_name")
           {
-            map<int, string>::const_iterator it1;
-            output << "M_.endo_partitions." << it->first << " = { ";
+            output << "M_.endo_partitions." << partition.first << " = { ";
             for (int id = 0; id < endo_nbr(); id++)
               {
                 output << "'";
-                it1 = it->second.find(endo_ids[id]);
-                if (it1 != it->second.end())
+                if (auto it1 = partition.second.find(endo_ids[id]); it1 != partition.second.end())
                   output << it1->second;
                 output << "' ";
               }
@@ -310,31 +343,36 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       output << "M_.param_names_long = cell(" << param_nbr() << ",1);" << endl;
       for (int id = 0; id < param_nbr(); id++)
         {
-          output << "M_.param_names(" << id+1 << ") = {'" << getName(param_ids[id]) << "'};" << endl
-                 << "M_.param_names_tex(" << id+1 << ") = {'" << getTeXName(param_ids[id]) << "'};" << endl
-                 << "M_.param_names_long(" << id+1 << ") = {'" << getLongName(param_ids[id]) << "'};" << endl;
+          output << "M_.param_names(" << id + 1 << ") = {'" << getName(param_ids[id]) << "'};"
+                 << endl
+                 << "M_.param_names_tex(" << id + 1 << ") = {'" << getTeXName(param_ids[id])
+                 << "'};" << endl
+                 << "M_.param_names_long(" << id + 1 << ") = {'" << getLongName(param_ids[id])
+                 << "'};" << endl;
           if (getName(param_ids[id]) == "dsge_prior_weight")
             output << "options_.dsge_var = 1;" << endl;
         }
       output << "M_.param_partitions = struct();" << endl;
-      map<string, map<int, string> > partitions = getPartitionsForType(eParameter);
-      for (map<string, map<int, string> >::const_iterator it = partitions.begin();
-           it != partitions.end(); it++)
-        if (it->first != "long_name")
+      for (auto& partition : getPartitionsForType(SymbolType::parameter))
+        if (partition.first != "long_name")
           {
-            map<int, string>::const_iterator it1;
-            output << "M_.param_partitions." << it->first << " = { ";
+            output << "M_.param_partitions." << partition.first << " = { ";
             for (int id = 0; id < param_nbr(); id++)
               {
                 output << "'";
-                it1 = it->second.find(param_ids[id]);
-                if (it1 != it->second.end())
+                if (auto it1 = partition.second.find(param_ids[id]); it1 != partition.second.end())
                   output << it1->second;
                 output << "' ";
               }
             output << "};" << endl;
           }
     }
+  else
+    {
+      output << "M_.param_names = {};" << endl;
+      output << "M_.param_names_tex = {};" << endl;
+      output << "M_.param_names_long = {};" << endl;
+    }
 
   output << "M_.exo_det_nbr = " << exo_det_nbr() << ";" << endl
          << "M_.exo_nbr = " << exo_nbr() << ";" << endl
@@ -346,445 +384,622 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
   if (aux_vars.size() == 0)
     output << "M_.aux_vars = [];" << endl;
   else
-    for (int i = 0; i < (int) aux_vars.size(); i++)
-      {
-        output << "M_.aux_vars(" << i+1 << ").endo_index = " << getTypeSpecificID(aux_vars[i].get_symb_id())+1 << ";" << endl
-               << "M_.aux_vars(" << i+1 << ").type = " << aux_vars[i].get_type() << ";" << endl;
-        switch (aux_vars[i].get_type())
+    for (int i = 0; i < static_cast<int>(aux_vars.size()); i++)
       {
-          case avEndoLead:
-          case avExoLead:
-            break;
-          case avEndoLag:
-          case avExoLag:
-            output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl
-                   << "M_.aux_vars(" << i+1 << ").orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
+        output << "M_.aux_vars(" << i + 1
+               << ").endo_index = " << getTypeSpecificID(aux_vars[i].symb_id) + 1 << ";" << endl
+               << "M_.aux_vars(" << i + 1 << ").type = " << aux_vars[i].get_type_id() << ";"
+               << endl;
+        switch (aux_vars[i].type)
+          {
+          case AuxVarType::endoLead:
+          case AuxVarType::exoLead:
+          case AuxVarType::expectation:
+          case AuxVarType::pacExpectation:
+          case AuxVarType::pacTargetNonstationary:
+          case AuxVarType::aggregationOp:
+          case AuxVarType::heterogeneousMultiplier:
             break;
-          case avMultiplier:
-            output << "M_.aux_vars(" << i+1 << ").eq_nbr = " << aux_vars[i].get_equation_number_for_multiplier() + 1 << ";" << endl;
+          case AuxVarType::endoLag:
+          case AuxVarType::exoLag:
+          case AuxVarType::logTransform:
+          case AuxVarType::diffLag:
+          case AuxVarType::diffLead:
+          case AuxVarType::diffForward:
+            output << "M_.aux_vars(" << i + 1
+                   << ").orig_index = " << getTypeSpecificID(aux_vars[i].orig_symb_id.value()) + 1
+                   << ";" << endl
+                   << "M_.aux_vars(" << i + 1
+                   << ").orig_lead_lag = " << aux_vars[i].orig_lead_lag.value() << ";" << endl;
             break;
-          case avDiffForward:
-            output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl;
+          case AuxVarType::unaryOp:
+            output << "M_.aux_vars(" << i + 1 << ").unary_op = '" << aux_vars[i].unary_op << "';"
+                   << endl;
+            [[fallthrough]];
+          case AuxVarType::diff:
+            if (aux_vars[i].orig_symb_id)
+              output << "M_.aux_vars(" << i + 1
+                     << ").orig_index = " << getTypeSpecificID(*aux_vars[i].orig_symb_id) + 1 << ";"
+                     << endl
+                     << "M_.aux_vars(" << i + 1
+                     << ").orig_lead_lag = " << aux_vars[i].orig_lead_lag.value() << ";" << endl;
             break;
-          case avExpectation:
-            output << "M_.aux_vars(" << i+1 << ").orig_expr = '\\mathbb{E}_{t"
-                   << (aux_vars[i].get_information_set() < 0 ? "" : "+")
-                   << aux_vars[i].get_information_set() << "}(";
-            aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
-            output << ")';" << endl;
+          case AuxVarType::multiplier:
+            output << "M_.aux_vars(" << i + 1
+                   << ").eq_nbr = " << aux_vars[i].equation_number_for_multiplier + 1 << ";"
+                   << endl;
             break;
           }
+
+        if (expr_t orig_expr = aux_vars[i].expr_node; orig_expr)
+          {
+            output << "M_.aux_vars(" << i + 1 << ").orig_expr = '";
+            orig_expr->writeJsonOutput(output, {}, {});
+            output << "';" << endl;
+          }
       }
 
   if (predeterminedNbr() > 0)
     {
       output << "M_.predetermined_variables = [ ";
-      for (set<int>::const_iterator it = predetermined_variables.begin();
-           it != predetermined_variables.end(); it++)
-        output << getTypeSpecificID(*it)+1 << " ";
+      for (int predetermined_variable : predetermined_variables)
+        output << getTypeSpecificID(predetermined_variable) + 1 << " ";
       output << "];" << endl;
     }
 
   if (observedVariablesNbr() > 0)
     {
-      int ic = 1;
       output << "options_.varobs = cell(" << observedVariablesNbr() << ", 1);" << endl;
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++, ic++)
-        output << "options_.varobs(" << ic << ")  = {'" << getName(*it) << "'};" << endl;
+      for (int ic {1}; int it : varobs)
+        output << "options_.varobs(" << ic++ << ")  = {'" << getName(it) << "'};" << endl;
 
       output << "options_.varobs_id = [ ";
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++)
-        output << getTypeSpecificID(*it)+1 << " ";
+      for (int varob : varobs)
+        output << getTypeSpecificID(varob) + 1 << " ";
       output << " ];" << endl;
     }
-}
 
-void
-SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
+  if (observedExogenousVariablesNbr() > 0)
     {
-  if (!frozen)
-    throw NotYetFrozenException();
+      output << "options_.varexobs = cell(1);" << endl;
+      for (int ic {1}; int it : varexobs)
+        output << "options_.varexobs(" << ic++ << ")  = {'" << getName(it) << "'};" << endl;
 
-  output << endl
-         << "int exo_nbr = " << exo_nbr() << ";" << endl;
-  if (exo_nbr() > 0)
-    {
-      output << "char *exo_names[" << exo_nbr() << "];" << endl;
-      for (int id = 0; id < exo_nbr(); id++)
-        output << "exo_names[" << id << "] = \"" << getName(exo_ids[id]) << "\";" << endl;
+      output << "options_.varexobs_id = [ ";
+      for (int varexob : varexobs)
+        output << getTypeSpecificID(varexob) + 1 << " ";
+      output << " ];" << endl;
     }
 
-  output << endl
-         << "int exo_det_nbr = " << exo_det_nbr() << ";" << endl;
-  if (exo_det_nbr() > 0)
+  // 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)
         {
-      output << "char *exo_det_names[" << exo_det_nbr() << "];" << endl;
-      for (int id = 0; id < exo_det_nbr(); id++)
-        output << "exo_det_names[" << id << "] = \"" << getName(exo_det_ids[id]) << "\";" << endl;
+          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 << endl
-         << "int endo_nbr = " << endo_nbr() << ";" << endl;
-  if (endo_nbr() > 0)
+      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 << "char *endo_names[" << endo_nbr() << "];" << endl;
-      for (int id = 0; id < endo_nbr(); id++)
-        output << "endo_names[" << id << "] = \"" << getName(endo_ids[id]) << "\";" << endl;
+            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;
+          }
+    }
 }
 
-  output << endl
-         << "int param_nbr = " << param_nbr() << ";" << endl;
-  if (param_nbr() > 0)
+int
+SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index, expr_t expr_arg) noexcept(false)
 {
-      output << "char *param_names[" << param_nbr() << "];" << endl;
-      for (int id = 0; id < param_nbr(); id++)
-        output << "param_names[" << id << "] = \"" << getName(param_ids[id]) << "\";" << endl;
+  string varname {(endo ? "AUX_ENDO_LEAD_" : "AUX_EXO_LEAD_") + to_string(index)};
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
     }
 
-  // Write the auxiliary variable table
-  output << "int aux_var_nbr = " << aux_vars.size() << ";" << endl;
-  if (aux_vars.size() > 0)
+  aux_vars.emplace_back(symb_id, (endo ? AuxVarType::endoLead : AuxVarType::exoLead), 0, 0, 0, 0,
+                        expr_arg, "");
+
+  return symb_id;
+}
+
+int
+SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag,
+                                        expr_t expr_arg) noexcept(false)
 {
-      output << "struct aux_vars_t *av[" << aux_vars.size() << "];" << endl;
-      for (int i = 0; i < (int) aux_vars.size(); i++)
+  string varname {(endo ? "AUX_ENDO_LAG_" : "AUX_EXO_LAG_") + to_string(orig_symb_id) + "_"
+                  + to_string(-orig_lead_lag)};
+  int symb_id;
+  try
     {
-          output << "av[" << i << "].endo_index = " << getTypeSpecificID(aux_vars[i].get_symb_id()) << ";" << endl
-                 << "av[" << i << "].type = " << aux_vars[i].get_type() << ";" << endl;
-          switch (aux_vars[i].get_type())
+      symb_id = addSymbol(varname, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException& e)
     {
-            case avEndoLead:
-            case avExoLead:
-            case avExpectation:
-            case avMultiplier:
-            case avDiffForward:
-              break;
-            case avEndoLag:
-            case avExoLag:
-              output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
-                     << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
-              break;
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
     }
+
+  aux_vars.emplace_back(symb_id, (endo ? AuxVarType::endoLag : AuxVarType::exoLag), orig_symb_id,
+                        orig_lead_lag, 0, 0, expr_arg, "");
+
+  return symb_id;
 }
+
+int
+SymbolTable::addEndoLeadAuxiliaryVar(int index, expr_t expr_arg) noexcept(false)
+{
+  return addLeadAuxiliaryVarInternal(true, index, expr_arg);
 }
 
-  output << "int predeterminedNbr = " << predeterminedNbr() << ";" << endl;
-  if (predeterminedNbr() > 0)
+int
+SymbolTable::addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag,
+                                    expr_t expr_arg) noexcept(false)
 {
-      output << "int predetermined_variables[" << predeterminedNbr() << "] = {";
-      for (set<int>::const_iterator it = predetermined_variables.begin();
-           it != predetermined_variables.end(); it++)
+  return addLagAuxiliaryVarInternal(true, orig_symb_id, orig_lead_lag, expr_arg);
+}
+
+int
+SymbolTable::addExoLeadAuxiliaryVar(int index, expr_t expr_arg) noexcept(false)
 {
-          if (it != predetermined_variables.begin())
-            output << ",";
-          output << getTypeSpecificID(*it);
+  return addLeadAuxiliaryVarInternal(false, index, expr_arg);
 }
-      output << "};" << endl;
+
+int
+SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag,
+                                   expr_t expr_arg) noexcept(false)
+{
+  return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lead_lag, expr_arg);
 }
 
-  output << "int observedVariablesNbr = " << observedVariablesNbr() << ";" << endl;
-  if (observedVariablesNbr() > 0)
+int
+SymbolTable::addExpectationAuxiliaryVar(int information_set, int index,
+                                        expr_t expr_arg) noexcept(false)
 {
-      output << "int varobs[" << observedVariablesNbr() << "] = {";
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++)
+  string varname {"AUX_EXPECT_"s + (information_set < 0 ? "LAG" : "LEAD") + "_"
+                  + to_string(abs(information_set)) + "_" + to_string(index)};
+  int symb_id;
+  try
     {
-          if (it != varobs.begin())
-            output << ",";
-          output << getTypeSpecificID(*it);
-        }
-      output << "};" << endl;
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
     }
 
-void
-SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
-{
-  if (!frozen)
-    throw NotYetFrozenException();
+  aux_vars.emplace_back(symb_id, AuxVarType::expectation, 0, 0, 0, information_set, expr_arg, "");
 
-  output << endl
-         << "exo_nbr = " << exo_nbr() << ";" << endl;
-  for (int id = 0; id < exo_nbr(); id++)
-    output << "exo_names[\"" << getName(exo_ids[id]) << "\"] = " << id << ";" << endl;
+  return symb_id;
+}
 
-  output << endl
-         << "exo_det_nbr = " << exo_det_nbr() << ";" << endl;
-  for (int id = 0; id < exo_det_nbr(); id++)
-    output << "exo_det_names[\"" << getName(exo_det_ids[id]) << "\"] = " << id << " ;" << endl;
+int
+SymbolTable::addLogTransformAuxiliaryVar(int orig_symb_id, int orig_lead_lag,
+                                         expr_t expr_arg) noexcept(false)
+{
+  string varname = "LOG_" + getName(orig_symb_id);
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", it conflicts with the auxiliary variable created for representing the log of "
+           << getName(orig_symb_id) << endl;
+      exit(EXIT_FAILURE);
+    }
 
-  output << endl
-         << "endo_nbr = " << endo_nbr() << ";" << endl;
-  for (int id = 0; id < endo_nbr(); id++)
-    output << "endo_names[\"" << getName(endo_ids[id]) << "\"] = " << id << ";" << endl;
+  aux_vars.emplace_back(symb_id, AuxVarType::logTransform, orig_symb_id, orig_lead_lag, 0, 0,
+                        expr_arg, "");
 
-  output << endl
-         << "param_nbr = " << param_nbr() << ";" << endl;
-  for (int id = 0; id < param_nbr(); id++)
-    output << "param_names[\"" << getName(param_ids[id]) << "\"] = " << id << ";" << endl;
+  return symb_id;
+}
 
-  // Write the auxiliary variable table
-  for (int i = 0; i < (int) aux_vars.size(); i++)
-    {
-      output << "aux_vars_t av" << i << ";" << endl;
-      output << "av" << i << ".endo_index = " << getTypeSpecificID(aux_vars[i].get_symb_id()) << ";" << endl
-             << "av" << i << ".type = " << aux_vars[i].get_type() << ";" << endl;
-      switch (aux_vars[i].get_type())
-        {
-        case avEndoLead:
-        case avExoLead:
-        case avExpectation:
-        case avMultiplier:
-        case avDiffForward:
-          break;
-        case avEndoLag:
-        case avExoLag:
-          output << "av" << i << ".orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
-                 << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
-          break;
+int
+SymbolTable::addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id,
+                                    int orig_lag) noexcept(false)
+{
+  string varname {"AUX_DIFF_LAG_" + to_string(index)};
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
-      output << "aux_vars.push_back(" << "av" << i << ");" << endl;
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
     }
 
-  for (set<int>::const_iterator it = predetermined_variables.begin();
-       it != predetermined_variables.end(); it++)
-    output << "predetermined_variables.push_back(" << getTypeSpecificID(*it) << ");" << endl;
+  aux_vars.emplace_back(symb_id, AuxVarType::diffLag, orig_symb_id, orig_lag, 0, 0, expr_arg, "");
 
-  for (vector<int>::const_iterator it = varobs.begin();
-       it != varobs.end(); it++)
-    output << "varobs.push_back(" << getTypeSpecificID(*it) << ");" << endl;
+  return symb_id;
 }
 
 int
-SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addDiffLeadAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id,
+                                     int orig_lead) noexcept(false)
 {
-  ostringstream varname;
-  if (endo)
-    varname << "AUX_ENDO_LEAD_";
-  else
-    varname << "AUX_EXO_LEAD_";
-  varname << index;
+  string varname {"AUX_DIFF_LEAD_" + to_string(index)};
   int symb_id;
   try
     {
-      symb_id = addSymbol(varname.str(), eEndogenous);
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
   catch (AlreadyDeclaredException& e)
     {
-      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLead : avExoLead), 0, 0, 0, 0, expr_arg));
+  aux_vars.emplace_back(symb_id, AuxVarType::diffLead, orig_symb_id, orig_lead, 0, 0, expr_arg, "");
 
   return symb_id;
 }
 
 int
-SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, const optional<int>& orig_symb_id,
+                                 const optional<int>& orig_lag) noexcept(false)
 {
-  ostringstream varname;
-  if (endo)
-    varname << "AUX_ENDO_LAG_";
-  else
-    varname << "AUX_EXO_LAG_";
-  varname << orig_symb_id << "_" << -orig_lead_lag;
-
+  string varname {"AUX_DIFF_" + to_string(index)};
   int symb_id;
   try
     {
-      symb_id = addSymbol(varname.str(), eEndogenous);
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
   catch (AlreadyDeclaredException& e)
     {
-      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLag : avExoLag), orig_symb_id, orig_lead_lag, 0, 0, expr_arg));
+  aux_vars.emplace_back(symb_id, AuxVarType::diff, orig_symb_id, orig_lag, 0, 0, expr_arg, "");
 
   return symb_id;
 }
 
 int
-SymbolTable::addEndoLeadAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op,
+                                    const optional<int>& orig_symb_id,
+                                    const optional<int>& orig_lag) noexcept(false)
 {
-  return addLeadAuxiliaryVarInternal(true, index, expr_arg);
+  string varname {"AUX_UOP_" + to_string(index)};
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
-
-int
-SymbolTable::addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
+  catch (AlreadyDeclaredException& e)
     {
-  return addLagAuxiliaryVarInternal(true, orig_symb_id, orig_lead_lag, expr_arg);
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
     }
 
-int
-SymbolTable::addExoLeadAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
-{
-  return addLeadAuxiliaryVarInternal(false, index, expr_arg);
+  aux_vars.emplace_back(symb_id, AuxVarType::unaryOp, orig_symb_id, orig_lag, 0, 0, expr_arg,
+                        move(unary_op));
+
+  return symb_id;
 }
 
 int
-SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addHeterogeneousMultiplierAuxiliaryVar(int het_dim, int index,
+                                                    const string& varname) noexcept(false)
 {
-  return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lead_lag, expr_arg);
+  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::addExpectationAuxiliaryVar(int information_set, int index, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addMultiplierAuxiliaryVar(int index) noexcept(false)
 {
-  ostringstream varname;
+  string varname {"MULT_" + to_string(index + 1)};
   int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException& e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
 
-  varname << "AUX_EXPECT_" << (information_set < 0 ? "LAG" : "LEAD") << "_"
-          << abs(information_set) << "_" << index;
+  aux_vars.emplace_back(symb_id, AuxVarType::multiplier, 0, 0, index, 0, nullptr, "");
+  return symb_id;
+}
 
+int
+SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id, int orig_lead_lag,
+                                        expr_t expr_arg) noexcept(false)
+{
+  string varname {"AUX_DIFF_FWRD_" + to_string(orig_symb_id + 1)};
+  int symb_id;
   try
     {
-      symb_id = addSymbol(varname.str(), eEndogenous);
+      symb_id = addSymbol(varname, SymbolType::endogenous);
     }
   catch (AlreadyDeclaredException& e)
     {
-      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      cerr << "ERROR: you should rename your variable called " << varname
+           << ", this name is internally used by Dynare" << endl;
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, avExpectation, 0, 0, 0, information_set, expr_arg));
-
+  aux_vars.emplace_back(symb_id, AuxVarType::diffForward, orig_symb_id, orig_lead_lag, 0, 0,
+                        expr_arg, "");
   return symb_id;
 }
 
 int
-SymbolTable::addMultiplierAuxiliaryVar(int index) throw (FrozenException)
+SymbolTable::addPacExpectationAuxiliaryVar(const string& name, expr_t expr_arg)
 {
-  ostringstream varname;
   int symb_id;
-  varname << "MULT_" << index+1;
-
   try
     {
-      symb_id = addSymbol(varname.str(), eEndogenous);
+      symb_id = addSymbol(name, SymbolType::endogenous);
     }
   catch (AlreadyDeclaredException& e)
     {
-      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      cerr << "ERROR: the variable/parameter '" << name
+           << "' conflicts with a variable that will be generated for a 'pac_expectation' "
+              "expression. Please rename it."
+           << endl;
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, avMultiplier, 0, 0, index, 0, NULL));
+  aux_vars.emplace_back(symb_id, AuxVarType::pacExpectation, 0, 0, 0, 0, expr_arg, "");
   return symb_id;
 }
 
 int
-SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id, expr_t expr_arg) throw (FrozenException)
+SymbolTable::addPacTargetNonstationaryAuxiliaryVar(const string& name, expr_t expr_arg)
 {
-  ostringstream varname;
   int symb_id;
-  varname << "AUX_DIFF_FWRD_" << orig_symb_id+1;
-
   try
     {
-      symb_id = addSymbol(varname.str(), eEndogenous);
+      symb_id = addSymbol(name, SymbolType::endogenous);
     }
   catch (AlreadyDeclaredException& e)
     {
-      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      cerr << "ERROR: the variable/parameter '" << name
+           << "' conflicts with a variable that will be generated for a 'pac_target_nonstationary' "
+              "expression. Please rename it."
+           << endl;
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, avDiffForward, orig_symb_id, 0, 0, 0, expr_arg));
+  aux_vars.emplace_back(symb_id, AuxVarType::pacTargetNonstationary, 0, 0, 0, 0, expr_arg, "");
   return symb_id;
 }
 
 int
-SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException)
+SymbolTable::addAggregationOpAuxiliaryVar(const string& name, expr_t expr_arg)
+{
+  int symb_id {[&] {
+    try
       {
-  for (size_t i = 0; i < aux_vars.size(); i++)
-    if ((aux_vars[i].get_type() == avEndoLag || aux_vars[i].get_type() == avExoLag)
-        && aux_vars[i].get_orig_symb_id() == orig_symb_id && aux_vars[i].get_orig_lead_lag() == orig_lead_lag)
-      return aux_vars[i].get_symb_id();
-  throw SearchFailedException(orig_symb_id, orig_lead_lag);
+        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);
+      }
+  }()};
 
-expr_t
-SymbolTable::getAuxiliaryVarsExprNode(int symb_id) const throw (SearchFailedException)
-// throw exception if it is a Lagrange multiplier
+  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)
 {
-  for (size_t i = 0; i < aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id)
+  for (const auto& aux_var : aux_vars)
+    if ((aux_var.type == AuxVarType::endoLag || aux_var.type == AuxVarType::exoLag)
+        && aux_var.orig_symb_id == orig_symb_id && aux_var.orig_lead_lag == orig_lead_lag)
+      return aux_var.symb_id;
+  throw SearchFailedException {orig_symb_id, orig_lead_lag};
+}
+
+int
+SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id_arg) const noexcept(false)
 {
-        expr_t expr_node = aux_vars[i].get_expr_node();
-        if (expr_node != NULL)
-          return expr_node;
+  for (const auto& aux_var : aux_vars)
+    if ((aux_var.type == AuxVarType::endoLag || aux_var.type == AuxVarType::exoLag
+         || aux_var.type == AuxVarType::diff || aux_var.type == AuxVarType::diffLag
+         || aux_var.type == AuxVarType::diffLead || aux_var.type == AuxVarType::diffForward
+         || aux_var.type == AuxVarType::unaryOp)
+        && aux_var.symb_id == aux_var_symb_id_arg)
+      {
+        if (optional<int> r = aux_var.orig_symb_id; r)
+          return *r;
         else
-          throw SearchFailedException(symb_id);
+          throw UnknownSymbolIDException {
+              aux_var_symb_id_arg}; // Some diff and unaryOp auxvars have orig_symb_id unset
+      }
+  throw UnknownSymbolIDException {aux_var_symb_id_arg};
 }
-  throw SearchFailedException(symb_id);
+
+pair<int, int>
+SymbolTable::unrollDiffLeadLagChain(int symb_id, int lag) const noexcept(false)
+{
+  for (const auto& aux_var : aux_vars)
+    if (aux_var.symb_id == symb_id)
+      if (aux_var.type == AuxVarType::diffLag || aux_var.type == AuxVarType::diffLead)
+        {
+          auto [orig_symb_id, orig_lag] = unrollDiffLeadLagChain(aux_var.orig_symb_id.value(), lag);
+          return {orig_symb_id, orig_lag + aux_var.orig_lead_lag.value()};
+        }
+  return {symb_id, lag};
 }
 
 void
-SymbolTable::markPredetermined(int symb_id) throw (UnknownSymbolIDException, FrozenException)
+SymbolTable::markPredetermined(int symb_id) noexcept(false)
 {
   validateSymbID(symb_id);
 
   if (frozen)
     throw FrozenException();
 
-  assert(getType(symb_id) == eEndogenous);
+  assert(getType(symb_id) == SymbolType::endogenous);
 
   predetermined_variables.insert(symb_id);
 }
 
+void
+SymbolTable::markWithLogTransform(int symb_id) noexcept(false)
+{
+  validateSymbID(symb_id);
+
+  if (frozen)
+    throw FrozenException();
+
+  assert(getType(symb_id) == SymbolType::endogenous);
+
+  with_log_transform.insert(symb_id);
+}
+
 bool
-SymbolTable::isPredetermined(int symb_id) const throw (UnknownSymbolIDException)
+SymbolTable::isPredetermined(int symb_id) const noexcept(false)
 {
   validateSymbID(symb_id);
-  return (predetermined_variables.find(symb_id) != predetermined_variables.end());
+  return predetermined_variables.contains(symb_id);
 }
 
 int
 SymbolTable::predeterminedNbr() const
 {
-  return (predetermined_variables.size());
+  return predetermined_variables.size();
 }
 
 void
-SymbolTable::addObservedVariable(int symb_id) throw (UnknownSymbolIDException)
+SymbolTable::addObservedVariable(int symb_id) noexcept(false)
 {
   validateSymbID(symb_id);
-  assert(getType(symb_id) == eEndogenous);
+  assert(getType(symb_id) == SymbolType::endogenous);
   varobs.push_back(symb_id);
 }
 
 int
 SymbolTable::observedVariablesNbr() const
 {
-  return (int) varobs.size();
+  return static_cast<int>(varobs.size());
 }
 
 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
 {
-  vector<int>::const_iterator it = find(varobs.begin(), varobs.end(), symb_id);
+  auto it = ranges::find(varobs, symb_id);
   assert(it != varobs.end());
-  return (int) (it - varobs.begin());
+  return static_cast<int>(it - varobs.begin());
+}
+
+void
+SymbolTable::addObservedExogenousVariable(int symb_id) noexcept(false)
+{
+  validateSymbID(symb_id);
+  assert(getType(symb_id) == SymbolType::exogenous);
+  varexobs.push_back(symb_id);
+}
+
+int
+SymbolTable::observedExogenousVariablesNbr() const
+{
+  return static_cast<int>(varexobs.size());
+}
+
+bool
+SymbolTable::isObservedExogenousVariable(int symb_id) const
+{
+  return ranges::find(varexobs, symb_id) != varexobs.end();
+}
+
+int
+SymbolTable::getObservedExogenousVariableIndex(int symb_id) const
+{
+  auto it = ranges::find(varexobs, symb_id);
+  assert(it != varexobs.end());
+  return static_cast<int>(it - varexobs.begin());
 }
 
 vector<int>
 SymbolTable::getTrendVarIds() const
 {
   vector<int> trendVars;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eTrend || getType(it->second) == eLogTrend)
-      trendVars.push_back(it->second);
+  for (const auto& it : symbol_table)
+    if (getType(it.second) == SymbolType::trend || getType(it.second) == SymbolType::logTrend)
+      trendVars.push_back(it.second);
   return trendVars;
 }
 
@@ -792,180 +1007,190 @@ set<int>
 SymbolTable::getExogenous() const
 {
   set<int> exogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eExogenous)
-      exogs.insert(it->second);
+  for (const auto& it : symbol_table)
+    if (getType(it.second) == SymbolType::exogenous)
+      exogs.insert(it.second);
   return exogs;
 }
 
+set<int>
+SymbolTable::getObservedExogenous() const
+{
+  set<int> oexogs;
+  for (const auto& it : symbol_table)
+    if (getType(it.second) == SymbolType::exogenous)
+      if (isObservedExogenousVariable(it.second))
+        oexogs.insert(it.second);
+  return oexogs;
+}
+
 set<int>
 SymbolTable::getEndogenous() const
 {
   set<int> endogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eEndogenous)
-      endogs.insert(it->second);
+  for (const auto& it : symbol_table)
+    if (getType(it.second) == SymbolType::endogenous)
+      endogs.insert(it.second);
   return endogs;
 }
 
 bool
 SymbolTable::isAuxiliaryVariable(int symb_id) const
 {
-  for (int i = 0; i < (int) aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id)
-      return true;
-  return false;
+  return ranges::any_of(aux_vars, [=](const auto& av) { return av.symb_id == symb_id; });
 }
 
 bool
-SymbolTable::isAuxiliaryVariableButNotMultiplier(int symb_id) const
+SymbolTable::isDiffAuxiliaryVariable(int symb_id) const
 {
-  for (int i = 0; i < (int) aux_vars.size(); i++)
-    if (aux_vars[i].get_symb_id() == symb_id && aux_vars[i].get_type() != avMultiplier)
-      return true;
-  return false;
+  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);
+  });
 }
 
 set<int>
 SymbolTable::getOrigEndogenous() const
 {
   set<int> origendogs;
-  for (symbol_table_type::const_iterator it = symbol_table.begin();
-       it != symbol_table.end(); it++)
-    if (getType(it->second) == eEndogenous && !isAuxiliaryVariable(it->second))
-      origendogs.insert(it->second);
+  for (const auto& it : symbol_table)
+    if (getType(it.second) == SymbolType::endogenous && !isAuxiliaryVariable(it.second))
+      origendogs.insert(it.second);
   return origendogs;
 }
 
 void
-SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenException)
+SymbolTable::writeJsonOutput(ostream& output) const
 {
-  if (!frozen)
-    throw NotYetFrozenException();
+  output << R"("endogenous": )";
+  writeJsonVarVector(output, endo_ids);
 
-  output << "# Endogenous Variables" << endl
-         << "model_.endo = [" << endl;
-  if (endo_nbr() > 0)
-    for (int id = 0; id < endo_nbr(); id++)
-      output << "              DynareModel.Endo(\""
-             << getName(endo_ids[id]) << "\", \""
-             << getTeXName(endo_ids[id]) << "\", \""
-             << getLongName(endo_ids[id]) << "\")" << endl;
-  output << "             ]" << endl;
+  output << R"(, "exogenous":)";
+  writeJsonVarVector(output, exo_ids);
 
-  output << "# Exogenous Variables" << endl
-         << "model_.exo = [" << endl;
-  if (exo_nbr() > 0)
-    for (int id = 0; id < exo_nbr(); id++)
-      output << "             DynareModel.Exo(\""
-             << getName(exo_ids[id]) << "\", \""
-             << getTeXName(exo_ids[id]) << "\", \""
-             << getLongName(exo_ids[id]) << "\")" << endl;
+  output << R"(, "exogenous_deterministic": )";
+  writeJsonVarVector(output, exo_det_ids);
+
+  output << R"(, "parameters": )";
+  writeJsonVarVector(output, param_ids);
+
+  if (observedVariablesNbr() > 0)
+    {
+      output << R"(, "varobs": [)";
+      for (size_t i = 0; i < varobs.size(); i++)
+        {
+          if (i != 0)
+            output << ", ";
+          output << R"(")" << getName(varobs[i]) << R"(")";
+        }
       output << "]" << endl;
 
-  if (exo_det_nbr() > 0)
+      output << R"(, "varobs_ids": [)";
+      for (size_t i = 0; i < varobs.size(); i++)
         {
-      output << "# Exogenous Deterministic Variables" << endl
-             << "model_.exo_det = [" << endl;
-      if (exo_det_nbr() > 0)
-        for (int id = 0; id < exo_det_nbr(); id++)
-          output << "                 DynareModel.ExoDet(\""
-                 << getName(exo_det_ids[id]) << "\", \""
-                 << getTeXName(exo_det_ids[id]) << "\", \""
-                 << getLongName(exo_det_ids[id]) << "\")" << endl;
+          if (i != 0)
+            output << ", ";
+          output << getTypeSpecificID(varobs[i]) + 1;
+        }
       output << "]" << endl;
     }
 
-  output << "# Parameters" << endl
-         << "model_.param = [" << endl;
-  if (param_nbr() > 0)
-    for (int id = 0; id < param_nbr(); id++)
-      output << "               DynareModel.Param(\""
-             << getName(param_ids[id]) << "\", \""
-             << getTeXName(param_ids[id]) << "\", \""
-             << getLongName(param_ids[id]) << "\")" << endl;
+  if (observedExogenousVariablesNbr() > 0)
+    {
+      output << R"(, "varexobs": [)";
+      for (size_t i = 0; i < varexobs.size(); i++)
+        {
+          if (i != 0)
+            output << ", ";
+          output << R"(")" << getName(varexobs[i]) << R"(")";
+        }
       output << "]" << endl;
 
-  output << "model_.orig_endo_nbr = " << orig_endo_nbr() << endl;
-
-  if (aux_vars.size() > 0)
+      output << R"(, "varexobs_ids": [)";
+      for (size_t i = 0; i < varexobs.size(); i++)
         {
-      output << "# Auxiliary Variables" << endl
-             << "model_.aux_vars = [" << endl;
-      for (int i = 0; i < (int) aux_vars.size(); i++)
+          if (i != 0)
+            output << ", ";
+          output << getTypeSpecificID(varexobs[i]) + 1;
+        }
+      output << "]" << endl;
+    }
+  // Write the auxiliary variable table
+  output << R"(, "orig_endo_nbr": )" << orig_endo_nbr() << endl;
+  if (aux_vars.size() == 0)
+    output << R"(, "aux_vars": [])";
+  else
     {
-          output << "                   DynareModel.AuxVars("
-                 << getTypeSpecificID(aux_vars[i].get_symb_id()) + 1 << ", "
-                 << aux_vars[i].get_type() << ", ";
-          switch (aux_vars[i].get_type())
+      output << R"(, "aux_vars": [)" << endl;
+      for (int i = 0; i < static_cast<int>(aux_vars.size()); i++)
         {
-            case avEndoLead:
-            case avExoLead:
-              break;
-            case avEndoLag:
-            case avExoLag:
-              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
-                     << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
+          if (i != 0)
+            output << ", ";
+          output << R"({"endo_index": )" << getTypeSpecificID(aux_vars[i].symb_id) + 1
+                 << R"(, "type": )" << aux_vars[i].get_type_id();
+          switch (aux_vars[i].type)
+            {
+            case AuxVarType::endoLead:
+            case AuxVarType::exoLead:
+            case AuxVarType::expectation:
+            case AuxVarType::pacExpectation:
+            case AuxVarType::pacTargetNonstationary:
+            case AuxVarType::aggregationOp:
+            case AuxVarType::heterogeneousMultiplier:
               break;
-            case avMultiplier:
-              output << "NaN, NaN, " << aux_vars[i].get_equation_number_for_multiplier() + 1
-                     << ", NaN";
+            case AuxVarType::endoLag:
+            case AuxVarType::exoLag:
+            case AuxVarType::logTransform:
+            case AuxVarType::diffLag:
+            case AuxVarType::diffLead:
+            case AuxVarType::diffForward:
+              output << R"(, "orig_index": )"
+                     << getTypeSpecificID(aux_vars[i].orig_symb_id.value()) + 1
+                     << R"(, "orig_lead_lag": )" << aux_vars[i].orig_lead_lag.value();
               break;
-            case avDiffForward:
-              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ", NaN, ";
+            case AuxVarType::unaryOp:
+              output << R"(, "unary_op": ")" << aux_vars[i].unary_op << R"(")";
+              [[fallthrough]];
+            case AuxVarType::diff:
+              if (aux_vars[i].orig_symb_id)
+                output << R"(, "orig_index": )" << getTypeSpecificID(*aux_vars[i].orig_symb_id) + 1
+                       << R"(, "orig_lead_lag": )" << aux_vars[i].orig_lead_lag.value();
               break;
-            case avExpectation:
-              output << "NaN, NaN, NaN, \"\\mathbb{E}_{t"
-                     << (aux_vars[i].get_information_set() < 0 ? "" : "+")
-                     << aux_vars[i].get_information_set() << "}(";
-              aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
-              output << ")\"";
+            case AuxVarType::multiplier:
+              output << R"(, "eq_nbr": )" << aux_vars[i].equation_number_for_multiplier + 1;
               break;
             }
-          output << ")" << endl;
-        }
-      output << "]" << endl;
-    }
 
-  if (predeterminedNbr() > 0)
+          if (expr_t orig_expr = aux_vars[i].expr_node; orig_expr)
             {
-      output << "# Predetermined Variables" << endl
-             << "model_.pred_vars = [ " << endl;
-      for (set<int>::const_iterator it = predetermined_variables.begin();
-           it != predetermined_variables.end(); it++)
-        output << "                   DynareModel.PredVars("
-               << getTypeSpecificID(*it)+1 << ")" << endl;
+              output << R"(, "orig_expr": ")";
+              orig_expr->writeJsonOutput(output, {}, {});
+              output << R"(")";
+            }
+          output << '}' << endl;
+        }
       output << "]" << endl;
     }
 
-  if (observedVariablesNbr() > 0)
+  if (!heterogeneity_table.empty())
     {
-      output << "# Observed Variables" << endl
-             << "options_.obs_vars = [" << endl;
-      for (vector<int>::const_iterator it = varobs.begin();
-           it != varobs.end(); it++)
-        output << "                    DynareModel.ObsVars("
-               << getTypeSpecificID(*it)+1 << ")" << endl;
-      output << "                   ]" << endl;
+      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
-SymbolTable::writeJsonOutput(ostream &output) const
-{
-  output << "\"endogenous\": ";
-  writeJsonVarVector(output, endo_ids);
-
-  output << ", \"exogenous\":";
-  writeJsonVarVector(output, exo_ids);
-
-  output << ", \"exogenous_deterministic\": ";
-  writeJsonVarVector(output, exo_det_ids);
-
-  output << ", \"parameters\": ";
-  writeJsonVarVector(output, param_ids);
 }
 
 void
@@ -977,10 +1202,62 @@ SymbolTable::writeJsonVarVector(ostream &output, const vector<int> &varvec) cons
       if (i != 0)
         output << ", ";
       output << "{"
-             << "\"name\":\"" << getName(varvec[i]) << "\", "
-             << "\"texName\":\"" << boost::replace_all_copy(getTeXName(varvec[i]), "\\", "\\\\") << "\", "
-             << "\"longName\":\"" << boost::replace_all_copy(getLongName(varvec[i]), "\\", "\\\\") << "\"}"
-             << endl;
+             << R"("name":")" << getName(varvec[i]) << R"(", )"
+             << R"("texName":")" << boost::replace_all_copy(getTeXName(varvec[i]), R"(\)", R"(\\)")
+             << R"(", )"
+             << R"("longName":")"
+             << boost::replace_all_copy(getLongName(varvec[i]), R"(\)", R"(\\)") << R"("})" << endl;
     }
   output << "]" << endl;
 }
+
+int
+SymbolTable::getUltimateOrigSymbID(int symb_id) const
+{
+  while (isAuxiliaryVariable(symb_id))
+    try
+      {
+        symb_id = getOrigSymbIdForAuxVar(symb_id);
+      }
+    catch (UnknownSymbolIDException&)
+      {
+        break;
+      }
+  return symb_id;
+}
+
+optional<int>
+SymbolTable::getEquationNumberForMultiplier(int symb_id) const
+{
+  for (const auto& aux_var : aux_vars)
+    if (aux_var.symb_id == symb_id && aux_var.type == AuxVarType::multiplier)
+      return aux_var.equation_number_for_multiplier;
+  return nullopt;
+}
+
+const set<int>&
+SymbolTable::getVariablesWithLogTransform() const
+{
+  return with_log_transform;
+}
+
+set<int>
+SymbolTable::getLagrangeMultipliers() const
+{
+  set<int> r;
+  for (const auto& aux_var : aux_vars)
+    if (aux_var.type == AuxVarType::multiplier)
+      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};
+}