diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index f5f2f8e9aea9be44b0b7d70b2cd73cc7a77bd62a..44bb929c38179ac31fae3d503fcdcb02ca4780b6 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3618,12 +3618,43 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   // Write Pac Model Consistent Expectation parameter info
   for (auto & it : pac_mce_alpha_symb_ids)
     {
-      output << modstruct << "pac." << it.first << ".mce_alpha_idxs = [";
+      output << modstruct << "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 << modstruct << "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 << modstruct << "pac." << it.first.first << ".equations." << it.second.first << ".max_lag = " << it.second.second << ";" << endl;
+
+  map <string, vector<pair<string, string>> > mymap;
+  mymap.insert(pair<string,vector<pair<string,string>> >("A", vector<pair<string,string>>()));
+  mymap["A"].push_back(make_pair("A","BB"));
+  mymap["A"].push_back(make_pair("C", "DD"));
+
+  // Write Pac equation tag info
+  map<string, vector<pair<string, string>>> for_writing;
+  for (auto & it : pac_eqtag_and_lag)
+    {
+      if (for_writing.find(it.first.first) == for_writing.end())
+        for_writing.insert(pair<string,vector<pair<string,string>> >(it.first.first, vector<pair<string,string>>()));
+      for_writing[it.first.first].push_back(make_pair(it.first.second, it.second.first));
+    }
+
+  for (auto & it : for_writing)
+    {
+      output << modstruct << "pac." << it.first << ".tag_map = [";
+      for (auto & it1 : it.second)
+        output << "{'" << it1.first << "', '" << it1.second << "'};" ;
+      output << "];" << endl;
+    }
+
   // Write PacExpectationInfo
   for (auto it : pac_expectation_info)
     it->ExprNode::writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
@@ -4250,8 +4281,9 @@ DynamicModel::getUndiffLHSForPac(const string &aux_model_name,
 }
 
 void
-DynamicModel::walkPacParameters()
+DynamicModel::walkPacParameters(const string &name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag)
 {
+  int i = 0;
   for (auto & equation : equations)
     {
       pair<int, int> lhs (-1, -1);
@@ -4296,16 +4328,33 @@ DynamicModel::walkPacParameters()
                   exit(EXIT_FAILURE);
                 }
             }
-          equation->addParamInfoToPac(lhs,
+
+          string eqtag = "";
+          for (auto & tag : equation_tags)
+            if (tag.first == (&equation - &equations[0]))
+              if (tag.second.first == "name")
+                {
+                  eqtag = tag.second.second;
+                  break;
+                }
+          if (eqtag == "")
+            {
+              cerr << "Every equation with a pac expectation must have been assigned an equation tag name" << endl;
+              exit(EXIT_FAILURE);
+            }
+          string eq = "eq" + to_string(i++);
+          equation->addParamInfoToPac(eq,
+                                      lhs,
                                       optim_share_index,
                                       ec_params_and_vars, ar_params_and_vars,
                                       non_optim_vars_params_and_constants);
+          eqtag_and_lag[make_pair(name, eqtag)] = make_pair(eq, 0);
         }
     }
 }
 
-int
-DynamicModel::getPacMaxLag(const string &pac_model_name) const
+void
+DynamicModel::getPacMaxLag(const string &pac_model_name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag) const
 {
   for (auto & equation : equations)
     if (equation->containsPacExpectation(pac_model_name))
@@ -4318,9 +4367,18 @@ DynamicModel::getPacMaxLag(const string &pac_model_name) const
                  << endl;
             exit(EXIT_FAILURE);
           }
-        return equation->PacMaxLag((*(endogs.begin())).first);
+
+        string eqtag = "";
+        for (auto & tag : equation_tags)
+          if (tag.first == (&equation - &equations[0]))
+            if (tag.second.first == "name")
+              {
+                eqtag = tag.second.second;
+                break;
+              }
+        string eq = eqtag_and_lag[make_pair(pac_model_name, eqtag)].first;
+        eqtag_and_lag[make_pair(pac_model_name, eqtag)] = make_pair(eq, equation->PacMaxLag((*(endogs.begin())).first));
       }
-  return 0;
 }
 
 int
@@ -4346,133 +4404,122 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
   return -1;
 }
 
-int
+void
 DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id,
-                                                       int discount_symb_id, int pac_max_lag_m,
+                                                       int discount_symb_id, const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
                                                        ExprNode::subst_table_t &diff_subst_table)
 {
-  int mce_symb_id;
-  try
-    {
-      mce_symb_id = symbol_table.addSymbol("mce_Z_" + name, SymbolType::endogenous);
-    }
-  catch (SymbolTable::AlreadyDeclaredException &e)
+  pac_eqtag_and_lag.insert(eqtag_and_lag.begin(), eqtag_and_lag.end());
+  int neqs = 0;
+  for (auto & it : eqtag_and_lag)
     {
-      cerr << "Variable name needed by PAC (mce_Z_" << name << endl;
-      exit(EXIT_FAILURE);
-    }
+      string eqtag = it.second.first;
+      int pac_max_lag_m = it.second.second;
+      string append_to_name = name + "_" + eqtag;
+      int mce_z1_symb_id;
+      try
+        {
+          mce_z1_symb_id = symbol_table.addSymbol("mce_Z1_" + append_to_name, SymbolType::endogenous);
+        }
+      catch (SymbolTable::AlreadyDeclaredException &e)
+        {
+          cerr << "Variable name needed by PAC (mce_Z1_" << append_to_name << endl;
+          exit(EXIT_FAILURE);
+        }
 
-  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_" + name + "_" + to_string(i),
-                                                     SymbolType::parameter);
-        pac_mce_alpha_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),
-                                       AddPower(beta, AddPossiblyNegativeConstant(i))),
-                              AddVariable(mce_symb_id, i)));
+      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[make_pair(name, 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_" << name << "_" << i << ")" << endl;
-        exit(EXIT_FAILURE);
-      }
+          }
+        catch (SymbolTable::AlreadyDeclaredException &e)
+          {
+            cerr << "Variable name needed by PAC (mce_alpha_" << append_to_name << "_" << i << ")" << endl;
+            exit(EXIT_FAILURE);
+          }
 
-  // Add diff nodes and eqs for pac_target_symb_id
-  int neqs = 0;
-  const VariableNode *target_base_diff_node;
-  expr_t diff_node_to_search = AddDiff(AddVariable(pac_target_symb_id));
-  ExprNode::subst_table_t::const_iterator sit = diff_subst_table.find(diff_node_to_search);
-  if (sit != diff_subst_table.end())
-    target_base_diff_node = sit->second;
-  else
-    {
-      int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
-      target_base_diff_node = AddVariable(symb_id);
-      addEquation(dynamic_cast<BinaryOpNode *>(AddEqual((expr_t) target_base_diff_node,
-                                                        AddMinus(AddVariable(pac_target_symb_id),
-                                                                 AddVariable(pac_target_symb_id, -1)))), -1);
-      neqs++;
-    }
+      // Add diff nodes and eqs for pac_target_symb_id
+      const VariableNode *target_base_diff_node;
+      expr_t diff_node_to_search = AddDiff(AddVariable(pac_target_symb_id));
+      ExprNode::subst_table_t::const_iterator sit = diff_subst_table.find(diff_node_to_search);
+      if (sit != diff_subst_table.end())
+        target_base_diff_node = sit->second;
+      else
+        {
+          int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
+          target_base_diff_node = AddVariable(symb_id);
+          addEquation(dynamic_cast<BinaryOpNode *>(AddEqual((expr_t) target_base_diff_node,
+                                                            AddMinus(AddVariable(pac_target_symb_id),
+                                                                     AddVariable(pac_target_symb_id, -1)))), -1);
+          neqs++;
+        }
 
-  map<int, VariableNode *> target_aux_var_to_add;
-  const VariableNode *last_aux_var = target_base_diff_node;
-  for (int i = 1; i <= pac_max_lag_m - 1; i++, neqs++)
-    {
-      expr_t this_diff_node = AddDiff(AddVariable(pac_target_symb_id, i));
-      int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node,
-                                                         last_aux_var->symb_id, last_aux_var->lag);
-      VariableNode *current_aux_var = AddVariable(symb_id);
-      addEquation(dynamic_cast<BinaryOpNode *>(AddEqual(current_aux_var,
-                                                        AddVariable(last_aux_var->symb_id, 1))), -1);
-      last_aux_var = current_aux_var;
-      target_aux_var_to_add[i] = current_aux_var;
-    }
+      map<int, VariableNode *> target_aux_var_to_add;
+      const VariableNode *last_aux_var = target_base_diff_node;
+      for (int i = 1; i <= pac_max_lag_m - 1; i++, neqs++)
+        {
+          expr_t this_diff_node = AddDiff(AddVariable(pac_target_symb_id, i));
+          int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node,
+                                                             last_aux_var->symb_id, last_aux_var->lag);
+          VariableNode *current_aux_var = AddVariable(symb_id);
+          addEquation(dynamic_cast<BinaryOpNode *>(AddEqual(current_aux_var,
+                                                            AddVariable(last_aux_var->symb_id, 1))), -1);
+          last_aux_var = current_aux_var;
+          target_aux_var_to_add[i] = current_aux_var;
+        }
 
-  expr_t fs = Zero;
-  for (int k = 1; k <= pac_max_lag_m - 1; k++)
-    {
-      expr_t ssum = Zero;
-      for (int j = k+1; j <= pac_max_lag_m; j++)
+      expr_t fs = Zero;
+      for (int k = 1; k <= pac_max_lag_m - 1; k++)
         {
-          int alpha_j_symb_id = -1;
-          string varname = "mce_alpha_" + name + "_" + to_string(j);
-          try
+          expr_t ssum = Zero;
+          for (int j = k+1; j <= pac_max_lag_m; j++)
             {
-              alpha_j_symb_id = symbol_table.getID(varname);
-            }
-          catch (SymbolTable::UnknownSymbolNameException &e)
-            {
-              alpha_j_symb_id = symbol_table.addSymbol(varname, SymbolType::parameter);
+              int alpha_j_symb_id = -1;
+              string varname = "mce_alpha_" + append_to_name + "_" + to_string(j);
+              try
+                {
+                  alpha_j_symb_id = symbol_table.getID(varname);
+                }
+              catch (SymbolTable::UnknownSymbolNameException &e)
+                {
+                  alpha_j_symb_id = symbol_table.addSymbol(varname, SymbolType::parameter);
+                }
+              ssum = AddPlus(ssum,
+                             AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j))));
             }
-          ssum = AddPlus(ssum,
-                         AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j))));
+          fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k]));
         }
-      fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k]));
+      addEquation(AddEqual(AddVariable(mce_z1_symb_id),
+                           AddMinus(AddTimes(A, AddMinus((expr_t) target_base_diff_node, fs)), fp)), -1);
+      neqs++;
+      pac_mce_z1_symb_ids[make_pair(name, eqtag)] = mce_z1_symb_id;
     }
-  addEquation(AddEqual(AddVariable(mce_symb_id),
-                       AddMinus(AddTimes(A, AddMinus((expr_t) target_base_diff_node, fs)), fp)), -1);
-  neqs++;
-  cout << "Pac Model Consistent Expectation: added " << neqs << " auxiliary variables and equations." << endl;
-  return mce_symb_id;
+    cout << "Pac Model Consistent Expectation: added " << neqs << " auxiliary variables and equations." << endl;
 }
 
 void
 DynamicModel::fillPacExpectationVarInfo(const string &pac_model_name,
                                         vector<int> &lhs,
                                         int max_lag,
-                                        int pac_max_lag,
-                                        vector<bool> &nonstationary,
+                                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                        const vector<bool> &nonstationary,
                                         int growth_symb_id, int growth_lag)
 {
+  pac_eqtag_and_lag.insert(eqtag_and_lag.begin(), eqtag_and_lag.end());
   for (size_t i = 0; i < equations.size(); i++)
     equations[i]->fillPacExpectationVarInfo(pac_model_name, lhs, max_lag,
-                                            pac_max_lag, nonstationary, growth_symb_id, growth_lag, i);
-}
-
-void
-DynamicModel::substitutePacExpectation(const string & name, int model_consistent_expectation_symb_id)
-{
-  map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
-  for (auto & it : local_variables_table)
-    it.second = it.second->substitutePacExpectation(name, model_consistent_expectation_symb_id, subst_table);
-
-  for (auto & equation : equations)
-    {
-      auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, model_consistent_expectation_symb_id, subst_table));
-      assert(substeq != nullptr);
-      equation = substeq;
-    }
-
-  for (map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator it = subst_table.begin();
-       it != subst_table.end(); it++)
-    pac_expectation_info.insert(const_cast<PacExpectationNode *>(it->first));
+                                            nonstationary, growth_symb_id, growth_lag, i);
 }
 
 void
@@ -4480,11 +4527,15 @@ DynamicModel::substitutePacExpectation(const string & name)
 {
   map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
   for (auto & it : local_variables_table)
-    it.second = it.second->substitutePacExpectation(name, subst_table);
+    it.second = pac_mce_z1_symb_ids.empty() ?
+      it.second->substitutePacExpectation(name, subst_table) :
+      it.second->substitutePacExpectation(name, pac_mce_z1_symb_ids, subst_table);
 
   for (auto & equation : equations)
     {
-      auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, subst_table));
+      auto *substeq = pac_mce_z1_symb_ids.empty() ?
+        dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, subst_table)) :
+        dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, pac_mce_z1_symb_ids, subst_table));
       assert(substeq != nullptr);
       equation = substeq;
     }
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 03b168998d7ed82f0ae24b1792870cbcd4aefa30..3b1d9f1aac3b17f898f742567b686a3b55a3c60f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -344,21 +344,18 @@ public:
   void addEquationsForVar();
 
   //! Get Pac equation parameter info
-  void walkPacParameters();
+  void walkPacParameters(const string &name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag);
   //! Add var_model info to pac_expectation nodes
   void fillPacExpectationVarInfo(const string &pac_model_name,
                                  vector<int> &lhs,
                                  int max_lag,
-                                 int pac_max_lag,
-                                 vector<bool> &nonstationary,
+                                 const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+                                 const vector<bool> &nonstationary,
                                  int growth_symb_id, int growth_lag);
 
   //! Substitutes pac_expectation operator with expectation based on auxiliary model
   void substitutePacExpectation(const string & name);
 
-  //! Substitute pac_expectation operator with model consistent expectation
-  void substitutePacExpectation(const string & name, int model_consistent_expectation_symb_id);
-
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &basename,
                                    const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const;
@@ -455,16 +452,20 @@ public:
   void substituteVarExpectation(const map<string, expr_t> &subst_table);
 
   //! Return max lag of pac equation
-  int getPacMaxLag(const string &pac_model_name) const;
+  void getPacMaxLag(const string &pac_model_name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag) const;
 
   //! Return target of the pac equation
   int getPacTargetSymbId(const string &pac_model_name) const;
 
   //! Add model consistent expectation equation for pac model
-  int addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id, int discount, int pac_max_lag_m, ExprNode::subst_table_t &diff_subst_table);
+  void addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id, int discount, const map<pair<string, string>, pair<string, int>> &eqtag_and_lag, ExprNode::subst_table_t &diff_subst_table);
 
   //! store symb_ids for alphas created in addPacModelConsistentExpectationEquation
-  map<string, vector<int>> pac_mce_alpha_symb_ids;
+  map<pair<string, string>, vector<int>> pac_mce_alpha_symb_ids;
+  //! store symb_ids for z1s created in addPacModelConsistentExpectationEquation
+  map<pair<string, string>, int> pac_mce_z1_symb_ids;
+  //! Store lag info for pac equations
+  map<pair<string, string>, pair<string, int>> pac_eqtag_and_lag;
 
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index feedee5644bd72053486cf974322e3fffa2e992f..48d983444f5baa274eda8dd45b26f4b485993d28 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -601,7 +601,7 @@ NumConstNode::substitutePacExpectation(const string & name, map<const PacExpecta
 }
 
 expr_t
-NumConstNode::substitutePacExpectation(const string & name, int mce_symb_id,
+NumConstNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<NumConstNode *>(this);
@@ -695,12 +695,12 @@ NumConstNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+NumConstNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
 }
 
 void
-NumConstNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+NumConstNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -1642,7 +1642,7 @@ VariableNode::substitutePacExpectation(const string & name, map<const PacExpecta
 }
 
 expr_t
-VariableNode::substitutePacExpectation(const string & name, int mce_symb_id,
+VariableNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VariableNode *>(this);
@@ -2013,12 +2013,12 @@ VariableNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+VariableNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
 }
 
 void
-VariableNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+VariableNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -3649,10 +3649,10 @@ UnaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectat
 }
 
 expr_t
-UnaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+UnaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t argsubst = arg->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t argsubst = arg->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
@@ -3851,15 +3851,15 @@ UnaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+UnaryOpNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
-  arg->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
 }
 
 void
-UnaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+UnaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
-  arg->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
 }
 
 bool
@@ -5475,11 +5475,11 @@ BinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpecta
 }
 
 expr_t
-BinaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+BinaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substitutePacExpectation(name, mce_symb_id, subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t arg1subst = arg1->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -6047,17 +6047,17 @@ BinaryOpNode::getPacLHS(pair<int, int> &lhs)
 }
 
 void
-BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+BinaryOpNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
-  arg1->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
-  arg2->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg1->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg2->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
 }
 
 void
-BinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+BinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
-  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
-  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
 }
 
 bool
@@ -6885,12 +6885,12 @@ TrinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpect
 }
 
 expr_t
-TrinaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+TrinaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                         map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substitutePacExpectation(name, mce_symb_id, subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
-  expr_t arg3subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t arg1subst = arg1->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
+  expr_t arg3subst = arg2->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -6994,19 +6994,19 @@ TrinaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+TrinaryOpNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
-  arg1->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
-  arg2->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
-  arg3->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg1->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg2->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+  arg3->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
 }
 
 void
-TrinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+TrinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
-  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
-  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
-  arg3->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+  arg3->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
 }
 
 bool
@@ -7419,12 +7419,12 @@ AbstractExternalFunctionNode::substitutePacExpectation(const string & name, map<
 }
 
 expr_t
-AbstractExternalFunctionNode::substitutePacExpectation(const string & name, int mce_symb_id,
+AbstractExternalFunctionNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
    vector<expr_t> arguments_subst;
   for (auto argument : arguments)
-    arguments_subst.push_back(argument->substitutePacExpectation(name, mce_symb_id, subst_table));
+    arguments_subst.push_back(argument->substitutePacExpectation(name, eqtag_and_symb_id, subst_table));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7586,17 +7586,17 @@ AbstractExternalFunctionNode::getPacOptimizingShareAndExprNodes(set<int> &optim_
 }
 
 void
-AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+AbstractExternalFunctionNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
   for (auto argument : arguments)
-    argument->addParamInfoToPac(lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
+    argument->addParamInfoToPac(eqtag_arg, lhs_arg, optim_share_arg, ec_params_and_vars_arg, ar_params_and_vars_arg, non_optim_vars_params_and_constants);
 }
 
 void
-AbstractExternalFunctionNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   for (auto argument : arguments)
-    argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
+    argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
 }
 
 bool
@@ -9126,7 +9126,7 @@ VarExpectationNode::substitutePacExpectation(const string & name, map<const PacE
 }
 
 expr_t
-VarExpectationNode::substitutePacExpectation(const string & name, int mce_symb_id,
+VarExpectationNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                              map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VarExpectationNode *>(this);
@@ -9239,12 +9239,12 @@ VarExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants)
+VarExpectationNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants)
 {
 }
 
 void
-VarExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+VarExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -9353,29 +9353,30 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       return;
     }
 
-  output << "M_.pac." << model_name << ".lhs_var = "
-         << datatree.symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl
-         << "M_.pac." << model_name << ".max_lag = " << pac_max_lag << ";" << endl;
+  string substruct = model_name + ".equations." + eqtag + ".";
+
+  output << "M_.pac." << substruct << "lhs_var = "
+         << datatree.symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl;
 
   if (growth_symb_id >= 0)
-    output << "M_.pac." << model_name << ".growth_neutrality_param_index = "
+    output << "M_.pac." << substruct << "growth_neutrality_param_index = "
            << datatree.symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl;
 
   if (optim_share_index >= 0)
-    output << "M_.pac." << model_name << ".share_of_optimizing_agents_index = "
+    output << "M_.pac." << substruct << "share_of_optimizing_agents_index = "
            << datatree.symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
 
-  output << "M_.pac." << model_name << ".ec.params = "
+  output << "M_.pac." << substruct << "ec.params = "
          << datatree.symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
-         << "M_.pac." << model_name << ".ec.vars = [";
+         << "M_.pac." << substruct << "ec.vars = [";
   for (auto it : ec_params_and_vars.second.first)
       output << datatree.symbol_table.getTypeSpecificID(it) + 1 << " ";
   output << "];" << endl
-         << "M_.pac." << model_name << ".ec.isendo = [";
+         << "M_.pac." << substruct << "ec.isendo = [";
   for (auto it : ec_params_and_vars.second.second)
     output << (it ? "true" : "false") << " ";
   output << "];" << endl
-         << "M_.pac." << model_name << ".ar.params = [";
+         << "M_.pac." << substruct << "ar.params = [";
   for (auto it = ar_params_and_vars.begin();
        it != ar_params_and_vars.end(); it++)
     {
@@ -9384,7 +9385,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       output << datatree.symbol_table.getTypeSpecificID(it->first) + 1;
     }
   output << "];" << endl
-         << "M_.pac." << model_name << ".ar.vars = [";
+         << "M_.pac." << substruct << "ar.vars = [";
   for (auto it = ar_params_and_vars.begin();
        it != ar_params_and_vars.end(); it++)
     {
@@ -9393,7 +9394,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       output << datatree.symbol_table.getTypeSpecificID(it->second.first) + 1;
     }
   output << "];" << endl
-         << "M_.pac." << model_name << ".ar.lags = [";
+         << "M_.pac." << substruct << "ar.lags = [";
   for (auto it = ar_params_and_vars.begin();
        it != ar_params_and_vars.end(); it++)
     {
@@ -9404,7 +9405,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
   output << "];" << endl;
   if (!non_optim_vars_params_and_constants.empty())
     {
-      output << "M_.pac." << model_name << ".non_optimizing_behaviour.params = [";
+      output << "M_.pac." << substruct << "non_optimizing_behaviour.params = [";
       for (auto it = non_optim_vars_params_and_constants.begin();
            it != non_optim_vars_params_and_constants.end(); ++it)
         {
@@ -9416,7 +9417,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
             output << "NaN";
         }
       output << "];"
-             << "M_.pac." << model_name << ".non_optimizing_behaviour.vars = [";
+             << "M_.pac." << substruct << "non_optimizing_behaviour.vars = [";
       for (auto it = non_optim_vars_params_and_constants.begin();
            it != non_optim_vars_params_and_constants.end(); ++it)
         {
@@ -9425,7 +9426,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           output << datatree.symbol_table.getTypeSpecificID(get<0>(*it)) + 1;
         }
       output << "];" << endl
-             << "M_.pac." << model_name << ".non_optimizing_behaviour.lags = [";
+             << "M_.pac." << substruct << "non_optimizing_behaviour.lags = [";
       for (auto it = non_optim_vars_params_and_constants.begin();
            it != non_optim_vars_params_and_constants.end(); ++it)
         {
@@ -9434,7 +9435,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           output << get<1>(*it);
         }
       output << "];" << endl
-             << "M_.pac." << model_name << ".non_optimizing_behaviour.scaling_factor = [";
+             << "M_.pac." << substruct << "non_optimizing_behaviour.scaling_factor = [";
       for (auto it = non_optim_vars_params_and_constants.begin();
            it != non_optim_vars_params_and_constants.end(); ++it)
         {
@@ -9444,7 +9445,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         }
       output << "];" << endl;
     }
-  output << "M_.pac." << model_name << ".h0_param_indices = [";
+  output << "M_.pac." << substruct << "h0_param_indices = [";
   for (auto it = h0_indices.begin();
        it != h0_indices.end(); it++)
     {
@@ -9453,7 +9454,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       output << datatree.symbol_table.getTypeSpecificID(*it) + 1;
     }
   output << "];" << endl
-         << "M_.pac." << model_name << ".h1_param_indices = [";
+         << "M_.pac." << substruct << "h1_param_indices = [";
   for (auto it = h1_indices.begin();
        it != h1_indices.end(); it++)
     {
@@ -9844,8 +9845,14 @@ PacExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 }
 
 void
-PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants_arg)
+PacExpectationNode::addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants_arg)
 {
+  if (eqtag_arg == "")
+    {
+      cerr << "Pac Expectation: every equation with a pac_expectation operator must have been assigned an equation tag name" << endl;
+      exit(EXIT_FAILURE);
+    }
+
   if (lhs_arg.first == -1)
     {
       cerr << "Pac Expectation: error in obtaining LHS varibale." << endl;
@@ -9858,23 +9865,23 @@ PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_a
       exit(EXIT_FAILURE);
     }
 
-  lhs_pac_var = lhs_arg;
-  optim_share_index = optim_share_arg;
-  ar_params_and_vars = ar_params_and_vars_arg;
-  ec_params_and_vars = ec_params_and_vars_arg;
-  non_optim_vars_params_and_constants = non_optim_vars_params_and_constants_arg;
+  eqtag = move(eqtag_arg);
+  lhs_pac_var = move(lhs_arg);
+  optim_share_index = move(optim_share_arg);
+  ar_params_and_vars = move(ar_params_and_vars_arg);
+  ec_params_and_vars = move(ec_params_and_vars_arg);
+  non_optim_vars_params_and_constants = move(non_optim_vars_params_and_constants_arg);
 }
 
 
 void
-PacExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+PacExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   if (model_name != model_name_arg)
     return;
 
-  lhs = lhs_arg;
+  lhs = move(lhs_arg);
   max_lag = max_lag_arg;
-  pac_max_lag = pac_max_lag_arg;
   growth_symb_id = growth_symb_id_arg;
   growth_lag = growth_lag_arg;
   equation_number = equation_number_arg;
@@ -9949,13 +9956,18 @@ PacExpectationNode::substitutePacExpectation(const string & name, map<const PacE
 }
 
 expr_t
-PacExpectationNode::substitutePacExpectation(const string & name, int mce_symb_id,
-                                             map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+PacExpectationNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   if (model_name != name)
     return const_cast<PacExpectationNode *>(this);
 
-  expr_t subExpr = datatree.AddVariable(mce_symb_id);
+  if (eqtag == "")
+    {
+      cerr << "Dynare internal error: eqtag should have been assigned already" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  expr_t subExpr = datatree.AddVariable(eqtag_and_symb_id.at(make_pair(model_name, eqtag)));
   subst_table[const_cast<PacExpectationNode *>(this)] = dynamic_cast<BinaryOpNode *>(subExpr);
   return subExpr;
 }
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index dcbdc144ad90668bcbf7e13e92fd65a83a67d358..8b29e59b9343ec5a5ee8774ff9c8016682faa315 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -569,7 +569,7 @@ class ExprNode
       virtual expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
 
       //! Substitute pac_expectation operator with model consistent expectation
-      virtual expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+      virtual expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                               map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
 
       //! Add ExprNodes to the provided datatree
@@ -608,10 +608,10 @@ class ExprNode
                                                      expr_t &optim_part,
                                                      expr_t &non_optim_part) const = 0;
       //! Adds PAC equation param info to pac_expectation
-      virtual void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) = 0;
+      virtual void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) = 0;
 
       //! Fills var_model info for pac_expectation node
-      virtual void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) = 0;
+      virtual void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) = 0;
 
       //! Fills the AR matrix structure
       virtual void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const = 0;
@@ -726,7 +726,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -740,8 +740,8 @@ public:
   expr_t clone(DataTree &datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -820,7 +820,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -834,8 +834,8 @@ public:
   expr_t clone(DataTree &datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -942,7 +942,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -956,8 +956,8 @@ public:
   expr_t clone(DataTree &datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1071,7 +1071,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1091,8 +1091,8 @@ public:
   //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
   expr_t getNonZeroPartofEquation() const;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> ar_params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRowHelper(expr_t arg1, expr_t arg2,
                                    int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
@@ -1201,7 +1201,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1215,8 +1215,8 @@ public:
   expr_t clone(DataTree &datatree) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1330,7 +1330,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
@@ -1346,8 +1346,8 @@ public:
   expr_t clone(DataTree &datatree) const override = 0;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1544,7 +1544,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<tuple<int, expr_t, expr_t>>  &List_of_Op_RHS) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
@@ -1565,8 +1565,8 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1590,12 +1590,12 @@ class PacExpectationNode : public ExprNode
 public:
   const string model_name;
 private:
-  string var_model_name;
+  string var_model_name, eqtag;
   int growth_symb_id, growth_lag;
   bool stationary_vars_present, nonstationary_vars_present;
   vector<int> lhs;
   pair<int, int> lhs_pac_var;
-  int max_lag, pac_max_lag;
+  int max_lag;
   vector<int> h0_indices, h1_indices;
   int growth_param_index, equation_number;
   int optim_share_index;
@@ -1649,7 +1649,7 @@ public:
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
-  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+  expr_t substitutePacExpectation(const string & name, const map<pair<string, string>, int> &eqtag_and_symb_id,
                                   map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<tuple<int, expr_t, expr_t>>  &List_of_Op_RHS) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
@@ -1670,8 +1670,8 @@ public:
   expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
-  void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void addParamInfoToPac(string eqtag_arg, pair<int, int> lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars_arg, set<pair<int, pair<int, int>>> params_and_vars_arg, vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> lhs_arg, int max_lag_arg, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &A0, map<tuple<int, int, int>, expr_t> &A0star) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 34ba4819620e3401f35f31129c2599cecbef4bef..804d984c1eef741d2b78f7a453852d33f334c5c4 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -400,7 +400,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     // substitute only those unary ops that appear in auxiliary model equations
     dynamic_model.substituteUnaryOps(diff_static_model, unary_ops_nodes, unary_ops_subst_table, eqtags);
 
-  // Create auxiliary variable and equations for Diff operators that appear in VAR equations
+  // Create auxiliary variable and equations for Diff operators
   diff_table_t diff_table;
   ExprNode::subst_table_t diff_subst_table;
   dynamic_model.substituteDiff(diff_static_model, diff_table, diff_subst_table);
@@ -443,23 +443,20 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
                exit(EXIT_FAILURE);
              }
            pms->fillUndiffedLHS(lhs);
-           dynamic_model.walkPacParameters();
-           int pac_max_lag_m = original_model.getPacMaxLag(pms->name);
+           map<pair<string, string>, pair<string, int>> eqtag_and_lag;
+           dynamic_model.walkPacParameters(pms->name, eqtag_and_lag);
+           original_model.getPacMaxLag(pms->name, eqtag_and_lag);
            if (pms->aux_model_name == "")
              {
                int pac_target_symb_id = dynamic_model.getPacTargetSymbId(pms->name);
-               int model_consistent_expectation_symb_id =
-                 dynamic_model.addPacModelConsistentExpectationEquation(pms->name, pac_target_symb_id,
-                                                                        symbol_table.getID(pms->discount), pac_max_lag_m,
-                                                                        diff_subst_table);
-               dynamic_model.substitutePacExpectation(pms->name, model_consistent_expectation_symb_id);
+               dynamic_model.addPacModelConsistentExpectationEquation(pms->name, pac_target_symb_id,
+                                                                      symbol_table.getID(pms->discount), eqtag_and_lag,
+                                                                      diff_subst_table);
              }
            else
-             {
-               dynamic_model.fillPacExpectationVarInfo(pms->name, lhs, max_lag,
-                                                       pac_max_lag_m, nonstationary, pms->growth_symb_id, pms->growth_lag);
-               dynamic_model.substitutePacExpectation(pms->name);
-             }
+             dynamic_model.fillPacExpectationVarInfo(pms->name, lhs, max_lag,
+                                                     eqtag_and_lag, nonstationary, pms->growth_symb_id, pms->growth_lag);
+           dynamic_model.substitutePacExpectation(pms->name);
          }
      }