diff --git a/parser.src/ComputingTasks.cc b/parser.src/ComputingTasks.cc
index 88020ee24199c76cfed0e076fd216d221c80ee33..85be18e6571ee912a8e055fab1a56aaa346c2ff8 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 cd62be38b057af91034b12e11d6db5e9f5a1528e..50eaddffb09c4134f3cd0e395b4dd0b37bd45bf9 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 4c0fc49c16cbc50250dc30d2b781d4c9401afd8d..d340717eaa8dfb17cfaadb441a652cab874b9b59 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<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
diff --git a/parser.src/Statement.cc b/parser.src/Statement.cc
index bdc4bb492b9c099b2549bf3f1f7f65578932214d..a92282d5b889dd7ddb11bb19fecfb2deccf167cc 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 edfdd8181cefc5b8637ebfc7fa4af3735c0aaeca..180cd216ccdf8c4c6dc228c5155e7067733eef03 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 c95157fef0c31a0f5244386d271128e93503b021..1940850edf2b4ba1427b38dbb9fc5826436cd309 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<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;
diff --git a/parser.src/include/Statement.hh b/parser.src/include/Statement.hh
index f71e6975caf2da5b4a922c11b8f7d08aee961027..e1f98fa75f34d21145f13b2457fdb09c2fe91a04 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