From b775f1812c9bae7c6db1e1453f983716ce9b5dc5 Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Mon, 20 Apr 2009 13:58:15 +0000
Subject: [PATCH] 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
---
preprocessor/ComputingTasks.cc | 15 +++++
preprocessor/ComputingTasks.hh | 9 +++
preprocessor/DynamicModel.cc | 100 ++++++++++++++++++++++++++++++++-
preprocessor/DynamicModel.hh | 19 ++++++-
preprocessor/DynareBison.yy | 25 +++++----
preprocessor/DynareFlex.ll | 1 +
preprocessor/ExprNode.cc | 9 +--
preprocessor/ModFile.cc | 6 +-
preprocessor/ModelTree.cc | 17 +++---
preprocessor/ModelTree.hh | 2 +-
preprocessor/ParsingDriver.cc | 6 ++
preprocessor/ParsingDriver.hh | 1 +
preprocessor/Statement.cc | 3 +-
preprocessor/Statement.hh | 2 +
preprocessor/StaticModel.cc | 2 +-
15 files changed, 183 insertions(+), 34 deletions(-)
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index ce800ac8da..782c38321d 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/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
+{
+}
diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh
index bbc09f69c8..cd39b4bc84 100644
--- a/preprocessor/ComputingTasks.hh
+++ b/preprocessor/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
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index a04e0ef106..3107b1b0a1 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/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();
+}
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index 871909a8b8..0e92c8ddfc 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/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;
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index ef8ed73900..c892520bed 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/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(); }
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 483a1b0890..e04af03ee6 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/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; }
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index a8e5a3a087..170ab263c4 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/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:
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 7d1d75fc36..0af95d2b69 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/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
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 75919b1cba..fbdef3ea68 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/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);
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index 278623ec2c..df156d5262 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/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;
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index e55055c846..9d823945df 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/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)
{
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 1f3864d13c..a272e9c0e9 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/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
diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc
index 0594bcc2f3..35df9360db 100644
--- a/preprocessor/Statement.cc
+++ b/preprocessor/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)
{
}
diff --git a/preprocessor/Statement.hh b/preprocessor/Statement.hh
index 47a57dbadb..4bcba3c544 100644
--- a/preprocessor/Statement.hh
+++ b/preprocessor/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
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 90aa477532..e4cec43b2a 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/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);
--
GitLab