diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index fa4cfb86dc3af2e01ba25bb7d84bd0a136aa2d2d..7beab03505ff96d0af8df881113eca6a250ea350 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 84fad95ed2819155ff5b5fd58a573c1a02550a9f..a2ad30bf2a3ec81053a3806bb4e905ea833a84df 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 cdade2eb6d7798fd5afcd3c74816ce5b798059fd..f4ce61ac70bc7b2f3377107870184044ef5d154e 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 4347d9947f40dc9d31201bcd1a86e55cc890bc0e..180b691973982e5b4b46c0f50c65fc184a6ea3e7 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 9982fcddc606963dbbed0f218315aa794d3dba7d..b1f77541b48f1786399b5ef6c30704a0b5ed67ba 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 53446c2f9f9f139b6c4d5654e6c59d66831563aa..766f00cfae873118b122ff38e5566f6d8f9901bd 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 bad563f31935522dcc08677bcea9fe8c84c69785..af3e4fcb9834ed8a0d49b6e3efa2a9f39302d2c4 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 77302c81a5f0adceeb48e4b3504f6a066b9d48d8..bd156cd2304d27e1f6bcd0c98092bf427e594d3a 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 c18bc17e3e7f74effa6ed9c05239ad3a960435d4..ec1f26615b7b80616e6c871e1b380078a0ba5959 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 627dbec6b829efe475c5893a1ab9f21b996c5995..37c35da73d46ad38aa83136bd16c871aac93293a 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;