diff --git a/DynamicModel.cc b/DynamicModel.cc
index fb23c71d0182026764b8194684d6ee7d7f371f84..dc79dcb888e66c6a6026f961ea775ac34c6e0704 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3344,7 +3344,8 @@ DynamicModel::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<
 void
 DynamicModel::substitutePacExpectation()
 {
-  map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > subst_table;
+  // maps PacExpectationNode to (subst node, (var_model_name, (growth_id, (h0_idxs, h1_idxs))))
+  map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > subst_table;
   for (map<int, expr_t>::iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     it->second = it->second->substitutePacExpectation(subst_table);
@@ -3356,7 +3357,7 @@ DynamicModel::substitutePacExpectation()
       equations[i] = substeq;
     }
 
-  for (map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::const_iterator it = subst_table.begin(); it != subst_table.end(); it++)
+  for (map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::const_iterator it = subst_table.begin(); it != subst_table.end(); it++)
     pac_expectation_info[it->second.second.first] = it->second.second.second;
 }
 
diff --git a/DynamicModel.hh b/DynamicModel.hh
index c2a3fa6466c4da7650032aea27ec42cd99a126ce..5edd6b1ddf0981e742639f8f58c28eb7a5520d1e 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -227,7 +227,7 @@ private:
   map<string, set<int> > var_expectation_functions_to_write;
 
   //! Used for pac_expectation operator
-  // maps model_name to (growth_idx, (h0_indices, h1_indices))
+  // maps var_model_name to (growth_idx, (h0_indices, h1_indices))
   map<string, pair<int, pair<vector<int>, vector<int> > > > pac_expectation_info;
 
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
diff --git a/ExprNode.cc b/ExprNode.cc
index 42254c75e2fd42fbfb7f3ff9d77cab3773fc443d..6d8a4b13ef75391452243909e12f42077fb01133 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -494,7 +494,7 @@ NumConstNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
 }
 
 expr_t
-NumConstNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+NumConstNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1273,7 +1273,7 @@ VariableNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
 }
 
 expr_t
-VariableNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+VariableNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   return const_cast<VariableNode *>(this);
 }
@@ -2821,7 +2821,7 @@ UnaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &
 }
 
 expr_t
-UnaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+UnaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   expr_t argsubst = arg->substitutePacExpectation(subst_table);
   return buildSimilarUnaryOpNode(argsubst, datatree);
@@ -4445,7 +4445,7 @@ BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
 }
 
 expr_t
-BinaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+BinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
   expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
@@ -5217,7 +5217,7 @@ TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
 }
 
 expr_t
-TrinaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+TrinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
   expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
@@ -5555,7 +5555,7 @@ AbstractExternalFunctionNode::substituteDiff(subst_table_t &subst_table, vector<
 }
 
 expr_t
-AbstractExternalFunctionNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+AbstractExternalFunctionNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
@@ -7046,7 +7046,7 @@ VarExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNo
 }
 
 expr_t
-VarExpectationNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -7509,9 +7509,9 @@ PacExpectationNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>,
 }
 
 expr_t
-PacExpectationNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+PacExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
 {
-  map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::iterator myit =
+  map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::const_iterator myit =
     subst_table.find(const_cast<PacExpectationNode *>(this));
   if (myit != subst_table.end())
     return const_cast<BinaryOpNode *>(myit->second.first);
diff --git a/ExprNode.hh b/ExprNode.hh
index 6ec139454c1663b1251057dedfc62b889a7c97c4..4f1d72972594a183e8d8e1af7ab2f1baf6acd455 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -35,6 +35,7 @@ using namespace std;
 class DataTree;
 class VariableNode;
 class BinaryOpNode;
+class PacExpectationNode;
 
 typedef class ExprNode *expr_t;
 
@@ -123,7 +124,7 @@ enum ExprNodeOutputType
 #define MIN_COST(is_matlab) ((is_matlab) ? MIN_COST_MATLAB : MIN_COST_C)
 
 //! Base class for expression nodes
-    class ExprNode
+class ExprNode
     {
       friend class DataTree;
       friend class DynamicModel;
@@ -465,7 +466,7 @@ enum ExprNodeOutputType
       virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
-      virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table) = 0;
+      virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table) = 0;
 
       //! Add ExprNodes to the provided datatree
       virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
@@ -542,7 +543,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -620,7 +621,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -721,7 +722,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -834,7 +835,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -927,7 +928,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -1020,7 +1021,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1204,7 +1205,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1270,7 +1271,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,