diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index b56685497b02c2832f03961c34ac8d0244d1ce09..ae7649105a0cfe86e8b44b55a792c1e7b5817865 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/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/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 94f5a1a59e22875fa67fe2158571e0718ab90ae6..2c44e3f2434c96d53e74db85772b701f7b4ad7ee 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/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/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index c78324ea87acf11b1d31dbb9a3839bf9c5b3aed9..d82efaeda73153ad657225ee293940cd08cafa74 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/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<int> 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<int>::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/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index f193c4c4dc7f7ce1dd75fdc1cd1ce1ea76f51bd9..e76a5f0fe3d1a0c1ac4a0b1c7ae9775c54c50f98 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/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,