From 934f5f21a7b8625584a02c5b6019f6ea7645893c Mon Sep 17 00:00:00 2001 From: sebastien Date: Thu, 8 Feb 2007 11:56:46 +0000 Subject: [PATCH] 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 | 25 ++++++ parser.src/ModFile.cc | 10 ++- parser.src/ModelTree.cc | 136 +++++++++++++++++++++++++------- parser.src/Statement.cc | 3 +- parser.src/include/ModFile.hh | 1 - parser.src/include/ModelTree.hh | 16 +++- parser.src/include/Statement.hh | 5 +- 7 files changed, 163 insertions(+), 33 deletions(-) diff --git a/parser.src/ComputingTasks.cc b/parser.src/ComputingTasks.cc index 88020ee24..85be18e65 100644 --- a/parser.src/ComputingTasks.cc +++ b/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 diff --git a/parser.src/ModFile.cc b/parser.src/ModFile.cc index cd62be38b..50eaddffb 100644 --- a/parser.src/ModFile.cc +++ b/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; - model_tree.computeHessian = 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(); diff --git a/parser.src/ModelTree.cc b/parser.src/ModelTree.cc index 4c0fc49c1..d340717ea 100644 --- a/parser.src/ModelTree.cc +++ b/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 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::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 diff --git a/parser.src/Statement.cc b/parser.src/Statement.cc index bdc4bb492..a92282d5b 100644 --- a/parser.src/Statement.cc +++ b/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) { } diff --git a/parser.src/include/ModFile.hh b/parser.src/include/ModFile.hh index edfdd8181..180cd216c 100644 --- a/parser.src/include/ModFile.hh +++ b/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; }; diff --git a/parser.src/include/ModelTree.hh b/parser.src/include/ModelTree.hh index c95157fef..1940850ed 100644 --- a/parser.src/include/ModelTree.hh +++ b/parser.src/include/ModelTree.hh @@ -36,6 +36,15 @@ private: */ second_derivatives_type second_derivatives; + typedef map > >, 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; diff --git a/parser.src/include/Statement.hh b/parser.src/include/Statement.hh index f71e6975c..e1f98fa75 100644 --- a/parser.src/include/Statement.hh +++ b/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