diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 1bc39dc199037007be8d3fac14f9494f757d0be9..77c8d60acd43e946eb77f486f507aa218133b021 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3789,7 +3789,7 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(const string &nam
                                                                expr_t growth_correction_term,
                                                                ExprNode::subst_table_t &diff_subst_table,
                                                                map<string, int> &pac_aux_var_symb_ids,
-                                                               map<string, vector<int>> &pac_mce_alpha_symb_ids,
+                                                               map<string, vector<int>> &pac_aux_param_symb_ids,
                                                                map<string, expr_t> &pac_expectation_substitution)
 {
   int pac_target_symb_id;
@@ -3827,7 +3827,7 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(const string &nam
       try
         {
           int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
-          pac_mce_alpha_symb_ids[name].push_back(alpha_i_symb_id);
+          pac_aux_param_symb_ids[name].push_back(alpha_i_symb_id);
           A = AddPlus(A, AddVariable(alpha_i_symb_id));
           fp = AddPlus(fp,
                        AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
@@ -3928,16 +3928,11 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
                                                         const vector<int> &lhs,
                                                         int max_lag,
                                                         const string &aux_model_type,
-                                                        const vector<bool> &nonstationary,
                                                         expr_t growth_correction_term,
                                                         map<string, int> &pac_aux_var_symb_ids,
-                                                        map<string, vector<int>> &pac_h0_indices,
-                                                        map<string, vector<int>> &pac_h1_indices,
+                                                        map<string, vector<int>> &pac_aux_param_symb_ids,
                                                         map<string, expr_t> &pac_expectation_substitution)
 {
-  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
-
   auto create_aux_param = [&](const string &param_name)
   {
     try
@@ -3952,51 +3947,25 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
   };
 
   expr_t subExpr = Zero;
-  if (stationary_vars_present)
-    {
-      if (aux_model_type == "var")
-        {
-          /* If the auxiliary model is a VAR, add a parameter corresponding
-             to the constant. */
-          int new_param_symb_id = create_aux_param("h0_" + name + "_constant");
-          pac_h0_indices[name].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)
-          {
-            int new_param_symb_id = create_aux_param("h0_" + name + "_var_"
-                                                     + symbol_table.getName(lhsit)
-                                                     + "_lag_" + to_string(i));
-            pac_h0_indices[name].push_back(new_param_symb_id);
-            subExpr = AddPlus(subExpr,
-                              AddTimes(AddVariable(new_param_symb_id),
-                                       AddVariable(lhsit, -i)));
-          }
-    }
-
-  if (nonstationary_vars_present)
+  if (aux_model_type == "var")
     {
-      if (aux_model_type == "var")
-        {
-          /* If the auxiliary model is a VAR, add a parameter corresponding
-             to the constant. */
-          int new_param_symb_id = create_aux_param("h1_" + name + "_constant");
-          pac_h1_indices[name].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)
-          {
-            int new_param_symb_id = create_aux_param("h1_" + name + "_var_"
-                                                     + symbol_table.getName(lhsit)
-                                                     + "_lag_" + to_string(i));
-            pac_h1_indices[name].push_back(new_param_symb_id);
-            subExpr = AddPlus(subExpr,
-                              AddTimes(AddVariable(new_param_symb_id),
-                                       AddVariable(lhsit, -i)));
-          }
+      /* If the auxiliary model is a VAR, add a parameter corresponding
+         to the constant. */
+      int new_param_symb_id = create_aux_param("h_" + name + "_constant");
+      pac_aux_param_symb_ids[name].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)
+      {
+        int new_param_symb_id = create_aux_param("h_" + name + "_var_"
+                                                 + symbol_table.getName(lhsit)
+                                                 + "_lag_" + to_string(i));
+        pac_aux_param_symb_ids[name].push_back(new_param_symb_id);
+        subExpr = AddPlus(subExpr,
+                          AddTimes(AddVariable(new_param_symb_id),
+                                   AddVariable(lhsit, -i)));
+      }
 
   int expect_var_id;
   string expect_var_name = "pac_expectation_" + name;
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 80579ceb65ad979cd28336c8592f7d6ca7be0367..7f0044752d3daa60c43e5cd3f2d8e016edb10414 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -539,7 +539,7 @@ public:
                                                         expr_t growth_correction_term,
                                                         ExprNode::subst_table_t &diff_subst_table,
                                                         map<string, int> &pac_aux_var_symb_ids,
-                                                        map<string, vector<int>> &pac_mce_alpha_symb_ids,
+                                                        map<string, vector<int>> &pac_aux_param_symb_ids,
                                                         map<string, expr_t> &pac_expectation_substitution);
 
 
@@ -550,11 +550,9 @@ public:
                                                  const vector<int> &lhs,
                                                  int max_lag,
                                                  const string &aux_model_type,
-                                                 const vector<bool> &nonstationary,
                                                  expr_t growth_correction_term,
                                                  map<string, int> &pac_aux_var_symb_ids,
-                                                 map<string, vector<int>> &pac_h0_indices,
-                                                 map<string, vector<int>> &pac_h1_indices,
+                                                 map<string, vector<int>> &pac_aux_param_symb_ids,
                                                  map<string, expr_t> &pac_expectation_substitution);
 
   //! Substitutes pac_expectation operator with expectation based on auxiliary model
diff --git a/src/SubModel.cc b/src/SubModel.cc
index f484b2afb973d8aaab4f09fe23a8fa27d9b4eca0..dbfaf0e1cc118f95109528efbd1a3b0b523e2eb7 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -842,22 +842,17 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
 
       // Collect some information about PAC models
       int max_lag;
-      vector<bool> nonstationary;
       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);
-          // 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[name] = "var";
           max_lag = var_model_table.getMaxLag(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]);
         }
       else if (aux_model_name[name].empty())
         max_lag = 0;
@@ -895,15 +890,14 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
                                                                        growth_correction_term,
                                                                        diff_subst_table,
                                                                        aux_var_symb_ids,
-                                                                       mce_alpha_symb_ids,
+                                                                       aux_param_symb_ids,
                                                                        pac_expectation_substitution);
       else
         dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag,
                                                                 aux_model_type[name],
-                                                                nonstationary,
                                                                 growth_correction_term,
                                                                 aux_var_symb_ids,
-                                                                h0_indices, h1_indices,
+                                                                aux_param_symb_ids,
                                                                 pac_expectation_substitution);
     }
 
@@ -975,10 +969,10 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
         }
     }
 
-  // Write PAC Model Consistent Expectation parameter info
-  for (auto &[name, ids] : mce_alpha_symb_ids)
+  // Write the auxiliary parameter IDs created for the pac_expectation operator
+  for (auto &[name, ids] : aux_param_symb_ids)
     {
-      output << "M_.pac." << name << ".mce.alpha = [";
+      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;
@@ -1171,25 +1165,6 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
             output << get<3>(it) << " ";
           output << "];" << endl;
         }
-      // Create empty h0 and h1 substructures that will be overwritten later if not empty
-      output << "M_.pac." << name << ".h0_param_indices = [];" << endl
-             << "M_.pac." << name << ".h1_param_indices = [];" << endl;
-    }
-
-  for (auto &[name, symb_ids] : h0_indices)
-    {
-      output << "M_.pac." << name << ".h0_param_indices = [";
-      for (auto it : symb_ids)
-        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
-      output << "];" << endl;
-    }
-
-  for (auto &[name, symb_ids] : h1_indices)
-    {
-      output << "M_.pac." << name << ".h1_param_indices = [";
-      for (auto it : symb_ids)
-        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
-      output << "];" << endl;
     }
 }
 
diff --git a/src/SubModel.hh b/src/SubModel.hh
index ae84b6855bd008ea01f433fd650b9e864729fa94..f3a7921e08d7a92c37c7e2ae7d990c883ee8eb1a 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -213,12 +213,13 @@ private:
      correction term)
      pac_model_name → symb_id */
   map<string, int> aux_var_symb_ids;
-  /* Stores symb_ids for alphas created by DynamicModel::addPacModelConsistentExpectationEquation()
-     pac_model_name → mce_alpha_symb_ids */
-  map<string, vector<int>> mce_alpha_symb_ids;
-  /* Stores symb_ids for h0, h1 parameters
-     pac_model_name → parameter symb_ids */
-  map<string, vector<int>> h0_indices, h1_indices;
+  /* Stores symb_ids for auxiliary parameters created for the expression
+     substituted to the pac_expectation operator (excluding the growth
+     neutrality correction):
+     - in the backward case, contains the “h” parameters
+     - in the MCE case, contains the “α” parameters
+     pac_model_name → symb_ids */
+  map<string, vector<int>> aux_param_symb_ids;
   /* Stores indices for growth neutrality parameters
      pac_model_name → growth_neutrality_param_index */
   map<string, int> growth_neutrality_params;