trunk preprocessor:

* added support for derivatives of Hessian w.r. to parameters
* added "identification" command which does nothing for the moment, except triggering those derivatives w.r. to params


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

ComputingTasks.cc

@@ -953,3 +953,18 @@ BVARForecastStatement::writeOutput(ostream &output, const string &basename) cons
   options_list.writeOutput(output);
   output << "bvar_forecast(" << nlags << ");" << endl;
 }
+
+IdentificationStatement::IdentificationStatement()
+{
+}
+
+void
+IdentificationStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.identification_present = true;
+}
+
+void
+IdentificationStatement::writeOutput(ostream &output, const string &basename) const
+{
+}

ComputingTasks.hh

@@ -429,4 +429,13 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
+class IdentificationStatement : public Statement
+{
+public:
+  IdentificationStatement();
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  /*! \todo add something inside this method when Matlab code is available */
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
 #endif

DynamicModel.cc

@@ -1619,7 +1619,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
 
   writeModelLocalVariables(model_output, output_type);
 
-  writeTemporaryTerms(model_output, output_type);
+  writeTemporaryTerms(temporary_terms, model_output, output_type);
 
   writeModelEquations(model_output, output_type);
 
@@ -2150,10 +2150,10 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock)
 }
 
 void
-DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives,
+DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives,
                             const eval_context_type &eval_context, bool no_tmp_terms)
 {
-  if (!jacobianExo && (hessian || thirdDerivatives))
+  if (!jacobianExo && (hessian || thirdDerivatives || paramsDerivatives))
     {
       cerr << "DynamicModel::computingPass: computing 2nd or 3rd order derivatives imply computing 1st derivatives w.r. to exogenous" << endl;
       exit(EXIT_FAILURE);
@@ -2182,6 +2182,15 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
       computeHessian(vars);
     }
 
+  if (paramsDerivatives)
+    {
+      cout << " - order 2 (derivatives of Jacobian w.r. to parameters)" << endl;
+      computeParamsDerivatives();
+
+      if (!no_tmp_terms)
+        computeParamsDerivativesTemporaryTerms();
+    }
+
   if (thirdDerivatives)
     {
       cout << " - order 3" << endl;
@@ -2291,6 +2300,9 @@ DynamicModel::computeDerivID(int symb_id, int lag)
       else if (-max_exo_det_lag > lag)
         max_exo_det_lag = -lag;
       break;
+    case eParameter:
+      // We wan't to compute a derivation ID for parameters
+      break;
     default:
       return -1;
     }
@@ -2384,8 +2396,12 @@ DynamicModel::computeDynJacobianCols(bool jacobianExo)
           if (jacobianExo)
             dyn_jacobian_cols_table[deriv_id] = dynJacobianColsNbr + symbol_table.exo_nbr() + tsid;
           break;
+        case eParameter:
+          // We don't assign a dynamic jacobian column to parameters
+          break;
         default:
           // Shut up GCC
+          cerr << "DynamicModel::computeDynJacobianCols: impossible case" << endl;
           exit(EXIT_FAILURE);
         }
     }
@@ -2411,3 +2427,81 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
     return it->second;
 }
 
+void
+DynamicModel::computeParamsDerivatives()
+{
+  for (int param = 0; param < getDerivIDNbr(); param++)
+    {
+      if (getTypeByDerivID(param) != eParameter)
+        continue;
+
+      for (first_derivatives_type::const_iterator it = first_derivatives.begin();
+           it != first_derivatives.end(); it++)
+        {
+          int eq = it->first.first;
+          int var = it->first.second;
+          NodeID d1 = it->second;
+
+          NodeID d2 = d1->getDerivative(param);
+          if (d2 == Zero)
+            continue;
+          params_derivatives[make_pair(eq, make_pair(var, param))] = d2;
+        }
+    }
+}
+
+void
+DynamicModel::computeParamsDerivativesTemporaryTerms()
+{
+  map<NodeID, int> reference_count;
+  params_derivs_temporary_terms.clear();
+
+  for (second_derivatives_type::iterator it = params_derivatives.begin();
+       it != params_derivatives.end(); it++)
+    it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
+}
+
+void
+DynamicModel::writeParamsDerivativesFile(const string &basename) const
+{
+  if (!params_derivatives.size())
+    return;
+
+  string filename = basename + "_params_derivs.m";
+
+  ofstream paramsDerivsFile;
+  paramsDerivsFile.open(filename.c_str(), ios::out | ios::binary);
+  if (!paramsDerivsFile.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  paramsDerivsFile << "function gp = " << basename << "_params_derivs(y, x, params, it_)" << endl
+                   << "%" << endl
+                   << "% Warning : this file is generated automatically by Dynare" << endl
+                   << "%           from model file (.mod)" << endl << endl;
+
+
+  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, oMatlabDynamicModel);
+
+  paramsDerivsFile << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", "
+                   << symbol_table.param_nbr() << ");" << endl;
+
+  for (second_derivatives_type::const_iterator it = params_derivatives.begin();
+       it != params_derivatives.end(); it++)
+    {
+      int eq = it->first.first;
+      int var = it->first.second.first;
+      int param = it->first.second.second;
+      NodeID d2 = it->second;
+
+      int var_col = getDynJacobianCol(var) + 1;
+      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
+
+      paramsDerivsFile << "gp(" << eq+1 << ", " << var_col << ", " << param_col << ") = ";
+      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms);
+      paramsDerivsFile << ";" << endl;
+    }
+
+  paramsDerivsFile.close();
+}

DynamicModel.hh

@@ -56,6 +56,16 @@ private:
   /*! Set by computeDerivID() and computeDynJacobianCols() */
   int dynJacobianColsNbr;
 
+  //! Derivatives of the jacobian w.r. to parameters
+  /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter.
+    Only non-null derivatives are stored in the map.
+    Variable and parameter indices are those of the getDerivID() method.
+  */
+  second_derivatives_type params_derivatives;
+
+  //! Temporary terms for the file containing parameters dervicatives
+  temporary_terms_type params_derivs_temporary_terms;
+
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &dynamic_basename) const;
   //! Writes dynamic model file (C version)
@@ -95,6 +105,10 @@ private:
   int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
   //! Compute the column indices of the dynamic Jacobian
   void computeDynJacobianCols(bool jacobianExo);
+  //! Computes derivatives of the Jacobian w.r. to parameters
+  void computeParamsDerivatives();
+  //! Computes temporary terms for the file containing parameters derivatives
+  void computeParamsDerivativesTemporaryTerms();
 
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
@@ -112,10 +126,11 @@ public:
     \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
     \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
     \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
+    \param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true)
     \param eval_context evaluation context for normalization
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
-  void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives,
+  void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives,
                      const eval_context_type &eval_context, bool no_tmp_terms);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output) const;
@@ -126,6 +141,8 @@ public:
                              const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
   //! Writes dynamic model file
   void writeDynamicFile(const string &basename) const;
+  //! Writes file containing parameters derivatives
+  void writeParamsDerivativesFile(const string &basename) const;
   //! Converts to static model (only the equations)
   /*! It assumes that the static model given in argument has just been allocated */
   void toStatic(StaticModel &static_model) const;

DynareBison.yy

@@ -100,7 +100,7 @@ class ParsingDriver;
 %token FORECAST
 %token GAMMA_PDF GAUSSIAN_ELIMINATION GMRES GRAPH
 %token HISTVAL HP_FILTER HP_NGRID
-%token INF_CONSTANT INITVAL INITVAL_FILE
+%token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE
 %token <string_val> INT_NUMBER
 %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF
 %token KALMAN_ALGO KALMAN_TOL
@@ -136,7 +136,7 @@ class ParsingDriver;
 %nonassoc POWER
 %token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH ASINH ACOSH ATANH SQRT
 /* GSA analysis */
-%token DYNARE_SENSITIVITY IDENTIFICATION MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU GLUE MORRIS_NLIV
+%token DYNARE_SENSITIVITY MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU GLUE MORRIS_NLIV
 %token MORRIS_NTRA NSAM LOAD_REDFORM LOAD_RMSE LOAD_STAB ALPHA2_STAB KSSTAT LOGTRANS_REDFORM THRESHOLD_REDFORM
 %token KSSTAT_REDFORM ALPHA2_REDFORM NAMENDO NAMLAGENDO NAMEXO RMSE LIK_ONLY VAR_RMSE PFILT_RMSE ISTART_RMSE
 %token ALPHA_RMSE ALPHA2_RMSE TRANS_IDENT
@@ -210,6 +210,7 @@ statement : parameters
           | forecast
           | load_params_and_steady_state
           | save_params_and_steady_state
+          | identification
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -1087,27 +1088,29 @@ calib : CALIB ';'
         { driver.run_calib(1); }
       ;
 
-dynatype : DYNATYPE '(' filename ')'';'
+dynatype : DYNATYPE '(' filename ')' ';'
            { driver.run_dynatype($3); }
          | DYNATYPE '(' filename ')' symbol_list ';'
            { driver.run_dynatype($3); }
          ;
 
-dynasave : DYNASAVE '(' filename ')'';'
+dynasave : DYNASAVE '(' filename ')' ';'
            { driver.run_dynasave($3); }
          | DYNASAVE '(' filename ')' symbol_list ';'
            { driver.run_dynasave($3); }
          ;
 
-load_params_and_steady_state: LOAD_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
-          {driver.run_load_params_and_steady_state($3);}
-          ;
-
-save_params_and_steady_state: SAVE_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
-          {driver.run_save_params_and_steady_state($3);}
-          ;
+load_params_and_steady_state : LOAD_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
+                               { driver.run_load_params_and_steady_state($3); }
+                             ;
 
+save_params_and_steady_state : SAVE_PARAMS_AND_STEADY_STATE '(' filename ')' ';'
+                               { driver.run_save_params_and_steady_state($3); }
+                             ;
 
+identification : IDENTIFICATION ';'
+                 { driver.run_identification(); }
+               ;
 
 model_comparison : MODEL_COMPARISON mc_filename_list ';'
                    { driver.run_model_comparison(); }

DynareFlex.ll

@@ -126,6 +126,7 @@ int sigma_e = 0;
 <INITIAL>calib {BEGIN DYNARE_STATEMENT; return token::CALIB;}
 <INITIAL>planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;}
 <INITIAL>ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;}
+<INITIAL>identification {BEGIN DYNARE_STATEMENT; return token::IDENTIFICATION;}
 
 <INITIAL>bvar_density {BEGIN DYNARE_STATEMENT; return token::BVAR_DENSITY; }
 <INITIAL>bvar_forecast {BEGIN DYNARE_STATEMENT; return token::BVAR_FORECAST; }

ExprNode.cc

@@ -198,11 +198,9 @@ VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg,
     case eEndogenous:
     case eExogenous:
     case eExogenousDet:
-      // For a variable, the only non-null derivative is with respect to itself
-      non_null_derivatives.insert(deriv_id);
-      break;
     case eParameter:
-      // All derivatives are null, do nothing
+      // For a variable or a parameter, the only non-null derivative is with respect to itself
+      non_null_derivatives.insert(deriv_id);
       break;
     case eModelLocalVariable:
       // Non null derivatives are those of the value of the local parameter
@@ -225,12 +223,11 @@ VariableNode::computeDerivative(int deriv_id_arg)
     case eEndogenous:
     case eExogenous:
     case eExogenousDet:
+    case eParameter:
       if (deriv_id == deriv_id_arg)
         return datatree.One;
       else
         return datatree.Zero;
-    case eParameter:
-      return datatree.Zero;
     case eModelLocalVariable:
       return datatree.local_variables_table[symb_id]->getDerivative(deriv_id_arg);
     case eModFileLocalVariable:

ModFile.cc

@@ -149,7 +149,7 @@ ModFile::computingPass(bool no_tmp_terms)
       // Set things to compute for dynamic model
       
       if (mod_file_struct.simul_present)
-        dynamic_model.computingPass(false, false, false, global_eval_context, no_tmp_terms);
+        dynamic_model.computingPass(false, false, false, false, global_eval_context, no_tmp_terms);
       else
         {
           if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
@@ -159,7 +159,8 @@ ModFile::computingPass(bool no_tmp_terms)
             }
           bool hessian = mod_file_struct.order_option >= 2;
           bool thirdDerivatives = mod_file_struct.order_option == 3;
-          dynamic_model.computingPass(true, hessian, thirdDerivatives, global_eval_context, no_tmp_terms);
+          bool paramsDerivatives = mod_file_struct.identification_present;
+          dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivatives, global_eval_context, no_tmp_terms);
         }
     }
 
@@ -254,6 +255,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
       dynamic_model.writeOutput(mOutputFile);
       static_model.writeStaticFile(basename);
       dynamic_model.writeDynamicFile(basename);
+      dynamic_model.writeParamsDerivativesFile(basename);
     }
 
   // Print statements

ModelTree.cc

@@ -143,21 +143,22 @@ ModelTree::computeTemporaryTerms()
 }
 
 void
-ModelTree::writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const
+ModelTree::writeTemporaryTerms(const temporary_terms_type &tt, ostream &output,
+                               ExprNodeOutputType output_type) const
 {
-  // A copy of temporary terms
+  // Local var used to keep track of temp nodes already written
   temporary_terms_type tt2;
 
-  if (temporary_terms.size() > 0 && (!OFFSET(output_type)))
-    output << "double\n";
+  if (tt.size() > 0 && (!OFFSET(output_type)))
+    output << "double" << endl;
 
-  for (temporary_terms_type::const_iterator it = temporary_terms.begin();
-       it != temporary_terms.end(); it++)
+  for (temporary_terms_type::const_iterator it = tt.begin();
+       it != tt.end(); it++)
     {
-      if (!OFFSET(output_type) && it != temporary_terms.begin())
+      if (!OFFSET(output_type) && it != tt.begin())
         output << "," << endl;
 
-      (*it)->writeOutput(output, output_type, temporary_terms);
+      (*it)->writeOutput(output, output_type, tt);
       output << " = ";
 
       (*it)->writeOutput(output, output_type, tt2);

ModelTree.hh

@@ -89,7 +89,7 @@ protected:
   //! Computes temporary terms (for all equations and derivatives)
   void computeTemporaryTerms();
   //! Writes temporary terms
-  void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const;
+  void writeTemporaryTerms(const temporary_terms_type &tt, ostream &output, ExprNodeOutputType output_type) const;
   //! Writes model local variables
   /*! No temporary term is used in the output, so that local parameters declarations can be safely put before temporary terms declaration in the output files */
   void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const;

ParsingDriver.cc

@@ -990,6 +990,12 @@ ParsingDriver::run_save_params_and_steady_state(string *filename)
   delete filename;
 }
 
+void
+ParsingDriver::run_identification()
+{
+  mod_file->addStatement(new IdentificationStatement());
+}
+
 void
 ParsingDriver::add_mc_filename(string *filename, string *prior)
 {

ParsingDriver.hh

@@ -329,6 +329,7 @@ public:
   void run_dynatype(string *filename);
   void run_load_params_and_steady_state(string *filename);
   void run_save_params_and_steady_state(string *filename);
+  void run_identification();
   void add_mc_filename(string *filename, string *prior = new string("1"));
   void run_model_comparison();
   //! Begin a planner_objective statement

Statement.cc

@@ -29,7 +29,8 @@ ModFileStructure::ModFileStructure() :
   ramsey_policy_present(false),
   order_option(0),
   bvar_density_present(false),
-  bvar_forecast_present(false)
+  bvar_forecast_present(false),
+  identification_present(false)
 {
 }
 

Statement.hh

@@ -54,6 +54,8 @@ public:
   bool bvar_density_present;
   //! Whether a bvar_forecast statement is present
   bool bvar_forecast_present;
+  //! Whether an identification statement is present
+  bool identification_present;
 };
 
 class Statement

StaticModel.cc

@@ -129,7 +129,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
 
   writeModelLocalVariables(model_output, output_type);
 
-  writeTemporaryTerms(model_output, output_type);
+  writeTemporaryTerms(temporary_terms, model_output, output_type);
 
   writeModelEquations(model_output, output_type);