From 7bfeef5d7f04e07237ea39339a0d3e5248c9368c Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Wed, 7 Oct 2009 16:34:42 +0000
Subject: [PATCH] preprocessor: * fixed substitution of endogenous with leads
 >= 2: take into account exogenous with leads in non-linear terms * fixed
 substitution of exogenous with leads: take into account other variables with
 leads in non-linear terms

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3027 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 DynamicModel.cc |  55 ++++++----
 DynamicModel.hh |  13 ++-
 ExprNode.cc     | 286 +++++++++++++++++++++++++++++++++++++-----------
 ExprNode.hh     |  49 +++++++--
 ModFile.cc      |   3 +-
 SymbolTable.cc  |  50 ++++-----
 SymbolTable.hh  |  10 +-
 7 files changed, 338 insertions(+), 128 deletions(-)

diff --git a/DynamicModel.cc b/DynamicModel.cc
index fb22c354..57a19564 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -2833,23 +2833,29 @@ DynamicModel::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOut
 void
 DynamicModel::substituteEndoLeadGreaterThanTwo()
 {
-  substituteLeadLagInternal(0);
+  substituteLeadLagInternal(avEndoLead);
 }
 
 void
 DynamicModel::substituteEndoLagGreaterThanTwo()
 {
-  substituteLeadLagInternal(1);
+  substituteLeadLagInternal(avEndoLag);
 }
 
 void
-DynamicModel::substituteExoLeadLag()
+DynamicModel::substituteExoLead()
 {
-  substituteLeadLagInternal(2);
+  substituteLeadLagInternal(avExoLead);
 }
 
 void
-DynamicModel::substituteLeadLagInternal(int vars)
+DynamicModel::substituteExoLag()
+{
+  substituteLeadLagInternal(avExoLag);
+}
+
+void
+DynamicModel::substituteLeadLagInternal(aux_var_t type)
 {
   ExprNode::subst_table_t subst_table;
   vector<BinaryOpNode *> neweqs;
@@ -2859,16 +2865,19 @@ DynamicModel::substituteLeadLagInternal(int vars)
       it != local_variables_table.end(); it++)
     {
       NodeID subst;
-      switch(vars)
+      switch(type)
         {
-        case 0:
+        case avEndoLead:
           subst = it->second->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
           break;
-        case 1:
+        case avEndoLag:
           subst = it->second->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
-        case 2:
-          subst = it->second->substituteExoLeadLag(subst_table, neweqs);
+        case avExoLead:
+          subst = it->second->substituteExoLead(subst_table, neweqs);
+          break;
+        case avExoLag:
+          subst = it->second->substituteExoLag(subst_table, neweqs);
           break;
         }
       it->second = subst;
@@ -2878,16 +2887,19 @@ DynamicModel::substituteLeadLagInternal(int vars)
   for(int i = 0; i < (int) equations.size(); i++)
     {
       NodeID subst;
-      switch(vars)
+      switch(type)
         {
-        case 0:
+        case avEndoLead:
           subst = equations[i]->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
           break;
-        case 1:
+        case avEndoLag:
           subst = equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
           break;
-        case 2:
-          subst = equations[i]->substituteExoLeadLag(subst_table, neweqs);
+        case avExoLead:
+          subst = equations[i]->substituteExoLead(subst_table, neweqs);
+          break;
+        case avExoLag:
+          subst = equations[i]->substituteExoLag(subst_table, neweqs);
           break;
         }
       BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(subst);
@@ -2905,16 +2917,19 @@ DynamicModel::substituteLeadLagInternal(int vars)
   if (neweqs.size() > 0)
     {
       cout << "Substitution of ";
-      switch(vars)
+      switch(type)
         {
-        case 0:
+        case avEndoLead:
           cout << "endo leads >= 2";
           break;
-        case 1:
+        case avEndoLag:
           cout << "endo lags >= 2";
           break;
-        case 2:
-          cout << "exo leads/lags >= 2";
+        case avExoLead:
+          cout << "exo leads";
+          break;
+        case avExoLag:
+          cout << "exo lags";
           break;
         }
       cout << ": added " << neweqs.size() << " auxiliary variables and equations." << endl;
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 04bf98b5..14f1aa12 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -145,8 +145,8 @@ private:
   void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
 
   //! Factorized code for substitutions of leads/lags
-  /*! \param[in] vars 0 for endo leads, 1 for endo lags, 2 for exo */
-  void substituteLeadLagInternal(int vars);
+  /*! \param[in] type determines which type of variables is concerned */
+  void substituteLeadLagInternal(aux_var_t type);
 
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
@@ -203,13 +203,18 @@ public:
   virtual bool isDynamic() const { return true; };
 
   //! 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();
 
   //! Transforms the model by removing all lags greater or equal than 2 on endos
   void substituteEndoLagGreaterThanTwo();
 
-  //! Transforms the model by removing all leads and lags on exos
-  void substituteExoLeadLag();
+  //! Transforms the model by removing all leads on exos
+  /*! Note that this can create new lags on endos and exos */
+  void substituteExoLead();
+
+  //! Transforms the model by removing all lags on exos
+  void substituteExoLag();
 
   //! Fills eval context with values of model local variables and auxiliary variables
   void fillEvalContext(eval_context_type &eval_context) const;
diff --git a/ExprNode.cc b/ExprNode.cc
index 089502e3..2dffa162 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -183,6 +183,42 @@ ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector
   return dynamic_cast<VariableNode *>(substexpr);
 }
 
+VariableNode *
+ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  int n = maxExoLead();
+  assert(n >= 1);
+
+  subst_table_t::const_iterator it = subst_table.find(this);
+  if (it != subst_table.end())
+    return const_cast<VariableNode *>(it->second);
+
+  NodeID substexpr = decreaseLeadsLags(n);
+  int lag = n-1;
+
+  // Each iteration tries to create an auxvar such that auxvar(+1)=expr(-lag)
+  // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to expr(-lag-1) (resp. expr(-lag))
+  while(lag >= 0)
+    {
+      NodeID orig_expr = decreaseLeadsLags(lag);
+      it = subst_table.find(orig_expr);
+      if (it == subst_table.end())
+        {
+          int symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(orig_expr->idx);
+          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
+          substexpr = datatree.AddVariable(symb_id, +1);
+          assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
+          subst_table[orig_expr] = dynamic_cast<VariableNode *>(substexpr);
+        }
+      else
+        substexpr = const_cast<VariableNode *>(it->second);
+
+      lag--;
+    }
+
+  return dynamic_cast<VariableNode *>(substexpr);
+}
+
 
 NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
     ExprNode(datatree_arg),
@@ -271,6 +307,12 @@ NumConstNode::maxEndoLead() const
   return 0;
 }
 
+int
+NumConstNode::maxExoLead() const
+{
+  return 0;
+}
+
 NodeID
 NumConstNode::decreaseLeadsLags(int n) const
 {
@@ -290,7 +332,13 @@ NumConstNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 NodeID
-NumConstNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+NodeID
+NumConstNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -755,6 +803,20 @@ VariableNode::maxEndoLead() const
     }
 }
 
+int
+VariableNode::maxExoLead() const
+{
+  switch(type)
+    {
+    case eExogenous:
+      return max(lag, 0);
+    case eModelLocalVariable:
+      return datatree.local_variables_table[symb_id]->maxExoLead();
+    default:
+      return 0;
+    }
+}
+
 NodeID
 VariableNode::decreaseLeadsLags(int n) const
 {
@@ -840,7 +902,29 @@ VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 NodeID
-VariableNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  NodeID value;
+  switch(type)
+    {
+    case eExogenous:
+      if (lag <= 0)
+        return const_cast<VariableNode *>(this);
+      else
+        return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+    case eModelLocalVariable:
+      value = datatree.local_variables_table[symb_id];
+      if (value->maxExoLead() == 0)
+        return const_cast<VariableNode *>(this);
+      else
+        return value->substituteExoLead(subst_table, neweqs);
+    default:
+      return const_cast<VariableNode *>(this);
+    }
+}
+
+NodeID
+VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   VariableNode *substexpr;
   subst_table_t::const_iterator it;
@@ -848,66 +932,38 @@ VariableNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNo
   switch(type)
     {
     case eExogenous:
-      if (lag == 0)
+      if (lag >= 0)
         return const_cast<VariableNode *>(this);
 
       it = subst_table.find(this);
       if (it != subst_table.end())
         return const_cast<VariableNode *>(it->second);
 
-      if (lag < 0)
-        {
-          substexpr = datatree.AddVariable(symb_id, 0);
-          cur_lag = -1;
-
-          // Each iteration tries to create an auxvar such that auxvar(-1)=curvar(cur_lag)
-          // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to curvar(cur_lag+1) (resp. curvar(cur_lag))
-          while(cur_lag >= lag)
-            {
-              VariableNode *orig_expr = datatree.AddVariable(symb_id, cur_lag);
-              it = subst_table.find(orig_expr);
-              if (it == subst_table.end())
-                {
-                  int aux_symb_id = datatree.symbol_table.addExoLagAuxiliaryVar(symb_id, cur_lag+1);
-                  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), substexpr)));
-                  substexpr = datatree.AddVariable(aux_symb_id, -1);
-                  subst_table[orig_expr] = substexpr;
-                }
-              else
-                substexpr = const_cast<VariableNode *>(it->second);
+      substexpr = datatree.AddVariable(symb_id, 0);
+      cur_lag = -1;
 
-              cur_lag--;
-            }
-          return substexpr;
-        }
-      else
+      // Each iteration tries to create an auxvar such that auxvar(-1)=curvar(cur_lag)
+      // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to curvar(cur_lag+1) (resp. curvar(cur_lag))
+      while(cur_lag >= lag)
         {
-          substexpr = datatree.AddVariable(symb_id, 0);
-          cur_lag = +1;
-
-          // Each iteration tries to create an auxvar such that auxvar(+1)=curvar(cur_lag)
-          // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to curvar(cur_lag-1) (resp. curvar(cur_lag))
-          while(cur_lag <= lag)
+          VariableNode *orig_expr = datatree.AddVariable(symb_id, cur_lag);
+          it = subst_table.find(orig_expr);
+          if (it == subst_table.end())
             {
-              VariableNode *orig_expr = datatree.AddVariable(symb_id, cur_lag);
-              it = subst_table.find(orig_expr);
-              if (it == subst_table.end())
-                {
-                  int aux_symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(symb_id, cur_lag-1);
-                  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), substexpr)));
-                  substexpr = datatree.AddVariable(aux_symb_id, +1);
-                  subst_table[orig_expr] = substexpr;
-                }
-              else
-                substexpr = const_cast<VariableNode *>(it->second);
-
-              cur_lag++;
+              int aux_symb_id = datatree.symbol_table.addExoLagAuxiliaryVar(symb_id, cur_lag+1);
+              neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), substexpr)));
+              substexpr = datatree.AddVariable(aux_symb_id, -1);
+              subst_table[orig_expr] = substexpr;
             }
-          return substexpr;
+          else
+            substexpr = const_cast<VariableNode *>(it->second);
+
+          cur_lag--;
         }
+      return substexpr;
 
     case eModelLocalVariable:
-      return datatree.local_variables_table[symb_id]->substituteExoLeadLag(subst_table, neweqs);
+      return datatree.local_variables_table[symb_id]->substituteExoLag(subst_table, neweqs);
     default:
       return const_cast<VariableNode *>(this);
     }
@@ -1543,6 +1599,12 @@ UnaryOpNode::maxEndoLead() const
   return arg->maxEndoLead();
 }
 
+int
+UnaryOpNode::maxExoLead() const
+{
+  return arg->maxExoLead();
+}
+
 NodeID
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -1575,9 +1637,26 @@ UnaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<
 }
 
 NodeID
-UnaryOpNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID argsubst = arg->substituteExoLeadLag(subst_table, neweqs);
+  if (op_code == oUminus)
+    {
+      NodeID argsubst = arg->substituteExoLead(subst_table, neweqs);
+      return buildSimilarUnaryOpNode(argsubst, datatree);
+    }
+  else
+    {
+      if (maxExoLead() >= 1)
+        return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+      else
+        return const_cast<UnaryOpNode *>(this);
+    }
+}
+
+NodeID
+UnaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  NodeID argsubst = arg->substituteExoLag(subst_table, neweqs);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
@@ -2443,6 +2522,12 @@ BinaryOpNode::maxEndoLead() const
   return max(arg1->maxEndoLead(), arg2->maxEndoLead());
 }
 
+int
+BinaryOpNode::maxExoLead() const
+{
+  return max(arg1->maxExoLead(), arg2->maxExoLead());
+}
+
 NodeID
 BinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -2455,9 +2540,9 @@ NodeID
 BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   NodeID arg1subst, arg2subst;
-  int maxlead1 = arg1->maxEndoLead(), maxlead2 = arg2->maxEndoLead();
+  int maxendolead1 = arg1->maxEndoLead(), maxendolead2 = arg2->maxEndoLead();
 
-  if (maxlead1 < 2 && maxlead2 < 2)
+  if (maxendolead1 < 2 && maxendolead2 < 2)
     return const_cast<BinaryOpNode *>(this);
 
   switch(op_code)
@@ -2465,17 +2550,18 @@ BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vecto
     case oPlus:
     case oMinus:
     case oEqual:
-      arg1subst = maxlead1 >= 2 ? arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg1;
-      arg2subst = maxlead2 >= 2 ? arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg2;
+      arg1subst = maxendolead1 >= 2 ? arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg1;
+      arg2subst = maxendolead2 >= 2 ? arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs) : arg2;
       return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
     case oTimes:
     case oDivide:
-      if (maxlead1 >= 2 && maxlead2 == 0)
+      if (maxendolead1 >= 2 && maxendolead2 == 0 && arg2->maxExoLead())
         {
           arg1subst = arg1->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
           return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
         }
-      if (maxlead1 == 0 && maxlead2 >= 2 && op_code == oTimes)
+      if (maxendolead1 == 0 && arg1->maxExoLead() == 0
+          && maxendolead2 >= 2 && op_code == oTimes)
         {
           arg2subst = arg2->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
           return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
@@ -2495,10 +2581,46 @@ BinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
 }
 
 NodeID
-BinaryOpNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  NodeID arg1subst, arg2subst;
+  int maxexolead1 = arg1->maxExoLead(), maxexolead2 = arg2->maxExoLead();
+
+  if (maxexolead1 < 1 && maxexolead2 < 1)
+    return const_cast<BinaryOpNode *>(this);
+
+  switch(op_code)
+    {
+    case oPlus:
+    case oMinus:
+    case oEqual:
+      arg1subst = maxexolead1 >= 1 ? arg1->substituteExoLead(subst_table, neweqs) : arg1;
+      arg2subst = maxexolead2 >= 1 ? arg2->substituteExoLead(subst_table, neweqs) : arg2;
+      return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+    case oTimes:
+    case oDivide:
+      if (maxexolead1 >= 1 && maxexolead2 == 0 && arg2->maxEndoLead())
+        {
+          arg1subst = arg1->substituteExoLead(subst_table, neweqs);
+          return buildSimilarBinaryOpNode(arg1subst, arg2, datatree);
+        }
+      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);
+    }
+}
+
+NodeID
+BinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteExoLeadLag(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteExoLeadLag(subst_table, neweqs);
+  NodeID arg1subst = arg1->substituteExoLag(subst_table, neweqs);
+  NodeID arg2subst = arg2->substituteExoLag(subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -2851,6 +2973,12 @@ TrinaryOpNode::maxEndoLead() const
   return max(arg1->maxEndoLead(), max(arg2->maxEndoLead(), arg3->maxEndoLead()));
 }
 
+int
+TrinaryOpNode::maxExoLead() const
+{
+  return max(arg1->maxExoLead(), max(arg2->maxExoLead(), arg3->maxExoLead()));
+}
+
 NodeID
 TrinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -2879,11 +3007,20 @@ TrinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vecto
 }
 
 NodeID
-TrinaryOpNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  if (maxExoLead() == 0)
+    return const_cast<TrinaryOpNode *>(this);
+  else
+    return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
+}
+
+NodeID
+TrinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteExoLeadLag(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteExoLeadLag(subst_table, neweqs);
-  NodeID arg3subst = arg3->substituteExoLeadLag(subst_table, neweqs);
+  NodeID arg1subst = arg1->substituteExoLag(subst_table, neweqs);
+  NodeID arg2subst = arg2->substituteExoLag(subst_table, neweqs);
+  NodeID arg3subst = arg3->substituteExoLag(subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -3028,6 +3165,16 @@ UnknownFunctionNode::maxEndoLead() const
   return val;
 }
 
+int
+UnknownFunctionNode::maxExoLead() const
+{
+  int val = 0;
+  for(vector<NodeID>::const_iterator it = arguments.begin();
+      it != arguments.end(); it++)
+    val = max(val, (*it)->maxExoLead());
+  return val;
+}
+
 NodeID
 UnknownFunctionNode::decreaseLeadsLags(int n) const
 {
@@ -3050,8 +3197,15 @@ UnknownFunctionNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table,
 }
 
 NodeID
-UnknownFunctionNode::substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnknownFunctionNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  cerr << "UnknownFunctionNode::substituteExoLead: not implemented!" << endl;
+  exit(EXIT_FAILURE);
+}
+
+NodeID
+UnknownFunctionNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  cerr << "UnknownFunctionNode::substituteExoLeadLag: not implemented!" << endl;
+  cerr << "UnknownFunctionNode::substituteExoLag: not implemented!" << endl;
   exit(EXIT_FAILURE);
 }
diff --git a/ExprNode.hh b/ExprNode.hh
index a18508c0..e6ae9c04 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -229,6 +229,10 @@ public:
   /*! Always returns a non-negative value */
   virtual int maxEndoLead() const = 0;
 
+  //! Returns the maximum lead of exogenous in this expression
+  /*! Always returns a non-negative value */
+  virtual int maxExoLead() const = 0;
+
   //! Returns a new expression where all the leads/lags have been shifted backwards by the same amount
   /*!
     Only acts on endogenous, exogenous, exogenous det
@@ -245,10 +249,21 @@ public:
     If maximum endogenous lead >= 3, this method will also create intermediary auxiliary var, and will add the equations of the form aux1 = aux2(+1) to the substitution table.
     \pre This expression is assumed to have maximum endogenous lead >= 2
     \param[in,out] subst_table The table to which new auxiliary variables and their correspondance will be added
+    \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
     \return The new variable node corresponding to the current expression
   */
   VariableNode *createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 
+  //! Creates auxiliary exo lead variables corresponding to this expression
+  /*! 
+    If maximum exogenous lead >= 2, this method will also create intermediary auxiliary var, and will add the equations of the form aux1 = aux2(+1) to the substitution table.
+    \pre This expression is assumed to have maximum exogenous lead >= 1
+    \param[in,out] subst_table The table to which new auxiliary variables and their correspondance will be added
+    \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
+    \return The new variable node corresponding to the current expression
+  */
+  VariableNode *createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+
   //! Constructs a new expression where sub-expressions with max endo lead >= 2 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.
@@ -269,12 +284,18 @@ public:
   */
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
-  //! Constructs a new expression where endo variables with max exo lag >= 2 or max exo lag >= 2 have been replaced by auxiliary variables
+  //! Constructs a new expression where exogenous variables with a lead 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.
+    \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
+  */
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) 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.
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -306,10 +327,12 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 //! Symbol or variable node
@@ -341,10 +364,12 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 //! Unary operator node
@@ -382,12 +407,14 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another UnaryOpNode with the same opcode, but with a possibly different datatree and argument
   NodeID buildSimilarUnaryOpNode(NodeID alt_arg, DataTree &alt_datatree) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 //! Binary operator node
@@ -430,12 +457,14 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another BinaryOpNode with the same opcode, but with a possibly different datatree and arguments
   NodeID buildSimilarBinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, DataTree &alt_datatree) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 //! Trinary operator node
@@ -472,12 +501,14 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another TrinaryOpNode with the same opcode, but with a possibly different datatree and arguments
   NodeID buildSimilarTrinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, NodeID alt_arg3, DataTree &alt_datatree) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 //! Unknown function node
@@ -508,10 +539,12 @@ public:
   virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
   virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
   virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
   virtual NodeID decreaseLeadsLags(int n) const;
   virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLeadLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
 };
 
 #endif
diff --git a/ModFile.cc b/ModFile.cc
index 28660a72..4bf811cd 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -142,8 +142,9 @@ ModFile::transformPass()
       || mod_file_struct.ramsey_policy_present)
     {
       dynamic_model.substituteEndoLeadGreaterThanTwo();
+      dynamic_model.substituteExoLead();
       dynamic_model.substituteEndoLagGreaterThanTwo();
-      dynamic_model.substituteExoLeadLag();
+      dynamic_model.substituteExoLag();
     }
 
   // Freeze the symbol table
diff --git a/SymbolTable.cc b/SymbolTable.cc
index ea72365e..8fc5195d 100644
--- a/SymbolTable.cc
+++ b/SymbolTable.cc
@@ -209,9 +209,9 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
       switch(aux_vars[i].type)
         {
         case avEndoLead:
+        case avExoLead:
           break;
         case avEndoLag:
-        case avExoLead:
         case avExoLag:
           output << "M_.aux_vars(" << i+1 << ").orig_endo_index = " << getTypeSpecificID(aux_vars[i].orig_symb_id)+1 << ";" << endl
                  << "M_.aux_vars(" << i+1 << ").orig_lag = " << aux_vars[i].orig_lag << ";" << endl;
@@ -221,10 +221,14 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
 }
 
 int
-SymbolTable::addEndoLeadAuxiliaryVar(int index) throw (FrozenException)
+SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenException)
 {
   ostringstream varname;
-  varname << "AUX_ENDO_LEAD_" << index;
+  if (endo)
+    varname << "AUX_ENDO_LEAD_";
+  else
+    varname << "AUX_EXO_LEAD_";
+  varname << index;
   int symb_id;
   try
     {
@@ -238,31 +242,21 @@ SymbolTable::addEndoLeadAuxiliaryVar(int index) throw (FrozenException)
 
   AuxVarInfo avi;
   avi.symb_id = symb_id;
-  avi.type = avEndoLead;
+  avi.type = (endo ? avEndoLead : avExoLead);
   aux_vars.push_back(avi);
 
   return symb_id;
 }
 
 int
-SymbolTable::addLeadLagAuxiliaryVarInternal(aux_var_t type, int orig_symb_id, int orig_lag) throw (FrozenException)
+SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lag) throw (FrozenException)
 {
   ostringstream varname;
-  switch(type)
-    {
-    case avEndoLag:
-      varname << "AUX_ENDO_LAG_" << orig_symb_id << "_" << -orig_lag;
-      break;
-    case avExoLead:
-      varname << "AUX_EXO_LEAD_" << orig_symb_id << "_" << orig_lag;
-      break;
-    case avExoLag:
-      varname << "AUX_EXO_LAG_" << orig_symb_id << "_" << -orig_lag;
-      break;
-    default:
-      cerr << "Impossible case!" << endl;
-      exit(EXIT_FAILURE);
-    }
+  if (endo)
+    varname << "AUX_ENDO_LAG_";
+  else
+    varname << "AUX_EXO_LAG_";
+  varname << orig_symb_id << "_" << -orig_lag;
 
   int symb_id;
   try
@@ -277,7 +271,7 @@ SymbolTable::addLeadLagAuxiliaryVarInternal(aux_var_t type, int orig_symb_id, in
 
   AuxVarInfo avi;
   avi.symb_id = symb_id;
-  avi.type = type;
+  avi.type = (endo ? avEndoLag : avExoLag);
   avi.orig_symb_id = orig_symb_id;
   avi.orig_lag = orig_lag;
   aux_vars.push_back(avi);
@@ -285,20 +279,26 @@ SymbolTable::addLeadLagAuxiliaryVarInternal(aux_var_t type, int orig_symb_id, in
   return symb_id;
 }
 
+int
+SymbolTable::addEndoLeadAuxiliaryVar(int index) throw (FrozenException)
+{
+  return addLeadAuxiliaryVarInternal(true, index);
+}
+
 int
 SymbolTable::addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException)
 {
-  return addLeadLagAuxiliaryVarInternal(avEndoLag, orig_symb_id, orig_lag);
+  return addLagAuxiliaryVarInternal(true, orig_symb_id, orig_lag);
 }
 
 int
-SymbolTable::addExoLeadAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException)
+SymbolTable::addExoLeadAuxiliaryVar(int index) throw (FrozenException)
 {
-  return addLeadLagAuxiliaryVarInternal(avExoLead, orig_symb_id, orig_lag);
+  return addLeadAuxiliaryVarInternal(false, index);
 }
 
 int
 SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException)
 {
-  return addLeadLagAuxiliaryVarInternal(avExoLag, orig_symb_id, orig_lag);
+  return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lag);
 }
diff --git a/SymbolTable.hh b/SymbolTable.hh
index 071ec3a5..89b11cac 100644
--- a/SymbolTable.hh
+++ b/SymbolTable.hh
@@ -137,7 +137,10 @@ public:
   };
 
 private:
-  int addLeadLagAuxiliaryVarInternal(aux_var_t type, int orig_symb_id, int orig_lag) throw (FrozenException);
+  //! Factorized code for adding aux lag variables
+  int addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lag) throw (FrozenException);
+  //! Factorized code for adding aux lead variables
+  int addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenException);
 
 public:
   //! Add a symbol
@@ -159,10 +162,9 @@ public:
   int addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException);
   //! Adds an auxiliary variable for endogenous with lead >= 1
   /*!
-    \param[in] orig_symb_id symbol ID of the exogenous declared by the user that this new variable will represent
-    \param[in] orig_lag lag value such that this new variable will be equivalent to orig_symb_id(orig_lag)
+    \param[in] index Used to construct the variable name
     \return the symbol ID of the new symbol */
-  int addExoLeadAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException);
+  int addExoLeadAuxiliaryVar(int index) throw (FrozenException);
   //! Adds an auxiliary variable for exogenous with lag >= 1
   /*!
     \param[in] orig_symb_id symbol ID of the exogenous declared by the user that this new variable will represent
-- 
GitLab