diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 07f7c893..f6fd184e 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -3648,100 +3648,99 @@ DynamicModel::getUndiffLHSForPac(const string &aux_model_name, } void -DynamicModel::analyzePacEquationStructure(const string &name, map, pair> &eqtag_and_lag, PacModelTable::equation_info_t &pac_equation_info) +DynamicModel::analyzePacEquationStructure(const string &name, map &pac_eq_name, PacModelTable::equation_info_t &pac_equation_info) { - int i = 0; for (auto &equation : equations) - { - pair lhs(-1, -1); - pair>> ec_params_and_vars; - vector> ar_params_and_vars; - vector> non_optim_vars_params_and_constants, optim_additive_vars_params_and_constants, additive_vars_params_and_constants; + if (equation->containsPacExpectation(name)) + { + if (!pac_eq_name[name].empty()) + { + cerr << "It is not possible to use 'pac_expectation(" << name << ")' in several equations." << endl; + exit(EXIT_FAILURE); + } + string eqn = equation_tags.getTagValueByEqnAndKey(&equation - &equations[0], "name"); + if (eqn.empty()) + { + cerr << "Every equation with a 'pac_expectation' operator must have been assigned an equation tag name" << endl; + exit(EXIT_FAILURE); + } + pac_eq_name[name] = eqn; - if (equation->containsPacExpectation(name)) - { - set> lhss; - equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss); - lhs = *lhss.begin(); - int lhs_symb_id = lhs.first; - int lhs_orig_symb_id = lhs_symb_id; - if (symbol_table.isAuxiliaryVariable(lhs_orig_symb_id)) + set> lhss; + equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss); + auto lhs = *lhss.begin(); + int lhs_symb_id = lhs.first; + int lhs_orig_symb_id = lhs_symb_id; + if (symbol_table.isAuxiliaryVariable(lhs_orig_symb_id)) + try + { + lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_orig_symb_id); + } + catch (...) + { + } + + auto arg2 = dynamic_cast(equation->arg2); + if (!arg2) + { + cerr << "Pac equation in incorrect format" << endl; + exit(EXIT_FAILURE); + } + auto [optim_share_index, optim_part, non_optim_part, additive_part] + = arg2->getPacOptimizingShareAndExprNodes(lhs_symb_id, lhs_orig_symb_id); + + pair>> ec_params_and_vars; + vector> ar_params_and_vars; + vector> non_optim_vars_params_and_constants, optim_additive_vars_params_and_constants, additive_vars_params_and_constants; + if (!optim_part) + { + auto bopn = dynamic_cast(equation->arg2); + if (!bopn) + { + cerr << "Error in PAC equation" << endl; + exit(EXIT_FAILURE); + } + bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars, additive_vars_params_and_constants); + } + else + { + auto bopn = dynamic_cast(optim_part); + if (!bopn) + { + cerr << "Error in PAC equation" << endl; + exit(EXIT_FAILURE); + } + bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars, optim_additive_vars_params_and_constants); try { - lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_orig_symb_id); + non_optim_vars_params_and_constants = non_optim_part->matchLinearCombinationOfVariables(); + if (additive_part) + additive_vars_params_and_constants = additive_part->matchLinearCombinationOfVariables(); } - catch (...) + catch (ExprNode::MatchFailureException &e) { + cerr << "Error in parsing non-optimizing agents or additive part of PAC equation: " + << e.message << endl; + exit(EXIT_FAILURE); } + } - auto arg2 = dynamic_cast(equation->arg2); - if (!arg2) - { - cerr << "Pac equation in incorrect format" << endl; - exit(EXIT_FAILURE); - } - auto [optim_share_index, optim_part, non_optim_part, additive_part] - = arg2->getPacOptimizingShareAndExprNodes(lhs_symb_id, lhs_orig_symb_id); - - if (!optim_part) - { - auto bopn = dynamic_cast(equation->arg2); - if (!bopn) - { - cerr << "Error in PAC equation" << endl; - exit(EXIT_FAILURE); - } - bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars, additive_vars_params_and_constants); - } - else - { - auto bopn = dynamic_cast(optim_part); - if (!bopn) - { - cerr << "Error in PAC equation" << endl; - exit(EXIT_FAILURE); - } - bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars, optim_additive_vars_params_and_constants); - try - { - non_optim_vars_params_and_constants = non_optim_part->matchLinearCombinationOfVariables(); - if (additive_part) - additive_vars_params_and_constants = additive_part->matchLinearCombinationOfVariables(); - } - catch (ExprNode::MatchFailureException &e) - { - cerr << "Error in parsing non-optimizing agents or additive part of PAC equation: " - << e.message << endl; - exit(EXIT_FAILURE); - } - } - - string eqtag = equation_tags.getTagValueByEqnAndKey(&equation - &equations[0], "name"); - if (eqtag.empty()) - { - cerr << "Every equation with a pac expectation must have been assigned an equation tag name" << endl; - exit(EXIT_FAILURE); - } - if (lhs.first == -1) - { - cerr << "analyzePacEquationStructure: error obtaining LHS variable." << endl; - exit(EXIT_FAILURE); - } - if (ec_params_and_vars.second.empty()) - { - cerr << "analyzePacEquationStructure: error obtaining RHS parameters." << endl; - exit(EXIT_FAILURE); - } - string eq = "eq" + to_string(i++); - pac_equation_info[{name, eq}] = {lhs, optim_share_index, - ar_params_and_vars, ec_params_and_vars, - non_optim_vars_params_and_constants, - additive_vars_params_and_constants, - optim_additive_vars_params_and_constants}; - int max_lag = ar_params_and_vars.size(); - eqtag_and_lag[{name, eqtag}] = {eq, max_lag}; - } - } + if (lhs.first == -1) + { + cerr << "analyzePacEquationStructure: error obtaining LHS variable." << endl; + exit(EXIT_FAILURE); + } + if (ec_params_and_vars.second.empty()) + { + cerr << "analyzePacEquationStructure: error obtaining RHS parameters." << endl; + exit(EXIT_FAILURE); + } + pac_equation_info[name] = {lhs, optim_share_index, + ar_params_and_vars, ec_params_and_vars, + non_optim_vars_params_and_constants, + additive_vars_params_and_constants, + optim_additive_vars_params_and_constants}; + } } int @@ -3785,10 +3784,10 @@ DynamicModel::getPacTargetSymbId(const string &pac_model_name) const void DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int discount_symb_id, - const map, pair> &eqtag_and_lag, + int pac_eq_max_lag, ExprNode::subst_table_t &diff_subst_table, - map, int> &pac_mce_z1_symb_ids, - map, vector> &pac_mce_alpha_symb_ids) + map &pac_mce_z1_symb_ids, + map> &pac_mce_alpha_symb_ids) { int pac_target_symb_id; try @@ -3801,152 +3800,131 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d exit(EXIT_FAILURE); } int neqs = 0; - for (const auto &[key, val] : eqtag_and_lag) + + // Create the endogenous representing Z₁ + int mce_z1_symb_id; + string varname = "mce_Z1_" + name; + try { - auto &[name2, eqtag] = key; - if (name2 != name) - continue; + mce_z1_symb_id = symbol_table.addSymbol(varname, SymbolType::endogenous); + } + catch (SymbolTable::AlreadyDeclaredException &e) + { + cerr << "The variable/parameter '" << varname << "' conflicts with the variable that will be generated for the Z_1 variable of the '" << name << "' PAC model. Please rename it." << endl; + exit(EXIT_FAILURE); + } + pac_mce_z1_symb_ids[name] = mce_z1_symb_id; - auto &standard_eqtag = val.first; - int pac_max_lag_m = val.second + 1; - string append_to_name = name + "_" + standard_eqtag; - - // Create the endogenous representing Z₁ - int mce_z1_symb_id; - string varname = "mce_Z1_" + append_to_name; + expr_t A = One; + expr_t fp = Zero; + expr_t beta = AddVariable(discount_symb_id); + for (int i = 1; i <= pac_eq_max_lag+1; i++) + { + string param_name = "mce_alpha_" + name + "_" + to_string(i); try { - mce_z1_symb_id = symbol_table.addSymbol(varname, SymbolType::endogenous); + int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter); + pac_mce_alpha_symb_ids[name].push_back(alpha_i_symb_id); + A = AddPlus(A, AddVariable(alpha_i_symb_id)); + fp = AddPlus(fp, + AddTimes(AddTimes(AddVariable(alpha_i_symb_id), + AddPower(beta, AddPossiblyNegativeConstant(i))), + AddVariable(mce_z1_symb_id, i))); } catch (SymbolTable::AlreadyDeclaredException &e) { - cerr << "The variable/parameter '" << varname << "' conflicts with the variable that will be generated for the Z_1 variable of the '" << name << "' PAC model. Please rename it." << endl; + cerr << "The variable/parameter '" << param_name << "' conflicts with a parameter that will be generated for the '" << name << "' PAC model. Please rename it." << endl; exit(EXIT_FAILURE); } - pac_mce_z1_symb_ids[{name, standard_eqtag}] = mce_z1_symb_id; + } - expr_t A = One; - expr_t fp = Zero; - expr_t beta = AddVariable(discount_symb_id); - for (int i = 1; i <= pac_max_lag_m; i++) - { - string param_name = "mce_alpha_" + append_to_name + "_" + to_string(i); - try - { - int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter); - pac_mce_alpha_symb_ids[{name, standard_eqtag}].push_back(alpha_i_symb_id); - A = AddPlus(A, AddVariable(alpha_i_symb_id)); - fp = AddPlus(fp, - AddTimes(AddTimes(AddVariable(alpha_i_symb_id), - AddPower(beta, AddPossiblyNegativeConstant(i))), - AddVariable(mce_z1_symb_id, i))); - } - catch (SymbolTable::AlreadyDeclaredException &e) - { - cerr << "The variable/parameter '" << param_name << "' conflicts with a parameter that will be generated for the '" << name << "' PAC model. Please rename it." << endl; - exit(EXIT_FAILURE); - } - } + // Add diff nodes and eqs for pac_target_symb_id + const VariableNode *target_base_diff_node; + auto create_target_lag = [&](int lag) + { + if (symbol_table.isAuxiliaryVariable(pac_target_symb_id)) + { + // We know it is a log, see ExprNode::matchParamTimesTargetMinusVariable() + /* We don’t use SymbolTable::getOrigSymbIdForAuxVar(), because it + does not work for unary ops, and changing this behaviour might + break stuff that relies on an exception in this case. */ + auto avi = symbol_table.getAuxVarInfo(pac_target_symb_id); + return AddLog(AddVariable(avi.get_orig_symb_id(), lag)); + } + else + return dynamic_cast(AddVariable(pac_target_symb_id, lag)); + }; - // Add diff nodes and eqs for pac_target_symb_id - const VariableNode *target_base_diff_node; - auto create_target_lag = [&](int lag) - { - if (symbol_table.isAuxiliaryVariable(pac_target_symb_id)) - { - // We know it is a log, see ExprNode::matchParamTimesTargetMinusVariable() - /* We don’t use SymbolTable::getOrigSymbIdForAuxVar(), because it - does not work for unary ops, and changing this behaviour might - break stuff that relies on an exception in this case. */ - auto avi = symbol_table.getAuxVarInfo(pac_target_symb_id); - return AddLog(AddVariable(avi.get_orig_symb_id(), lag)); - } - else - return dynamic_cast(AddVariable(pac_target_symb_id, lag)); - }; - - expr_t diff_node_to_search = AddDiff(create_target_lag(0)); - if (auto sit = diff_subst_table.find(diff_node_to_search); - sit != diff_subst_table.end()) - target_base_diff_node = sit->second; - else - { - int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search); - target_base_diff_node = AddVariable(symb_id); - auto neweq = AddEqual(const_cast(target_base_diff_node), - AddMinus(create_target_lag(0), - create_target_lag(-1))); - addEquation(neweq, -1); - addAuxEquation(neweq); - neqs++; - } - - map target_aux_var_to_add; - const VariableNode *last_aux_var = target_base_diff_node; - for (int i = 1; i <= pac_max_lag_m - 1; i++, neqs++) - { - expr_t this_diff_node = AddDiff(create_target_lag(i)); - int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node, - last_aux_var->symb_id, last_aux_var->lag); - VariableNode *current_aux_var = AddVariable(symb_id); - auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1)); - addEquation(neweq, -1); - addAuxEquation(neweq); - last_aux_var = current_aux_var; - target_aux_var_to_add[i] = current_aux_var; - } - - expr_t fs = Zero; - for (int k = 1; k <= pac_max_lag_m - 1; k++) - { - expr_t ssum = Zero; - for (int j = k+1; j <= pac_max_lag_m; j++) - { - int alpha_j_symb_id = -1; - string param_name = "mce_alpha_" + append_to_name + "_" + to_string(j); - try - { - alpha_j_symb_id = symbol_table.getID(param_name); - } - catch (SymbolTable::UnknownSymbolNameException &e) - { - alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter); - } - ssum = AddPlus(ssum, - AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j)))); - } - fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k])); - } - auto neweq = AddEqual(AddVariable(mce_z1_symb_id), - AddMinus(AddTimes(A, AddMinus(const_cast(target_base_diff_node), fs)), fp)); + expr_t diff_node_to_search = AddDiff(create_target_lag(0)); + if (auto sit = diff_subst_table.find(diff_node_to_search); + sit != diff_subst_table.end()) + target_base_diff_node = sit->second; + else + { + int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search); + target_base_diff_node = AddVariable(symb_id); + auto neweq = AddEqual(const_cast(target_base_diff_node), + AddMinus(create_target_lag(0), + create_target_lag(-1))); addEquation(neweq, -1); + addAuxEquation(neweq); neqs++; } + + map target_aux_var_to_add; + const VariableNode *last_aux_var = target_base_diff_node; + for (int i = 1; i <= pac_eq_max_lag; i++, neqs++) + { + expr_t this_diff_node = AddDiff(create_target_lag(i)); + int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node, + last_aux_var->symb_id, last_aux_var->lag); + VariableNode *current_aux_var = AddVariable(symb_id); + auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1)); + addEquation(neweq, -1); + addAuxEquation(neweq); + last_aux_var = current_aux_var; + target_aux_var_to_add[i] = current_aux_var; + } + + expr_t fs = Zero; + for (int k = 1; k <= pac_eq_max_lag; k++) + { + expr_t ssum = Zero; + for (int j = k+1; j <= pac_eq_max_lag+1; j++) + { + int alpha_j_symb_id = -1; + string param_name = "mce_alpha_" + name + "_" + to_string(j); + try + { + alpha_j_symb_id = symbol_table.getID(param_name); + } + catch (SymbolTable::UnknownSymbolNameException &e) + { + alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter); + } + ssum = AddPlus(ssum, + AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j)))); + } + fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k])); + } + auto neweq = AddEqual(AddVariable(mce_z1_symb_id), + AddMinus(AddTimes(A, AddMinus(const_cast(target_base_diff_node), fs)), fp)); + addEquation(neweq, -1); + neqs++; + cout << "PAC Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl; } void DynamicModel::computePacModelConsistentExpectationSubstitution(const string &name, - const map, pair> &eqtag_and_lag, - expr_t growth, - const map &pac_growth_neutrality_params, - const map, int> &pac_mce_z1_symb_ids, - map, expr_t> &pac_expectation_substitution) + expr_t growth_correction_term, + int pac_mce_z1_symb_id, + map &pac_expectation_substitution) { - for (const auto &[key, val] : eqtag_and_lag) - { - auto &[name2, eqtag] = key; - if (name2 != name) - continue; - - int mce_z1_symb_id = pac_mce_z1_symb_ids.at({name, val.first}); - - /* The growth correction term is not added to the definition of Z₁ - because the latter is recursive. Rather put it at the level of the - substition of pac_expectation operator. */ - expr_t growth_correction = growth ? AddTimes(AddVariable(pac_growth_neutrality_params.at(name)), growth) : Zero; - pac_expectation_substitution[{name, eqtag}] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction); - } + /* The growth correction term is not added to the definition of Z₁ + because the latter is recursive. Rather put it at the level of the + substition of pac_expectation operator. */ + pac_expectation_substitution[name] = AddPlus(AddVariable(pac_mce_z1_symb_id), growth_correction_term); } void @@ -3954,13 +3932,11 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name, const vector &lhs, int max_lag, const string &aux_model_type, - const map, pair> &eqtag_and_lag, const vector &nonstationary, - expr_t growth, - const map &pac_growth_neutrality_params, - map, vector> &pac_h0_indices, - map, vector> &pac_h1_indices, - map, expr_t> &pac_expectation_substitution) + expr_t growth_correction_term, + map> &pac_h0_indices, + map> &pac_h1_indices, + map &pac_expectation_substitution) { bool stationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), logical_not()); bool nonstationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), [](bool b) { return b; }); // FIXME: use std::identity instead of an anonymous function when we upgrade to C++20 @@ -3978,77 +3954,63 @@ DynamicModel::computePacBackwardExpectationSubstitution(const string &name, } }; - for (const auto &[key, val] : eqtag_and_lag) + expr_t subExpr = Zero; + if (stationary_vars_present) { - auto &[name2, eqtag] = key; - if (name2 != name) - continue; - auto &standard_eqtag = val.first; - expr_t subExpr = Zero; - if (stationary_vars_present) + if (aux_model_type == "var") { - if (aux_model_type == "var") - { - /* If the auxiliary model is a VAR, add a parameter corresponding - to the constant. */ - int new_param_symb_id = create_aux_param("h0_" + name + "_" + standard_eqtag + "_constant"); - pac_h0_indices[{name, standard_eqtag}].push_back(new_param_symb_id); - subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id)); - } - for (int i = 1; i < max_lag + 1; i++) - for (auto lhsit : lhs) - { - int new_param_symb_id = create_aux_param("h0_" + name + "_" + standard_eqtag - + "_var_" + symbol_table.getName(lhsit) - + "_lag_" + to_string(i)); - pac_h0_indices[{name, standard_eqtag}].push_back(new_param_symb_id); - subExpr = AddPlus(subExpr, - AddTimes(AddVariable(new_param_symb_id), - AddVariable(lhsit, -i))); - } + /* If the auxiliary model is a VAR, add a parameter corresponding + to the constant. */ + int new_param_symb_id = create_aux_param("h0_" + name + "_constant"); + pac_h0_indices[name].push_back(new_param_symb_id); + subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id)); } - - if (nonstationary_vars_present) - { - if (aux_model_type == "var") - { - /* If the auxiliary model is a VAR, add a parameter corresponding - to the constant. */ - int new_param_symb_id = create_aux_param("h1_" + name + "_" + standard_eqtag + "_constant"); - pac_h1_indices[{name, standard_eqtag}].push_back(new_param_symb_id); - subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id)); - } - for (int i = 1; i < max_lag + 1; i++) - for (auto lhsit : lhs) - { - int new_param_symb_id = create_aux_param("h1_" + name + "_" + standard_eqtag - + "_var_" + symbol_table.getName(lhsit) - + "_lag_" + to_string(i)); - pac_h1_indices[{name, standard_eqtag}].push_back(new_param_symb_id); - subExpr = AddPlus(subExpr, - AddTimes(AddVariable(new_param_symb_id), - AddVariable(lhsit, -i))); - } - } - - if (growth) - { - int growth_param_index = pac_growth_neutrality_params.at(name); - subExpr = AddPlus(subExpr, - AddTimes(AddVariable(growth_param_index), growth)); - } - - pac_expectation_substitution[{name, eqtag}] = subExpr; + for (int i = 1; i < max_lag + 1; i++) + for (auto lhsit : lhs) + { + int new_param_symb_id = create_aux_param("h0_" + name + "_var_" + + symbol_table.getName(lhsit) + + "_lag_" + to_string(i)); + pac_h0_indices[name].push_back(new_param_symb_id); + subExpr = AddPlus(subExpr, + AddTimes(AddVariable(new_param_symb_id), + AddVariable(lhsit, -i))); + } } + + if (nonstationary_vars_present) + { + if (aux_model_type == "var") + { + /* If the auxiliary model is a VAR, add a parameter corresponding + to the constant. */ + int new_param_symb_id = create_aux_param("h1_" + name + "_constant"); + pac_h1_indices[name].push_back(new_param_symb_id); + subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id)); + } + for (int i = 1; i < max_lag + 1; i++) + for (auto lhsit : lhs) + { + int new_param_symb_id = create_aux_param("h1_" + name + "_var_" + + symbol_table.getName(lhsit) + + "_lag_" + to_string(i)); + pac_h1_indices[name].push_back(new_param_symb_id); + subExpr = AddPlus(subExpr, + AddTimes(AddVariable(new_param_symb_id), + AddVariable(lhsit, -i))); + } + } + + pac_expectation_substitution[name] = AddPlus(subExpr, growth_correction_term); } void -DynamicModel::substitutePacExpectation(const map, expr_t> &pac_expectation_substitution) +DynamicModel::substitutePacExpectation(const map &pac_expectation_substitution, + const map &pac_eq_name) { - for (auto &[key, substexpr] : pac_expectation_substitution) + for (auto &[model_name, substexpr] : pac_expectation_substitution) { - auto &[model_name, eq_name] = key; - int eq = equation_tags.getEqnByTag("name", eq_name); + int eq = equation_tags.getEqnByTag("name", pac_eq_name.at(model_name)); auto substeq = dynamic_cast(equations[eq]->substitutePacExpectation(model_name, substexpr)); assert(substeq); equations[eq] = substeq; diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index d493837e..387f0e8b 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -513,7 +513,7 @@ public: //! Substitute VarExpectation operators void substituteVarExpectation(const map &subst_table); - void analyzePacEquationStructure(const string &name, map, pair> &eqtag_and_lag, PacModelTable::equation_info_t &pac_equation_info); + void analyzePacEquationStructure(const string &name, map &pac_eq_name, PacModelTable::equation_info_t &pac_equation_info); // Exception thrown by getPacTargetSymbId() class PacTargetNotIdentifiedException @@ -532,19 +532,17 @@ public: The symbol IDs of the new endogenous and parameters are also added to pac_mce_{z1,alpha}_symb_ids. */ void addPacModelConsistentExpectationEquation(const string &name, int discount, - const map, pair> &eqtag_and_lag, + int pac_eq_max_lag, ExprNode::subst_table_t &diff_subst_table, - map, int> &pac_mce_z1_symb_ids, - map, vector> &pac_mce_alpha_symb_ids); + map &pac_mce_z1_symb_ids, + map> &pac_mce_alpha_symb_ids); /* For a PAC MCE model, fill pac_expectation_substitution with the expression that will be substituted for the pac_expectation operator */ void computePacModelConsistentExpectationSubstitution(const string &name, - const map, pair> &eqtag_and_lag, - expr_t growth, - const map &pac_growth_neutrality_params, - const map, int> &pac_mce_z1_symb_ids, - map, expr_t> &pac_expectation_substitution); + expr_t growth_correction_term, + int pac_mce_z1_symb_id, + map &pac_expectation_substitution); /* For a PAC backward model, fill pac_expectation_substitution with the @@ -554,16 +552,15 @@ public: const vector &lhs, int max_lag, const string &aux_model_type, - const map, pair> &eqtag_and_lag, const vector &nonstationary, - expr_t growth, - const map &pac_growth_neutrality_params, - map, vector> &pac_h0_indices, - map, vector> &pac_h1_indices, - map, expr_t> &pac_expectation_substitution); + expr_t growth_correction_term, + map> &pac_h0_indices, + map> &pac_h1_indices, + map &pac_expectation_substitution); //! Substitutes pac_expectation operator with expectation based on auxiliary model - void substitutePacExpectation(const map, expr_t> &pac_expectation_substitution); + void substitutePacExpectation(const map &pac_expectation_substitution, + const map &pac_eq_name); //! Table to undiff LHS variables for pac vector z vector getUndiffLHSForPac(const string &aux_model_name, diff --git a/src/SubModel.cc b/src/SubModel.cc index f6f9ee7b..3c911437 100644 --- a/src/SubModel.cc +++ b/src/SubModel.cc @@ -815,8 +815,8 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table, DynamicModel &dynamic_model, const VarModelTable &var_model_table, const TrendComponentModelTable &trend_component_model_table) { - // (model name, eq_name) → expression for pac_expectation - map, expr_t> pac_expectation_substitution; + // model name → expression for pac_expectation + map pac_expectation_substitution; for (const auto &name : names) { @@ -866,7 +866,7 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table, cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl; exit(EXIT_FAILURE); } - dynamic_model.analyzePacEquationStructure(name, eqtag_and_lag, equation_info); + dynamic_model.analyzePacEquationStructure(name, eq_name, equation_info); // Declare parameter for growth neutrality correction if (growth[name]) @@ -887,28 +887,30 @@ PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table, // In the MCE case, add the variable and the equation defining Z₁ if (aux_model_name[name].empty()) dynamic_model.addPacModelConsistentExpectationEquation(name, symbol_table.getID(discount[name]), - eqtag_and_lag, diff_subst_table, + pacEquationMaxLag(name), + diff_subst_table, mce_z1_symb_ids, mce_alpha_symb_ids); // Compute the expressions that will be substituted for the pac_expectation operators + expr_t growth_correction_term = dynamic_model.Zero; + if (growth[name]) + growth_correction_term = dynamic_model.AddTimes(growth[name], dynamic_model.AddVariable(growth_neutrality_params[name])); if (aux_model_name[name].empty()) - dynamic_model.computePacModelConsistentExpectationSubstitution(name, eqtag_and_lag, - growth[name], - growth_neutrality_params, - mce_z1_symb_ids, + dynamic_model.computePacModelConsistentExpectationSubstitution(name, + growth_correction_term, + mce_z1_symb_ids[name], pac_expectation_substitution); else dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag, aux_model_type[name], - eqtag_and_lag, nonstationary, - growth[name], - growth_neutrality_params, + nonstationary, + growth_correction_term, h0_indices, h1_indices, pac_expectation_substitution); } // Actually perform the substitution of pac_expectation - dynamic_model.substitutePacExpectation(pac_expectation_substitution); + dynamic_model.substitutePacExpectation(pac_expectation_substitution, eq_name); dynamic_model.checkNoRemainingPacExpectation(); } @@ -976,35 +978,22 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const } // Write PAC Model Consistent Expectation parameter info - for (auto &[key, ids] : mce_alpha_symb_ids) + for (auto &[name, ids] : mce_alpha_symb_ids) { - output << "M_.pac." << key.first << ".equations." << key.second << ".mce.alpha = ["; + output << "M_.pac." << name << ".mce.alpha = ["; for (auto id : ids) output << symbol_table.getTypeSpecificID(id) + 1 << " "; output << "];" << endl; } // Write PAC Model Consistent Expectation Z1 info - for (auto &[key, id] : mce_z1_symb_ids) - output << "M_.pac." << key.first << ".equations." << key.second << ".mce.z1 = " + for (auto &[name, id] : mce_z1_symb_ids) + output << "M_.pac." << name << ".mce.z1 = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl; - // Write PAC lag info - for (auto &[key, val] : eqtag_and_lag) - output << "M_.pac." << key.first << ".equations." << val.first << ".max_lag = " << val.second << ";" << endl; - - // Write PAC equation tag info - map>> for_writing; - for (auto &[key, val] : eqtag_and_lag) - for_writing[key.first].emplace_back(key.second, val.first); - - for (auto &[key, val] : for_writing) - { - output << "M_.pac." << key << ".tag_map = ["; - for (auto &[eqtag, standard_eqtag] : val) - output << "{'" << eqtag << "', '" << standard_eqtag << "'};"; - output << "];" << endl; - } + // Write PAC equation name info + for (auto &[name, eq] : eq_name) + output << "M_.pac." << name << ".eq_name = '" << eq << "';" << endl; for (auto &[model, growth_neutrality_param_index] : growth_neutrality_params) output << "M_.pac." << model << ".growth_neutrality_param_index = " @@ -1021,33 +1010,31 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const for (auto &[model, type] : aux_model_type) output << "M_.pac." << model << ".auxiliary_model_type = '" << type << "';" << endl; - for (auto &[key, val] : equation_info) + for (auto &[name, val] : equation_info) { auto [lhs_pac_var, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants] = val; - string substruct = key.first + ".equations." + key.second + "."; - - output << "M_.pac." << substruct << "lhs_var = " + output << "M_.pac." << name << ".lhs_var = " << symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl; if (optim_share_index >= 0) - output << "M_.pac." << substruct << "share_of_optimizing_agents_index = " + output << "M_.pac." << name << ".share_of_optimizing_agents_index = " << symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl; - output << "M_.pac." << substruct << "ec.params = " + output << "M_.pac." << name << ".ec.params = " << symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl - << "M_.pac." << substruct << "ec.vars = ["; + << "M_.pac." << name << ".ec.vars = ["; for (auto it : ec_params_and_vars.second) output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " "; output << "];" << endl - << "M_.pac." << substruct << "ec.istarget = ["; + << "M_.pac." << name << ".ec.istarget = ["; for (auto it : ec_params_and_vars.second) output << (get<1>(it) ? "true " : "false "); output << "];" << endl - << "M_.pac." << substruct << "ec.scale = ["; + << "M_.pac." << name << ".ec.scale = ["; for (auto it : ec_params_and_vars.second) output << get<2>(it) << " "; output << "];" << endl - << "M_.pac." << substruct << "ec.isendo = ["; + << "M_.pac." << name << ".ec.isendo = ["; for (auto it : ec_params_and_vars.second) switch (symbol_table.getType(get<0>(it))) { @@ -1062,32 +1049,33 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const exit(EXIT_FAILURE); } output << "];" << endl - << "M_.pac." << substruct << "ar.params = ["; + << "M_.pac." << name << ".ar.params = ["; for (auto &[pid, vid, vlag] : ar_params_and_vars) output << (pid != -1 ? symbol_table.getTypeSpecificID(pid) + 1 : -1) << " "; output << "];" << endl - << "M_.pac." << substruct << "ar.vars = ["; + << "M_.pac." << name << ".ar.vars = ["; for (auto &[pid, vid, vlag] : ar_params_and_vars) output << (vid != -1 ? symbol_table.getTypeSpecificID(vid) + 1 : -1) << " "; output << "];" << endl - << "M_.pac." << substruct << "ar.lags = ["; + << "M_.pac." << name << ".ar.lags = ["; for (auto &[pid, vid, vlag] : ar_params_and_vars) output << vlag << " "; - output << "];" << endl; + output << "];" << endl + << "M_.pac." << name << ".max_lag = " << pacEquationMaxLag(name) << ";" << endl; if (!non_optim_vars_params_and_constants.empty()) { - output << "M_.pac." << substruct << "non_optimizing_behaviour.params = ["; + output << "M_.pac." << name << ".non_optimizing_behaviour.params = ["; for (auto &it : non_optim_vars_params_and_constants) if (get<2>(it) >= 0) output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " "; else output << "NaN "; output << "];" << endl - << "M_.pac." << substruct << "non_optimizing_behaviour.vars = ["; + << "M_.pac." << name << ".non_optimizing_behaviour.vars = ["; for (auto &it : non_optim_vars_params_and_constants) output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " "; output << "];" << endl - << "M_.pac." << substruct << "non_optimizing_behaviour.isendo = ["; + << "M_.pac." << name << ".non_optimizing_behaviour.isendo = ["; for (auto &it : non_optim_vars_params_and_constants) switch (symbol_table.getType(get<0>(it))) { @@ -1102,29 +1090,29 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const exit(EXIT_FAILURE); } output << "];" << endl - << "M_.pac." << substruct << "non_optimizing_behaviour.lags = ["; + << "M_.pac." << name << ".non_optimizing_behaviour.lags = ["; for (auto &it : non_optim_vars_params_and_constants) output << get<1>(it) << " "; output << "];" << endl - << "M_.pac." << substruct << "non_optimizing_behaviour.scaling_factor = ["; + << "M_.pac." << name << ".non_optimizing_behaviour.scaling_factor = ["; for (auto &it : non_optim_vars_params_and_constants) output << get<3>(it) << " "; output << "];" << endl; } if (!additive_vars_params_and_constants.empty()) { - output << "M_.pac." << substruct << "additive.params = ["; + output << "M_.pac." << name << ".additive.params = ["; for (auto &it : additive_vars_params_and_constants) if (get<2>(it) >= 0) output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " "; else output << "NaN "; output << "];" << endl - << "M_.pac." << substruct << "additive.vars = ["; + << "M_.pac." << name << ".additive.vars = ["; for (auto &it : additive_vars_params_and_constants) output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " "; output << "];" << endl - << "M_.pac." << substruct << "additive.isendo = ["; + << "M_.pac." << name << ".additive.isendo = ["; for (auto &it : additive_vars_params_and_constants) switch (symbol_table.getType(get<0>(it))) { @@ -1139,29 +1127,29 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const exit(EXIT_FAILURE); } output << "];" << endl - << "M_.pac." << substruct << "additive.lags = ["; + << "M_.pac." << name << ".additive.lags = ["; for (auto &it : additive_vars_params_and_constants) output << get<1>(it) << " "; output << "];" << endl - << "M_.pac." << substruct << "additive.scaling_factor = ["; + << "M_.pac." << name << ".additive.scaling_factor = ["; for (auto &it : additive_vars_params_and_constants) output << get<3>(it) << " "; output << "];" << endl; } if (!optim_additive_vars_params_and_constants.empty()) { - output << "M_.pac." << substruct << "optim_additive.params = ["; + output << "M_.pac." << name << ".optim_additive.params = ["; for (auto &it : optim_additive_vars_params_and_constants) if (get<2>(it) >= 0) output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " "; else output << "NaN "; output << "];" << endl - << "M_.pac." << substruct << "optim_additive.vars = ["; + << "M_.pac." << name << ".optim_additive.vars = ["; for (auto &it : optim_additive_vars_params_and_constants) output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " "; output << "];" << endl - << "M_.pac." << substruct << "optim_additive.isendo = ["; + << "M_.pac." << name << ".optim_additive.isendo = ["; for (auto &it : optim_additive_vars_params_and_constants) switch (symbol_table.getType(get<0>(it))) { @@ -1176,31 +1164,31 @@ PacModelTable::writeOutput(const string &basename, ostream &output) const exit(EXIT_FAILURE); } output << "];" << endl - << "M_.pac." << substruct << "optim_additive.lags = ["; + << "M_.pac." << name << ".optim_additive.lags = ["; for (auto &it : optim_additive_vars_params_and_constants) output << get<1>(it) << " "; output << "];" << endl - << "M_.pac." << substruct << "optim_additive.scaling_factor = ["; + << "M_.pac." << name << ".optim_additive.scaling_factor = ["; for (auto &it : optim_additive_vars_params_and_constants) output << get<3>(it) << " "; output << "];" << endl; } // Create empty h0 and h1 substructures that will be overwritten later if not empty - output << "M_.pac." << substruct << "h0_param_indices = [];" << endl - << "M_.pac." << substruct << "h1_param_indices = [];" << endl; + output << "M_.pac." << name << ".h0_param_indices = [];" << endl + << "M_.pac." << name << ".h1_param_indices = [];" << endl; } - for (auto &[model_eqtag, symb_ids] : h0_indices) + for (auto &[name, symb_ids] : h0_indices) { - output << "M_.pac." << model_eqtag.first << ".equations." << model_eqtag.second << ".h0_param_indices = ["; + output << "M_.pac." << name << ".h0_param_indices = ["; for (auto it : symb_ids) output << symbol_table.getTypeSpecificID(it) + 1 << " "; output << "];" << endl; } - for (auto &[model_eqtag, symb_ids] : h1_indices) + for (auto &[name, symb_ids] : h1_indices) { - output << "M_.pac." << model_eqtag.first << ".equations." << model_eqtag.second << ".h1_param_indices = ["; + output << "M_.pac." << name << ".h1_param_indices = ["; for (auto it : symb_ids) output << symbol_table.getTypeSpecificID(it) + 1 << " "; output << "];" << endl; @@ -1228,3 +1216,8 @@ PacModelTable::writeJsonOutput(ostream &output) const } } +int +PacModelTable::pacEquationMaxLag(const string &name_arg) const +{ + return get<2>(equation_info.at(name_arg)).size(); +} diff --git a/src/SubModel.hh b/src/SubModel.hh index a4fe1906..00997796 100644 --- a/src/SubModel.hh +++ b/src/SubModel.hh @@ -201,17 +201,21 @@ private: map growth, original_growth; map>> growth_info; + /* Stores the name of the PAC equation associated to the model. + pac_model_name → eq_name */ + map eq_name; + /* Stores symb_ids for alphas created by DynamicModel::addPacModelConsistentExpectationEquation() - (pac_model_name, standardized_eqtag) -> mce_alpha_symb_ids */ - map, vector> mce_alpha_symb_ids; + pac_model_name → mce_alpha_symb_ids */ + map> mce_alpha_symb_ids; /* Stores symb_ids for z1s created by DynamicModel::addPacModelConsistentExpectationEquation() - (pac_model_name, standardized_eqtag) -> mce_z1_symb_id */ - map, int> mce_z1_symb_ids; + pac_model_name → mce_z1_symb_id */ + map mce_z1_symb_ids; /* Stores symb_ids for h0, h1 parameters - (pac_model_name, standardized_eqtag) -> parameter symb_ids */ - map, vector> h0_indices, h1_indices; + pac_model_name → parameter symb_ids */ + map> h0_indices, h1_indices; /* Stores indices for growth neutrality parameters - pac_model_name -> growth_neutrality_param_index */ + pac_model_name → growth_neutrality_param_index */ map growth_neutrality_params; // Stores LHS vars (only for backward PAC models) @@ -220,19 +224,18 @@ private: // Stores auxiliary model type (only for backward PAC models) map aux_model_type; - /* Stores lag info for pac equations - (pac_model_name, equation_tag) -> (standardized_eqtag, lag) */ - map, pair> eqtag_and_lag; public: - /* Stores info about PAC models specific to the equation they appear in - (pac_model_name, standardized_eqtag) -> + /* Stores info about PAC equations + pac_model_name → (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants) */ - using equation_info_t = map, + using equation_info_t = map, int, vector>, pair>>, vector>, vector>, vector>>>; private: equation_info_t equation_info; + int pacEquationMaxLag(const string &name_arg) const; + public: explicit PacModelTable(SymbolTable &symbol_table_arg); void addPacModel(string name_arg, string aux_model_name_arg, string discount_arg, expr_t growth_arg);