diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index f6fd184ee678c9385d53abe061ae85cf96dce7ec..1bc39dc199037007be8d3fac14f9494f757d0be9 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3783,11 +3783,14 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
 }
 
 void
-DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int discount_symb_id,
-                                                       int pac_eq_max_lag,
-                                                       ExprNode::subst_table_t &diff_subst_table,
-                                                       map<string, int> &pac_mce_z1_symb_ids,
-                                                       map<string, vector<int>> &pac_mce_alpha_symb_ids)
+DynamicModel::computePacModelConsistentExpectationSubstitution(const string &name,
+                                                               int discount_symb_id,
+                                                               int pac_eq_max_lag,
+                                                               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, expr_t> &pac_expectation_substitution)
 {
   int pac_target_symb_id;
   try
@@ -3813,7 +3816,7 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
       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);
     }
-  pac_mce_z1_symb_ids[name] = mce_z1_symb_id;
+  pac_aux_var_symb_ids[name] = mce_z1_symb_id;
 
   expr_t A = One;
   expr_t fp = Zero;
@@ -3913,18 +3916,11 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
   neqs++;
 
   cout << "PAC Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl;
-}
 
-void
-DynamicModel::computePacModelConsistentExpectationSubstitution(const string &name,
-                                                               expr_t growth_correction_term,
-                                                               int pac_mce_z1_symb_id,
-                                                               map<string, expr_t> &pac_expectation_substitution)
-{
   /* 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. */
-  pac_expectation_substitution[name] = AddPlus(AddVariable(pac_mce_z1_symb_id), growth_correction_term);
+  pac_expectation_substitution[name] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction_term);
 }
 
 void
@@ -3934,6 +3930,7 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
                                                         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, expr_t> &pac_expectation_substitution)
@@ -4001,7 +3998,21 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
           }
     }
 
-  pac_expectation_substitution[name] = AddPlus(subExpr, growth_correction_term);
+  int expect_var_id;
+  string expect_var_name = "pac_expectation_" + name;
+  try
+    {
+      expect_var_id = symbol_table.addSymbol(expect_var_name, SymbolType::endogenous);
+    }
+  catch (SymbolTable::AlreadyDeclaredException &e)
+    {
+      cerr << "The variable/parameter '" << expect_var_name << "' conflicts with the variable that will be generated for the 'pac_expectation' expression of the '" << name << "' PAC model. Please rename it." << endl;
+      exit(EXIT_FAILURE);
+    }
+  addEquation(AddEqual(AddVariable(expect_var_id), AddPlus(subExpr, growth_correction_term)), -1);
+
+  pac_aux_var_symb_ids[name] = expect_var_id;
+  pac_expectation_substitution[name] = AddVariable(expect_var_id);
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 387f0e8b989bd07a72b3acbc0d15e20f940af787..80579ceb65ad979cd28336c8592f7d6ca7be0367 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -528,20 +528,18 @@ public:
   //! Return target of the pac equation
   int getPacTargetSymbId(const string &pac_model_name) const;
 
-  /* 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,
-                                                int pac_eq_max_lag,
-                                                ExprNode::subst_table_t &diff_subst_table,
-                                                map<string, int> &pac_mce_z1_symb_ids,
-                                                map<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 */
+     expression that will be substituted for the pac_expectation operator.
+     In the process, add the variable and the equation defining Z₁.
+     The symbol IDs of the new endogenous are added to pac_aux_var_symb_ids,
+     and the new auxiliary parameters to pac_mce_alpha_symb_ids.
+  */
   void computePacModelConsistentExpectationSubstitution(const string &name,
+                                                        int discount_symb_id, int pac_eq_max_lag,
                                                         expr_t growth_correction_term,
-                                                        int pac_mce_z1_symb_id,
+                                                        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, expr_t> &pac_expectation_substitution);
 
 
@@ -554,6 +552,7 @@ public:
                                                  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, expr_t> &pac_expectation_substitution);
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 3c9114376680b0272be586225707fa6b89f28a4c..f484b2afb973d8aaab4f09fe23a8fa27d9b4eca0 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -884,27 +884,25 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
             }
         }
 
-      // 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]),
-                                                               pacEquationMaxLag(name),
-                                                               diff_subst_table,
-                                                               mce_z1_symb_ids, mce_alpha_symb_ids);
-
       // Compute the expressions that will be substituted for the pac_expectation operators
       expr_t growth_correction_term = dynamic_model.Zero;
       if (growth[name])
         growth_correction_term = dynamic_model.AddTimes(growth[name], dynamic_model.AddVariable(growth_neutrality_params[name]));
       if (aux_model_name[name].empty())
         dynamic_model.computePacModelConsistentExpectationSubstitution(name,
+                                                                       symbol_table.getID(discount[name]),
+                                                                       pacEquationMaxLag(name),
                                                                        growth_correction_term,
-                                                                       mce_z1_symb_ids[name],
+                                                                       diff_subst_table,
+                                                                       aux_var_symb_ids,
+                                                                       mce_alpha_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,
                                                                 pac_expectation_substitution);
     }
@@ -986,10 +984,10 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const
       output << "];" << endl;
     }
 
-  // Write PAC Model Consistent Expectation Z1 info
-  for (auto &[name, id] : mce_z1_symb_ids)
-    output << "M_.pac." << name << ".mce.z1 = "
-           << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
+  // Write the auxiliary variable IDs created for the pac_expectation operator
+  for (auto &[name, id] : aux_var_symb_ids)
+    output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z1" : "aux_id")
+           << " = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
 
   // Write PAC equation name info
   for (auto &[name, eq] : eq_name)
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 009977968d59c4b660f6cdd7d1d66eefeb468991..ae84b6855bd008ea01f433fd650b9e864729fa94 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -205,12 +205,17 @@ private:
      pac_model_name → eq_name */
   map<string, string> eq_name;
 
+  /* Stores symb_ids for auxiliary endogenous created for the expression
+     substituted to the pac_expectation operator:
+     - in the backward case, this auxiliary contains exactly the
+     pac_expectation value
+     - in the MCE case, this auxiliary represents Z₁ (i.e. without the growth
+     correction term)
+     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 z1s created by DynamicModel::addPacModelConsistentExpectationEquation()
-     pac_model_name → mce_z1_symb_id */
-  map<string, int> mce_z1_symb_ids;
   /* Stores symb_ids for h0, h1 parameters
      pac_model_name → parameter symb_ids */
   map<string, vector<int>> h0_indices, h1_indices;