v4 parser:

* added third order derivatives in dynamic file (triggered by order=3 option of stoch_simul/olr/osr/ramsey_policy)
* don't output hessian in dynamic file if order=1 is specified


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1180 ac1d8469-bf42-47a9-8791-bf33cf982152

parser.src/ComputingTasks.cc

@@ -71,6 +71,11 @@ void
 StochSimulStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = atoi(it->second.c_str());
 }
 
 void
@@ -92,6 +97,11 @@ void
 RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = atoi(it->second.c_str());
 }
 
 void
@@ -113,6 +123,11 @@ void
 EstimationStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = atoi(it->second.c_str());
 }
 
 void
@@ -543,6 +558,11 @@ void
 OsrStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = atoi(it->second.c_str());
 }
 
 void
@@ -576,6 +596,11 @@ void
 OlrStatement::checkPass(ModFileStructure &mod_file_struct)
 {
   mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = atoi(it->second.c_str());
 }
 
 void

parser.src/ModFile.cc

@@ -47,8 +47,16 @@ ModFile::computingPass()
     model_tree.computeJacobian = true;
   else
     {
+      if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
+        {
+          cerr << "Incorrect order option..." << endl;
+          exit(-1);
+        }
       model_tree.computeJacobianExo = true;
+      if (mod_file_struct.order_option >= 2)
         model_tree.computeHessian = true;
+      if (mod_file_struct.order_option == 3)
+        model_tree.computeThirdDerivatives = true;
     }
 
   model_tree.computingPass();

parser.src/ModelTree.cc

@@ -11,7 +11,8 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
   computeJacobian(false),
   computeJacobianExo(false),
   computeHessian(false),
-  computeStaticHessian(false)
+  computeStaticHessian(false),
+  computeThirdDerivatives(false)
 {
 }
 
@@ -31,7 +32,7 @@ ModelTree::derive(int order)
       }
   cout << "done" << endl;
 
-  if (order == 2)
+  if (order >= 2)
     {
       cout << "  Processing Order 2... ";
       for(first_derivatives_type::const_iterator it = first_derivatives.begin();
@@ -52,6 +53,32 @@ ModelTree::derive(int order)
         }
       cout << "done" << endl;
     }
+
+  if (order >= 3)
+    {
+      cout << "  Processing Order 3... ";
+      for(second_derivatives_type::const_iterator it = second_derivatives.begin();
+          it != second_derivatives.end(); it++)
+        {
+          int eq = it->first.first;
+
+          int var1 = it->first.second.first;
+          int var2 = it->first.second.second;
+          // By construction, var2 <= var1
+
+          NodeID d2 = it->second;
+
+          // Store only third derivatives such that var3 <= var2 <= var1
+          for(int var3 = 0; var3 <= var2; var3++)
+            {
+              NodeID d3 = d2->getDerivative(var3);
+              if (d3 == Zero)
+                continue;
+              third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
+            }
+        }
+      cout << "done" << endl;
+    }
 }
 
 void
@@ -68,10 +95,15 @@ ModelTree::computeTemporaryTerms(int order)
       it != first_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count, temporary_terms);
 
-  if (order == 2)
+  if (order >= 2)
     for(second_derivatives_type::iterator it = second_derivatives.begin();
         it != second_derivatives.end(); it++)
       it->second->computeTemporaryTerms(reference_count, temporary_terms);
+
+  if (order >= 3)
+    for(third_derivatives_type::iterator it = third_derivatives.begin();
+        it != third_derivatives.end(); it++)
+      it->second->computeTemporaryTerms(reference_count, temporary_terms);
 }
 
 void
@@ -170,7 +202,7 @@ ModelTree::writeDynamicMFile(const string &dynamic_basename) const
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(-1);
     }
-  mDynamicModelFile << "function [residual, g1, g2] = " << dynamic_basename << "(y, x)\n";
+  mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x)\n";
   mDynamicModelFile << interfaces::comment()+"\n"+interfaces::comment();
   mDynamicModelFile << "Status : Computes dynamic model for Dynare\n" << interfaces::comment() << "\n";
   mDynamicModelFile << interfaces::comment();
@@ -476,6 +508,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
   ostringstream model_output;    // Used for storing model equations
   ostringstream jacobian_output; // Used for storing jacobian equations
   ostringstream hessian_output;  // Used for storing Hessian equations
+  ostringstream third_derivatives_output;
 
   writeTemporaryTerms(model_output, true);
 
@@ -483,6 +516,14 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
 
   writeModelEquations(model_output, true);
 
+  int nrows = equations.size();
+  int nvars;
+  if (computeJacobianExo)
+    nvars = variable_table.get_dyn_var_nbr();
+  else
+    nvars = variable_table.var_endo_nbr;
+  int nvars_sq = nvars * nvars;
+
   // Writing Jacobian
   if (computeJacobian || computeJacobianExo)
     for(first_derivatives_type::const_iterator it = first_derivatives.begin();
@@ -516,8 +557,8 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
         int id1 = variable_table.getSortID(var1);
         int id2 = variable_table.getSortID(var2);
 
-        int col_nb = id1*variable_table.get_dyn_var_nbr()+id2+1;
-        int col_nb_sym = id2*variable_table.get_dyn_var_nbr()+id1+1;
+        int col_nb = id1*nvars+id2+1;
+        int col_nb_sym = id2*nvars+id1+1;
 
         hessian_output << "  g2" << lpar << eq+1 << ", " << col_nb << rpar << " = ";
         d2->writeOutput(hessian_output, true, temporary_terms);
@@ -529,10 +570,42 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
                     <<  "g2" << lpar << eq+1 << ", " << col_nb << rpar << ";" << endl;
       }
 
-  int nrows = equations.size();
-  int nvars = variable_table.var_endo_nbr;
-  if (computeJacobianExo)
-    nvars += symbol_table.exo_nbr + symbol_table.exo_det_nbr;
+  // Writing third derivatives
+  if (computeThirdDerivatives)
+    for(third_derivatives_type::const_iterator it = third_derivatives.begin();
+        it != third_derivatives.end(); it++)
+      {
+        int eq = it->first.first;
+        int var1 = it->first.second.first;
+        int var2 = it->first.second.second.first;
+        int var3 = it->first.second.second.second;
+        NodeID d3 = it->second;
+
+        int id1 = variable_table.getSortID(var1);
+        int id2 = variable_table.getSortID(var2);
+        int id3 = variable_table.getSortID(var3);
+
+        // Reference column number for the g3 matrix
+        int ref_col = id1 * nvars_sq + id2 * nvars + id3 + 1;
+
+        third_derivatives_output << "  g3" << lpar << eq+1 << ", " << ref_col << rpar << " = ";
+        d3->writeOutput(third_derivatives_output, true, temporary_terms);
+        third_derivatives_output << ";" << endl;
+
+        // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
+        set<int> cols;
+        cols.insert(id1 * nvars_sq + id3 * nvars + id2 + 1);
+        cols.insert(id2 * nvars_sq + id1 * nvars + id3 + 1);
+        cols.insert(id2 * nvars_sq + id3 * nvars + id1 + 1);
+        cols.insert(id3 * nvars_sq + id1 * nvars + id2 + 1);
+        cols.insert(id3 * nvars_sq + id2 * nvars + id1 + 1);
+
+        for(set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
+          if (*it2 != ref_col)
+            third_derivatives_output << "  g3" << lpar << eq+1 << ", " << *it2 << rpar << " = "
+                                     << "g3" << lpar << eq+1 << ", " << ref_col << rpar
+                                     << ";" << endl;
+      }
 
   if (offset == 1)
     {
@@ -561,16 +634,25 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
         {
           DynamicOutput << "if nargout >= 3,\n";
           // Writing initialization instruction for matrix g2
-          int ncols = nvars*nvars;
-          DynamicOutput << "  g2 = " <<
-            "sparse([],[],[]," << nrows << ", " << ncols << ", " <<
-            5*ncols << ");\n";
-          DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
-          DynamicOutput << "Hessian matrix\n\t";
-          DynamicOutput << interfaces::comment() + "\n\n";
+          int ncols = nvars_sq;
+          DynamicOutput << "  g2 = sparse([],[],[]," << nrows << ", " << ncols << ", "
+                        << 5*ncols << ");\n";
+          DynamicOutput << "\n\t"+interfaces::comment() << "\n\t" << interfaces::comment();
+          DynamicOutput << "Hessian matrix\n\t" << interfaces::comment() << "\n\n";
           DynamicOutput << hessian_output.str() << lsymetric.str();
           DynamicOutput << "end;\n";
         }
+      if (computeThirdDerivatives)
+        {
+          DynamicOutput << "if nargout >= 4,\n";
+          int ncols = nvars_sq * nvars;
+          DynamicOutput << "  g3 = sparse([],[],[]," << nrows << ", " << ncols << ", "
+                        << 5*ncols << ");\n";
+          DynamicOutput << "\n\t" + interfaces::comment() + "\n\t" + interfaces::comment();
+          DynamicOutput << "Third order derivatives\n\t" << interfaces::comment() << "\n\n";
+          DynamicOutput << third_derivatives_output.str();
+          DynamicOutput << "end;\n";
+        }
     }
   else
     {
@@ -709,16 +791,16 @@ ModelTree::computingPass()
       rpar = ']';
     }
 
-  if (computeHessian || computeStaticHessian)
-    {
-      derive(2);
-      computeTemporaryTerms(2);
-    }
-  else
-    {
-      derive(1);
-      computeTemporaryTerms(1);
-    }
+  // Determine derivation order
+  int order = 1;
+  if (computeThirdDerivatives)
+    order = 3;
+  else if (computeHessian || computeStaticHessian)
+    order = 2;
+
+  // Launch computations
+  derive(order);
+  computeTemporaryTerms(order);
 }
 
 void

parser.src/Statement.cc

@@ -3,7 +3,8 @@
 ModFileStructure::ModFileStructure() :
   check_present(false),
   simul_present(false),
-  stoch_simul_or_similar_present(false)
+  stoch_simul_or_similar_present(false),
+  order_option(2)
 {
 }
 

parser.src/include/ModFile.hh

@@ -43,7 +43,6 @@ public:
   /*!
     \param basename The base name used for writing output files. Should be the name of the mod file without its extension
     \param clear_all Should a "clear all" instruction be written to output ?
-    \todo make this method const
   */
   void writeOutputFiles(const string &basename, bool clear_all) const;
 };

parser.src/include/ModelTree.hh

@@ -36,6 +36,15 @@ private:
   */
   second_derivatives_type second_derivatives;
 
+  typedef map<pair<int, pair<int, pair<int, int> > >, NodeID> third_derivatives_type;
+  //! Third order derivatives
+  /*! First index is equation number, second, third and fourth are variables w.r. to which is computed the derivative.
+    Only non-null derivatives are stored in the map.
+    Contains only third order derivatives where var1 >= var2 >= var3 (for obvious symmetry reasons).
+    Variable indexes used are those of the variable_table, before sorting.
+  */
+  third_derivatives_type third_derivatives;
+
   //! Temporary terms (those which will be noted Txxxx)
   temporary_terms_type temporary_terms;
 
@@ -54,6 +63,7 @@ private:
   /*! \todo handle hessian in C output */
   void writeStaticModel(ostream &StaticOutput) const;
   //! Writes the dynamic model equations and its derivatives
+  /*! \todo add third derivatives handling in C output */
   void writeDynamicModel(ostream &DynamicOutput) const;
   //! Writes static model file (Matlab version)
   void writeStaticMFile(const string &static_basename) const;
@@ -62,6 +72,7 @@ private:
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &dynamic_basename) const;
   //! Writes dynamic model file (C version)
+  /*! \todo add third derivatives handling */
   void writeDynamicCFile(const string &dynamic_basename) const;
 
 public:
@@ -78,9 +89,10 @@ public:
   bool computeHessian;
   //! Whether static Hessian (w.r. to endogenous only) should be written
   bool computeStaticHessian;
+  //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written
+  bool computeThirdDerivatives;
   //! Execute computations (variable sorting + derivation)
-  /*! You must set computeJacobian, computeJacobianExo and computeHessian to correct values before
-    calling this function */
+  /*! You must set computeJacobian, computeJacobianExo, computeHessian, computeStaticHessian and computeThirdDerivatives to correct values before calling this function */
   void computingPass();
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output) const;

parser.src/include/Statement.hh

@@ -15,8 +15,11 @@ public:
   bool check_present;
   //! Whether a simul statement is present
   bool simul_present;
-  //! Whether a stoch_simul, estimation, olr, osr statement is present
+  //! Whether a stoch_simul, estimation, olr, osr, ramsey_policy statement is present
   bool stoch_simul_or_similar_present;
+  //! The value of the "order" option of stoch_simul, estimation, olr, osr, ramsey_policy
+  /*! Defaults to 2 */
+  int order_option;
 };
 
 class Statement