diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 61f91926f05f4dab1bc730aee554d53ff6ed8964..63aae410efb2127db6865da7854309661cdc1dfa 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -300,22 +300,17 @@ PacModelStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolida
     }
 }
 
-void
-PacModelStatement::fillUndiffedLHS(vector<int> &lhs_arg)
-{
-  lhs = lhs_arg;
-}
-
 void
 PacModelStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
   output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name << "';" << endl
-         << "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount) + 1 << ";" << endl
-         << "M_.pac." << name << ".steady_state_growth_rate = ";
-  if (steady_state_growth_rate_symb_id < 0)
-    output << steady_state_growth_rate_number << ";" << endl;
-  else
-    output << symbol_table.getTypeSpecificID(steady_state_growth_rate_symb_id) + 1 << ";" << endl;
+         << "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount) + 1 << ";" << endl;
+  if (steady_state_growth_rate_symb_id < 0 && steady_state_growth_rate_number > 0)
+    output << "M_.pac." << name << ".steady_state_growth_rate = "
+           << steady_state_growth_rate_number << ";" << endl;
+  else if (steady_state_growth_rate_symb_id >= 0)
+    output << "M_.pac." << name << ".steady_state_growth_rate = "
+           << symbol_table.getTypeSpecificID(steady_state_growth_rate_symb_id) + 1 << ";" << endl;
 
   if (growth_symb_id >= 0)
     {
@@ -364,14 +359,6 @@ PacModelStatement::writeOutput(ostream &output, const string &basename, bool min
             }
         }
     }
-  output << "M_.pac." << name << ".lhs = [";
-  for (auto it = lhs.begin(); it !=lhs.end(); it++)
-    {
-      if (it != lhs.begin())
-        output << " ";
-      output << *it + 1;
-    }
-  output << "];" << endl;
 }
 
 void
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 5a37b8b7f8b07843a99b665572e71b9f87b8fc34..47c4446e492e23c545a8b03de4df4a2b513868fb 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -156,7 +156,6 @@ public:
   void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
   void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream &output) const override;
-  void fillUndiffedLHS(vector<int> &lhs);
 };
 
 class VarRestrictionsStatement : public Statement
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 44bb929c38179ac31fae3d503fcdcb02ca4780b6..ae512e56304c418af0c210601ec7a43e04192000 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -85,9 +85,6 @@ DynamicModel::copyHelper(const DynamicModel &m)
     derivative_exo.push_back(convert_derivative_t(it));
   for (const auto &it : m.derivative_exo_det)
     derivative_exo_det.push_back(convert_derivative_t(it));
-
-  for (const auto &it : m.pac_expectation_info)
-    pac_expectation_info.insert(dynamic_cast<const PacExpectationNode *>(f(it)));
 }
 
 
@@ -231,8 +228,6 @@ DynamicModel::operator=(const DynamicModel &m)
   equation_block = m.equation_block;
   var_expectation_functions_to_write = m.var_expectation_functions_to_write;
 
-  pac_expectation_info.clear();
-
   endo_max_leadlag_block = m.endo_max_leadlag_block;
   other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
   exo_max_leadlag_block = m.exo_max_leadlag_block;
@@ -3633,11 +3628,6 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   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)
@@ -3655,9 +3645,98 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
       output << "];" << endl;
     }
 
-  // Write PacExpectationInfo
-  for (auto it : pac_expectation_info)
-    it->ExprNode::writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
+  for (auto & it : pac_model_info)
+    {
+      int growth_param_index = get<6>(it.second);
+      if (growth_param_index >= 0)
+        output << modstruct << "pac." << it.first << ".growth_neutrality_param_index = "
+               << symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl;
+
+      vector<int> lhs = get<0>(it.second);
+      output << modstruct << "pac." << it.first << ".lhs = [";
+      for (auto it : lhs)
+        output << it + 1 << " ";
+      output << "];" << endl;
+    }
+
+  for (auto & pit : pac_equation_info)
+    {
+      pair<int, int> lhs_pac_var;
+      int optim_share_index;
+      set<pair<int, pair<int, int>>> ar_params_and_vars;
+      pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars;
+      vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants;
+      tie(lhs_pac_var, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants) = pit.second;
+      string substruct = pit.first.first + ".equations." + pit.first.second + ".";
+
+      output << "M_.pac." << substruct << "lhs_var = "
+             << symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl;
+
+      if (optim_share_index >= 0)
+        output << "M_.pac." << substruct << "share_of_optimizing_agents_index = "
+               << symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
+
+      output << "M_.pac." << substruct << "ec.params = "
+             << symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
+             << "M_.pac." << substruct << "ec.vars = [";
+      for (auto it : ec_params_and_vars.second.first)
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ec.isendo = [";
+      for (auto it : ec_params_and_vars.second.second)
+        output << (it ? "true" : "false") << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ar.params = [";
+      for (auto & it : ar_params_and_vars)
+        output << symbol_table.getTypeSpecificID(it.first) + 1 << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ar.vars = [";
+      for (auto & it : ar_params_and_vars)
+        output << symbol_table.getTypeSpecificID(it.second.first) + 1 << " ";
+      output << "];" << endl
+             << "M_.pac." << substruct << "ar.lags = [";
+      for (auto & it : ar_params_and_vars)
+        output << it.second.second << " ";
+      output << "];" << endl;
+      if (!non_optim_vars_params_and_constants.empty())
+        {
+          output << "M_.pac." << substruct << "non_optimizing_behavior.params = [";
+          for (auto & it : non_optim_vars_params_and_constants)
+            if (get<2>(it) >= 0)
+              output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
+            else
+              output << "NaN ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "non_optimizing_behavior.vars = [";
+          for (auto & it : non_optim_vars_params_and_constants)
+            output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "non_optimizing_behavior.lags = [";
+          for (auto & it : non_optim_vars_params_and_constants)
+            output << get<1>(it) << " ";
+          output << "];" << endl
+                 << "M_.pac." << substruct << "non_optimizing_behavior.scaling_factor = [";
+          for (auto & it : non_optim_vars_params_and_constants)
+            output << get<3>(it) << " ";
+          output << "];" << endl;
+        }
+    }
+
+  for (auto & it : pac_h0_indices)
+    {
+      output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".h0_param_indices = [";
+      for (auto it1 : it.second)
+        output << symbol_table.getTypeSpecificID(it1) + 1 << " ";
+      output << "];" << endl;
+    }
+
+  for (auto & it : pac_h1_indices)
+    {
+      output << "M_.pac." << it.first.first << ".equations." << it.first.second << ".h1_param_indices = [";
+      for (auto it1 : it.second)
+        output << symbol_table.getTypeSpecificID(it1) + 1 << " ";
+      output << "];" << endl;
+    }
 }
 
 map<tuple<int, int, int>, expr_t>
@@ -4297,7 +4376,9 @@ DynamicModel::walkPacParameters(const string &name, map<pair<string, string>, pa
           set<int> optim_share;
           expr_t optim_part = nullptr;
           expr_t non_optim_part = nullptr;
-          equation->getPacLHS(lhs);
+          set<pair<int, int>> lhss;
+          equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss);
+          lhs = *(lhss.begin());
           int lhs_orig_symb_id = lhs.first;
           if (symbol_table.isAuxiliaryVariable(lhs_orig_symb_id))
             try
@@ -4342,12 +4423,20 @@ DynamicModel::walkPacParameters(const string &name, map<pair<string, string>, pa
               cerr << "Every equation with a pac expectation must have been assigned an equation tag name" << endl;
               exit(EXIT_FAILURE);
             }
+          if (lhs.first == -1)
+            {
+              cerr << "walkPacParameters: error obtaining LHS varibale." << endl;
+              exit(EXIT_FAILURE);
+            }
+          if (ec_params_and_vars.second.first.empty() || ar_params_and_vars.empty())
+            {
+              cerr << "walkPacParameters: error obtaining RHS parameters." << 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);
+          pac_equation_info[make_pair(name, eq)] = make_tuple(lhs, optim_share_index,
+                                                              ar_params_and_vars, ec_params_and_vars,
+                                                              non_optim_vars_params_and_constants);
           eqtag_and_lag[make_pair(name, eqtag)] = make_pair(eq, 0);
         }
     }
@@ -4388,7 +4477,9 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
     if (equation->containsPacExpectation(pac_model_name))
       {
         pair<int, int> lhs (-1, -1);
-        equation->getPacLHS(lhs);
+        set<pair<int, int>> lhss;
+        equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss);
+        lhs = *(lhss.begin());
         int lhs_symb_id = lhs.first;
         int lhs_orig_symb_id = lhs_symb_id;
         if (symbol_table.isAuxiliaryVariable(lhs_symb_id))
@@ -4405,17 +4496,19 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
 }
 
 void
-DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id,
-                                                       int discount_symb_id, const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
+DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int discount_symb_id,
+                                                       const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
                                                        ExprNode::subst_table_t &diff_subst_table)
 {
+  int pac_target_symb_id = getPacTargetSymbId(name);
   pac_eqtag_and_lag.insert(eqtag_and_lag.begin(), eqtag_and_lag.end());
   int neqs = 0;
   for (auto & it : eqtag_and_lag)
     {
-      string eqtag = it.second.first;
+      string eqtag = it.first.second;
+      string standard_eqtag = it.second.first;
       int pac_max_lag_m = it.second.second;
-      string append_to_name = name + "_" + eqtag;
+      string append_to_name = name + "_" + standard_eqtag;
       int mce_z1_symb_id;
       try
         {
@@ -4435,7 +4528,7 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int
           {
             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);
+            pac_mce_alpha_symb_ids[make_pair(name, standard_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),
@@ -4503,46 +4596,106 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int
       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;
+      pac_mce_z1_symb_ids[make_pair(name, standard_eqtag)] = mce_z1_symb_id;
+      pac_expectation_substitution[make_pair(name, eqtag)] = AddVariable(mce_z1_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,
-                                        const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
-                                        const vector<bool> &nonstationary,
-                                        int growth_symb_id, int growth_lag)
+DynamicModel::fillPacModelInfo(const string &pac_model_name,
+                               vector<int> &lhs,
+                               int max_lag,
+                               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,
-                                            nonstationary, growth_symb_id, growth_lag, i);
-}
 
-void
-DynamicModel::substitutePacExpectation(const string & name)
-{
-  map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
-  for (auto & it : local_variables_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);
+  bool stationary_vars_present = false;
+  bool nonstationary_vars_present = false;
+  for (auto it : nonstationary)
+    if (nonstationary_vars_present && stationary_vars_present)
+      break;
+    else
+      if (it)
+        nonstationary_vars_present = true;
+      else
+        stationary_vars_present = true;
 
-  for (auto & equation : equations)
+  int growth_param_index = -1;
+   if (growth_symb_id >= 0)
+     growth_param_index = symbol_table.addSymbol(pac_model_name +
+                                                 "_pac_growth_neutrality_correction",
+                                                 SymbolType::parameter);
+
+   for (auto pac_models_and_eqtags : pac_eqtag_and_lag)
     {
-      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;
+      if (pac_models_and_eqtags.first.first != pac_model_name)
+        continue;
+      string eqtag = pac_models_and_eqtags.first.second;
+      string standard_eqtag = pac_models_and_eqtags.second.first;
+      expr_t subExpr = Zero;
+      if (stationary_vars_present)
+        for (int i = 1; i < max_lag + 1; i++)
+          for (auto lhsit : lhs)
+            {
+              stringstream param_name_h0;
+              param_name_h0 << "h0_" << pac_model_name
+                            << "_" << standard_eqtag
+                            << "_var_" << symbol_table.getName(lhsit)
+                            << "_lag_" << i;
+              int new_param_symb_id = symbol_table.addSymbol(param_name_h0.str(), SymbolType::parameter);
+              pac_h0_indices[make_pair(pac_model_name, standard_eqtag)].push_back(new_param_symb_id);
+              subExpr = AddPlus(subExpr,
+                                AddTimes(AddVariable(new_param_symb_id),
+                                         AddVariable(lhsit, -i)));
+            }
+
+      if (nonstationary_vars_present)
+        for (int i = 1; i < max_lag + 1; i++)
+          for (auto lhsit : lhs)
+            {
+              stringstream param_name_h1;
+              param_name_h1 << "h1_" << pac_model_name
+                            << "_" << standard_eqtag
+                            << "_var_" << symbol_table.getName(lhsit)
+                            << "_lag_" << i;
+              int new_param_symb_id = symbol_table.addSymbol(param_name_h1.str(), SymbolType::parameter);
+              pac_h1_indices[make_pair(pac_model_name, standard_eqtag)].push_back(new_param_symb_id);
+              subExpr = AddPlus(subExpr,
+                                AddTimes(AddVariable(new_param_symb_id),
+                                         AddVariable(lhsit, -i)));
+            }
+
+      if (growth_symb_id >= 0)
+        subExpr = AddPlus(subExpr,
+                          AddTimes(AddVariable(growth_param_index),
+                                   AddVariable(growth_symb_id, growth_lag)));
+
+      pac_expectation_substitution[make_pair(pac_model_name, eqtag)] = subExpr;
     }
 
-  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));
+  pac_model_info[pac_model_name] = make_tuple(lhs, max_lag,
+                                              nonstationary_vars_present, stationary_vars_present,
+                                              growth_symb_id, growth_lag, growth_param_index);
+}
+
+void
+DynamicModel::substitutePacExpectation(const string & pac_model_name)
+{
+  for (auto & it : pac_expectation_substitution)
+    if (it.first.first == pac_model_name)
+      for (auto & equation : equations)
+        for (auto & tag : equation_tags)
+          if (tag.first == (&equation - &equations[0]))
+            if (tag.second.first == "name" && tag.second.second == it.first.second)
+              {
+                auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(pac_model_name, it.second));
+                assert(substeq != nullptr);
+                equation = substeq;
+                break;
+              }
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 3b1d9f1aac3b17f898f742567b686a3b55a3c60f..10e4ab59b0df29427c7cd5b32e14ff2f5c23a884 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -233,9 +233,6 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int>> var_expectation_functions_to_write;
 
-  //! Used for pac_expectation operator
-  set<const PacExpectationNode *> pac_expectation_info; // PacExpectationNode pointers
-
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int>> endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
 
@@ -346,15 +343,15 @@ public:
   //! Get Pac equation parameter info
   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,
-                                 const map<pair<string, string>, pair<string, int>> &eqtag_and_lag,
-                                 const vector<bool> &nonstationary,
-                                 int growth_symb_id, int growth_lag);
+  void fillPacModelInfo(const string &pac_model_name,
+                        vector<int> &lhs,
+                        int max_lag,
+                        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);
+  void substitutePacExpectation(const string & pac_model_name);
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &basename,
@@ -458,15 +455,36 @@ public:
   int getPacTargetSymbId(const string &pac_model_name) const;
 
   //! Add model consistent expectation equation for pac model
-  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);
+  void addPacModelConsistentExpectationEquation(const string & name, 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
+  //! (pac_model_name, standardized_eqtag) -> mce_alpha_symb_id
   map<pair<string, string>, vector<int>> pac_mce_alpha_symb_ids;
+  //! store symb_ids for h0, h1 parameters
+  //! (pac_model_name, standardized_eqtag) -> parameter symb_ids
+  map<pair<string, string>, vector<int>> pac_h0_indices, pac_h1_indices;
   //! store symb_ids for z1s created in addPacModelConsistentExpectationEquation
+  //! (pac_model_name, standardized_eqtag) -> mce_z1_symb_id
   map<pair<string, string>, int> pac_mce_z1_symb_ids;
   //! Store lag info for pac equations
+  //! (pac_model_name, equation_tag) -> (standardized_eqtag, lag)
   map<pair<string, string>, pair<string, int>> pac_eqtag_and_lag;
 
+  //! (pac_model_name, equation_tag) -> expr_t
+  map<pair<string, string>, expr_t> pac_expectation_substitution;
+
+  //! Store info about pac models:
+  //! pac_model_name -> (lhsvars, max_lag, nonstationary_vars_present, stationary_vars_present, growth_symb_id, growth_lag, growth_param_index)
+  map<string, tuple<vector<int>, int, bool, bool, int, int, int>> pac_model_info;
+
+  //! Store info about pac models specific to the equation they appear in
+  //! (pac_model_name, standardized_eqtag) ->
+  //!     (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants)
+  map<pair<string, string>,
+      tuple<pair<int, int>, int, set<pair<int, pair<int, int>>>, pair<int, pair<vector<int>, vector<bool>>>, vector<tuple<int, int, int, double>>>> pac_equation_info;
+
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
                                  ExprNode::subst_table_t &diff_subst_table) const;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 48d983444f5baa274eda8dd45b26f4b485993d28..39c1fdc461b13d6af77f53522872daebc3c7bb87 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -595,14 +595,7 @@ NumConstNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &no
 }
 
 expr_t
-NumConstNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  return const_cast<NumConstNode *>(this);
-}
-
-expr_t
-NumConstNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+NumConstNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -694,16 +687,6 @@ NumConstNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 {
 }
 
-void
-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, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
-{
-}
-
 bool
 NumConstNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -1636,14 +1619,7 @@ VariableNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &no
 }
 
 expr_t
-VariableNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  return const_cast<VariableNode *>(this);
-}
-
-expr_t
-VariableNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+VariableNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
   return const_cast<VariableNode *>(this);
 }
@@ -2012,16 +1988,6 @@ VariableNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 {
 }
 
-void
-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, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
-{
-}
-
 bool
 VariableNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -3642,17 +3608,9 @@ UnaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nod
 }
 
 expr_t
-UnaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  expr_t argsubst = arg->substitutePacExpectation(name, subst_table);
-  return buildSimilarUnaryOpNode(argsubst, datatree);
-}
-
-expr_t
-UnaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                      map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+UnaryOpNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
-  expr_t argsubst = arg->substitutePacExpectation(name, eqtag_and_symb_id, subst_table);
+  expr_t argsubst = arg->substitutePacExpectation(name, subexpr);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
@@ -3850,18 +3808,6 @@ UnaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
   arg->getPacOptimizingShareAndExprNodes(optim_share, optim_part, non_optim_part);
 }
 
-void
-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(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, 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, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
-}
-
 bool
 UnaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -5467,19 +5413,10 @@ BinaryOpNode::countDiffs() const
 }
 
 expr_t
-BinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  expr_t arg1subst = arg1->substitutePacExpectation(name, subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(name, subst_table);
-  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
-}
-
-expr_t
-BinaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+BinaryOpNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
-  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 arg1subst = arg1->substitutePacExpectation(name, subexpr);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, subexpr);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -6037,29 +5974,6 @@ BinaryOpNode::fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, c
   arg2->fillErrorCorrectionRow(eqn, nontrend_lhs, trend_lhs, A0, A0star);
 }
 
-void
-BinaryOpNode::getPacLHS(pair<int, int> &lhs)
-{
-  set<pair<int, int>> general_lhs;
-  arg1->collectDynamicVariables(SymbolType::endogenous, general_lhs);
-  if (general_lhs.size() == 1)
-    lhs = *(general_lhs.begin());
-}
-
-void
-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(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, 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, 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
 BinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -6876,21 +6790,11 @@ TrinaryOpNode::countDiffs() const
 }
 
 expr_t
-TrinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  expr_t arg1subst = arg1->substitutePacExpectation(name, subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(name, subst_table);
-  expr_t arg3subst = arg3->substitutePacExpectation(name, subst_table);
-  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
-}
-
-expr_t
-TrinaryOpNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+TrinaryOpNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
-  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);
+  expr_t arg1subst = arg1->substitutePacExpectation(name, subexpr);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, subexpr);
+  expr_t arg3subst = arg3->substitutePacExpectation(name, subexpr);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -6993,22 +6897,6 @@ TrinaryOpNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
   arg3->getPacOptimizingShareAndExprNodes(optim_share, optim_part, non_optim_part);
 }
 
-void
-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(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, 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, 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
 TrinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -7410,21 +7298,11 @@ AbstractExternalFunctionNode::countDiffs() const
 }
 
 expr_t
-AbstractExternalFunctionNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+AbstractExternalFunctionNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
   vector<expr_t> arguments_subst;
   for (auto argument : arguments)
-    arguments_subst.push_back(argument->substitutePacExpectation(name, subst_table));
-  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
-}
-
-expr_t
-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, eqtag_and_symb_id, subst_table));
+    arguments_subst.push_back(argument->substitutePacExpectation(name, subexpr));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7452,7 +7330,6 @@ AbstractExternalFunctionNode::getIndxInTefTerms(int the_symb_id, const deriv_nod
   auto it = tef_terms.find({ the_symb_id, arguments });
   if (it != tef_terms.end())
     return it->second;
-  cout << endl << endl << tef_terms.size() <<  "." <<the_symb_id << endl << endl;
   throw UnknownFunctionNameAndArgs();
 }
 
@@ -7585,20 +7462,6 @@ AbstractExternalFunctionNode::getPacOptimizingShareAndExprNodes(set<int> &optim_
     argument->getPacOptimizingShareAndExprNodes(optim_share, optim_part, non_optim_part);
 }
 
-void
-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(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, 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, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
-}
-
 bool
 AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -9120,14 +8983,7 @@ VarExpectationNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table
 }
 
 expr_t
-VarExpectationNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
-{
-  return const_cast<VarExpectationNode *>(this);
-}
-
-expr_t
-VarExpectationNode::substitutePacExpectation(const string &name, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                             map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+VarExpectationNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -9238,16 +9094,6 @@ VarExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 {
 }
 
-void
-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, const vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
-{
-}
-
 expr_t
 VarExpectationNode::substituteStaticAuxiliaryVariable() const
 {
@@ -9346,123 +9192,11 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                 const deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != ExprNodeOutputType::matlabOutsideModel);
-
   if (isLatexOutput(output_type))
     {
       output << "PAC_EXPECTATION" << LEFT_PAR(output_type) << model_name << RIGHT_PAR(output_type);
       return;
     }
-
-  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." << substruct << "growth_neutrality_param_index = "
-           << datatree.symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl;
-
-  if (optim_share_index >= 0)
-    output << "M_.pac." << substruct << "share_of_optimizing_agents_index = "
-           << datatree.symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
-
-  output << "M_.pac." << substruct << "ec.params = "
-         << datatree.symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
-         << "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." << substruct << "ec.isendo = [";
-  for (auto it : ec_params_and_vars.second.second)
-    output << (it ? "true" : "false") << " ";
-  output << "];" << endl
-         << "M_.pac." << substruct << "ar.params = [";
-  for (auto it = ar_params_and_vars.begin();
-       it != ar_params_and_vars.end(); it++)
-    {
-      if (it != ar_params_and_vars.begin())
-        output << " ";
-      output << datatree.symbol_table.getTypeSpecificID(it->first) + 1;
-    }
-  output << "];" << endl
-         << "M_.pac." << substruct << "ar.vars = [";
-  for (auto it = ar_params_and_vars.begin();
-       it != ar_params_and_vars.end(); it++)
-    {
-      if (it != ar_params_and_vars.begin())
-        output << " ";
-      output << datatree.symbol_table.getTypeSpecificID(it->second.first) + 1;
-    }
-  output << "];" << endl
-         << "M_.pac." << substruct << "ar.lags = [";
-  for (auto it = ar_params_and_vars.begin();
-       it != ar_params_and_vars.end(); it++)
-    {
-      if (it != ar_params_and_vars.begin())
-        output << " ";
-      output << it->second.second;
-    }
-  output << "];" << endl;
-  if (!non_optim_vars_params_and_constants.empty())
-    {
-      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)
-        {
-          if (it != non_optim_vars_params_and_constants.begin())
-            output << " ";
-          if (get<2>(*it) >= 0)
-            output << datatree.symbol_table.getTypeSpecificID(get<2>(*it)) + 1;
-          else
-            output << "NaN";
-        }
-      output << "];"
-             << "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)
-        {
-          if (it != non_optim_vars_params_and_constants.begin())
-            output << " ";
-          output << datatree.symbol_table.getTypeSpecificID(get<0>(*it)) + 1;
-        }
-      output << "];" << endl
-             << "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)
-        {
-          if (it != non_optim_vars_params_and_constants.begin())
-            output << " ";
-          output << get<1>(*it);
-        }
-      output << "];" << endl
-             << "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)
-        {
-          if (it != non_optim_vars_params_and_constants.begin())
-            output << " ";
-          output << get<3>(*it);
-        }
-      output << "];" << endl;
-    }
-  output << "M_.pac." << substruct << "h0_param_indices = [";
-  for (auto it = h0_indices.begin();
-       it != h0_indices.end(); it++)
-    {
-      if (it != h0_indices.begin())
-        output << " ";
-      output << datatree.symbol_table.getTypeSpecificID(*it) + 1;
-    }
-  output << "];" << endl
-         << "M_.pac." << substruct << "h1_param_indices = [";
-  for (auto it = h1_indices.begin();
-       it != h1_indices.end(); it++)
-    {
-      if (it != h1_indices.begin())
-        output << " ";
-      output << datatree.symbol_table.getTypeSpecificID(*it) + 1;
-    }
-  output << "];" << endl;
 }
 
 int
@@ -9844,132 +9578,12 @@ PacExpectationNode::getPacOptimizingShareAndExprNodes(set<int> &optim_share,
 {
 }
 
-void
-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;
-      exit(EXIT_FAILURE);
-    }
-
-  if (ec_params_and_vars_arg.second.first.empty() || ar_params_and_vars_arg.empty())
-    {
-      cerr << "Pac Expectation: error in obtaining RHS parameters." << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  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, 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 = move(lhs_arg);
-  max_lag = max_lag_arg;
-  growth_symb_id = growth_symb_id_arg;
-  growth_lag = growth_lag_arg;
-  equation_number = equation_number_arg;
-
-  for (vector<bool>::const_iterator it = nonstationary_arg.begin();
-       it != nonstationary_arg.end(); it++)
-    {
-      if (*it)
-        nonstationary_vars_present = true;
-      else
-        stationary_vars_present = true;
-      if (nonstationary_vars_present && stationary_vars_present)
-        break;
-    }
-}
-
 expr_t
-PacExpectationNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+PacExpectationNode::substitutePacExpectation(const string & name, expr_t subexpr)
 {
   if (model_name != name)
     return const_cast<PacExpectationNode *>(this);
-
-  map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator myit =
-    subst_table.find(const_cast<PacExpectationNode *>(this));
-  if (myit != subst_table.end())
-    return const_cast<BinaryOpNode *>(myit->second);
-
-  expr_t subExpr = datatree.AddNonNegativeConstant("0");
-  if (stationary_vars_present)
-    for (int i = 1; i < max_lag + 1; i++)
-      for (vector<int>::const_iterator it = lhs.begin(); it != lhs.end(); it++)
-        {
-          stringstream param_name_h0;
-          param_name_h0 << "h0_" << model_name
-                        << "_var_" << datatree.symbol_table.getName(*it)
-                        << "_lag_" << i;
-          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h0.str(), SymbolType::parameter);
-          h0_indices.push_back(new_param_symb_id);
-          subExpr = datatree.AddPlus(subExpr,
-                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
-                                                       datatree.AddVariable(*it, -i)));
-        }
-
-  if (nonstationary_vars_present)
-    for (int i = 1; i < max_lag + 1; i++)
-      for (vector<int>::const_iterator it = lhs.begin(); it != lhs.end(); it++)
-        {
-          stringstream param_name_h1;
-          param_name_h1 << "h1_" << model_name
-                        << "_var_" << datatree.symbol_table.getName(*it)
-                        << "_lag_" << i;
-          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h1.str(), SymbolType::parameter);
-          h1_indices.push_back(new_param_symb_id);
-          subExpr = datatree.AddPlus(subExpr,
-                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
-                                                       datatree.AddVariable(*it, -i)));
-        }
-
-  if (growth_symb_id >= 0)
-    {
-      growth_param_index = datatree.symbol_table.addSymbol(model_name +
-                                                           "_pac_growth_neutrality_correction",
-                                                           SymbolType::parameter);
-      subExpr = datatree.AddPlus(subExpr,
-                                 datatree.AddTimes(datatree.AddVariable(growth_param_index),
-                                                   datatree.AddVariable(growth_symb_id, growth_lag)));
-    }
-
-  subst_table[const_cast<PacExpectationNode *>(this)] = dynamic_cast<BinaryOpNode *>(subExpr);
-
-  return subExpr;
-}
-
-expr_t
-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);
-
-  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;
+  return subexpr;
 }
 
 void
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 8b29e59b9343ec5a5ee8774ff9c8016682faa315..eaee71aa80e62e20c4fdbd726a7c760739f5aec8 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -566,11 +566,7 @@ class ExprNode
       virtual expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
-      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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                              map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
+      virtual expr_t substitutePacExpectation(const string & name, expr_t subexpr) = 0;
 
       //! Add ExprNodes to the provided datatree
       virtual expr_t clone(DataTree &datatree) const = 0;
@@ -607,11 +603,6 @@ class ExprNode
       virtual void getPacOptimizingShareAndExprNodes(set<int> &optim_share,
                                                      expr_t &optim_part,
                                                      expr_t &non_optim_part) const = 0;
-      //! Adds PAC equation param info to pac_expectation
-      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, 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;
@@ -725,9 +716,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -740,8 +729,6 @@ 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(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;
@@ -819,9 +806,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -834,8 +819,6 @@ 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(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;
@@ -941,9 +924,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -956,8 +937,6 @@ 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(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;
@@ -1036,7 +1015,6 @@ public:
                                   int lhs_orig_symb_id,
                                   pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars,
                                   set<pair<int, pair<int, int>>> &ar_params_and_vars) const;
-  void getPacLHS(pair<int, int> &lhs);
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<tuple<int, expr_t, expr_t>>  &List_of_Op_RHS) const override;
@@ -1070,9 +1048,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1091,8 +1067,6 @@ 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(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;
@@ -1200,9 +1174,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1215,8 +1187,6 @@ 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(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;
@@ -1329,9 +1299,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) override;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1346,8 +1314,6 @@ 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(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;
@@ -1543,9 +1509,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) 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,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1565,8 +1529,6 @@ 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(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;
@@ -1589,19 +1551,6 @@ class PacExpectationNode : public ExprNode
 {
 public:
   const string model_name;
-private:
-  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;
-  vector<int> h0_indices, h1_indices;
-  int growth_param_index, equation_number;
-  int optim_share_index;
-  pair<int, pair<vector<int>, vector<bool>>> ec_params_and_vars;
-  set<pair<int, pair<int, int>>> ar_params_and_vars;
-  vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants;
 public:
   PacExpectationNode(DataTree &datatree_arg, int idx_arg, string model_name);
   void computeTemporaryTerms(const pair<int, int> &derivOrder,
@@ -1648,9 +1597,7 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   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, const map<pair<string, string>, int> &eqtag_and_symb_id,
-                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, expr_t subexpr) 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,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1670,8 +1617,6 @@ 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(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 804d984c1eef741d2b78f7a453852d33f334c5c4..6cdea487ad96a0dfbecf89c1007f74384294e0a5 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -442,20 +442,15 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
                cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
                exit(EXIT_FAILURE);
              }
-           pms->fillUndiffedLHS(lhs);
            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);
-               dynamic_model.addPacModelConsistentExpectationEquation(pms->name, pac_target_symb_id,
-                                                                      symbol_table.getID(pms->discount), eqtag_and_lag,
-                                                                      diff_subst_table);
-             }
+             dynamic_model.addPacModelConsistentExpectationEquation(pms->name, symbol_table.getID(pms->discount),
+                                                                    eqtag_and_lag, diff_subst_table);
            else
-             dynamic_model.fillPacExpectationVarInfo(pms->name, lhs, max_lag,
-                                                     eqtag_and_lag, nonstationary, pms->growth_symb_id, pms->growth_lag);
+             dynamic_model.fillPacModelInfo(pms->name, lhs, max_lag,
+                                            eqtag_and_lag, nonstationary, pms->growth_symb_id, pms->growth_lag);
            dynamic_model.substitutePacExpectation(pms->name);
          }
      }