From 365fb27f3d48ed208332a52f0f96e255f54c2dd4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 7 Jul 2021 10:54:04 +0200
Subject: [PATCH] =?UTF-8?q?New=20=E2=80=9Cstructural=E2=80=9D=20option=20t?=
=?UTF-8?q?o=20=E2=80=9Cvar=5Fmodel=E2=80=9D?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
As the name implies, this option allows contemporaneous variables on the RHS.
The A0 matrix for contemporaneous variables is added as a second (optional)
output to the generated var_ar.m file. Note that for reduced-form VAR, this
matrix will be the identity.
Also, the user is now allowed to write the VAR models in a more flexible form:
the LHS must still be a single variable, but the RHS can be an arbitrary
expression (as long as it is linear, obviously). Internally, the preprocessor
now uses derivation to compute the coefficients of the AR and A0. This change
applies to both reduced-form and structural VAR models.
Ref. dynare#1785
---
src/DynamicModel.cc | 132 +++++++++++++++++++++++++++++++++++++++----
src/DynamicModel.hh | 18 +++++-
src/DynareBison.yy | 4 +-
src/DynareFlex.ll | 1 +
src/ExprNode.cc | 10 ++++
src/ExprNode.hh | 4 ++
src/ModFile.cc | 3 +
src/ParsingDriver.cc | 7 ++-
src/SubModel.cc | 37 ++++++++++--
src/SubModel.hh | 12 +++-
10 files changed, 205 insertions(+), 23 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index fa4cfb86..7beab035 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3402,9 +3402,6 @@ DynamicModel::fillVarModelTable() const
var_model_table.setLhs(lhsr);
var_model_table.setRhs(rhsr);
var_model_table.setLhsExprT(lhs_expr_tr);
-
- // Fill AR Matrix
- var_model_table.setAR(computeAutoregressiveMatrices(true));
}
void
@@ -3424,12 +3421,19 @@ DynamicModel::fillVarModelTableFromOrigModel() const
set<pair<int, int>> rhs_endo_set, rhs_exo_set;
equations[eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_endo_set);
for (const auto &[symb_id, lag] : rhs_endo_set)
- if (lag >= 0)
+ if (lag > 0)
+ {
+ cerr << "ERROR: in Equation " << eqtag
+ << ". A VAR model may not have leaded endogenous variables on the RHS. " << endl;
+ exit(EXIT_FAILURE);
+ }
+ else if (!var_model_table.getStructural().at(model_name) && lag == 0)
{
cerr << "ERROR: in Equation " << eqtag
- << ". A VAR model may not have leaded or contemporaneous endogenous variables on the RHS. " << endl;
+ << ". A non-structural VAR model may not have contemporaneous endogenous variables on the RHS. " << endl;
exit(EXIT_FAILURE);
}
+
equations[eqn]->arg2->collectDynamicVariables(SymbolType::exogenous, rhs_exo_set);
for (const auto &[symb_id, lag] : rhs_exo_set)
if (lag != 0)
@@ -3487,18 +3491,122 @@ DynamicModel::fillVarModelTableFromOrigModel() const
var_model_table.setOrigDiffVar(orig_diff_var);
}
+vector<int>
+DynamicModel::getVARDerivIDs(int lhs_symb_id, int lead_lag) const
+{
+ vector<int> deriv_ids;
+
+ // First directly look for the variable itself
+ if (auto it = deriv_id_table.find({ lhs_symb_id, lead_lag });
+ it != deriv_id_table.end())
+ deriv_ids.push_back(it->second);
+
+ // Then go through auxiliary variables
+ for (auto &[key, deriv_id2] : deriv_id_table)
+ {
+ auto [symb_id2, lead_lag2] = key;
+ const AuxVarInfo *avi;
+ try
+ {
+ avi = &symbol_table.getAuxVarInfo(symb_id2);
+ }
+ catch (SymbolTable::UnknownSymbolIDException)
+ {
+ continue;
+ }
+
+ if (avi->get_type() == AuxVarType::endoLag && avi->get_orig_symb_id() == lhs_symb_id
+ && avi->get_orig_lead_lag() + lead_lag2 == lead_lag)
+ deriv_ids.push_back(deriv_id2);
+
+ // Handle diff lag auxvar, possibly nested several times
+ int diff_lag_depth = 0;
+ while (avi->get_type() == AuxVarType::diffLag)
+ {
+ diff_lag_depth++;
+ if (avi->get_orig_symb_id() == lhs_symb_id && lead_lag2 - diff_lag_depth == lead_lag)
+ {
+ deriv_ids.push_back(deriv_id2);
+ break;
+ }
+ try
+ {
+ avi = &symbol_table.getAuxVarInfo(avi->get_orig_symb_id());
+ }
+ catch (SymbolTable::UnknownSymbolIDException)
+ {
+ break;
+ }
+ }
+ }
+
+ return deriv_ids;
+}
+
+void
+DynamicModel::fillVarModelTableMatrices()
+{
+ map<string, map<tuple<int, int, int>, expr_t>> AR;
+ map<string, map<tuple<int, int>, expr_t>> A0;
+ for (const auto &[model_name, eqns] : var_model_table.getEqNums())
+ {
+ const vector<int> &lhs = var_model_table.getLhs(model_name);
+ int max_lag = var_model_table.getMaxLag(model_name);
+ for (auto lhs_symb_id : lhs)
+ {
+ // Fill autoregressive matrix (AR)
+ for (int lag = 1; lag <= max_lag; lag++)
+ {
+ vector<int> deriv_ids = getVARDerivIDs(lhs_symb_id, -lag);;
+ for (size_t i = 0; i < eqns.size(); i++)
+ {
+ expr_t d = Zero;
+ for (int deriv_id : deriv_ids)
+ d = AddPlus(d, equations[eqns[i]]->getDerivative(deriv_id));
+ if (d != Zero)
+ {
+ if (!d->isConstant())
+ {
+ cerr << "ERROR: Equation '" << equation_tags.getTagValueByEqnAndKey(eqns[i], "name") << "' is not linear" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ AR[model_name][{ i, lag, lhs_symb_id }] = AddUMinus(d);
+ }
+ }
+ }
+
+ // Fill A0 matrix (for contemporaneous variables)
+ int lhs_deriv_id = getDerivID(lhs_symb_id, 0);
+ for (size_t i = 0; i < eqns.size(); i++)
+ {
+ expr_t d = equations[eqns[i]]->getDerivative(lhs_deriv_id);
+ if (d != Zero)
+ {
+ if (!d->isConstant())
+ {
+ cerr << "ERROR: Equation '" << equation_tags.getTagValueByEqnAndKey(eqns[i], "name") << "' is not linear" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ A0[model_name][{ i, lhs_symb_id }] = d;
+ }
+ }
+ }
+ }
+ var_model_table.setAR(AR);
+ var_model_table.setA0(A0);
+}
+
map<string, map<tuple<int, int, int>, expr_t>>
-DynamicModel::computeAutoregressiveMatrices(bool is_var) const
+DynamicModel::computeAutoregressiveMatrices() const
{
map<string, map<tuple<int, int, int>, expr_t>> ARr;
- const auto &all_eqnums = is_var ?
- var_model_table.getEqNums() : trend_component_model_table.getNonTargetEqNums();
- for (const auto &[model_name, eqns] : all_eqnums)
+ for (const auto &[model_name, eqns] : trend_component_model_table.getNonTargetEqNums())
{
int i = 0;
map<tuple<int, int, int>, expr_t> AR;
- const vector<int> &lhs = is_var ? var_model_table.getLhsOrigIds(model_name)
- : trend_component_model_table.getNonTargetLhs(model_name);
+ const vector<int> &lhs = trend_component_model_table.getNonTargetLhs(model_name);
for (auto eqn : eqns)
{
auto bopn = dynamic_cast<BinaryOpNode *>(equations[eqn]->arg2);
@@ -3710,7 +3818,7 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
void
DynamicModel::fillTrendComponentModelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const
{
- auto ARr = computeAutoregressiveMatrices(false);
+ auto ARr = computeAutoregressiveMatrices();
trend_component_model_table.setAR(ARr);
auto [A0r, A0starr] = computeErrorComponentMatrices(diff_subst_table);
trend_component_model_table.setA0(A0r, A0starr);
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 84fad95e..a2ad30bf 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -271,13 +271,23 @@ private:
str.erase(str.find_last_not_of("\t\n\v\f\r ") + 1);
}
- //! Compute autoregressive matrices of VAR or trend component models
- map<string, map<tuple<int, int, int>, expr_t>> computeAutoregressiveMatrices(bool is_var) const;
+ //! Compute autoregressive matrices of trend component models
+ /* The algorithm uses matching rules over expression trees. It cannot handle
+ arbitrarily-written expressions. */
+ map<string, map<tuple<int, int, int>, expr_t>> computeAutoregressiveMatrices() const;
//! Compute error component matrices of trend component_models
/*! Returns a pair (A0r, A0starr) */
pair<map<string, map<tuple<int, int, int>, expr_t>>, map<string, map<tuple<int, int, int>, expr_t>>> computeErrorComponentMatrices(const ExprNode::subst_table_t &diff_subst_table) const;
+ /* For a VAR model, given the symbol ID of a LHS variable, and a (negative)
+ lag, returns all the corresponding deriv_ids (by properly dealing with two
+ types of auxiliary variables: endo lags and diff lags). It returns a
+ vector because in some cases there may be sereval corresponding deriv_ids
+ (for example, in the deriv_id table, AUX_DIFF_nn(-1) may appear as itself
+ (with a lag), and also as a contemporaneous diff lag auxvar). */
+ vector<int> getVARDerivIDs(int lhs_symb_id, int lead_lag) const;
+
public:
DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
@@ -363,9 +373,13 @@ public:
void fillTrendComponentModelTableAREC(const ExprNode::subst_table_t &diff_subst_table) const;
//! Fill the VAR model table with information available from the transformed model
+ // NB: Does not fill the AR and A0 matrices
void fillVarModelTable() const;
//! Fill the VAR model table with information available from the original model
void fillVarModelTableFromOrigModel() const;
+ //! Fill the AR and A0 matrices of the VAR model table
+ // Uses derivatives, hence must be called after computingPass()
+ void fillVarModelTableMatrices();
//! Update the rhs references in the var model and trend component tables
//! after substitution of auxiliary variables and find the trend variables
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index cdade2eb..f4ce61ac 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -166,7 +166,7 @@ class ParsingDriver;
%token PARAMETER_CONVERGENCE_CRITERION NUMBER_OF_LARGE_PERTURBATIONS NUMBER_OF_SMALL_PERTURBATIONS
%token NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION MAX_NUMBER_OF_STAGES
%token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION NO_INIT_ESTIMATION_CHECK_FIRST_OBS
-%token HETEROSKEDASTIC_FILTER TIME_SHIFT
+%token HETEROSKEDASTIC_FILTER TIME_SHIFT STRUCTURAL
/* Method of Moments */
%token METHOD_OF_MOMENTS MOM_METHOD
%token BARTLETT_KERNEL_LAG WEIGHTING_MATRIX WEIGHTING_MATRIX_SCALING_FACTOR ANALYTIC_STANDARD_ERRORS ANALYTIC_JACOBIAN PENALIZED_ESTIMATOR VERBOSE
@@ -386,6 +386,7 @@ var_model_options_list : var_model_options_list COMMA var_model_options
var_model_options : o_var_name
| o_var_eq_tags
+ | o_var_structural
;
trend_component_model : TREND_COMPONENT_MODEL '(' trend_component_model_options_list ')' ';' { driver.trend_component_model(); }
@@ -3230,6 +3231,7 @@ o_series : SERIES EQUAL symbol { driver.option_str("series", $3); };
o_datafile : DATAFILE EQUAL filename { driver.option_str("datafile", $3); };
o_filename : FILENAME EQUAL filename { driver.option_str("filename", $3); };
o_var_eq_tags : EQTAGS EQUAL vec_str { driver.option_vec_str("var.eqtags", $3); }
+o_var_structural : STRUCTURAL { driver.option_num("var.structural", "true"); }
o_dirname : DIRNAME EQUAL filename { driver.option_str("dirname", $3); };
o_huge_number : HUGE_NUMBER EQUAL non_negative_number { driver.option_num("huge_number", $3); };
o_nobs : NOBS EQUAL vec_int
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index 4347d994..180b6919 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -815,6 +815,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
<DYNARE_STATEMENT>variances {return token::VARIANCES;}
<DYNARE_STATEMENT>equations {return token::EQUATIONS;}
<DYNARE_STATEMENT>time_shift {return token::TIME_SHIFT;}
+<DYNARE_STATEMENT>structural {return token::STRUCTURAL;}
<DYNARE_STATEMENT>true {
yylval->build<string>(yytext);
return token::TRUE;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 9982fcdd..b1f77541 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -379,6 +379,16 @@ ExprNode::matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<in
throw MatchFailureException{"Unsupported expression"};
}
+bool
+ExprNode::isConstant() const
+{
+ set<pair<int, int>> symbs_lags;
+ collectDynamicVariables(SymbolType::endogenous, symbs_lags);
+ collectDynamicVariables(SymbolType::exogenous, symbs_lags);
+ collectDynamicVariables(SymbolType::exogenousDet, symbs_lags);
+ return symbs_lags.empty();
+}
+
NumConstNode::NumConstNode(DataTree &datatree_arg, int idx_arg, int id_arg) :
ExprNode{datatree_arg, idx_arg},
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 53446c2f..766f00cf 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -703,6 +703,10 @@ public:
the first, lag in the second, exponent in the third).
Throws a MatchFailureException if not of the right form. */
virtual void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const;
+
+ /* Returns true if the expression contains no endogenous, no exogenous and no
+ exogenous deterministic */
+ bool isConstant() const;
};
//! Object used to compare two nodes (using their indexes)
diff --git a/src/ModFile.cc b/src/ModFile.cc
index bad563f3..af3e4fcb 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -837,6 +837,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 &statement : statements)
statement->computingPass(mod_file_struct);
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 77302c81..bd156cd2 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -1391,7 +1391,12 @@ ParsingDriver::var_model()
error("You must pass the eqtags option to the var_model statement.");
auto eqtags = itvs->second;
- mod_file->var_model_table.addVarModel(name, eqtags);
+ bool structural = false;
+ if (auto itn = options_list.num_options.find("var.structural");
+ itn != options_list.num_options.end() && itn->second == "true")
+ structural = true;
+
+ mod_file->var_model_table.addVarModel(name, structural, eqtags);
options_list.clear();
}
diff --git a/src/SubModel.cc b/src/SubModel.cc
index c18bc17e..ec1f2661 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -428,7 +428,7 @@ VarModelTable::VarModelTable(SymbolTable &symbol_table_arg) :
}
void
-VarModelTable::addVarModel(string name_arg, vector<string> eqtags_arg)
+VarModelTable::addVarModel(string name_arg, bool structural_arg, vector<string> eqtags_arg)
{
if (isExistingVarModelName(name_arg))
{
@@ -436,6 +436,7 @@ VarModelTable::addVarModel(string name_arg, vector<string> eqtags_arg)
exit(EXIT_FAILURE);
}
+ structural[name_arg] = structural_arg;
eqtags[name_arg] = move(eqtags_arg);
names.insert(move(name_arg));
}
@@ -454,14 +455,15 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
- ar_output << "function ar = var_ar(model_name, params)" << endl
- << "%function ar = var_ar(model_name, params)" << endl
+ ar_output << "function [ar, a0] = var_ar(model_name, params)" << endl
+ << "%function [ar, a0] = var_ar(model_name, params)" << endl
<< "% File automatically generated by the Dynare preprocessor" << endl << endl;
for (const auto &name : names)
{
- output << "M_.var." << name << ".model_name = '" << name << "';" << endl;
- output << "M_.var." << name << ".eqtags = {";
+ output << "M_.var." << name << ".model_name = '" << name << "';" << endl
+ << "M_.var." << name << ".structural = " << (structural.at(name) ? "true" : "false") << ";" << endl
+ << "M_.var." << name << ".eqtags = {";
for (const auto &it : eqtags.at(name))
output << "'" << it << "'; ";
output << "};" << endl
@@ -512,7 +514,18 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
ar_output << ";" << endl;
}
- ar_output << " return" << endl
+ ar_output << " if nargout > 1" << endl
+ << " a0 = zeros(" << lhs.size() << ", " << lhs.size() << ");" << endl;
+ for (const auto &[key, expr] : A0.at(name))
+ {
+ auto [eqn, lhs_symb_id] = key;
+ int colidx = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
+ ar_output << " a0(" << eqn + 1 << ", " << colidx + 1 << ") = ";
+ expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
+ ar_output << ";" << endl;
+ }
+ ar_output << " end" << endl
+ << " return" << endl
<< "end" << endl << endl;
}
ar_output << "error([model_name ' is not a valid var_model name'])" << endl
@@ -540,6 +553,12 @@ VarModelTable::writeJsonOutput(ostream &output) const
}
}
+const map<string, bool> &
+VarModelTable::getStructural() const
+{
+ return structural;
+}
+
const map<string, vector<string>> &
VarModelTable::getEqTags() const
{
@@ -652,6 +671,12 @@ VarModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
AR = move(AR_arg);
}
+void
+VarModelTable::setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg)
+{
+ A0 = move(A0_arg);
+}
+
const vector<int> &
VarModelTable::getMaxLags(const string &name_arg) const
{
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 627dbec6..37c35da7 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -47,6 +47,8 @@ private:
map<string, vector<expr_t>> lhs_expr_t;
map<string, vector<int>> target_vars;
map<string, map<tuple<int, int, int>, expr_t>> AR; // name -> (eqn, lag, lhs_symb_id) -> expr_t
+ /* Note that A0 in the trend-component model context is not the same thing as
+ in the structural VAR context. */
map<string, map<tuple<int, int, int>, expr_t>> A0, A0star; // name -> (eqn, lag, col) -> expr_t
public:
explicit TrendComponentModelTable(SymbolTable &symbol_table_arg);
@@ -116,21 +118,28 @@ class VarModelTable
private:
SymbolTable &symbol_table;
set<string> names;
+ map<string, bool> structural; // Whether VARs are structural or reduced-form
map<string, vector<string>> eqtags;
map<string, vector<int>> eqnums, max_lags, lhs, lhs_orig_symb_ids, orig_diff_var;
map<string, vector<set<pair<int, int>>>> rhs; // name -> for each equation: set of pairs (var, lag)
map<string, vector<bool>> diff;
map<string, vector<expr_t>> lhs_expr_t;
map<string, map<tuple<int, int, int>, expr_t>> AR; // name -> (eqn, lag, lhs_symb_id) -> param_expr_t
+ /* The A0 matrix is mainly for structural VARs. For reduced-form VARs, it
+ will be equal to the identity matrix. Also note that A0 in the structural
+ VAR context is not the same thing as in the trend-component model
+ context. */
+ map<string, map<tuple<int, int>, expr_t>> A0; // name -> (eqn, lhs_symb_id) -> param_expr_t
public:
explicit VarModelTable(SymbolTable &symbol_table_arg);
//! Add a VAR model
- void addVarModel(string name, vector<string> eqtags);
+ void addVarModel(string name, bool structural_arg, vector<string> eqtags);
inline bool isExistingVarModelName(const string &name_arg) const;
inline bool empty() const;
+ const map<string, bool> &getStructural() const;
const map<string, vector<string>> &getEqTags() const;
const vector<string> &getEqTags(const string &name_arg) const;
const map<string, vector<int>> &getEqNums() const;
@@ -151,6 +160,7 @@ public:
void setMaxLags(map<string, vector<int>> max_lags_arg);
void setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg);
void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
+ void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg);
//! Write output of this class
void writeOutput(const string &basename, ostream &output) const;
--
GitLab