diff --git a/DynamicModel.cc b/DynamicModel.cc
index 054482adc0a247ee91002af316c8c028a9c4a23b..67b28c81542501be1e901e96b5f4cd66dae736ee 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -4750,10 +4750,37 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
 }
 
 void
-DynamicModel::substituteAdlAndDiff()
+DynamicModel::substituteAdl()
 {
   for (int i = 0; i < (int) equations.size(); i++)
-    equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdlAndDiff());
+    equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdl());
+}
+
+void
+DynamicModel::substituteDiff()
+{
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode *> neweqs;
+
+  // Substitute in model local variables
+  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);
+
+  // Substitute in equations
+  for (int i = 0; i < (int) equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteDiff(subst_table, neweqs));
+      assert(substeq != NULL);
+      equations[i] = substeq;
+    }
+
+  // Add new equations
+  for (int i = 0; i < (int) neweqs.size(); i++)
+    addEquation(neweqs[i], -1);
+
+  if (subst_table.size() > 0)
+    cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
 }
 
 void
diff --git a/DynamicModel.hh b/DynamicModel.hh
index edd506d25fe71f539de65cf2f1bf2849eb86971a..0d5c94b7e822f1be9bf9c92ef334b529e9ddbdeb 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -366,8 +366,11 @@ public:
   //! Transforms the model by removing trends specified by the user
   void detrendEquations();
 
-  //! Substitutes adl and diff operators
-  void substituteAdlAndDiff();
+  //! Substitutes adl operator
+  void substituteAdl();
+
+  //! Substitutes diff operator
+  void substituteDiff();
 
   //! Fill var_expectation_functions_to_write
   void fillVarExpectationFunctionsToWrite();
diff --git a/ExprNode.cc b/ExprNode.cc
index fbdf400da0b4305f5bb443fd7ceeba96d9f6637e..2b149036d125cd0a7a7480efaa05f1d1655585dd 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2017 Dynare Team
+ * Copyright (C) 2007-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -482,7 +482,13 @@ NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-NumConstNode::substituteAdlAndDiff() const
+NumConstNode::substituteAdl() const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+expr_t
+NumConstNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1244,7 +1250,13 @@ VariableNode::maxLead() const
 }
 
 expr_t
-VariableNode::substituteAdlAndDiff() const
+VariableNode::substituteAdl() const
+{
+  return const_cast<VariableNode *>(this);
+}
+
+expr_t
+VariableNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VariableNode *>(this);
 }
@@ -2709,22 +2721,15 @@ UnaryOpNode::maxLead() const
 }
 
 expr_t
-UnaryOpNode::substituteAdlAndDiff() const
+UnaryOpNode::substituteAdl() const
 {
-  if (op_code != oDiff && op_code != oAdl)
+  if (op_code != oAdl)
     {
-      expr_t argsubst = arg->substituteAdlAndDiff();
+      expr_t argsubst = arg->substituteAdl();
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 
-  if (op_code == oDiff)
-    {
-      expr_t argsubst = arg->substituteAdlAndDiff();
-      return datatree.AddMinus(argsubst,
-                               argsubst->decreaseLeadsLags(1));
-    }
-
-  expr_t arg1subst = arg->substituteAdlAndDiff();
+  expr_t arg1subst = arg->substituteAdl();
   expr_t retval;
   ostringstream inttostr;
   if (adl_param_lag >= 0)
@@ -2765,6 +2770,34 @@ UnaryOpNode::substituteAdlAndDiff() const
   return retval;
 }
 
+expr_t
+UnaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  if (op_code != oDiff)
+    {
+      expr_t argsubst = arg->substituteDiff(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);
+
+  int 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,
+                                                                  datatree.AddMinus(argsubst, argsubst->decreaseLeadsLags(1)))));
+
+  assert(dynamic_cast<VariableNode *>(newAuxVar) != NULL);
+  subst_table[this] = dynamic_cast<VariableNode *>(newAuxVar);
+  return newAuxVar;
+}
+
 expr_t
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -4361,10 +4394,18 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-BinaryOpNode::substituteAdlAndDiff() const
+BinaryOpNode::substituteAdl() const
+{
+  expr_t arg1subst = arg1->substituteAdl();
+  expr_t arg2subst = arg2->substituteAdl();
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
+expr_t
+BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteAdlAndDiff();
-  expr_t arg2subst = arg2->substituteAdlAndDiff();
+  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -5108,11 +5149,20 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
 }
 
 expr_t
-TrinaryOpNode::substituteAdlAndDiff() const
+TrinaryOpNode::substituteAdl() const
+{
+  expr_t arg1subst = arg1->substituteAdl();
+  expr_t arg2subst = arg2->substituteAdl();
+  expr_t arg3subst = arg3->substituteAdl();
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
+expr_t
+TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->substituteAdlAndDiff();
-  expr_t arg2subst = arg2->substituteAdlAndDiff();
-  expr_t arg3subst = arg3->substituteAdlAndDiff();
+  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiff(subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -5420,11 +5470,20 @@ AbstractExternalFunctionNode::substituteExpectation(subst_table_t &subst_table,
 }
 
 expr_t
-AbstractExternalFunctionNode::substituteAdlAndDiff() const
+AbstractExternalFunctionNode::substituteAdl() const
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteAdlAndDiff());
+    arguments_subst.push_back((*it)->substituteAdl());
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
+expr_t
+AbstractExternalFunctionNode::substituteDiff(subst_table_t &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));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -6892,7 +6951,13 @@ VarExpectationNode::substituteExpectation(subst_table_t &subst_table, vector<Bin
 }
 
 expr_t
-VarExpectationNode::substituteAdlAndDiff() const
+VarExpectationNode::substituteAdl() const
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
+expr_t
+VarExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VarExpectationNode *>(this);
 }
diff --git a/ExprNode.hh b/ExprNode.hh
index 9fa28e05c513dd4f363f1cf7fba59b6966a9984e..e2e20d6ae8b2d847a9664dbb1a9c4503055c2131 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2017 Dynare Team
+ * Copyright (C) 2007-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -457,8 +457,11 @@ enum ExprNodeOutputType
       */
       virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const = 0;
 
-      //! Substitute adl and diff operators
-      virtual expr_t substituteAdlAndDiff() const = 0;
+      //! Substitute adl operator
+      virtual expr_t substituteAdl() const = 0;
+
+      //! Substitute diff operator
+      virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Add ExprNodes to the provided datatree
       virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
@@ -530,7 +533,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -605,7 +609,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -703,7 +708,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -813,7 +819,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -903,7 +910,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
@@ -993,7 +1001,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1174,7 +1183,8 @@ public:
   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 substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   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,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
diff --git a/ModFile.cc b/ModFile.cc
index 1de813fcde3c47fc60359e16bfad9dd21930651c..f81cb0963af1896aa687046adfbd6ecf1065a9a3 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -346,7 +346,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   // - except adl and diff which we always want expanded
-  dynamic_model.substituteAdlAndDiff();
+  dynamic_model.substituteAdl();
   dynamic_model.setLeadsLagsOrig();
   dynamic_model.cloneDynamic(original_model);
 
@@ -381,6 +381,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   if (symbol_table.predeterminedNbr() > 0)
     dynamic_model.transformPredeterminedVariables();
 
+  // Create auxiliary variable and equations for Diff operator
+  dynamic_model.substituteDiff();
+
   // Create auxiliary vars for Expectation operator
   dynamic_model.substituteExpectation(mod_file_struct.partial_information);
 
diff --git a/SymbolTable.cc b/SymbolTable.cc
index dde3d141bc047ed7d1202920b61538721ab0cd92..4b18946c0faa0fd21ca3193a6c0ff6c04327eddd 100644
--- a/SymbolTable.cc
+++ b/SymbolTable.cc
@@ -354,6 +354,7 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
           {
           case avEndoLead:
           case avExoLead:
+          case avDiff:
             break;
           case avEndoLag:
           case avExoLag:
@@ -475,6 +476,7 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
             case avExpectation:
             case avMultiplier:
             case avDiffForward:
+            case avDiff:
               break;
             case avEndoLag:
             case avExoLag:
@@ -568,6 +570,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
         case avExpectation:
         case avMultiplier:
         case avDiffForward:
+        case avDiff:
           break;
         case avEndoLag:
         case avExoLag:
@@ -691,6 +694,29 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
   return symb_id;
 }
 
+int
+SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_DIFF_" << index;
+
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.push_back(AuxVarInfo(symb_id, avDiff, 0, 0, 0, 0, expr_arg));
+
+  return symb_id;
+}
+
 int
 SymbolTable::addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException)
 {
@@ -1001,6 +1027,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
             {
             case avEndoLead:
             case avExoLead:
+            case avDiff:
               break;
             case avEndoLag:
             case avExoLag:
diff --git a/SymbolTable.hh b/SymbolTable.hh
index 76aa3896141782b1533d21bad56da4adaa58953b..53ba84d4990c0096f48edb24433c637f0f8ce223 100644
--- a/SymbolTable.hh
+++ b/SymbolTable.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -43,7 +43,8 @@ enum aux_var_t
     avExpectation = 4,    //!< Substitute for Expectation Operator
     avDiffForward = 5,    //!< Substitute for the differentiate of a forward variable
     avMultiplier = 6,     //!< Multipliers for FOC of Ramsey Problem
-    avVarModel = 7        //!< Variable for var_model with order > abs(min_lag()) present in model
+    avVarModel = 7,       //!< Variable for var_model with order > abs(min_lag()) present in model
+    avDiff = 8            //!< Variable for Diff operator
   };
 
 //! Information on some auxiliary variables
@@ -283,6 +284,8 @@ public:
   //! Adds an auxiliary variable when var_model is used with an order that is greater in absolute value
   //! than the largest lag present in the model.
   int addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException);
+  //! Adds an auxiliary variable when the diff operator is encountered
+  int addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const