From 1e68b2fcda6834919b61660e0ba7fb4d20f27dc9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 23 May 2024 20:45:36 +0200
Subject: [PATCH] Initial HANK support
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
– New “heterogeneity_dimension†statement
– New option “heterogeneity†to “varâ€, “varexoâ€, “parametersâ€, “model†and
“shocks†statements
– New “SUM()†operator in “model†block
---
src/CommonEnums.hh | 31 ++++--
src/DataTree.cc | 12 +-
src/DataTree.hh | 9 +-
src/DynamicModel.cc | 36 +++++-
src/DynamicModel.hh | 5 +
src/DynareBison.yy | 53 +++++++--
src/DynareFlex.ll | 3 +
src/ExprNode.cc | 210 +++++++++++++++++++++++++++++++++-
src/ExprNode.hh | 50 ++++++---
src/HeterogeneityTable.cc | 139 +++++++++++++++++++++++
src/HeterogeneityTable.hh | 100 +++++++++++++++++
src/HeterogeneousModel.cc | 229 ++++++++++++++++++++++++++++++++++++++
src/HeterogeneousModel.hh | 101 +++++++++++++++++
src/ModFile.cc | 207 +++++++++++++++++++++++++++++++---
src/ModFile.hh | 6 +
src/ModelEquationBlock.cc | 30 +++--
src/ModelEquationBlock.hh | 6 +-
src/ModelTree.cc | 38 ++++++-
src/ModelTree.hh | 143 ++++++++++++------------
src/ParsingDriver.cc | 165 +++++++++++++++++++++++----
src/ParsingDriver.hh | 29 +++--
src/Shocks.cc | 209 +++++++++++++++++++++++++++++++++-
src/Shocks.hh | 44 +++++++-
src/StaticModel.cc | 10 +-
src/StaticModel.hh | 3 +-
src/SymbolList.cc | 9 ++
src/SymbolTable.cc | 152 +++++++++++++++++++++++--
src/SymbolTable.hh | 84 +++++++++++++-
src/meson.build | 2 +
29 files changed, 1921 insertions(+), 194 deletions(-)
create mode 100644 src/HeterogeneityTable.cc
create mode 100644 src/HeterogeneityTable.hh
create mode 100644 src/HeterogeneousModel.cc
create mode 100644 src/HeterogeneousModel.hh
diff --git a/src/CommonEnums.hh b/src/CommonEnums.hh
index 2b20225b..ceadbafb 100644
--- a/src/CommonEnums.hh
+++ b/src/CommonEnums.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -25,14 +25,17 @@
* command */
enum class SymbolType
{
- endogenous = 0, //!< Endogenous
- exogenous = 1, //!< Exogenous
- exogenousDet = 2, //!< Exogenous deterministic
- parameter = 4, //!< Parameter
- modelLocalVariable = 10, //!< Local variable whose scope is model (pound expression)
- modFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
- externalFunction = 12, //!< External (user-defined) function
- trend = 13, //!< Trend variable
+ endogenous = 0, // Endogenous (non-heterogeneous)
+ exogenous = 1, // Exogenous (non-heterogeneous)
+ exogenousDet = 2, // Exogenous deterministic (non-heterogeneous)
+ parameter = 4, // Parameter (non-heterogeneous)
+ heterogeneousEndogenous = 5, // Endogenous that is heterogeneous across some dimension
+ heterogeneousExogenous = 6, // Exogenous that is heterogeneous across some dimension
+ heterogeneousParameter = 7, // Parameter that is heterogeneous across some dimension
+ modelLocalVariable = 10, // Local variable whose scope is model (pound expression)
+ modFileLocalVariable = 11, // Local variable whose scope is mod file (model excluded)
+ externalFunction = 12, // External (user-defined) function
+ trend = 13, // Trend variable
statementDeclaredVariable
= 14, //!< Local variable assigned within a Statement (see subsample statement for example)
logTrend = 15, //!< Log-trend variable
@@ -45,6 +48,13 @@ enum class SymbolType
excludedVariable = 19 //!< Variable excluded via model_remove/var_remove/include_eqs/exclude_eqs
};
+constexpr bool
+isHeterogeneous(SymbolType type)
+{
+ return type == SymbolType::heterogeneousEndogenous || type == SymbolType::heterogeneousExogenous
+ || type == SymbolType::heterogeneousParameter;
+}
+
enum class UnaryOpcode
{
uminus,
@@ -75,7 +85,8 @@ enum class UnaryOpcode
erf,
erfc,
diff,
- adl
+ adl,
+ sum
};
enum class BinaryOpcode
diff --git a/src/DataTree.cc b/src/DataTree.cc
index 5f77077a..acfd28aa 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -47,10 +47,12 @@ DataTree::initConstants()
}
DataTree::DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
symbol_table {symbol_table_arg},
num_constants {num_constants_arg},
external_functions_table {external_functions_table_arg},
+ heterogeneity_table {heterogeneity_table_arg},
is_dynamic {is_dynamic_arg}
{
initConstants();
@@ -60,6 +62,7 @@ DataTree::DataTree(const DataTree& d) :
symbol_table {d.symbol_table},
num_constants {d.num_constants},
external_functions_table {d.external_functions_table},
+ heterogeneity_table {d.heterogeneity_table},
is_dynamic {d.is_dynamic},
local_variables_vector {d.local_variables_vector}
{
@@ -81,6 +84,7 @@ DataTree::operator=(const DataTree& d)
assert(&symbol_table == &d.symbol_table);
assert(&num_constants == &d.num_constants);
assert(&external_functions_table == &d.external_functions_table);
+ assert(&heterogeneity_table == &d.heterogeneity_table);
assert(is_dynamic == d.is_dynamic);
num_const_node_map.clear();
@@ -791,6 +795,12 @@ DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<exp
return p;
}
+expr_t
+DataTree::AddSum(expr_t arg)
+{
+ return AddUnaryOp(UnaryOpcode::sum, arg);
+}
+
bool
DataTree::isSymbolUsed(int symb_id) const
{
diff --git a/src/DataTree.hh b/src/DataTree.hh
index 4cdd06fe..a0fe9185 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -33,6 +33,7 @@
#include "ExprNode.hh"
#include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
#include "NumericalConstants.hh"
#include "SubModel.hh"
#include "SymbolTable.hh"
@@ -48,6 +49,8 @@ public:
NumericalConstants& num_constants;
//! A reference to the external functions table
ExternalFunctionsTable& external_functions_table;
+ // A reference to the heterogeneity table
+ HeterogeneityTable& heterogeneity_table;
//! Is it possible to use leads/lags on variable nodes?
/* NB: This data member cannot be replaced by a virtual method, because this information is needed
in AddVariable(), which itself can be called from the copy constructor. */
@@ -139,7 +142,8 @@ private:
public:
DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
virtual ~DataTree() = default;
@@ -275,6 +279,9 @@ public:
//! Adds an external function node for the second derivative of an external function
expr_t AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t>& arguments,
int input_index1, int input_index2);
+ // Adds "SUM(arg)" to model tree
+ expr_t AddSum(expr_t arg);
+
//! Checks if a given symbol is used somewhere in the data tree
[[nodiscard]] bool isSymbolUsed(int symb_id) const;
//! Checks if a given unary op is used somewhere in the data tree
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e06ba481..d0e2e15f 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -50,9 +50,11 @@ DynamicModel::copyHelper(const DynamicModel& m)
DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg) :
- ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
+ ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg, true},
trend_component_model_table {trend_component_model_table_arg},
var_model_table {var_model_table_arg}
{
@@ -1133,7 +1135,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
output << "];" << endl;
- writeDriverSparseIndicesHelper<true, false>(output);
+ writeDriverSparseIndicesHelper("dynamic", output);
// Write LHS of each equation in text form
output << "M_.lhs = {" << endl;
@@ -3697,6 +3699,36 @@ DynamicModel::substituteLogTransform()
}
}
+void
+DynamicModel::substituteAggregationOperators()
+{
+ ExprNode::subst_table_t subst_table;
+ vector<BinaryOpNode*> neweqs;
+
+ for (auto& [symb_id, expr] : local_variables_table)
+ expr = expr->substituteAggregationOperators(subst_table, neweqs);
+
+ for (auto& equation : equations)
+ {
+ equation = dynamic_cast<BinaryOpNode*>(
+ equation->substituteAggregationOperators(subst_table, neweqs));
+ assert(equation);
+ }
+
+ for (auto& equation : static_only_equations)
+ {
+ equation = dynamic_cast<BinaryOpNode*>(
+ equation->substituteAggregationOperators(subst_table, neweqs));
+ assert(equation);
+ }
+
+ for (auto neweq : neweqs)
+ {
+ addEquation(neweq, nullopt);
+ addAuxEquation(neweq);
+ }
+}
+
void
DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
{
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 699e3001..2cbcc2cf 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -325,6 +325,7 @@ protected:
public:
DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg);
@@ -566,6 +567,10 @@ public:
// Performs the transformations associated to variables declared with “var(log)â€
void substituteLogTransform();
+ /* Performs the transformations associated to aggregation operators in heterogeneous models, such
+ as SUM(…) */
+ void substituteAggregationOperators();
+
// Check that no variable was declared with “var(log)†in the given equations
void checkNoWithLogTransform(const set<int>& eqnumbers);
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 83cd7cbf..4ca1954f 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -216,6 +216,7 @@ str_tolower(string s)
%token HOMOTOPY_MAX_COMPLETION_SHARE HOMOTOPY_MIN_STEP_SIZE HOMOTOPY_INITIAL_STEP_SIZE HOMOTOPY_STEP_SIZE_INCREASE_SUCCESS_COUNT
%token HOMOTOPY_LINEARIZATION_FALLBACK HOMOTOPY_MARGINAL_LINEARIZATION_FALLBACK HOMOTOPY_EXCLUDE_VAREXO FROM_INITVAL_TO_ENDVAL
%token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS PERPENDICULAR
+%token HETEROGENEITY HETEROGENEITY_DIMENSION SUM
%token <vector<string>> SYMBOL_VEC
@@ -272,6 +273,7 @@ statement : parameters
| predetermined_variables
| model_local_variable
| change_type
+ | heterogeneity_dimension
| model
| initval
| initval_file
@@ -531,9 +533,9 @@ log_trend_var : LOG_TREND_VAR '(' LOG_GROWTH_FACTOR EQUAL { driver.begin_model()
;
var : VAR symbol_list_with_tex_and_partition ';'
- { driver.var($2, false); }
+ { driver.var($2, {}, false); }
| VAR '(' LOG ')' symbol_list_with_tex_and_partition ';'
- { driver.var($5, true); }
+ { driver.var($5, {}, true); }
| VAR '(' DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
{ driver.end_nonstationary_var(false, $6, $8, false); }
| VAR '(' LOG COMMA DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
@@ -542,6 +544,8 @@ var : VAR symbol_list_with_tex_and_partition ';'
{ driver.end_nonstationary_var(true, $6, $8, false); }
/* The case LOG + LOG_DEFLATOR is omitted, because it does not make much sense
from an economic point of view (amounts to taking the log two times) */
+ | VAR '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+ { driver.var($7, $5, false); }
;
var_remove : VAR_REMOVE symbol_list ';' { driver.var_remove($2); };
@@ -615,7 +619,9 @@ var_expectation_model_option : VARIABLE EQUAL symbol
;
varexo : VAREXO symbol_list_with_tex_and_partition ';'
- { driver.varexo($2); }
+ { driver.varexo($2, {}); }
+ | VAREXO '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+ { driver.varexo($7, $5); }
;
varexo_det : VAREXO_DET symbol_list_with_tex_and_partition ';'
@@ -627,7 +633,9 @@ predetermined_variables : PREDETERMINED_VARIABLES symbol_list ';'
;
parameters : PARAMETERS symbol_list_with_tex_and_partition ';'
- { driver.parameters($2); }
+ { driver.parameters($2, {}); }
+ | PARAMETERS '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+ { driver.parameters($7, $5); }
;
model_local_variable : MODEL_LOCAL_VARIABLE symbol_list_with_tex ';'
@@ -648,6 +656,10 @@ change_type_arg : PARAMETERS
{ $$ = SymbolType::exogenousDet; }
;
+heterogeneity_dimension : HETEROGENEITY_DIMENSION symbol_list ';'
+ { driver.heterogeneity_dimension($2); }
+ ;
+
init_param : symbol EQUAL expression ';' { driver.init_param($1, $3); };
expression : '(' expression ')'
@@ -995,6 +1007,8 @@ model : MODEL ';' { driver.begin_model(); }
equation_list END ';' { driver.end_model(); }
| MODEL '(' model_options_list ')' ';' { driver.begin_model(); }
equation_list END ';' { driver.end_model(); }
+ | MODEL '(' HETEROGENEITY EQUAL symbol ')' { driver.begin_heterogeneous_model($5); }';'
+ equation_list END ';' { driver.end_model(); }
;
equation_list : equation_list equation
@@ -1152,6 +1166,8 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_erfc($3); }
| STEADY_STATE '(' hand_side ')'
{ $$ = driver.add_steady_state($3); }
+ | SUM '(' hand_side ')'
+ { $$ = driver.add_sum($3); }
;
comma_hand_side : hand_side
@@ -1205,6 +1221,12 @@ shocks : SHOCKS ';' shock_list END ';' { driver.end_shocks(false); }
| SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, false); }
| SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER COMMA OVERWRITE ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, true); }
| SHOCKS '(' OVERWRITE COMMA LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($7, true); }
+ | SHOCKS '(' HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+ { driver.end_heterogeneous_shocks($5, false); }
+ | SHOCKS '(' HETEROGENEITY EQUAL symbol COMMA OVERWRITE ')' ';' stoch_shock_list END ';'
+ { driver.end_heterogeneous_shocks($5, true); }
+ | SHOCKS '(' OVERWRITE COMMA HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+ { driver.end_heterogeneous_shocks($7, true); }
;
shock_list : shock_list shock_elem
@@ -1212,16 +1234,23 @@ shock_list : shock_list shock_elem
;
shock_elem : det_shock_elem
- | VAR symbol ';' STDERR expression ';'
- { driver.add_stderr_shock($2, $5); }
- | VAR symbol EQUAL expression ';'
- { driver.add_var_shock($2, $4); }
- | VAR symbol COMMA symbol EQUAL expression ';'
- { driver.add_covar_shock($2, $4, $6); }
- | CORR symbol COMMA symbol EQUAL expression ';'
- { driver.add_correl_shock($2, $4, $6); }
+ | stoch_shock_elem
;
+stoch_shock_elem : VAR symbol ';' STDERR expression ';'
+ { driver.add_stderr_shock($2, $5); }
+ | VAR symbol EQUAL expression ';'
+ { driver.add_var_shock($2, $4); }
+ | VAR symbol COMMA symbol EQUAL expression ';'
+ { driver.add_covar_shock($2, $4, $6); }
+ | CORR symbol COMMA symbol EQUAL expression ';'
+ { driver.add_correl_shock($2, $4, $6); }
+ ;
+
+stoch_shock_list : stoch_shock_list stoch_shock_elem
+ | stoch_shock_elem
+ ;
+
det_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
{ driver.add_det_shock($2, $5, $8, ParsingDriver::DetShockType::standard); }
| VAR symbol ';' PERIODS period_list ';' ADD value_list ';'
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 233cfcc3..f826c8b2 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -111,6 +111,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
<INITIAL>predetermined_variables {BEGIN DYNARE_STATEMENT; return token::PREDETERMINED_VARIABLES;}
<INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;}
<INITIAL>model_local_variable {BEGIN DYNARE_STATEMENT; return token::MODEL_LOCAL_VARIABLE;}
+<INITIAL>heterogeneity_dimension {BEGIN DYNARE_STATEMENT; return token::HETEROGENEITY_DIMENSION;}
<INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
<INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
<INITIAL>set_time {BEGIN DYNARE_STATEMENT; return token::SET_TIME;}
@@ -811,6 +812,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
}
<DYNARE_STATEMENT,DYNARE_BLOCK>static_mfs {return token::STATIC_MFS;}
<DYNARE_STATEMENT,DYNARE_BLOCK>balanced_growth_test_tol {return token::BALANCED_GROWTH_TEST_TOL;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>heterogeneity {return token::HETEROGENEITY;}
<DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
<DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}
<DYNARE_BLOCK>normal_pdf {return token::NORMAL_PDF;}
@@ -1011,6 +1013,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
<DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
<DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
<DYNARE_BLOCK>pac_target_nonstationary {return token::PAC_TARGET_NONSTATIONARY;}
+<DYNARE_BLOCK>sum {return token::SUM;}
<DYNARE_STATEMENT>discount {return token::DISCOUNT;}
<DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
<DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 69533fa1..a890a357 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -870,6 +870,13 @@ NumConstNode::substituteLogTransform([[maybe_unused]] int orig_symb_id,
return const_cast<NumConstNode*>(this);
}
+expr_t
+NumConstNode::substituteAggregationOperators([[maybe_unused]] subst_table_t& subst_table,
+ [[maybe_unused]] vector<BinaryOpNode*>& neweqs) const
+{
+ return const_cast<NumConstNode*>(this);
+}
+
VariableNode::VariableNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg, int lag_arg) :
ExprNode {datatree_arg, idx_arg}, symb_id {symb_id_arg}, lag {lag_arg}
{
@@ -900,6 +907,7 @@ VariableNode::prepareForDerivation()
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
+ case SymbolType::heterogeneousEndogenous:
non_null_derivatives.insert(getDerivID());
break;
case SymbolType::modelLocalVariable:
@@ -907,6 +915,8 @@ VariableNode::prepareForDerivation()
// Non null derivatives are those of the value of the local parameter
non_null_derivatives = datatree.getLocalVariable(symb_id, lag)->non_null_derivatives;
break;
+ case SymbolType::heterogeneousExogenous:
+ case SymbolType::heterogeneousParameter:
case SymbolType::modFileLocalVariable:
case SymbolType::statementDeclaredVariable:
case SymbolType::unusedEndogenous:
@@ -955,6 +965,9 @@ VariableNode::prepareForChainRuleDerivation(
case SymbolType::modFileLocalVariable:
case SymbolType::statementDeclaredVariable:
case SymbolType::unusedEndogenous:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
+ case SymbolType::heterogeneousParameter:
// Those variables are never derived using chain rule
non_null_chain_rule_derivatives.try_emplace(const_cast<VariableNode*>(this));
break;
@@ -990,10 +1003,14 @@ VariableNode::computeDerivative(int deriv_id)
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
+ case SymbolType::heterogeneousEndogenous:
if (deriv_id == getDerivID())
return datatree.One;
else
return datatree.Zero;
+ case SymbolType::heterogeneousExogenous:
+ case SymbolType::heterogeneousParameter:
+ return datatree.Zero;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)->getDerivative(deriv_id);
case SymbolType::modFileLocalVariable:
@@ -1069,8 +1086,23 @@ VariableNode::writeJsonAST(ostream& output) const
case SymbolType::excludedVariable:
cerr << "VariableNode::computeDerivative: Impossible case!" << endl;
exit(EXIT_FAILURE);
+ case SymbolType::heterogeneousEndogenous:
+ output << "heterogeneousEndogenous";
+ break;
+ case SymbolType::heterogeneousExogenous:
+ output << "heterogeneousExogenous";
+ break;
+ case SymbolType::heterogeneousParameter:
+ output << "heterogeneousParameter";
+ break;
}
- output << R"(", "lag" : )" << lag << "}";
+ output << '"';
+ if (isHeterogeneous(get_type()))
+ output << R"(, "heterogeneity_dimension" : ")"
+ << datatree.heterogeneity_table.getName(
+ datatree.symbol_table.getHeterogeneityDimension(symb_id))
+ << '"';
+ output << R"(, "lag" : )" << lag << "}";
}
void
@@ -1411,6 +1443,64 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
case SymbolType::excludedVariable:
cerr << "VariableNode::writeOutput: Impossible case" << endl;
exit(EXIT_FAILURE);
+
+ case SymbolType::heterogeneousEndogenous:
+ switch (int tsid = datatree.symbol_table.getTypeSpecificID(symb_id); output_type)
+ {
+ case ExprNodeOutputType::juliaSparseDynamicModel:
+ case ExprNodeOutputType::matlabSparseDynamicModel:
+ case ExprNodeOutputType::CSparseDynamicModel:
+ assert(lag >= -1 && lag <= 1);
+ i = tsid
+ + (lag + 1)
+ * datatree.symbol_table.het_endo_nbr(
+ datatree.symbol_table.getHeterogeneityDimension(symb_id))
+ + ARRAY_SUBSCRIPT_OFFSET(output_type);
+ output << "yh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
+ << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
+ default:
+ cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousEndogenous "
+ "symbol"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ break;
+ case SymbolType::heterogeneousExogenous:
+ i = datatree.symbol_table.getTypeSpecificID(symb_id) + ARRAY_SUBSCRIPT_OFFSET(output_type);
+ switch (output_type)
+ {
+ case ExprNodeOutputType::juliaSparseDynamicModel:
+ case ExprNodeOutputType::matlabSparseDynamicModel:
+ case ExprNodeOutputType::CSparseDynamicModel:
+ assert(lag == 0);
+ output << "xh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
+ << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
+ default:
+ cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousExogenous "
+ "symbol"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ break;
+ case SymbolType::heterogeneousParameter:
+ i = datatree.symbol_table.getTypeSpecificID(symb_id) + ARRAY_SUBSCRIPT_OFFSET(output_type);
+ switch (output_type)
+ {
+ case ExprNodeOutputType::juliaSparseDynamicModel:
+ case ExprNodeOutputType::matlabSparseDynamicModel:
+ case ExprNodeOutputType::CSparseDynamicModel:
+ output << "paramsh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
+ << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
+ default:
+ cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousParameter "
+ "symbol"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ break;
}
}
@@ -1525,6 +1615,7 @@ VariableNode::computeChainRuleDerivative(
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
+ case SymbolType::heterogeneousEndogenous:
if (deriv_id == getDerivID())
return datatree.One;
// If there is in the equation a recursive variable we could use a chaine rule derivation
@@ -1533,6 +1624,9 @@ VariableNode::computeChainRuleDerivative(
non_null_chain_rule_derivatives, cache);
else
return datatree.Zero;
+ case SymbolType::heterogeneousExogenous:
+ case SymbolType::heterogeneousParameter:
+ return datatree.Zero;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)
@@ -1590,6 +1684,9 @@ VariableNode::computeXrefs(EquationInfo& ei) const
case SymbolType::externalFunction:
case SymbolType::epilogue:
case SymbolType::excludedVariable:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
+ case SymbolType::heterogeneousParameter:
break;
}
}
@@ -1688,6 +1785,8 @@ VariableNode::maxLead() const
case SymbolType::endogenous:
case SymbolType::exogenous:
case SymbolType::exogenousDet:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
return lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)->maxLead();
@@ -1704,6 +1803,8 @@ VariableNode::maxLag() const
case SymbolType::endogenous:
case SymbolType::exogenous:
case SymbolType::exogenousDet:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
return -lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)->maxLag();
@@ -1721,6 +1822,8 @@ VariableNode::maxLagWithDiffsExpanded() const
case SymbolType::exogenous:
case SymbolType::exogenousDet:
case SymbolType::epilogue:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
return -lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)->maxLagWithDiffsExpanded();
@@ -1844,6 +1947,8 @@ VariableNode::decreaseLeadsLags(int n) const
case SymbolType::exogenousDet:
case SymbolType::trend:
case SymbolType::logTrend:
+ case SymbolType::heterogeneousEndogenous:
+ case SymbolType::heterogeneousExogenous:
return datatree.AddVariable(symb_id, lag - n);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id, lag)->decreaseLeadsLags(n);
@@ -2240,6 +2345,17 @@ VariableNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return const_cast<VariableNode*>(this);
}
+expr_t
+VariableNode::substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+{
+ if (get_type() == SymbolType::modelLocalVariable)
+ return datatree.getLocalVariable(symb_id, lag)
+ ->substituteAggregationOperators(subst_table, neweqs);
+
+ return const_cast<VariableNode*>(this);
+}
+
UnaryOpNode::UnaryOpNode(DataTree& datatree_arg, int idx_arg, UnaryOpcode op_code_arg,
const expr_t arg_arg, int expectation_information_set_arg,
int param1_symb_id_arg, int param2_symb_id_arg, string adl_param_name_arg,
@@ -2265,12 +2381,12 @@ UnaryOpNode::prepareForDerivation()
/* Non-null derivatives are those of the argument (except for STEADY_STATE in
a dynamic context, in which case the potentially non-null derivatives are
- all the parameters) */
+ all the parameters; and for SUM, which will be substituted out before deriving) */
if ((op_code == UnaryOpcode::steadyState || op_code == UnaryOpcode::steadyStateParamDeriv
|| op_code == UnaryOpcode::steadyStateParam2ndDeriv)
&& datatree.is_dynamic)
datatree.addAllParamDerivId(non_null_derivatives);
- else
+ else if (op_code != UnaryOpcode::sum)
{
arg->prepareForDerivation();
non_null_derivatives = arg->non_null_derivatives;
@@ -2447,6 +2563,9 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
case UnaryOpcode::adl:
cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
+ case UnaryOpcode::sum:
+ cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::sum" << endl;
+ exit(EXIT_FAILURE);
}
__builtin_unreachable(); // Silence GCC warning
}
@@ -2541,6 +2660,8 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
case UnaryOpcode::adl:
cerr << "UnaryOpNode::cost: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
+ case UnaryOpcode::sum:
+ return 0; // In the generated files, the SUM() operator behaves like a variable
}
else
// Cost for C files
@@ -2591,6 +2712,8 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
case UnaryOpcode::adl:
cerr << "UnaryOpNode::cost: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
+ case UnaryOpcode::sum:
+ return 0; // In the generated files, the SUM() operator behaves like a variable
}
__builtin_unreachable(); // Silence GCC warning
}
@@ -2735,6 +2858,9 @@ UnaryOpNode::writeJsonAST(ostream& output) const
case UnaryOpcode::erfc:
output << "erfc";
break;
+ case UnaryOpcode::sum:
+ output << "sum";
+ break;
}
output << R"(", "arg" : )";
arg->writeJsonAST(output);
@@ -2887,6 +3013,9 @@ UnaryOpNode::writeJsonOutput(ostream& output, const temporary_terms_t& temporary
case UnaryOpcode::erfc:
output << "erfc";
break;
+ case UnaryOpcode::sum:
+ output << "sum";
+ break;
}
bool close_parenthesis = false;
@@ -3117,6 +3246,20 @@ UnaryOpNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
case UnaryOpcode::adl:
output << "adl";
break;
+ case UnaryOpcode::sum:
+ if (isLatexOutput(output_type))
+ output << "sum";
+ else
+ {
+ auto varg {dynamic_cast<VariableNode*>(arg)};
+ assert(varg && varg->lag == 0);
+ int idx {
+ datatree.heterogeneity_table.getSummedHeterogenousEndogenousIndex(varg->symb_id)};
+ output << "yagg" << LEFT_ARRAY_SUBSCRIPT(output_type) << idx + 1
+ << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ return;
+ }
+ break;
}
if (output_type == ExprNodeOutputType::juliaTimeDataFrame && op_code != UnaryOpcode::uminus)
@@ -3246,6 +3389,9 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) noexcept(false)
case UnaryOpcode::adl:
cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
+ case UnaryOpcode::sum:
+ cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::sum" << endl;
+ exit(EXIT_FAILURE);
}
__builtin_unreachable(); // Silence GCC warning
}
@@ -3503,6 +3649,8 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree& alt_datatree) con
return alt_datatree.AddDiff(alt_arg);
case UnaryOpcode::adl:
return alt_datatree.AddAdl(alt_arg, adl_param_name, adl_lags);
+ case UnaryOpcode::sum:
+ return alt_datatree.AddSum(alt_arg);
}
__builtin_unreachable(); // Silence GCC warning
}
@@ -4095,6 +4243,34 @@ UnaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
+expr_t
+UnaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+{
+ if (op_code == UnaryOpcode::sum)
+ {
+ if (auto it = subst_table.find(const_cast<UnaryOpNode*>(this)); it != subst_table.end())
+ return const_cast<VariableNode*>(it->second);
+
+ // Arriving here, we need to create an auxiliary variable for this operator
+ const VariableNode* varg {dynamic_cast<const VariableNode*>(arg)};
+ assert(varg && varg->lag == 0);
+ const string varname {"SUM_" + varg->getName()};
+ int symb_id {datatree.symbol_table.addAggregationOpAuxiliaryVar(
+ varname, const_cast<UnaryOpNode*>(this))};
+ expr_t newAuxE = datatree.AddVariable(symb_id, 0);
+
+ subst_table[this] = dynamic_cast<VariableNode*>(newAuxE);
+ neweqs.push_back(datatree.AddEqual(newAuxE, const_cast<UnaryOpNode*>(this)));
+
+ datatree.heterogeneity_table.addSummedHeterogeneousEndogenous(varg->symb_id);
+
+ return newAuxE;
+ }
+ else
+ return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
+}
+
BinaryOpNode::BinaryOpNode(DataTree& datatree_arg, int idx_arg, const expr_t arg1_arg,
BinaryOpcode op_code_arg, const expr_t arg2_arg,
int powerDerivOrder_arg) :
@@ -6020,6 +6196,13 @@ BinaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
+expr_t
+BinaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+{
+ return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
+}
+
TrinaryOpNode::TrinaryOpNode(DataTree& datatree_arg, int idx_arg, const expr_t arg1_arg,
TrinaryOpcode op_code_arg, const expr_t arg2_arg,
const expr_t arg3_arg) :
@@ -6859,6 +7042,13 @@ TrinaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
+expr_t
+TrinaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+{
+ return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
+}
+
AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree& datatree_arg, int idx_arg,
int symb_id_arg,
vector<expr_t> arguments_arg) :
@@ -7408,6 +7598,13 @@ AbstractExternalFunctionNode::substituteLogTransform(int orig_symb_id, int aux_s
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
+expr_t
+AbstractExternalFunctionNode::substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+{
+ return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
+}
+
expr_t
AbstractExternalFunctionNode::toStatic(DataTree& static_datatree) const
{
@@ -8766,6 +8963,13 @@ SubModelNode::substituteLogTransform([[maybe_unused]] int orig_symb_id,
return const_cast<SubModelNode*>(this);
}
+expr_t
+SubModelNode::substituteAggregationOperators([[maybe_unused]] subst_table_t& subst_table,
+ [[maybe_unused]] vector<BinaryOpNode*>& neweqs) const
+{
+ return const_cast<SubModelNode*>(this);
+}
+
VarExpectationNode::VarExpectationNode(DataTree& datatree_arg, int idx_arg, string model_name_arg) :
SubModelNode {datatree_arg, idx_arg, move(model_name_arg)}
{
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index aa503a70..dfdb04bb 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -572,38 +572,35 @@ public:
{
};
- //! Returns the maximum lead of endogenous in this expression
+ //! Returns the maximum lead of endogenous in this expression (not incl. heterogeneous endo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxEndoLead() const = 0;
- //! Returns the maximum lead of exogenous in this expression
+ //! Returns the maximum lead of exogenous in this expression (not incl. heterogeneous exo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxExoLead() const = 0;
- //! Returns the maximum lag of endogenous in this expression
+ //! Returns the maximum lag of endogenous in this expression (not incl. heterogeneous endo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxEndoLag() const = 0;
- //! Returns the maximum lag of exogenous in this expression
+ //! Returns the maximum lag of exogenous in this expression (not incl. heterogeneous exo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxExoLag() const = 0;
- //! Returns the maximum lead of endo/exo/exodet in this expression
- /*! A negative value means that the expression contains only lagged
- variables. A value of numeric_limits<int>::min() means that there is
- no variable. */
+ /* Returns the maximum lead of endo/exo/exodet in this expression (including heterogeneous
+ endo/exo). A negative value means that the expression contains only lagged variables. A value
+ of numeric_limits<int>::min() means that there is no variable. */
[[nodiscard]] virtual int maxLead() const = 0;
- //! Returns the maximum lag of endo/exo/exodet in this expression
- /*! A negative value means that the expression contains only leaded
- variables. A value of numeric_limits<int>::min() means that there is
- no variable. */
+ /* Returns the maximum lag of endo/exo/exodet in this expression (including heterogeneous
+ endo/exo). A negative value means that the expression contains only leaded variables. A value
+ of numeric_limits<int>::min() means that there is no variable. */
[[nodiscard]] virtual int maxLag() const = 0;
- //! Returns the maximum lag of endo/exo/exodet, as if diffs were expanded
- /*! This function behaves as maxLag(), except that it treats diff()
- differently. For e.g., on diff(diff(x(-1))), maxLag() returns 1 while
- maxLagWithDiffsExpanded() returns 3. */
+ /* Returns the maximum lag of endo/exo/exodet (including heterogeneous endo/exo), as if diffs were
+ expanded. This function behaves as maxLag(), except that it treats diff() differently. For
+ e.g., on diff(diff(x(-1))), maxLag() returns 1 while maxLagWithDiffsExpanded() returns 3. */
[[nodiscard]] virtual int maxLagWithDiffsExpanded() const = 0;
[[nodiscard]] virtual expr_t undiff() const = 0;
@@ -947,6 +944,13 @@ public:
If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
Otherwise, throws a MatchFailureException. */
[[nodiscard]] virtual tuple<int, expr_t, expr_t> matchComplementarityCondition() const;
+
+ /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
+ Also declares those aggregation operators in the HeterogeneityTable, so as to
+ compute their index in the dedicated vector in argument of the dynamic/static files. */
+ [[nodiscard]] virtual expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const
+ = 0;
};
//! Object used to compare two nodes (using their indexes)
@@ -1058,6 +1062,8 @@ public:
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
};
//! Symbol or variable node
@@ -1165,6 +1171,8 @@ public:
vector<int>& powers) const override;
[[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
};
//! Unary operator node
@@ -1312,6 +1320,8 @@ public:
[[nodiscard]] bool isParamTimesEndogExpr() const override;
void decomposeAdditiveTerms(vector<pair<expr_t, int>>& terms, int current_sign) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
};
//! Binary operator node
@@ -1507,6 +1517,8 @@ public:
vector<int>& powers) const override;
[[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
[[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
};
@@ -1651,6 +1663,8 @@ public:
[[nodiscard]] bool containsPacTargetNonstationary(const string& pac_model_name
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
};
@@ -1826,6 +1840,8 @@ public:
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
};
class ExternalFunctionNode : public AbstractExternalFunctionNode
@@ -2028,6 +2044,8 @@ public:
expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
[[nodiscard]] expr_t removeTrendLeadLag(const map<int, expr_t>& trend_symbols_map) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+ [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+ vector<BinaryOpNode*>& neweqs) const override;
protected:
void prepareForDerivation() override;
diff --git a/src/HeterogeneityTable.cc b/src/HeterogeneityTable.cc
new file mode 100644
index 00000000..1edb63cd
--- /dev/null
+++ b/src/HeterogeneityTable.cc
@@ -0,0 +1,139 @@
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cassert>
+#include <utility>
+
+#include "HeterogeneityTable.hh"
+#include "SymbolTable.hh"
+
+void
+HeterogeneityTable::setSymbolTable(SymbolTable* symbol_table_arg)
+{
+ symbol_table = symbol_table_arg;
+}
+
+int
+HeterogeneityTable::addDimension(string name)
+{
+ if (name_to_id.contains(name))
+ throw AlreadyDeclaredDimensionException {move(name)};
+
+ int id {static_cast<int>(id_to_name.size())};
+ name_to_id.emplace(name, id);
+ id_to_name.push_back(move(name));
+ return id;
+}
+
+bool
+HeterogeneityTable::exists(const string& name) const
+{
+ return name_to_id.contains(name);
+}
+
+int
+HeterogeneityTable::getID(const string& name) const
+{
+ if (auto it = name_to_id.find(name); it != name_to_id.end())
+ return it->second;
+ else
+ throw UnknownDimensionNameException {name};
+}
+
+string
+HeterogeneityTable::getName(int id) const
+{
+ if (id < 0 || id >= static_cast<int>(id_to_name.size()))
+ throw UnknownDimensionIDException {id};
+ else
+ return id_to_name[id];
+}
+
+bool
+HeterogeneityTable::empty() const
+{
+ return name_to_id.empty();
+}
+
+vector<string>
+HeterogeneityTable::getDimensions() const
+{
+ return id_to_name;
+}
+
+int
+HeterogeneityTable::size() const
+{
+ return static_cast<int>(name_to_id.size());
+}
+
+void
+HeterogeneityTable::addSummedHeterogeneousEndogenous(int symb_id)
+{
+ assert(symbol_table->getType(symb_id) == SymbolType::heterogeneousEndogenous);
+ if (summed_het_endo_to_index.contains(symb_id))
+ throw AlreadyDeclaredSummedHeterogeneousEndogenousException {symb_id};
+
+ int index {static_cast<int>(index_to_summed_het_endo.size())};
+ summed_het_endo_to_index.emplace(symb_id, index);
+ index_to_summed_het_endo.push_back(symb_id);
+}
+
+int
+HeterogeneityTable::getSummedHeterogenousEndogenousIndex(int symb_id) const
+{
+ if (auto it = summed_het_endo_to_index.find(symb_id); it != summed_het_endo_to_index.end())
+ return it->second;
+ else
+ throw UnknownSummedHeterogeneousEndogenousException {symb_id};
+}
+
+int
+HeterogeneityTable::aggregateEndoSize() const
+{
+ return index_to_summed_het_endo.size();
+}
+
+void
+HeterogeneityTable::writeOutput(ostream& output) const
+{
+ for (size_t id {0}; id < id_to_name.size(); id++)
+ output << "M_.heterogeneity(" << id + 1 << ").dimension_name = '" << id_to_name[id] << "';"
+ << endl;
+
+ output << "M_.heterogeneity_aggregates = {" << endl;
+ for (int symb_id : index_to_summed_het_endo)
+ output << "'sum', " << symbol_table->getHeterogeneityDimension(symb_id) + 1 << ", "
+ << symbol_table->getTypeSpecificID(symb_id) + 1 << ";" << endl;
+ output << "};" << endl;
+}
+
+void
+HeterogeneityTable::writeJsonOutput(ostream& output) const
+{
+ assert(!empty());
+ output << R"("heterogeneity_dimension": [)";
+ for (bool first_written {false}; const auto& dim : id_to_name)
+ {
+ if (exchange(first_written, true))
+ output << ", ";
+ output << '"' << dim << '"';
+ }
+ output << "]" << endl;
+}
diff --git a/src/HeterogeneityTable.hh b/src/HeterogeneityTable.hh
new file mode 100644
index 00000000..2b76f41d
--- /dev/null
+++ b/src/HeterogeneityTable.hh
@@ -0,0 +1,100 @@
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef HETEROGENEITY_TABLE_HH
+#define HETEROGENEITY_TABLE_HH
+
+#include <map>
+#include <ostream>
+#include <string>
+#include <vector>
+
+using namespace std;
+
+class SymbolTable; // Forward declaration, to avoid circularity
+
+/*
+ There is a guarantee that heterogeneity IDs are increasing, i.e. if dimension A is added after
+ dimension B, then the ID of A is greater than the ID of B.
+ Moreover, the IDs form a contiguous interval starting at 0.
+*/
+class HeterogeneityTable
+{
+private:
+ // Maps dimension names to IDs
+ map<string, int> name_to_id;
+ // Maps dimension IDs to names
+ vector<string> id_to_name;
+
+ SymbolTable* symbol_table {nullptr}; // Cannot be a reference, because of circularity
+
+ /* Keeps track of the SUM() operator instances.
+ Maps a symbol ID that appears inside a SUM() operator into an index in
+ M_.heterogeneity_aggregates */
+ map<int, int> summed_het_endo_to_index;
+
+ // Maps an index in M_.heterogeneity_aggregates into a symbol ID
+ vector<int> index_to_summed_het_endo;
+
+public:
+ void setSymbolTable(SymbolTable* symbol_table_arg);
+
+ struct AlreadyDeclaredDimensionException
+ {
+ // Dimension name
+ const string name;
+ };
+ struct UnknownDimensionNameException
+ {
+ // Dimension name
+ const string name;
+ };
+ struct UnknownDimensionIDException
+ {
+ // Dimension ID
+ const int id;
+ };
+
+ // Returns the dimension ID
+ int addDimension(string name);
+ [[nodiscard]] bool exists(const string& name) const;
+ [[nodiscard]] int getID(const string& name) const;
+ [[nodiscard]] string getName(int id) const;
+ [[nodiscard]] bool empty() const;
+ [[nodiscard]] vector<string> getDimensions() const;
+ [[nodiscard]] int size() const;
+
+ struct AlreadyDeclaredSummedHeterogeneousEndogenousException
+ {
+ const int symb_id;
+ };
+ struct UnknownSummedHeterogeneousEndogenousException
+ {
+ const int symb_id;
+ };
+
+ void addSummedHeterogeneousEndogenous(int symb_id);
+ int getSummedHeterogenousEndogenousIndex(int symb_id) const;
+ [[nodiscard]] int aggregateEndoSize() const;
+
+ void writeOutput(ostream& output) const;
+ void writeJsonOutput(ostream& output) const;
+};
+
+#endif
diff --git a/src/HeterogeneousModel.cc b/src/HeterogeneousModel.cc
new file mode 100644
index 00000000..a8edec13
--- /dev/null
+++ b/src/HeterogeneousModel.cc
@@ -0,0 +1,229 @@
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cstdlib>
+#include <iostream>
+
+#include "HeterogeneousModel.hh"
+
+HeterogeneousModel::HeterogeneousModel(SymbolTable& symbol_table_arg,
+ NumericalConstants& num_constants_arg,
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
+ int heterogeneity_dimension_arg) :
+ ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg, true},
+ heterogeneity_dimension {heterogeneity_dimension_arg}
+{
+}
+
+HeterogeneousModel::HeterogeneousModel(const HeterogeneousModel& m) :
+ ModelTree {m},
+ heterogeneity_dimension {m.heterogeneity_dimension},
+ deriv_id_table {m.deriv_id_table},
+ inv_deriv_id_table {m.inv_deriv_id_table}
+{
+}
+
+HeterogeneousModel&
+HeterogeneousModel::operator=(const HeterogeneousModel& m)
+{
+ ModelTree::operator=(m);
+
+ assert(heterogeneity_dimension == m.heterogeneity_dimension);
+
+ deriv_id_table = m.deriv_id_table;
+ inv_deriv_id_table = m.inv_deriv_id_table;
+
+ return *this;
+}
+
+void
+HeterogeneousModel::computeChainRuleJacobian()
+{
+ cerr << "Heterogeneous::computeChainRuleJacobian(): unimplemented" << endl;
+ exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getBlockJacobianEndoCol([[maybe_unused]] int blk, [[maybe_unused]] int var,
+ [[maybe_unused]] int lead_lag) const
+{
+ cerr << "Heterogeneous::getBlockJacobianEndoCol(): unimplemented" << endl;
+ exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getMFS() const
+{
+ cerr << "Heterogeneous::getMFS(): unimplemented" << endl;
+ exit(EXIT_FAILURE);
+}
+
+void
+HeterogeneousModel::computeDerivIDs()
+{
+ set<pair<int, int>> dynvars;
+
+ for (auto& equation : equations)
+ {
+ equation->collectDynamicVariables(SymbolType::heterogeneousEndogenous, dynvars);
+ equation->collectDynamicVariables(SymbolType::heterogeneousExogenous, dynvars);
+ equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
+ equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
+ equation->collectDynamicVariables(SymbolType::parameter, dynvars);
+ }
+
+ for (const auto& [symb_id, lead_lag] : dynvars)
+ {
+ auto type {symbol_table.getType(symb_id)};
+ if (isHeterogeneous(type))
+ assert(symbol_table.getHeterogeneityDimension(symb_id) == heterogeneity_dimension);
+ if (type == SymbolType::heterogeneousEndogenous || type == SymbolType::endogenous)
+ assert(abs(lead_lag) <= 1);
+ if (type == SymbolType::heterogeneousExogenous || type == SymbolType::exogenous)
+ assert(lead_lag == 0);
+ int deriv_id {static_cast<int>(deriv_id_table.size())};
+ deriv_id_table.emplace(pair {symb_id, lead_lag}, deriv_id);
+ inv_deriv_id_table.emplace_back(symb_id, lead_lag);
+ }
+}
+
+void
+HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll)
+{
+ assert(!use_dll); // Not yet implemented
+
+ computeDerivIDs();
+
+ set<int> vars;
+ for (auto& [symb_lag, deriv_id] : deriv_id_table)
+ if (symbol_table.getType(symb_lag.first) != SymbolType::parameter)
+ vars.insert(deriv_id);
+
+ cout << "Computing " << modelClassName() << " derivatives (order " << derivsOrder << ")." << endl;
+
+ computeDerivatives(derivsOrder, vars);
+
+ computeTemporaryTerms(!use_dll, no_tmp_terms);
+
+ if (ranges::any_of(complementarity_conditions, [](const auto& x) { return x.has_value(); }))
+ {
+ // Implementing it requires modifications in ModelTree::computeMCPEquationsReordering()
+ cerr << "ERROR: Complementarity conditions are not yet implemented in "
+ "model(heterogeneity=...) blocks"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+}
+
+void
+HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
+{
+ assert(!julia); // Not yet implemented
+ writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+}
+
+int
+HeterogeneousModel::getJacobianCol(int deriv_id, bool sparse) const
+{
+ assert(sparse);
+
+ SymbolType type {getTypeByDerivID(deriv_id)};
+ int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+ int lag {getLagByDerivID(deriv_id)};
+
+ if (type == SymbolType::heterogeneousEndogenous)
+ return tsid + (lag + 1) * symbol_table.het_endo_nbr(heterogeneity_dimension);
+
+ int shift {3 * symbol_table.het_endo_nbr(heterogeneity_dimension)};
+
+ if (type == SymbolType::heterogeneousExogenous)
+ return shift + tsid;
+
+ shift += symbol_table.het_exo_nbr(heterogeneity_dimension);
+
+ if (type == SymbolType::endogenous)
+ return shift + tsid + (lag + 1) * symbol_table.endo_nbr();
+
+ shift += symbol_table.endo_nbr();
+
+ if (type == SymbolType::exogenous)
+ return shift + tsid;
+
+ throw UnknownDerivIDException();
+}
+
+int
+HeterogeneousModel::getJacobianColsNbr(bool sparse) const
+{
+ assert(sparse);
+ return 3 * (symbol_table.het_endo_nbr(heterogeneity_dimension) + symbol_table.endo_nbr())
+ + symbol_table.het_exo_nbr(heterogeneity_dimension) + symbol_table.exo_nbr();
+}
+
+SymbolType
+HeterogeneousModel::getTypeByDerivID(int deriv_id) const noexcept(false)
+{
+ return symbol_table.getType(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getLagByDerivID(int deriv_id) const noexcept(false)
+{
+ if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+ throw UnknownDerivIDException();
+
+ return inv_deriv_id_table[deriv_id].second;
+}
+
+int
+HeterogeneousModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
+{
+ if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+ throw UnknownDerivIDException();
+
+ return inv_deriv_id_table[deriv_id].first;
+}
+
+int
+HeterogeneousModel::getTypeSpecificIDByDerivID(int deriv_id) const
+{
+ return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getDerivID(int symb_id, int lead_lag) const noexcept(false)
+{
+ if (auto it = deriv_id_table.find({symb_id, lead_lag}); it == deriv_id_table.end())
+ throw UnknownDerivIDException();
+ else
+ return it->second;
+}
+
+void
+HeterogeneousModel::writeDriverOutput(ostream& output) const
+{
+ output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").dynamic_tmp_nbr = [";
+ for (const auto& it : temporary_terms_derivatives)
+ output << it.size() << "; ";
+ output << "];" << endl;
+ writeDriverSparseIndicesHelper(
+ "heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
+}
diff --git a/src/HeterogeneousModel.hh b/src/HeterogeneousModel.hh
new file mode 100644
index 00000000..d4eb7d88
--- /dev/null
+++ b/src/HeterogeneousModel.hh
@@ -0,0 +1,101 @@
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef HETEROGENEOUS_MODEL_HH
+#define HETEROGENEOUS_MODEL_HH
+
+#include <string>
+
+#include "ModelTree.hh"
+
+using namespace std;
+
+class HeterogeneousModel : public ModelTree
+{
+public:
+ const int heterogeneity_dimension;
+
+ HeterogeneousModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg, int heterogeneity_dimension_arg);
+
+ HeterogeneousModel(const HeterogeneousModel& m);
+ HeterogeneousModel& operator=(const HeterogeneousModel& m);
+
+ void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);
+
+ void writeModelFiles(const string& basename, bool julia) const;
+ void writeDriverOutput(ostream& output) const;
+
+ [[nodiscard]] int getJacobianCol(int deriv_id, bool sparse) const override;
+ [[nodiscard]] int getJacobianColsNbr(bool sparse) const override;
+
+#if 0
+ void substituteEndoLeadGreaterThanTwo();
+
+ //! Transforms the model by removing all lags greater or equal than 2 on endos
+ void substituteEndoLagGreaterThanTwo();
+
+ //! Transforms the model by removing all leads on exos
+ /*! Note that this can create new lags on endos and exos */
+ void substituteExoLead();
+
+ //! Transforms the model by removing all lags on exos
+ void substituteExoLag();
+
+ //! Transforms the model by removing all UnaryOpcode::expectation
+ void substituteExpectation(bool partial_information_model);
+
+ //! Transforms the model by decreasing the lead/lag of predetermined variables in model equations
+ //! by one
+ void transformPredeterminedVariables();
+
+ //! Substitutes out all model-local variables
+ void substituteModelLocalVariables();
+#endif
+
+ // FIXME: the following 5 functions are identical to those in DynamicModel. Factorization?
+ [[nodiscard]] int getDerivID(int symb_id, int lead_lag) const noexcept(false) override;
+ [[nodiscard]] SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
+ [[nodiscard]] int getLagByDerivID(int deriv_id) const noexcept(false) override;
+ [[nodiscard]] int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
+ [[nodiscard]] int getTypeSpecificIDByDerivID(int deriv_id) const override;
+
+protected:
+ void computeChainRuleJacobian() override;
+ int getBlockJacobianEndoCol(int blk, int var, int lead_lag) const override;
+ string
+ modelClassName() const override
+ {
+ return "dynamic model for heterogeneity dimension '"
+ + heterogeneity_table.getName(heterogeneity_dimension) + "'";
+ }
+ int getMFS() const override;
+
+private:
+ // Maps a pair (symbol ID, lead/lag) to a deriv ID
+ map<pair<int, int>, int> deriv_id_table;
+ // Maps a deriv ID to a pair (symbol ID, lead/lag)
+ vector<pair<int, int>> inv_deriv_id_table;
+
+ // Allocates the derivation IDs for all endogenous variables for this heterogeneity dimension
+ void computeDerivIDs();
+};
+
+#endif
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 4e70c100..9b8211f4 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -31,25 +31,48 @@
#include "Shocks.hh"
ModFile::ModFile(WarningConsolidation& warnings_arg) :
+ symbol_table {heterogeneity_table},
var_model_table {symbol_table},
trend_component_model_table {symbol_table},
var_expectation_model_table {symbol_table},
pac_model_table {symbol_table},
- expressions_tree {symbol_table, num_constants, external_functions_table},
- original_model {symbol_table, num_constants, external_functions_table,
- trend_component_model_table, var_model_table},
- dynamic_model {symbol_table, num_constants, external_functions_table,
- trend_component_model_table, var_model_table},
- trend_dynamic_model {symbol_table, num_constants, external_functions_table,
- trend_component_model_table, var_model_table},
- orig_ramsey_dynamic_model {symbol_table, num_constants, external_functions_table,
- trend_component_model_table, var_model_table},
- epilogue {symbol_table, num_constants, external_functions_table, trend_component_model_table,
+ expressions_tree {symbol_table, num_constants, external_functions_table, heterogeneity_table},
+ original_model {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
+ var_model_table},
+ dynamic_model {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
+ var_model_table},
+ trend_dynamic_model {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
+ var_model_table},
+ orig_ramsey_dynamic_model {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
+ var_model_table},
+ epilogue {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
var_model_table},
- static_model {symbol_table, num_constants, external_functions_table},
- steady_state_model {symbol_table, num_constants, external_functions_table, static_model},
+ static_model {symbol_table, num_constants, external_functions_table, heterogeneity_table},
+ steady_state_model {symbol_table, num_constants, external_functions_table, heterogeneity_table,
+ static_model},
warnings {warnings_arg}
{
+ heterogeneity_table.setSymbolTable(&symbol_table);
}
void
@@ -396,6 +419,120 @@ ModFile::checkPass(bool nostrict, bool stochastic)
exit(EXIT_FAILURE);
}
}
+
+ if (!heterogeneity_table.empty())
+ {
+ if (block)
+ {
+ cerr << "ERROR: the 'block' option of the 'model' block is not supported for "
+ "heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.check_present)
+ {
+ cerr << "ERROR: The 'check' command is not supported for heterogeneous models" << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.steady_present)
+ {
+ cerr << "ERROR: The 'steady' command is not supported for heterogeneous models" << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.perfect_foresight_solver_present)
+ {
+ cerr << "ERROR: The 'perfect_foresight_solver' command is not supported for "
+ "heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.perfect_foresight_with_expectation_errors_solver_present)
+ {
+ cerr << "ERROR: The 'perfect_foresight_with_expectation_errors_solver' command is not "
+ "supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.stoch_simul_present)
+ {
+ cerr << "ERROR: The 'stoch_simul' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.estimation_present)
+ {
+ cerr << "ERROR: The 'estimation' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.osr_present)
+ {
+ cerr << "ERROR: The 'osr' command is not supported for heterogeneous models" << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.osr_params_present)
+ {
+ cerr << "ERROR: The 'osr_params' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.optim_weights_present)
+ {
+ cerr << "ERROR: The 'optim_weights' block is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.ramsey_model_present)
+ {
+ cerr << "ERROR: The 'ramsey_model' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.discretionary_policy_present)
+ {
+ cerr << "ERROR: The 'discretionary_policy' command is not supported for heterogeneous "
+ "models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.planner_objective_present)
+ {
+ cerr << "ERROR: The 'planner_objective' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.extended_path_present)
+ {
+ cerr << "ERROR: The 'extended_path' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.identification_present)
+ {
+ cerr << "ERROR: The 'identification' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.sensitivity_present)
+ {
+ cerr << "ERROR: The 'sensitivity' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.mom_estimation_present)
+ {
+ cerr
+ << "ERROR: The 'methods_of_moments' command is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ if (mod_file_struct.occbin_constraints_present)
+ {
+ cerr << "ERROR: The 'occbin_constraints' block is not supported for heterogeneous models"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ }
}
void
@@ -529,9 +666,12 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
*/
if (linear)
orig_ramsey_dynamic_model = dynamic_model;
- DynamicModel ramsey_FOC_equations_dynamic_model {
- symbol_table, num_constants, external_functions_table, trend_component_model_table,
- var_model_table};
+ DynamicModel ramsey_FOC_equations_dynamic_model {symbol_table,
+ num_constants,
+ external_functions_table,
+ heterogeneity_table,
+ trend_component_model_table,
+ var_model_table};
ramsey_FOC_equations_dynamic_model = dynamic_model;
auto clone_if_not_null
= [&](expr_t e) { return e ? e->clone(ramsey_FOC_equations_dynamic_model) : nullptr; };
@@ -550,6 +690,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
// Must come after detrending of variables and Ramsey policy transformation
dynamic_model.substituteLogTransform();
+ if (!heterogeneity_table.empty())
+ dynamic_model.substituteAggregationOperators();
+
/* Create auxiliary vars for leads and lags greater than 2, on both endos and
exos. The transformation is not exactly the same on stochastic and
deterministic models, because there is no need to take into account the
@@ -702,6 +845,16 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
exit(EXIT_FAILURE);
}
}
+
+ for (int dim {0}; dim < heterogeneity_table.size(); dim++)
+ if (heterogeneous_models.at(dim).equation_number() != symbol_table.het_endo_nbr(dim))
+ {
+ cerr << "ERROR: There are " << heterogeneous_models.at(dim).equation_number()
+ << " equations but " << symbol_table.het_endo_nbr(dim)
+ << " endogenous variables in the model for heterogeneity dimension '"
+ << heterogeneity_table.getName(dim) << "'!" << endl;
+ exit(EXIT_FAILURE);
+ }
}
void
@@ -822,6 +975,9 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
// Those matrices can only be filled here, because we use derivatives
dynamic_model.fillVarModelTableMatrices();
+ for (auto& hm : heterogeneous_models)
+ hm.computingPass(mod_file_struct.order_option, no_tmp_terms, use_dll);
+
for (auto& statement : statements)
statement->computingPass(mod_file_struct);
@@ -864,6 +1020,11 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
auto plusfolder {DataTree::packageDir(basename)};
+ if (check_model_changes && !heterogeneity_table.empty())
+ {
+ cerr << "ERROR: the 'fast' option is not supported for heterogeneous models" << endl;
+ exit(EXIT_FAILURE);
+ }
bool hasModelChanged = !dynamic_model.isChecksumMatching(basename) || !check_model_changes;
if (hasModelChanged)
{
@@ -939,6 +1100,7 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
mOutputFile << "M_.parameter_used_with_lead_lag = true;" << endl;
symbol_table.writeOutput(mOutputFile);
+ heterogeneity_table.writeOutput(mOutputFile);
var_model_table.writeOutput(basename, mOutputFile);
trend_component_model_table.writeOutput(basename, mOutputFile);
@@ -950,6 +1112,10 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
<< ");" << endl
<< "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", "
<< symbol_table.exo_nbr() << ");" << endl;
+ for (int hd {0}; hd < heterogeneity_table.size(); hd++)
+ mOutputFile << "M_.heterogeneity(" << hd + 1 << ").Sigma_e = zeros("
+ << symbol_table.het_exo_nbr(hd) << ", " << symbol_table.het_exo_nbr(hd) << ");"
+ << endl;
if (mod_file_struct.calibrated_measurement_errors)
mOutputFile << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", "
@@ -1038,6 +1204,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
}
}
+ for (const auto& hm : heterogeneous_models)
+ hm.writeDriverOutput(mOutputFile);
+
if (onlymodel || gui)
for (const auto& statement : statements)
{
@@ -1186,6 +1355,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
epilogue.writeEpilogueFile(basename);
pac_model_table.writeTargetCoefficientsFile(basename);
+
+ for (const auto& hm : heterogeneous_models)
+ hm.writeModelFiles(basename, false);
}
}
@@ -1263,6 +1435,11 @@ ModFile::writeJsonOutputParsingCheck(const string& basename, JsonFileOutputType
symbol_table.writeJsonOutput(output);
output << ", ";
+ if (!heterogeneity_table.empty())
+ {
+ heterogeneity_table.writeJsonOutput(output);
+ output << ", ";
+ }
dynamic_model.writeJsonOutput(output);
output << ", ";
static_model.writeJsonOutput(output);
diff --git a/src/ModFile.hh b/src/ModFile.hh
index 22d50f9a..97423f16 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -30,6 +30,8 @@
#include "DynamicModel.hh"
#include "ExtendedPreprocessorTypes.hh"
#include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
+#include "HeterogeneousModel.hh"
#include "ModelEquationBlock.hh"
#include "NumericalConstants.hh"
#include "NumericalInitialization.hh"
@@ -46,6 +48,8 @@ class ModFile
{
public:
explicit ModFile(WarningConsolidation& warnings_arg);
+ // For heterogeneity dimensions
+ HeterogeneityTable heterogeneity_table;
//! Symbol table
SymbolTable symbol_table;
//! External Functions table
@@ -76,6 +80,8 @@ public:
StaticModel static_model;
//! Static model, as declared in the "steady_state_model" block if present
SteadyStateModel steady_state_model;
+ // Heterogeneous model blocks, ordered per heterogeneity dimension ID
+ vector<HeterogeneousModel> heterogeneous_models;
//! Option linear
bool linear {false};
diff --git a/src/ModelEquationBlock.cc b/src/ModelEquationBlock.cc
index 224d6db4..a2cd7070 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -26,8 +26,10 @@
PlannerObjective::PlannerObjective(SymbolTable& symbol_table_arg,
NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg) :
- StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg}
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg) :
+ StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg}
{
}
@@ -38,7 +40,7 @@ PlannerObjective::writeDriverOutput(ostream& output) const
for (const auto& it : temporary_terms_derivatives)
output << it.size() << "; ";
output << "];" << endl;
- writeDriverSparseIndicesHelper<false, true>(output);
+ writeDriverSparseIndicesHelper("objective", output);
}
void
@@ -51,9 +53,14 @@ PlannerObjective::computingPassBlock([[maybe_unused]] const eval_context_t& eval
OrigRamseyDynamicModel::OrigRamseyDynamicModel(
SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg, VarModelTable& var_model_table_arg) :
- DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
- trend_component_model_table_arg, var_model_table_arg}
+ DynamicModel {symbol_table_arg,
+ num_constants_arg,
+ external_functions_table_arg,
+ heterogeneity_table_arg,
+ trend_component_model_table_arg,
+ var_model_table_arg}
{
}
@@ -67,8 +74,10 @@ OrigRamseyDynamicModel::operator=(const DynamicModel& m)
SteadyStateModel::SteadyStateModel(SymbolTable& symbol_table_arg,
NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
const StaticModel& static_model_arg) :
- DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg},
+ DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg},
static_model {static_model_arg}
{
}
@@ -327,10 +336,15 @@ SteadyStateModel::writeJsonSteadyStateFile(ostream& output, bool transformComput
Epilogue::Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg) :
- DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
- trend_component_model_table_arg, var_model_table_arg}
+ DynamicModel {symbol_table_arg,
+ num_constants_arg,
+ external_functions_table_arg,
+ heterogeneity_table_arg,
+ trend_component_model_table_arg,
+ var_model_table_arg}
{
}
diff --git a/src/ModelEquationBlock.hh b/src/ModelEquationBlock.hh
index cd6844b1..9b943b24 100644
--- a/src/ModelEquationBlock.hh
+++ b/src/ModelEquationBlock.hh
@@ -30,7 +30,8 @@ class PlannerObjective : public StaticModel
{
public:
PlannerObjective(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg);
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg);
// NB: masks the method with the same name in StaticModel (not in a virtual fashion)
void writeDriverOutput(ostream& output) const;
@@ -50,6 +51,7 @@ class OrigRamseyDynamicModel : public DynamicModel
public:
OrigRamseyDynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg);
OrigRamseyDynamicModel& operator=(const DynamicModel& m);
@@ -75,6 +77,7 @@ private:
public:
SteadyStateModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
const StaticModel& static_model_arg);
SteadyStateModel(const SteadyStateModel& m);
@@ -108,6 +111,7 @@ private:
public:
Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg);
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 59413577..1f823ec1 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -144,8 +144,10 @@ ModelTree::copyHelper(const ModelTree& m)
}
ModelTree::ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
- DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
+ DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg, is_dynamic_arg},
derivatives(4),
NNZDerivatives(4, 0),
temporary_terms_derivatives(4)
@@ -2072,6 +2074,8 @@ ModelTree::waitForMEXCompilationWorkers()
void
ModelTree::computingPassBlock(const eval_context_t& eval_context, bool no_tmp_terms)
{
+ if (!heterogeneity_table.empty())
+ return;
if (!computeNonSingularNormalization(eval_context))
return;
auto [prologue, epilogue] = computePrologueAndEpilogue();
@@ -2148,3 +2152,33 @@ ModelTree::computeMCPEquationsReordering()
swap(mcp_equations_reordering[endo_id], *it);
}
}
+
+void
+ModelTree::writeDriverSparseIndicesHelper(const string& prefix, ostream& output) const
+{
+ // Write indices for the sparse Jacobian (both naive and CSC storage)
+ output << "M_." << prefix << "_g1_sparse_rowval = int32([";
+ for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
+ output << indices.first + 1 << ' ';
+ output << "]);" << endl << "M_." << prefix << "_g1_sparse_colval = int32([";
+ for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
+ output << indices.second + 1 << ' ';
+ output << "]);" << endl << "M_." << prefix << "_g1_sparse_colptr = int32([";
+ for (int it : jacobian_sparse_colptr)
+ output << it + 1 << ' ';
+ output << "]);" << endl;
+
+ // Write indices for the sparse higher-order derivatives
+ for (int i {2}; i <= computed_derivs_order; i++)
+ {
+ output << "M_." << prefix << "_g" << i << "_sparse_indices = int32([";
+ for (const auto& [vidx, d] : derivatives[i])
+ {
+ for (bool row_number {true}; // First element of vidx is row number
+ int it : vidx)
+ output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 1 << ' ';
+ output << ';' << endl;
+ }
+ output << "]);" << endl;
+ }
+}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index f6247744..4198c299 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -366,9 +366,8 @@ protected:
temporary_terms_t& temporary_terms_union) const;
/* Helper for writing sparse derivatives indices in MATLAB/Octave driver file.
- Also supports the planner objective through the corresponding boolean. */
- template<bool dynamic, bool objective>
- void writeDriverSparseIndicesHelper(ostream& output) const;
+ The “prefix†will be prepended to variable names to construct objects under M_. */
+ void writeDriverSparseIndicesHelper(const string& prefix, ostream& output) const;
// Helper for writing sparse derivatives indices in JSON
template<bool dynamic>
@@ -409,7 +408,8 @@ protected:
// Writes the sparse representation of the model in MATLAB/Octave
template<bool dynamic>
- void writeSparseModelMFiles(const string& basename) const;
+ void writeSparseModelMFiles(const string& basename,
+ const optional<int>& heterogeneous_dimension = nullopt) const;
// Writes and compiles the sparse representation of the model in C
template<bool dynamic>
@@ -659,7 +659,8 @@ private:
public:
ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
protected:
ModelTree(const ModelTree& m);
@@ -2264,41 +2265,6 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output)};
}
-template<bool dynamic, bool objective>
-void
-ModelTree::writeDriverSparseIndicesHelper(ostream& output) const
-{
- static_assert(!(objective && dynamic), "There is no such thing as a dynamic planner objective");
- // TODO: when C++20 support is complete, mark this constexpr
- const string model_name {objective ? "objective" : (dynamic ? "dynamic" : "static")};
-
- // Write indices for the sparse Jacobian (both naive and CSC storage)
- output << "M_." << model_name << "_g1_sparse_rowval = int32([";
- for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
- output << indices.first + 1 << ' ';
- output << "]);" << endl << "M_." << model_name << "_g1_sparse_colval = int32([";
- for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
- output << indices.second + 1 << ' ';
- output << "]);" << endl << "M_." << model_name << "_g1_sparse_colptr = int32([";
- for (int it : jacobian_sparse_colptr)
- output << it + 1 << ' ';
- output << "]);" << endl;
-
- // Write indices for the sparse higher-order derivatives
- for (int i {2}; i <= computed_derivs_order; i++)
- {
- output << "M_." << model_name << "_g" << i << "_sparse_indices = int32([";
- for (const auto& [vidx, d] : derivatives[i])
- {
- for (bool row_number {true}; // First element of vidx is row number
- int it : vidx)
- output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 1 << ' ';
- output << ';' << endl;
- }
- output << "]);" << endl;
- }
-}
-
template<bool dynamic>
void
ModelTree::writeJsonSparseIndicesHelper(ostream& output) const
@@ -2408,6 +2374,8 @@ template<bool dynamic>
void
ModelTree::writeSparseModelJuliaFiles(const string& basename) const
{
+ assert(heterogeneity_table.empty());
+
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper < dynamic
? ExprNodeOutputType::juliaSparseDynamicModel
: ExprNodeOutputType::juliaSparseStaticModel > ();
@@ -2518,7 +2486,8 @@ ModelTree::writeSparseModelJuliaFiles(const string& basename) const
template<bool dynamic>
void
-ModelTree::writeSparseModelMFiles(const string& basename) const
+ModelTree::writeSparseModelMFiles(const string& basename,
+ const optional<int>& heterogeneous_dimension) const
{
constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::matlabSparseDynamicModel
: ExprNodeOutputType::matlabSparseStaticModel};
@@ -2526,9 +2495,17 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
const filesystem::path m_dir {packageDir(basename) / "+sparse"};
// TODO: when C++20 support is complete, mark the following strings constexpr
- const string prefix {dynamic ? "dynamic_" : "static_"};
+ const string prefix {
+ (dynamic ? "dynamic_"s : "static_"s)
+ + (heterogeneous_dimension ? "het"s + to_string(*heterogeneous_dimension + 1) + "_"s : ""s)};
const string full_prefix {basename + ".sparse." + prefix};
- const string ss_arg {dynamic ? ", steady_state" : ""};
+ const string extra_args {(dynamic ? ", steady_state"s : ""s)
+ + (heterogeneous_dimension
+ ? ", yh, xh, paramsh"s
+ : (heterogeneity_table.empty() ? ""s : ", yagg"s))};
+ const int nextra_args {
+ static_cast<int>(dynamic)
+ + (heterogeneous_dimension ? 3 : static_cast<int>(!heterogeneity_table.empty()))};
size_t ttlen {temporary_terms_derivatives[0].size()};
@@ -2544,7 +2521,7 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
// Residuals (non-block)
open_file(m_dir / (prefix + "resid_tt.m"));
- output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << ss_arg
+ output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if T_order >= 0" << endl
<< " return" << endl
@@ -2557,13 +2534,13 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
output.close();
open_file(m_dir / (prefix + "resid.m"));
- output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << ss_arg
+ output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << extra_args
<< ", T_order, T)" << endl
- << "if nargin < " << 5 + static_cast<int>(dynamic) << endl
+ << "if nargin < " << 5 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
- << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
+ << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "residual = NaN(" << equations.size() << ", 1);" << endl
<< d_sparse_output[0].str();
@@ -2574,12 +2551,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
ttlen += temporary_terms_derivatives[1].size();
open_file(m_dir / (prefix + "g1_tt.m"));
- output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << ss_arg
+ output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if T_order >= 1" << endl
<< " return" << endl
<< "end" << endl
- << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
+ << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "T_order = 1;" << endl
<< "if size(T, 1) < " << ttlen << endl
@@ -2590,14 +2567,14 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
open_file(m_dir / (prefix + "g1.m"));
// NB: At first order, sparse indices are passed as extra arguments
- output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << ss_arg
+ output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << extra_args
<< ", sparse_rowval, sparse_colval, sparse_colptr, T_order, T)" << endl
- << "if nargin < " << 8 + static_cast<int>(dynamic) << endl
+ << "if nargin < " << 8 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
- << "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T);"
- << endl
+ << "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << extra_args
+ << ", T_order, T);" << endl
<< "g1_v = NaN(" << jacobian_sparse_column_major_order.size() << ", 1);" << endl
<< d_sparse_output[1].str();
// On MATLAB < R2020a, sparse() does not accept int32 indices
@@ -2616,12 +2593,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
ttlen += temporary_terms_derivatives[i].size();
open_file(m_dir / (prefix + "g" + to_string(i) + "_tt.m"));
- output << "function [T_order, T] = " << prefix << "g" << i << "_tt(y, x, params" << ss_arg
+ output << "function [T_order, T] = " << prefix << "g" << i << "_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if T_order >= " << i << endl
<< " return" << endl
<< "end" << endl
- << "[T_order, T] = " << full_prefix << "g" << i - 1 << "_tt(y, x, params" << ss_arg
+ << "[T_order, T] = " << full_prefix << "g" << i - 1 << "_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "T_order = " << i << ";" << endl
<< "if size(T, 1) < " << ttlen << endl
@@ -2632,12 +2609,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
open_file(m_dir / (prefix + "g" + to_string(i) + ".m"));
output << "function [g" << i << "_v, T_order, T] = " << prefix << "g" << i << "(y, x, params"
- << ss_arg << ", T_order, T)" << endl
- << "if nargin < " << 5 + static_cast<int>(dynamic) << endl
+ << extra_args << ", T_order, T)" << endl
+ << "if nargin < " << 5 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
- << "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << ss_arg
+ << "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "g" << i << "_v = NaN(" << derivatives[i].size() << ", 1);" << endl
<< d_sparse_output[i].str() << "end" << endl;
@@ -2658,7 +2635,7 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
const string resid_g1_arg {evaluate ? "" : ", residual, g1"};
open_file(block_dir / (funcname + ".m"));
output << "function [y, T" << resid_g1_arg << "] = " << funcname << "(y, x, params"
- << ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
+ << extra_args << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
if (!evaluate)
output << "residual=NaN(" << blocks[blk].mfs_size << ", 1);" << endl;
@@ -2706,8 +2683,13 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src" / "sparse"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
- const string ss_argin {dynamic ? ", const double *restrict steady_state" : ""};
- const string ss_argout {dynamic ? ", steady_state" : ""};
+ const string extra_argin {
+ (dynamic ? ", const double *restrict steady_state"s : ""s)
+ + (heterogeneity_table.empty() ? ""s : ", const double *restrict yagg"s)};
+ const string extra_argout {(dynamic ? ", steady_state"s : ""s)
+ + (heterogeneity_table.empty() ? ""s : ", yagg"s)};
+ const int nextra_args {static_cast<int>(dynamic)
+ + static_cast<int>(!heterogeneity_table.empty())};
const int ylen {(dynamic ? 3 : 1) * symbol_table.endo_nbr()};
const int xlen {symbol_table.exo_nbr() + symbol_table.exo_det_nbr()};
@@ -2725,8 +2707,8 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
size_t ttlen {0};
- // Helper for dealing with y, x, params and steady_state inputs (shared with block case)
- auto y_x_params_ss_inputs = [&](bool assign_y) {
+ // Helper for dealing with y, x, params, steady_state and yagg inputs (shared with block case)
+ auto y_x_params_ss_yagg_inputs = [&](bool assign_y) {
output << " if (!(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0]) && !mxIsSparse(prhs[0]) && "
"mxGetNumberOfElements(prhs[0]) == "
<< ylen << "))" << endl
@@ -2753,12 +2735,22 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
<< R"( mexErrMsgTxt("steady_state must be a real dense numeric array with )"
<< symbol_table.endo_nbr() << R"( elements");)" << endl
<< " const double *restrict steady_state = mxGetPr(prhs[3]);" << endl;
+ if (!heterogeneity_table.empty())
+ {
+ const int idx {3 + static_cast<int>(dynamic)};
+ output << " if (!(mxIsDouble(prhs[" << idx << "]) && !mxIsComplex(prhs[" << idx
+ << "]) && !mxIsSparse(prhs[" << idx << "]) && mxGetNumberOfElements(prhs[" << idx
+ << "]) == " << heterogeneity_table.aggregateEndoSize() << "))" << endl
+ << R"( mexErrMsgTxt("yagg must be a real dense numeric array with )"
+ << heterogeneity_table.aggregateEndoSize() << R"( elements");)" << endl
+ << " const double *restrict yagg = mxGetPr(prhs[" << idx << "]);" << endl;
+ }
};
// Helper for dealing with sparse_rowval and sparse_colptr inputs (shared with block case)
auto sparse_indices_inputs = [&](int ncols, int nzval) {
// We use sparse_rowval and sparse_colptr (sparse_colval is unused)
- const int row_idx {3 + static_cast<int>(dynamic)}, col_idx {row_idx + 2};
+ const int row_idx {3 + nextra_args}, col_idx {row_idx + 2};
output << " if (!(mxIsInt32(prhs[" << row_idx << "]) && mxGetNumberOfElements(prhs[" << row_idx
<< "]) == " << nzval << "))" << endl
<< R"( mexErrMsgTxt("sparse_rowval must be an int32 array with )" << nzval
@@ -2800,7 +2792,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const string prototype_tt {
"void " + funcname
+ "_tt(const double *restrict y, const double *restrict x, const double *restrict params"
- + ss_argin + ", double *restrict T)"};
+ + extra_argin + ", double *restrict T)"};
const filesystem::path header_tt {model_src_dir / (funcname + "_tt.h")};
open_file(header_tt);
@@ -2825,7 +2817,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const string prototype_main {
"void " + funcname
+ "(const double *restrict y, const double *restrict x, const double *restrict params"
- + ss_argin + ", const double *restrict T, double *restrict "
+ + extra_argin + ", const double *restrict T, double *restrict "
+ (i == 0 ? "residual" : "g" + to_string(i) + "_v") + ")"};
const filesystem::path header_main {model_src_dir / (funcname + ".h")};
@@ -2850,7 +2842,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
compileMEX(model_src_dir, funcname, mexext, {source_main}, matlabroot, false)};
const filesystem::path source_mex {model_src_dir / (funcname + "_mex.c")};
- int nargin {5 + static_cast<int>(dynamic) + 3 * static_cast<int>(i == 1)};
+ int nargin {5 + nextra_args + 3 * static_cast<int>(i == 1)};
open_file(source_mex);
output << "#include <string.h>" << endl // For memcpy()
<< R"(#include "mex.h")" << endl
@@ -2870,7 +2862,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
<< " if (nlhs != 1 && nlhs != 3)" << endl
<< R"( mexErrMsgTxt("Accepts exactly 1 or 3 output arguments");)" << endl;
- y_x_params_ss_inputs(true);
+ y_x_params_ss_yagg_inputs(true);
if (i == 1)
sparse_indices_inputs(getJacobianColsNbr(true), jacobian_sparse_column_major_order.size());
@@ -2918,7 +2910,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
output << " default:" << endl << " " << prefix << "resid";
else
output << " case " << j - 1 << ":" << endl << " " << prefix << "g" << j;
- output << "_tt(y, x, params" << ss_argout << ", T);" << endl;
+ output << "_tt(y, x, params" << extra_argout << ", T);" << endl;
}
output << " }" << endl;
if (i == 1)
@@ -2928,7 +2920,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
output << " plhs[0] = mxCreateDoubleMatrix("
<< (i == 0 ? equations.size() : derivatives[i].size()) << ", 1, mxREAL);" << endl;
output << " " << prefix << (i == 0 ? "resid" : "g" + to_string(i)) << "(y, x, params"
- << ss_argout << ", T, mxGetPr(plhs[0]));" << endl
+ << extra_argout << ", T, mxGetPr(plhs[0]));" << endl
<< " if (nlhs == 3)" << endl
<< " {" << endl
<< " plhs[1] = T_order_mx;" << endl
@@ -2981,7 +2973,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
evaluate the recursive variables. */
output << " if (nlhs < 2 || nlhs > 4)" << endl
<< R"( mexErrMsgTxt("Accepts 2 to 4 output arguments");)" << endl;
- y_x_params_ss_inputs(false);
+ y_x_params_ss_yagg_inputs(false);
/* We’d like to avoid copying y if this is a “solve†block without
recursive variables. Unfortunately “plhs[0]=prhs[0]†leads to
@@ -3054,7 +3046,6 @@ ModelTree::writeDebugModelMFiles(const string& basename) const
const filesystem::path m_dir {packageDir(basename) / "+debug"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
- const string ss_arg {dynamic ? ", steady_state" : ""};
const filesystem::path resid_filename {m_dir / (prefix + "resid.m")};
ofstream output {resid_filename, ios::out | ios::binary};
@@ -3064,7 +3055,13 @@ ModelTree::writeDebugModelMFiles(const string& basename) const
exit(EXIT_FAILURE);
}
- output << "function [lhs, rhs] = " << prefix << "resid(y, x, params" << ss_arg << ")" << endl
+ output << "function [lhs, rhs] = " << prefix << "resid(y, x, params";
+ if (dynamic)
+ output << ", steady_state";
+ if (!heterogeneity_table.empty())
+ output << ", yagg";
+
+ output << ")" << endl
<< "T = NaN(" << temporary_terms_derivatives[0].size() << ", 1);" << endl
<< "lhs = NaN(" << equations.size() << ", 1);" << endl
<< "rhs = NaN(" << equations.size() << ", 1);" << endl;
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 4cf55d3e..7f29ae96 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -72,6 +72,7 @@ ParsingDriver::set_current_data_tree(DataTree* data_tree_arg)
data_tree = data_tree_arg;
model_tree = dynamic_cast<ModelTree*>(data_tree_arg);
dynamic_model = dynamic_cast<DynamicModel*>(data_tree_arg);
+ heterogeneous_model = dynamic_cast<HeterogeneousModel*>(data_tree_arg);
}
void
@@ -184,12 +185,14 @@ ParsingDriver::warning(const string& m)
int
ParsingDriver::declare_symbol(const string& name, SymbolType type, const string& tex_name,
- const vector<pair<string, string>>& partition_value)
+ const vector<pair<string, string>>& partition_value,
+ const optional<int>& heterogeneity_dimension)
{
int symb_id;
try
{
- symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value);
+ symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value,
+ heterogeneity_dimension);
}
catch (SymbolTable::AlreadyDeclaredException& e)
{
@@ -208,16 +211,30 @@ int
ParsingDriver::declare_endogenous(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
- return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value);
+ return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value, {});
}
void
ParsingDriver::var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
- bool log_option)
+ const optional<string>& heterogeneity_dimension, bool log_option)
{
for (auto& [name, tex_name, partition] : symbol_list)
{
- int symb_id = declare_endogenous(name, tex_name, partition);
+ int symb_id {[&] {
+ if (heterogeneity_dimension)
+ try
+ {
+ return declare_symbol(name, SymbolType::heterogeneousEndogenous, tex_name, partition,
+ mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+ }
+ catch (HeterogeneityTable::UnknownDimensionNameException&)
+ {
+ error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+ }
+ else
+ return declare_endogenous(name, tex_name, partition);
+ }()};
+
if (log_option)
mod_file->symbol_table.markWithLogTransform(symb_id);
}
@@ -227,15 +244,27 @@ int
ParsingDriver::declare_exogenous(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
- return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value);
+ return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value, {});
}
void
ParsingDriver::varexo(
- const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
+ const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+ const optional<string>& heterogeneity_dimension)
{
for (auto& [name, tex_name, partition] : symbol_list)
- declare_exogenous(name, tex_name, partition);
+ if (heterogeneity_dimension)
+ try
+ {
+ declare_symbol(name, SymbolType::heterogeneousExogenous, tex_name, partition,
+ mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+ }
+ catch (HeterogeneityTable::UnknownDimensionNameException&)
+ {
+ error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+ }
+ else
+ declare_exogenous(name, tex_name, partition);
}
void
@@ -243,22 +272,34 @@ ParsingDriver::varexo_det(
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
{
for (auto& [name, tex_name, partition] : symbol_list)
- declare_symbol(name, SymbolType::exogenousDet, tex_name, partition);
+ declare_symbol(name, SymbolType::exogenousDet, tex_name, partition, {});
}
int
ParsingDriver::declare_parameter(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
- return declare_symbol(name, SymbolType::parameter, tex_name, partition_value);
+ return declare_symbol(name, SymbolType::parameter, tex_name, partition_value, {});
}
void
ParsingDriver::parameters(
- const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
+ const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+ const optional<string>& heterogeneity_dimension)
{
for (auto& [name, tex_name, partition] : symbol_list)
- declare_parameter(name, tex_name, partition);
+ if (heterogeneity_dimension)
+ try
+ {
+ declare_symbol(name, SymbolType::heterogeneousParameter, tex_name, partition,
+ mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
+ }
+ catch (HeterogeneityTable::UnknownDimensionNameException&)
+ {
+ error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
+ }
+ else
+ declare_parameter(name, tex_name, partition);
}
void
@@ -267,7 +308,7 @@ ParsingDriver::declare_statement_local_variable(const string& name)
if (mod_file->symbol_table.exists(name))
error("Symbol " + name + " cannot be assigned within a statement "
+ "while being assigned elsewhere in the modfile");
- declare_symbol(name, SymbolType::statementDeclaredVariable, "", {});
+ declare_symbol(name, SymbolType::statementDeclaredVariable, "", {}, {});
}
void
@@ -294,7 +335,7 @@ ParsingDriver::end_trend_var(bool log_trend, expr_t growth_factor,
for (auto& [name, tex_name] : symbol_list)
{
int symb_id = declare_symbol(name, log_trend ? SymbolType::logTrend : SymbolType::trend,
- tex_name, {});
+ tex_name, {}, {});
declared_trend_vars.push_back(symb_id);
}
@@ -837,7 +878,7 @@ ParsingDriver::end_epilogue()
void
ParsingDriver::add_epilogue_variable(const string& name)
{
- declare_symbol(name, SymbolType::epilogue, "", {});
+ declare_symbol(name, SymbolType::epilogue, "", {}, {});
}
void
@@ -852,6 +893,22 @@ ParsingDriver::begin_model()
set_current_data_tree(&mod_file->dynamic_model);
}
+void
+ParsingDriver::begin_heterogeneous_model(const string& heterogeneity_dimension)
+{
+ int het_dim_id {[&] {
+ try
+ {
+ return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
+ }
+ catch (HeterogeneityTable::UnknownDimensionNameException&)
+ {
+ error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
+ }
+ }()};
+ set_current_data_tree(&mod_file->heterogeneous_models.at(het_dim_id));
+}
+
void
ParsingDriver::end_model()
{
@@ -901,6 +958,28 @@ ParsingDriver::end_shocks(bool overwrite)
corr_shocks.clear();
}
+void
+ParsingDriver::end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite)
+{
+ int het_dim_id {[&] {
+ try
+ {
+ return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
+ }
+ catch (HeterogeneityTable::UnknownDimensionNameException&)
+ {
+ error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
+ }
+ }()};
+ mod_file->addStatement(make_unique<HeterogeneousShocksStatement>(
+ het_dim_id, overwrite, move(var_shocks), move(std_shocks), move(covar_shocks),
+ move(corr_shocks), mod_file->symbol_table, mod_file->heterogeneity_table));
+ var_shocks.clear();
+ std_shocks.clear();
+ covar_shocks.clear();
+ corr_shocks.clear();
+}
+
void
ParsingDriver::end_mshocks(bool overwrite, bool relative_to_initval)
{
@@ -2221,7 +2300,8 @@ void
ParsingDriver::begin_planner_objective()
{
planner_objective = make_unique<PlannerObjective>(mod_file->symbol_table, mod_file->num_constants,
- mod_file->external_functions_table);
+ mod_file->external_functions_table,
+ mod_file->heterogeneity_table);
set_current_data_tree(planner_objective.get());
}
@@ -2798,7 +2878,7 @@ void
ParsingDriver::model_local_variable(const vector<pair<string, string>>& symbol_list)
{
for (auto& [name, tex_name] : symbol_list)
- declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {});
+ declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {}, {});
}
void
@@ -3217,6 +3297,23 @@ ParsingDriver::add_steady_state(expr_t arg1)
return data_tree->AddSteadyState(arg1);
}
+expr_t
+ParsingDriver::add_sum(expr_t arg)
+{
+ if (heterogeneous_model)
+ error("The SUM() operator cannot be used inside a model(heterogeneity=...) block");
+
+ VariableNode* varg {dynamic_cast<VariableNode*>(arg)};
+ if (!varg)
+ error("The argument to the SUM() operator must be a single variable");
+ if (varg->lag != 0)
+ error("The argument to the SUM() operator must not have a lead or lag");
+ if (mod_file->symbol_table.getType(varg->symb_id) != SymbolType::heterogeneousEndogenous)
+ error("The argument to the SUM() operator must be a heterogeneous endogenous variable");
+
+ return data_tree->AddSum(arg);
+}
+
void
ParsingDriver::external_function_option(const string& name_option, const string& opt)
{
@@ -3225,7 +3322,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
if (opt.empty())
error("An argument must be passed to the 'name' option of the external_function() "
"statement.");
- declare_symbol(opt, SymbolType::externalFunction, "", {});
+ declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_id = mod_file->symbol_table.getID(opt);
}
else if (name_option == "first_deriv_provided")
@@ -3235,7 +3332,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
= ExternalFunctionsTable::IDSetButNoNameProvided;
else
{
- int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
+ int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_options.firstDerivSymbID = symb_id;
}
}
@@ -3246,7 +3343,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
= ExternalFunctionsTable::IDSetButNoNameProvided;
else
{
- int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
+ int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_options.secondDerivSymbID = symb_id;
}
}
@@ -3412,7 +3509,7 @@ ParsingDriver::add_model_var_or_external_function(const string& function_name, b
+ ") within the model block, you must first declare it via the "
"external_function() statement.");
}
- int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {});
+ int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {}, {});
current_external_function_options.nargs = stack_external_function_args.top().size();
mod_file->external_functions_table.addExternalFunction(
symb_id, current_external_function_options, in_model_block);
@@ -3838,7 +3935,8 @@ ParsingDriver::begin_occbin_constraints()
if added to the main DynamicModel tree. It also simplifies the
enforcement of various constraints at parsing time. */
occbin_constraints_tree = make_unique<DataTree>(mod_file->symbol_table, mod_file->num_constants,
- mod_file->external_functions_table, false);
+ mod_file->external_functions_table,
+ mod_file->heterogeneity_table, false);
set_current_data_tree(occbin_constraints_tree.get());
}
@@ -3973,3 +4071,26 @@ ParsingDriver::matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_we
{
mod_file->addStatement(make_unique<MatchedIrfsWeightsStatement>(move(weights), overwrite));
}
+
+void
+ParsingDriver::heterogeneity_dimension(const vector<string>& dims)
+{
+ for (const auto& dim : dims)
+ {
+ int het_dim_id {[&] {
+ try
+ {
+ return mod_file->heterogeneity_table.addDimension(dim);
+ }
+ catch (HeterogeneityTable::AlreadyDeclaredDimensionException&)
+ {
+ error("Heterogeneity dimension '" + dim + "' already declared");
+ }
+ }()};
+
+ assert(static_cast<int>(mod_file->heterogeneous_models.size()) == het_dim_id);
+ mod_file->heterogeneous_models.emplace_back(mod_file->symbol_table, mod_file->num_constants,
+ mod_file->external_functions_table,
+ mod_file->heterogeneity_table, het_dim_id);
+ }
+}
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index a8fba462..628cd5a0 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -114,7 +114,8 @@ private:
//! Helper to add a symbol declaration (returns its symbol ID)
int declare_symbol(const string& name, SymbolType type, const string& tex_name,
- const vector<pair<string, string>>& partition_value);
+ const vector<pair<string, string>>& partition_value,
+ const optional<int>& heterogeneity_dimension);
//! Temporary store for the planner objective
unique_ptr<PlannerObjective> planner_objective;
@@ -137,6 +138,11 @@ private:
* DynamicModel instance */
DynamicModel* dynamic_model;
+ //! The heterogeneous model tree in which to add expressions currently parsed
+ /*! It is only a dynamic cast of data_tree pointer, and is therefore null if data_tree is not a
+ * HeterogeneousModel instance */
+ HeterogeneousModel* heterogeneous_model;
+
//! Sets data_tree and model_tree pointers
void set_current_data_tree(DataTree* data_tree_arg);
@@ -371,19 +377,21 @@ public:
const vector<pair<string, string>>& partition_value = {});
// Handles a “var†or “var(log)†statement (without “deflator†or “log_deflator†options)
void var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
- bool log_option);
+ const optional<string>& heterogeneity_dimension, bool log_option);
//! Declares an exogenous variable (and returns its symbol ID)
int declare_exogenous(const string& name, const string& tex_name = "",
const vector<pair<string, string>>& partition_value = {});
// Handles a “varexo†statement
- void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
+ void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+ const optional<string>& heterogeneity_dimension);
// Handles a “varexo_det†statement
void varexo_det(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
//! Declares a parameter (and returns its symbol ID)
int declare_parameter(const string& name, const string& tex_name = "",
const vector<pair<string, string>>& partition_value = {});
// Handles a “parameters†statement
- void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
+ void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
+ const optional<string>& heterogeneity_dimension);
// Handles a “model_local_variable†statement
void model_local_variable(const vector<pair<string, string>>& symbol_list);
//! Declares a statement local variable
@@ -454,10 +462,11 @@ public:
void add_epilogue_variable(const string& varname);
//! Add equation in epilogue block
void add_epilogue_equal(const string& varname, expr_t expr);
- /* Begin a model or model_replace block, or an expression as an option value
- of some statement.
- Must be followed by a call to reset_data_tree(). */
+ /* Begin a model block (without heterogeneity option), or an expression as an option value of some
+ statement. Must be followed by a call to reset_data_tree(). */
void begin_model();
+ // Begin a model(heterogeneity=…) block
+ void begin_heterogeneous_model(const string& heterogeneity_dimension);
//! End a model or model_replace block, printing errors that were encountered in parsing
void end_model();
//! Writes a shocks statement
@@ -468,6 +477,8 @@ public:
void end_shocks_surprise(bool overwrite);
//! Writes a shocks(learnt_in=…) block
void end_shocks_learnt_in(const string& learnt_in_period, bool overwrite);
+ // For a shocks(heterogeneity=…) block
+ void end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite);
//! Writes a mshocks(learnt_in=…) block
void end_mshocks_learnt_in(const string& learnt_in_period, bool overwrite,
bool relative_to_initval);
@@ -844,6 +855,8 @@ public:
expr_t add_erfc(expr_t arg);
//! Writes token "steadyState(arg1)" to model tree
expr_t add_steady_state(expr_t arg1);
+ // Add a “sum(arg)†node to model tree
+ expr_t add_sum(expr_t arg);
//! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)
void push_external_function_arg_vector_onto_stack();
//! Adds an external function argument
@@ -952,6 +965,8 @@ public:
// Add a matched_irfs_weights block
void matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_weights_t weights,
bool overwrite);
+ void heterogeneity_dimension(const vector<string>& dims);
+
// Returns true iff the string is a legal symbol identifier (see NAME token in lexer)
static bool isSymbolIdentifier(const string& str);
// Given an Occbin regime name, returns the corresponding auxiliary parameter
diff --git a/src/Shocks.cc b/src/Shocks.cc
index 1b7f54db..963dfd33 100644
--- a/src/Shocks.cc
+++ b/src/Shocks.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -592,6 +592,213 @@ ShocksLearntInStatement::writeJsonOutput(ostream& output) const
output << "]}";
}
+HeterogeneousShocksStatement::HeterogeneousShocksStatement(
+ int heterogeneity_dimension_arg, bool overwrite_arg, var_and_std_shocks_t var_shocks_arg,
+ var_and_std_shocks_t std_shocks_arg, covar_and_corr_shocks_t covar_shocks_arg,
+ covar_and_corr_shocks_t corr_shocks_arg, const SymbolTable& symbol_table_arg,
+ const HeterogeneityTable& heterogeneity_table_arg) :
+ heterogeneity_dimension {heterogeneity_dimension_arg},
+ overwrite {overwrite_arg},
+ var_shocks {move(var_shocks_arg)},
+ std_shocks {move(std_shocks_arg)},
+ covar_shocks {move(covar_shocks_arg)},
+ corr_shocks {move(corr_shocks_arg)},
+ symbol_table {symbol_table_arg},
+ heterogeneity_table {heterogeneity_table_arg}
+{
+}
+
+void
+HeterogeneousShocksStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
+ [[maybe_unused]] bool minimal_workspace) const
+{
+ if (overwrite)
+ output << sigmaeName() << " = zeros(" << symbol_table.het_exo_nbr(heterogeneity_dimension)
+ << ", " << symbol_table.het_exo_nbr(heterogeneity_dimension) << ");" << endl;
+
+ writeVarAndStdShocks(output);
+ writeCovarAndCorrShocks(output);
+}
+
+void
+HeterogeneousShocksStatement::writeJsonOutput(ostream& output) const
+{
+ output << R"({"statementName": "shocks")"
+ << R"(, "heterogeneity": ")" << heterogeneity_table.getName(heterogeneity_dimension)
+ << R"(", "overwrite": )" << boolalpha << overwrite;
+ output << R"(, "variance": [)";
+ for (bool printed_something {false}; auto& [id, value] : var_shocks)
+ {
+ if (exchange(printed_something, true))
+ output << ", ";
+ output << R"({"name": ")" << symbol_table.getName(id) << R"(", )"
+ << R"("variance": ")";
+ value->writeJsonOutput(output, {}, {});
+ output << R"("})";
+ }
+ output << "]"
+ << R"(, "stderr": [)";
+ for (bool printed_something {false}; auto& [id, value] : std_shocks)
+ {
+ if (exchange(printed_something, true))
+ output << ", ";
+ output << R"({"name": ")" << symbol_table.getName(id) << R"(", )"
+ << R"("stderr": ")";
+ value->writeJsonOutput(output, {}, {});
+ output << R"("})";
+ }
+ output << "]"
+ << R"(, "covariance": [)";
+ for (bool printed_something {false}; auto& [ids, value] : covar_shocks)
+ {
+ if (exchange(printed_something, true))
+ output << ", ";
+ output << "{"
+ << R"("name": ")" << symbol_table.getName(ids.first) << R"(", )"
+ << R"("name2": ")" << symbol_table.getName(ids.second) << R"(", )"
+ << R"("covariance": ")";
+ value->writeJsonOutput(output, {}, {});
+ output << R"("})";
+ }
+ output << "]"
+ << R"(, "correlation": [)";
+ for (bool printed_something {false}; auto& [ids, value] : corr_shocks)
+ {
+ if (exchange(printed_something, true))
+ output << ", ";
+ output << "{"
+ << R"("name": ")" << symbol_table.getName(ids.first) << R"(", )"
+ << R"("name2": ")" << symbol_table.getName(ids.second) << R"(", )"
+ << R"("correlation": ")";
+ value->writeJsonOutput(output, {}, {});
+ output << R"("})";
+ }
+ output << "]"
+ << "}";
+}
+
+void
+HeterogeneousShocksStatement::writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it,
+ bool stddev) const
+{
+ SymbolType type = symbol_table.getType(it.first);
+ assert(type == SymbolType::heterogeneousExogenous);
+
+ int id {symbol_table.getTypeSpecificID(it.first) + 1};
+ output << sigmaeName() << "(" << id << ", " << id << ") = ";
+ if (stddev)
+ output << "(";
+ it.second->writeOutput(output);
+ if (stddev)
+ output << ")^2";
+ output << ";" << endl;
+}
+
+void
+HeterogeneousShocksStatement::writeVarAndStdShocks(ostream& output) const
+{
+ for (const auto& it : var_shocks)
+ writeVarOrStdShock(output, it, false);
+
+ for (const auto& it : std_shocks)
+ writeVarOrStdShock(output, it, true);
+}
+
+void
+HeterogeneousShocksStatement::writeCovarOrCorrShock(ostream& output,
+ const pair<pair<int, int>, expr_t>& it,
+ bool corr) const
+{
+ assert(symbol_table.getType(it.first.first) == SymbolType::heterogeneousExogenous
+ && symbol_table.getType(it.first.second) == SymbolType::heterogeneousExogenous);
+ int id1 {symbol_table.getTypeSpecificID(it.first.first) + 1},
+ id2 {symbol_table.getTypeSpecificID(it.first.second) + 1};
+
+ output << sigmaeName() << "(" << id1 << ", " << id2 << ") = ";
+ it.second->writeOutput(output);
+ if (corr)
+ output << "*sqrt(" << sigmaeName() << "(" << id1 << ", " << id1 << ")*" << sigmaeName() << "("
+ << id2 << ", " << id2 << "))";
+ output << ";" << endl
+ << sigmaeName() << "(" << id2 << ", " << id1 << ") = " << sigmaeName() << "(" << id1
+ << ", " << id2 << ");" << endl;
+}
+
+void
+HeterogeneousShocksStatement::writeCovarAndCorrShocks(ostream& output) const
+{
+ for (const auto& it : covar_shocks)
+ writeCovarOrCorrShock(output, it, false);
+
+ for (const auto& it : corr_shocks)
+ writeCovarOrCorrShock(output, it, true);
+}
+
+void
+HeterogeneousShocksStatement::checkPass(ModFileStructure& mod_file_struct,
+ [[maybe_unused]] WarningConsolidation& warnings)
+{
+ /* Error out if variables are not of the right type. This must be done here
+ and not at parsing time (see #448). */
+ for (auto [id, val] : var_shocks)
+ if (symbol_table.getType(id) != SymbolType::heterogeneousExogenous)
+ {
+ cerr << "shocks: setting a variance on '" << symbol_table.getName(id)
+ << "' is not allowed, because it is not a heterogeneous exogenous variable" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ for (auto [id, val] : std_shocks)
+ if (symbol_table.getType(id) != SymbolType::heterogeneousExogenous)
+ {
+ cerr << "shocks: setting a standard error on '" << symbol_table.getName(id)
+ << "' is not allowed, because it is not a heterogeneous exogenous variable" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ for (const auto& [ids, val] : covar_shocks)
+ {
+ auto& [symb_id1, symb_id2] = ids;
+
+ if (!(symbol_table.getType(symb_id1) == SymbolType::heterogeneousExogenous
+ && symbol_table.getType(symb_id2) == SymbolType::heterogeneousExogenous))
+ {
+ cerr << "shocks: setting a covariance between '" << symbol_table.getName(symb_id1)
+ << "' and '" << symbol_table.getName(symb_id2)
+ << "'is not allowed; covariances can only be specified for heterogeneous exogenous "
+ "variables"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ for (const auto& [ids, val] : corr_shocks)
+ {
+ auto& [symb_id1, symb_id2] = ids;
+
+ if (!(symbol_table.getType(symb_id1) == SymbolType::heterogeneousExogenous
+ && symbol_table.getType(symb_id2) == SymbolType::heterogeneousExogenous))
+ {
+ cerr << "shocks: setting a correlation between '" << symbol_table.getName(symb_id1)
+ << "' and '" << symbol_table.getName(symb_id2)
+ << "'is not allowed; covariances can only be specified for heterogeneous exogenous "
+ "variables"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ // Fill in mod_file_struct.parameters_with_shocks_values (related to #469)
+ for (auto [id, val] : var_shocks)
+ val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
+ for (auto [id, val] : std_shocks)
+ val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
+ for (const auto& [ids, val] : covar_shocks)
+ val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
+ for (const auto& [ids, val] : corr_shocks)
+ val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
+}
+
ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(
AbstractShocksStatement::det_shocks_t paths_arg, const SymbolTable& symbol_table_arg) :
paths {move(paths_arg)}, symbol_table {symbol_table_arg}, path_length {computePathLength(paths)}
diff --git a/src/Shocks.hh b/src/Shocks.hh
index 9378f246..a92ef733 100644
--- a/src/Shocks.hh
+++ b/src/Shocks.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -25,6 +25,7 @@
#include <vector>
#include "ExprNode.hh"
+#include "HeterogeneityTable.hh"
#include "Statement.hh"
#include "SymbolTable.hh"
@@ -155,6 +156,47 @@ public:
void writeJsonOutput(ostream& output) const override;
};
+class HeterogeneousShocksStatement : public Statement
+{
+public:
+ const int heterogeneity_dimension;
+ const bool overwrite;
+
+ using var_and_std_shocks_t = map<int, expr_t>;
+ using covar_and_corr_shocks_t = map<pair<int, int>, expr_t>;
+
+ const var_and_std_shocks_t var_shocks, std_shocks;
+ const covar_and_corr_shocks_t covar_shocks, corr_shocks;
+
+private:
+ const SymbolTable& symbol_table;
+ const HeterogeneityTable& heterogeneity_table;
+
+ void writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it, bool stddev) const;
+ void writeVarAndStdShocks(ostream& output) const;
+ void writeCovarOrCorrShock(ostream& output, const pair<pair<int, int>, expr_t>& it,
+ bool corr) const;
+ void writeCovarAndCorrShocks(ostream& output) const;
+
+ string
+ sigmaeName() const
+ {
+ return "M_.heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").Sigma_e"s;
+ }
+
+public:
+ HeterogeneousShocksStatement(int heterogeneity_dimension_arg, bool overwrite_arg,
+ var_and_std_shocks_t var_shocks_arg,
+ var_and_std_shocks_t std_shocks_arg,
+ covar_and_corr_shocks_t covar_shocks_arg,
+ covar_and_corr_shocks_t corr_shocks_arg,
+ const SymbolTable& symbol_table_arg,
+ const HeterogeneityTable& heterogeneity_table_arg);
+ void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
+ void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
+ void writeJsonOutput(ostream& output) const override;
+};
+
class ConditionalForecastPathsStatement : public Statement
{
private:
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index f11b23c9..780e3ed8 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -31,8 +31,10 @@
#include "StaticModel.hh"
StaticModel::StaticModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
- ExternalFunctionsTable& external_functions_table_arg) :
- ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg}
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg) :
+ ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+ heterogeneity_table_arg}
{
}
@@ -72,7 +74,7 @@ StaticModel::operator=(const StaticModel& m)
}
StaticModel::StaticModel(const DynamicModel& m) :
- ModelTree {m.symbol_table, m.num_constants, m.external_functions_table}
+ ModelTree {m.symbol_table, m.num_constants, m.external_functions_table, m.heterogeneity_table}
{
// Convert model local variables (need to be done first)
for (int it : m.local_variables_vector)
@@ -335,7 +337,7 @@ StaticModel::writeDriverOutput(ostream& output) const
if (block_decomposed)
writeBlockDriverOutput(output);
- writeDriverSparseIndicesHelper<false, false>(output);
+ writeDriverSparseIndicesHelper("static", output);
output << "M_.static_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index acadbb56..eb6ebe8b 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -121,7 +121,8 @@ protected:
public:
StaticModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants,
- ExternalFunctionsTable& external_functions_table_arg);
+ ExternalFunctionsTable& external_functions_table_arg,
+ HeterogeneityTable& heterogeneity_table_arg);
StaticModel(const StaticModel& m);
StaticModel& operator=(const StaticModel& m);
diff --git a/src/SymbolList.cc b/src/SymbolList.cc
index 40293ce7..98032ead 100644
--- a/src/SymbolList.cc
+++ b/src/SymbolList.cc
@@ -103,6 +103,15 @@ SymbolList::checkPass(WarningConsolidation& warnings, const vector<SymbolType>&
case SymbolType::excludedVariable:
valid_types += "excludedVariable, ";
break;
+ case SymbolType::heterogeneousEndogenous:
+ valid_types += "heterogeneousEndogenous, ";
+ break;
+ case SymbolType::heterogeneousExogenous:
+ valid_types += "heterogeneousExogenous, ";
+ break;
+ case SymbolType::heterogeneousParameter:
+ valid_types += "heterogeneousParameter, ";
+ break;
}
valid_types = valid_types.erase(valid_types.size() - 2, 2);
throw SymbolListException {"Variable " + symbol + " is not one of {" + valid_types + "}"};
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 0c265f1d..304b4faf 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -31,7 +31,8 @@
int
SymbolTable::addSymbol(const string& name, SymbolType type, const string& tex_name,
- const vector<pair<string, string>>& partition_value) noexcept(false)
+ const vector<pair<string, string>>& partition_value,
+ const optional<int>& heterogeneity_dimension) noexcept(false)
{
if (frozen)
throw FrozenException();
@@ -78,13 +79,20 @@ SymbolTable::addSymbol(const string& name, SymbolType type, const string& tex_na
pmv[it.first] = it.second;
partition_value_map[id] = pmv;
}
+
+ assert(!isHeterogeneous(type)
+ || (heterogeneity_dimension.has_value() && *heterogeneity_dimension >= 0
+ && *heterogeneity_dimension < heterogeneity_table.size()));
+ if (isHeterogeneous(type))
+ heterogeneity_dimensions.emplace(id, *heterogeneity_dimension);
+
return id;
}
int
SymbolTable::addSymbol(const string& name, SymbolType type) noexcept(false)
{
- return addSymbol(name, type, "", {});
+ return addSymbol(name, type, "", {}, {});
}
void
@@ -95,6 +103,10 @@ SymbolTable::freeze() noexcept(false)
frozen = true;
+ het_endo_ids.resize(heterogeneity_table.size());
+ het_exo_ids.resize(heterogeneity_table.size());
+ het_param_ids.resize(heterogeneity_table.size());
+
for (int i = 0; i < static_cast<int>(symbol_table.size()); i++)
{
int tsi;
@@ -116,6 +128,18 @@ SymbolTable::freeze() noexcept(false)
tsi = param_ids.size();
param_ids.push_back(i);
break;
+ case SymbolType::heterogeneousEndogenous:
+ tsi = het_endo_ids.at(heterogeneity_dimensions.at(i)).size();
+ het_endo_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+ break;
+ case SymbolType::heterogeneousExogenous:
+ tsi = het_exo_ids.at(heterogeneity_dimensions.at(i)).size();
+ het_exo_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+ break;
+ case SymbolType::heterogeneousParameter:
+ tsi = het_param_ids.at(heterogeneity_dimensions.at(i)).size();
+ het_param_ids.at(heterogeneity_dimensions.at(i)).push_back(i);
+ break;
default:
continue;
}
@@ -131,12 +155,18 @@ SymbolTable::unfreeze()
exo_ids.clear();
exo_det_ids.clear();
param_ids.clear();
+ het_endo_ids.clear();
+ het_exo_ids.clear();
+ het_param_ids.clear();
type_specific_ids.clear();
}
void
SymbolTable::changeType(int id, SymbolType newtype) noexcept(false)
{
+ // FIXME: implement switch to heterogeneous variable; dimension will have to be provided
+ assert(!isHeterogeneous(newtype));
+
if (frozen)
throw FrozenException();
@@ -146,7 +176,8 @@ SymbolTable::changeType(int id, SymbolType newtype) noexcept(false)
}
int
-SymbolTable::getID(SymbolType type, int tsid) const noexcept(false)
+SymbolTable::getID(SymbolType type, int tsid, const optional<int>& heterogeneity_dimension) const
+ noexcept(false)
{
if (!frozen)
throw NotYetFrozenException();
@@ -155,26 +186,44 @@ SymbolTable::getID(SymbolType type, int tsid) const noexcept(false)
{
case SymbolType::endogenous:
if (tsid < 0 || tsid >= static_cast<int>(endo_ids.size()))
- throw UnknownTypeSpecificIDException {tsid, type};
+ throw UnknownTypeSpecificIDException {tsid, type, {}};
else
return endo_ids[tsid];
case SymbolType::exogenous:
if (tsid < 0 || tsid >= static_cast<int>(exo_ids.size()))
- throw UnknownTypeSpecificIDException {tsid, type};
+ throw UnknownTypeSpecificIDException {tsid, type, {}};
else
return exo_ids[tsid];
case SymbolType::exogenousDet:
if (tsid < 0 || tsid >= static_cast<int>(exo_det_ids.size()))
- throw UnknownTypeSpecificIDException {tsid, type};
+ throw UnknownTypeSpecificIDException {tsid, type, {}};
else
return exo_det_ids[tsid];
case SymbolType::parameter:
if (tsid < 0 || tsid >= static_cast<int>(param_ids.size()))
- throw UnknownTypeSpecificIDException {tsid, type};
+ throw UnknownTypeSpecificIDException {tsid, type, {}};
else
return param_ids[tsid];
+ case SymbolType::heterogeneousEndogenous:
+ assert(heterogeneity_dimension.has_value());
+ if (tsid < 0 || tsid >= static_cast<int>(het_endo_ids.at(*heterogeneity_dimension).size()))
+ throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+ else
+ return het_endo_ids.at(*heterogeneity_dimension).at(tsid);
+ case SymbolType::heterogeneousExogenous:
+ assert(heterogeneity_dimension.has_value());
+ if (tsid < 0 || tsid >= static_cast<int>(het_exo_ids.at(*heterogeneity_dimension).size()))
+ throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+ else
+ return het_exo_ids.at(*heterogeneity_dimension).at(tsid);
+ case SymbolType::heterogeneousParameter:
+ assert(heterogeneity_dimension.has_value());
+ if (tsid < 0 || tsid >= static_cast<int>(het_param_ids.at(*heterogeneity_dimension).size()))
+ throw UnknownTypeSpecificIDException {tsid, type, *heterogeneity_dimension};
+ else
+ return het_param_ids.at(*heterogeneity_dimension).at(tsid);
default:
- throw UnknownTypeSpecificIDException {tsid, type};
+ throw UnknownTypeSpecificIDException {tsid, type, {}};
}
}
@@ -348,6 +397,7 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
case AuxVarType::expectation:
case AuxVarType::pacExpectation:
case AuxVarType::pacTargetNonstationary:
+ case AuxVarType::aggregationOp:
break;
case AuxVarType::endoLag:
case AuxVarType::exoLag:
@@ -419,6 +469,39 @@ SymbolTable::writeOutput(ostream& output) const noexcept(false)
output << getTypeSpecificID(varexob) + 1 << " ";
output << " ];" << endl;
}
+
+ // Heterogeneous symbols
+ // FIXME: the following helper could be used to simplify non-heterogenous variables
+ auto print_symb_names = [this, &output](const string& field, const auto& symb_ids) {
+ auto helper = [this, &output, &symb_ids](auto nameMethod) {
+ for (bool first_printed {false}; int symb_id : symb_ids)
+ {
+ if (exchange(first_printed, true))
+ output << "; ";
+ output << "'" << (this->*nameMethod)(symb_id) << "'";
+ }
+ };
+ output << field << " = {";
+ helper(&SymbolTable::getName);
+ output << "};" << endl << field << "_tex = {";
+ helper(&SymbolTable::getTeXName);
+ output << "};" << endl << field << "_long = {";
+ helper(&SymbolTable::getLongName);
+ output << "};" << endl;
+ };
+ for (int het_dim {0}; het_dim < heterogeneity_table.size(); het_dim++)
+ {
+ const string basefield {"M_.heterogeneity(" + to_string(het_dim + 1) + ")."};
+
+ output << basefield << "endo_nbr = " << het_endo_nbr(het_dim) << ";" << endl;
+ print_symb_names(basefield + "endo_names", het_endo_ids.at(het_dim));
+
+ 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));
+ }
}
int
@@ -714,6 +797,28 @@ SymbolTable::addPacTargetNonstationaryAuxiliaryVar(const string& name, expr_t ex
return symb_id;
}
+int
+SymbolTable::addAggregationOpAuxiliaryVar(const string& name, expr_t expr_arg)
+{
+ int symb_id {[&] {
+ try
+ {
+ return addSymbol(name, SymbolType::endogenous);
+ }
+ catch (AlreadyDeclaredException&)
+ {
+ cerr << "ERROR: the variable/parameter '" << name
+ << "' conflicts with a variable that will be generated for an aggregation operator. "
+ "Please rename it."
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+ }()};
+
+ aux_vars.emplace_back(symb_id, AuxVarType::aggregationOp, 0, 0, 0, 0, expr_arg, "");
+ return symb_id;
+}
+
int
SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const noexcept(false)
{
@@ -994,6 +1099,7 @@ SymbolTable::writeJsonOutput(ostream& output) const
case AuxVarType::expectation:
case AuxVarType::pacExpectation:
case AuxVarType::pacTargetNonstationary:
+ case AuxVarType::aggregationOp:
break;
case AuxVarType::endoLag:
case AuxVarType::exoLag:
@@ -1028,6 +1134,25 @@ SymbolTable::writeJsonOutput(ostream& output) const
}
output << "]" << endl;
}
+
+ if (!heterogeneity_table.empty())
+ {
+ output << R"(, "heterogeneous_symbols": [)";
+ for (int i {0}; i < heterogeneity_table.size(); i++)
+ {
+ if (i != 0)
+ output << ", ";
+ output << R"({ "dimension": ")" << heterogeneity_table.getName(i)
+ << R"(", "endogenous": )";
+ writeJsonVarVector(output, het_endo_ids.at(i));
+ output << R"(, "exogenous": )";
+ writeJsonVarVector(output, het_exo_ids.at(i));
+ output << R"(, "parameters": )";
+ writeJsonVarVector(output, het_param_ids.at(i));
+ output << "}";
+ }
+ output << "]" << endl;
+ }
}
void
@@ -1087,3 +1212,14 @@ SymbolTable::getLagrangeMultipliers() const
r.insert(aux_var.symb_id);
return r;
}
+
+int
+SymbolTable::getHeterogeneityDimension(int symb_id) const
+{
+ validateSymbID(symb_id);
+ auto it = heterogeneity_dimensions.find(symb_id);
+ if (it != heterogeneity_dimensions.end())
+ return it->second;
+ else
+ throw NonHeteregeneousSymbolException {symb_id};
+}
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 01bbf91f..3a3a1a3e 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -1,5 +1,5 @@
/*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
@@ -30,6 +30,7 @@
#include "CommonEnums.hh"
#include "ExprNode.hh"
+#include "HeterogeneityTable.hh"
using namespace std;
@@ -53,7 +54,9 @@ enum class AuxVarType
diffLead = 11, //!< Variable for timing between Diff operators (lead)
pacExpectation = 12, //!< Variable created for the substitution of the pac_expectation operator
pacTargetNonstationary
- = 13 //!< Variable created for the substitution of the pac_target_nonstationary operator
+ = 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()
};
//! Information on some auxiliary variables
@@ -106,6 +109,8 @@ struct AuxVarInfo
class SymbolTable
{
private:
+ HeterogeneityTable& heterogeneity_table;
+
//! Has method freeze() been called?
bool frozen {false};
@@ -123,6 +128,8 @@ private:
map<int, map<string, string>> partition_value_map;
//! Maps IDs to types
vector<SymbolType> type_table;
+ // Maps IDs of heterogenous symbols to heterogeneity dimension IDs
+ map<int, int> heterogeneity_dimensions;
//! Maps symbol IDs to type specific IDs
map<int, int> type_specific_ids;
@@ -135,6 +142,16 @@ private:
vector<int> exo_det_ids;
//! Maps type specific IDs of parameters to symbol IDs
vector<int> param_ids;
+ /* Maps type specific IDs of heterogeneous endogenous to symbol IDs (outer vector is for
+ heterogeneity dimensions) */
+ vector<vector<int>> het_endo_ids;
+ /* Maps type specific IDs of heterogeneous exogenous to symbol IDs (outer vector is for
+ heterogeneity dimensions) */
+ vector<vector<int>> het_exo_ids;
+ /* Maps type specific IDs of heterogeneous parameters to symbol IDs (outer vector is for
+ heterogeneity dimensions) */
+ vector<vector<int>> het_param_ids;
+
//! Information about auxiliary variables
vector<AuxVarInfo> aux_vars;
@@ -168,6 +185,7 @@ public:
{
const int tsid;
const SymbolType type;
+ const optional<int> heterogeneity_dimension;
};
/* Thrown when requesting the type specific ID of a symbol which doesn’t
have one */
@@ -202,6 +220,12 @@ public:
}
};
+ // Thrown by getHeterogeneityDimension() on non-heterogeneous symbols
+ struct NonHeteregeneousSymbolException
+ {
+ const int id;
+ };
+
private:
//! Factorized code for adding aux lag variables
int addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag,
@@ -214,11 +238,19 @@ private:
inline void validateSymbID(int symb_id) const noexcept(false);
public:
+ SymbolTable(HeterogeneityTable& heterogeneity_table_arg) :
+ heterogeneity_table {heterogeneity_table_arg}
+ {
+ }
+
//! Add a symbol
- /*! Returns the symbol ID */
+ /* Returns the symbol ID.
+ heterogeneity_dimension must be defined if this is a heterogeneous symbol (otherwise it is
+ ignored) */
int addSymbol(const string& name, SymbolType type, const string& tex_name,
- const vector<pair<string, string>>& partition_value) noexcept(false);
- //! Add a symbol without its TeX name (will be equal to its name)
+ const vector<pair<string, string>>& partition_value,
+ const optional<int>& heterogeneity_dimension) noexcept(false);
+ //! Add a (non-heterogenous) symbol without its TeX name (will be equal to its name)
/*! Returns the symbol ID */
int addSymbol(const string& name, SymbolType type) noexcept(false);
//! Adds an auxiliary variable for endogenous with lead >= 2
@@ -313,6 +345,8 @@ public:
int addPacExpectationAuxiliaryVar(const string& name, expr_t expr_arg);
//! An auxiliary variable for a pac_target_nonstationary operator
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);
//! Returns the number of auxiliary variables
[[nodiscard]] int
AuxVarsSize() const
@@ -340,7 +374,9 @@ public:
//! Get ID (by name)
[[nodiscard]] inline int getID(const string& name) const noexcept(false);
//! Get ID (by type specific ID)
- [[nodiscard]] int getID(SymbolType type, int tsid) const noexcept(false);
+ [[nodiscard]] int getID(SymbolType type, int tsid,
+ const optional<int>& heterogeneity_dimension = nullopt) const
+ noexcept(false);
//! Freeze symbol table
void freeze() noexcept(false);
//! unreeze symbol table
@@ -360,6 +396,12 @@ public:
[[nodiscard]] inline int exo_det_nbr() const noexcept(false);
//! Get number of parameters
[[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 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
+ [[nodiscard]] inline int het_param_nbr(int het_dim) const noexcept(false);
//! Returns the greatest symbol ID (the smallest is zero)
[[nodiscard]] inline int maxID() const;
//! Get number of user-declared endogenous variables (without the auxiliary variables)
@@ -423,6 +465,9 @@ public:
[[nodiscard]] const set<int>& getVariablesWithLogTransform() const;
// Returns all Lagrange multipliers
[[nodiscard]] set<int> getLagrangeMultipliers() const;
+ /* Get heterogeneity dimension of a given symbol. Throws NonHeterogeneousSymbolException
+ if there is no such dimension. */
+ [[nodiscard]] int getHeterogeneityDimension(int symb_id) const;
};
inline void
@@ -537,6 +582,33 @@ SymbolTable::param_nbr() const noexcept(false)
return param_ids.size();
}
+inline int
+SymbolTable::het_endo_nbr(int het_dim) const noexcept(false)
+{
+ if (!frozen)
+ throw NotYetFrozenException();
+
+ return het_endo_ids.at(het_dim).size();
+}
+
+inline int
+SymbolTable::het_exo_nbr(int het_dim) const noexcept(false)
+{
+ if (!frozen)
+ throw NotYetFrozenException();
+
+ return het_exo_ids.at(het_dim).size();
+}
+
+inline int
+SymbolTable::het_param_nbr(int het_dim) const noexcept(false)
+{
+ if (!frozen)
+ throw NotYetFrozenException();
+
+ return het_param_ids.at(het_dim).size();
+}
+
inline int
SymbolTable::maxID() const
{
diff --git a/src/meson.build b/src/meson.build
index 4dd87161..efebcfe7 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -53,6 +53,8 @@ preprocessor_src = [ 'ComputingTasks.cc',
'ModelEquationBlock.cc',
'WarningConsolidation.cc',
'SubModel.cc',
+ 'HeterogeneityTable.cc',
+ 'HeterogeneousModel.cc',
'macro/Driver.cc',
'macro/Environment.cc',
'macro/Expressions.cc',
--
GitLab