/*
 * 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.
 *
 * 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 (find(target_eqnums[it.first].begin(), target_eqnums[it.first].end(), 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(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), nontrend_it))));
      nontarget_lhs[name] = nontarget_lhs_vec;

      for (int trend_it : getTargetEqNums(name))
        target_lhs_vec.push_back(
            lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), trend_it))));
      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 (find(target_eqnums.at(name).begin(), target_eqnums.at(name).end(), 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 << distance(eqtags_vec.begin(), find(eqtags_vec.begin(), eqtags_vec.end(), 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>(
              distance(nontarget_lhs_vec.begin(),
                       find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), 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>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), 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>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), 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 = " << symbol_table.getTypeSpecificID(disc_var->symb_id) + 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 = find_if(
                    components.begin(), components.end(),
                    [&](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))
            {
              cerr << "ERROR: the block 'pac_target_info(" << name
                   << ")' is not supported in the context of a PAC model with model-consistent "
                      "expectations (MCE)."
                   << endl;
              exit(EXIT_FAILURE);
            }
          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.z1" : "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 = "
               << symbol_table.getTypeSpecificID(dynamic_cast<VariableNode*>(component)->symb_id)
                      + 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++;
      }
}

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();
}