Select Git revision
SubModel.cc
Sébastien Villemot authored
There is no reason to keep this distinction. Additionally, since the data structure is now symmetric with the MCE case, unify this “h” parameter vector with the “α” parameter vector.
SubModel.cc 41.52 KiB
/*
* Copyright © 2018-2021 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
*/
#include <algorithm>
#include "SubModel.hh"
#include "DynamicModel.hh"
TrendComponentModelTable::TrendComponentModelTable(SymbolTable &symbol_table_arg) :
symbol_table{symbol_table_arg}
{
}
void
TrendComponentModelTable::addTrendComponentModel(string name_arg,
vector<string> eqtags_arg,
vector<string> target_eqtags_arg)
{
if (isExistingTrendComponentModelName(name_arg))
{
cerr << "Error: a trend component model already exists with the name " << name_arg << endl;
exit(EXIT_FAILURE);
}
eqtags[name_arg] = move(eqtags_arg);
target_eqtags[name_arg] = move(target_eqtags_arg);
names.insert(move(name_arg));
}
void
TrendComponentModelTable::setVals(map<string, vector<int>> eqnums_arg, map<string, vector<int>> target_eqnums_arg,
map<string, vector<int>> lhs_arg, map<string, vector<expr_t>> lhs_expr_t_arg)
{
eqnums = move(eqnums_arg);
target_eqnums = move(target_eqnums_arg);
lhs = move(lhs_arg);
lhs_expr_t = move(lhs_expr_t_arg);
for (const auto &it : eqnums)
{
vector<int> nontrend_vec;
for (auto eq : it.second)
if (find(target_eqnums[it.first].begin(), target_eqnums[it.first].end(), eq) == target_eqnums[it.first].end())
nontrend_vec.push_back(eq);
nontarget_eqnums[it.first] = nontrend_vec;
}
for (const auto &name : names)
{
vector<int> nontarget_lhs_vec, target_lhs_vec;
vector<int> lhsv = getLhs(name);
vector<int> eqnumsv = getEqNums(name);
for (int nontrend_it : getNonTargetEqNums(name))
nontarget_lhs_vec.push_back(lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), nontrend_it))));
nontarget_lhs[name] = nontarget_lhs_vec;
for (int trend_it : getTargetEqNums(name))
target_lhs_vec.push_back(lhsv.at(distance(eqnumsv.begin(), find(eqnumsv.begin(), eqnumsv.end(), trend_it))));
target_lhs[name] = target_lhs_vec;
}
}
void
TrendComponentModelTable::setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg)
{
rhs = move(rhs_arg);
}
void
TrendComponentModelTable::setTargetVar(map<string, vector<int>> target_vars_arg)
{
target_vars = move(target_vars_arg);
}
void
TrendComponentModelTable::setMaxLags(map<string, vector<int>> max_lags_arg)
{
max_lags = move(max_lags_arg);
}
void
TrendComponentModelTable::setDiff(map<string, vector<bool>> diff_arg)
{
diff = move(diff_arg);
}
void
TrendComponentModelTable::setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg)
{
orig_diff_var = move(orig_diff_var_arg);
}
void
TrendComponentModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
{
AR = move(AR_arg);
}
void
TrendComponentModelTable::setA0(map<string, map<tuple<int, int, int>, expr_t>> A0_arg,
map<string, map<tuple<int, int, int>, expr_t>> A0star_arg)
{
A0 = move(A0_arg);
A0star = move(A0star_arg);
}
const map<string, vector<string>> &
TrendComponentModelTable::getEqTags() const
{
return eqtags;
}
const vector<string> &
TrendComponentModelTable::getEqTags(const string &name_arg) const
{
checkModelName(name_arg);
return eqtags.find(name_arg)->second;
}
void
TrendComponentModelTable::checkModelName(const string &name_arg) const
{
if (!isExistingTrendComponentModelName(name_arg))
{
cerr << name_arg
<< " is not a recognized equation tag of a trend component model equation" << endl;
exit(EXIT_FAILURE);
}
}
const vector<int> &
TrendComponentModelTable::getNonTargetLhs(const string &name_arg) const
{
checkModelName(name_arg);
return nontarget_lhs.find(name_arg)->second;
}
const vector<int> &
TrendComponentModelTable::getTargetLhs(const string &name_arg) const
{
checkModelName(name_arg);
return target_lhs.find(name_arg)->second;
}
const vector<int> &
TrendComponentModelTable::getLhs(const string &name_arg) const
{
checkModelName(name_arg);
return lhs.find(name_arg)->second;
}
const vector<expr_t> &
TrendComponentModelTable::getLhsExprT(const string &name_arg) const
{
checkModelName(name_arg);
return lhs_expr_t.find(name_arg)->second;
}
const map<string, vector<string>> &
TrendComponentModelTable::getTargetEqTags() const
{
return target_eqtags;
}
const map<string, vector<int>> &
TrendComponentModelTable::getEqNums() const
{
return eqnums;
}
const map<string, vector<int>> &
TrendComponentModelTable::getTargetEqNums() const
{
return target_eqnums;
}
const vector<int> &
TrendComponentModelTable::getTargetEqNums(const string &name_arg) const
{
checkModelName(name_arg);
return target_eqnums.find(name_arg)->second;
}
const map<string, vector<int>> &
TrendComponentModelTable::getNonTargetEqNums() const
{
return nontarget_eqnums;
}
const vector<int> &
TrendComponentModelTable::getNonTargetEqNums(const string &name_arg) const
{
checkModelName(name_arg);
return nontarget_eqnums.find(name_arg)->second;
}
const vector<int> &
TrendComponentModelTable::getEqNums(const string &name_arg) const
{
checkModelName(name_arg);
return eqnums.find(name_arg)->second;
}
const vector<int> &
TrendComponentModelTable::getMaxLags(const string &name_arg) const
{
checkModelName(name_arg);
return max_lags.find(name_arg)->second;
}
int
TrendComponentModelTable::getMaxLag(const string &name_arg) const
{
int max_lag_int = 0;
for (auto it : getMaxLags(name_arg))
max_lag_int = max(max_lag_int, it);
return max_lag_int;
}
const vector<bool> &
TrendComponentModelTable::getDiff(const string &name_arg) const
{
checkModelName(name_arg);
return diff.find(name_arg)->second;
}
const vector<int> &
TrendComponentModelTable::getOrigDiffVar(const string &name_arg) const
{
checkModelName(name_arg);
return orig_diff_var.find(name_arg)->second;
}
void
TrendComponentModelTable::writeOutput(const string &basename, ostream &output) const
{
if (names.empty())
return;
string filename = "+" + basename + "/trend_component_ar_a0.m";
ofstream ar_ec_output;
ar_ec_output.open(filename, ios::out | ios::binary);
if (!ar_ec_output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
ar_ec_output << "function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
<< "%function [AR, A0, A0star] = trend_component_ar_a0(model_name, params)" << endl
<< "% File automatically generated by the Dynare preprocessor" << endl << endl;
for (const auto &name : names)
{
output << "M_.trend_component." << name << ".model_name = '" << name << "';" << endl
<< "M_.trend_component." << name << ".eqtags = {";
for (const auto &it : eqtags.at(name))
output << "'" << it << "'; ";
output << "};" << endl
<< "M_.trend_component." << name << ".eqn = [";
for (auto it : eqnums.at(name))
output << it + 1 << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".targets = [";
for (auto it : eqnums.at(name))
if (find(target_eqnums.at(name).begin(), target_eqnums.at(name).end(), it)
== target_eqnums.at(name).end())
output << "false ";
else
output << "true ";
output << "];" << endl
<< "M_.trend_component." << name << ".lhs = [";
for (auto it : lhs.at(name))
output << symbol_table.getTypeSpecificID(it) + 1 << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".max_lag = [";
for (auto it : max_lags.at(name))
output << it << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".diff = [";
for (const auto &it : diff.at(name))
output << (it ? "true" : "false") << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".orig_diff_var = [";
for (auto it : orig_diff_var.at(name))
output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".nonstationary = [";
for (size_t i = 0; i < diff.at(name).size(); i++)
output << "true ";
output << "];" << endl;
int i = 1;
for (const auto &it : rhs.at(name))
{
output << "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
for (auto [var, lag] : it)
output << symbol_table.getTypeSpecificID(var) + 1 << " ";
output << "];" << endl
<< "M_.trend_component." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
for (auto [var, lag] : it)
output << lag << " ";
output << "];" << endl;
i++;
}
output << "M_.trend_component." << name << ".target_vars = [";
for (auto it : target_vars.at(name))
output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
output << "];" << endl;
vector<string> target_eqtags_vec = target_eqtags.at(name);
output << "M_.trend_component." << name << ".target_eqtags = {";
for (auto it : target_eqtags_vec)
output << "'" << it << "';";
output << "};" << endl;
vector<string> eqtags_vec = eqtags.at(name);
output << "M_.trend_component." << name << ".target_eqn = [";
for (auto it : target_eqtags_vec)
{
int i = 0;
for (auto it1 : eqtags_vec)
{
i++;
if (it == it1)
{
output << i << " ";
break;
}
}
}
output << "];" << endl;
vector<int> target_lhs_vec = getTargetLhs(name);
vector<int> nontarget_lhs_vec = getNonTargetLhs(name);
ar_ec_output << "if strcmp(model_name, '" << name << "')" << endl
<< " % AR" << endl
<< " AR = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << getMaxLag(name) << ");" << endl;
for (const auto &[key, expr] : AR.at(name))
{
auto [eqn, lag, lhs_symb_id] = key;
int colidx = static_cast<int>(distance(nontarget_lhs_vec.begin(), find(nontarget_lhs_vec.begin(), nontarget_lhs_vec.end(), lhs_symb_id)));
ar_ec_output << " AR(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
ar_ec_output << ";" << endl;
}
int a0_lag = 0;
for (const auto &[key, expr] : A0.at(name))
a0_lag = max(a0_lag, get<1>(key));
ar_ec_output << endl
<< " % A0" << endl
<< " A0 = zeros(" << nontarget_lhs_vec.size() << ", " << nontarget_lhs_vec.size() << ", " << a0_lag << ");" << endl;
for (const auto &[key, expr] : A0.at(name))
{
auto [eqn, lag, colidx] = key;
ar_ec_output << " A0(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
ar_ec_output << ";" << endl;
}
int a0star_lag = 0;
for (const auto &[key, expr] : A0star.at(name))
a0star_lag = max(a0star_lag, get<1>(key));
ar_ec_output << endl
<< " % A0star" << endl
<< " A0star = zeros(" << nontarget_lhs_vec.size() << ", " << target_lhs_vec.size() << ", " << a0star_lag << ");" << endl;
for (const auto &[key, expr] : A0star.at(name))
{
auto [eqn, lag, colidx] = key;
ar_ec_output << " A0star(" << eqn + 1 << ", " << colidx + 1 << ", " << lag << ") = ";
expr->writeOutput(ar_ec_output, ExprNodeOutputType::matlabDynamicModel);
ar_ec_output << ";" << endl;
}
ar_ec_output << " return" << endl
<< "end" << endl << endl;
}
ar_ec_output << "error([model_name ' is not a valid trend_component_model name'])" << endl
<< "end" << endl;
ar_ec_output.close();
}
void
TrendComponentModelTable::writeJsonOutput(ostream &output) const
{
for (const auto &name : names)
{
if (name != *names.begin())
output << ", ";
output << R"({"statementName": "trend_component_model",)"
<< R"("model_name": ")" << name << R"(",)"
<< R"("eqtags": [)";
for (const auto &it : eqtags.at(name))
{
output << R"(")" << it << R"(")";
if (&it != &eqtags.at(name).back())
output << ", ";
}
output << R"(], "target_eqtags": [)";
for (const auto &it : target_eqtags.at(name))
{
output << R"(")" << it << R"(")";
if (&it != &target_eqtags.at(name).back())
output << ", ";
}
output << "]}";
}
}
VarModelTable::VarModelTable(SymbolTable &symbol_table_arg) :
symbol_table{symbol_table_arg}
{
}
void
VarModelTable::addVarModel(string name_arg, bool structural_arg, vector<string> eqtags_arg)
{
if (isExistingVarModelName(name_arg))
{
cerr << "Error: a VAR model already exists with the name " << name_arg << endl;
exit(EXIT_FAILURE);
}
structural[name_arg] = structural_arg;
eqtags[name_arg] = move(eqtags_arg);
names.insert(move(name_arg));
}
void
VarModelTable::writeOutput(const string &basename, ostream &output) const
{
if (names.empty())
return;
string filename = "+" + basename + "/varmatrices.m";
ofstream ar_output;
ar_output.open(filename, ios::out | ios::binary);
if (!ar_output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
ar_output << "function [ar, a0, constants] = varmatrices(model_name, params, reducedform)" << endl
<< "% File automatically generated by the Dynare preprocessor" << endl << endl
<< "if nargin<3" << endl
<< " reducedform = false;" << endl
<< "end" << endl << endl;
for (const auto &name : names)
{
output << "M_.var." << name << ".model_name = '" << name << "';" << endl
<< "M_.var." << name << ".structural = " << (structural.at(name) ? "true" : "false") << ";" << endl
<< "M_.var." << name << ".eqtags = {";
for (const auto &it : eqtags.at(name))
output << "'" << it << "'; ";
output << "};" << endl
<< "M_.var." << name << ".eqn = [";
for (auto it : eqnums.at(name))
output << it + 1 << " ";
output << "];" << endl
<< "M_.var." << name << ".lhs = [";
for (auto it : lhs.at(name))
output << symbol_table.getTypeSpecificID(it) + 1 << " ";
output << "];" << endl
<< "M_.var." << name << ".max_lag = [";
for (auto it : max_lags.at(name))
output << it << " ";
output << "];" << endl
<< "M_.var." << name << ".diff = [";
for (const auto &it : diff.at(name))
output << (it ? "true" : "false") << " ";
output << "];" << endl
<< "M_.var." << name << ".nonstationary = M_.var." << name << ".diff;" << endl
<< "M_.var." << name << ".orig_diff_var = [";
for (auto it : orig_diff_var.at(name))
output << (it >= 0 ? symbol_table.getTypeSpecificID(it) + 1 : -1) << " ";
output << "];" << endl;
int i = 1;
for (const auto &it : rhs.at(name))
{
output << "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.var = [";
for (auto [var, lag] : it)
output << symbol_table.getTypeSpecificID(var) + 1 << " ";
output << "];" << endl
<< "M_.var." << name << ".rhs.vars_at_eq{" << i << "}.lag = [";
for (auto [var, lag] : it)
output << lag << " ";
output << "];" << endl;
i++;
}
vector<int> lhs = getLhsOrigIds(name);
ar_output << "if strcmp(model_name, '" << name << "')" << endl
<< " ar = zeros(" << lhs.size() << ", " << lhs.size() << ", " << getMaxLag(name) << ");" << endl;
for (const auto &[key, expr] : AR.at(name))
{
auto [eqn, lag, lhs_symb_id] = key;
int colidx = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
ar_output << " ar(" << eqn + 1 << "," << colidx + 1 << "," << lag << ") = ";
expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
ar_output << ";" << endl;
}
ar_output << " if nargout>1" << endl
<< " a0 = eye(" << lhs.size() << ");" << endl;
for (const auto &[key, expr] : A0.at(name))
{
auto [eqn, lhs_symb_id] = key;
int colidx = static_cast<int>(distance(lhs.begin(), find(lhs.begin(), lhs.end(), lhs_symb_id)));
if (eqn!=colidx)
{
ar_output << " a0(" << eqn + 1 << "," << colidx + 1 << ") = ";
expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
ar_output << ";" << endl;
}
}
ar_output << " if reducedform" << endl
<< " for i=1:" << getMaxLag(name) << endl
<< " ar(:,:,i) = a0\\ar(:,:,i);" << endl
<< " end" << endl
<< " if nargout<3" << endl
<< " a0 = eye(" << lhs.size() << ");" << endl
<< " end" << endl
<< " end" << endl
<< " if nargout>2" << endl
<< " constants = zeros(" << lhs.size() << ",1);" << endl;
for (auto [eqn, expr] : constants.at(name))
{
ar_output << " constants(" << eqn + 1 << ") = ";
expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
ar_output << ";" << endl;
}
ar_output << " end" << endl
<< " if reducedform" << endl
<< " constants = a0\\constants;" << endl
<< " a0 = eye(" << lhs.size() << ");" << endl
<< " end" << endl
<< " end" << endl
<< " return" << endl
<< "end" << endl << endl;
}
ar_output << "error('%s is not a valid var_model name', model_name)" << endl;
ar_output.close();
}
void
VarModelTable::writeJsonOutput(ostream &output) const
{
for (const auto &name : names)
{
if (name != *names.begin())
output << ", ";
output << R"({"statementName": "var_model",)"
<< R"("model_name": ")" << name << R"(",)"
<< R"("eqtags": [)";
for (const auto &it : eqtags.at(name))
{
output << R"(")" << it << R"(")";
if (&it != &eqtags.at(name).back())
output << ", ";
}
output << "]}";
}
}
const map<string, bool> &
VarModelTable::getStructural() const
{
return structural;
}
const map<string, vector<string>> &
VarModelTable::getEqTags() const
{
return eqtags;
}
const vector<string> &
VarModelTable::getEqTags(const string &name_arg) const
{
checkModelName(name_arg);
return eqtags.find(name_arg)->second;
}
void
VarModelTable::checkModelName(const string &name_arg) const
{
if (!isExistingVarModelName(name_arg))
{
cerr << name_arg
<< " is not a recognized equation tag of a VAR model equation" << endl;
exit(EXIT_FAILURE);
}
}
void
VarModelTable::setEqNums(map<string, vector<int>> eqnums_arg)
{
eqnums = move(eqnums_arg);
}
void
VarModelTable::setLhs(map<string, vector<int>> lhs_arg)
{
lhs = move(lhs_arg);
for (auto it : lhs)
{
vector<int> lhsvec;
for (auto ids : it.second)
{
int lhs_last_orig_symb_id = ids;
int lhs_orig_symb_id = ids;
if (symbol_table.isAuxiliaryVariable(lhs_orig_symb_id))
try
{
lhs_last_orig_symb_id = lhs_orig_symb_id;
lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_orig_symb_id);
}
catch (...)
{
}
lhsvec.emplace_back(lhs_last_orig_symb_id);
}
lhs_orig_symb_ids[it.first] = lhsvec;
}
}
void
VarModelTable::setRhs(map<string, vector<set<pair<int, int>>>> rhs_arg)
{
rhs = move(rhs_arg);
}
void
VarModelTable::setLhsExprT(map<string, vector<expr_t>> lhs_expr_t_arg)
{
lhs_expr_t = move(lhs_expr_t_arg);
}
const map<string, vector<int>> &
VarModelTable::getEqNums() const
{
return eqnums;
}
const vector<bool> &
VarModelTable::getDiff(const string &name_arg) const
{
checkModelName(name_arg);
return diff.find(name_arg)->second;
}
const vector<int> &
VarModelTable::getEqNums(const string &name_arg) const
{
checkModelName(name_arg);
return eqnums.find(name_arg)->second;
}
void
VarModelTable::setMaxLags(map<string, vector<int>> max_lags_arg)
{
max_lags = move(max_lags_arg);
}
void
VarModelTable::setDiff(map<string, vector<bool>> diff_arg)
{
diff = move(diff_arg);
}
void
VarModelTable::setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg)
{
orig_diff_var = move(orig_diff_var_arg);
}
void
VarModelTable::setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg)
{
AR = move(AR_arg);
}
void
VarModelTable::setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg)
{
A0 = move(A0_arg);
}
void
VarModelTable::setConstants(map<string, map<int, expr_t>> constants_arg)
{
constants = move(constants_arg);
}
const vector<int> &
VarModelTable::getMaxLags(const string &name_arg) const
{
checkModelName(name_arg);
return max_lags.find(name_arg)->second;
}
int
VarModelTable::getMaxLag(const string &name_arg) const
{
int max_lag_int = 0;
for (auto it : getMaxLags(name_arg))
max_lag_int = max(max_lag_int, it);
return max_lag_int;
}
const vector<int> &
VarModelTable::getLhs(const string &name_arg) const
{
checkModelName(name_arg);
return lhs.find(name_arg)->second;
}
const vector<int> &
VarModelTable::getLhsOrigIds(const string &name_arg) const
{
checkModelName(name_arg);
return lhs_orig_symb_ids.find(name_arg)->second;
}
const vector<set<pair<int, int>>> &
VarModelTable::getRhs(const string &name_arg) const
{
checkModelName(name_arg);
return rhs.find(name_arg)->second;
}
const vector<expr_t> &
VarModelTable::getLhsExprT(const string &name_arg) const
{
checkModelName(name_arg);
return lhs_expr_t.find(name_arg)->second;
}
PacModelTable::PacModelTable(SymbolTable &symbol_table_arg) :
symbol_table{symbol_table_arg}
{
}
void
PacModelTable::addPacModel(string name_arg, string aux_model_name_arg, string discount_arg, expr_t growth_arg)
{
if (isExistingPacModelName(name_arg))
{
cerr << "Error: a PAC model already exists with the name " << name_arg << endl;
exit(EXIT_FAILURE);
}
aux_model_name[name_arg] = move(aux_model_name_arg);
discount[name_arg] = move(discount_arg);
growth[name_arg] = growth_arg;
original_growth[name_arg] = growth_arg;
names.insert(move(name_arg));
}
bool
PacModelTable::isExistingPacModelName(const string &name_arg) const
{
return names.find(name_arg) != names.end();
}
bool
PacModelTable::empty() const
{
return names.empty();
}
void
PacModelTable::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
for (auto &[name, gv] : growth)
if (gv)
gv->collectVariables(SymbolType::exogenous, mod_file_struct.pac_params);
}
void
PacModelTable::findDiffNodesInGrowth(lag_equivalence_table_t &diff_nodes) const
{
for (auto &[name, gv] : growth)
if (gv)
gv->findDiffNodes(diff_nodes);
}
void
PacModelTable::substituteDiffNodesInGrowth(const lag_equivalence_table_t &diff_nodes, ExprNode::subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs)
{
for (auto &[name, gv] : growth)
if (gv)
gv = gv->substituteDiff(diff_nodes, diff_subst_table, neweqs);
}
void
PacModelTable::transformPass(ExprNode::subst_table_t &diff_subst_table,
DynamicModel &dynamic_model, const VarModelTable &var_model_table,
const TrendComponentModelTable &trend_component_model_table)
{
// model name → expression for pac_expectation
map<string, expr_t> pac_expectation_substitution;
for (const auto &name : names)
{
/* Fill the growth_info structure.
Cannot be done in an earlier pass since growth terms can be
transformed by DynamicModel::substituteDiff(). */
if (growth[name])
try
{
growth_info[name] = growth[name]->matchLinearCombinationOfVariables(false);
}
catch (ExprNode::MatchFailureException &e)
{
auto gv = dynamic_cast<const VariableNode *>(growth[name]);
if (gv)
growth_info[name].emplace_back(gv->symb_id, gv->lag, -1, 1);
else
{
cerr << "Pac growth must be a linear combination of variables" << endl;
exit(EXIT_FAILURE);
}
}
// Collect some information about PAC models
int max_lag;
if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name[name]))
{
aux_model_type[name] = "trend_component";
max_lag = trend_component_model_table.getMaxLag(aux_model_name[name]) + 1;
lhs[name] = dynamic_model.getUndiffLHSForPac(aux_model_name[name], diff_subst_table);
}
else if (var_model_table.isExistingVarModelName(aux_model_name[name]))
{
aux_model_type[name] = "var";
max_lag = var_model_table.getMaxLag(aux_model_name[name]);
lhs[name] = var_model_table.getLhs(aux_model_name[name]);
}
else if (aux_model_name[name].empty())
max_lag = 0;
else
{
cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
exit(EXIT_FAILURE);
}
dynamic_model.analyzePacEquationStructure(name, eq_name, equation_info);
// Declare parameter for growth neutrality correction
if (growth[name])
{
string param_name = name + "_pac_growth_neutrality_correction";
try
{
int param_idx = symbol_table.addSymbol(param_name, SymbolType::parameter);
growth_neutrality_params[name] = param_idx;
}
catch (SymbolTable::AlreadyDeclaredException)
{
cerr << "The variable/parameter '" << param_name << "' conflicts with the auxiliary parameter that will be generated for the growth neutrality correction of the '" << name << "' PAC model. Please rename that parameter." << endl;
exit(EXIT_FAILURE);
}
}
// Compute the expressions that will be substituted for the pac_expectation operators
expr_t growth_correction_term = dynamic_model.Zero;
if (growth[name])
growth_correction_term = dynamic_model.AddTimes(growth[name], dynamic_model.AddVariable(growth_neutrality_params[name]));
if (aux_model_name[name].empty())
dynamic_model.computePacModelConsistentExpectationSubstitution(name,
symbol_table.getID(discount[name]),
pacEquationMaxLag(name),
growth_correction_term,
diff_subst_table,
aux_var_symb_ids,
aux_param_symb_ids,
pac_expectation_substitution);
else
dynamic_model.computePacBackwardExpectationSubstitution(name, lhs[name], max_lag,
aux_model_type[name],
growth_correction_term,
aux_var_symb_ids,
aux_param_symb_ids,
pac_expectation_substitution);
}
// Actually perform the substitution of pac_expectation
dynamic_model.substitutePacExpectation(pac_expectation_substitution, eq_name);
dynamic_model.checkNoRemainingPacExpectation();
}
void
PacModelTable::writeOutput(const string &basename, ostream &output) const
{
for (const auto &name : names)
{
output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name.at(name) << "';" << endl
<< "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount.at(name)) + 1 << ";" << endl;
if (growth.at(name))
{
output << "M_.pac." << name << ".growth_str = '";
original_growth.at(name)->writeJsonOutput(output, {}, {}, true);
output << "';" << endl;
int i = 0;
for (auto [growth_symb_id, growth_lag, param_id, constant] : growth_info.at(name))
{
string structname = "M_.pac." + name + ".growth_linear_comb(" + to_string(++i) + ").";
if (growth_symb_id >= 0)
{
string var_field = "endo_id";
if (symbol_table.getType(growth_symb_id) == SymbolType::exogenous)
{
var_field = "exo_id";
output << structname << "endo_id = 0;" << endl;
}
else
output << structname << "exo_id = 0;" << endl;
try
{
// case when this is not the highest lag of the growth variable
int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, growth_lag);
output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
<< structname << "lag = 0;" << endl;
}
catch (...)
{
try
{
// case when this is the highest lag of the growth variable
int tmp_growth_lag = growth_lag + 1;
int aux_symb_id = symbol_table.searchAuxiliaryVars(growth_symb_id, tmp_growth_lag);
output << structname << var_field << " = " << symbol_table.getTypeSpecificID(aux_symb_id) + 1 << ";" << endl
<< structname << "lag = -1;" << endl;
}
catch (...)
{
// case when there is no aux var for the variable
output << structname << var_field << " = "<< symbol_table.getTypeSpecificID(growth_symb_id) + 1 << ";" << endl
<< structname << "lag = " << growth_lag << ";" << endl;
}
}
}
else
output << structname << "endo_id = 0;" << endl
<< structname << "exo_id = 0;" << endl
<< structname << "lag = 0;" << endl;
output << structname << "param_id = "
<< (param_id == -1 ? 0 : symbol_table.getTypeSpecificID(param_id) + 1) << ";" << endl
<< structname << "constant = " << constant << ";" << endl;
}
}
}
// Write the auxiliary parameter IDs created for the pac_expectation operator
for (auto &[name, ids] : aux_param_symb_ids)
{
output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.alpha" : "h_param_indices") << " = [";
for (auto id : ids)
output << symbol_table.getTypeSpecificID(id) + 1 << " ";
output << "];" << endl;
}
// Write the auxiliary variable IDs created for the pac_expectation operator
for (auto &[name, id] : aux_var_symb_ids)
output << "M_.pac." << name << "." << (aux_model_name.at(name).empty() ? "mce.z1" : "aux_id")
<< " = " << symbol_table.getTypeSpecificID(id) + 1 << ";" << endl;
// Write PAC equation name info
for (auto &[name, eq] : eq_name)
output << "M_.pac." << name << ".eq_name = '" << eq << "';" << endl;
for (auto &[model, growth_neutrality_param_index] : growth_neutrality_params)
output << "M_.pac." << model << ".growth_neutrality_param_index = "
<< symbol_table.getTypeSpecificID(growth_neutrality_param_index) + 1 << ";" << endl;
for (auto &[model, lhs] : lhs)
{
output << "M_.pac." << model << ".lhs = [";
for (auto id : lhs)
output << id + 1 << " ";
output << "];" << endl;
}
for (auto &[model, type] : aux_model_type)
output << "M_.pac." << model << ".auxiliary_model_type = '" << type << "';" << endl;
for (auto &[name, val] : equation_info)
{
auto [lhs_pac_var, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants] = val;
output << "M_.pac." << name << ".lhs_var = "
<< symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl;
if (optim_share_index >= 0)
output << "M_.pac." << name << ".share_of_optimizing_agents_index = "
<< symbol_table.getTypeSpecificID(optim_share_index) + 1 << ";" << endl;
output << "M_.pac." << name << ".ec.params = "
<< symbol_table.getTypeSpecificID(ec_params_and_vars.first) + 1 << ";" << endl
<< "M_.pac." << name << ".ec.vars = [";
for (auto it : ec_params_and_vars.second)
output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
output << "];" << endl
<< "M_.pac." << name << ".ec.istarget = [";
for (auto it : ec_params_and_vars.second)
output << (get<1>(it) ? "true " : "false ");
output << "];" << endl
<< "M_.pac." << name << ".ec.scale = [";
for (auto it : ec_params_and_vars.second)
output << get<2>(it) << " ";
output << "];" << endl
<< "M_.pac." << name << ".ec.isendo = [";
for (auto it : ec_params_and_vars.second)
switch (symbol_table.getType(get<0>(it)))
{
case SymbolType::endogenous:
output << "true ";
break;
case SymbolType::exogenous:
output << "false ";
break;
default:
cerr << "expecting endogenous or exogenous" << endl;
exit(EXIT_FAILURE);
}
output << "];" << endl
<< "M_.pac." << name << ".ar.params = [";
for (auto &[pid, vid, vlag] : ar_params_and_vars)
output << (pid != -1 ? symbol_table.getTypeSpecificID(pid) + 1 : -1) << " ";
output << "];" << endl
<< "M_.pac." << name << ".ar.vars = [";
for (auto &[pid, vid, vlag] : ar_params_and_vars)
output << (vid != -1 ? symbol_table.getTypeSpecificID(vid) + 1 : -1) << " ";
output << "];" << endl
<< "M_.pac." << name << ".ar.lags = [";
for (auto &[pid, vid, vlag] : ar_params_and_vars)
output << vlag << " ";
output << "];" << endl
<< "M_.pac." << name << ".max_lag = " << pacEquationMaxLag(name) << ";" << endl;
if (!non_optim_vars_params_and_constants.empty())
{
output << "M_.pac." << name << ".non_optimizing_behaviour.params = [";
for (auto &it : non_optim_vars_params_and_constants)
if (get<2>(it) >= 0)
output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
else
output << "NaN ";
output << "];" << endl
<< "M_.pac." << name << ".non_optimizing_behaviour.vars = [";
for (auto &it : non_optim_vars_params_and_constants)
output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
output << "];" << endl
<< "M_.pac." << name << ".non_optimizing_behaviour.isendo = [";
for (auto &it : non_optim_vars_params_and_constants)
switch (symbol_table.getType(get<0>(it)))
{
case SymbolType::endogenous:
output << "true ";
break;
case SymbolType::exogenous:
output << "false ";
break;
default:
cerr << "expecting endogenous or exogenous" << endl;
exit(EXIT_FAILURE);
}
output << "];" << endl
<< "M_.pac." << name << ".non_optimizing_behaviour.lags = [";
for (auto &it : non_optim_vars_params_and_constants)
output << get<1>(it) << " ";
output << "];" << endl
<< "M_.pac." << name << ".non_optimizing_behaviour.scaling_factor = [";
for (auto &it : non_optim_vars_params_and_constants)
output << get<3>(it) << " ";
output << "];" << endl;
}
if (!additive_vars_params_and_constants.empty())
{
output << "M_.pac." << name << ".additive.params = [";
for (auto &it : additive_vars_params_and_constants)
if (get<2>(it) >= 0)
output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
else
output << "NaN ";
output << "];" << endl
<< "M_.pac." << name << ".additive.vars = [";
for (auto &it : additive_vars_params_and_constants)
output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
output << "];" << endl
<< "M_.pac." << name << ".additive.isendo = [";
for (auto &it : additive_vars_params_and_constants)
switch (symbol_table.getType(get<0>(it)))
{
case SymbolType::endogenous:
output << "true ";
break;
case SymbolType::exogenous:
output << "false ";
break;
default:
cerr << "expecting endogenous or exogenous" << endl;
exit(EXIT_FAILURE);
}
output << "];" << endl
<< "M_.pac." << name << ".additive.lags = [";
for (auto &it : additive_vars_params_and_constants)
output << get<1>(it) << " ";
output << "];" << endl
<< "M_.pac." << name << ".additive.scaling_factor = [";
for (auto &it : additive_vars_params_and_constants)
output << get<3>(it) << " ";
output << "];" << endl;
}
if (!optim_additive_vars_params_and_constants.empty())
{
output << "M_.pac." << name << ".optim_additive.params = [";
for (auto &it : optim_additive_vars_params_and_constants)
if (get<2>(it) >= 0)
output << symbol_table.getTypeSpecificID(get<2>(it)) + 1 << " ";
else
output << "NaN ";
output << "];" << endl
<< "M_.pac." << name << ".optim_additive.vars = [";
for (auto &it : optim_additive_vars_params_and_constants)
output << symbol_table.getTypeSpecificID(get<0>(it)) + 1 << " ";
output << "];" << endl
<< "M_.pac." << name << ".optim_additive.isendo = [";
for (auto &it : optim_additive_vars_params_and_constants)
switch (symbol_table.getType(get<0>(it)))
{
case SymbolType::endogenous:
output << "true ";
break;
case SymbolType::exogenous:
output << "false ";
break;
default:
cerr << "expecting endogenous or exogenous" << endl;
exit(EXIT_FAILURE);
}
output << "];" << endl
<< "M_.pac." << name << ".optim_additive.lags = [";
for (auto &it : optim_additive_vars_params_and_constants)
output << get<1>(it) << " ";
output << "];" << endl
<< "M_.pac." << name << ".optim_additive.scaling_factor = [";
for (auto &it : optim_additive_vars_params_and_constants)
output << get<3>(it) << " ";
output << "];" << endl;
}
}
}
void
PacModelTable::writeJsonOutput(ostream &output) const
{
for (const auto &name : names)
{
if (name != *names.begin())
output << ", ";
output << R"({"statementName": "pac_model",)"
<< R"("model_name": ")" << name << R"(",)"
<< R"("auxiliary_model_name": ")" << aux_model_name.at(name) << R"(",)"
<< R"("discount_index": )" << symbol_table.getTypeSpecificID(discount.at(name)) + 1;
if (growth.at(name))
{
output << R"(,"growth_str": ")";
original_growth.at(name)->writeJsonOutput(output, {}, {}, true);
output << R"(")";
}
output << "}" << endl;
}
}
int
PacModelTable::pacEquationMaxLag(const string &name_arg) const
{
return get<2>(equation_info.at(name_arg)).size();
}