diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index ce800ac8d..782c38321 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 bbc09f69c..cd39b4bc8 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 a04e0ef10..3107b1b0a 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 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 871909a8b..0e92c8ddf 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 ef8ed7390..c892520be 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 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 483a1b089..e04af03ee 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -126,6 +126,7 @@ int sigma_e = 0; calib {BEGIN DYNARE_STATEMENT; return token::CALIB;} planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;} ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;} +identification {BEGIN DYNARE_STATEMENT; return token::IDENTIFICATION;} bvar_density {BEGIN DYNARE_STATEMENT; return token::BVAR_DENSITY; } bvar_forecast {BEGIN DYNARE_STATEMENT; return token::BVAR_FORECAST; } diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index a8e5a3a08..170ab263c 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 7d1d75fc3..0af95d2b6 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 75919b1cb..fbdef3ea6 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 278623ec2..df156d526 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 e55055c84..9d823945d 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 1f3864d13..a272e9c0e 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 0594bcc2f..35df9360d 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 47a57dbad..4bcba3c54 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 90aa47753..e4cec43b2 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);