diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 0534695afacee1dbeecb66a01a8a9423961abc82..0ed0c9fb3880501a1f49297d3f4289e5ea9d4b0e 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -333,129 +333,6 @@ PriorPosteriorFunctionStatement::writeJsonOutput(ostream &output) const
   output << "}";
 }
 
-PacModelStatement::PacModelStatement(string name_arg,
-                                     string aux_model_name_arg,
-                                     string discount_arg,
-                                     expr_t growth_arg,
-                                     const SymbolTable &symbol_table_arg) :
-  name{move(name_arg)},
-  aux_model_name{move(aux_model_name_arg)},
-  discount{move(discount_arg)},
-  growth{growth_arg},
-  original_growth{growth_arg},
-  symbol_table{symbol_table_arg}
-{
-}
-
-void
-PacModelStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
-{
-  if (growth)
-    growth->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
-}
-
-void
-PacModelStatement::overwriteGrowth(expr_t new_growth)
-{
-  if (!new_growth || !growth)
-    return;
-
-  growth = new_growth;
-
-  try
-    {
-      growth_info = growth->matchLinearCombinationOfVariables(false);
-    }
-  catch (ExprNode::MatchFailureException &e)
-    {
-      auto gv = dynamic_cast<const VariableNode *>(growth);
-      if (gv)
-        growth_info.emplace_back(gv->symb_id, gv->lag, -1, 1);
-      else
-        {
-          cerr << "Pac growth must be a linear combination of variables" << endl;
-          exit(EXIT_FAILURE);
-        }
-    }
-}
-
-void
-PacModelStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
-{
-  output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name << "';" << endl
-         << "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount) + 1 << ";" << endl;
-
-  if (growth)
-    {
-      output << "M_.pac." << name << ".growth_str = '";
-      original_growth->writeJsonOutput(output, {}, {}, true);
-      output << "';" << endl;
-      int i = 0;
-      for (auto [growth_symb_id, growth_lag, param_id, constant] : growth_info)
-        {
-          string structname = "M_.pac." + name + ".growth_linear_comb(" + to_string(++i) + ").";
-          if (growth_symb_id >= 0)
-            {
-              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 == -1 ? 0 : symbol_table.getTypeSpecificID(param_id) + 1) << ";" << endl
-                 << structname << "constant = " << constant << ";" << endl;
-        }
-    }
-}
-
-void
-PacModelStatement::writeJsonOutput(ostream &output) const
-{
-  output << R"({"statementName": "pac_model",)"
-         << R"("model_name": ")" << name << R"(",)"
-         << R"("auxiliary_model_name": ")" << aux_model_name << R"(",)"
-         << R"("discount_index": )" << symbol_table.getTypeSpecificID(discount) + 1;
-  if (growth)
-    {
-      output << R"(,"growth_str": ")";
-      original_growth->writeJsonOutput(output, {}, {}, true);
-      output << R"(")";
-    }
-  output << "}" << endl;
-}
-
 StochSimulStatement::StochSimulStatement(SymbolList symbol_list_arg,
                                          OptionsList options_list_arg) :
   symbol_list{move(symbol_list_arg)},
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 0a79bc149e1a566f20a266d4b55cb5d4691b4df1..47cdb503c2eb1bddae4db76f3c849e1b07909715 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -140,26 +140,6 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
-class PacModelStatement : public Statement
-{
-public:
-  const string name, aux_model_name, discount;
-  expr_t growth, original_growth;
-private:
-  const SymbolTable &symbol_table;
-  vector<tuple<int, int, int, double>> growth_info;
-public:
-  PacModelStatement(string name_arg,
-                    string aux_model_name_arg,
-                    string discount_arg,
-                    expr_t growth_arg,
-                    const SymbolTable &symbol_table_arg);
-  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
-  void overwriteGrowth(expr_t new_growth);
-  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
-  void writeJsonOutput(ostream &output) const override;
-};
-
 class ForecastStatement : public Statement
 {
 private:
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 8717e739e1d989aa92feeaff5aacb889db500d22..d2c0b009115b787ffaf535e9c54fb944223a45aa 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5777,7 +5777,7 @@ DynamicModel::substituteUnaryOps(const vector<int> &eqnumbers)
 }
 
 pair<lag_equivalence_table_t, ExprNode::subst_table_t>
-DynamicModel::substituteDiff(vector<expr_t> &pac_growth)
+DynamicModel::substituteDiff(PacModelTable &pac_model_table)
 {
   /* Note: at this point, we know that there is no diff operator with a lead,
      because they have been expanded by DataTree::AddDiff().
@@ -5799,9 +5799,7 @@ DynamicModel::substituteDiff(vector<expr_t> &pac_growth)
   for (const auto &equation : equations)
     equation->findDiffNodes(diff_nodes);
 
-  for (const auto &gv : pac_growth)
-    if (gv)
-      gv->findDiffNodes(diff_nodes);
+  pac_model_table.findDiffNodesInGrowth(diff_nodes);
 
   // Substitute in model local variables
   vector<BinaryOpNode *> neweqs;
@@ -5817,9 +5815,7 @@ DynamicModel::substituteDiff(vector<expr_t> &pac_growth)
       equation = substeq;
     }
 
-  for (auto &it : pac_growth)
-    if (it)
-      it = it->substituteDiff(diff_nodes, diff_subst_table, neweqs);
+  pac_model_table.substituteDiffNodesInGrowth(diff_nodes, diff_subst_table, neweqs);
 
   // Add new equations
   for (auto &neweq : neweqs)
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index ee08a698d7d51d8ded6d1814362af52834695fed..2edbe0f96b73a0bcdb3e2b7edb53e52f2588114c 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -558,7 +558,7 @@ public:
   pair<lag_equivalence_table_t, ExprNode::subst_table_t> substituteUnaryOps(const vector<int> &eqnumbers);
 
   //! Substitutes diff operator
-  pair<lag_equivalence_table_t, ExprNode::subst_table_t> substituteDiff(vector<expr_t> &pac_growth);
+  pair<lag_equivalence_table_t, ExprNode::subst_table_t> substituteDiff(PacModelTable &pac_model_table);
 
   //! Substitute VarExpectation operators
   void substituteVarExpectation(const map<string, expr_t> &subst_table);
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 511f29c018f3368b512bbfca10d27e7c36c70a57..8c4f43b85375df5fdb1046b1dff7121c4246811d 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -34,6 +34,7 @@
 ModFile::ModFile(WarningConsolidation &warnings_arg)
   : var_model_table{symbol_table},
     trend_component_model_table{symbol_table},
+    pac_model_table{symbol_table},
     expressions_tree{symbol_table, num_constants, external_functions_table},
     original_model{symbol_table, num_constants, external_functions_table,
                    trend_component_model_table, var_model_table},
@@ -114,6 +115,8 @@ ModFile::checkPass(bool nostrict, bool stochastic)
   // Check epilogue block
   epilogue.checkPass(mod_file_struct, warnings);
 
+  pac_model_table.checkPass(mod_file_struct, warnings);
+
   if (mod_file_struct.write_latex_steady_state_model_present
       && !mod_file_struct.steady_state_model_present)
     {
@@ -409,17 +412,6 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
   if (unusedEndogsIsErr)
     exit(EXIT_FAILURE);
 
-  // Declare endogenous used for PAC model-consistent expectations
-  for (auto &statement : statements)
-    if (auto pms = dynamic_cast<PacModelStatement *>(statement.get());
-        pms)
-      {
-        if (pms->growth)
-          pac_growth.push_back(pms->growth);
-        if (pms->aux_model_name.empty())
-          dynamic_model.declarePacModelConsistentExpectationEndogs(pms->name);
-      }
-
   // Get all equation tags associated with VARs and Trend Component Models
   set<string> eqtags;
   for (const auto &it : trend_component_model_table.getEqTags())
@@ -440,7 +432,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
     tie(unary_ops_nodes, unary_ops_subst_table) = dynamic_model.substituteUnaryOps(eqtags);
 
   // Create auxiliary variable and equations for Diff operators
-  auto [diff_nodes, diff_subst_table] = dynamic_model.substituteDiff(pac_growth);
+  auto [diff_nodes, diff_subst_table] = dynamic_model.substituteDiff(pac_model_table);
 
   // Fill trend component and VAR model tables
   dynamic_model.fillTrendComponentModelTable();
@@ -449,55 +441,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
   dynamic_model.fillVarModelTable();
   original_model.fillVarModelTableFromOrigModel();
 
-  // Pac Model
-  int i = 0;
-  for (auto &statement : statements)
-    if (auto pms = dynamic_cast<PacModelStatement *>(statement.get()); pms)
-      {
-        if (pms->growth)
-          pms->overwriteGrowth(pac_growth.at(i++));
-
-        int max_lag;
-        vector<int> lhs;
-        vector<bool> nonstationary;
-        string aux_model_type;
-        if (trend_component_model_table.isExistingTrendComponentModelName(pms->aux_model_name))
-          {
-            aux_model_type = "trend_component";
-            max_lag = trend_component_model_table.getMaxLag(pms->aux_model_name) + 1;
-            lhs = dynamic_model.getUndiffLHSForPac(pms->aux_model_name, diff_subst_table);
-            // All lhs variables in a trend component model are nonstationary
-            nonstationary.insert(nonstationary.end(), trend_component_model_table.getDiff(pms->aux_model_name).size(), true);
-          }
-        else if (var_model_table.isExistingVarModelName(pms->aux_model_name))
-          {
-            aux_model_type = "var";
-            max_lag = var_model_table.getMaxLag(pms->aux_model_name);
-            lhs = var_model_table.getLhs(pms->aux_model_name);
-            // nonstationary variables in a VAR are those that are in diff
-            nonstationary = var_model_table.getDiff(pms->aux_model_name);
-          }
-        else if (pms->aux_model_name.empty())
-          max_lag = 0;
-        else
-          {
-            cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
-            exit(EXIT_FAILURE);
-          }
-        if (pms->growth)
-          dynamic_model.createPacGrowthNeutralityParameter(pms->name);
-
-        auto eqtag_and_lag = dynamic_model.walkPacParameters(pms->name);
-        if (pms->aux_model_name.empty())
-          dynamic_model.addPacModelConsistentExpectationEquation(pms->name, symbol_table.getID(pms->discount),
-                                                                 eqtag_and_lag, diff_subst_table,
-                                                                 pms->growth);
-        else
-          dynamic_model.fillPacModelInfo(pms->name, lhs, max_lag, aux_model_type,
-                                         eqtag_and_lag, nonstationary, pms->growth);
-        dynamic_model.substitutePacExpectation(pms->name);
-      }
-  dynamic_model.checkNoRemainingPacExpectation();
+  // PAC model
+  pac_model_table.transformPass(diff_subst_table, dynamic_model, var_model_table,
+                                trend_component_model_table);
 
   // Create auxiliary vars for Expectation operator
   dynamic_model.substituteExpectation(mod_file_struct.partial_information);
@@ -979,6 +925,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
 
   var_model_table.writeOutput(basename, mOutputFile);
   trend_component_model_table.writeOutput(basename, mOutputFile);
+  pac_model_table.writeOutput(basename, mOutputFile);
 
   // Initialize M_.Sigma_e, M_.Correlation_matrix, M_.H, and M_.Correlation_matrix_ME
   mOutputFile << "M_.Sigma_e = zeros(" << symbol_table.exo_nbr() << ", "
@@ -1270,6 +1217,12 @@ ModFile::writeJsonOutputParsingCheck(const string &basename, JsonFileOutputType
           output << ", ";
         }
 
+      if (!pac_model_table.empty())
+        {
+          pac_model_table.writeJsonOutput(output);
+          output << ", ";
+        }
+
       for (auto it = statements.begin(); it != statements.end(); ++it)
         {
           if (it != statements.begin())
@@ -1305,6 +1258,11 @@ ModFile::writeJsonOutputParsingCheck(const string &basename, JsonFileOutputType
               trend_component_model_table.writeJsonOutput(original_model_output);
               original_model_output << ", ";
             }
+          if (!pac_model_table.empty())
+            {
+              pac_model_table.writeJsonOutput(original_model_output);
+              original_model_output << ", ";
+            }
           int i = 0;
           for (const auto &it : statements)
             {
diff --git a/src/ModFile.hh b/src/ModFile.hh
index e2f31264202db9a44855941b69552cb3d94b52af..d5942c8447e6b173cc5030e8d8d8e015ab127546 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -55,6 +55,8 @@ public:
   VarModelTable var_model_table;
   //! Trend Component Model Table used for storing info about trend component models
   TrendComponentModelTable trend_component_model_table;
+  //! PAC Model Table used for storing info about pac models
+  PacModelTable pac_model_table;
   //! Expressions outside model block
   DataTree expressions_tree;
   //! Original model, as declared in the "model" block, that won't be modified by the preprocessor
@@ -124,7 +126,6 @@ private:
   void writeJsonOutputParsingCheck(const string &basename, JsonFileOutputType json_output_mode, bool transformpass, bool computingpass) const;
   void writeJsonComputingPassOutput(const string &basename, JsonFileOutputType json_output_mode, bool jsonderivsimple) const;
   void writeJsonFileHelper(const string &fname, ostringstream &output) const;
-  vector<expr_t> pac_growth;
   /* Generate a random temporary path, in the current directory. Equivalent to
      boost::filesystem::unique_path(). Both are insecure, but currently there
      is no better portable solution. Maybe in a later C++ standard? */
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 86b066b3ea1070a9d548ef86fcecaf0184c8327b..04e3541d358fadc484f55f2f45c3efc86a46ff62 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2621,9 +2621,7 @@ ParsingDriver::pac_model()
   auto discount = it->second;
   check_symbol_is_parameter(discount);
 
-  mod_file->addStatement(make_unique<PacModelStatement>(name, aux_model_name, discount,
-                                                        pac_growth,
-                                                        mod_file->symbol_table));
+  mod_file->pac_model_table.addPacModel(name, aux_model_name, discount, pac_growth);
   parsing_pac_model = false;
 }
 
diff --git a/src/SubModel.cc b/src/SubModel.cc
index a84ddfaa24b3a37ff3ca41b2d9c9ac15299d8aba..0866b4f2b19c3e10fcdb1282122424d18b23a93c 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -20,6 +20,7 @@
 #include <algorithm>
 
 #include "SubModel.hh"
+#include "DynamicModel.hh"
 
 TrendComponentModelTable::TrendComponentModelTable(SymbolTable &symbol_table_arg) :
   symbol_table{symbol_table_arg}
@@ -751,3 +752,225 @@ VarModelTable::getLhsExprT(const string &name_arg) const
   checkModelName(name_arg);
   return lhs_expr_t.find(name_arg)->second;
 }
+
+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)
+{
+  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;
+  names.insert(move(name_arg));
+}
+
+bool
+PacModelTable::isExistingPacModelName(const string &name_arg) const
+{
+  return names.find(name_arg) != names.end();
+}
+
+bool
+PacModelTable::empty() const
+{
+  return names.empty();
+}
+
+void
+PacModelTable::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+{
+  for (auto &[name, gv] : growth)
+    if (gv)
+      gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
+}
+
+void
+PacModelTable::findDiffNodesInGrowth(lag_equivalence_table_t &diff_nodes) const
+{
+  for (auto &[name, gv] : growth)
+    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);
+}
+
+void
+PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
+                             DynamicModel &dynamic_model, const VarModelTable &var_model_table,
+                             const TrendComponentModelTable &trend_component_model_table)
+{
+  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]->matchLinearCombinationOfVariables(false);
+          }
+        catch (ExprNode::MatchFailureException &e)
+          {
+            auto gv = dynamic_cast<const VariableNode *>(growth[name]);
+            if (gv)
+              growth_info[name].emplace_back(gv->symb_id, gv->lag, -1, 1);
+            else
+              {
+                cerr << "Pac growth must be a linear combination of variables" << endl;
+                exit(EXIT_FAILURE);
+              }
+          }
+
+      // Declare endogenous used for PAC model-consistent expectations
+      if (aux_model_name[name].empty())
+        dynamic_model.declarePacModelConsistentExpectationEndogs(name);
+
+      // Declare parameter for growth neutrality correction
+      if (growth[name])
+        dynamic_model.createPacGrowthNeutralityParameter(name);
+
+      // Substitute pac_expectation operators
+      int max_lag;
+      vector<int> lhs;
+      vector<bool> nonstationary;
+      string aux_model_type;
+      if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name[name]))
+        {
+          aux_model_type = "trend_component";
+          max_lag = trend_component_model_table.getMaxLag(aux_model_name[name]) + 1;
+          lhs = dynamic_model.getUndiffLHSForPac(aux_model_name[name], diff_subst_table);
+          // All lhs variables in a trend component model are nonstationary
+          nonstationary.insert(nonstationary.end(), trend_component_model_table.getDiff(aux_model_name[name]).size(), true);
+        }
+      else if (var_model_table.isExistingVarModelName(aux_model_name[name]))
+        {
+          aux_model_type = "var";
+          max_lag = var_model_table.getMaxLag(aux_model_name[name]);
+          lhs = var_model_table.getLhs(aux_model_name[name]);
+          // nonstationary variables in a VAR are those that are in diff
+          nonstationary = var_model_table.getDiff(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);
+        }
+
+      auto eqtag_and_lag = dynamic_model.walkPacParameters(name);
+      if (aux_model_name[name].empty())
+        dynamic_model.addPacModelConsistentExpectationEquation(name, symbol_table.getID(discount[name]),
+                                                               eqtag_and_lag, diff_subst_table,
+                                                               growth[name]);
+      else
+        dynamic_model.fillPacModelInfo(name, lhs, max_lag, aux_model_type,
+                                       eqtag_and_lag, nonstationary, growth[name]);
+
+      dynamic_model.substitutePacExpectation(name);
+    }
+
+  dynamic_model.checkNoRemainingPacExpectation();
+}
+
+void
+PacModelTable::writeOutput(const string &basename, ostream &output) const
+{
+  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;
+          int i = 0;
+          for (auto [growth_symb_id, growth_lag, param_id, constant] : growth_info.at(name))
+            {
+              string structname = "M_.pac." + name + ".growth_linear_comb(" + to_string(++i) + ").";
+              if (growth_symb_id >= 0)
+                {
+                  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 == -1 ? 0 : symbol_table.getTypeSpecificID(param_id) + 1) << ";" << endl
+                     << structname << "constant = " << constant << ";" << endl;
+            }
+        }
+    }
+}
+
+void
+PacModelTable::writeJsonOutput(ostream &output) const
+{
+  for (const auto &name : names)
+    {
+      if (name != *names.begin())
+        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;
+    }
+}
+
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 06af8284a1198c110f734451cd39ae4927a1fe3d..bb1b798c01d8b59f2c6ad6be32f35bbe182bcc80 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -28,6 +28,10 @@
 #include "ExprNode.hh"
 #include "SymbolTable.hh"
 #include "SymbolList.hh"
+#include "Statement.hh"
+
+// DynamicModel.hh can’t be included here, otherwise it would be a circular dependency
+class DynamicModel;
 
 using namespace std;
 
@@ -186,4 +190,32 @@ VarModelTable::empty() const
   return names.empty();
 }
 
+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
+  map<string, expr_t> growth, original_growth;
+  map<string, vector<tuple<int, int, int, double>>> growth_info;
+public:
+  explicit PacModelTable(SymbolTable &symbol_table_arg);
+  void addPacModel(string name_arg, string aux_model_name_arg, string discount_arg, expr_t growth_arg);
+  bool isExistingPacModelName(const string &name_arg) const;
+  bool empty() const;
+  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings);
+  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()
+  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(const string &basename, ostream &output) const;
+  void writeJsonOutput(ostream &output) const;
+};
+
+
 #endif