diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 142a2afd2d3579b177510ed69115ff1323d03bee..10ca5aa2fcb854ce5600ebcfac5affe04cb40c1f 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -3975,6 +3975,9 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         case avExoLag:
           subst = value->substituteExoLag(subst_table, neweqs);
           break;
+        case avDiffForward:
+          subst = value->differentiateForwardVars(subst_table, neweqs);
+          break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
           exit(EXIT_FAILURE);
@@ -4000,6 +4003,9 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         case avExoLag:
           subst = equations[i]->substituteExoLag(subst_table, neweqs);
           break;
+        case avDiffForward:
+          subst = equations[i]->differentiateForwardVars(subst_table, neweqs);
+          break;
         default:
           cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
           exit(EXIT_FAILURE);
@@ -4039,6 +4045,9 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         case avExpectation:
           cout << "expectation";
           break;
+        case avDiffForward:
+          cout << "forward vars";
+          break;
         case avMultiplier:
           cerr << "avMultiplier encountered: impossible case" << endl;
           exit(EXIT_FAILURE);
@@ -4124,6 +4133,12 @@ DynamicModel::removeTrendVariableFromEquations()
     }
 }
 
+void
+DynamicModel::differentiateForwardVars()
+{
+  substituteLeadLagInternal(avDiffForward, true);
+}
+
 void
 DynamicModel::fillEvalContext(eval_context_t &eval_context) const
 {
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index 57bd30a79374fed0b3d5e8885b53da4f3d6e432b..b3ca75c9bd45de71bfba4ef6951d4ecada9d9262 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -280,6 +280,9 @@ public:
   //! Transforms the model by replacing trend variables with a 1
   void removeTrendVariableFromEquations();
 
+  //! Transforms the model by creating aux vars for the diff of forward vars
+  void differentiateForwardVars();
+
   //! Fills eval context with values of model local variables and auxiliary variables
   void fillEvalContext(eval_context_t &eval_context) const;
 
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 05a766efb9628672942b9e68976195b997abc826..92d4927555a3b5023729fddb32a6d3b4c3425ae6 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -96,7 +96,7 @@ class ParsingDriver;
 %token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
 %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
-%token DATAFILE FILE DETERMINISTIC DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
+%token DATAFILE FILE DETERMINISTIC DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
 %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
 %token <string_val> FLOAT_NUMBER
@@ -568,6 +568,7 @@ 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(); }
               | o_linear
               ;
 
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 7f5904f27096cc584d55bcbac60d68c934f0b2f7..8a8f1f4ae32743a6f0d0780c6bb173dd3c7ac194 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -566,6 +566,7 @@ string eofbuff;
 <DYNARE_BLOCK>bytecode {return token::BYTECODE;}
 <DYNARE_BLOCK>all_values_required {return token::ALL_VALUES_REQUIRED;}
 <DYNARE_BLOCK>no_static {return token::NO_STATIC;}
+<DYNARE_BLOCK>differentiate_forward_vars {return token::DIFFERENTIATE_FORWARD_VARS;}
 
 <DYNARE_STATEMENT,DYNARE_BLOCK>linear {return token::LINEAR;}
 
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index abe47a2ad7be7011f41dcb415080e9e5c5908c3c..59d26fa264274ffa9679d0fc79cb33eee81a11c6 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -430,6 +430,12 @@ NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 bool
 NumConstNode::isNumConstNodeEqualTo(double value) const
 {
@@ -1218,6 +1224,43 @@ VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t value;
+  switch (type)
+    {
+    case eEndogenous:
+      assert(lag <= 1);
+      if (lag <= 0)
+        return const_cast<VariableNode *>(this);
+      else
+        {
+          subst_table_t::iterator it = subst_table.find(this);
+          VariableNode *diffvar;
+          if (it != subst_table.end())
+            diffvar = const_cast<VariableNode *>(it->second);
+          else
+            {
+              int aux_symb_id = datatree.symbol_table.addDiffForwardAuxiliaryVar(symb_id);
+              neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), datatree.AddMinus(datatree.AddVariable(symb_id, 0),
+                                                                                                                                      datatree.AddVariable(symb_id, -1)))));
+              diffvar = datatree.AddVariable(aux_symb_id, 1);
+              subst_table[this] = diffvar;
+            }
+          return datatree.AddPlus(datatree.AddVariable(symb_id, 0), diffvar);
+        }
+    case eModelLocalVariable:
+      value = datatree.local_variables_table[symb_id];
+      if (value->maxEndoLead() <= 0)
+        return const_cast<VariableNode *>(this);
+      else
+        return value->differentiateForwardVars(subst_table, neweqs);
+    default:
+      return const_cast<VariableNode *>(this);
+    }
+}
+
 bool
 VariableNode::isNumConstNodeEqualTo(double value) const
 {
@@ -2295,6 +2338,13 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
     }
 }
 
+expr_t
+UnaryOpNode::differentiateForwardVars(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t argsubst = arg->differentiateForwardVars(subst_table, neweqs);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 bool
 UnaryOpNode::isNumConstNodeEqualTo(double value) const
 {
@@ -3531,6 +3581,15 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+
+expr_t
+BinaryOpNode::differentiateForwardVars(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);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 expr_t
 BinaryOpNode::addMultipliersToConstraints(int i)
 {
@@ -4139,6 +4198,16 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+
+expr_t
+TrinaryOpNode::differentiateForwardVars(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);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 bool
 TrinaryOpNode::isNumConstNodeEqualTo(double value) const
 {
@@ -4695,6 +4764,15 @@ ExternalFunctionNode::substituteExpectation(subst_table_t &subst_table, vector<B
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+ExternalFunctionNode::differentiateForwardVars(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));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 expr_t
 ExternalFunctionNode::buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const
 {
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index 526a26b168e289a5f9a8dfc12ce758ea7ed4b6cc..7ba4ccc2cf36d8e8dec3705bdb66f6521fea827f 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -359,6 +359,14 @@ public:
 
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const = 0;
 
+  //! 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,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;
+
   //! Return true if the nodeID is a numerical constant equal to value and false otherwise
   /*!
     \param[in] value of the numerical constante
@@ -444,6 +452,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 bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -505,6 +514,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 bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -581,6 +591,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 bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -670,6 +681,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 bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -739,6 +751,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 bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
@@ -812,6 +825,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 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/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 73236f142a1a7eb88534a47f9abadb3e0d2d4ef6..503c9e20a4d3856547290b393093e7b19ffc55dc 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -38,6 +38,7 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
     static_model(symbol_table, num_constants, external_functions_table),
     steady_state_model(symbol_table, num_constants, external_functions_table, static_model),
     linear(false), block(false), byte_code(false), use_dll(false), no_static(false), 
+    differentiate_forward_vars(false),
     nonstationary_variables(false), ramsey_policy_orig_eqn_nbr(0),
     warnings(warnings_arg)
 {
@@ -307,6 +308,9 @@ ModFile::transformPass()
       dynamic_model.substituteEndoLagGreaterThanTwo(true);
     }
 
+  if (differentiate_forward_vars)
+    dynamic_model.differentiateForwardVars();
+
   if (mod_file_struct.dsge_var_estimated || !mod_file_struct.dsge_var_calibrated.empty())
     try
       {
diff --git a/preprocessor/ModFile.hh b/preprocessor/ModFile.hh
index 6e030d22074737f4fb0d7f7a5d394ad0f9ab199d..ceba97b1ccbd0b03d414e111957ccc16cad7b197 100644
--- a/preprocessor/ModFile.hh
+++ b/preprocessor/ModFile.hh
@@ -75,6 +75,9 @@ public:
   //! Is the static model have to computed (no_static=false) or not (no_static=true). Option of 'model'
   bool no_static;
 
+  //! Is the 'differentiate_forward_vars' option used?
+  bool differentiate_forward_vars;
+
   //! Are nonstationary variables present ?
   bool nonstationary_variables;
 
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index ca7c515b4da7536bbfc176e345c7ac8627679700..2d8d4fa9c93e8faa6515990e457780a5a6d8919f 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -493,6 +493,12 @@ ParsingDriver::byte_code()
   mod_file->byte_code = true;
 }
 
+void
+ParsingDriver::differentiate_forward_vars()
+{
+  mod_file->differentiate_forward_vars = true;
+}
+
 void
 ParsingDriver::cutoff(string *value)
 {
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 279ba13d95d8601f8d8b893b3fd7758e19c01347..7da190a25875958ff27c4f29df2af3406ca982a1 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -245,6 +245,8 @@ public:
   void byte_code();
   //! the static model is not computed
   void no_static();
+  //! the differentiate_forward_vars option is enabled
+  void differentiate_forward_vars();
   //! cutoff option of model block
   void cutoff(string *value);
   //! mfs option of model block
diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc
index 9d1184c5649bac4c34f0a3eaad59a78632e76c7a..b90b1ce59de9c83189b6f98f536f282e7fef3037 100644
--- a/preprocessor/SymbolTable.cc
+++ b/preprocessor/SymbolTable.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -236,6 +236,9 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
           case avMultiplier:
             output << "M_.aux_vars(" << i+1 << ").eq_nbr = '" << aux_vars[i].get_equation_number_for_multiplier() << "';" << endl;
             break;
+          case avDiffForward:
+            output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl;
+            break;
           }
       }
 
@@ -386,6 +389,27 @@ SymbolTable::addMultiplierAuxiliaryVar(int index) throw (FrozenException)
   return symb_id;
 }
 
+int
+SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+  varname << "AUX_DIFF_FWRD_" << orig_symb_id+1;
+
+  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, avDiffForward, orig_symb_id, 0, 0));
+  return symb_id;
+}
+
 int
 SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException)
 {
diff --git a/preprocessor/SymbolTable.hh b/preprocessor/SymbolTable.hh
index 3a3f6ebb3eb6238ca79f5aa2eb6b04e8481e62de..a8a032f5a97c0d1847e163ebbbae6105c8fd0d80 100644
--- a/preprocessor/SymbolTable.hh
+++ b/preprocessor/SymbolTable.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -38,7 +38,7 @@ enum aux_var_t
     avExoLead = 2,        //!< Substitute for exo leads >= 2
     avExoLag = 3,         //!< Substitute for exo lags >= 2
     avExpectation = 4,    //!< Substitute for Expectation Operator
-    // Type 5 now unused
+    avDiffForward = 5,    //!< Substitute for the differentiate of a forward variable
     avMultiplier = 6      //!< Multipliers for FOC of Ramsey Problem
   };
 
@@ -221,6 +221,12 @@ public:
     \return the symbol ID of the new symbol
   */
   int addMultiplierAuxiliaryVar(int index) throw (FrozenException);
+  //! Adds an auxiliary variable for the (time) differentiate of a forward var
+  /*!
+    \param[in] orig_symb_id The symb_id of the forward variable
+    \return the symbol ID of the new symbol
+  */
+  int addDiffForwardAuxiliaryVar(int orig_symb_id) throw (FrozenException);
   //! Searches auxiliary variables which are substitutes for a given symbol_id and lead/lag
   /*!
     The search is only performed among auxiliary variables of endo/exo lag.