From a724cd0c578ffead357af67c10279075c58827f7 Mon Sep 17 00:00:00 2001
From: Normann Rion <normann@dynare.org>
Date: Mon, 8 Apr 2024 11:51:04 +0200
Subject: [PATCH] Handle PAC models with MCE and a composite target
(pac_target_info block)
---
src/DynamicModel.cc | 202 ++++++++++++++++++++++++++++++++++++++++++++
src/DynamicModel.hh | 12 +++
src/SubModel.cc | 52 ++++++++++--
src/SymbolTable.cc | 2 +-
4 files changed, 262 insertions(+), 6 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 4f884b8b..3d9f2df4 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -2084,6 +2084,208 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(
pac_expectation_substitution[name] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction_term);
}
+void
+DynamicModel::computePacModelConsistentExpectationSubstitutionWithComponents(
+ const string& name, int discount_symb_id, int pac_eq_max_lag,
+ ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+ vector<PacModelTable::target_component_t>& pac_target_components,
+ map<string, expr_t>& pac_expectation_substitution)
+{
+ auto create_aux_param = [&](const string& param_name) {
+ try
+ {
+ return symbol_table.addSymbol(param_name, SymbolType::parameter);
+ }
+ catch (SymbolTable::AlreadyDeclaredException)
+ {
+ cerr << "ERROR: the variable/parameter '" << param_name
+ << "' conflicts with some auxiliary parameter that will be generated for the '" << name
+ << "' PAC model. Please rename that parameter." << endl;
+ exit(EXIT_FAILURE);
+ }
+ };
+ int neqs = 0;
+
+ // At the end of this loop:
+ // ₘ ₘ
+ // A_1 = 1+∑ αᵢ = A(1) and A_beta = 1+∑ αᵢβⁱ
+ // ᵢ₌₁ ᵢ₌₁
+ expr_t A_1 = One;
+ expr_t A_beta = One;
+ 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
+ {
+ int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+ pac_aux_param_symb_ids[name].push_back(alpha_i_symb_id);
+ A_1 = AddPlus(A_1, AddVariable(alpha_i_symb_id));
+ A_beta = AddPlus(A_beta, AddTimes(AddVariable(alpha_i_symb_id),
+ AddPower(beta, AddPossiblyNegativeConstant(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);
+ }
+ }
+
+ auto create_target_lag = [&](int variable, int lag) {
+ if (symbol_table.isAuxiliaryVariable(variable))
+ return AddVariable(symbol_table.getOrigSymbIdForAuxVar(variable), lag);
+ else
+ return AddVariable(variable, lag);
+ };
+
+ expr_t substexpr = Zero;
+ for (int component_idx {1};
+ auto& [component, growth, auxname, kind, coeff, growth_neutrality_param, h_indices,
+ original_growth, growth_info] : pac_target_components)
+ {
+ // Create the auxiliary variable for this component
+ int aux_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, component);
+ expr_t aux_var = AddVariable(aux_id);
+
+ // Add the equation defining the auxiliary variable for this component
+ expr_t neweq = AddEqual(aux_var, component);
+ addEquation(neweq, nullopt);
+ addAuxEquation(neweq);
+ neqs++;
+
+ // ₘ
+ // fp = ∑ αᵢβⁱZₜ₊ᵢ
+ // ᵢ₌₁
+ expr_t fp = Zero;
+ int mce_z_symb_id = symbol_table.addPacExpectationAuxiliaryVar("mce_Z_" + auxname, nullptr);
+ for (int i = 1; i <= pac_eq_max_lag + 1; i++)
+ {
+ int alpha_i_symb_id = -1;
+ string param_name = "mce_alpha_" + name + "_" + to_string(i);
+ try
+ {
+ alpha_i_symb_id = symbol_table.getID(param_name);
+ }
+ catch (SymbolTable::UnknownSymbolNameException& e)
+ {
+ alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+ }
+ fp = AddPlus(fp, AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
+ AddPower(beta, AddPossiblyNegativeConstant(i))),
+ AddVariable(mce_z_symb_id, i)));
+ }
+ if (kind != PacTargetKind::ll) // non-stationary component y¹ₜ
+ {
+ // Add diff nodes and eqs for the non-stationary component
+ const VariableNode* aux_target_ns_var_diff_node;
+ expr_t diff_node_to_search = AddDiff(create_target_lag(aux_id, 0));
+ if (auto sit = diff_subst_table.find(diff_node_to_search); sit != diff_subst_table.end())
+ aux_target_ns_var_diff_node = sit->second;
+ else
+ {
+ int symb_id
+ = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
+ aux_target_ns_var_diff_node = AddVariable(symb_id);
+ auto neweq
+ = AddEqual(const_cast<VariableNode*>(aux_target_ns_var_diff_node),
+ AddMinus(create_target_lag(aux_id, 0), create_target_lag(aux_id, -1)));
+ addEquation(neweq, nullopt);
+ addAuxEquation(neweq);
+ neqs++;
+ }
+
+ map<int, VariableNode*> target_ns_aux_var_to_add;
+ const VariableNode* last_aux_var = aux_target_ns_var_diff_node;
+ for (int i = 1; i <= pac_eq_max_lag; i++, neqs++)
+ {
+ expr_t this_diff_node = AddDiff(create_target_lag(aux_id, i));
+ int symb_id = symbol_table.addDiffLeadAuxiliaryVar(
+ this_diff_node->idx, this_diff_node, last_aux_var->symb_id, 1);
+ VariableNode* current_aux_var = AddVariable(symb_id);
+ auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1));
+ addEquation(neweq, nullopt);
+ addAuxEquation(neweq);
+ last_aux_var = current_aux_var;
+ target_ns_aux_var_to_add[i] = current_aux_var;
+ }
+
+ // At the end of this loop,
+ // ₘ₋₁ ⎛ ₘ₋₁ ⎞
+ // fs = ∑ ⎢ ∑ αᵢβⁱ⎥ Δy¹ₜ₊ᵢ
+ // ₖ₌₁ ⎝ ᵢ₌ₖ ⎠
+ 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_ns_aux_var_to_add[k]));
+ }
+
+ // Assembling the equation Zₜ¹ = A_1 ( Δyₜ¹ - fs ) - fp_1
+ auto neweq_1 = AddEqual(
+ AddVariable(mce_z_symb_id),
+ AddMinus(
+ AddTimes(A_1,
+ AddMinus(const_cast<VariableNode*>(aux_target_ns_var_diff_node), fs)),
+ fp));
+ addEquation(neweq_1, nullopt);
+ neqs++;
+ /* This equation is not added to the list of auxiliary equations, because it
+ is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+ }
+ else // Stationary component yₜ⁰
+ {
+ // Assembling the equation Zₜ⁰ = A_beta A_1 yₜ⁰ - fp
+ auto neweq_0 = AddEqual(AddVariable(mce_z_symb_id),
+ AddMinus(AddTimes(A_beta, AddTimes(A_1, component)), fp));
+ addEquation(neweq_0, nullopt);
+ neqs++;
+ /* This equation is not added to the list of auxiliary equations, because it
+ is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+ }
+
+ // If needed, add the growth neutrality correction for this component
+ expr_t growth_correction_term = Zero;
+ string name_component = name + "_component" + to_string(component_idx);
+ if (growth)
+ {
+ growth_neutrality_param
+ = create_aux_param(name_component + "_pac_growth_neutrality_correction");
+ growth_correction_term = AddTimes(growth, AddVariable(growth_neutrality_param));
+ }
+ else
+ growth_neutrality_param = -1;
+
+ substexpr = AddPlus(
+ substexpr, AddTimes(coeff, AddPlus(AddVariable(mce_z_symb_id), growth_correction_term)));
+ component_idx++;
+ }
+
+ cout << "PAC Model Consistent Expectation: added " << neqs
+ << " auxiliary variables and equations for model " << name << "." << endl;
+
+ /* 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] = substexpr;
+}
+
void
DynamicModel::computePacBackwardExpectationSubstitution(
const string& name, const vector<int>& lhs, int max_lag, const string& aux_model_type,
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 05d5b444..262d7599 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -638,6 +638,18 @@ public:
map<string, int>& pac_aux_var_symb_ids, map<string, vector<int>>& pac_aux_param_symb_ids,
map<string, expr_t>& pac_expectation_substitution);
+ /* For a PAC MCE model with an associated pac_target_info, fill pac_expectation_substitution with
+ the expression that will be substituted for the pac_expectation operator. In the process, add
+ the variables and the equations defining Z₁ and Z₀ for each component. The new auxiliary
+ parameters are added to pac_mce_alpha_symb_ids. The routine also creates the auxiliary
+ variables for the components, and adds the corresponding equations.
+ */
+ void computePacModelConsistentExpectationSubstitutionWithComponents(
+ const string& name, int discount_symb_id, int pac_eq_max_lag,
+ ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+ vector<PacModelTable::target_component_t>& pac_target_components,
+ map<string, expr_t>& pac_expectation_substitution);
+
/* For a PAC backward model, fill pac_expectation_substitution with the
expression that will be substituted for the pac_expectation operator.
The symbol IDs of the new parameters are also added to pac_aux_param_symb_ids.
diff --git a/src/SubModel.cc b/src/SubModel.cc
index 38b729ba..6f504424 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -1384,11 +1384,11 @@ PacModelTable::transformPass(const lag_equivalence_table_t& unary_ops_nodes,
{
if (target_info.contains(name))
{
- cerr << "ERROR: the block 'pac_target_info(" << name
- << ")' is not supported in the context of a PAC model with model-consistent "
- "expectations (MCE)."
- << endl;
- exit(EXIT_FAILURE);
+ assert(growth_correction_term == dynamic_model.Zero);
+ dynamic_model.computePacModelConsistentExpectationSubstitutionWithComponents(
+ name, symbol_table.getID(discount[name]), pacEquationMaxLag(name),
+ diff_subst_table, aux_param_symb_ids, get<2>(target_info[name]),
+ pac_expectation_substitution);
}
else
dynamic_model.computePacModelConsistentExpectationSubstitution(
@@ -1709,6 +1709,48 @@ PacModelTable::writeOutput(ostream& output) const
}
component_idx++;
}
+
+ for (auto& [name, val] : target_info)
+ if (aux_model_name.at(name).empty())
+ {
+ string pac_model_name = "M_.pac." + name + ".mce.";
+ output << pac_model_name << "z = NaN(" << get<2>(val).size() << ",1);" << endl
+ << pac_model_name << "w = NaN(" << get<2>(val).size() << ",1);" << endl
+ << pac_model_name << "b = NaN(" << get<2>(val).size() << ",1);" << endl;
+ for (int component_idx {1};
+ auto& [component, growth_component, auxname, kind, coeff, growth_neutrality_param,
+ h_indices, original_growth_component, growth_component_info] : get<2>(val))
+ {
+ output << pac_model_name << "z(" << component_idx
+ << ") = " << symbol_table.getTypeSpecificID("mce_Z_" + auxname) + 1 << ";"
+ << endl;
+ set<int> params;
+ coeff->collectVariables(SymbolType::parameter, params);
+ /* If the coefficient expression uses numerical constants only, set b to 0 and w to the
+ computed numerical weight
+ */
+ if (params.empty())
+ {
+ ostringstream mce_w;
+ coeff->writeOutput(mce_w);
+ output << pac_model_name << "b(" << component_idx << ") = 0;" << endl
+ << pac_model_name << "w(" << component_idx << ") = " << mce_w.str() << ";"
+ << endl;
+ }
+ /* If the coefficient expression is simply a parameter, set b to 1 and w to
+ the parameter index
+ */
+ if ((params.size() == 1) && (typeid(*coeff) == typeid(VariableNode)))
+ output << pac_model_name << "b(" << component_idx << ") = 1;" << endl
+ << pac_model_name << "w(" << component_idx
+ << ") = " << symbol_table.getTypeSpecificID(*params.begin()) + 1 << ";"
+ << endl;
+ /* If the coefficient expression mixes numerical constants and parameters, the
+ associated values in w and b are NaN
+ */
+ component_idx++;
+ }
+ }
}
void
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index d836566d..83e46953 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -731,7 +731,7 @@ SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id_arg) const noexcept(fals
if ((aux_var.type == AuxVarType::endoLag || aux_var.type == AuxVarType::exoLag
|| aux_var.type == AuxVarType::diff || aux_var.type == AuxVarType::diffLag
|| aux_var.type == AuxVarType::diffLead || aux_var.type == AuxVarType::diffForward
- || aux_var.type == AuxVarType::unaryOp)
+ || aux_var.type == AuxVarType::unaryOp || aux_var.type == AuxVarType::pacExpectation)
&& aux_var.symb_id == aux_var_symb_id_arg)
{
if (optional<int> r = aux_var.orig_symb_id; r)
--
GitLab