diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index f43021d7f240baa5a4cda1e7f64431bc641cab66..73d34bacfbdf85a8ab5013fb60c5a04460976261 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5188,6 +5188,44 @@ DynamicModel::substituteAdl()
     equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdl());
 }
 
+void
+DynamicModel::substituteUnaryOps(StaticModel &static_model)
+{
+  diff_table_t nodes;
+  // Find matching unary ops that may be outside of diffs (i.e., those with different lags)
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
+
+  for (int i = 0; i < (int) equations.size(); i++)
+    equations[i]->findUnaryOpNodesForAuxVarCreation(static_model, nodes);
+
+  // Substitute in model local variables
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode *> neweqs;
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second = it->second->substituteUnaryOpNodes(static_model, nodes, subst_table, neweqs);
+
+  // Substitute in equations
+  for (int i = 0; i < (int) equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->
+                                                           substituteUnaryOpNodes(static_model, nodes, 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);
+
+  copy(neweqs.begin(), neweqs.end(), back_inserter(diff_aux_equations));
+
+  if (subst_table.size() > 0)
+    cout << "Substitution of Unary Ops: added " << neweqs.size() << " auxiliary variables and equations." << endl;
+}
+
 void
 DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table)
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index ee9f981c80f6e928deefc84b4fb5e381e76f61e0..726a45a33295dea0e7e29d7d7003f53bd0212d06 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -409,6 +409,9 @@ public:
   //! Substitutes adl operator
   void substituteAdl();
 
+  //! Creates aux vars for certain unary operators: originally implemented for support of VARs
+  void substituteUnaryOps(StaticModel &static_model);
+
   //! Substitutes diff operator
   void substituteDiff(StaticModel &static_model, ExprNode::subst_table_t &diff_subst_table);
 
diff --git a/src/DynareMain.cc b/src/DynareMain.cc
index 30a3aa16d1ddb0214b755fd9882bf56bc69bc55d..00974428e90c73f9d388890811d1d96324ec4628 100644
--- a/src/DynareMain.cc
+++ b/src/DynareMain.cc
@@ -41,7 +41,7 @@ void main2(stringstream &in, string &basename, bool debug, bool clear_all, bool
            bool nograph, bool nointeractive, bool parallel, ConfigFile &config_file,
            WarningConsolidation &warnings_arg, bool nostrict, bool stochastic, bool check_model_changes,
            bool minimal_workspace, bool compute_xrefs, FileOutputType output_mode,
-           LanguageOutputType lang, int params_derivs_order
+           LanguageOutputType lang, int params_derivs_order, bool transform_unary_ops
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
            , bool cygwin, bool msvc, bool mingw
 #endif
@@ -58,7 +58,7 @@ usage()
   cerr << "Dynare usage: dynare mod_file [debug] [noclearall] [onlyclearglobals] [savemacro[=macro_file]] [onlymacro] [nolinemacro] [noemptylinemacro] [notmpterms] [nolog] [warn_uninit]"
        << " [console] [nograph] [nointeractive] [parallel[=cluster_name]] [conffile=parallel_config_path_and_filename] [parallel_slave_open_mode] [parallel_test]"
        << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] [compute_xrefs] [output=dynamic|first|second|third] [language=C|C++|julia]"
-       << " [params_derivs_order=0|1|2]"
+       << " [params_derivs_order=0|1|2] [transform_unary_ops]"
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
        << " [cygwin] [msvc] [mingw]"
 #endif
@@ -114,6 +114,7 @@ main(int argc, char **argv)
   bool check_model_changes = false;
   bool minimal_workspace = false;
   bool compute_xrefs = false;
+  bool transform_unary_ops = false;
   map<string, string> defines;
   vector<string> path;
   FileOutputType output_mode = none;
@@ -210,6 +211,8 @@ main(int argc, char **argv)
         minimal_workspace = true;
       else if (!strcmp(argv[arg], "compute_xrefs"))
         compute_xrefs = true;
+      else if (!strcmp(argv[arg], "transform_unary_ops"))
+        transform_unary_ops = true;
       else if (strlen(argv[arg]) >= 8 && !strncmp(argv[arg], "parallel", 8))
         {
           parallel = true;
@@ -401,7 +404,7 @@ main(int argc, char **argv)
   main2(macro_output, basename, debug, clear_all, clear_global,
         no_tmp_terms, no_log, no_warn, warn_uninit, console, nograph, nointeractive,
         parallel, config_file, warnings, nostrict, stochastic, check_model_changes, minimal_workspace,
-        compute_xrefs, output_mode, language, params_derivs_order
+        compute_xrefs, output_mode, language, params_derivs_order, transform_unary_ops
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
         , cygwin, msvc, mingw
 #endif
diff --git a/src/DynareMain2.cc b/src/DynareMain2.cc
index d61df45b5039433290c4a67cc0e50dffd0c42276..b5ae53a150395cb75b08e4dbf13f084f47f27600 100644
--- a/src/DynareMain2.cc
+++ b/src/DynareMain2.cc
@@ -30,7 +30,7 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
       bool nograph, bool nointeractive, bool parallel, ConfigFile &config_file,
       WarningConsolidation &warnings, bool nostrict, bool stochastic, bool check_model_changes,
       bool minimal_workspace, bool compute_xrefs, FileOutputType output_mode,
-      LanguageOutputType language, int params_derivs_order
+      LanguageOutputType language, int params_derivs_order, bool transform_unary_ops
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
       , bool cygwin, bool msvc, bool mingw
 #endif
@@ -51,7 +51,7 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
     mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput);
 
   // Perform transformations on the model (creation of auxiliary vars and equations)
-  mod_file->transformPass(nostrict, stochastic, compute_xrefs || json == transformpass, nopreprocessoroutput);
+  mod_file->transformPass(nostrict, stochastic, compute_xrefs || json == transformpass, nopreprocessoroutput, transform_unary_ops);
   if (json == transformpass)
     mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput);
 
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index f6018f937b9efb76e313449943337ccba80a206e..c38154e00973d00f16ec6d01b7f70399ed6fc45f 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -538,12 +538,23 @@ NumConstNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table)
 {
 }
 
+void
+NumConstNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+}
+
 expr_t
 NumConstNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 expr_t
 NumConstNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -1413,6 +1424,11 @@ VariableNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table)
 {
 }
 
+void
+VariableNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+}
+
 expr_t
 VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                              vector<BinaryOpNode *> &neweqs) const
@@ -1420,6 +1436,12 @@ VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VariableNode *>(this);
+}
+
 expr_t
 VariableNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -3013,6 +3035,60 @@ UnaryOpNode::isDiffPresent() const
   return arg->isDiffPresent();
 }
 
+bool
+UnaryOpNode::createAuxVarForUnaryOpNodeInDiffOp() const
+{
+  switch (op_code)
+    {
+    case oExp:
+    case oLog:
+    case oLog10:
+    case oCos:
+    case oSin:
+    case oTan:
+    case oAcos:
+    case oAsin:
+    case oAtan:
+    case oCosh:
+    case oSinh:
+    case oTanh:
+    case oAcosh:
+    case oAsinh:
+    case oAtanh:
+    case oSqrt:
+    case oAbs:
+    case oSign:
+    case oErf:
+      return true;
+    default:
+      return false;
+    }
+}
+
+void
+UnaryOpNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+  if (!this->createAuxVarForUnaryOpNodeInDiffOp())
+    {
+      arg->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+      return;
+    }
+
+  expr_t sthis = this->toStatic(static_datatree);
+  int arg_max_lag = -arg->maxLag();
+  diff_table_t::iterator it = nodes.find(sthis);
+  if (it != nodes.end())
+    {
+      for (map<int, expr_t>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        if (arg == it1->second)
+          return;
+      it->second[arg_max_lag] = const_cast<UnaryOpNode *>(this);
+    }
+  else
+    nodes[sthis][arg_max_lag] = const_cast<UnaryOpNode *>(this);
+}
+
 void
 UnaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -3197,6 +3273,45 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
   return const_cast<VariableNode *>(subst_table[this]);
 }
 
+expr_t
+UnaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  subst_table_t::const_iterator sit = subst_table.find(this);
+  if (sit != subst_table.end())
+    return const_cast<VariableNode *>(sit->second);
+
+  UnaryOpNode *sthis = dynamic_cast<UnaryOpNode *>(this->toStatic(static_datatree));
+  diff_table_t::iterator it = nodes.find(sthis);
+  if (it == nodes.end())
+    {
+      expr_t argsubst = arg->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+      return buildSimilarUnaryOpNode(argsubst, datatree);
+    }
+
+  VariableNode *aux_var = NULL;
+  for (map<int, expr_t>::reverse_iterator rit = it->second.rbegin();
+       rit != it->second.rend(); rit++)
+    {
+      if (rit == it->second.rbegin())
+        {
+          VariableNode *vn = dynamic_cast<VariableNode *>(const_cast<UnaryOpNode *>(this)->get_arg());
+          int symb_id = datatree.symbol_table.addUnaryOpAuxiliaryVar(this->idx, const_cast<UnaryOpNode *>(this),
+                                                                               vn->get_symb_id(), vn->get_lag());
+          aux_var = datatree.AddVariable(symb_id, 0);
+          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(aux_var, const_cast<UnaryOpNode *>(this))));
+          subst_table[rit->second] = dynamic_cast<VariableNode *>(aux_var);
+        }
+      else
+        {
+          VariableNode *vn = dynamic_cast<VariableNode *>(dynamic_cast<UnaryOpNode *>(rit->second)->get_arg());
+          subst_table[rit->second] = dynamic_cast<VariableNode *>(aux_var->decreaseLeadsLags(-vn->get_lag()));
+        }
+    }
+
+  sit = subst_table.find(this);
+  return const_cast<VariableNode *>(sit->second);
+}
+
 expr_t
 UnaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -4865,6 +4980,13 @@ BinaryOpNode::substituteAdl() const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+void
+BinaryOpNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+  arg1->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+  arg2->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+}
+
 void
 BinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
@@ -4881,6 +5003,14 @@ BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+expr_t
+BinaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 bool
 BinaryOpNode::isDiffPresent() const
 {
@@ -5766,6 +5896,14 @@ TrinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table
   arg3->findDiffNodes(static_datatree, diff_table);
 }
 
+void
+TrinaryOpNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+  arg1->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+  arg2->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+  arg3->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+}
+
 expr_t
 TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                               vector<BinaryOpNode *> &neweqs) const
@@ -5776,6 +5914,15 @@ TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_tabl
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+expr_t
+TrinaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 bool
 TrinaryOpNode::isDiffPresent() const
 {
@@ -6189,6 +6336,13 @@ AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_tabl
     (*it)->findDiffNodes(static_datatree, diff_table);
 }
 
+void
+AbstractExternalFunctionNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->findUnaryOpNodesForAuxVarCreation(static_datatree, nodes);
+}
+
 expr_t
 AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                              vector<BinaryOpNode *> &neweqs) const
@@ -6199,6 +6353,15 @@ AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_tab
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+AbstractExternalFunctionNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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)->substituteUnaryOpNodes(static_datatree, nodes, subst_table, neweqs));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 bool
 AbstractExternalFunctionNode::isDiffPresent() const
 {
@@ -7752,6 +7915,11 @@ VarExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_
 {
 }
 
+void
+VarExpectationNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+}
+
 expr_t
 VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                    vector<BinaryOpNode *> &neweqs) const
@@ -7759,6 +7927,12 @@ VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff
   return const_cast<VarExpectationNode *>(this);
 }
 
+expr_t
+VarExpectationNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
 expr_t
 VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
@@ -8186,6 +8360,11 @@ PacExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_
 {
 }
 
+void
+PacExpectationNode::findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const
+{
+}
+
 expr_t
 PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table,
                                    vector<BinaryOpNode *> &neweqs) const
@@ -8193,6 +8372,12 @@ PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff
   return const_cast<PacExpectationNode *>(this);
 }
 
+expr_t
+PacExpectationNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
 expr_t
 PacExpectationNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index f1c1772f7e0b5092e861f32280d2def8695ca85e..221b85aa055f5d9e44b1e5419af0d78b7d9f4f42 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -492,7 +492,9 @@ class ExprNode
 
       //! Substitute diff operator
       virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const = 0;
+      virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) const = 0;
       virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+      virtual expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
       virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
@@ -584,7 +586,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -672,7 +676,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -782,8 +788,11 @@ 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;
+  bool createAuxVarForUnaryOpNodeInDiffOp() const;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -910,7 +919,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1013,7 +1024,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, 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;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1116,7 +1129,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &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;
@@ -1309,7 +1324,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &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,
@@ -1392,7 +1409,9 @@ 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;
+  virtual void findUnaryOpNodesForAuxVarCreation(DataTree &static_datatree, diff_table_t &nodes) 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 substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &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 0ad90e612ee85a95426ea9c64c9690584d842bfd..540a59a4e671be2cb532eb29fb9ec896f75fe291 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -347,7 +347,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
 }
 
 void
-ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput)
+ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput, const bool transform_unary_ops)
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   // - except adl and diff which we always want expanded
@@ -366,6 +366,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
         }
     }
 
+  if (transform_unary_ops)
+    dynamic_model.substituteUnaryOps(diff_static_model);
+
   // Create auxiliary variable and equations for Diff operator
   ExprNode::subst_table_t diff_subst_table;
   dynamic_model.substituteDiff(diff_static_model, diff_subst_table);
diff --git a/src/ModFile.hh b/src/ModFile.hh
index 08ee6c4213b839a854b705ab1b0acf70e43ae3e1..6664c88b2d1185fef3868a756185879426e20b2b 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -136,7 +136,7 @@ public:
   void checkPass(bool nostrict, bool stochastic);
   //! Perform some transformations on the model (creation of auxiliary vars and equations)
   /*! \param compute_xrefs if true, equation cross references will be computed */
-  void transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput);
+  void transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput, const bool transform_unary_ops);
   //! Execute computations
   /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
   /*! \param params_derivs_order compute this order of derivs wrt parameters */
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 2ec4060a3096c38e5474b0e77f838d468452afc1..79ebe86c7bbe7607db048877dc141238ead7dbab 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -373,6 +373,7 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
           case avEndoLag:
           case avExoLag:
           case avVarModel:
+          case avUnaryOp:
             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;
             break;
@@ -506,6 +507,7 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
             case avEndoLag:
             case avExoLag:
             case avVarModel:
+            case avUnaryOp:
               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;
               break;
@@ -605,6 +607,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
         case avEndoLag:
         case avExoLag:
         case avVarModel:
+        case avUnaryOp:
           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;
           break;
@@ -805,6 +808,28 @@ SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenExcept
   return addDiffAuxiliaryVar(index, expr_arg, -1, 0);
 }
 
+int
+SymbolTable::addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_UOP_IN_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, avUnaryOp, orig_symb_id, orig_lag, 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)
 {
@@ -1133,6 +1158,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
             case avEndoLag:
             case avExoLag:
             case avVarModel:
+            case avUnaryOp:
               output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
                      << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
               break;
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 58dee978321161bfe82f8ec728908f29fd250621..1649f96349cc70f3a2eb1858db3ce26ae5f5b0cf 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -36,16 +36,17 @@ typedef class ExprNode *expr_t;
 //! Types of auxiliary variables
 enum aux_var_t
   {
-    avEndoLead = 0,       //!< Substitute for endo leads >= 2
-    avEndoLag = 1,        //!< Substitute for endo lags >= 2
-    avExoLead = 2,        //!< Substitute for exo leads >= 1
-    avExoLag = 3,         //!< Substitute for exo lags >= 1
-    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
-    avDiff = 8,           //!< Variable for Diff operator
-    avDiffLag = 9         //!< Variable for timing between Diff operators
+    avEndoLead = 0,    //!< Substitute for endo leads >= 2
+    avEndoLag = 1,     //!< Substitute for endo lags >= 2
+    avExoLead = 2,     //!< Substitute for exo leads >= 1
+    avExoLag = 3,      //!< Substitute for exo lags >= 1
+    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
+    avDiff = 8,        //!< Variable for Diff operator
+    avDiffLag = 9,     //!< Variable for timing between Diff operators
+    avUnaryOp = 10     //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
   };
 
 //! Information on some auxiliary variables
@@ -300,6 +301,8 @@ public:
   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);
+  //! An Auxiliary variable for a unary op
+  int addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const