From 75b5f1d18a4d47a341508de2c4c7a18c20b44aea Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 25 Apr 2013 18:09:31 +0200
Subject: [PATCH] Add new model option 'differentiate_forward_vars' (ref. #351)

---
 preprocessor/DynamicModel.cc  | 15 +++++++
 preprocessor/DynamicModel.hh  |  3 ++
 preprocessor/DynareBison.yy   |  3 +-
 preprocessor/DynareFlex.ll    |  1 +
 preprocessor/ExprNode.cc      | 78 +++++++++++++++++++++++++++++++++++
 preprocessor/ExprNode.hh      | 14 +++++++
 preprocessor/ModFile.cc       |  4 ++
 preprocessor/ModFile.hh       |  3 ++
 preprocessor/ParsingDriver.cc |  6 +++
 preprocessor/ParsingDriver.hh |  2 +
 preprocessor/SymbolTable.cc   | 26 +++++++++++-
 preprocessor/SymbolTable.hh   | 10 ++++-
 12 files changed, 161 insertions(+), 4 deletions(-)

diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 142a2afd2d..10ca5aa2fc 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 57bd30a793..b3ca75c9bd 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 05a766efb9..92d4927555 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 7f5904f270..8a8f1f4ae3 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 abe47a2ad7..59d26fa264 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 526a26b168..7ba4ccc2cf 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 73236f142a..503c9e20a4 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 6e030d2207..ceba97b1cc 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 ca7c515b4d..2d8d4fa9c9 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 279ba13d95..7da190a258 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 9d1184c564..b90b1ce59d 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 3a3f6ebb3e..a8a032f5a9 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.
-- 
GitLab