diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 921bbc364d590b4493396a1623b50aa6aea36dc2..c4812ebde64ee1aeadd09b74ab2fbf9cb20c4934 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -1367,14 +1367,8 @@ VariableNode::undiff() const
 void
 VariableNode::VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs, int &max_lag) const
 {
-  for (set<expr_t>::const_iterator it = static_lhs.begin();
-       it != static_lhs.end(); it++)
-    if (*it == this->toStatic(static_datatree))
-      {
-        if (-lag > max_lag)
-          max_lag = -lag;
-        return;
-      }
+  if (-lag > max_lag)
+    max_lag = -lag;
 }
 
 int
@@ -3011,6 +3005,74 @@ UnaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table)
     diff_table[sthis][arg_max_lag] = const_cast<UnaryOpNode *>(this);
 }
 
+void
+UnaryOpNode::getDiffArgUnaryOperatorIfAny(string &op_handle) const
+{
+  switch (op_code)
+    {
+    case oExp:
+      op_handle = "@exp";
+      break;
+    case oLog:
+      op_handle = "@log";
+      break;
+    case oLog10:
+      op_handle = "@log10";
+      break;
+    case oCos:
+      op_handle = "@cos";
+      break;
+    case oSin:
+      op_handle = "@sin";
+      break;
+    case oTan:
+      op_handle = "@tan";
+      break;
+    case oAcos:
+      op_handle = "@acos";
+      break;
+    case oAsin:
+      op_handle = "@asin";
+      break;
+    case oAtan:
+      op_handle = "@atan";
+      break;
+    case oCosh:
+      op_handle = "@cosh";
+      break;
+    case oSinh:
+      op_handle = "@sinh";
+      break;
+    case oTanh:
+      op_handle = "@tanh";
+      break;
+    case oAcosh:
+      op_handle = "@acosh";
+      break;
+    case oAsinh:
+      op_handle = "@asinh";
+      break;
+    case oAtanh:
+      op_handle = "@atanh";
+      break;
+    case oSqrt:
+      op_handle = "@sqrt";
+      break;
+    case oAbs:
+      op_handle = "@abs";
+      break;
+    case oSign:
+      op_handle = "@sign";
+      break;
+    case oErf:
+      op_handle = "@erf";
+      break;
+    default:
+      op_handle = "";
+      break;
+    }
+}
+
 expr_t
 UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                             vector<BinaryOpNode *> &neweqs) const
@@ -3047,8 +3109,29 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
           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);
-
+            {
+              UnaryOpNode *diffarg = dynamic_cast<UnaryOpNode *>(argsubst);
+              if (diffarg != NULL)
+                {
+                  string op;
+                  diffarg->getDiffArgUnaryOperatorIfAny(op);
+                  VariableNode *vnarg = dynamic_cast<VariableNode *>(diffarg->get_arg());
+                  if (vnarg != NULL)
+                    symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst,
+                                                                        vnarg->get_symb_id(), vnarg->get_lag(), op);
+                  else
+                    {
+                      // The case where we have diff(log(exp(x))) for example
+                      cerr << "diffs of nested non-diff expressions are not yet supported" << endl;
+                      exit(EXIT_FAILURE);
+                    }
+                }
+              else
+                {
+                  cerr << "diffs of non unary expressions are not yet supported" << endl;
+                  exit(EXIT_FAILURE);
+                }
+            }
           // make originating aux var & equation
           last_arg_max_lag = rit->first;
           last_aux_var = datatree.AddVariable(symb_id, 0);
@@ -3064,18 +3147,12 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
           VariableNode *new_aux_var = NULL;
           for (int i = last_arg_max_lag; i > rit->first; i--)
             {
-              if (vn == NULL)
-                if (i == last_arg_max_lag)
-                  symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
-                else
-                  symb_id = datatree.symbol_table.addDiffAuxiliaryVar(new_aux_var->idx, new_aux_var);
+              if (i == last_arg_max_lag)
+                symb_id = datatree.symbol_table.addDiffLagAuxiliaryVar(argsubst->idx, argsubst,
+                                                                       last_aux_var->get_symb_id(), last_aux_var->get_lag());
               else
-                if (i == last_arg_max_lag)
-                  symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst,
-                                                                      vn->get_symb_id(), i - 1);
-                else
-                  symb_id = datatree.symbol_table.addDiffAuxiliaryVar(new_aux_var->idx, new_aux_var,
-                                                                      vn->get_symb_id(), i - 1);
+                symb_id = datatree.symbol_table.addDiffLagAuxiliaryVar(new_aux_var->idx, new_aux_var,
+                                                                       last_aux_var->get_symb_id(), last_aux_var->get_lag());
 
               new_aux_var = datatree.AddVariable(symb_id, 0);
               neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 7d77ddb1776ac469990519aa019c2acf69291462..59f4ebcb099e96072d2efd542b6d21c3fe17318d 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -59,7 +59,7 @@ typedef map<int, double> eval_context_t;
 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
-//! diff_table[static_expr_t][lag] -> dynamic_expr_t
+//! diff_table[static_expr_t][lag] -> [dynamic_expr_t]
 typedef map<expr_t, map<int, expr_t> > diff_table_t;
 
 //! Possible types of output when writing ExprNode(s)
@@ -776,6 +776,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
   virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  void getDiffArgUnaryOperatorIfAny(string &op_handle) const;
   virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 6a037ac1735990b474df27cba8f99feee85396d8..65e7fbfe98c9e115d4920a44977ae741430a8b37 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -518,6 +518,28 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 
   dynamic_model.combineDiffAuxEquations();
 
+  for (vector<Statement *>::const_iterator it = statements.begin();
+       it != statements.end(); it++)
+    {
+      VarModelStatement *vms = dynamic_cast<VarModelStatement *>(*it);
+      if (vms != NULL)
+        {
+          string var_model_name;
+          vms->getVarModelInfo(var_model_name, var_model_info_var_expectation, var_model_eq_tags);
+          vector<expr_t> lhs_expr_t;
+          vector<int> lhs, eqnumber, orig_diff_var;
+          vector<set<pair<int, int> > > rhs;
+          vector<bool> nonstationary, diff;
+          vector<string> eqtags = var_model_eq_tags[var_model_name];
+          dynamic_model.getVarModelVariablesFromEqTags(eqtags,
+                                                       eqnumber, lhs, lhs_expr_t, rhs, nonstationary);
+          int max_lag = original_model.getVarMaxLag(diff_static_model, eqnumber);
+          original_model.getVarLhsDiffAndInfo(eqnumber, diff, orig_diff_var);
+          vms->fillVarModelInfoFromEquations(eqnumber, lhs, rhs, nonstationary,
+                                             diff, orig_diff_var, max_lag);
+        }
+    }
+
   if (differentiate_forward_vars)
     dynamic_model.differentiateForwardVars(differentiate_forward_vars_subset);
 
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 2742e21bb1243a38865562cd0be6ab068932caf3..2ec4060a3096c38e5474b0e77f838d468452afc1 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -25,6 +25,20 @@
 
 #include "SymbolTable.hh"
 
+AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg,
+                       int equation_number_for_multiplier_arg, int information_set_arg,
+                       expr_t expr_node_arg, string &unary_op_handle_arg) :
+  symb_id(symb_id_arg),
+  type(type_arg),
+  orig_symb_id(orig_symb_id_arg),
+  orig_lead_lag(orig_lead_lag_arg),
+  equation_number_for_multiplier(equation_number_for_multiplier_arg),
+  information_set(information_set_arg),
+  expr_node(expr_node_arg),
+  unary_op_handle(unary_op_handle_arg)
+{
+}
+
 AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg,
                        int equation_number_for_multiplier_arg, int information_set_arg,
                        expr_t expr_node_arg) :
@@ -34,7 +48,8 @@ AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg
   orig_lead_lag(orig_lead_lag_arg),
   equation_number_for_multiplier(equation_number_for_multiplier_arg),
   information_set(information_set_arg),
-  expr_node(expr_node_arg)
+  expr_node(expr_node_arg),
+  unary_op_handle("")
 {
 }
 
@@ -375,9 +390,16 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             output << ")';" << endl;
             break;
           case avDiff:
+          case avDiffLag:
             if (aux_vars[i].get_orig_symb_id() >= 0)
               output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl
                      << "M_.aux_vars(" << i+1 << ").orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
+
+            output << "M_.aux_vars(" << i+1 << ").unary_op_handle = ";
+            if (!aux_vars[i].get_unary_op_handle().empty())
+              output << aux_vars[i].get_unary_op_handle() << ";" << endl;
+            else
+              output << "'';" << endl;
             break;
           }
       }
@@ -488,6 +510,7 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
                      << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
               break;
             case avDiff:
+            case avDiffLag:
               if (aux_vars[i].get_orig_symb_id() >= 0)
                 output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                        << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -586,6 +609,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
                  << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
           break;
         case avDiff:
+        case avDiffLag:
           if (aux_vars[i].get_orig_symb_id() >= 0)
             output << "av" << i << ".orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                    << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -706,6 +730,29 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
   return symb_id;
 }
 
+int
+SymbolTable::addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_DIFF_LAG_" << 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, avDiffLag, orig_symb_id, orig_lag, 0, 0, expr_arg));
+
+  return symb_id;
+}
+
 int
 SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException)
 {
@@ -729,6 +776,29 @@ SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, i
   return symb_id;
 }
 
+int
+SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag, string &unary_op_handle) 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, orig_symb_id, orig_lag, 0, 0, expr_arg, unary_op_handle));
+
+  return symb_id;
+}
+
 int
 SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
 {
@@ -1067,9 +1137,10 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
                      << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
               break;
             case avDiff:
+            case avDiffLag:
               if (aux_vars[i].get_orig_symb_id() >= 0)
                 output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
-                       << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
+                       << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN," << aux_vars[i].get_unary_op_handle();
               break;
             case avMultiplier:
               output << "NaN, NaN, " << aux_vars[i].get_equation_number_for_multiplier() + 1
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 7a0d89a45ff7606ebdc215b22b2a07856e3801da..58dee978321161bfe82f8ec728908f29fd250621 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -44,7 +44,8 @@ enum aux_var_t
     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
-    avDiff = 8            //!< Variable for Diff operator
+    avDiff = 8,           //!< Variable for Diff operator
+    avDiffLag = 9         //!< Variable for timing between Diff operators
   };
 
 //! Information on some auxiliary variables
@@ -58,8 +59,10 @@ private:
   int equation_number_for_multiplier; //!< Stores the original constraint equation number associated with this aux var. Only used for avMultiplier.
   int information_set; //! Argument of expectation operator. Only used for avExpectation.
   expr_t expr_node; //! Auxiliary variable definition
+  string unary_op_handle; //!Unary op potentially opplied to aux vars of type avDiff
 public:
   AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, int equation_number_for_multiplier_arg, int information_set_arg, expr_t expr_node_arg);
+  AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, int equation_number_for_multiplier_arg, int information_set_arg, expr_t expr_node_arg, string &unary_op_handle);
   int
   get_symb_id() const
   {
@@ -95,6 +98,11 @@ public:
   {
     return expr_node;
   };
+  string
+  get_unary_op_handle() const
+  {
+    return unary_op_handle;
+  }
 };
 
 //! Stores the symbol table
@@ -289,6 +297,9 @@ public:
   //! Adds an auxiliary variable when the diff operator is encountered
   int addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException);
   int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
+  int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag, string &unary_op_handle) throw (FrozenException);
+  //! Takes care of timing between diff statements
+  int addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const