diff --git a/ComputingTasks.cc b/ComputingTasks.cc index 70d20bf5..42c9af2e 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -2411,6 +2411,12 @@ CalibSmootherStatement::CalibSmootherStatement(const SymbolList &symbol_list_arg { } +void +CalibSmootherStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) +{ + mod_file_struct.calib_smoother_present = true; +} + void CalibSmootherStatement::writeOutput(ostream &output, const string &basename) const { diff --git a/ComputingTasks.hh b/ComputingTasks.hh index 926e082d..59c6a264 100644 --- a/ComputingTasks.hh +++ b/ComputingTasks.hh @@ -516,6 +516,7 @@ private: public: CalibSmootherStatement(const SymbolList &symbol_list_arg, const OptionsList &options_list_arg); + virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings); virtual void writeOutput(ostream &output, const string &basename) const; }; diff --git a/DynamicModel.cc b/DynamicModel.cc index d19727fa..8c47907c 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -2913,26 +2913,23 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl << "M_.maximum_lag = " << max_lag << ";" << endl << "M_.maximum_lead = " << max_lead << ";" << endl; - if (symbol_table.endo_nbr()) - { - output << "M_.maximum_endo_lag = " << max_endo_lag << ";" << endl - << "M_.maximum_endo_lead = " << max_endo_lead << ";" << endl - << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl; - } - if (symbol_table.exo_nbr()) - { - output << "M_.maximum_exo_lag = " << max_exo_lag << ";" << endl - << "M_.maximum_exo_lead = " << max_exo_lead << ";" << endl - << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl; - } + + output << "M_.maximum_endo_lag = " << max_endo_lag << ";" << endl + << "M_.maximum_endo_lead = " << max_endo_lead << ";" << endl + << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl; + + output << "M_.maximum_exo_lag = " << max_exo_lag << ";" << endl + << "M_.maximum_exo_lead = " << max_exo_lead << ";" << endl + << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl; + if (symbol_table.exo_det_nbr()) { output << "M_.maximum_exo_det_lag = " << max_exo_det_lag << ";" << endl << "M_.maximum_exo_det_lead = " << max_exo_det_lead << ";" << endl << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);" << endl; } - if (symbol_table.param_nbr()) - output << "M_.params = NaN(" << symbol_table.param_nbr() << ", 1);" << endl; + + output << "M_.params = NaN(" << symbol_table.param_nbr() << ", 1);" << endl; // Write number of non-zero derivatives // Use -1 if the derivatives have not been computed @@ -3757,22 +3754,28 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context || symbol_table.getType(it->first.first) == eLogTrend) for (int eq = 0; eq < (int) equations.size(); eq++) { - expr_t testeq = AddLog(AddMinus(equations[eq]->get_arg1(), // F: a = b -> ln(a - b) - equations[eq]->get_arg2())); - testeq = testeq->getDerivative(it->second); // d F / d Trend - for (deriv_id_table_t::const_iterator endogit = deriv_id_table.begin(); - endogit != deriv_id_table.end(); endogit++) - if (symbol_table.getType(endogit->first.first) == eEndogenous) - { - double nearZero = testeq->getDerivative(endogit->second)->eval(eval_context); // eval d F / d Trend d Endog - if (fabs(nearZero) > ZERO_BAND) + expr_t homogeneq = AddMinus(equations[eq]->get_arg1(), + equations[eq]->get_arg2()); + + // Do not run the test if the term inside the log is zero + if (fabs(homogeneq->eval(eval_context)) > ZERO_BAND) + { + expr_t testeq = AddLog(homogeneq); // F = log(lhs-rhs) + testeq = testeq->getDerivative(it->second); // d F / d Trend + for (deriv_id_table_t::const_iterator endogit = deriv_id_table.begin(); + endogit != deriv_id_table.end(); endogit++) + if (symbol_table.getType(endogit->first.first) == eEndogenous) { - cerr << "ERROR: trends not compatible with balanced growth path; the second-order cross partial of equation " << eq + 1 << " w.r.t. trend variable " - << symbol_table.getName(it->first.first) << " and endogenous variable " - << symbol_table.getName(endogit->first.first) << " is not null. " << endl; - exit(EXIT_FAILURE); + double nearZero = testeq->getDerivative(endogit->second)->eval(eval_context); // eval d F / d Trend d Endog + if (fabs(nearZero) > ZERO_BAND) + { + cerr << "ERROR: trends not compatible with balanced growth path; the second-order cross partial of equation " << eq + 1 << " w.r.t. trend variable " + << symbol_table.getName(it->first.first) << " and endogenous variable " + << symbol_table.getName(endogit->first.first) << " is not null. " << endl; + exit(EXIT_FAILURE); + } } - } + } } } @@ -4134,8 +4137,9 @@ DynamicModel::transformPredeterminedVariables() void DynamicModel::detrendEquations() { - for (nonstationary_symbols_map_t::const_iterator it = nonstationary_symbols_map.begin(); - it != nonstationary_symbols_map.end(); it++) + // We go backwards in the list of trend_vars, to deal correctly with I(2) processes + for (nonstationary_symbols_map_t::const_reverse_iterator it = nonstationary_symbols_map.rbegin(); + it != nonstationary_symbols_map.rend(); ++it) for (int i = 0; i < (int) equations.size(); i++) { BinaryOpNode *substeq = dynamic_cast(equations[i]->detrend(it->first, it->second.first, it->second.second)); diff --git a/DynareBison.yy b/DynareBison.yy index 81e3f29d..53c1c767 100644 --- a/DynareBison.yy +++ b/DynareBison.yy @@ -101,7 +101,7 @@ class ParsingDriver; %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME %token FLOAT_NUMBER %token DEFAULT FIXED_POINT -%token FORECAST K_ORDER_SOLVER INSTRUMENTS PRIOR SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN +%token FORECAST K_ORDER_SOLVER INSTRUMENTS SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN %token GAMMA_PDF GRAPH GRAPH_FORMAT CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK STD %token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HOMOTOPY_FORCE_CONTINUE HP_FILTER HP_NGRID HYBRID %token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE BOUNDS JSCALE INIT @@ -113,7 +113,7 @@ class ParsingDriver; %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT %token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN -%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL +%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE %token NAME %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS @@ -151,7 +151,7 @@ class ParsingDriver; %token VLISTLOG VLISTPER %token RESTRICTION RESTRICTION_FNAME CROSS_RESTRICTIONS NLAGS CONTEMP_REDUCED_FORM REAL_PSEUDO_FORECAST %token DUMMY_OBS NSTATES INDXSCALESSTATES NO_BAYESIAN_PRIOR SPECIFICATION SIMS_ZHA -%token ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE +%token ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX %token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD %token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT %token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION @@ -1554,6 +1554,7 @@ estimation_options : o_datafile | o_qz_zero_threshold | o_taper_steps | o_geweke_interval + | o_mcmc_jumping_covariance ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -2673,6 +2674,14 @@ o_discretionary_tol: DISCRETIONARY_TOL EQUAL non_negative_number { driver.option o_analytic_derivation : ANALYTIC_DERIVATION { driver.option_num("analytic_derivation", "1"); } o_endogenous_prior : ENDOGENOUS_PRIOR { driver.option_num("endogenous_prior", "1"); } o_use_univariate_filters_if_singularity_is_detected : USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED EQUAL INT_NUMBER { driver.option_num("use_univariate_filters_if_singularity_is_detected", $3); } +o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN + { driver.option_str("MCMC_jumping_covariance", $3); } | MCMC_JUMPING_COVARIANCE EQUAL PRIOR_VARIANCE + { driver.option_str("MCMC_jumping_covariance", $3); } + | MCMC_JUMPING_COVARIANCE EQUAL IDENTITY_MATRIX + { driver.option_str("MCMC_jumping_covariance", $3); } + | MCMC_JUMPING_COVARIANCE EQUAL filename + { driver.option_str("MCMC_jumping_covariance", $3); } + ; range : symbol ':' symbol { @@ -2791,6 +2800,7 @@ symbol : NAME | FIG | NONE | DR + | PRIOR ; %% diff --git a/DynareFlex.ll b/DynareFlex.ll index f7ea3a9a..66f4edca 100644 --- a/DynareFlex.ll +++ b/DynareFlex.ll @@ -196,7 +196,10 @@ string eofbuff; subsamples {return token::SUBSAMPLES;} options {return token::OPTIONS;} -prior {return token::PRIOR;} +prior { + yylval->string_val = new string(yytext); + return token::PRIOR; +} std {BEGIN DYNARE_STATEMENT; return token::STD;} corr {BEGIN DYNARE_STATEMENT; return token::CORR;} @@ -436,6 +439,19 @@ string eofbuff; random_parameter_convergence_criterion {return token::RANDOM_PARAMETER_CONVERGENCE_CRITERION;} tolf {return token::TOLF;} instruments {return token::INSTRUMENTS;} +hessian { + yylval->string_val = new string(yytext); + return token::HESSIAN; +} +prior_variance { + yylval->string_val = new string(yytext); + return token::PRIOR_VARIANCE; +} +identity_matrix { + yylval->string_val = new string(yytext); + return token::IDENTITY_MATRIX; +} +mcmc_jumping_covariance {return token::MCMC_JUMPING_COVARIANCE;} /* These four (var, varexo, varexo_det, parameters) are for change_type */ var { return token::VAR; } diff --git a/ExprNode.cc b/ExprNode.cc index e3a5bad8..45247ec1 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -21,15 +21,11 @@ #include #include -// For select1st() -#ifdef __GNUC__ -# include -using namespace __gnu_cxx; -#endif - #include #include +#include + #include "ExprNode.hh" #include "DataTree.hh" #include "ModFile.hh" @@ -110,7 +106,7 @@ ExprNode::collectModelLocalVariables(set &result) const set > symb_ids; collectVariables(eModelLocalVariable, symb_ids); transform(symb_ids.begin(), symb_ids.end(), inserter(result, result.begin()), - select1st >()); + boost::bind(&pair::first,_1)); } void diff --git a/ModFile.cc b/ModFile.cc index 4ef0ab43..1abb76ce 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -123,7 +123,8 @@ ModFile::checkPass() || mod_file_struct.estimation_present || mod_file_struct.osr_present || mod_file_struct.ramsey_policy_present - || mod_file_struct.discretionary_policy_present; + || mod_file_struct.discretionary_policy_present + || mod_file_struct.calib_smoother_present; // Allow empty model only when doing a standalone BVAR estimation if (dynamic_model.equation_number() == 0 @@ -318,7 +319,8 @@ ModFile::transformPass(bool nostrict) || mod_file_struct.estimation_present || mod_file_struct.osr_present || mod_file_struct.ramsey_policy_present - || mod_file_struct.discretionary_policy_present) + || mod_file_struct.discretionary_policy_present + || mod_file_struct.calib_smoother_present) { // In stochastic models, create auxiliary vars for leads and lags greater than 2, on both endos and exos dynamic_model.substituteEndoLeadGreaterThanTwo(false); @@ -423,7 +425,8 @@ ModFile::computingPass(bool no_tmp_terms) { if (mod_file_struct.stoch_simul_present || mod_file_struct.estimation_present || mod_file_struct.osr_present - || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present) + || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present + || mod_file_struct.calib_smoother_present) static_model.set_cutoff_to_zero(); const bool static_hessian = mod_file_struct.identification_present @@ -437,7 +440,8 @@ ModFile::computingPass(bool no_tmp_terms) if (mod_file_struct.simul_present || mod_file_struct.check_present || mod_file_struct.stoch_simul_present || mod_file_struct.estimation_present || mod_file_struct.osr_present - || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present) + || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present + || mod_file_struct.calib_smoother_present) { if (mod_file_struct.simul_present) dynamic_model.computingPass(true, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code); @@ -445,7 +449,8 @@ ModFile::computingPass(bool no_tmp_terms) { if (mod_file_struct.stoch_simul_present || mod_file_struct.estimation_present || mod_file_struct.osr_present - || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present) + || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present + || mod_file_struct.calib_smoother_present) dynamic_model.set_cutoff_to_zero(); if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3) { @@ -532,15 +537,20 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b symbol_table.writeOutput(mOutputFile); - // Initialize M_.Sigma_e and M_.H + // Initialize M_.Sigma_e, M_.Correlation_matrix, M_.H, and M_.Correlation_matrix_ME mOutputFile << "M_.Sigma_e = zeros(" << symbol_table.exo_nbr() << ", " + << symbol_table.exo_nbr() << ");" << endl + << "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", " << symbol_table.exo_nbr() << ");" << endl; if (mod_file_struct.calibrated_measurement_errors) mOutputFile << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", " + << symbol_table.observedVariablesNbr() << ");" << endl + << "M_.Correlation_matrix_ME = eye(" << symbol_table.observedVariablesNbr() << ", " << symbol_table.observedVariablesNbr() << ");" << endl; else - mOutputFile << "M_.H = 0;" << endl; + mOutputFile << "M_.H = 0;" << endl + << "M_.Correlation_matrix_ME = 1;" << endl; if (linear == 1) mOutputFile << "options_.linear = 1;" << endl; diff --git a/ModelTree.cc b/ModelTree.cc index 2a798017..362a5622 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -1567,3 +1567,10 @@ ModelTree::computeParamsDerivativesTemporaryTerms() it != hessian_params_derivatives.end(); ++it) it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true); } + +bool ModelTree::isNonstationary(int symb_id) const +{ + return (nonstationary_symbols_map.find(symb_id) + != nonstationary_symbols_map.end()); +} + diff --git a/ModelTree.hh b/ModelTree.hh index 3a29aa90..30afd27d 100644 --- a/ModelTree.hh +++ b/ModelTree.hh @@ -307,6 +307,8 @@ public: void addTrendVariables(vector trend_vars, expr_t growth_factor) throw (TrendException); //! Adds a nonstationary variables with their (common) deflator void addNonstationaryVariables(vector nonstationary_vars, bool log_deflator, expr_t deflator) throw (TrendException); + //! Is a given variable non-stationary? + bool isNonstationary(int symb_id) const; void set_cutoff_to_zero(); //! Helper for writing the Jacobian elements in MATLAB and C /*! Writes either (i+1,j+1) or [i+j*no_eq] */ diff --git a/ParsingDriver.cc b/ParsingDriver.cc index 31d5fbfa..645546df 100644 --- a/ParsingDriver.cc +++ b/ParsingDriver.cc @@ -357,6 +357,13 @@ ParsingDriver::end_nonstationary_var(bool log_deflator, expr_t deflator) { error("Variable " + e.name + " was listed more than once as following a trend."); } + + set > r; + deflator->collectVariables(eEndogenous, r); + for (set >::const_iterator it = r.begin(); it != r.end(); ++it) + if (dynamic_model->isNonstationary(it->first)) + error("The deflator contains a non-stationary endogenous variable. This is not allowed. Please use only stationary endogenous and/or {log_}trend_vars."); + declared_nonstationary_vars.clear(); reset_data_tree(); } diff --git a/Shocks.cc b/Shocks.cc index ab5b36df..1a943d2c 100644 --- a/Shocks.cc +++ b/Shocks.cc @@ -155,17 +155,19 @@ ShocksStatement::writeCovarOrCorrShock(ostream &output, covar_and_corr_shocks_t: SymbolType type2 = symbol_table.getType(it->first.second); assert((type1 == eExogenous && type2 == eExogenous) || (symbol_table.isObservedVariable(it->first.first) && symbol_table.isObservedVariable(it->first.second))); - string matrix; + string matrix, corr_matrix; int id1, id2; if (type1 == eExogenous) { matrix = "M_.Sigma_e"; + corr_matrix = "M_.Correlation_matrix"; id1 = symbol_table.getTypeSpecificID(it->first.first) + 1; id2 = symbol_table.getTypeSpecificID(it->first.second) + 1; } else { matrix = "M_.H"; + corr_matrix = "M_.Correlation_matrix_ME"; id1 = symbol_table.getObservedVariableIndex(it->first.first) + 1; id2 = symbol_table.getObservedVariableIndex(it->first.second) + 1; } @@ -178,6 +180,15 @@ ShocksStatement::writeCovarOrCorrShock(ostream &output, covar_and_corr_shocks_t: output << ";" << endl << matrix << "(" << id2 << ", " << id1 << ") = " << matrix << "(" << id1 << ", " << id2 << ");" << endl; + + if (corr) + { + output << corr_matrix << "(" << id1 << ", " << id2 << ") = "; + it->second->writeOutput(output); + output << ";" << endl + << corr_matrix << "(" << id2 << ", " << id1 << ") = " + << corr_matrix << "(" << id1 << ", " << id2 << ");" << endl; + } } void diff --git a/Statement.cc b/Statement.cc index 76b600b7..53d06911 100644 --- a/Statement.cc +++ b/Statement.cc @@ -47,7 +47,8 @@ ModFileStructure::ModFileStructure() : dsge_var_estimated(false), bayesian_irf_present(false), estimation_data_statement_present(false), - last_markov_switching_chain(0) + last_markov_switching_chain(0), + calib_smoother_present(false) { } diff --git a/Statement.hh b/Statement.hh index cd3fd5cf..ff0a53e0 100644 --- a/Statement.hh +++ b/Statement.hh @@ -97,6 +97,8 @@ public: bool estimation_data_statement_present; //! Last chain number for Markov Switching statement int last_markov_switching_chain; + //! Whether a calib_smoother statement is present + bool calib_smoother_present; }; class Statement