/*
* 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;
}
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