diff --git a/DynamicModel.cc b/DynamicModel.cc
index 10ca5aa2fcb854ce5600ebcfac5affe04cb40c1f..2dfba144b54e0d72468f87df2eae5a282f1a1a5f 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3924,29 +3924,29 @@ DynamicModel::writeLatexFile(const string &basename) const
 void
 DynamicModel::substituteEndoLeadGreaterThanTwo(bool deterministic_model)
 {
-  substituteLeadLagInternal(avEndoLead, deterministic_model);
+  substituteLeadLagInternal(avEndoLead, deterministic_model, vector<string>());
 }
 
 void
 DynamicModel::substituteEndoLagGreaterThanTwo(bool deterministic_model)
 {
-  substituteLeadLagInternal(avEndoLag, deterministic_model);
+  substituteLeadLagInternal(avEndoLag, deterministic_model, vector<string>());
 }
 
 void
 DynamicModel::substituteExoLead(bool deterministic_model)
 {
-  substituteLeadLagInternal(avExoLead, deterministic_model);
+  substituteLeadLagInternal(avExoLead, deterministic_model, vector<string>());
 }
 
 void
 DynamicModel::substituteExoLag(bool deterministic_model)
 {
-  substituteLeadLagInternal(avExoLag, deterministic_model);
+  substituteLeadLagInternal(avExoLag, deterministic_model, vector<string>());
 }
 
 void
-DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model)
+DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model, const vector<string> &subset)
 {
   ExprNode::subst_table_t subst_table;
   vector<BinaryOpNode *> neweqs;
@@ -3976,7 +3976,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           subst = value->substituteExoLag(subst_table, neweqs);
           break;
         case avDiffForward:
-          subst = value->differentiateForwardVars(subst_table, neweqs);
+          subst = value->differentiateForwardVars(subset, subst_table, neweqs);
           break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
@@ -4004,7 +4004,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           subst = equations[i]->substituteExoLag(subst_table, neweqs);
           break;
         case avDiffForward:
-          subst = equations[i]->differentiateForwardVars(subst_table, neweqs);
+          subst = equations[i]->differentiateForwardVars(subset, subst_table, neweqs);
           break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
@@ -4134,9 +4134,9 @@ DynamicModel::removeTrendVariableFromEquations()
 }
 
 void
-DynamicModel::differentiateForwardVars()
+DynamicModel::differentiateForwardVars(const vector<string> &subset)
 {
-  substituteLeadLagInternal(avDiffForward, true);
+  substituteLeadLagInternal(avDiffForward, true, subset);
 }
 
 void
diff --git a/DynamicModel.hh b/DynamicModel.hh
index b3ca75c9bd45de71bfba4ef6951d4ecada9d9262..63ac665f4807d79d7079d034045fb24b860c2fcb 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -139,8 +139,11 @@ private:
   void collectBlockVariables();
 
   //! Factorized code for substitutions of leads/lags
-  /*! \param[in] type determines which type of variables is concerned */
-  void substituteLeadLagInternal(aux_var_t type, bool deterministic_model);
+  /*! \param[in] type determines which type of variables is concerned
+    \param[in] deterministic_model whether we are in a deterministic model (only for exogenous leads/lags)
+    \param[in] subset variables to which to apply the transformation (only for diff of forward vars)
+  */
+  void substituteLeadLagInternal(aux_var_t type, bool deterministic_model, const vector<string> &subset);
 
 private:
   //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
@@ -281,7 +284,9 @@ public:
   void removeTrendVariableFromEquations();
 
   //! Transforms the model by creating aux vars for the diff of forward vars
-  void differentiateForwardVars();
+  /*! If subset is empty, does the transformation for all fwrd vars; otherwise
+      restrict it to the vars in subset */
+  void differentiateForwardVars(const vector<string> &subset);
 
   //! Fills eval context with values of model local variables and auxiliary variables
   void fillEvalContext(eval_context_t &eval_context) const;
diff --git a/DynareBison.yy b/DynareBison.yy
index 92d4927555a3b5023729fddb32a6d3b4c3425ae6..7bc5698f193aeb083fba553259eccbc538b9f255 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -568,7 +568,8 @@ model_options : BLOCK { driver.block(); }
               | BYTECODE { driver.byte_code(); }
               | USE_DLL { driver.use_dll(); }
               | NO_STATIC { driver.no_static();}
-              | DIFFERENTIATE_FORWARD_VARS { driver.differentiate_forward_vars(); }
+              | DIFFERENTIATE_FORWARD_VARS { driver.differentiate_forward_vars_all(); }
+              | DIFFERENTIATE_FORWARD_VARS EQUAL '(' symbol_list ')' { driver.differentiate_forward_vars_some(); }
               | o_linear
               ;
 
diff --git a/ExprNode.cc b/ExprNode.cc
index 59d26fa264274ffa9679d0fc79cb33eee81a11c6..930ee3e0d85aea5fad74bde54ca26a4e4755bc8c 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -431,7 +431,7 @@ NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-NumConstNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1225,14 +1225,16 @@ VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-VariableNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   expr_t value;
   switch (type)
     {
     case eEndogenous:
       assert(lag <= 1);
-      if (lag <= 0)
+      if (lag <= 0
+          || (subset.size() > 0
+              && find(subset.begin(), subset.end(), datatree.symbol_table.getName(symb_id)) == subset.end()))
         return const_cast<VariableNode *>(this);
       else
         {
@@ -1255,7 +1257,7 @@ VariableNode::differentiateForwardVars(subst_table_t &subst_table, vector<Binary
       if (value->maxEndoLead() <= 0)
         return const_cast<VariableNode *>(this);
       else
-        return value->differentiateForwardVars(subst_table, neweqs);
+        return value->differentiateForwardVars(subset, subst_table, neweqs);
     default:
       return const_cast<VariableNode *>(this);
     }
@@ -2339,9 +2341,9 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
 }
 
 expr_t
-UnaryOpNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t argsubst = arg->differentiateForwardVars(subst_table, neweqs);
+  expr_t argsubst = arg->differentiateForwardVars(subset, subst_table, neweqs);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
@@ -3583,10 +3585,10 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 
 
 expr_t
-BinaryOpNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->differentiateForwardVars(subst_table, neweqs);
-  expr_t arg2subst = arg2->differentiateForwardVars(subst_table, neweqs);
+  expr_t arg1subst = arg1->differentiateForwardVars(subset, subst_table, neweqs);
+  expr_t arg2subst = arg2->differentiateForwardVars(subset, subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -4200,11 +4202,11 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
 
 
 expr_t
-TrinaryOpNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  expr_t arg1subst = arg1->differentiateForwardVars(subst_table, neweqs);
-  expr_t arg2subst = arg2->differentiateForwardVars(subst_table, neweqs);
-  expr_t arg3subst = arg3->differentiateForwardVars(subst_table, neweqs);
+  expr_t arg1subst = arg1->differentiateForwardVars(subset, subst_table, neweqs);
+  expr_t arg2subst = arg2->differentiateForwardVars(subset, subst_table, neweqs);
+  expr_t arg3subst = arg3->differentiateForwardVars(subset, subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -4765,11 +4767,11 @@ ExternalFunctionNode::substituteExpectation(subst_table_t &subst_table, vector<B
 }
 
 expr_t
-ExternalFunctionNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+ExternalFunctionNode::differentiateForwardVars(const vector<string> &subset, 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)->differentiateForwardVars(subst_table, neweqs));
+    arguments_subst.push_back((*it)->differentiateForwardVars(subset, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
diff --git a/ExprNode.hh b/ExprNode.hh
index 7ba4ccc2cf36d8e8dec3705bdb66f6521fea827f..917643f1f6e831cb89fbf8c01f859a3165b06b54 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -361,11 +361,13 @@ public:
 
   //! Constructs a new expression where forward variables (supposed to be at most in t+1) have been replaced by themselves at t, plus a new aux var representing their (time) differentiate
   /*!
+    \param[in] subset variables to which to limit the transformation; transform
+    all fwrd vars if empty
     \param[in,out] subst_table Map used to store mapping between a given
     forward variable and the aux var that contains its differentiate
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
   //! Return true if the nodeID is a numerical constant equal to value and false otherwise
   /*!
@@ -452,7 +454,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 decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -514,7 +516,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 decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -591,7 +593,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 decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -681,7 +683,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 decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -751,7 +753,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 decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -825,7 +827,7 @@ public:
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
-  virtual expr_t differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
diff --git a/ModFile.cc b/ModFile.cc
index a966c353f99bf4d8ea8da6268c8ae5df6e884eaf..86296273fccf48d5c178e049a4f213f193cc55b0 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -322,7 +322,7 @@ ModFile::transformPass()
     }
 
   if (differentiate_forward_vars)
-    dynamic_model.differentiateForwardVars();
+    dynamic_model.differentiateForwardVars(differentiate_forward_vars_subset);
 
   if (mod_file_struct.dsge_var_estimated || !mod_file_struct.dsge_var_calibrated.empty())
     try
diff --git a/ModFile.hh b/ModFile.hh
index ceba97b1ccbd0b03d414e111957ccc16cad7b197..d2784d1a436e9cd4280027aef4987afca8a21379 100644
--- a/ModFile.hh
+++ b/ModFile.hh
@@ -78,6 +78,12 @@ public:
   //! Is the 'differentiate_forward_vars' option used?
   bool differentiate_forward_vars;
 
+  /*! If the 'differentiate_forward_vars' option is used, contains the set of
+      endogenous with respect to which to do the transformation;
+      if empty, means that the transformation must be applied to all endos
+      with a lead */
+  vector<string> differentiate_forward_vars_subset;
+
   //! Are nonstationary variables present ?
   bool nonstationary_variables;
 
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 2d8d4fa9c93e8faa6515990e457780a5a6d8919f..9747b6408c7fea571314d43f34c1e4d774f42602 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -494,11 +494,26 @@ ParsingDriver::byte_code()
 }
 
 void
-ParsingDriver::differentiate_forward_vars()
+ParsingDriver::differentiate_forward_vars_all()
 {
   mod_file->differentiate_forward_vars = true;
 }
 
+void
+ParsingDriver::differentiate_forward_vars_some()
+{
+  mod_file->differentiate_forward_vars = true;
+  mod_file->differentiate_forward_vars_subset = symbol_list.get_symbols();
+  for (vector<string>::const_iterator it = mod_file->differentiate_forward_vars_subset.begin();
+       it != mod_file->differentiate_forward_vars_subset.end(); ++it)
+    {
+      check_symbol_existence(*it);
+      if (mod_file->symbol_table.getType(*it) != eEndogenous)
+        error("Symbol " + *it + " is not an endogenous");
+    }
+  symbol_list.clear();
+}
+
 void
 ParsingDriver::cutoff(string *value)
 {
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index 7da190a25875958ff27c4f29df2af3406ca982a1..6d3fda9f09309234fbb20ac3b2df0eaf23dcca89 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -245,8 +245,10 @@ public:
   void byte_code();
   //! the static model is not computed
   void no_static();
-  //! the differentiate_forward_vars option is enabled
-  void differentiate_forward_vars();
+  //! the differentiate_forward_vars option is enabled (for all vars)
+  void differentiate_forward_vars_all();
+  //! the differentiate_forward_vars option is enabled (for a subset of vars)
+  void differentiate_forward_vars_some();
   //! cutoff option of model block
   void cutoff(string *value);
   //! mfs option of model block