diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 535e7a08e4239db15a4b5446b68c9d184eaccfe2..aa703ff558aa8c1b14ef9e6f82c503a3a8a1b191 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5378,18 +5378,38 @@ DynamicModel::substituteAdl()
 }
 
 void
-DynamicModel::substituteUnaryOps(StaticModel &static_model, set<string> &var_model_eqtags)
+DynamicModel::getEquationNumbersFromTags(vector<int> &eqnumbers, set<string> &eqtags) const
 {
-  diff_table_t nodes;
-  vector<int> eqnumber;
-  for (auto & eqtag : var_model_eqtags)
+  for (auto & eqtag : eqtags)
     for (const auto & equation_tag : equation_tags)
       if (equation_tag.second.first == "name"
           && equation_tag.second.second == eqtag)
         {
-          eqnumber.push_back(equation_tag.first);
+          eqnumbers.push_back(equation_tag.first);
           break;
         }
+}
+
+void
+DynamicModel::findPacExpectationEquationNumbers(vector<int> &eqnumbers) const
+{
+  int i = 0;
+  for (auto & equation : equations)
+    {
+      if (equation->containsPacExpectation())
+        if (find(eqnumbers.begin(), eqnumbers.end(), i) == eqnumbers.end())
+          eqnumbers.push_back(i);
+      i++;
+    }
+}
+
+void
+DynamicModel::substituteUnaryOps(StaticModel &static_model, set<string> &var_model_eqtags)
+{
+  diff_table_t nodes;
+  vector<int> eqnumber;
+  getEquationNumbersFromTags(eqnumber, var_model_eqtags);
+  findPacExpectationEquationNumbers(eqnumber);
 
   // Find matching unary ops that may be outside of diffs (i.e., those with different lags)
   set<int> used_local_vars;
@@ -5430,15 +5450,24 @@ DynamicModel::substituteUnaryOps(StaticModel &static_model, set<string> &var_mod
 }
 
 void
-DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table)
+DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table, set<string> &var_model_eqtags)
 {
-  // Find diff Nodes
+  vector<int> eqnumbers;
+  getEquationNumbersFromTags(eqnumbers, var_model_eqtags);
+  findPacExpectationEquationNumbers(eqnumbers);
+
+  set<int> used_local_vars;
+  for (int eqnumber : eqnumbers)
+    equations[eqnumber]->collectVariables(eModelLocalVariable, used_local_vars);
+
+  // Only substitute diffs in model local variables that appear in VAR equations
   diff_table_t diff_table;
   for (auto & it : local_variables_table)
-    it.second->findDiffNodes(static_model, diff_table);
+    if (used_local_vars.find(it.first) != used_local_vars.end())
+      it.second->findDiffNodes(static_model, diff_table);
 
-  for (auto & equation : equations)
-    equation->findDiffNodes(static_model, diff_table);
+  for (int eqnumber : eqnumbers)
+    equations[eqnumber]->findDiffNodes(static_model, diff_table);
 
   // Substitute in model local variables
   vector<BinaryOpNode *> neweqs;
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 36141a2e90f13df3152dc2627f4dcd9fbcbacfd6..578a9ad523d2c027f7f32a79585b19a33252e0ea 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -245,6 +245,10 @@ private:
   //! Create a legacy *_dynamic.m file for Matlab/Octave not yet using the temporary terms array interface
   void writeDynamicMatlabCompatLayer(const string &name) const;
 
+  void getEquationNumbersFromTags(vector<int> &eqnumber, set<string> &eqtags) const;
+
+  void findPacExpectationEquationNumbers(vector<int> &eqnumber) const;
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_argx);
   //! Adds a variable node
@@ -423,7 +427,7 @@ public:
   void substituteUnaryOps(StaticModel &static_model, set<string> &eq_tags);
 
   //! Substitutes diff operator
-  void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table);
+  void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table, set<string> &var_model_eqtags);
 
   //! 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,
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 9ab518f11e982469b9ca0aed87cebd18f108dace..8619492facc33708e2e1b3f7dc76cb9282e588cd 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -600,6 +600,12 @@ NumConstNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 {
 }
 
+bool
+NumConstNode::containsPacExpectation() const
+{
+  return false;
+}
+
 bool
 NumConstNode::containsEndogenous() const
 {
@@ -1688,6 +1694,12 @@ VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
     return false;
 }
 
+bool
+VariableNode::containsPacExpectation() const
+{
+  return false;
+}
+
 bool
 VariableNode::containsEndogenous() const
 {
@@ -3127,8 +3139,11 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
   auto it = diff_table.find(sthis);
   if (it == diff_table.end() || it->second[-arg->maxLag()] != this)
     {
-      cerr << "Internal error encountered. Please report" << endl;
-      exit(EXIT_FAILURE);
+      // diff does not appear in VAR equations
+      // so simply substitute diff(x) with x-x(-1)
+      expr_t argsubst = arg->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+      return dynamic_cast<BinaryOpNode *>(datatree.AddMinus(argsubst,
+                                                            argsubst->decreaseLeadsLags(1)));
     }
 
   int last_arg_max_lag = 0;
@@ -3358,6 +3373,12 @@ UnaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag
   return false;
 }
 
+bool
+UnaryOpNode::containsPacExpectation() const
+{
+  return arg->containsPacExpectation();
+}
+
 bool
 UnaryOpNode::containsEndogenous() const
 {
@@ -4955,6 +4976,12 @@ BinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
   return false;
 }
 
+bool
+BinaryOpNode::containsPacExpectation() const
+{
+  return (arg1->containsPacExpectation() || arg2->containsPacExpectation());
+}
+
 bool
 BinaryOpNode::containsEndogenous() const
 {
@@ -5856,6 +5883,12 @@ TrinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int l
   return false;
 }
 
+bool
+TrinaryOpNode::containsPacExpectation() const
+{
+  return (arg1->containsPacExpectation() || arg2->containsPacExpectation() || arg3->containsPacExpectation());
+}
+
 bool
 TrinaryOpNode::containsEndogenous() const
 {
@@ -6332,6 +6365,16 @@ AbstractExternalFunctionNode::isVariableNodeEqualTo(SymbolType type_arg, int var
   return false;
 }
 
+
+bool
+AbstractExternalFunctionNode::containsPacExpectation() const
+{
+  bool result = false;
+  for (auto argument : arguments)
+    result = result || argument->containsPacExpectation();
+  return result;
+}
+
 bool
 AbstractExternalFunctionNode::containsEndogenous() const
 {
@@ -7838,6 +7881,12 @@ VarExpectationNode::differentiateForwardVars(const vector<string> &subset, subst
   return const_cast<VarExpectationNode *>(this);
 }
 
+bool
+VarExpectationNode::containsPacExpectation() const
+{
+  return false;
+}
+
 bool
 VarExpectationNode::containsEndogenous() const
 {
@@ -8280,6 +8329,12 @@ PacExpectationNode::differentiateForwardVars(const vector<string> &subset, subst
   return const_cast<PacExpectationNode *>(this);
 }
 
+bool
+PacExpectationNode::containsPacExpectation() const
+{
+  return true;
+}
+
 bool
 PacExpectationNode::containsEndogenous() const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 6eae1e05446efa5e8591346e1e0d3d1c558ef9e2..7b2b62808bfa3f54a6c3023d0e86ca20092b1d98 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -540,6 +540,9 @@ class ExprNode
       //! 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;
 
+      //! Returns true if PacExpectationNode encountered
+      virtual bool containsPacExpectation() const = 0;
+
       //! Fills map
       virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
     };
@@ -623,6 +626,7 @@ public:
   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;
   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;
@@ -713,6 +717,7 @@ public:
   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;
   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
@@ -827,6 +832,7 @@ public:
   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;
   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
@@ -963,6 +969,7 @@ public:
   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;
   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
@@ -1063,6 +1070,7 @@ public:
   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;
   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
@@ -1175,6 +1183,7 @@ public:
   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;
   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
@@ -1377,6 +1386,7 @@ public:
   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;
   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;
@@ -1462,6 +1472,7 @@ public:
   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;
   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 4aceca53bcd295dc8e2efc97ac3b4d769b2be67c..3332cf479ac966a9d2a3b88639ef0a0062fdd856 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -364,6 +364,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
         }
     }
 
+  // Get all equation tags associated with VARs and Pac Models
   string var_model_name;
   set<string> eqtags;
   map<string, vector<string>> var_model_eq_tags;
@@ -383,9 +384,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     // substitute only those unary ops that appear in VAR equations
     dynamic_model.substituteUnaryOps(diff_static_model, eqtags);
 
-  // Create auxiliary variable and equations for Diff operator
+  // Create auxiliary variable and equations for Diff operators that appear in VAR equations
   ExprNode::subst_table_t diff_subst_table;
-  dynamic_model.substituteDiff(diff_static_model, diff_subst_table);
+  dynamic_model.substituteDiff(diff_static_model, diff_subst_table, eqtags);
 
   // Var Model
   map<string, tuple<vector<int>, vector<expr_t>, vector<bool>, vector<int>, int, vector<bool>, vector<int>>>