From 3c79e231c7dc2226e580466cd288ae3574f31090 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:41:56 +0100 Subject: [PATCH 1/8] Removed penalty_hessian routine. + Code factorization. + Added an option for using the penalized objective when computing numerically the hessian at the mode. Previous behaviour (introduced with penalty_hessian routine) was to compute the hessian matrix at the mode with the penalized objective function (instead of the original objective function). This behaviour hides problematic situations, where the computed hessian (using the original objective) would not be full rank. For instance, if the estimation ends up with a parameter on (or very close to) the bounds of its possible values (which is often not a desirable outcome), the estimated posterior variance would be zero for this parameter (with the original objective) because the hessian is not finite in this direction, while the posterior variance would be positive if the penalized objective is used instead. But this estimate would not be reliable by construction of the penalty which is quite ad-hoc (more fundamentally I do not think that there exists any rational for approximating the covariance matrix with the inverse of the hessian matrix if the mode is on the border of the set of possible values). This commit restore the behaviour previous to 2446ab02ba4b3ed88c9c5021aced076078d96007. An option is available for computing the hessian with the penalized objective function. --- DynareBison.yy | 3 +++ DynareFlex.ll | 1 + 2 files changed, 4 insertions(+) diff --git a/DynareBison.yy b/DynareBison.yy index bbae9875..7f6d2883 100644 --- a/DynareBison.yy +++ b/DynareBison.yy @@ -113,6 +113,7 @@ class ParsingDriver; %token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW %token FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION %token NAME +%token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY %token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS %token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE @@ -1825,6 +1826,7 @@ estimation_options : o_datafile | o_posterior_sampling_method | o_posterior_sampler_options | o_keep_kalman_algo_if_singularity_is_detected + | o_use_penalized_objective_for_hessian ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -3232,6 +3234,7 @@ o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN | MCMC_JUMPING_COVARIANCE EQUAL filename { driver.option_str("MCMC_jumping_covariance", $3); } ; +o_use_penalized_objective_for_hessian : USE_PENALIZED_OBJECTIVE_FOR_HESSIAN { driver.option_num("hessian.use_penalized_objective","true"); }; o_irf_plot_threshold : IRF_PLOT_THRESHOLD EQUAL non_negative_number { driver.option_num("impulse_responses.plot_threshold", $3); }; o_dr_display_tol : DR_DISPLAY_TOL EQUAL non_negative_number { driver.option_num("dr_display_tol", $3); }; o_consider_all_endogenous : CONSIDER_ALL_ENDOGENOUS { driver.option_str("endo_vars_for_moment_computations_in_estimation", "all_endogenous_variables"); }; diff --git a/DynareFlex.ll b/DynareFlex.ll index 8d00b287..e1e0aa4c 100644 --- a/DynareFlex.ll +++ b/DynareFlex.ll @@ -400,6 +400,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;} proposal_distribution {return token::PROPOSAL_DISTRIBUTION;} no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;} +use_penalized_objective_for_hessian {return token::USE_PENALIZED_OBJECTIVE_FOR_HESSIAN;} alpha { yylval->string_val = new string(yytext); From 5e1d20c8d94193c7e6a2e0bbe4692cb90e025072 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 12:46:37 +0100 Subject: [PATCH 2/8] preprocessor: replace exit(1) with exit(EXIT_FAILURE) --- ComputingTasks.cc | 4 ++-- DynamicModel.cc | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ComputingTasks.cc b/ComputingTasks.cc index f12c4086..ca83e94e 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2016 Dynare Team + * Copyright (C) 2003-2017 Dynare Team * * This file is part of Dynare. * @@ -336,7 +336,7 @@ RamseyConstraintsStatement::writeOutput(ostream &output, const string &basename, break; default: cerr << "Ramsey constraints: this shouldn't happen." << endl; - exit(1); + exit(EXIT_FAILURE); } output << "', '"; it->expression->writeOutput(output); diff --git a/DynamicModel.cc b/DynamicModel.cc index 08dd2bc4..21b18d7a 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -5012,7 +5012,7 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons break; default: std::cerr << "This case shouldn't happen" << std::endl; - exit(1); + exit(EXIT_FAILURE); } derivative deriv(col_id + eq*cols_nbr,col_id,eq,it->second); D.push_back(deriv); From d6e801cac80409b05b1563be9d15fcc9e5c18303 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 16:21:18 +0100 Subject: [PATCH 3/8] =?UTF-8?q?preprocessor:=20change=20error=20to=20warni?= =?UTF-8?q?ng=20because=20we=20don=E2=80=99t=20end=20processing=20in=20thi?= =?UTF-8?q?s=20situation;=20fix=20error=20message?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ComputingTasks.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ComputingTasks.cc b/ComputingTasks.cc index ca83e94e..9843ba7d 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -1058,7 +1058,7 @@ ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, output << "';" << endl; } else - cout << "Error : Non-variable symbol used in TREND_COEFF: " << it->first << endl; + cerr << "Warning : Non-variable symbol used in observation_trends: " << it->first << endl; } } From 4dc640c635d75491e6c541d98848184a6f0d6632 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 16:22:00 +0100 Subject: [PATCH 4/8] preprocessor: aesthetic fix --- ComputingTasks.cc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/ComputingTasks.cc b/ComputingTasks.cc index 9843ba7d..ae580372 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -1044,10 +1044,7 @@ void ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const { output << "options_.trend_coeff = {};" << endl; - - trend_elements_t::const_iterator it; - - for (it = trend_elements.begin(); it != trend_elements.end(); it++) + for (trend_elements_t::const_iterator it = trend_elements.begin(); it != trend_elements.end(); it++) { SymbolType type = symbol_table.getType(it->first); if (type == eEndogenous) From 83ea2691556bfb8271745b10039296e940495fc9 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 10 Feb 2017 11:07:45 +0100 Subject: [PATCH 5/8] preprocessor: aesthetic fix --- ComputingTasks.cc | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/ComputingTasks.cc b/ComputingTasks.cc index ae580372..a05727d3 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -821,9 +821,7 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename, b << "estim_params_.corrn = [];" << endl << "estim_params_.param_vals = [];" << endl; - vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); @@ -892,9 +890,7 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam if (use_calibration) output << "options_.use_calibration_initialization = 1;" << endl; - vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); @@ -955,9 +951,7 @@ EstimatedParamsBoundsStatement::EstimatedParamsBoundsStatement(const vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); From ad849141186d97513b11df0714bf7143e90c2c75 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 24 Feb 2017 11:20:54 +0100 Subject: [PATCH 6/8] preprocessor: for consistency, use equations.size() instead of equation_number() in classes that have equations as a field --- DynamicModel.cc | 12 ++++++------ ModelTree.cc | 20 ++++++++++---------- StaticModel.cc | 8 ++++---- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/DynamicModel.cc b/DynamicModel.cc index 21b18d7a..ebf6b9b0 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -3252,8 +3252,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative if (!no_tmp_terms) computeTemporaryTermsOrdered(); int k = 0; - equation_block = vector(equation_number()); - variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number()); + equation_block = vector(equations.size()); + variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equations.size()); for (unsigned int i = 0; i < getNbBlocks(); i++) { for (unsigned int j = 0; j < getBlockSize(i); j++) @@ -4210,10 +4210,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "% from model file (.mod)" << endl << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "if nargout >= 3" << endl << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl @@ -4238,10 +4238,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "ss_param_deriv, ss_param_2nd_deriv)" << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl << hessian1_output.str() diff --git a/ModelTree.cc b/ModelTree.cc index 70db28e0..1faac8ff 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -36,7 +36,7 @@ using namespace MFS; bool ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose) { - const int n = equation_number(); + const int n = equations.size(); assert(n == symbol_table.endo_nbr()); @@ -96,7 +96,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo #endif // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them - endo2eq.resize(equation_number()); + endo2eq.resize(equations.size()); transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus(), n)); #ifdef DEBUG @@ -143,7 +143,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian cout << "Normalizing the model..." << endl; - int n = equation_number(); + int n = equations.size(); // compute the maximum value of each row of the contemporaneous Jacobian matrix //jacob_map normalized_contemporaneous_jacobian; @@ -329,7 +329,7 @@ ModelTree::writeRevXrefs(ostream &output, const map > &xrefmap, co void ModelTree::computeNormalizedEquations(multimap &endo2eqs) const { - for (int i = 0; i < equation_number(); i++) + for (int i = 0; i < equations.size(); i++) { VariableNode *lhs = dynamic_cast(equations[i]->get_arg1()); if (lhs == NULL) @@ -417,11 +417,11 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m void ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector &equation_reordered, vector &variable_reordered) { - vector eq2endo(equation_number(), 0); - equation_reordered.resize(equation_number()); - variable_reordered.resize(equation_number()); + vector eq2endo(equations.size(), 0); + equation_reordered.resize(equations.size()); + variable_reordered.resize(equations.size()); bool *IM; - int n = equation_number(); + int n = equations.size(); IM = (bool *) calloc(n*n, sizeof(bool)); int i = 0; for (vector::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++) @@ -1697,7 +1697,7 @@ ModelTree::addEquation(expr_t eq, int lineno) void ModelTree::addEquation(expr_t eq, int lineno, vector > &eq_tags) { - int n = equation_number(); + int n = equations.size(); for (size_t i = 0; i < eq_tags.size(); i++) equation_tags.push_back(make_pair(n, eq_tags[i])); addEquation(eq, lineno); @@ -1741,7 +1741,7 @@ ModelTree::addNonstationaryVariables(vector nonstationary_vars, bool log_de void ModelTree::initializeVariablesAndEquations() { - for (int j = 0; j < equation_number(); j++) + for (int j = 0; j < equations.size(); j++) { equation_reordered.push_back(j); variable_reordered.push_back(j); diff --git a/StaticModel.cc b/StaticModel.cc index a7652f9b..2cdb2f0a 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -2368,10 +2368,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons << "% from model file (.mod)" << endl << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", " + << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "if nargout >= 3" << endl @@ -2396,10 +2396,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons << "function params_derivs(y, x, params)" << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", " + << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl From 9eb7c00da10f965e44ef51340981cda1c884fb7c Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 1 Mar 2017 12:20:55 +0100 Subject: [PATCH 7/8] preprocessor: fix bug that caused temporary terms not to be printed for the parameter derivatives of the static model. closes #1397 --- StaticModel.hh | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/StaticModel.hh b/StaticModel.hh index bf6cad6f..8535a714 100644 --- a/StaticModel.hh +++ b/StaticModel.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2016 Dynare Team + * Copyright (C) 2003-2017 Dynare Team * * This file is part of Dynare. * @@ -30,9 +30,6 @@ using namespace std; class StaticModel : public ModelTree { private: - //! Temporary terms for the file containing parameters dervicatives - temporary_terms_t params_derivs_temporary_terms; - //! global temporary terms for block decomposed models vector > v_temporary_terms; From fc48bf26b9b5c6f231e3dc1cd9c7451568984369 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 1 Mar 2017 12:44:31 +0100 Subject: [PATCH 8/8] preprocessor: write error messages to cerr instead of cout, replace \n with endl --- DynamicModel.cc | 6 +++--- ModelTree.cc | 2 +- StaticModel.cc | 6 +++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/DynamicModel.cc b/DynamicModel.cc index ebf6b9b0..f3c90c25 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -800,7 +800,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } @@ -1076,7 +1076,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } //Temporary variables declaration @@ -1768,7 +1768,7 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const SaveCode.open((bin_basename + "_dynamic.bin").c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0; diff --git a/ModelTree.cc b/ModelTree.cc index 1faac8ff..f295423e 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -1594,7 +1594,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &basename, SaveCode.open(bin_basename.c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0; diff --git a/StaticModel.cc b/StaticModel.cc index 2cdb2f0a..1d176027 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -412,7 +412,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } int count_u; @@ -596,7 +596,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } //Temporary variables declaration @@ -990,7 +990,7 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const st SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0;