From dfd316b945234b0b81e7a6d74d4a8df4cd682ef2 Mon Sep 17 00:00:00 2001
From: Normann Rion <normann@dynare.org>
Date: Fri, 20 Jun 2025 15:25:49 +0200
Subject: [PATCH] Unfold complementarity conditions in heterogeneous models
using Lagrange multipliers
---
src/HeterogeneousModel.cc | 66 +++++++++++++++++++++++++++++++++++++++
src/HeterogeneousModel.hh | 2 ++
src/ModFile.cc | 4 +++
src/SymbolTable.cc | 38 ++++++++++++++++++++++
src/SymbolTable.hh | 35 ++++++++++++++++++++-
5 files changed, 144 insertions(+), 1 deletion(-)
diff --git a/src/HeterogeneousModel.cc b/src/HeterogeneousModel.cc
index d3cc1967..8a967e36 100644
--- a/src/HeterogeneousModel.cc
+++ b/src/HeterogeneousModel.cc
@@ -97,6 +97,72 @@ HeterogeneousModel::computeDerivIDs()
}
}
+/*
+ * Unfold complementarity conditions: (i) declare the multipliers associated
+ * with each bound constraint μ_l and μ_u ; (ii) add or substract the
+ * multiplier into the associated condition; (iii) add the the complementarity
+ * slackness conditions into the set of equations. For example,
+ * households choose {cₜ, aₜ₊₁} to maximize expected lifetime utility:
+ * max 𝐸ₜ [∑ₛ₌₀^∞ βˢ · u(cₜ₊ₛ)]
+ *
+ * Subject to:
+ * 1. Budget constraint: cₜ + aₜ₊₁ = yₜ + (1 + rₜ) · aₜ
+ * 2. Borrowing constraint: aₜ₊₁ ≥ aₘᵢₙ
+ *
+ * Let u'(cₜ) denote the marginal utility of consumption.
+ * Let μₜ ≥ 0 be the Lagrange multiplier on the borrowing constraint.
+ *
+ * Then, the Euler equation becomes:
+ * u′(cₜ) = β · (1 + rₜ₊₁) · u′(cₜ₊₁) − μₜ
+ *
+ * Together with:
+ * aₜ₊₁ ≥ aₘᵢₙ [primal feasibility]
+ * μₜ ≥ 0 [dual feasibility]
+ * μₜ · (aₜ₊₁ − aₘᵢₙ) = 0 [complementarity slackness]
+ * Note that the primal feasibility and dual feasibility constraints are not
+ * introduced here, but Bhandari et al. (2023) show in Appendix B.1 that they
+ * are redundant.
+ */
+void
+HeterogeneousModel::transformPass()
+{
+ for (int i = 0; i < static_cast<int>(equations.size()); ++i)
+ {
+ if (!complementarity_conditions[i])
+ continue;
+
+ /*
+ * `const auto& [symb_id, lb, ub] = *complementarity_conditions[i];` was not used here because
+ * the call to `addEquation` may eventually lead to a resize of the
+ * `complementarity_conditions` vector, which may invalidate the reference to its element. We
+ * take a copy instead for safety.
+ */
+ auto [symb_id, lb, ub] = *complementarity_conditions[i];
+
+ VariableNode* var = getVariable(symb_id);
+ if (lb)
+ {
+ int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
+ heterogeneity_dimension, i, "MULT_L_" + symbol_table.getName(symb_id));
+ expr_t mu_L = AddVariable(mu_id);
+ auto substeq = AddEqual(AddPlus(equations[i]->arg1, mu_L), equations[i]->arg2);
+ assert(substeq);
+ equations[i] = substeq;
+ addEquation(AddEqual(AddTimes(mu_L, AddMinus(var, lb)), Zero), nullopt);
+ }
+ if (ub)
+ {
+ int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
+ heterogeneity_dimension, i, "MULT_U_" + symbol_table.getName(symb_id));
+ auto mu_U = AddVariable(mu_id);
+ auto substeq = AddEqual(AddMinus(equations[i]->arg1, mu_U), equations[i]->arg2);
+ assert(substeq);
+ equations[i] = substeq;
+ addEquation(AddEqual(AddTimes(mu_U, AddMinus(ub, var)), Zero), nullopt);
+ }
+ }
+}
+
void
HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll)
{
diff --git a/src/HeterogeneousModel.hh b/src/HeterogeneousModel.hh
index ed8986ec..15169de7 100644
--- a/src/HeterogeneousModel.hh
+++ b/src/HeterogeneousModel.hh
@@ -38,6 +38,8 @@ public:
HeterogeneousModel(const HeterogeneousModel& m) = default;
HeterogeneousModel& operator=(const HeterogeneousModel& m);
+ void transformPass();
+
void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);
void writeModelFiles(const string& basename, bool julia) const;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index eb2de0b8..fe88bd1a 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -731,6 +731,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
dynamic_model.reorderAuxiliaryEquations();
+ symbol_table.resizeHetAuxVars();
+ for (auto& hm : heterogeneous_models)
+ hm.transformPass();
+
// Freeze the symbol table
symbol_table.freeze();
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 304b4faf..191d1c3c 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -398,6 +398,7 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
case AuxVarType::pacExpectation:
case AuxVarType::pacTargetNonstationary:
case AuxVarType::aggregationOp:
+ case AuxVarType::heterogeneousMultiplier:
break;
case AuxVarType::endoLag:
case AuxVarType::exoLag:
@@ -496,11 +497,26 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
output << basefield << "endo_nbr = " << het_endo_nbr(het_dim) << ";" << endl;
print_symb_names(basefield + "endo_names", het_endo_ids.at(het_dim));
+ output << basefield << "orig_endo_nbr = " << het_orig_endo_nbr(het_dim) << ";" << endl;
+
output << basefield << "exo_nbr = " << het_exo_nbr(het_dim) << ";" << endl;
print_symb_names(basefield + "exo_names", het_exo_ids.at(het_dim));
output << basefield << "param_nbr = " << het_param_nbr(het_dim) << ";" << endl;
print_symb_names(basefield + "param_names", het_param_ids.at(het_dim));
+
+ const vector<AuxVarInfo>& hav = het_aux_vars[het_dim];
+ if (hav.size() == 0)
+ output << basefield << "aux_vars = [];";
+ else
+ for (int i = 0; i < static_cast<int>(hav.size()); i++)
+ {
+ output << basefield << "aux_vars(" << i + 1
+ << ").endo_index = " << getTypeSpecificID(hav[i].symb_id) + 1 << ";" << endl
+ << basefield << "aux_vars(" << i + 1 << ").type = " << hav[i].get_type_id()
+ << ";" << basefield << "aux_vars(" << i + 1
+ << ").eq_nbr = " << hav[i].equation_number_for_multiplier + 1 << ";" << endl;
+ }
}
}
@@ -713,6 +729,27 @@ SymbolTable::addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op,
return symb_id;
}
+int
+SymbolTable::addHeterogeneousMultiplierAuxiliaryVar(int het_dim, int index,
+ const string& varname) noexcept(false)
+{
+ int symb_id;
+ try
+ {
+ symb_id = addSymbol(varname, SymbolType::heterogeneousEndogenous, "", {}, het_dim);
+ }
+ catch (AlreadyDeclaredException& e)
+ {
+ cerr << "ERROR: you should rename your variable called " << varname
+ << ", this name is internally used by Dynare" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ het_aux_vars[het_dim].emplace_back(symb_id, AuxVarType::heterogeneousMultiplier, 0, 0, index, 0,
+ nullptr, "");
+ return symb_id;
+}
+
int
SymbolTable::addMultiplierAuxiliaryVar(int index) noexcept(false)
{
@@ -1100,6 +1137,7 @@ SymbolTable::writeJsonOutput(ostream& output) const
case AuxVarType::pacExpectation:
case AuxVarType::pacTargetNonstationary:
case AuxVarType::aggregationOp:
+ case AuxVarType::heterogeneousMultiplier:
break;
case AuxVarType::endoLag:
case AuxVarType::exoLag:
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 14eb6ca5..e9b9d98a 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -56,7 +56,9 @@ enum class AuxVarType
pacTargetNonstationary
= 13, //!< Variable created for the substitution of the pac_target_nonstationary operator
aggregationOp
- = 14 // Substitute for an aggregation operator in a heterogeneous setup, such as SUM()
+ = 14, // Substitute for an aggregation operator in a heterogeneous setup, such as SUM()
+ heterogeneousMultiplier = 15 /* Multiplier for bound conditions of complementarity conditions in
+ the heterogeneous equations */
};
//! Information on some auxiliary variables
@@ -155,6 +157,9 @@ private:
//! Information about auxiliary variables
vector<AuxVarInfo> aux_vars;
+ //! Information about heterogeneous auxiliary variables
+ vector<vector<AuxVarInfo>> het_aux_vars;
+
//! Stores the predetermined variables (by symbol IDs)
set<int> predetermined_variables;
@@ -282,6 +287,14 @@ public:
\return the symbol ID of the new symbol
*/
int addExpectationAuxiliaryVar(int information_set, int index, expr_t arg) noexcept(false);
+ /* Adds an auxiliary variable for the bound-specific multiplier of the MCPs and returns its
+ symbolID.
+ – het_dim is the heterogeneity dimension of the model
+ – index is the equation's index
+ – varname is the multiplier name
+ */
+ int addHeterogeneousMultiplierAuxiliaryVar(int het_dim, int index,
+ const string& varname) noexcept(false);
//! Adds an auxiliary variable for the multiplier for the FOCs of the Ramsey Problem
/*!
\param[in] index Used to construct the variable name
@@ -347,6 +360,8 @@ public:
int addPacTargetNonstationaryAuxiliaryVar(const string& name, expr_t expr_arg);
// An auxiliary variable for an aggregation operator (e.g. SUM(yh) where yh is heterogeneous)
int addAggregationOpAuxiliaryVar(const string& name, expr_t expr_arg);
+ //! Set the size of het_aux_vars
+ inline void resizeHetAuxVars();
//! Returns the number of auxiliary variables
[[nodiscard]] int
AuxVarsSize() const
@@ -398,6 +413,9 @@ public:
[[nodiscard]] inline int param_nbr() const noexcept(false);
//! Get number of heterogeneous endogenous variables along a given dimension
[[nodiscard]] inline int het_endo_nbr(int het_dim) const noexcept(false);
+ //! Get number of user-declared heterogeneous endogenous variables (without
+ //! the auxiliary variables)
+ [[nodiscard]] inline int het_orig_endo_nbr(int het_dim) const noexcept(false);
//! Get number of heterogeneous exogenous variables along a given dimension
[[nodiscard]] inline int het_exo_nbr(int het_dim) const noexcept(false);
//! Get number of heterogeneous parameters along a given dimension
@@ -477,6 +495,12 @@ SymbolTable::validateSymbID(int symb_id) const noexcept(false)
throw UnknownSymbolIDException {symb_id};
}
+inline void
+SymbolTable::resizeHetAuxVars()
+{
+ het_aux_vars.resize(heterogeneity_table.size());
+}
+
inline bool
SymbolTable::exists(const string& name) const
{
@@ -591,6 +615,15 @@ SymbolTable::het_endo_nbr(int het_dim) const noexcept(false)
return het_endo_ids.at(het_dim).size();
}
+inline int
+SymbolTable::het_orig_endo_nbr(int het_dim) const noexcept(false)
+{
+ if (!frozen)
+ throw NotYetFrozenException();
+
+ return het_endo_ids.at(het_dim).size() - het_aux_vars.at(het_dim).size();
+}
+
inline int
SymbolTable::het_exo_nbr(int het_dim) const noexcept(false)
{
--
GitLab