From 9dfcf897f7c2162782e5b6d94b136ea76f1d1072 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Fri, 22 Nov 2013 21:09:04 +0100 Subject: [PATCH] adding 3rd order derivatives to Static Model for evaluation of Ramsey policy computed at order = 2 --- ComputingTasks.cc | 6 +-- ModFile.cc | 2 +- StaticModel.cc | 98 ++++++++++++++++++++++++++++++++++++++++++++++- StaticModel.hh | 2 +- 4 files changed, 101 insertions(+), 7 deletions(-) diff --git a/ComputingTasks.cc b/ComputingTasks.cc index b5668549..ae764910 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -179,9 +179,9 @@ RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct, WarningConso if (it != options_list.num_options.end()) { int order = atoi(it->second.c_str()); - if (order > 1) + if (order > 2) { - cerr << "ERROR: ramsey_policy: order > 1 is not yet implemented" << endl; + cerr << "ERROR: ramsey_policy: order > 2 is not implemented" << endl; exit(EXIT_FAILURE); } mod_file_struct.order_option = max(mod_file_struct.order_option, order + 1); @@ -968,7 +968,7 @@ PlannerObjectiveStatement::getPlannerObjective() const void PlannerObjectiveStatement::computingPass() { - model_tree->computingPass(eval_context_t(), false, true, false, false, false); + model_tree->computingPass(eval_context_t(), false, true, true, false, false, false); } void diff --git a/ModFile.cc b/ModFile.cc index 94f5a1a5..2c44e3f2 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -470,7 +470,7 @@ ModFile::computingPass(bool no_tmp_terms) const bool paramsDerivatives = mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation; static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian, - paramsDerivatives, block, byte_code); + false, paramsDerivatives, block, byte_code); } // Set things to compute for dynamic model if (mod_file_struct.simul_present || mod_file_struct.check_present diff --git a/StaticModel.cc b/StaticModel.cc index c78324ea..d82efaed 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -1047,7 +1047,7 @@ StaticModel::collect_first_order_derivatives_endogenous() } void -StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool paramsDerivatives, bool block, bool bytecode) +StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatives, bool paramsDerivatives, bool block, bool bytecode) { initializeVariablesAndEquations(); @@ -1070,6 +1070,12 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms computeHessian(vars); } + if (thirdDerivatives) + { + cout << " - order 3" << endl; + computeThirdDerivatives(vars); + } + if (paramsDerivatives) { cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl; @@ -1142,7 +1148,7 @@ StaticModel::writeStaticMFile(const string &func_name) const exit(EXIT_FAILURE); } - output << "function [residual, g1, g2] = " << func_name + "_static(y, x, params)" << endl + output << "function [residual, g1, g2, g3] = " << func_name + "_static(y, x, params)" << endl << "%" << endl << "% Status : Computes static model for Dynare" << endl << "%" << endl @@ -1160,6 +1166,9 @@ StaticModel::writeStaticMFile(const string &func_name) const << "% g2 [M_.endo_nbr by (M_.endo_nbr)^2] double Hessian matrix of the static model equations;" << endl << "% columns: equations in order of declaration" << endl << "% rows: variables in declaration order" << endl + << "% g3 [M_.endo_nbr by (M_.endo_nbr)^3] double Third derivatives matrix of the static model equations;" << endl + << "% columns: variables in declaration order" << endl + << "% rows: equations in order of declaration" << endl << "%" << endl << "%" << endl << "% Warning : this file is generated automatically by Dynare" << endl @@ -1176,6 +1185,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) 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; // Used for storing third order derivatives equations ExprNodeOutputType output_type = (use_dll ? oCStaticModel : oMatlabStaticModel); deriv_node_temp_terms_t tef_terms; @@ -1185,6 +1195,10 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const writeModelEquations(model_output, output_type); + int nrows = equations.size(); + int JacobianColsNbr = symbol_table.endo_nbr(); + int hessianColsNbr = JacobianColsNbr*JacobianColsNbr; + // Write Jacobian w.r. to endogenous only for (first_derivatives_t::const_iterator it = first_derivatives.begin(); it != first_derivatives.end(); it++) @@ -1247,6 +1261,64 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const } } + // Writing third derivatives + k = 0; // Keep the line of a 3rd derivative in v3 + for (third_derivatives_t::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; + expr_t d3 = it->second; + + int id1 = getSymbIDByDerivID(var1); + int id2 = getSymbIDByDerivID(var2); + int id3 = getSymbIDByDerivID(var3); + + + // Reference column number for the g3 matrix + int ref_col = id1 * hessianColsNbr + id2 * JacobianColsNbr + id3; + + sparseHelper(3, third_derivatives_output, k, 0, output_type); + third_derivatives_output << "=" << eq + 1 << ";" << endl; + + sparseHelper(3, third_derivatives_output, k, 1, output_type); + third_derivatives_output << "=" << ref_col + 1 << ";" << endl; + + sparseHelper(3, third_derivatives_output, k, 2, output_type); + third_derivatives_output << "="; + d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_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 * hessianColsNbr + id3 * JacobianColsNbr + id2); + cols.insert(id2 * hessianColsNbr + id1 * JacobianColsNbr + id3); + cols.insert(id2 * hessianColsNbr + id3 * JacobianColsNbr + id1); + cols.insert(id3 * hessianColsNbr + id1 * JacobianColsNbr + id2); + cols.insert(id3 * hessianColsNbr + id2 * JacobianColsNbr + id1); + + int k2 = 1; // Keeps the offset of the permutation relative to k + for (set::iterator it2 = cols.begin(); it2 != cols.end(); it2++) + if (*it2 != ref_col) + { + sparseHelper(3, third_derivatives_output, k+k2, 0, output_type); + third_derivatives_output << "=" << eq + 1 << ";" << endl; + + sparseHelper(3, third_derivatives_output, k+k2, 1, output_type); + third_derivatives_output << "=" << *it2 + 1 << ";" << endl; + + sparseHelper(3, third_derivatives_output, k+k2, 2, output_type); + third_derivatives_output << "="; + sparseHelper(3, third_derivatives_output, k, 2, output_type); + third_derivatives_output << ";" << endl; + + k2++; + } + k += k2; + } + if (!use_dll) { StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl @@ -1280,6 +1352,20 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const else StaticOutput << " g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl; StaticOutput << "end" << endl; + // Initialize g3 matrix + StaticOutput << "if nargout >= 4," << endl + << " %" << endl + << " % Third order derivatives" << endl + << " %" << endl + << endl; + int ncols = hessianColsNbr * JacobianColsNbr; + if (third_derivatives.size()) + StaticOutput << " v3 = zeros(" << NNZDerivatives[2] << ",3);" << endl + << third_derivatives_output.str() + << " g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");" << endl; + else // Either 3rd derivatives is all zero, or we didn't compute it + StaticOutput << " g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl; + } else { @@ -1305,6 +1391,14 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const << " {" << endl << hessian_output.str() << " }" << endl; + if (third_derivatives.size()) + StaticOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl + << " if (v3 == NULL)" << endl + << " return;" << endl + << " else" << endl + << " {" << endl + << third_derivatives_output.str() + << " }" << endl; } } diff --git a/StaticModel.hh b/StaticModel.hh index f193c4c4..e76a5f0f 100644 --- a/StaticModel.hh +++ b/StaticModel.hh @@ -161,7 +161,7 @@ public: \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed \param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed */ - void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian,bool paramsDerivatives, bool block, bool bytecode); + void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatices, bool paramsDerivatives, bool block, bool bytecode); //! Adds informations for simulation in a binary file for a block decomposed model void Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num,