diff --git a/DynamicModel.cc b/DynamicModel.cc
index cee51df6..e282f1a8 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3296,31 +3296,31 @@ DynamicModel::sparseHelper(int order, ostream &output, int row_nb, int col_nb, E
 }
 
 void
-DynamicModel::substituteEndoLeadGreaterThanTwo()
+DynamicModel::substituteEndoLeadGreaterThanTwo(bool deterministic_model)
 {
-  substituteLeadLagInternal(avEndoLead);
+  substituteLeadLagInternal(avEndoLead, deterministic_model);
 }
 
 void
-DynamicModel::substituteEndoLagGreaterThanTwo()
+DynamicModel::substituteEndoLagGreaterThanTwo(bool deterministic_model)
 {
-  substituteLeadLagInternal(avEndoLag);
+  substituteLeadLagInternal(avEndoLag, deterministic_model);
 }
 
 void
-DynamicModel::substituteExoLead()
+DynamicModel::substituteExoLead(bool deterministic_model)
 {
-  substituteLeadLagInternal(avExoLead);
+  substituteLeadLagInternal(avExoLead, deterministic_model);
 }
 
 void
-DynamicModel::substituteExoLag()
+DynamicModel::substituteExoLag(bool deterministic_model)
 {
-  substituteLeadLagInternal(avExoLag);
+  substituteLeadLagInternal(avExoLag, deterministic_model);
 }
 
 void
-DynamicModel::substituteLeadLagInternal(aux_var_t type)
+DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model)
 {
   ExprNode::subst_table_t subst_table;
   vector<BinaryOpNode *> neweqs;
@@ -3333,13 +3333,13 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type)
       switch (type)
         {
         case avEndoLead:
-          subst = it->second->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
+          subst = it->second->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
           break;
         case avEndoLag:
           subst = it->second->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
         case avExoLead:
-          subst = it->second->substituteExoLead(subst_table, neweqs);
+          subst = it->second->substituteExoLead(subst_table, neweqs, deterministic_model);
           break;
         case avExoLag:
           subst = it->second->substituteExoLag(subst_table, neweqs);
@@ -3358,13 +3358,13 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type)
       switch (type)
         {
         case avEndoLead:
-          subst = equations[i]->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
+          subst = equations[i]->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
           break;
         case avEndoLag:
           subst = equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
         case avExoLead:
-          subst = equations[i]->substituteExoLead(subst_table, neweqs);
+          subst = equations[i]->substituteExoLead(subst_table, neweqs, deterministic_model);
           break;
         case avExoLag:
           subst = equations[i]->substituteExoLag(subst_table, neweqs);
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 1d42f281..05189ff9 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -185,7 +185,7 @@ private:
 
   //! Factorized code for substitutions of leads/lags
   /*! \param[in] type determines which type of variables is concerned */
-  void substituteLeadLagInternal(aux_var_t type);
+  void substituteLeadLagInternal(aux_var_t type, bool deterministic_model);
 
 private:
   //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
@@ -296,17 +296,17 @@ public:
 
   //! Transforms the model by removing all leads greater or equal than 2 on endos
   /*! Note that this can create new lags on endos and exos */
-  void substituteEndoLeadGreaterThanTwo();
+  void substituteEndoLeadGreaterThanTwo(bool deterministic_model);
 
   //! Transforms the model by removing all lags greater or equal than 2 on endos
-  void substituteEndoLagGreaterThanTwo();
+  void substituteEndoLagGreaterThanTwo(bool deterministic_model);
 
   //! Transforms the model by removing all leads on exos
   /*! Note that this can create new lags on endos and exos */
-  void substituteExoLead();
+  void substituteExoLead(bool deterministic_model);
 
   //! Transforms the model by removing all lags on exos
-  void substituteExoLag();
+  void substituteExoLag(bool deterministic_model);
 
   //! Transforms the model by removing all oExpectation
   void substituteExpectation(bool partial_information_model);
diff --git a/ExprNode.cc b/ExprNode.cc
index 9e0ea561..43c9afff 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -366,7 +366,7 @@ NumConstNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-NumConstNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -378,7 +378,7 @@ NumConstNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 expr_t
-NumConstNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -934,7 +934,7 @@ VariableNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   expr_t value;
   switch (type)
@@ -949,7 +949,7 @@ VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vecto
       if (value->maxEndoLead() <= 1)
         return const_cast<VariableNode *>(this);
       else
-        return value->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
+        return value->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
     default:
       return const_cast<VariableNode *>(this);
     }
@@ -1007,7 +1007,7 @@ VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 expr_t
-VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   expr_t value;
   switch (type)
@@ -1022,7 +1022,7 @@ VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode
       if (value->maxExoLead() == 0)
         return const_cast<VariableNode *>(this);
       else
-        return value->substituteExoLead(subst_table, neweqs);
+        return value->substituteExoLead(subst_table, neweqs, deterministic_model);
     default:
       return const_cast<VariableNode *>(this);
     }
@@ -1814,11 +1814,11 @@ UnaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-UnaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
-  if (op_code == oUminus)
+  if (op_code == oUminus || deterministic_model)
     {
-      expr_t argsubst = arg->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
+      expr_t argsubst = arg->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
   else
@@ -1838,11 +1838,11 @@ UnaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<
 }
 
 expr_t
-UnaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
-  if (op_code == oUminus)
+  if (op_code == oUminus || deterministic_model)
     {
-      expr_t argsubst = arg->substituteExoLead(subst_table, neweqs);
+      expr_t argsubst = arg->substituteExoLead(subst_table, neweqs, deterministic_model);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
   else
@@ -2854,39 +2854,48 @@ BinaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   expr_t arg1subst, arg2subst;
   int maxendolead1 = arg1->maxEndoLead(), maxendolead2 = arg2->maxEndoLead();
 
   if (maxendolead1 < 2 && maxendolead2 < 2)
     return const_cast<BinaryOpNode *>(this);
-
-  switch (op_code)
+  if (deterministic_model)
     {
-    case oPlus:
-    case oMinus:
-    case oEqual:
-      arg1subst = maxendolead1 >= 2 ? arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg1;
-      arg2subst = maxendolead2 >= 2 ? arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg2;
+      arg1subst = maxendolead1 >= 2 ? arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model) : arg1;
+      arg2subst = maxendolead2 >= 2 ? arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model) : arg2;
       return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
-    case oTimes:
-    case oDivide:
-      if (maxendolead1 >= 2 && maxendolead2 == 0 && arg2->maxExoLead() == 0)
-        {
-          arg1subst = arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
-          return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
-        }
-      if (maxendolead1 == 0 && arg1->maxExoLead() == 0
-          && maxendolead2 >= 2 && op_code == oTimes)
-        {
-          arg2subst = arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
-          return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
-        }
-      return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
-    default:
-      return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
     }
+  else
+    {
+      switch (op_code)
+        {
+        case oPlus:
+        case oMinus:
+        case oEqual:
+          arg1subst = maxendolead1 >= 2 ? arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model) : arg1;
+          arg2subst = maxendolead2 >= 2 ? arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model) : arg2;
+          return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+        case oTimes:
+        case oDivide:
+          if (maxendolead1 >= 2 && maxendolead2 == 0 && arg2->maxExoLead() == 0)
+            {
+              arg1subst = arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+              return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
+            }
+          if (maxendolead1 == 0 && arg1->maxExoLead() == 0
+              && maxendolead2 >= 2 && op_code == oTimes)
+            {
+              arg2subst = arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+              return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
+            }
+          return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+        default:
+          return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+        }
+    }
+
 }
 
 expr_t
@@ -2898,38 +2907,46 @@ BinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 expr_t
-BinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   expr_t arg1subst, arg2subst;
   int maxexolead1 = arg1->maxExoLead(), maxexolead2 = arg2->maxExoLead();
 
   if (maxexolead1 < 1 && maxexolead2 < 1)
     return const_cast<BinaryOpNode *>(this);
-
-  switch (op_code)
+  if (deterministic_model)
     {
-    case oPlus:
-    case oMinus:
-    case oEqual:
-      arg1subst = maxexolead1 >= 1 ? arg1->substituteExoLead(subst_table, neweqs) : arg1;
-      arg2subst = maxexolead2 >= 1 ? arg2->substituteExoLead(subst_table, neweqs) : arg2;
+      arg1subst = maxexolead1 >= 1 ? arg1->substituteExoLead(subst_table, neweqs, deterministic_model) : arg1;
+      arg2subst = maxexolead2 >= 1 ? arg2->substituteExoLead(subst_table, neweqs, deterministic_model) : arg2;
       return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
-    case oTimes:
-    case oDivide:
-      if (maxexolead1 >= 1 && maxexolead2 == 0 && arg2->maxEndoLead() == 0)
+    }
+  else
+    {
+      switch (op_code)
         {
-          arg1subst = arg1->substituteExoLead(subst_table, neweqs);
-          return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
+        case oPlus:
+        case oMinus:
+        case oEqual:
+          arg1subst = maxexolead1 >= 1 ? arg1->substituteExoLead(subst_table, neweqs, deterministic_model) : arg1;
+          arg2subst = maxexolead2 >= 1 ? arg2->substituteExoLead(subst_table, neweqs, deterministic_model) : arg2;
+          return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+        case oTimes:
+        case oDivide:
+          if (maxexolead1 >= 1 && maxexolead2 == 0 && arg2->maxEndoLead() == 0)
+            {
+              arg1subst = arg1->substituteExoLead(subst_table, neweqs, deterministic_model);
+              return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
+            }
+          if (maxexolead1 == 0 && arg1->maxEndoLead() == 0
+              && maxexolead2 >= 1 && op_code == oTimes)
+            {
+              arg2subst = arg2->substituteExoLead(subst_table, neweqs, deterministic_model);
+              return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
+            }
+          return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+        default:
+          return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
         }
-      if (maxexolead1 == 0 && arg1->maxEndoLead() == 0
-          && maxexolead2 >= 1 && op_code == oTimes)
-        {
-          arg2subst = arg2->substituteExoLead(subst_table, neweqs);
-          return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
-        }
-      return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
-    default:
-      return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
     }
 }
 
@@ -3425,10 +3442,17 @@ TrinaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-TrinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   if (maxEndoLead() < 2)
     return const_cast<TrinaryOpNode *>(this);
+  else if (deterministic_model)
+    {
+      expr_t arg1subst = arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+      expr_t arg2subst = arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+      expr_t arg3subst = arg3->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
+      return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+    }
   else
     return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
 }
@@ -3443,10 +3467,17 @@ TrinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vecto
 }
 
 expr_t
-TrinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   if (maxExoLead() == 0)
     return const_cast<TrinaryOpNode *>(this);
+  else if (deterministic_model)
+    {
+      expr_t arg1subst = arg1->substituteExoLead(subst_table, neweqs, deterministic_model);
+      expr_t arg2subst = arg2->substituteExoLead(subst_table, neweqs, deterministic_model);
+      expr_t arg3subst = arg3->substituteExoLead(subst_table, neweqs, deterministic_model);
+      return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+    }
   else
     return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
 }
@@ -3759,11 +3790,11 @@ ExternalFunctionNode::decreaseLeadsLagsPredeterminedVariables() const
 }
 
 expr_t
-ExternalFunctionNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+ExternalFunctionNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteEndoLeadGreaterThanTwo(subst_table, neweqs));
+    arguments_subst.push_back((*it)->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -3777,11 +3808,11 @@ ExternalFunctionNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table
 }
 
 expr_t
-ExternalFunctionNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+ExternalFunctionNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteExoLead(subst_table, neweqs));
+    arguments_subst.push_back((*it)->substituteExoLead(subst_table, neweqs, deterministic_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
diff --git a/ExprNode.hh b/ExprNode.hh
index 825e4343..e57d52e1 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -309,7 +309,7 @@ public:
 
     \return A new equivalent expression where sub-expressions with max endo lead >= 2 have been replaced by auxiliary variables
   */
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const = 0;
 
   //! Constructs a new expression where endo variables with max endo lag >= 2 have been replaced by auxiliary variables
   /*!
@@ -323,7 +323,7 @@ public:
     \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const = 0;
   //! Constructs a new expression where exogenous variables with a lag have been replaced by auxiliary variables
   /*!
     \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
@@ -396,9 +396,9 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -449,9 +449,9 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -512,11 +512,11 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   //! Creates another UnaryOpNode with the same opcode, but with a possibly different datatree and argument
   expr_t buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -582,11 +582,11 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   //! Creates another BinaryOpNode with the same opcode, but with a possibly different datatree and arguments
   expr_t buildSimilarBinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, DataTree &alt_datatree) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -634,11 +634,11 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   //! Creates another TrinaryOpNode with the same opcode, but with a possibly different datatree and arguments
   expr_t buildSimilarTrinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, expr_t alt_arg3, DataTree &alt_datatree) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -692,9 +692,9 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
-  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const;
diff --git a/ModFile.cc b/ModFile.cc
index c3009e2a..c346f784 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -220,11 +220,17 @@ ModFile::transformPass()
       || mod_file_struct.osr_present
       || mod_file_struct.ramsey_policy_present)
     {
-      // In stochastic models, create auxiliary vars for leads and lags greater than 2
-      dynamic_model.substituteEndoLeadGreaterThanTwo();
-      dynamic_model.substituteExoLead();
-      dynamic_model.substituteEndoLagGreaterThanTwo();
-      dynamic_model.substituteExoLag();
+      // In stochastic models, create auxiliary vars for leads and lags greater than 2, on both endos and exos
+      dynamic_model.substituteEndoLeadGreaterThanTwo(false);
+      dynamic_model.substituteExoLead(false);
+      dynamic_model.substituteEndoLagGreaterThanTwo(false);
+      dynamic_model.substituteExoLag(false);
+    }
+  else
+    {
+      // In deterministic models, create auxiliary vars for leads and lags endogenous greater than 2, only on endos (useless on exos)
+      dynamic_model.substituteEndoLeadGreaterThanTwo(true);
+      dynamic_model.substituteEndoLagGreaterThanTwo(true);
     }
 
   if (mod_file_struct.dsge_var_estimated || !mod_file_struct.dsge_var_calibrated.empty())