/*
* Copyright © 2003-2023 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 .
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "DynamicModel.hh"
#include "ParsingDriver.hh"
void
DynamicModel::copyHelper(const DynamicModel &m)
{
auto f = [this](const ExprNode *e) { return e->clone(*this); };
for (const auto &it : m.static_only_equations)
static_only_equations.push_back(dynamic_cast(f(it)));
}
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_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},
trend_component_model_table{trend_component_model_table_arg},
var_model_table{var_model_table_arg}
{
}
DynamicModel::DynamicModel(const DynamicModel &m) :
ModelTree{m},
trend_component_model_table{m.trend_component_model_table},
var_model_table{m.var_model_table},
balanced_growth_test_tol{m.balanced_growth_test_tol},
static_only_equations_lineno{m.static_only_equations_lineno},
static_only_equations_equation_tags{m.static_only_equations_equation_tags},
deriv_id_table{m.deriv_id_table},
inv_deriv_id_table{m.inv_deriv_id_table},
dyn_jacobian_cols_table{m.dyn_jacobian_cols_table},
dyn_jacobian_ncols{m.dyn_jacobian_ncols},
max_lag{m.max_lag},
max_lead{m.max_lead},
max_endo_lag{m.max_endo_lag},
max_endo_lead{m.max_endo_lead},
max_exo_lag{m.max_exo_lag},
max_exo_lead{m.max_exo_lead},
max_exo_det_lag{m.max_exo_det_lag},
max_exo_det_lead{m.max_exo_det_lead},
max_lag_orig{m.max_lag_orig},
max_lead_orig{m.max_lead_orig},
max_lag_with_diffs_expanded_orig{m.max_lag_with_diffs_expanded_orig},
max_endo_lag_orig{m.max_endo_lag_orig},
max_endo_lead_orig{m.max_endo_lead_orig},
max_exo_lag_orig{m.max_exo_lag_orig},
max_exo_lead_orig{m.max_exo_lead_orig},
max_exo_det_lag_orig{m.max_exo_det_lag_orig},
max_exo_det_lead_orig{m.max_exo_det_lead_orig},
xrefs{m.xrefs},
xref_param{m.xref_param},
xref_endo{m.xref_endo},
xref_exo{m.xref_exo},
xref_exo_det{m.xref_exo_det},
nonzero_hessian_eqs{m.nonzero_hessian_eqs},
variableMapping{m.variableMapping},
blocks_jacob_cols_endo{m.blocks_jacob_cols_endo},
var_expectation_functions_to_write{m.var_expectation_functions_to_write},
mfs{m.mfs},
static_mfs{m.static_mfs}
{
copyHelper(m);
}
DynamicModel &
DynamicModel::operator=(const DynamicModel &m)
{
ModelTree::operator=(m);
assert(&trend_component_model_table == &m.trend_component_model_table);
assert(&var_model_table == &m.var_model_table);
balanced_growth_test_tol = m.balanced_growth_test_tol;
static_only_equations_lineno = m.static_only_equations_lineno;
static_only_equations_equation_tags = m.static_only_equations_equation_tags;
deriv_id_table = m.deriv_id_table;
inv_deriv_id_table = m.inv_deriv_id_table;
dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
dyn_jacobian_ncols = m.dyn_jacobian_ncols;
max_lag = m.max_lag;
max_lead = m.max_lead;
max_endo_lag = m.max_endo_lag;
max_endo_lead = m.max_endo_lead;
max_exo_lag = m.max_exo_lag;
max_exo_lead = m.max_exo_lead;
max_exo_det_lag = m.max_exo_det_lag;
max_exo_det_lead = m.max_exo_det_lead;
max_lag_orig = m.max_lag_orig;
max_lead_orig = m.max_lead_orig;
max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
max_endo_lag_orig = m.max_endo_lag_orig;
max_endo_lead_orig = m.max_endo_lead_orig;
max_exo_lag_orig = m.max_exo_lag_orig;
max_exo_lead_orig = m.max_exo_lead_orig;
max_exo_det_lag_orig = m.max_exo_det_lag_orig;
max_exo_det_lead_orig = m.max_exo_det_lead_orig;
xrefs = m.xrefs;
xref_param = m.xref_param;
xref_endo = m.xref_endo;
xref_exo = m.xref_exo;
xref_exo_det = m.xref_exo_det;
nonzero_hessian_eqs = m.nonzero_hessian_eqs;
variableMapping = m.variableMapping;
blocks_jacob_cols_endo = m.blocks_jacob_cols_endo;
var_expectation_functions_to_write = m.var_expectation_functions_to_write;
mfs = m.mfs;
static_mfs = m.static_mfs;
copyHelper(m);
return *this;
}
void
DynamicModel::writeDynamicBytecode(const string &basename) const
{
/* Bytecode only works when there are with as many endogenous as equations.
(e.g. the constructor of FBEGINBLOCK_ makes this assumption) */
assert(static_cast(equations.size()) == symbol_table.endo_nbr());
// Determine the type of model (used for typing the single block)
BlockSimulationType simulation_type;
if (max_endo_lag > 0 && max_endo_lead > 0)
simulation_type = BlockSimulationType::solveTwoBoundariesComplete;
else if (max_endo_lag >= 0 && max_endo_lead == 0)
simulation_type = BlockSimulationType::solveForwardComplete;
else
simulation_type = BlockSimulationType::solveBackwardComplete;
// First write the .bin file
int u_count_int { writeBytecodeBinFile(basename + "/model/bytecode/dynamic.bin",
simulation_type == BlockSimulationType::solveTwoBoundariesComplete) };
BytecodeWriter code_file {basename + "/model/bytecode/dynamic.cod"};
// Declare temporary terms
code_file << FDIMT_{static_cast(temporary_terms_derivatives[0].size()
+ temporary_terms_derivatives[1].size())};
// Declare the (single) block
vector exo(symbol_table.exo_nbr()), exo_det(symbol_table.exo_det_nbr());
iota(exo.begin(), exo.end(), 0);
iota(exo_det.begin(), exo_det.end(), 0);
int jacobian_ncols_endo
{ static_cast(count_if(dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
[this](const auto &v)
{ return getTypeByDerivID(v.first) == SymbolType::endogenous; }))
};
vector eq_idx(equations.size());
iota(eq_idx.begin(), eq_idx.end(), 0);
vector endo_idx(symbol_table.endo_nbr());
iota(endo_idx.begin(), endo_idx.end(), 0);
code_file << FBEGINBLOCK_{symbol_table.endo_nbr(),
simulation_type,
0,
symbol_table.endo_nbr(),
endo_idx,
eq_idx,
false,
u_count_int,
jacobian_ncols_endo,
symbol_table.exo_det_nbr(),
symbol_table.exo_nbr(),
exo_det,
exo};
writeBytecodeHelper(code_file);
}
void
DynamicModel::writeDynamicBlockBytecode(const string &basename) const
{
BytecodeWriter code_file {basename + "/model/bytecode/block/dynamic.cod"};
const filesystem::path bin_filename {basename + "/model/bytecode/block/dynamic.bin"};
ofstream bin_file {bin_filename, ios::out | ios::binary};
if (!bin_file.is_open())
{
cerr << R"(Error : Can't open file ")" << bin_filename.string() << R"(" for writing)" << endl;
exit(EXIT_FAILURE);
}
// Temporary variables declaration
code_file << FDIMT_{static_cast(blocks_temporary_terms_idxs.size())};
temporary_terms_t temporary_terms_written;
for (int block {0}; block < static_cast(blocks.size()); block++)
{
const BlockSimulationType simulation_type {blocks[block].simulation_type};
// Write section of .bin file except for evaluate blocks and solve simple blocks
const int u_count {simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete
? writeBlockBytecodeBinFile(bin_file, block)
: 0};
code_file << FBEGINBLOCK_{blocks[block].mfs_size,
simulation_type,
blocks[block].first_equation,
blocks[block].size,
endo_idx_block2orig,
eq_idx_block2orig,
blocks[block].linear,
u_count,
static_cast(blocks_jacob_cols_endo[block].size())};
writeBlockBytecodeHelper(code_file, block, temporary_terms_written);
}
code_file << FEND_{};
}
void
DynamicModel::writeDynamicMFile(const string &basename) const
{
auto [d_output, tt_output] = writeModelFileHelper();
ostringstream init_output, end_output;
init_output << "residual = zeros(" << equations.size() << ", 1);";
writeDynamicMFileHelper(basename, "dynamic_resid", "residual", "dynamic_resid_tt",
temporary_terms_derivatives[0].size(),
"", init_output, end_output, d_output[0], tt_output[0]);
init_output.str("");
init_output << "g1 = zeros(" << equations.size() << ", " << getJacobianColsNbr(false) << ");";
writeDynamicMFileHelper(basename, "dynamic_g1", "g1", "dynamic_g1_tt",
temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
"dynamic_resid_tt", init_output, end_output, d_output[1], tt_output[1]);
writeDynamicMWrapperFunction(basename, "g1");
// For order ≥ 2
int ncols{getJacobianColsNbr(false)};
int ntt { static_cast(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())};
for (size_t i{2}; i < derivatives.size(); i++)
{
ncols *= getJacobianColsNbr(false);
ntt += temporary_terms_derivatives[i].size();
string gname{"g" + to_string(i)};
string gprevname{"g" + to_string(i-1)};
init_output.str("");
end_output.str("");
if (derivatives[i].size())
{
init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
<< gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
<< gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
end_output << gname << " = sparse("
<< gname << "_i," << gname << "_j," << gname << "_v,"
<< equations.size() << "," << ncols << ");";
}
else
init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");";
writeDynamicMFileHelper(basename, "dynamic_" + gname, gname, "dynamic_" + gname + "_tt", ntt,
"dynamic_" + gprevname + "_tt", init_output, end_output,
d_output[i], tt_output[i]);
if (i <= 3)
writeDynamicMWrapperFunction(basename, gname);
}
writeDynamicMCompatFile(basename);
}
string
DynamicModel::reform(const string &name1) const
{
string name = name1;
int pos = name.find(R"(\)", 0);
while (pos >= 0)
{
if (name.substr(pos + 1, 1) != R"(\)")
{
name = name.insert(pos, R"(\)");
pos++;
}
pos++;
pos = name.find(R"(\)", pos);
}
return name;
}
void
DynamicModel::printNonZeroHessianEquations(ostream &output) const
{
if (nonzero_hessian_eqs.size() != 1)
output << "[";
for (bool printed_something{false};
int it : nonzero_hessian_eqs)
{
if (exchange(printed_something, true))
output << " ";
output << it + 1;
}
if (nonzero_hessian_eqs.size() != 1)
output << "]";
}
void
DynamicModel::writeDynamicMWrapperFunction(const string &basename, const string &ending) const
{
string name;
if (ending == "g1")
name = "dynamic_resid_g1";
else if (ending == "g2")
name = "dynamic_resid_g1_g2";
else if (ending == "g3")
name = "dynamic_resid_g1_g2_g3";
filesystem::path filename {packageDir(basename) / (name + ".m")};
ofstream output{filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
if (ending == "g1")
output << "function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl;
else if (ending == "g2")
output << "function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl;
else if (ending == "g3")
output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl;
output << "%" << endl
<< "% Wrapper function automatically created by Dynare" << endl
<< "%" << endl
<< endl
<< " if T_flag" << endl
<< " T = " << basename << ".dynamic_" << ending << "_tt(T, y, x, params, steady_state, it_);" << endl
<< " end" << endl;
if (ending == "g1")
output << " residual = " << basename << ".dynamic_resid(T, y, x, params, steady_state, it_, false);" << endl
<< " g1 = " << basename << ".dynamic_g1(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g2")
output << " [residual, g1] = " << basename << ".dynamic_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
<< " g2 = " << basename << ".dynamic_g2(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g3")
output << " [residual, g1, g2] = " << basename << ".dynamic_resid_g1_g2(T, y, x, params, steady_state, it_, false);" << endl
<< " g3 = " << basename << ".dynamic_g3(T, y, x, params, steady_state, it_, false);" << endl;
output << endl << "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicMFileHelper(const string &basename,
const string &name, const string &retvalname,
const string &name_tt, size_t ttlen,
const string &previous_tt_name,
const ostringstream &init_s, const ostringstream &end_s,
const ostringstream &s, const ostringstream &s_tt) const
{
filesystem::path filename {packageDir(basename) / (name_tt + ".m")};
ofstream output{filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "% function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables in the order stored" << endl
<< "% in M_.lead_lag_incidence; see the Manual" << endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables (in declaration order)" << endl
<< "% for all simulation periods" << endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values" << endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
<< "% it_ scalar double time period for exogenous variables for which" << endl
<< "% to evaluate the model" << endl
<< "%" << endl
<< "% Output:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms" << endl
<< "%" << endl << endl
<< "assert(length(T) >= " << ttlen << ");" << endl
<< endl;
if (!previous_tt_name.empty())
output << "T = " << basename << "." << previous_tt_name << "(T, y, x, params, steady_state, it_);" << endl << endl;
output << s_tt.str() << endl
<< "end" << endl;
output.close();
filename = packageDir(basename) / (name + ".m");
output.open(filename, ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables in the order stored" << endl
<< "% in M_.lead_lag_incidence; see the Manual" << endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables (in declaration order)" << endl
<< "% for all simulation periods" << endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values" << endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
<< "% it_ scalar double time period for exogenous variables for which" << endl
<< "% to evaluate the model" << endl
<< "% T_flag boolean boolean flag saying whether or not to calculate temporary terms" << endl
<< "%" << endl
<< "% Output:" << endl
<< "% " << retvalname << endl
<< "%" << endl << endl;
if (!name_tt.empty())
output << "if T_flag" << endl
<< " T = " << basename << "." << name_tt << "(T, y, x, params, steady_state, it_);" << endl
<< "end" << endl;
output << init_s.str() << endl
<< s.str()
<< end_s.str() << endl
<< "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicMCompatFile(const string &basename) const
{
filesystem::path filename {packageDir(basename) / "dynamic.m"};
ofstream output{filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
int ntt { static_cast(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size()) };
output << "function [residual, g1, g2, g3] = dynamic(y, x, params, steady_state, it_)" << endl
<< " T = NaN(" << ntt << ", 1);" << endl
<< " if nargout <= 1" << endl
<< " residual = " << basename << ".dynamic_resid(T, y, x, params, steady_state, it_, true);" << endl
<< " elseif nargout == 2" << endl
<< " [residual, g1] = " << basename << ".dynamic_resid_g1(T, y, x, params, steady_state, it_, true);" << endl
<< " elseif nargout == 3" << endl
<< " [residual, g1, g2] = " << basename << ".dynamic_resid_g1_g2(T, y, x, params, steady_state, it_, true);" << endl
<< " else" << endl
<< " [residual, g1, g2, g3] = " << basename << ".dynamic_resid_g1_g2_g3(T, y, x, params, steady_state, it_, true);" << endl
<< " end" << endl
<< "end" << endl;
output.close();
}
vector>
DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs)
{
auto removeLeadingTrailingWhitespace = [](string &str)
{
str.erase(0, str.find_first_not_of("\t\n\v\f\r "));
str.erase(str.find_last_not_of("\t\n\v\f\r ") + 1);
};
string tags;
if (filesystem::exists(inc_exc_option_value))
{
ifstream exclude_file;
exclude_file.open(inc_exc_option_value, ifstream::in);
if (!exclude_file.is_open())
{
cerr << "ERROR: Could not open " << inc_exc_option_value << endl;
exit(EXIT_FAILURE);
}
string line;
bool tagname_on_first_line = false;
while (getline(exclude_file, line))
{
removeLeadingTrailingWhitespace(line);
if (!line.empty())
{
if (tags.empty() && line.find("=") != string::npos)
{
tagname_on_first_line = true;
tags += line + "(";
}
else if (line.find("'") != string::npos)
tags += line + ",";
else
tags += "'" + line + "',";
}
}
if (!tags.empty())
{
tags = tags.substr(0, tags.size()-1);
if (tagname_on_first_line)
tags += ")";
}
}
else
tags = inc_exc_option_value;
removeLeadingTrailingWhitespace(tags);
if (tags.front() == '[' && tags.back() != ']')
{
cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
<< ": if the first character is '[' the last must be ']'" << endl;
exit(EXIT_FAILURE);
}
if (tags.front() == '[' && tags.back() == ']')
tags = tags.substr(1, tags.length() - 2);
removeLeadingTrailingWhitespace(tags);
regex q(R"(^\w+\s*=)");
smatch matches;
string tagname = "name";
if (regex_search(tags, matches, q))
{
tagname = matches[0].str();
tags = tags.substr(tagname.size(), tags.length() - tagname.size() + 1);
removeLeadingTrailingWhitespace(tags);
if (tags.front() == '(' && tags.back() == ')')
{
tags = tags.substr(1, tags.length() - 2);
removeLeadingTrailingWhitespace(tags);
}
tagname = tagname.substr(0, tagname.size()-1);
removeLeadingTrailingWhitespace(tagname);
}
string quote_regex = "'[^']+'";
string non_quote_regex = R"([^,\s]+)";
regex r(R"((\s*)" + quote_regex + "|" + non_quote_regex + R"(\s*)(,\s*()" + quote_regex + "|" + non_quote_regex + R"()\s*)*)");
if (!regex_match(tags, r))
{
cerr << "ERROR: " << (exclude_eqs ? "exclude_eqs" : "include_eqs")
<< ": argument is of incorrect format." << endl;
exit(EXIT_FAILURE);
}
vector> eq_tag_set;
regex s(quote_regex + "|" + non_quote_regex);
for (auto it = sregex_iterator(tags.begin(), tags.end(), s);
it != sregex_iterator(); ++it)
{
string_view str {it->str()};
if (str.front() == '\'' && str.back() == '\'')
{
str.remove_prefix(1);
str.remove_suffix(1);
}
eq_tag_set.emplace_back(tagname, str);
}
return eq_tag_set;
}
vector
DynamicModel::removeEquationsHelper(set> &listed_eqs_by_tag, bool exclude_eqs,
bool excluded_vars_change_type,
vector &all_equations,
vector> &all_equations_lineno,
EquationTags &all_equation_tags, bool static_equations) const
{
if (all_equations.empty())
return {};
/* Try to convert the list of equations by tags into a list of equation
numbers.
The tag pairs that match an equation are removed from the list, so that
the caller knows which tag pairs have not been handled. */
set listed_eqs_by_number;
for (auto it = listed_eqs_by_tag.begin(); it != listed_eqs_by_tag.end();)
if (auto tmp = all_equation_tags.getEqnsByTag(it->first, it->second);
!tmp.empty())
{
listed_eqs_by_number.insert(tmp.begin(), tmp.end());
it = listed_eqs_by_tag.erase(it);
}
else
++it;
// Compute the indices of equations to be actually deleted
set eqs_to_delete_by_number;
if (exclude_eqs)
eqs_to_delete_by_number = listed_eqs_by_number;
else
for (size_t i = 0; i < all_equations.size(); i++)
if (!listed_eqs_by_number.contains(i))
eqs_to_delete_by_number.insert(i);
// remove from equations, equations_lineno, equation_tags
vector new_equations;
vector> new_equations_lineno;
map old_eqn_num_2_new;
vector excluded_vars;
for (size_t i = 0; i < all_equations.size(); i++)
if (eqs_to_delete_by_number.contains(i))
{
if (excluded_vars_change_type)
{
if (auto tmp = all_equation_tags.getTagValueByEqnAndKey(i, "endogenous"); tmp)
excluded_vars.push_back(symbol_table.getID(*tmp));
else
{
set result;
all_equations[i]->arg1->collectVariables(SymbolType::endogenous, result);
if (result.size() == 1)
excluded_vars.push_back(*result.begin());
else
{
cerr << "ERROR: Equation " << i+1
<< " has been excluded but it does not have a single variable on its left-hand side or an `endogenous` tag" << endl;
exit(EXIT_FAILURE);
}
}
}
}
else
{
new_equations.emplace_back(all_equations[i]);
old_eqn_num_2_new[i] = new_equations.size() - 1;
new_equations_lineno.emplace_back(all_equations_lineno[i]);
}
int n_excl = all_equations.size() - new_equations.size();
all_equations = new_equations;
all_equations_lineno = new_equations_lineno;
all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
if (!static_equations)
for (size_t i = 0; i < excluded_vars.size(); i++)
for (size_t j = i+1; j < excluded_vars.size(); j++)
if (excluded_vars[i] == excluded_vars[j])
{
cerr << "ERROR: Variable " << symbol_table.getName(i) << " was excluded twice"
<< " via a model_remove or model_replace statement, or via the include_eqs or exclude_eqs option" << endl;
exit(EXIT_FAILURE);
}
cout << "Excluded " << n_excl << (static_equations ? " static " : " dynamic ")
<< "equation" << (n_excl > 1 ? "s" : "") << " via model_remove or model_replace statement, or via include_eqs or exclude_eqs option" << endl;
return excluded_vars;
}
void
DynamicModel::removeEquations(const vector> &listed_eqs_by_tag, bool exclude_eqs,
bool excluded_vars_change_type)
{
/* Convert the const vector to a (mutable) set */
set listed_eqs_by_tag2(listed_eqs_by_tag.begin(), listed_eqs_by_tag.end());
vector excluded_vars = removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs,
excluded_vars_change_type,
equations, equations_lineno,
equation_tags, false);
// Ignore output because variables are not excluded when equations marked 'static' are excluded
removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
static_only_equations, static_only_equations_lineno,
static_only_equations_equation_tags, true);
if (!listed_eqs_by_tag2.empty())
{
cerr << "ERROR: model_remove/model_replace/exclude_eqs/include_eqs: The equations specified by" << endl;
for (const auto &[tagname, tagvalue] : listed_eqs_by_tag)
cerr << " " << tagname << "=" << tagvalue << endl;
cerr << "were not found." << endl;
exit(EXIT_FAILURE);
}
if (excluded_vars_change_type)
{
// Collect list of used variables in updated list of equations
set eqn_vars;
for (auto eqn : equations)
eqn->collectVariables(SymbolType::endogenous, eqn_vars);
for (auto eqn : static_only_equations)
eqn->collectVariables(SymbolType::endogenous, eqn_vars);
/* Change type of endogenous variables determined by excluded equations.
They become exogenous if they are still used somewhere, otherwise they are
completely excluded from the model. */
for (auto ev : excluded_vars)
if (eqn_vars.contains(ev))
{
symbol_table.changeType(ev, SymbolType::exogenous);
cerr << "Variable '" << symbol_table.getName(ev) << "' turned into an exogenous, as its defining equation has been removed (but it still appears in an equation)" << endl;
}
else
{
symbol_table.changeType(ev, SymbolType::excludedVariable);
cerr << "Variable '" << symbol_table.getName(ev) << "' has been excluded from the model, as its defining equation has been removed and it appears nowhere else" << endl;
}
}
}
void
DynamicModel::includeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs)
{
if (inc_exc_option_value.empty())
return;
auto listed_eqs_by_tag = parseIncludeExcludeEquations(inc_exc_option_value, exclude_eqs);
removeEquations(listed_eqs_by_tag, exclude_eqs, true);
/* There is already a check about #static and #dynamic in
ModFile::checkPass(), but the present method is called from
ModFile::transformPass(), so we must do the check again */
if (staticOnlyEquationsNbr() != dynamicOnlyEquationsNbr())
{
cerr << "ERROR: exclude_eqs/include_eqs: You must remove the same number of equations marked `static` as equations marked `dynamic`." << endl;
exit(EXIT_FAILURE);
}
}
void
DynamicModel::writeBlockDriverOutput(ostream &output) const
{
output << "M_.block_structure.time_recursive = " << boolalpha << time_recursive_block_decomposition << ";" << endl;
for (int blk = 0; blk < static_cast(blocks.size()); blk++)
{
int block_size = blocks[blk].size;
output << "M_.block_structure.block(" << blk+1 << ").Simulation_Type = " << static_cast(blocks[blk].simulation_type) << ";" << endl
<< "M_.block_structure.block(" << blk+1 << ").endo_nbr = " << block_size << ";" << endl
<< "M_.block_structure.block(" << blk+1 << ").mfs = " << blocks[blk].mfs_size << ";" << endl
<< "M_.block_structure.block(" << blk+1 << ").equation = [";
for (int eq = 0; eq < block_size; eq++)
output << " " << getBlockEquationID(blk, eq)+1;
output << "];" << endl
<< "M_.block_structure.block(" << blk+1 << ").variable = [";
for (int var = 0; var < block_size; var++)
output << " " << getBlockVariableID(blk, var)+1;
output << "];" << endl
<< "M_.block_structure.block(" << blk+1 << ").is_linear = " << boolalpha << blocks[blk].linear << ';' << endl
<< "M_.block_structure.block(" << blk+1 << ").NNZDerivatives = " << blocks_derivatives[blk].size() << ';' << endl
<< "M_.block_structure.block(" << blk+1 << ").bytecode_jacob_cols_to_sparse = [";
const bool one_boundary {blocks[blk].simulation_type == BlockSimulationType::solveBackwardSimple
|| blocks[blk].simulation_type == BlockSimulationType::solveForwardSimple
|| blocks[blk].simulation_type == BlockSimulationType::solveBackwardComplete
|| blocks[blk].simulation_type == BlockSimulationType::solveForwardComplete};
/* bytecode_jacob_cols_to_sparse is an array that maps column indices in
the Jacobian returned by bytecode MEX in evaluate mode (which only
contains nonzero columns, but also incorporate recursive variables),
into column indices in the sparse representiation (which has 3×n
columns for two-boundaries blocks and n columns for one-boundary
blocks). Columns unused in the sparse representation are indicated by
a zero. */
vector bytecode_jacob_cols_to_sparse(blocks_jacob_cols_endo[blk].size());
for (auto &[key, index] : blocks_jacob_cols_endo[blk])
{
auto [var, lag] {key};
if (var >= blocks[blk].getRecursiveSize() // NB: this check can be removed once Jacobian no longer contains columns for derivatives w.r.t. recursive variables
&& !(one_boundary && lag != 0))
bytecode_jacob_cols_to_sparse[index] = static_cast(!one_boundary)*(lag+1)*blocks[blk].mfs_size + var - blocks[blk].getRecursiveSize();
else
bytecode_jacob_cols_to_sparse[index] = -1;
}
for (int i : bytecode_jacob_cols_to_sparse)
output << i+1 << " ";
output << "];" << endl;
}
writeBlockDriverSparseIndicesHelper(output);
output << "M_.block_structure.variable_reordered = [";
for (int i = 0; i < symbol_table.endo_nbr(); i++)
output << " " << endo_idx_block2orig[i]+1;
output << "];" << endl
<< "M_.block_structure.equation_reordered = [";
for (int i = 0; i < symbol_table.endo_nbr(); i++)
output << " " << eq_idx_block2orig[i]+1;
output << "];" << endl;
map>> lag_row_incidence;
for (const auto &[indices, d1] : derivatives[1])
if (int deriv_id = indices[1];
getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
int eq = indices[0];
int var { getTypeSpecificIDByDerivID(deriv_id) };
int lag = getLagByDerivID(deriv_id);
lag_row_incidence[lag].insert({ eq, var });
}
for (auto [lag, eq_var_set] : lag_row_incidence)
{
output << "M_.block_structure.incidence(" << max_endo_lag+lag+1 << ").lead_lag = " << lag << ";" << endl
<< "M_.block_structure.incidence(" << max_endo_lag+lag+1 << ").sparse_IM = [" << endl;
for (auto [eq, var] : eq_var_set)
output << " " << eq+1 << " " << var+1 << ";" << endl;
output << "];" << endl;
}
output << "M_.block_structure.dyn_tmp_nbr = " << blocks_temporary_terms_idxs.size() << ';' << endl;
}
void
DynamicModel::writeDriverOutput(ostream &output, bool compute_xrefs) const
{
/* Writing initialisation for M_.lead_lag_incidence matrix
M_.lead_lag_incidence is a matrix with as many columns as there are
endogenous variables and as many rows as there are periods in the
models (nbr of rows = M_.max_lag+M_.max_lead+1)
The matrix elements are equal to zero if a variable isn't present in the
model at a given period.
*/
output << "M_.orig_maximum_endo_lag = " << max_endo_lag_orig << ";" << endl
<< "M_.orig_maximum_endo_lead = " << max_endo_lead_orig << ";" << endl
<< "M_.orig_maximum_exo_lag = " << max_exo_lag_orig << ";" << endl
<< "M_.orig_maximum_exo_lead = " << max_exo_lead_orig << ";" << endl
<< "M_.orig_maximum_exo_det_lag = " << max_exo_det_lag_orig << ";" << endl
<< "M_.orig_maximum_exo_det_lead = " << max_exo_det_lead_orig << ";" << endl
<< "M_.orig_maximum_lag = " << max_lag_orig << ";" << endl
<< "M_.orig_maximum_lead = " << max_lead_orig << ";" << endl
<< "M_.orig_maximum_lag_with_diffs_expanded = " << max_lag_with_diffs_expanded_orig << ";" << endl
<< "M_.lead_lag_incidence = [";
// Loop on endogenous variables
int nstatic = 0,
nfwrd = 0,
npred = 0,
nboth = 0;
for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
{
output << endl;
int sstatic = 1,
sfwrd = 0,
spred = 0,
sboth = 0;
// Loop on periods
for (int lag = -max_endo_lag; lag <= max_endo_lead; lag++)
{
// Print variableID if exists with current period, otherwise print 0
try
{
int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag);
output << " " << getJacobianCol(varID, false) + 1;
if (lag == -1)
{
sstatic = 0;
spred = 1;
}
else if (lag == 1)
{
if (spred == 1)
{
sboth = 1;
spred = 0;
}
else
{
sstatic = 0;
sfwrd = 1;
}
}
}
catch (UnknownDerivIDException &e)
{
output << " 0";
}
}
nstatic += sstatic;
nfwrd += sfwrd;
npred += spred;
nboth += sboth;
output << ";";
}
output << "]';" << endl
<< "M_.nstatic = " << nstatic << ";" << endl
<< "M_.nfwrd = " << nfwrd << ";" << endl
<< "M_.npred = " << npred << ";" << endl
<< "M_.nboth = " << nboth << ";" << endl
<< "M_.nsfwrd = " << nfwrd+nboth << ";" << endl
<< "M_.nspred = " << npred+nboth << ";" << endl
<< "M_.ndynamic = " << npred+nboth+nfwrd << ";" << endl
<< "M_.dynamic_tmp_nbr = [";
for (const auto &tts : temporary_terms_derivatives)
output << tts.size() << "; ";
output << "];" << endl;
// Write equation tags
equation_tags.writeOutput(output);
// Write mapping for variables and equations they are present in
for (const auto &variable : variableMapping)
{
output << "M_.mapping." << symbol_table.getName(variable.first) << ".eqidx = [";
for (auto equation : variable.second)
output << equation + 1 << " ";
output << "];" << endl;
}
/* Say if static and dynamic models differ (because of [static] and [dynamic]
equation tags) */
output << "M_.static_and_dynamic_models_differ = "
<< boolalpha << (static_only_equations.size() > 0)
<< ";" << endl;
// Say if model contains an external function call
bool has_external_function = false;
for (auto equation : equations)
if (equation->containsExternalFunction())
{
has_external_function = true;
break;
}
output << "M_.has_external_function = " << boolalpha << has_external_function
<< ';' << endl;
// Compute list of state variables, ordered in block-order
vector state_var;
for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
// Loop on negative lags
for (int lag = -max_endo_lag; lag < 0; lag++)
try
{
getDerivID(symbol_table.getID(SymbolType::endogenous, endo_idx_block2orig[endoID]), lag);
if (find(state_var.begin(), state_var.end(), endo_idx_block2orig[endoID]) == state_var.end())
state_var.push_back(endo_idx_block2orig[endoID]);
}
catch (UnknownDerivIDException &e)
{
}
// Write the block structure of the model
if (block_decomposed)
writeBlockDriverOutput(output);
output << "M_.state_var = [";
for (int it : state_var)
output << it+1 << " ";
output << "];" << endl;
// Writing initialization for some other variables
output << "M_.maximum_lag = " << max_lag << ";" << endl
<< "M_.maximum_lead = " << max_lead << ";" << endl;
output << "M_.maximum_endo_lag = " << max_endo_lag << ";" << endl
<< "M_.maximum_endo_lead = " << max_endo_lead << ";" << endl
<< "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl;
output << "M_.maximum_exo_lag = " << max_exo_lag << ";" << endl
<< "M_.maximum_exo_lead = " << max_exo_lead << ";" << endl
<< "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl;
if (symbol_table.exo_det_nbr())
{
output << "M_.maximum_exo_det_lag = " << max_exo_det_lag << ";" << endl
<< "M_.maximum_exo_det_lead = " << max_exo_det_lead << ";" << endl
<< "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);" << endl;
}
output << "M_.params = " << "NaN(" << symbol_table.param_nbr() << ", 1);" << endl;
string empty_cell = "cell(" + to_string(symbol_table.endo_nbr()) + ", 1)";
output << "M_.endo_trends = struct('deflator', " << empty_cell
<< ", 'log_deflator', " << empty_cell << ", 'growth_factor', " << empty_cell
<< ", 'log_growth_factor', " << empty_cell << ");" << endl;
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
int symb_id = symbol_table.getID(SymbolType::endogenous, i);
if (auto it = nonstationary_symbols_map.find(symb_id); it != nonstationary_symbols_map.end())
{
auto [is_log, deflator] = it->second;
output << "M_.endo_trends(" << i+1 << ")."
<< (is_log ? "log_deflator" : "deflator") << " = '";
deflator->writeJsonOutput(output, {}, {});
output << "';" << endl;
auto growth_factor = const_cast(this)->AddDivide(deflator, deflator->decreaseLeadsLags(1))->removeTrendLeadLag(trend_symbols_map)->replaceTrendVar();
output << "M_.endo_trends(" << i+1 << ")."
<< (is_log ? "log_growth_factor" : "growth_factor") << " = '";
growth_factor->writeJsonOutput(output, {}, {});
output << "';" << endl;
}
}
if (compute_xrefs)
writeXrefs(output);
// Write number of non-zero derivatives
// Use -1 if the derivatives have not been computed
output << "M_.NNZDerivatives = [";
for (int i = 1; i < static_cast(NNZDerivatives.size()); i++)
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
output << "];" << endl;
writeDriverSparseIndicesHelper(output);
// Write LHS of each equation in text form
output << "M_.lhs = {" << endl;
for (auto eq : equations)
{
output << "'";
eq->arg1->writeJsonOutput(output, {}, {});
output << "'; " << endl;
}
output << "};" << endl;
}
void
DynamicModel::runTrendTest(const eval_context_t &eval_context)
{
computeDerivIDs();
testTrendDerivativesEqualToZero(eval_context);
}
void
DynamicModel::updateVarAndTrendModel() const
{
for (bool var : { true, false })
{
map>> trend_varr;
map>>> rhsr;
for (const auto &[model_name, eqns] : (var ? var_model_table.getEqNums()
: trend_component_model_table.getEqNums()))
{
vector lhs, trend_lhs;
vector> trend_var;
vector>> rhs;
if (!var)
{
lhs = trend_component_model_table.getLhs(model_name);
for (auto teqn : trend_component_model_table.getTargetEqNums().at(model_name))
{
int eqnidx = 0;
for (auto eqn : eqns)
{
if (eqn == teqn)
trend_lhs.push_back(lhs[eqnidx]);
eqnidx++;
}
}
}
int lhs_idx = 0;
for (auto eqn : eqns)
{
set> rhs_set;
equations[eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_set);
rhs.push_back(rhs_set);
if (!var)
{
int lhs_symb_id = lhs[lhs_idx++];
if (symbol_table.isDiffAuxiliaryVariable(lhs_symb_id))
try
{
lhs_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_symb_id);
}
catch (...)
{
}
optional trend_var_symb_id = equations[eqn]->arg2->findTargetVariable(lhs_symb_id);
if (trend_var_symb_id)
{
if (symbol_table.isDiffAuxiliaryVariable(*trend_var_symb_id))
try
{
trend_var_symb_id = symbol_table.getOrigSymbIdForAuxVar(*trend_var_symb_id);
}
catch (...)
{
}
if (find(trend_lhs.begin(), trend_lhs.end(), *trend_var_symb_id) == trend_lhs.end())
{
cerr << "ERROR: trend found in trend_component equation #" << eqn << " ("
<< symbol_table.getName(*trend_var_symb_id) << ") does not correspond to a trend equation" << endl;
exit(EXIT_FAILURE);
}
}
trend_var.push_back(move(trend_var_symb_id));
}
}
rhsr[model_name] = rhs;
if (!var)
trend_varr[model_name] = trend_var;
}
if (var)
var_model_table.setRhs(move(rhsr));
else
{
trend_component_model_table.setRhs(move(rhsr));
trend_component_model_table.setTargetVar(move(trend_varr));
}
}
}
void
DynamicModel::fillVarModelTable() const
{
map> eqnums, lhsr;
map> lhs_expr_tr;
map>>> rhsr;
for (const auto &[model_name, eqtags] : var_model_table.getEqTags())
{
vector eqnumber, lhs;
vector lhs_expr_t;
vector>> rhs;
for (const auto &eqtag : eqtags)
{
set> lhs_set, lhs_tmp_set, rhs_set;
optional eqn { equation_tags.getEqnByTag("name", eqtag) };
if (!eqn)
{
cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
exit(EXIT_FAILURE);
}
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::endogenous, lhs_set);
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::exogenous, lhs_tmp_set);
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::parameter, lhs_tmp_set);
if (lhs_set.size() != 1 || !lhs_tmp_set.empty())
{
cerr << "ERROR: in Equation " << eqtag
<< ". A VAR may only have one endogenous variable on the LHS. " << endl;
exit(EXIT_FAILURE);
}
auto itlhs = lhs_set.begin();
if (itlhs->second != 0)
{
cerr << "ERROR: in Equation " << eqtag
<< ". The variable on the LHS of a VAR may not appear with a lead or a lag. "
<< endl;
exit(EXIT_FAILURE);
}
eqnumber.push_back(*eqn);
lhs.push_back(itlhs->first);
lhs_set.clear();
set lhs_expr_t_set;
equations[*eqn]->arg1->collectVARLHSVariable(lhs_expr_t_set);
lhs_expr_t.push_back(*(lhs_expr_t_set.begin()));
equations[*eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_set);
rhs.push_back(rhs_set);
}
eqnums[model_name] = eqnumber;
lhsr[model_name] = lhs;
lhs_expr_tr[model_name] = lhs_expr_t;
rhsr[model_name] = rhs;
}
var_model_table.setEqNums(move(eqnums));
var_model_table.setLhs(move(lhsr));
var_model_table.setRhs(move(rhsr));
var_model_table.setLhsExprT(move(lhs_expr_tr));
}
void
DynamicModel::fillVarModelTableFromOrigModel() const
{
map> lags;
map>> orig_diff_var;
map> diff;
for (const auto &[model_name, eqns] : var_model_table.getEqNums())
{
set lhs;
vector> orig_diff_var_vec;
vector diff_vec;
for (auto eqn : eqns)
{
// Perform some sanity checks on the RHS
optional eqtag { equation_tags.getTagValueByEqnAndKey(eqn, "name") };
set> 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)
{
cerr << "ERROR: in Equation " << eqtag.value_or(to_string(eqn+1))
<< ". 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.value_or(to_string(eqn+1))
<< ". 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)
{
cerr << "ERROR: in Equation " << eqtag.value_or(to_string(eqn+1))
<< ". A VAR model may not have lagged or leaded exogenous variables on the RHS. " << endl;
exit(EXIT_FAILURE);
}
// save lhs variables
equations[eqn]->arg1->collectVARLHSVariable(lhs);
equations[eqn]->arg1->countDiffs() > 0 ?
diff_vec.push_back(true) : diff_vec.push_back(false);
if (diff_vec.back())
{
set> diff_set;
equations[eqn]->arg1->collectDynamicVariables(SymbolType::endogenous, diff_set);
if (diff_set.size() != 1)
{
cerr << "ERROR: problem getting variable for LHS diff operator in equation "
<< eqn << endl;
exit(EXIT_FAILURE);
}
orig_diff_var_vec.push_back(diff_set.begin()->first);
}
else
orig_diff_var_vec.push_back(nullopt);
}
if (eqns.size() != lhs.size())
{
cerr << "ERROR: The LHS variables of the VAR model are not unique" << endl;
exit(EXIT_FAILURE);
}
set lhs_lag_equiv;
for (const auto &lh : lhs)
{
auto [lag_equiv_repr, index] = lh->getLagEquivalenceClass();
lhs_lag_equiv.insert(lag_equiv_repr);
}
vector max_lag;
for (auto eqn : eqns)
max_lag.push_back(equations[eqn]->arg2->VarMaxLag(lhs_lag_equiv));
lags[model_name] = max_lag;
diff[model_name] = diff_vec;
orig_diff_var[model_name] = orig_diff_var_vec;
}
var_model_table.setDiff(move(diff));
var_model_table.setMaxLags(move(lags));
var_model_table.setOrigDiffVar(move(orig_diff_var));
}
vector
DynamicModel::getVARDerivIDs(int lhs_symb_id, int lead_lag) const
{
vector 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->type == AuxVarType::endoLag && avi->orig_symb_id.value() == lhs_symb_id
&& avi->orig_lead_lag.value() + 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->type == AuxVarType::diffLag)
{
diff_lag_depth++;
if (avi->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->orig_symb_id.value());
}
catch (SymbolTable::UnknownSymbolIDException)
{
break;
}
}
}
return deriv_ids;
}
void
DynamicModel::fillVarModelTableMatrices()
{
map, expr_t>> AR;
map, expr_t>> A0;
map> constants;
for (const auto &[model_name, eqns] : var_model_table.getEqNums())
{
const vector &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 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").value_or(to_string(eqns[i]+1)) << " 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").value_or(to_string(eqns[i]+1)) << " is not linear" << endl;
exit(EXIT_FAILURE);
}
A0[model_name][{ i, lhs_symb_id }] = d;
}
}
// Fill constants vector
// Constants are computed by replacing all (transformed) endos and exos by zero
constants[model_name] = {}; // Ensure that the map exists, even if constants are all zero
for (size_t i = 0; i < eqns.size(); i++)
{
auto rhs = equations[eqns[i]]->arg2;
map subst_table;
auto rhs_vars = var_model_table.getRhs(model_name)[i]; // All the (transformed) endogenous on RHS, as computed by updateVarAndTrendModel()
rhs->collectDynamicVariables(SymbolType::exogenous, rhs_vars); // Add exos
for (auto [symb_id, lag] : rhs_vars)
subst_table[AddVariable(symb_id, lag)] = Zero;
expr_t c = rhs->replaceVarsInEquation(subst_table);
if (c != Zero)
constants[model_name][i] = c;
}
}
}
var_model_table.setAR(move(AR));
var_model_table.setA0(move(A0));
var_model_table.setConstants(move(constants));
}
map, expr_t>>
DynamicModel::computeAutoregressiveMatrices() const
{
map, expr_t>> ARr;
for (const auto &[model_name, eqns] : trend_component_model_table.getNonTargetEqNums())
{
map, expr_t> AR;
const vector &lhs = trend_component_model_table.getNonTargetLhs(model_name);
for (int i{0};
auto eqn : eqns)
{
auto bopn = dynamic_cast(equations[eqn]->arg2);
bopn->fillAutoregressiveRow(i++, lhs, AR);
}
ARr[model_name] = AR;
}
return ARr;
}
void
DynamicModel::fillTrendComponentModelTable() const
{
map> eqnums, trend_eqnums, lhsr;
map> lhs_expr_tr;
map>>> rhsr;
for (const auto &[model_name, eqtags] : trend_component_model_table.getTargetEqTags())
{
vector trend_eqnumber;
for (const auto &eqtag : eqtags)
{
optional eqn { equation_tags.getEqnByTag("name", eqtag) };
if (!eqn)
{
cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
exit(EXIT_FAILURE);
}
trend_eqnumber.push_back(*eqn);
}
trend_eqnums[model_name] = trend_eqnumber;
}
for (const auto &[model_name, eqtags] : trend_component_model_table.getEqTags())
{
vector eqnumber, lhs;
vector lhs_expr_t;
vector>> rhs;
for (const auto &eqtag : eqtags)
{
set> lhs_set, lhs_tmp_set, rhs_set;
optional eqn { equation_tags.getEqnByTag("name", eqtag) };
if (!eqn)
{
cerr << "ERROR: no equation is named '" << eqtag << "'" << endl;
exit(EXIT_FAILURE);
}
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::endogenous, lhs_set);
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::exogenous, lhs_tmp_set);
equations[*eqn]->arg1->collectDynamicVariables(SymbolType::parameter, lhs_tmp_set);
if (lhs_set.size() != 1 || !lhs_tmp_set.empty())
{
cerr << "ERROR: in Equation " << eqtag
<< ". A trend component model may only have one endogenous variable on the LHS. " << endl;
exit(EXIT_FAILURE);
}
auto itlhs = lhs_set.begin();
if (itlhs->second != 0)
{
cerr << "ERROR: in Equation " << eqtag
<< ". The variable on the LHS of a trend component model may not appear with a lead or a lag. "
<< endl;
exit(EXIT_FAILURE);
}
eqnumber.push_back(*eqn);
lhs.push_back(itlhs->first);
lhs_set.clear();
set lhs_expr_t_set;
equations[*eqn]->arg1->collectVARLHSVariable(lhs_expr_t_set);
lhs_expr_t.push_back(*(lhs_expr_t_set.begin()));
equations[*eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_set);
rhs.push_back(rhs_set);
}
eqnums[model_name] = eqnumber;
lhsr[model_name] = lhs;
lhs_expr_tr[model_name] = lhs_expr_t;
rhsr[model_name] = rhs;
}
trend_component_model_table.setRhs(move(rhsr));
trend_component_model_table.setVals(move(eqnums), move(trend_eqnums), move(lhsr), move(lhs_expr_tr));
}
pair