diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 23d3030617f1a8a501901ad4cae911d94547d5d2..e81c10cc1d4c05755224b205062d5e3d710f2697 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3759,15 +3759,35 @@ DynamicModel::walkPacParameters()
     }
 }
 
+int
+DynamicModel::getPacMaxLag(const string &pac_model_name) const
+{
+  for (auto & equation : equations)
+    if (equation->containsPacExpectation(pac_model_name))
+      {
+        set<pair<int, int>> endogs;
+        equation->get_arg1()->collectDynamicVariables(SymbolType::endogenous, endogs);
+        if (endogs.size() != 1)
+          {
+            cerr << "The LHS of the PAC equation may only be comprised of one endogenous variable"
+                 << endl;
+            exit(EXIT_FAILURE);
+          }
+        return equation->PacMaxLag((*(endogs.begin())).first);
+      }
+  return 0;
+}
+
 void
 DynamicModel::fillPacExpectationVarInfo(string &pac_model_name,
                                         vector<int> &lhs,
                                         int max_lag,
+                                        int pac_max_lag,
                                         vector<bool> &nonstationary,
                                         int growth_symb_id)
 {
   for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->fillPacExpectationVarInfo(pac_model_name, lhs, max_lag, nonstationary, growth_symb_id, i);
+    equations[i]->fillPacExpectationVarInfo(pac_model_name, lhs, max_lag, pac_max_lag, nonstationary, growth_symb_id, i);
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 9d014c492bd463211fb71496b6eccc81e9a59e66..f87e68962f9cd6e5d6f3e72787241218f6c65cd2 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -326,6 +326,7 @@ public:
   void fillPacExpectationVarInfo(string &pac_model_name,
                                  vector<int> &lhs,
                                  int max_lag,
+                                 int pac_max_lag,
                                  vector<bool> &nonstationary,
                                  int growth_symb_id);
 
@@ -434,6 +435,9 @@ public:
   //! Substitute VarExpectation operators
   void substituteVarExpectation(const map<string, expr_t> &subst_table);
 
+  //! Return max lag of pac equation
+  int getPacMaxLag(const string &pac_model_name) const;
+
   //! Table to undiff LHS variables for pac vector z
   void getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, vector<bool> &diff, vector<int> &orig_diff_var,
                           vector<int> &eqnumber, map<string, int> &undiff, ExprNode::subst_table_t &diff_subst_table);
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 8a09162b42152a7f5d14df10e41e047212c3368f..0c2b089e8a0c3577dfac04fad26dc77ac2a112e0 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -494,7 +494,7 @@ NumConstNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) cons
 }
 
 int
-NumConstNode::PacMaxLag(vector<int> &lhs) const
+NumConstNode::PacMaxLag(int lhs_symb_id) const
 {
   return 0;
 }
@@ -608,7 +608,7 @@ NumConstNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 }
 
 bool
-NumConstNode::containsPacExpectation() const
+NumConstNode::containsPacExpectation(const string &pac_model_name) const
 {
   return false;
 }
@@ -660,7 +660,7 @@ NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int,
 }
 
 void
-NumConstNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+NumConstNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
 }
 
@@ -1409,9 +1409,11 @@ VariableNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) cons
 }
 
 int
-VariableNode::PacMaxLag(vector<int> &lhs) const
+VariableNode::PacMaxLag(int lhs_symb_id) const
 {
-  return -lag;
+  if (lhs_symb_id == symb_id)
+    return -lag;
+  return 0;
 }
 
 expr_t
@@ -1690,7 +1692,7 @@ VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
 }
 
 bool
-VariableNode::containsPacExpectation() const
+VariableNode::containsPacExpectation(const string &pac_model_name) const
 {
   return false;
 }
@@ -1811,7 +1813,7 @@ VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int,
 }
 
 void
-VariableNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+VariableNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
 }
 
@@ -2972,10 +2974,10 @@ UnaryOpNode::VarMinLag() const
 }
 
 int
-UnaryOpNode::PacMaxLag(vector<int> &lhs) const
+UnaryOpNode::PacMaxLag(int lhs_symb_id) const
 {
   //This will never be an UnaryOpcode::diff node
-  return arg->PacMaxLag(lhs);
+  return arg->PacMaxLag(lhs_symb_id);
 }
 
 expr_t
@@ -3350,9 +3352,9 @@ UnaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag
 }
 
 bool
-UnaryOpNode::containsPacExpectation() const
+UnaryOpNode::containsPacExpectation(const string &pac_model_name) const
 {
-  return arg->containsPacExpectation();
+  return arg->containsPacExpectation(pac_model_name);
 }
 
 bool
@@ -3412,9 +3414,9 @@ UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int,
 }
 
 void
-UnaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+UnaryOpNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
-  arg->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
 }
 
 bool
@@ -4732,9 +4734,9 @@ BinaryOpNode::undiff() const
 }
 
 int
-BinaryOpNode::PacMaxLag(vector<int> &lhs) const
+BinaryOpNode::PacMaxLag(int lhs_symb_id) const
 {
-  return max(arg1->PacMaxLag(lhs), arg2->PacMaxLag(lhs));
+  return max(arg1->PacMaxLag(lhs_symb_id), arg2->PacMaxLag(lhs_symb_id));
 }
 
 expr_t
@@ -4956,9 +4958,9 @@ BinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
 }
 
 bool
-BinaryOpNode::containsPacExpectation() const
+BinaryOpNode::containsPacExpectation(const string &pac_model_name) const
 {
-  return (arg1->containsPacExpectation() || arg2->containsPacExpectation());
+  return (arg1->containsPacExpectation(pac_model_name) || arg2->containsPacExpectation(pac_model_name));
 }
 
 bool
@@ -5069,10 +5071,10 @@ BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int,
 }
 
 void
-BinaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+BinaryOpNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
-  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
-  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
 }
 
 bool
@@ -5693,9 +5695,9 @@ TrinaryOpNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) con
 }
 
 int
-TrinaryOpNode::PacMaxLag(vector<int> &lhs) const
+TrinaryOpNode::PacMaxLag(int lhs_symb_id) const
 {
-  return max(arg1->PacMaxLag(lhs), max(arg2->PacMaxLag(lhs), arg3->PacMaxLag(lhs)));
+  return max(arg1->PacMaxLag(lhs_symb_id), max(arg2->PacMaxLag(lhs_symb_id), arg3->PacMaxLag(lhs_symb_id)));
 }
 
 expr_t
@@ -5866,9 +5868,9 @@ TrinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int l
 }
 
 bool
-TrinaryOpNode::containsPacExpectation() const
+TrinaryOpNode::containsPacExpectation(const string &pac_model_name) const
 {
-  return (arg1->containsPacExpectation() || arg2->containsPacExpectation() || arg3->containsPacExpectation());
+  return (arg1->containsPacExpectation(pac_model_name) || arg2->containsPacExpectation(pac_model_name) || arg3->containsPacExpectation(pac_model_name));
 }
 
 bool
@@ -5933,11 +5935,11 @@ TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int
 }
 
 void
-TrinaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+TrinaryOpNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
-  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
-  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
-  arg3->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+  arg3->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
 }
 
 bool
@@ -6141,11 +6143,11 @@ AbstractExternalFunctionNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &
 }
 
 int
-AbstractExternalFunctionNode::PacMaxLag(vector<int> &lhs) const
+AbstractExternalFunctionNode::PacMaxLag(int lhs_symb_id) const
 {
   int val = 0;
   for (auto argument : arguments)
-    val = max(val, argument->PacMaxLag(lhs));
+    val = max(val, argument->PacMaxLag(lhs_symb_id));
   return val;
 }
 
@@ -6352,11 +6354,11 @@ AbstractExternalFunctionNode::isVariableNodeEqualTo(SymbolType type_arg, int var
 
 
 bool
-AbstractExternalFunctionNode::containsPacExpectation() const
+AbstractExternalFunctionNode::containsPacExpectation(const string &pac_model_name) const
 {
   bool result = false;
   for (auto argument : arguments)
-    result = result || argument->containsPacExpectation();
+    result = result || argument->containsPacExpectation(pac_model_name);
   return result;
 }
 
@@ -6429,10 +6431,10 @@ AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pai
 }
 
 void
-AbstractExternalFunctionNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
   for (auto argument : arguments)
-    argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
+    argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, equation_number_arg);
 }
 
 bool
@@ -7696,7 +7698,7 @@ VarExpectationNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs
 }
 
 int
-VarExpectationNode::PacMaxLag(vector<int> &lhs) const
+VarExpectationNode::PacMaxLag(int lhs_symb_id) const
 {
   cerr << "VarExpectationNode::PacMaxLag not implemented." << endl;
   exit(EXIT_FAILURE);
@@ -7879,7 +7881,7 @@ VarExpectationNode::differentiateForwardVars(const vector<string> &subset, subst
 }
 
 bool
-VarExpectationNode::containsPacExpectation() const
+VarExpectationNode::containsPacExpectation(const string &pac_model_name) const
 {
   return false;
 }
@@ -7969,7 +7971,7 @@ VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pai
 }
 
 void
-VarExpectationNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+VarExpectationNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
 }
 
@@ -8046,7 +8048,7 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
 
   output << "M_.pac." << model_name << ".lhs_var = "
          << datatree.symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl
-         << "M_.pac." << model_name << ".max_lag = " << max_lag << ";" << endl;
+         << "M_.pac." << model_name << ".max_lag = " << pac_max_lag << ";" << endl;
 
   if (growth_symb_id >= 0)
     output << "M_.pac." << model_name << ".growth_neutrality_param_index = "
@@ -8167,7 +8169,7 @@ PacExpectationNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs
 }
 
 int
-PacExpectationNode::PacMaxLag(vector<int> &lhs) const
+PacExpectationNode::PacMaxLag(int lhs_symb_id) const
 {
   return 0;
 }
@@ -8331,9 +8333,12 @@ PacExpectationNode::differentiateForwardVars(const vector<string> &subset, subst
 }
 
 bool
-PacExpectationNode::containsPacExpectation() const
+PacExpectationNode::containsPacExpectation(const string &pac_model_name) const
 {
-  return true;
+  if (pac_model_name.empty())
+    return true;
+  else
+    return pac_model_name == model_name;
 }
 
 bool
@@ -8446,13 +8451,14 @@ PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pai
 
 
 void
-PacExpectationNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg)
+PacExpectationNode::fillPacExpectationVarInfo(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 equation_number_arg)
 {
   if (model_name != model_name_arg)
     return;
 
   lhs = lhs_arg;
   max_lag = max_lag_arg;
+  pac_max_lag = pac_max_lag_arg;
   growth_symb_id = growth_symb_id_arg;
   equation_number = equation_number_arg;
 
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index e2710b7522ad39c007c39430d47c60aa789c9fdb..fd8320692371604de0187ee685389cedd0250a67 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -370,7 +370,7 @@ class ExprNode
 
       //! Get Max lag of var associated with Pac model
       //! Takes account of undiffed LHS variables in calculating the max lag
-      virtual int PacMaxLag(vector<int> &lhs) const = 0;
+      virtual int PacMaxLag(int lhs_symb_id) const = 0;
 
       virtual expr_t undiff() const = 0;
 
@@ -533,10 +533,10 @@ class ExprNode
       virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) = 0;
 
       //! Fills var_model info for pac_expectation node
-      virtual void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) = 0;
+      virtual void fillPacExpectationVarInfo(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 equation_number_arg) = 0;
 
       //! Returns true if PacExpectationNode encountered
-      virtual bool containsPacExpectation() const = 0;
+      virtual bool containsPacExpectation(const string &pac_model_name = "") const = 0;
 
       //! Fills map
       virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
@@ -591,7 +591,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -620,8 +620,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  virtual bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  virtual bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   expr_t substituteStaticAuxiliaryVariable() const override;
@@ -682,7 +682,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -711,8 +711,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
@@ -794,7 +794,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -826,8 +826,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
@@ -926,7 +926,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -963,8 +963,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
@@ -1033,7 +1033,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1064,8 +1064,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
@@ -1146,7 +1146,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1177,8 +1177,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
@@ -1335,7 +1335,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
@@ -1377,8 +1377,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   expr_t substituteStaticAuxiliaryVariable() const override;
@@ -1394,7 +1394,7 @@ private:
   bool stationary_vars_present, nonstationary_vars_present;
   vector<int> lhs;
   pair<int, int> lhs_pac_var;
-  int max_lag;
+  int max_lag, pac_max_lag;
   vector<int> h0_indices, h1_indices;
   int growth_param_index, equation_number;
   set<pair<int, pair<int, int>>> ec_params_and_vars;
@@ -1421,7 +1421,7 @@ public:
   int maxLag() const override;
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
-  int PacMaxLag(vector<int> &lhs) const override;
+  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
@@ -1463,8 +1463,8 @@ public:
   bool isInStaticForm() const override;
   void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int>>> &ec_params_and_vars, set<pair<int, pair<int, int>>> &params_and_vars) const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int equation_number_arg) override;
-  bool containsPacExpectation() const override;
+  void fillPacExpectationVarInfo(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 equation_number_arg) override;
+  bool containsPacExpectation(const string &pac_model_name = "") const override;
   bool isVarModelReferenced(const string &model_info_name) const override;
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   expr_t substituteStaticAuxiliaryVariable() const override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index b5fdf32233fa68ff90f3e93034f0ab9bdb120078..815a7d99043ded7c1a3e3e0a4acc9aa1b15e9d7a 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -437,7 +437,8 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
                  }
                pms->fillUndiffedLHS(lhs);
                dynamic_model.walkPacParameters();
-               dynamic_model.fillPacExpectationVarInfo(pac_model_name, lhs, max_lag, nonstationary, growth_symb_id);
+               int pac_max_lag = original_model.getPacMaxLag(pac_model_name);
+               dynamic_model.fillPacExpectationVarInfo(pac_model_name, lhs, max_lag, pac_max_lag, nonstationary, growth_symb_id);
                dynamic_model.substitutePacExpectation();
              }
          }