From 1e20717f5874d7b488cc9dda096659aa47ac94ce Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Wed, 28 Feb 2018 17:33:00 +0100
Subject: [PATCH] fix bug in substitution of diff operator

---
 src/DynamicModel.cc |  13 +++--
 src/DynamicModel.hh |   2 +-
 src/ExprNode.cc     | 118 +++++++++++++++++++++++++++++++++-----------
 src/ExprNode.hh     |  36 +++++++++-----
 src/ModFile.cc      |   3 +-
 src/ModFile.hh      |   2 +
 6 files changed, 126 insertions(+), 48 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 10024fc6..99fc9cf1 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4897,20 +4897,19 @@ DynamicModel::substituteAdl()
 }
 
 void
-DynamicModel::substituteDiff()
+DynamicModel::substituteDiff(StaticModel &static_model)
 {
-  ExprNode::subst_table_t subst_table;
-  vector<BinaryOpNode *> neweqs;
-
   // Substitute in model local variables
+  vector<BinaryOpNode *> neweqs;
+  diff_subst_table_t diff_subst_table;
   for (map<int, expr_t>::iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
-    it->second = it->second->substituteDiff(subst_table, neweqs);
+    it->second = it->second->substituteDiff(static_model, diff_subst_table, neweqs);
 
   // Substitute in equations
   for (int i = 0; i < (int) equations.size(); i++)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteDiff(subst_table, neweqs));
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteDiff(static_model, diff_subst_table, neweqs));
       assert(substeq != NULL);
       equations[i] = substeq;
     }
@@ -4919,7 +4918,7 @@ DynamicModel::substituteDiff()
   for (int i = 0; i < (int) neweqs.size(); i++)
     addEquation(neweqs[i], -1);
 
-  if (subst_table.size() > 0)
+  if (diff_subst_table.size() > 0)
     cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
 }
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index afcd72b9..dda2a2aa 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -392,7 +392,7 @@ public:
   void substituteAdl();
 
   //! Substitutes diff operator
-  void substituteDiff();
+  void substituteDiff(StaticModel &static_model);
 
   //! Fill var_expectation_functions_to_write
   void fillVarExpectationFunctionsToWrite();
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 43095bcf..3ae0a2a7 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -451,6 +451,12 @@ NumConstNode::maxLead() const
   return 0;
 }
 
+int
+NumConstNode::maxLag() const
+{
+  return 0;
+}
+
 expr_t
 NumConstNode::decreaseLeadsLags(int n) const
 {
@@ -500,7 +506,7 @@ NumConstNode::substituteAdl() const
 }
 
 expr_t
-NumConstNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1282,6 +1288,22 @@ VariableNode::maxLead() const
     }
 }
 
+int
+VariableNode::maxLag() const
+{
+  switch (type)
+    {
+    case eEndogenous:
+      return -lag;
+    case eExogenous:
+      return -lag;
+    case eModelLocalVariable:
+      return datatree.local_variables_table[symb_id]->maxLag();
+    default:
+      return 0;
+    }
+}
+
 expr_t
 VariableNode::substituteAdl() const
 {
@@ -1289,7 +1311,7 @@ VariableNode::substituteAdl() const
 }
 
 expr_t
-VariableNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VariableNode *>(this);
 }
@@ -2788,6 +2810,12 @@ UnaryOpNode::maxLead() const
   return arg->maxLead();
 }
 
+int
+UnaryOpNode::maxLag() const
+{
+  return arg->maxLag();
+}
+
 expr_t
 UnaryOpNode::substituteAdl() const
 {
@@ -2830,37 +2858,37 @@ UnaryOpNode::isDiffPresent() const
 }
 
 expr_t
-UnaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (op_code != oDiff)
     {
-      expr_t argsubst = arg->substituteDiff(subst_table, neweqs);
+      expr_t argsubst = arg->substituteDiff(static_datatree, diff_subst_table, neweqs);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 
-  subst_table_t::iterator it = subst_table.find(const_cast<UnaryOpNode *>(this));
-  if (it != subst_table.end())
-    return const_cast<VariableNode *>(it->second);
-
-  expr_t argsubst = arg->substituteDiff(subst_table, neweqs);
-  assert(argsubst != NULL);
+  expr_t sarg = arg->toStatic(static_datatree);
+  int arg_max_lag = -arg->maxLag();
+  diff_subst_table_t::iterator it = diff_subst_table.find(sarg);
+  if (it != diff_subst_table.end())
+    if (it->second.first.first == arg)
+      return const_cast<VariableNode *>(it->second.second);
+    else
+      return it->second.second->decreaseLeadsLags(it->second.first.second - arg_max_lag);
 
   int symb_id = -1;
+  expr_t argsubst = arg->substituteDiff(static_datatree, diff_subst_table, neweqs);
   VariableNode *vn = dynamic_cast<VariableNode *>(argsubst);
   if (vn != NULL)
     symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
   else
     symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
 
-  expr_t newAuxVar = datatree.AddVariable(symb_id, 0);
-
   //AUX_DIFF_IDX = argsubst - argsubst(-1)
-  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(newAuxVar,
+  VariableNode *new_aux_var = datatree.AddVariable(symb_id, 0);
+  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
                                                                   datatree.AddMinus(argsubst, argsubst->decreaseLeadsLags(1)))));
-
-  assert(dynamic_cast<VariableNode *>(newAuxVar) != NULL);
-  subst_table[this] = dynamic_cast<VariableNode *>(newAuxVar);
-  return newAuxVar;
+  diff_subst_table[sarg] = make_pair(make_pair(arg, arg_max_lag), new_aux_var);
+  return new_aux_var;
 }
 
 expr_t
@@ -4355,6 +4383,12 @@ BinaryOpNode::maxLead() const
   return max(arg1->maxLead(), arg2->maxLead());
 }
 
+int
+BinaryOpNode::maxLag() const
+{
+  return max(arg1->maxLag(), arg2->maxLag());
+}
+
 expr_t
 BinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -4492,10 +4526,10 @@ BinaryOpNode::substituteAdl() const
 }
 
 expr_t
-BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
-  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -5260,6 +5294,12 @@ TrinaryOpNode::maxLead() const
   return max(arg1->maxLead(), max(arg2->maxLead(), arg3->maxLead()));
 }
 
+int
+TrinaryOpNode::maxLag() const
+{
+  return max(arg1->maxLag(), max(arg2->maxLag(), arg3->maxLag()));
+}
+
 expr_t
 TrinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -5347,11 +5387,11 @@ TrinaryOpNode::substituteAdl() const
 }
 
 expr_t
-TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
-  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
-  expr_t arg3subst = arg3->substituteDiff(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiff(static_datatree, diff_subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -5634,6 +5674,16 @@ AbstractExternalFunctionNode::maxLead() const
   return val;
 }
 
+int
+AbstractExternalFunctionNode::maxLag() const
+{
+  int val = 0;
+  for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+    val = max(val, (*it)->maxLag());
+  return val;
+}
+
 expr_t
 AbstractExternalFunctionNode::decreaseLeadsLags(int n) const
 {
@@ -5707,11 +5757,11 @@ AbstractExternalFunctionNode::substituteAdl() const
 }
 
 expr_t
-AbstractExternalFunctionNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteDiff(subst_table, neweqs));
+    arguments_subst.push_back((*it)->substituteDiff(static_datatree, diff_subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7112,6 +7162,12 @@ VarExpectationNode::maxLead() const
   return 0;
 }
 
+int
+VarExpectationNode::maxLag() const
+{
+  return 0;
+}
+
 expr_t
 VarExpectationNode::decreaseLeadsLags(int n) const
 {
@@ -7230,7 +7286,7 @@ VarExpectationNode::substituteAdl() const
 }
 
 expr_t
-VarExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -7530,6 +7586,12 @@ PacExpectationNode::maxLead() const
   return 0;
 }
 
+int
+PacExpectationNode::maxLag() const
+{
+  return 0;
+}
+
 expr_t
 PacExpectationNode::decreaseLeadsLags(int n) const
 {
@@ -7647,7 +7709,7 @@ PacExpectationNode::substituteAdl() const
 }
 
 expr_t
-PacExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<PacExpectationNode *>(this);
 }
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 7fe371cc..c3bbcb14 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -57,6 +57,9 @@ typedef map<int, double> eval_context_t;
 //! Type for tracking first/second derivative functions that have already been written as temporary terms
 typedef map<pair<int, vector<expr_t> >, int> deriv_node_temp_terms_t;
 
+//! Type for the substitution map used in the process of substitutitng diff expressions
+typedef map<expr_t, pair<pair<expr_t, int>, VariableNode *> > diff_subst_table_t;
+
 //! Possible types of output when writing ExprNode(s)
 enum ExprNodeOutputType
   {
@@ -339,6 +342,10 @@ class ExprNode
       /*! A negative value means that the expression contains only lagged variables */
       virtual int maxLead() const = 0;
 
+      //! Returns the relative period of the most backward term in this expression
+      /*! A negative value means that the expression contains only leaded variables */
+      virtual int maxLag() 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
@@ -350,8 +357,7 @@ class ExprNode
       //! Type for the substitution map used in the process of creating auxiliary vars for leads >= 2
       typedef map<const ExprNode *, const VariableNode *> subst_table_t;
 
-      //! Type for the substitution map used in the process of substituting adl/diff expressions
-      //      typedef map<const ExprNode *, const UnaryOpNode *> subst_table_diff_t;
+      //! Type for the substitution map used in the process of substituting adl expressions
       typedef map<const ExprNode *, const expr_t> subst_table_adl_t;
 
       //! Creates auxiliary endo lead variables corresponding to this expression
@@ -466,7 +472,7 @@ class ExprNode
       virtual expr_t substituteAdl() const = 0;
 
       //! Substitute diff operator
-      virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+      virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
       virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
@@ -544,6 +550,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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;
@@ -551,7 +558,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -625,6 +632,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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;
@@ -632,7 +640,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -726,6 +734,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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
@@ -735,7 +744,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -845,6 +854,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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
@@ -854,7 +864,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -941,6 +951,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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
@@ -950,7 +961,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1039,6 +1050,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) 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;
@@ -1046,7 +1058,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -1220,6 +1232,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
   virtual void prepareForDerivation();
   virtual expr_t computeDerivative(int deriv_id);
@@ -1233,7 +1246,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
@@ -1295,6 +1308,7 @@ public:
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
   virtual int maxLead() const;
+  virtual int maxLag() const;
   virtual expr_t decreaseLeadsLags(int n) const;
   virtual void prepareForDerivation();
   virtual expr_t computeDerivative(int deriv_id);
@@ -1308,7 +1322,7 @@ public:
   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 substituteAdl() const;
-  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
diff --git a/src/ModFile.cc b/src/ModFile.cc
index eb4a374b..df36512c 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -39,6 +39,7 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
     orig_ramsey_dynamic_model(symbol_table, num_constants, external_functions_table),
     static_model(symbol_table, num_constants, external_functions_table),
     steady_state_model(symbol_table, num_constants, external_functions_table, static_model),
+    diff_static_model(symbol_table, num_constants, external_functions_table),
     linear(false), block(false), byte_code(false), use_dll(false), no_static(false),
     differentiate_forward_vars(false), nonstationary_variables(false),
     param_used_with_lead_lag(false), warnings(warnings_arg)
@@ -362,7 +363,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     }
 
   // Create auxiliary variable and equations for Diff operator
-  dynamic_model.substituteDiff();
+  dynamic_model.substituteDiff(diff_static_model);
 
   // Var Model
   map<string, pair<SymbolList, int> > var_model_info_var_expectation;
diff --git a/src/ModFile.hh b/src/ModFile.hh
index ac678b5b..08ee6c42 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -72,6 +72,8 @@ public:
   StaticModel static_model;
   //! Static model, as declared in the "steady_state_model" block if present
   SteadyStateModel steady_state_model;
+  //! Static model used for mapping arguments of diff operator
+  StaticModel diff_static_model;
   //! Option linear
   bool linear;
 
-- 
GitLab