diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 77c8d60acd43e946eb77f486f507aa218133b021..9dfec6b170f9f21a44e1feedd63eb8197f2a9479 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3804,18 +3804,8 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(const string &nam
     }
   int neqs = 0;
 
-  // Create the endogenous representing Z₁
-  int mce_z1_symb_id;
-  string varname = "mce_Z1_" + name;
-  try
-    {
-      mce_z1_symb_id = symbol_table.addSymbol(varname, SymbolType::endogenous);
-    }
-  catch (SymbolTable::AlreadyDeclaredException &e)
-    {
-      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);
-    }
+  // Create the endogenous representing Z₁ (no orig_expr is given since its definition is recursive)
+  int mce_z1_symb_id = symbol_table.addPacExpectationAuxiliaryVar("mce_Z1_" + name, nullptr);
   pac_aux_var_symb_ids[name] = mce_z1_symb_id;
 
   expr_t A = One;
@@ -3913,6 +3903,8 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(const string &nam
   auto neweq = AddEqual(AddVariable(mce_z1_symb_id),
                         AddMinus(AddTimes(A, AddMinus(const_cast<VariableNode *>(target_base_diff_node), fs)), fp));
   addEquation(neweq, -1);
+  /* This equation is not added to the list of auxiliary equations, because it
+     is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
   neqs++;
 
   cout << "PAC Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl;
@@ -3967,19 +3959,12 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name,
                                    AddVariable(lhsit, -i)));
       }
 
-  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);
+  subExpr = AddPlus(subExpr, growth_correction_term);
 
+  int expect_var_id = symbol_table.addPacExpectationAuxiliaryVar("pac_expectation_" + name, subExpr);
+  expr_t neweq = AddEqual(AddVariable(expect_var_id), subExpr);
+  addEquation(neweq, -1);
+  addAuxEquation(neweq);
   pac_aux_var_symb_ids[name] = expect_var_id;
   pac_expectation_substitution[name] = AddVariable(expect_var_id);
 }
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 76ad7718e712fba65bdec008616f5ea3762ace98..23c028ea585f9e979c63a9fc7a34b43457e8cbb8 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -376,6 +376,7 @@ SymbolTable::writeOutput(ostream &output) const noexcept(false)
             output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl;
             break;
           case AuxVarType::expectation:
+          case AuxVarType::pacExpectation:
             break;
           case AuxVarType::diff:
           case AuxVarType::diffLag:
@@ -668,6 +669,24 @@ SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id, expr_t expr_arg) noexc
   return symb_id;
 }
 
+int
+SymbolTable::addPacExpectationAuxiliaryVar(const string &name, expr_t expr_arg)
+{
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(name, SymbolType::endogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: the variable/parameter '" << name << "' conflicts with a variable that will be generated for a 'pac_expectation' expression. Please rename it." << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.emplace_back(symb_id, AuxVarType::pacExpectation, 0, 0, 0, 0, expr_arg, "");
+  return symb_id;
+}
+
 int
 SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const noexcept(false)
 {
@@ -984,6 +1003,7 @@ SymbolTable::writeJsonOutput(ostream &output) const
 	      output << R"(, orig_index": )" << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1;
 	      break;
 	    case AuxVarType::expectation:
+	    case AuxVarType::pacExpectation:
 	      break;
 	    case AuxVarType::diff:
 	    case AuxVarType::diffLag:
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index c63719eed0c9710fd7e158a788d80383c02718a9..b3fd7be48db353386486e8e1bc5487669fd62c0e 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -52,7 +52,8 @@ enum class AuxVarType
    diff = 8, //!< Variable for Diff operator
    diffLag = 9, //!< Variable for timing between Diff operators (lag)
    unaryOp = 10, //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
-   diffLead = 11 //!< Variable for timing between Diff operators (lead)
+   diffLead = 11, //!< Variable for timing between Diff operators (lead)
+   pacExpectation = 12 //!< Variable created for the substitution of the pac_expectation operator
   };
 
 //! Information on some auxiliary variables
@@ -334,6 +335,8 @@ public:
   int addDiffLeadAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lead) noexcept(false);
   //! An Auxiliary variable for a unary op
   int addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op, int orig_symb_id = -1, int orig_lag = 0) noexcept(false);
+  //! An auxiliary variable for a pac_expectation operator
+  int addPacExpectationAuxiliaryVar(const string &name, expr_t expr_arg);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const