diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index d2c0b009115b787ffaf535e9c54fb944223a45aa..a7393d1e081365d0fb16c83b643596547ec2ca2e 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -50,9 +50,6 @@ DynamicModel::copyHelper(const DynamicModel &m)
     blocks_derivatives_exo.emplace_back(convert_block_derivative(it));
   for (const auto &it : m.blocks_derivatives_exo_det)
     blocks_derivatives_exo_det.emplace_back(convert_block_derivative(it));
-
-  for (const auto &[key, expr] : m.pac_expectation_substitution)
-    pac_expectation_substitution.emplace(key, f(expr));
 }
 
 DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
@@ -108,15 +105,7 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
   blocks_jacob_cols_other_endo{m.blocks_jacob_cols_other_endo},
   blocks_jacob_cols_exo{m.blocks_jacob_cols_exo},
   blocks_jacob_cols_exo_det{m.blocks_jacob_cols_exo_det},
-  var_expectation_functions_to_write{m.var_expectation_functions_to_write},
-  pac_mce_alpha_symb_ids{m.pac_mce_alpha_symb_ids},
-  pac_h0_indices{m.pac_h0_indices},
-  pac_h1_indices{m.pac_h1_indices},
-  pac_mce_z1_symb_ids{m.pac_mce_z1_symb_ids},
-  pac_eqtag_and_lag{m.pac_eqtag_and_lag},
-  pac_growth_neutrality_params{m.pac_growth_neutrality_params},
-  pac_model_info{m.pac_model_info},
-  pac_equation_info{m.pac_equation_info}
+  var_expectation_functions_to_write{m.var_expectation_functions_to_write}
 {
   copyHelper(m);
 }
@@ -173,16 +162,6 @@ DynamicModel::operator=(const DynamicModel &m)
 
   var_expectation_functions_to_write = m.var_expectation_functions_to_write;
 
-  pac_mce_alpha_symb_ids = m.pac_mce_alpha_symb_ids;
-  pac_h0_indices = m.pac_h0_indices;
-  pac_h1_indices = m.pac_h1_indices;
-  pac_mce_z1_symb_ids = m.pac_mce_z1_symb_ids;
-  pac_eqtag_and_lag = m.pac_eqtag_and_lag;
-  pac_expectation_substitution.clear();
-  pac_growth_neutrality_params = m.pac_growth_neutrality_params;
-  pac_model_info = m.pac_model_info;
-  pac_equation_info = m.pac_equation_info;
-
   copyHelper(m);
 
   return *this;
@@ -3003,237 +2982,6 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl
   for (int i = 1; i < static_cast<int>(NNZDerivatives.size()); i++)
     output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
   output << "];" << endl;
-
-  // Write Pac Model Consistent Expectation parameter info
-  for (auto &it : pac_mce_alpha_symb_ids)
-    {
-      output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".mce.alpha = [";
-      for (auto it : it.second)
-        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
-      output << "];" << endl;
-    }
-
-  // Write Pac Model Consistent Expectation Z1 info
-  for (auto &it : pac_mce_z1_symb_ids)
-    output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".mce.z1 = "
-           << symbol_table.getTypeSpecificID(it.second) + 1 << ";" << endl;
-
-  // Write Pac lag info
-  for (auto &it : pac_eqtag_and_lag)
-    output << "M_.pac." << it.first.first << ".equations." << it.second.first << ".max_lag = " << it.second.second << ";" << endl;
-
-  // Write Pac equation tag info
-  map<string, vector<pair<string, string>>> for_writing;
-  for (auto &it : pac_eqtag_and_lag)
-    for_writing[it.first.first].emplace_back(it.first.second, it.second.first);
-
-  for (auto &it : for_writing)
-    {
-      output << "M_.pac." << it.first << ".tag_map = [";
-      for (auto &it1 : it.second)
-        output << "{'" << it1.first << "', '" << it1.second << "'};";
-      output << "];" << endl;
-    }
-
-  for (auto &[model, growth_neutrality_param_index] : pac_growth_neutrality_params)
-    output << "M_.pac." << model << ".growth_neutrality_param_index = "
-           << symbol_table.getTypeSpecificID(growth_neutrality_param_index) + 1 << ";" << endl;
-
-  for (auto &[model, values] : pac_model_info)
-    {
-      auto &[lhs, aux_model_type] = values;
-      output << "M_.pac." << model << ".lhs = [";
-      for (auto it : lhs)
-        output << it + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.pac." << model << ".auxiliary_model_type = '" << aux_model_type << "';" << endl;
-    }
-
-  for (auto &pit : pac_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] = pit.second;
-      string substruct = pit.first.first + ".equations." + pit.first.second + ".";
-
-      output << "M_.pac." << substruct << "lhs_var = "
-             << symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl;
-
-      if (optim_share_index >= 0)
-        output << "M_.pac." << substruct << "share_of_optimizing_agents_index = "
-               << symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
-
-      output << "M_.pac." << substruct << "ec.params = "
-             << symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
-             << "M_.pac." << substruct << "ec.vars = [";
-      for (auto it : ec_params_and_vars.second)
-        output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
-      output << "];" << endl
-             << "M_.pac." << substruct << "ec.istarget = [";
-      for (auto it : ec_params_and_vars.second)
-        output << (get<1>(it) ? "true " : "false ");
-      output << "];" << endl
-             << "M_.pac." << substruct << "ec.scale = [";
-      for (auto it : ec_params_and_vars.second)
-        output << get<2>(it) << " ";
-      output << "];" << endl
-             << "M_.pac." << substruct << "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." << substruct << "ar.params = [";
-      for (auto &[pid, vid, vlag] : ar_params_and_vars)
-        output << (pid != -1 ? symbol_table.getTypeSpecificID(pid) + 1 : -1) << " ";
-      output << "];" << endl
-             << "M_.pac." << substruct << "ar.vars = [";
-      for (auto &[pid, vid, vlag] : ar_params_and_vars)
-        output << (vid != -1 ? symbol_table.getTypeSpecificID(vid) + 1 : -1) << " ";
-      output << "];" << endl
-             << "M_.pac." << substruct << "ar.lags = [";
-      for (auto &[pid, vid, vlag] : ar_params_and_vars)
-        output << vlag << " ";
-      output << "];" << endl;
-      if (!non_optim_vars_params_and_constants.empty())
-        {
-          output << "M_.pac." << substruct << "non_optimizing_behaviour.params = [";
-          for (auto &it : non_optim_vars_params_and_constants)
-            if (get<2>(it) >= 0)
-              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
-            else
-              output << "NaN ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "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." << substruct << "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." << substruct << "non_optimizing_behaviour.lags = [";
-          for (auto &it : non_optim_vars_params_and_constants)
-            output << get<1>(it) << " ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "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." << substruct << "additive.params = [";
-          for (auto &it : additive_vars_params_and_constants)
-            if (get<2>(it) >= 0)
-              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
-            else
-              output << "NaN ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "additive.vars = [";
-          for (auto &it : additive_vars_params_and_constants)
-            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "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." << substruct << "additive.lags = [";
-          for (auto &it : additive_vars_params_and_constants)
-            output << get<1>(it) << " ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "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." << substruct << "optim_additive.params = [";
-          for (auto &it : optim_additive_vars_params_and_constants)
-            if (get<2>(it) >= 0)
-              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
-            else
-              output << "NaN ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "optim_additive.vars = [";
-          for (auto &it : optim_additive_vars_params_and_constants)
-            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "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." << substruct << "optim_additive.lags = [";
-          for (auto &it : optim_additive_vars_params_and_constants)
-            output << get<1>(it) << " ";
-          output << "];" << endl
-                 << "M_.pac." << substruct << "optim_additive.scaling_factor = [";
-          for (auto &it : optim_additive_vars_params_and_constants)
-            output << get<3>(it) << " ";
-          output << "];" << endl;
-        }
-      // Create empty h0 and h1 substructures that will be overwritten later if not empty
-      output << "M_.pac." << substruct << "h0_param_indices = [];" << endl
-             << "M_.pac." << substruct << "h1_param_indices = [];" << endl;
-    }
-
-  for (auto &it : pac_h0_indices)
-    {
-      output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".h0_param_indices = [";
-      for (auto it1 : it.second)
-        output << symbol_table.getTypeSpecificID(it1) + 1 << " ";
-      output << "];" << endl;
-    }
-
-  for (auto &it : pac_h1_indices)
-    {
-      output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".h1_param_indices = [";
-      for (auto it1 : it.second)
-        output << symbol_table.getTypeSpecificID(it1) + 1 << " ";
-      output << "];" << endl;
-    }
 }
 
 void
@@ -3899,11 +3647,9 @@ DynamicModel::getUndiffLHSForPac(const string &aux_model_name,
   return lhs;
 }
 
-map<pair<string, string>, pair<string, int>>
-DynamicModel::walkPacParameters(const string &name)
+void
+DynamicModel::analyzePacEquationStructure(const string &name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag, PacModelTable::equation_info_t &pac_equation_info)
 {
-  map<pair<string, string>, pair<string, int>> eqtag_and_lag;
-
   int i = 0;
   for (auto &equation : equations)
     {
@@ -3978,12 +3724,12 @@ DynamicModel::walkPacParameters(const string &name)
             }
           if (lhs.first == -1)
             {
-              cerr << "walkPacParameters: error obtaining LHS variable." << endl;
+              cerr << "analyzePacEquationStructure: error obtaining LHS variable." << endl;
               exit(EXIT_FAILURE);
             }
           if (ec_params_and_vars.second.empty())
             {
-              cerr << "walkPacParameters: error obtaining RHS parameters." << endl;
+              cerr << "analyzePacEquationStructure: error obtaining RHS parameters." << endl;
               exit(EXIT_FAILURE);
             }
           string eq = "eq" + to_string(i++);
@@ -3996,7 +3742,6 @@ DynamicModel::walkPacParameters(const string &name)
           eqtag_and_lag[{name, eqtag}] = {eq, max_lag};
         }
     }
-  return eqtag_and_lag;
 }
 
 int
@@ -4038,37 +3783,12 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
   throw PacTargetNotIdentifiedException{pac_model_name, "No equation with the corresponding pac_expectation operator"};
 }
 
-void
-DynamicModel::declarePacModelConsistentExpectationEndogs(const string &name)
-{
-  int i = 0;
-  for (auto &equation : equations)
-    if (equation->containsPacExpectation(name))
-      {
-        if (!equation_tags.exists(&equation - &equations[0], "name"))
-          {
-            cerr << "Every equation with a pac expectation must have been assigned an equation tag name" << endl;
-            exit(EXIT_FAILURE);
-          }
-        string standard_eqtag = "eq" + to_string(i++);
-        try
-          {
-            pac_mce_z1_symb_ids[{name, standard_eqtag}]
-              = symbol_table.addSymbol("mce_Z1_" + name + "_" + standard_eqtag, SymbolType::endogenous);
-          }
-        catch (SymbolTable::AlreadyDeclaredException &e)
-          {
-            cerr << "Variable name needed by PAC (mce_Z1_" << name << "_" << standard_eqtag << endl;
-            exit(EXIT_FAILURE);
-          }
-      }
-}
-
 void
 DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int discount_symb_id,
                                                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
                                                        ExprNode::subst_table_t &diff_subst_table,
-                                                       expr_t growth)
+                                                       map<pair<string, string>, int> &pac_mce_z1_symb_ids,
+                                                       map<pair<string, string>, vector<int>> &pac_mce_alpha_symb_ids)
 {
   int pac_target_symb_id;
   try
@@ -4080,43 +3800,53 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
       cerr << "Can't identify target for PAC model " << name << ": " << e.message;
       exit(EXIT_FAILURE);
     }
-  pac_eqtag_and_lag.insert(eqtag_and_lag.begin(), eqtag_and_lag.end());
   int neqs = 0;
-  for (auto &it : eqtag_and_lag)
+  for (const auto &[key, val] : eqtag_and_lag)
     {
-      assert(it.first.first == name);
-      string eqtag = it.first.second;
-      string standard_eqtag = it.second.first;
-      int pac_max_lag_m = it.second.second + 1;
+      auto &[name2, eqtag] = key;
+      if (name2 != name)
+        continue;
+
+      auto &standard_eqtag = val.first;
+      int pac_max_lag_m = val.second + 1;
       string append_to_name = name + "_" + standard_eqtag;
-      if (pac_mce_z1_symb_ids.find({name, standard_eqtag}) == pac_mce_z1_symb_ids.end())
+
+      // Create the endogenous representing Z₁
+      int mce_z1_symb_id;
+      string varname = "mce_Z1_" + append_to_name;
+      try
+        {
+          mce_z1_symb_id = symbol_table.addSymbol(varname, SymbolType::endogenous);
+        }
+      catch (SymbolTable::AlreadyDeclaredException &e)
         {
-          cerr << "Error finding pac MCE Z1 symb id" << endl;
+          cerr << "The variable/parameter '" << varname << "' conflicts with the variable that will be generated for the Z_1 variable of the '" << name << "' PAC model. Please rename it." << endl;
           exit(EXIT_FAILURE);
         }
-      int mce_z1_symb_id = pac_mce_z1_symb_ids[{name, standard_eqtag}];
+      pac_mce_z1_symb_ids[{name, standard_eqtag}] = mce_z1_symb_id;
 
       expr_t A = One;
       expr_t fp = Zero;
       expr_t beta = AddVariable(discount_symb_id);
       for (int i = 1; i <= pac_max_lag_m; i++)
-        try
-          {
-            int alpha_i_symb_id = symbol_table.addSymbol("mce_alpha_" + append_to_name + "_" + to_string(i),
-                                                         SymbolType::parameter);
-            pac_mce_alpha_symb_ids[{name, standard_eqtag}].push_back(alpha_i_symb_id);
-            A = AddPlus(A, AddVariable(alpha_i_symb_id));
-            fp = AddPlus(fp,
-                         AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
-                                           AddPower(beta, AddPossiblyNegativeConstant(i))),
-                                  AddVariable(mce_z1_symb_id, i)));
-
-          }
-        catch (SymbolTable::AlreadyDeclaredException &e)
-          {
-            cerr << "Variable name needed by PAC (mce_alpha_" << append_to_name << "_" << i << ")" << endl;
-            exit(EXIT_FAILURE);
-          }
+        {
+          string param_name = "mce_alpha_" + append_to_name + "_" + to_string(i);
+          try
+            {
+              int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+              pac_mce_alpha_symb_ids[{name, standard_eqtag}].push_back(alpha_i_symb_id);
+              A = AddPlus(A, AddVariable(alpha_i_symb_id));
+              fp = AddPlus(fp,
+                           AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
+                                             AddPower(beta, AddPossiblyNegativeConstant(i))),
+                                    AddVariable(mce_z1_symb_id, i)));
+            }
+          catch (SymbolTable::AlreadyDeclaredException &e)
+            {
+              cerr << "The variable/parameter '" << param_name << "' conflicts with a parameter that will be generated for the '" << name << "' PAC model. Please rename it." << endl;
+              exit(EXIT_FAILURE);
+            }
+        }
 
       // Add diff nodes and eqs for pac_target_symb_id
       const VariableNode *target_base_diff_node;
@@ -4173,14 +3903,14 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
           for (int j = k+1; j <= pac_max_lag_m; j++)
             {
               int alpha_j_symb_id = -1;
-              string varname = "mce_alpha_" + append_to_name + "_" + to_string(j);
+              string param_name = "mce_alpha_" + append_to_name + "_" + to_string(j);
               try
                 {
-                  alpha_j_symb_id = symbol_table.getID(varname);
+                  alpha_j_symb_id = symbol_table.getID(param_name);
                 }
               catch (SymbolTable::UnknownSymbolNameException &e)
                 {
-                  alpha_j_symb_id = symbol_table.addSymbol(varname, SymbolType::parameter);
+                  alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
                 }
               ssum = AddPlus(ssum,
                              AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j))));
@@ -4191,51 +3921,56 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
                             AddMinus(AddTimes(A, AddMinus(const_cast<VariableNode *>(target_base_diff_node), fs)), fp));
       addEquation(neweq, -1);
       neqs++;
-      /* The growth correction term is not added to the definition of Z₁
-         because the latter is recursive. Rather put it at the level of the
-         substition of pac_expectation operator. */
-      expr_t growth_correction = growth ? AddTimes(AddVariable(pac_growth_neutrality_params.at(name)), growth) : Zero;
-      pac_expectation_substitution[{name, eqtag}] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction);
     }
-  cout << "Pac Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl;
+  cout << "PAC Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl;
 }
 
 void
-DynamicModel::createPacGrowthNeutralityParameter(const string &pac_model_name)
+DynamicModel::computePacModelConsistentExpectationSubstitution(const string &name,
+                                                               const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                                               expr_t growth,
+                                                               const map<string, int> &pac_growth_neutrality_params,
+                                                               const map<pair<string, string>, int> &pac_mce_z1_symb_ids,
+                                                               map<pair<string, string>, expr_t> &pac_expectation_substitution)
 {
-  string param_name = pac_model_name + "_pac_growth_neutrality_correction";
-  try
+  for (const auto &[key, val] : eqtag_and_lag)
     {
-      int param_idx = symbol_table.addSymbol(param_name, SymbolType::parameter);
-      pac_growth_neutrality_params[pac_model_name] = param_idx;
-    }
-  catch (SymbolTable::AlreadyDeclaredException)
-    {
-      cerr << "The parameter '" << param_name << "' conflicts with the auxiliary parameter that will be generated for the growth neutrality correction of the '" << pac_model_name << "' PAC model. Please rename that parameter." << endl;
-      exit(EXIT_FAILURE);
+      auto &[name2, eqtag] = key;
+      if (name2 != name)
+        continue;
+
+      int mce_z1_symb_id = pac_mce_z1_symb_ids.at({name, val.first});
+
+      /* The growth correction term is not added to the definition of Z₁
+         because the latter is recursive. Rather put it at the level of the
+         substition of pac_expectation operator. */
+      expr_t growth_correction = growth ? AddTimes(AddVariable(pac_growth_neutrality_params.at(name)), growth) : Zero;
+      pac_expectation_substitution[{name, eqtag}] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction);
     }
 }
 
 void
-DynamicModel::fillPacModelInfo(const string &pac_model_name,
-                               vector<int> lhs,
-                               int max_lag,
-                               string aux_model_type,
-                               const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
-                               const vector<bool> &nonstationary,
-                               expr_t growth)
+DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
+                                                        const vector<int> &lhs,
+                                                        int max_lag,
+                                                        const string &aux_model_type,
+                                                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                                        const vector<bool> &nonstationary,
+                                                        expr_t growth,
+                                                        const map<string, int> &pac_growth_neutrality_params,
+                                                        map<pair<string, string>, vector<int>> &pac_h0_indices,
+                                                        map<pair<string, string>, vector<int>> &pac_h1_indices,
+                                                        map<pair<string, string>, expr_t> &pac_expectation_substitution)
 {
-  pac_eqtag_and_lag.insert(eqtag_and_lag.begin(), eqtag_and_lag.end());
-
   bool stationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), logical_not<bool>());
   bool nonstationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), [](bool b) { return b; }); // FIXME: use std::identity instead of an anonymous function when we upgrade to C++20
 
-  for (auto pac_models_and_eqtags : pac_eqtag_and_lag)
+  for (const auto &[key, val] : eqtag_and_lag)
     {
-      if (pac_models_and_eqtags.first.first != pac_model_name)
+      auto &[name2, eqtag] = key;
+      if (name2 != name)
         continue;
-      string eqtag = pac_models_and_eqtags.first.second;
-      string standard_eqtag = pac_models_and_eqtags.second.first;
+      auto &standard_eqtag = val.first;
       expr_t subExpr = Zero;
       if (stationary_vars_present)
         {
@@ -4243,21 +3978,21 @@ DynamicModel::fillPacModelInfo(const string &pac_model_name,
             {
               /* If the auxiliary model is a VAR, add a parameter corresponding
                  to the constant. */
-              string param_name_h0 = "h0_" + pac_model_name + "_" + standard_eqtag + "_constant";
+              string param_name_h0 = "h0_" + name + "_" + standard_eqtag + "_constant";
               int new_param_symb_id = symbol_table.addSymbol(param_name_h0, SymbolType::parameter);
-              pac_h0_indices[{pac_model_name, standard_eqtag}].push_back(new_param_symb_id);
+              pac_h0_indices[{name, standard_eqtag}].push_back(new_param_symb_id);
               subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id));
             }
           for (int i = 1; i < max_lag + 1; i++)
             for (auto lhsit : lhs)
               {
                 stringstream param_name_h0;
-                param_name_h0 << "h0_" << pac_model_name
+                param_name_h0 << "h0_" << name
                               << "_" << standard_eqtag
                               << "_var_" << symbol_table.getName(lhsit)
                               << "_lag_" << i;
                 int new_param_symb_id = symbol_table.addSymbol(param_name_h0.str(), SymbolType::parameter);
-                pac_h0_indices[{pac_model_name, standard_eqtag}].push_back(new_param_symb_id);
+                pac_h0_indices[{name, standard_eqtag}].push_back(new_param_symb_id);
                 subExpr = AddPlus(subExpr,
                                   AddTimes(AddVariable(new_param_symb_id),
                                            AddVariable(lhsit, -i)));
@@ -4270,21 +4005,21 @@ DynamicModel::fillPacModelInfo(const string &pac_model_name,
             {
               /* If the auxiliary model is a VAR, add a parameter corresponding
                  to the constant. */
-              string param_name_h1 = "h1_" + pac_model_name + "_" + standard_eqtag + "_constant";
+              string param_name_h1 = "h1_" + name + "_" + standard_eqtag + "_constant";
               int new_param_symb_id = symbol_table.addSymbol(param_name_h1, SymbolType::parameter);
-              pac_h1_indices[{pac_model_name, standard_eqtag}].push_back(new_param_symb_id);
+              pac_h1_indices[{name, standard_eqtag}].push_back(new_param_symb_id);
               subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id));
             }
           for (int i = 1; i < max_lag + 1; i++)
             for (auto lhsit : lhs)
               {
                 stringstream param_name_h1;
-                param_name_h1 << "h1_" << pac_model_name
+                param_name_h1 << "h1_" << name
                               << "_" << standard_eqtag
                               << "_var_" << symbol_table.getName(lhsit)
                               << "_lag_" << i;
                 int new_param_symb_id = symbol_table.addSymbol(param_name_h1.str(), SymbolType::parameter);
-                pac_h1_indices[{pac_model_name, standard_eqtag}].push_back(new_param_symb_id);
+                pac_h1_indices[{name, standard_eqtag}].push_back(new_param_symb_id);
                 subExpr = AddPlus(subExpr,
                                   AddTimes(AddVariable(new_param_symb_id),
                                            AddVariable(lhsit, -i)));
@@ -4293,29 +4028,26 @@ DynamicModel::fillPacModelInfo(const string &pac_model_name,
 
       if (growth)
         {
-          int growth_param_index = pac_growth_neutrality_params.at(pac_model_name);
+          int growth_param_index = pac_growth_neutrality_params.at(name);
           subExpr = AddPlus(subExpr,
                             AddTimes(AddVariable(growth_param_index), growth));
         }
 
-      pac_expectation_substitution[{pac_model_name, eqtag}] = subExpr;
+      pac_expectation_substitution[{name, eqtag}] = subExpr;
     }
-  pac_model_info[pac_model_name] = {move(lhs), move(aux_model_type)};
 }
 
 void
-DynamicModel::substitutePacExpectation(const string &pac_model_name)
+DynamicModel::substitutePacExpectation(const map<pair<string, string>, expr_t> &pac_expectation_substitution)
 {
-  for (auto &it : pac_expectation_substitution)
-    if (it.first.first == pac_model_name)
-      for (auto &equation : equations)
-        if (equation_tags.exists(&equation - &equations[0], "name", it.first.second))
-          {
-            auto substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(pac_model_name, it.second));
-            assert(substeq);
-            equation = substeq;
-            break;
-          }
+  for (auto &[key, substexpr] : pac_expectation_substitution)
+    {
+      auto &[model_name, eq_name] = key;
+      int eq = equation_tags.getEqnByTag("name", eq_name);
+      auto substeq = dynamic_cast<BinaryOpNode *>(equations[eq]->substitutePacExpectation(model_name, substexpr));
+      assert(substeq);
+      equations[eq] = substeq;
+    }
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 2edbe0f96b73a0bcdb3e2b7edb53e52f2588114c..d493837efc6c9eb470ecb8950dcf4b653f811809 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -122,35 +122,6 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int>> var_expectation_functions_to_write;
 
-  //! store symb_ids for alphas created in addPacModelConsistentExpectationEquation
-  //! (pac_model_name, standardized_eqtag) -> mce_alpha_symb_id
-  map<pair<string, string>, vector<int>> pac_mce_alpha_symb_ids;
-  //! store symb_ids for h0, h1 parameters
-  //! (pac_model_name, standardized_eqtag) -> parameter symb_ids
-  map<pair<string, string>, vector<int>> pac_h0_indices, pac_h1_indices;
-  //! store symb_ids for z1s created in addPacModelConsistentExpectationEquation
-  //! (pac_model_name, standardized_eqtag) -> mce_z1_symb_id
-  map<pair<string, string>, int> pac_mce_z1_symb_ids;
-  //! Store lag info for pac equations
-  //! (pac_model_name, equation_tag) -> (standardized_eqtag, lag)
-  map<pair<string, string>, pair<string, int>> pac_eqtag_and_lag;
-  // Store indices for growth neutrality parameters
-  // pac_model_name -> growth_param_index
-  map<string, int> pac_growth_neutrality_params;
-
-  //! (pac_model_name, equation_tag) -> expr_t
-  map<pair<string, string>, expr_t> pac_expectation_substitution;
-
-  //! Store info about backward PAC models:
-  //! pac_model_name -> (lhsvars, aux_model_type)
-  map<string, pair<vector<int>, string>> pac_model_info;
-
-  //! Store info about pac models specific to the equation they appear in
-  //! (pac_model_name, standardized_eqtag) ->
-  //!     (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants)
-  map<pair<string, string>,
-      tuple<pair<int, int>, int, vector<tuple<int, int, int>>, pair<int, vector<tuple<int, bool, int>>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>>> pac_equation_info;
-
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &basename) const;
   //! Writes dynamic model file (Julia version)
@@ -415,27 +386,6 @@ public:
   //! in the trend_component model
   void updateVarAndTrendModel() const;
 
-  //! Get Pac equation parameter info
-  map<pair<string, string>, pair<string, int>> walkPacParameters(const string &name);
-
-  // Create the growth neutrality parameter of a given PAC model (when the
-  // “growth” option has been passed)
-  void createPacGrowthNeutralityParameter(const string &pac_model_name);
-
-  //! Add var_model info to pac_expectation nodes
-  // Must be called after createPacGrowthNeutralityParameter() has been called
-  // on the relevant models
-  void fillPacModelInfo(const string &pac_model_name,
-                        vector<int> lhs,
-                        int max_lag,
-                        string aux_model_type,
-                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
-                        const vector<bool> &nonstationary,
-                        expr_t growth);
-
-  //! Substitutes pac_expectation operator with expectation based on auxiliary model
-  void substitutePacExpectation(const string &pac_model_name);
-
   //! Writes dynamic model file
   void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const;
   //! Writes file containing parameters derivatives
@@ -563,6 +513,8 @@ public:
   //! Substitute VarExpectation operators
   void substituteVarExpectation(const map<string, expr_t> &subst_table);
 
+  void analyzePacEquationStructure(const string &name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag, PacModelTable::equation_info_t &pac_equation_info);
+
   // Exception thrown by getPacTargetSymbId()
   class PacTargetNotIdentifiedException
   {
@@ -576,13 +528,42 @@ public:
   //! Return target of the pac equation
   int getPacTargetSymbId(const string &pac_model_name) const;
 
-  //! Declare Z1 variables before creating aux vars so it comes right after the endo names declared by the user
-  void declarePacModelConsistentExpectationEndogs(const string &name);
-
-  //! Add model consistent expectation equation for pac model
+  /* For a PAC MCE model, add the variable and the equation defining Z₁.
+     The symbol IDs of the new endogenous and parameters are also added to
+     pac_mce_{z1,alpha}_symb_ids. */
   void addPacModelConsistentExpectationEquation(const string &name, int discount,
                                                 const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
-                                                ExprNode::subst_table_t &diff_subst_table, expr_t growth);
+                                                ExprNode::subst_table_t &diff_subst_table,
+                                                map<pair<string, string>, int> &pac_mce_z1_symb_ids,
+                                                map<pair<string, string>, vector<int>> &pac_mce_alpha_symb_ids);
+
+  /* For a PAC MCE model, fill pac_expectation_substitution with the
+     expression that will be substituted for the pac_expectation operator */
+  void computePacModelConsistentExpectationSubstitution(const string &name,
+                                                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                                        expr_t growth,
+                                                        const map<string, int> &pac_growth_neutrality_params,
+                                                        const map<pair<string, string>, int> &pac_mce_z1_symb_ids,
+                                                        map<pair<string, string>, expr_t> &pac_expectation_substitution);
+
+
+  /* For a PAC backward model, fill pac_expectation_substitution with the
+     expression that will be substituted for the pac_expectation operator.
+     The symbol IDs of the new parameters are also added to pac_{h0,h1}_indices. */
+  void computePacBackwardExpectationSubstitution(const string &name,
+                                                 const vector<int> &lhs,
+                                                 int max_lag,
+                                                 const string &aux_model_type,
+                                                 const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                                 const vector<bool> &nonstationary,
+                                                 expr_t growth,
+                                                 const map<string, int> &pac_growth_neutrality_params,
+                                                 map<pair<string, string>, vector<int>> &pac_h0_indices,
+                                                 map<pair<string, string>, vector<int>> &pac_h1_indices,
+                                                 map<pair<string, string>, expr_t> &pac_expectation_substitution);
+
+  //! Substitutes pac_expectation operator with expectation based on auxiliary model
+  void substitutePacExpectation(const map<pair<string, string>, expr_t> &pac_expectation_substitution);
 
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 0866b4f2b19c3e10fcdb1282122424d18b23a93c..f6f9ee7b6b961e06756529625427103fc3b1012e 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -815,6 +815,9 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
                              DynamicModel &dynamic_model, const VarModelTable &var_model_table,
                              const TrendComponentModelTable &trend_component_model_table)
 {
+  // (model name, eq_name) → expression for pac_expectation
+  map<pair<string, string>, expr_t> pac_expectation_substitution;
+
   for (const auto &name : names)
     {
       /* Fill the growth_info structure.
@@ -837,32 +840,22 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
               }
           }
 
-      // 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
+      // Collect some information about PAC models
       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";
+          aux_model_type[name] = "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);
+          lhs[name] = 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";
+          aux_model_type[name] = "var";
           max_lag = var_model_table.getMaxLag(aux_model_name[name]);
-          lhs = var_model_table.getLhs(aux_model_name[name]);
+          lhs[name] = 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]);
         }
@@ -873,19 +866,49 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
           cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
           exit(EXIT_FAILURE);
         }
+      dynamic_model.analyzePacEquationStructure(name, eqtag_and_lag, 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);
+            }
+        }
 
-      auto eqtag_and_lag = dynamic_model.walkPacParameters(name);
+      // In the MCE case, add the variable and the equation defining Z₁
       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]);
+                                                               mce_z1_symb_ids, mce_alpha_symb_ids);
 
-      dynamic_model.substitutePacExpectation(name);
+      // Compute the expressions that will be substituted for the pac_expectation operators
+      if (aux_model_name[name].empty())
+        dynamic_model.computePacModelConsistentExpectationSubstitution(name, eqtag_and_lag,
+                                                                       growth[name],
+                                                                       growth_neutrality_params,
+                                                                       mce_z1_symb_ids,
+                                                                       pac_expectation_substitution);
+      else
+        dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag,
+                                                                aux_model_type[name],
+                                                                eqtag_and_lag, nonstationary,
+                                                                growth[name],
+                                                                growth_neutrality_params,
+                                                                h0_indices, h1_indices,
+                                                                pac_expectation_substitution);
     }
 
+  // Actually perform the substitution of pac_expectation
+  dynamic_model.substitutePacExpectation(pac_expectation_substitution);
   dynamic_model.checkNoRemainingPacExpectation();
 }
 
@@ -951,6 +974,237 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
             }
         }
     }
+
+  // Write PAC Model Consistent Expectation parameter info
+  for (auto &[key, ids] : mce_alpha_symb_ids)
+    {
+      output << "M_.pac." << key.first << ".equations." << key.second << ".mce.alpha = [";
+      for (auto id : ids)
+        output << symbol_table.getTypeSpecificID(id) + 1 << " ";
+      output << "];" << endl;
+    }
+
+  // Write PAC Model Consistent Expectation Z1 info
+  for (auto &[key, id] : mce_z1_symb_ids)
+    output << "M_.pac." << key.first << ".equations." << key.second << ".mce.z1 = "
+           << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
+
+  // Write PAC lag info
+  for (auto &[key, val] : eqtag_and_lag)
+    output << "M_.pac." << key.first << ".equations." << val.first << ".max_lag = " << val.second << ";" << endl;
+
+  // Write PAC equation tag info
+  map<string, vector<pair<string, string>>> for_writing;
+  for (auto &[key, val] : eqtag_and_lag)
+    for_writing[key.first].emplace_back(key.second, val.first);
+
+  for (auto &[key, val] : for_writing)
+    {
+      output << "M_.pac." << key << ".tag_map = [";
+      for (auto &[eqtag, standard_eqtag] : val)
+        output << "{'" << eqtag << "', '" << standard_eqtag << "'};";
+      output << "];" << 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, lhs] : lhs)
+    {
+      output << "M_.pac." << model << ".lhs = [";
+      for (auto id : lhs)
+        output << id + 1 << " ";
+      output << "];" << endl;
+    }
+
+  for (auto &[model, type] : aux_model_type)
+      output << "M_.pac." << model << ".auxiliary_model_type = '" << type << "';" << endl;
+
+  for (auto &[key, 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;
+      string substruct = key.first + ".equations." + key.second + ".";
+
+      output << "M_.pac." << substruct << "lhs_var = "
+             << symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl;
+
+      if (optim_share_index >= 0)
+        output << "M_.pac." << substruct << "share_of_optimizing_agents_index = "
+               << symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
+
+      output << "M_.pac." << substruct << "ec.params = "
+             << symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
+             << "M_.pac." << substruct << "ec.vars = [";
+      for (auto it : ec_params_and_vars.second)
+        output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ec.istarget = [";
+      for (auto it : ec_params_and_vars.second)
+        output << (get<1>(it) ? "true " : "false ");
+      output << "];" << endl
+             << "M_.pac." << substruct << "ec.scale = [";
+      for (auto it : ec_params_and_vars.second)
+        output << get<2>(it) << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "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." << substruct << "ar.params = [";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << (pid != -1 ? symbol_table.getTypeSpecificID(pid) + 1 : -1) << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ar.vars = [";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << (vid != -1 ? symbol_table.getTypeSpecificID(vid) + 1 : -1) << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ar.lags = [";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << vlag << " ";
+      output << "];" << endl;
+      if (!non_optim_vars_params_and_constants.empty())
+        {
+          output << "M_.pac." << substruct << "non_optimizing_behaviour.params = [";
+          for (auto &it : non_optim_vars_params_and_constants)
+            if (get<2>(it) >= 0)
+              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "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." << substruct << "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." << substruct << "non_optimizing_behaviour.lags = [";
+          for (auto &it : non_optim_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "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." << substruct << "additive.params = [";
+          for (auto &it : additive_vars_params_and_constants)
+            if (get<2>(it) >= 0)
+              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "additive.vars = [";
+          for (auto &it : additive_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "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." << substruct << "additive.lags = [";
+          for (auto &it : additive_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "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." << substruct << "optim_additive.params = [";
+          for (auto &it : optim_additive_vars_params_and_constants)
+            if (get<2>(it) >= 0)
+              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "optim_additive.vars = [";
+          for (auto &it : optim_additive_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "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." << substruct << "optim_additive.lags = [";
+          for (auto &it : optim_additive_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "optim_additive.scaling_factor = [";
+          for (auto &it : optim_additive_vars_params_and_constants)
+            output << get<3>(it) << " ";
+          output << "];" << endl;
+        }
+      // Create empty h0 and h1 substructures that will be overwritten later if not empty
+      output << "M_.pac." << substruct << "h0_param_indices = [];" << endl
+             << "M_.pac." << substruct << "h1_param_indices = [];" << endl;
+    }
+
+  for (auto &[model_eqtag, symb_ids] : h0_indices)
+    {
+      output << "M_.pac." << model_eqtag.first << ".equations." << model_eqtag.second << ".h0_param_indices = [";
+      for (auto it : symb_ids)
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl;
+    }
+
+  for (auto &[model_eqtag, symb_ids] : h1_indices)
+    {
+      output << "M_.pac." << model_eqtag.first << ".equations." << model_eqtag.second << ".h1_param_indices = [";
+      for (auto it : symb_ids)
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl;
+    }
 }
 
 void
diff --git a/src/SubModel.hh b/src/SubModel.hh
index bb1b798c01d8b59f2c6ad6be32f35bbe182bcc80..a4fe19063cb8a35a83c393eaccc37fa0b17ff990 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -200,6 +200,39 @@ private:
   // 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;
+
+  /* Stores symb_ids for alphas created by DynamicModel::addPacModelConsistentExpectationEquation()
+     (pac_model_name, standardized_eqtag) -> mce_alpha_symb_ids */
+  map<pair<string, string>, vector<int>> mce_alpha_symb_ids;
+  /* Stores symb_ids for z1s created by DynamicModel::addPacModelConsistentExpectationEquation()
+     (pac_model_name, standardized_eqtag) -> mce_z1_symb_id */
+  map<pair<string, string>, int> mce_z1_symb_ids;
+  /* Stores symb_ids for h0, h1 parameters
+     (pac_model_name, standardized_eqtag) -> parameter symb_ids */
+  map<pair<string, string>, vector<int>> h0_indices, h1_indices;
+  /* Stores indices for growth neutrality parameters
+     pac_model_name -> growth_neutrality_param_index */
+  map<string, int> growth_neutrality_params;
+
+  // Stores LHS vars (only for backward PAC models)
+  map<string, vector<int>> lhs;
+
+  // Stores auxiliary model type (only for backward PAC models)
+  map<string, string> aux_model_type;
+
+  /* Stores lag info for pac equations
+     (pac_model_name, equation_tag) -> (standardized_eqtag, lag) */
+  map<pair<string, string>, pair<string, int>> eqtag_and_lag;
+public:
+  /* Stores info about PAC models specific to the equation they appear in
+     (pac_model_name, standardized_eqtag) ->
+         (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants)
+  */
+  using equation_info_t = map<pair<string, string>,
+                              tuple<pair<int, int>, int, vector<tuple<int, int, int>>, pair<int, vector<tuple<int, bool, int>>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>>>;
+private:
+  equation_info_t equation_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);