4024 lines
147 KiB
C++
4024 lines
147 KiB
C++
/*
|
||
* 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 <https://www.gnu.org/licenses/>.
|
||
*/
|
||
|
||
#include <algorithm>
|
||
#include <cassert>
|
||
#include <cmath>
|
||
#include <cstdlib>
|
||
#include <iostream>
|
||
#include <numeric>
|
||
#include <ranges>
|
||
#include <regex>
|
||
#include <sstream>
|
||
#include <string_view>
|
||
|
||
#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<BinaryOpNode*>(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<int>(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)};
|
||
|
||
Bytecode::Writer code_file {basename + "/model/bytecode/dynamic.cod"};
|
||
|
||
// Declare temporary terms
|
||
code_file << Bytecode::FDIMT {static_cast<int>(temporary_terms_derivatives[0].size()
|
||
+ temporary_terms_derivatives[1].size())};
|
||
|
||
// Declare the (single) block
|
||
vector<int> 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<int>(count_if(
|
||
dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
|
||
[this](const auto& v) { return getTypeByDerivID(v.first) == SymbolType::endogenous; }))};
|
||
vector<int> eq_idx(equations.size());
|
||
iota(eq_idx.begin(), eq_idx.end(), 0);
|
||
vector<int> endo_idx(symbol_table.endo_nbr());
|
||
iota(endo_idx.begin(), endo_idx.end(), 0);
|
||
|
||
code_file << Bytecode::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<true>(code_file);
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeDynamicBlockBytecode(const string& basename) const
|
||
{
|
||
Bytecode::Writer 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 << Bytecode::FDIMT {static_cast<int>(blocks_temporary_terms_idxs.size())};
|
||
|
||
temporary_terms_t temporary_terms_written;
|
||
|
||
for (int block {0}; block < static_cast<int>(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 << Bytecode::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<int>(blocks_jacob_cols_endo[block].size())};
|
||
|
||
writeBlockBytecodeHelper<true>(code_file, block, temporary_terms_written);
|
||
}
|
||
code_file << Bytecode::FEND {};
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeDynamicMFile(const string& basename) const
|
||
{
|
||
auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::matlabDynamicModel>();
|
||
|
||
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<int>(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('\\', 0);
|
||
while (pos >= 0)
|
||
{
|
||
if (name.substr(pos + 1, 1) != R"(\)")
|
||
{
|
||
name = name.insert(pos, R"(\)");
|
||
pos++;
|
||
}
|
||
pos++;
|
||
pos = name.find('\\', 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<int>(
|
||
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<map<string, string>>
|
||
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<map<string, string>> 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.push_back({{tagname, string {str}}});
|
||
}
|
||
return eq_tag_set;
|
||
}
|
||
|
||
vector<int>
|
||
DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
|
||
bool excluded_vars_change_type,
|
||
vector<BinaryOpNode*>& all_equations,
|
||
vector<optional<int>>& 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<int> listed_eqs_by_number;
|
||
for (auto it = listed_eqs_by_tag.begin(); it != listed_eqs_by_tag.end();)
|
||
if (auto tmp = all_equation_tags.getEqnsByTags(*it); !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<int> 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<BinaryOpNode*> new_equations;
|
||
vector<optional<int>> new_equations_lineno;
|
||
map<int, int> old_eqn_num_2_new;
|
||
vector<int> 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<int> 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<map<string, string>>& 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<int> 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& m : listed_eqs_by_tag)
|
||
{
|
||
cerr << " ";
|
||
if (m.size() > 1)
|
||
cerr << "[ ";
|
||
bool first_printed {false};
|
||
for (const auto& [tagname, tagvalue] : m)
|
||
{
|
||
if (exchange(first_printed, true))
|
||
cerr << ", ";
|
||
cerr << tagname << "=" << tagvalue;
|
||
}
|
||
if (m.size() > 1)
|
||
cerr << " ]";
|
||
cerr << 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<int> 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<int>(blocks.size()); blk++)
|
||
{
|
||
int block_size = blocks[blk].size;
|
||
output << "M_.block_structure.block(" << blk + 1
|
||
<< ").Simulation_Type = " << static_cast<int>(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<int> 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<int>(!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<true>(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<int, set<pair<int, int>>> 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<int> 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<DynamicModel*>(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<int>(NNZDerivatives.size()); i++)
|
||
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
|
||
output << "];" << endl;
|
||
|
||
writeDriverSparseIndicesHelper<true>(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<string, vector<optional<int>>> trend_varr;
|
||
map<string, vector<set<pair<int, int>>>> rhsr;
|
||
for (const auto& [model_name, eqns] :
|
||
(var ? var_model_table.getEqNums() : trend_component_model_table.getEqNums()))
|
||
{
|
||
vector<int> lhs, trend_lhs;
|
||
vector<optional<int>> trend_var;
|
||
vector<set<pair<int, int>>> 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<pair<int, int>> 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<int> 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(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<string, vector<int>> eqnums, lhsr;
|
||
map<string, vector<expr_t>> lhs_expr_tr;
|
||
map<string, vector<set<pair<int, int>>>> rhsr;
|
||
|
||
for (const auto& [model_name, eqtags] : var_model_table.getEqTags())
|
||
{
|
||
vector<int> eqnumber, lhs;
|
||
vector<expr_t> lhs_expr_t;
|
||
vector<set<pair<int, int>>> rhs;
|
||
|
||
for (const auto& eqtag : eqtags)
|
||
{
|
||
set<pair<int, int>> lhs_set, lhs_tmp_set, rhs_set;
|
||
optional<int> 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<expr_t> 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<string, vector<int>> lags;
|
||
map<string, vector<optional<int>>> orig_diff_var;
|
||
map<string, vector<bool>> diff;
|
||
for (const auto& [model_name, eqns] : var_model_table.getEqNums())
|
||
{
|
||
set<expr_t> lhs;
|
||
vector<optional<int>> orig_diff_var_vec;
|
||
vector<bool> diff_vec;
|
||
for (auto eqn : eqns)
|
||
{
|
||
// Perform some sanity checks on the RHS
|
||
optional<string> eqtag {equation_tags.getTagValueByEqnAndKey(eqn, "name")};
|
||
set<pair<int, int>> rhs_endo_set, rhs_exo_set;
|
||
equations[eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_endo_set);
|
||
for (const auto& [symb_id, lag] : rhs_endo_set)
|
||
if (lag > 0)
|
||
{
|
||
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<pair<int, int>> 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.emplace_back(diff_set.begin()->first);
|
||
}
|
||
else
|
||
orig_diff_var_vec.emplace_back(nullopt);
|
||
}
|
||
|
||
if (eqns.size() != lhs.size())
|
||
{
|
||
cerr << "ERROR: The LHS variables of the VAR model are not unique" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
set<expr_t> lhs_lag_equiv;
|
||
for (const auto& lh : lhs)
|
||
{
|
||
auto [lag_equiv_repr, index] = lh->getLagEquivalenceClass();
|
||
lhs_lag_equiv.insert(lag_equiv_repr);
|
||
}
|
||
|
||
vector<int> 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<int>
|
||
DynamicModel::getVARDerivIDs(int lhs_symb_id, int lead_lag) const
|
||
{
|
||
vector<int> deriv_ids;
|
||
|
||
// First directly look for the variable itself
|
||
if (auto it = deriv_id_table.find({lhs_symb_id, lead_lag}); it != deriv_id_table.end())
|
||
deriv_ids.push_back(it->second);
|
||
|
||
// Then go through auxiliary variables
|
||
for (auto& [key, deriv_id2] : deriv_id_table)
|
||
{
|
||
auto [symb_id2, lead_lag2] = key;
|
||
const AuxVarInfo* avi;
|
||
try
|
||
{
|
||
avi = &symbol_table.getAuxVarInfo(symb_id2);
|
||
}
|
||
catch (SymbolTable::UnknownSymbolIDException)
|
||
{
|
||
continue;
|
||
}
|
||
|
||
if (avi->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<string, map<tuple<int, int, int>, expr_t>> AR;
|
||
map<string, map<tuple<int, int>, expr_t>> A0;
|
||
map<string, map<int, expr_t>> constants;
|
||
for (const auto& [model_name, eqns] : var_model_table.getEqNums())
|
||
{
|
||
const vector<int>& lhs = var_model_table.getLhs(model_name);
|
||
int max_lag = var_model_table.getMaxLag(model_name);
|
||
for (auto lhs_symb_id : lhs)
|
||
{
|
||
// Fill autoregressive matrix (AR)
|
||
for (int lag = 1; lag <= max_lag; lag++)
|
||
{
|
||
vector<int> deriv_ids = getVARDerivIDs(lhs_symb_id, -lag);
|
||
for (size_t i = 0; i < eqns.size(); i++)
|
||
{
|
||
expr_t d = Zero;
|
||
for (int deriv_id : deriv_ids)
|
||
d = AddPlus(d, equations[eqns[i]]->getDerivative(deriv_id));
|
||
if (d != Zero)
|
||
{
|
||
if (!d->isConstant())
|
||
{
|
||
cerr << "ERROR: Equation "
|
||
<< equation_tags.getTagValueByEqnAndKey(eqns[i], "name")
|
||
.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<VariableNode*, NumConstNode*> 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<string, map<tuple<int, int, int>, expr_t>>
|
||
DynamicModel::computeAutoregressiveMatrices() const
|
||
{
|
||
map<string, map<tuple<int, int, int>, expr_t>> ARr;
|
||
for (const auto& [model_name, eqns] : trend_component_model_table.getNonTargetEqNums())
|
||
{
|
||
map<tuple<int, int, int>, expr_t> AR;
|
||
const vector<int>& lhs = trend_component_model_table.getNonTargetLhs(model_name);
|
||
for (int i {0}; auto eqn : eqns)
|
||
{
|
||
auto bopn = dynamic_cast<BinaryOpNode*>(equations[eqn]->arg2);
|
||
bopn->fillAutoregressiveRow(i++, lhs, AR);
|
||
}
|
||
ARr[model_name] = AR;
|
||
}
|
||
return ARr;
|
||
}
|
||
|
||
void
|
||
DynamicModel::fillTrendComponentModelTable() const
|
||
{
|
||
map<string, vector<int>> eqnums, trend_eqnums, lhsr;
|
||
map<string, vector<expr_t>> lhs_expr_tr;
|
||
map<string, vector<set<pair<int, int>>>> rhsr;
|
||
for (const auto& [model_name, eqtags] : trend_component_model_table.getTargetEqTags())
|
||
{
|
||
vector<int> trend_eqnumber;
|
||
for (const auto& eqtag : eqtags)
|
||
{
|
||
optional<int> 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<int> eqnumber, lhs;
|
||
vector<expr_t> lhs_expr_t;
|
||
vector<set<pair<int, int>>> rhs;
|
||
|
||
for (const auto& eqtag : eqtags)
|
||
{
|
||
set<pair<int, int>> lhs_set, lhs_tmp_set, rhs_set;
|
||
optional<int> 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<expr_t> 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<map<string, map<tuple<int, int>, expr_t>>, map<string, map<tuple<int, int>, expr_t>>>
|
||
DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t& diff_subst_table) const
|
||
{
|
||
map<string, map<tuple<int, int>, expr_t>> A0r, A0starr;
|
||
|
||
for (const auto& [model_name, eqns] : trend_component_model_table.getEqNums())
|
||
{
|
||
map<tuple<int, int>, expr_t> A0, A0star;
|
||
vector<int> target_lhs = trend_component_model_table.getTargetLhs(model_name);
|
||
vector<int> nontarget_eqnums = trend_component_model_table.getNonTargetEqNums(model_name);
|
||
vector<int> undiff_nontarget_lhs = getUndiffLHSForPac(model_name, diff_subst_table);
|
||
vector<int> parsed_undiff_nontarget_lhs;
|
||
|
||
for (int i {0}; auto eqn : eqns)
|
||
{
|
||
if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
|
||
parsed_undiff_nontarget_lhs.push_back(undiff_nontarget_lhs.at(i));
|
||
i++;
|
||
}
|
||
|
||
for (int i {0}; auto eqn : eqns)
|
||
if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
|
||
equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs,
|
||
A0, A0star);
|
||
A0r[model_name] = A0;
|
||
A0starr[model_name] = A0star;
|
||
}
|
||
|
||
return {A0r, A0starr};
|
||
}
|
||
|
||
void
|
||
DynamicModel::fillTrendComponentModelTableFromOrigModel() const
|
||
{
|
||
map<string, vector<int>> lags;
|
||
map<string, vector<optional<int>>> orig_diff_var;
|
||
map<string, vector<bool>> diff;
|
||
for (const auto& [model_name, eqns] : trend_component_model_table.getEqNums())
|
||
{
|
||
set<expr_t> lhs;
|
||
vector<optional<int>> orig_diff_var_vec;
|
||
vector<bool> diff_vec;
|
||
for (auto eqn : eqns)
|
||
{
|
||
// Perform some sanity checks on the RHS
|
||
optional<string> eqtag {equation_tags.getTagValueByEqnAndKey(eqn, "name")};
|
||
set<pair<int, int>> rhs_endo_set, rhs_exo_set;
|
||
equations[eqn]->arg2->collectDynamicVariables(SymbolType::endogenous, rhs_endo_set);
|
||
for (const auto& [symb_id, lag] : rhs_endo_set)
|
||
if (lag >= 0)
|
||
{
|
||
cerr << "ERROR: in Equation " << eqtag.value_or(to_string(eqn + 1))
|
||
<< ". A trend component model may not have leaded or 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 trend component model may not have lagged or leaded exogenous "
|
||
"variables on the RHS. "
|
||
<< endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
// save lhs variables
|
||
equations[eqn]->arg1->collectVARLHSVariable(lhs);
|
||
|
||
if (equations[eqn]->arg1->countDiffs() > 0)
|
||
diff_vec.push_back(true);
|
||
else
|
||
diff_vec.push_back(false);
|
||
if (diff_vec.back())
|
||
{
|
||
set<pair<int, int>> 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.emplace_back(diff_set.begin()->first);
|
||
}
|
||
else
|
||
orig_diff_var_vec.emplace_back(nullopt);
|
||
}
|
||
|
||
if (eqns.size() != lhs.size())
|
||
{
|
||
cerr << "ERROR: The LHS variables of the trend component model are not unique" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
set<expr_t> lhs_lag_equiv;
|
||
for (const auto& lh : lhs)
|
||
{
|
||
auto [lag_equiv_repr, index] = lh->getLagEquivalenceClass();
|
||
lhs_lag_equiv.insert(lag_equiv_repr);
|
||
}
|
||
|
||
vector<int> 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;
|
||
}
|
||
trend_component_model_table.setDiff(move(diff));
|
||
trend_component_model_table.setMaxLags(move(lags));
|
||
trend_component_model_table.setOrigDiffVar(move(orig_diff_var));
|
||
}
|
||
|
||
void
|
||
DynamicModel::fillTrendComponentModelTableAREC(
|
||
const ExprNode::subst_table_t& diff_subst_table) const
|
||
{
|
||
auto ARr = computeAutoregressiveMatrices();
|
||
trend_component_model_table.setAR(move(ARr));
|
||
auto [A0r, A0starr] = computeErrorComponentMatrices(diff_subst_table);
|
||
trend_component_model_table.setA0(move(A0r), move(A0starr));
|
||
}
|
||
|
||
vector<int>
|
||
DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
|
||
const ExprNode::subst_table_t& diff_subst_table) const
|
||
{
|
||
vector<expr_t> lhs_expr_t = trend_component_model_table.getLhsExprT(aux_model_name);
|
||
vector<int> lhs = trend_component_model_table.getLhs(aux_model_name);
|
||
vector<bool> diff = trend_component_model_table.getDiff(aux_model_name);
|
||
vector<int> eqnumber = trend_component_model_table.getEqNums(aux_model_name);
|
||
vector<int> nontrend_eqnums = trend_component_model_table.getNonTargetEqNums(aux_model_name);
|
||
|
||
for (auto eqn : nontrend_eqnums)
|
||
{
|
||
auto i = distance(eqnumber.begin(), find(eqnumber.begin(), eqnumber.end(), eqn));
|
||
|
||
if (eqnumber[i] != eqn)
|
||
{
|
||
cerr << "ERROR: equation " << eqn << " not found in VAR" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
if (diff.at(i) != true)
|
||
{
|
||
cerr << "ERROR: the variable on the LHS of equation #" << eqn
|
||
<< " does not have the diff operator applied to it yet you are trying to undiff it."
|
||
<< endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
bool printerr = false;
|
||
expr_t node = nullptr;
|
||
expr_t aux_var = lhs_expr_t.at(i);
|
||
for (const auto& it : diff_subst_table)
|
||
if (it.second == aux_var)
|
||
{
|
||
node = const_cast<expr_t>(it.first);
|
||
break;
|
||
}
|
||
|
||
if (!node)
|
||
{
|
||
cerr << "Unexpected error encountered." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
node = node->undiff();
|
||
auto it1 = diff_subst_table.find(node);
|
||
if (it1 == diff_subst_table.end())
|
||
printerr = true;
|
||
|
||
if (printerr)
|
||
{ // we have undiffed something like diff(x), hence x is not in diff_subst_table
|
||
lhs_expr_t.at(i) = node;
|
||
lhs.at(i) = dynamic_cast<VariableNode*>(node)->symb_id;
|
||
}
|
||
else
|
||
{
|
||
lhs_expr_t.at(i) = const_cast<expr_t>(it1->first);
|
||
lhs.at(i) = const_cast<VariableNode*>(it1->second)->symb_id;
|
||
}
|
||
}
|
||
return lhs;
|
||
}
|
||
|
||
void
|
||
DynamicModel::analyzePacEquationStructure(const string& name, map<string, string>& pac_eq_name,
|
||
PacModelTable::equation_info_t& pac_equation_info)
|
||
{
|
||
for (auto& equation : equations)
|
||
if (equation->containsPacExpectation(name))
|
||
{
|
||
if (!pac_eq_name[name].empty())
|
||
{
|
||
cerr << "It is not possible to use 'pac_expectation(" << name
|
||
<< ")' in several equations." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
optional<string> eqn {
|
||
equation_tags.getTagValueByEqnAndKey(&equation - &equations[0], "name")};
|
||
if (!eqn)
|
||
{
|
||
cerr << "Every equation with a 'pac_expectation' operator must have been assigned an "
|
||
"equation tag name"
|
||
<< endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
pac_eq_name[name] = *eqn;
|
||
|
||
set<pair<int, int>> lhss;
|
||
equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss);
|
||
auto lhs = *lhss.begin();
|
||
int lhs_symb_id = lhs.first;
|
||
int lhs_orig_symb_id = lhs_symb_id;
|
||
if (symbol_table.isDiffAuxiliaryVariable(lhs_orig_symb_id))
|
||
try
|
||
{
|
||
lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_orig_symb_id);
|
||
}
|
||
catch (...)
|
||
{
|
||
}
|
||
|
||
auto arg2 = dynamic_cast<BinaryOpNode*>(equation->arg2);
|
||
if (!arg2)
|
||
{
|
||
cerr << "Pac equation in incorrect format" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
auto [optim_share_index, optim_part, non_optim_part, additive_part]
|
||
= arg2->getPacOptimizingShareAndExprNodes(lhs_orig_symb_id);
|
||
|
||
pair<int, vector<tuple<int, bool, int>>> ec_params_and_vars;
|
||
vector<tuple<optional<int>, optional<int>, int>> ar_params_and_vars;
|
||
vector<tuple<int, int, optional<int>, double>> non_optim_vars_params_and_constants,
|
||
optim_additive_vars_params_and_constants, additive_vars_params_and_constants;
|
||
if (!optim_part)
|
||
{
|
||
auto bopn = dynamic_cast<BinaryOpNode*>(equation->arg2);
|
||
if (!bopn)
|
||
{
|
||
cerr << "Error in PAC equation" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars,
|
||
additive_vars_params_and_constants);
|
||
}
|
||
else
|
||
{
|
||
auto bopn = dynamic_cast<BinaryOpNode*>(optim_part);
|
||
if (!bopn)
|
||
{
|
||
cerr << "Error in PAC equation" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
bopn->getPacAREC(lhs_symb_id, lhs_orig_symb_id, ec_params_and_vars, ar_params_and_vars,
|
||
optim_additive_vars_params_and_constants);
|
||
try
|
||
{
|
||
non_optim_vars_params_and_constants
|
||
= non_optim_part->matchLinearCombinationOfVariables();
|
||
if (additive_part)
|
||
additive_vars_params_and_constants
|
||
= additive_part->matchLinearCombinationOfVariables();
|
||
}
|
||
catch (ExprNode::MatchFailureException& e)
|
||
{
|
||
cerr << "Error in parsing non-optimizing agents or additive part of PAC equation: "
|
||
<< e.message << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
if (lhs.first == -1)
|
||
{
|
||
cerr << "analyzePacEquationStructure: error obtaining LHS variable." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
if (ec_params_and_vars.second.empty())
|
||
{
|
||
cerr << "analyzePacEquationStructure: error obtaining RHS parameters." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
pac_equation_info[name] = {lhs,
|
||
optim_share_index,
|
||
move(ar_params_and_vars),
|
||
move(ec_params_and_vars),
|
||
move(non_optim_vars_params_and_constants),
|
||
move(additive_vars_params_and_constants),
|
||
move(optim_additive_vars_params_and_constants)};
|
||
}
|
||
}
|
||
|
||
int
|
||
DynamicModel::getPacTargetSymbId(const string& pac_model_name) const
|
||
{
|
||
for (auto equation : equations)
|
||
if (equation->containsPacExpectation(pac_model_name))
|
||
{
|
||
set<pair<int, int>> lhss;
|
||
equation->arg1->collectDynamicVariables(SymbolType::endogenous, lhss);
|
||
if (lhss.size() != 1)
|
||
throw PacTargetNotIdentifiedException {pac_model_name,
|
||
"LHS must contain a single endogenous"};
|
||
int lhs_symb_id = lhss.begin()->first;
|
||
if (!symbol_table.isDiffAuxiliaryVariable(lhs_symb_id))
|
||
throw PacTargetNotIdentifiedException {pac_model_name, "LHS must be a diff operator"};
|
||
int undiff_lhs_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_symb_id);
|
||
auto barg2 = dynamic_cast<BinaryOpNode*>(equation->arg2);
|
||
if (!barg2)
|
||
throw PacTargetNotIdentifiedException {pac_model_name, "RHS must be a binary operator"};
|
||
auto [optim_share_index, optim_part, non_optim_part, additive_part]
|
||
= barg2->getPacOptimizingShareAndExprNodes(undiff_lhs_symb_id);
|
||
/* If there is an optimization part, restrict the search to that part,
|
||
since it contains the MCE . */
|
||
expr_t mce = optim_part ? optim_part : equation->arg2;
|
||
|
||
vector<pair<expr_t, int>> terms;
|
||
mce->decomposeAdditiveTerms(terms);
|
||
for (auto [term, sign] : terms)
|
||
try
|
||
{
|
||
auto [param_id, target_id]
|
||
= term->matchParamTimesTargetMinusVariable(undiff_lhs_symb_id);
|
||
return target_id;
|
||
}
|
||
catch (ExprNode::MatchFailureException&)
|
||
{
|
||
}
|
||
throw PacTargetNotIdentifiedException {pac_model_name,
|
||
"No term of the form parameter*(target-LHS_level)"};
|
||
}
|
||
throw PacTargetNotIdentifiedException {
|
||
pac_model_name, "No equation with the corresponding pac_expectation operator"};
|
||
}
|
||
|
||
void
|
||
DynamicModel::computePacModelConsistentExpectationSubstitution(
|
||
const string& name, int discount_symb_id, int pac_eq_max_lag, expr_t growth_correction_term,
|
||
string auxname, ExprNode::subst_table_t& diff_subst_table,
|
||
map<string, int>& pac_aux_var_symb_ids, map<string, vector<int>>& pac_aux_param_symb_ids,
|
||
map<string, expr_t>& pac_expectation_substitution)
|
||
{
|
||
int pac_target_symb_id;
|
||
try
|
||
{
|
||
pac_target_symb_id = getPacTargetSymbId(name);
|
||
}
|
||
catch (PacTargetNotIdentifiedException& e)
|
||
{
|
||
cerr << "Can't identify target for PAC model " << name << ": " << e.message;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
int neqs = 0;
|
||
|
||
// Create the endogenous representing Z₁ (no orig_expr is given since its definition is recursive)
|
||
if (auxname.empty())
|
||
auxname = "mce_Z1_" + name;
|
||
int mce_z1_symb_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, nullptr);
|
||
pac_aux_var_symb_ids[name] = mce_z1_symb_id;
|
||
|
||
expr_t A = One;
|
||
expr_t fp = Zero;
|
||
expr_t beta = AddVariable(discount_symb_id);
|
||
for (int i = 1; i <= pac_eq_max_lag + 1; i++)
|
||
{
|
||
string param_name = "mce_alpha_" + name + "_" + to_string(i);
|
||
try
|
||
{
|
||
int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
|
||
pac_aux_param_symb_ids[name].push_back(alpha_i_symb_id);
|
||
A = AddPlus(A, AddVariable(alpha_i_symb_id));
|
||
fp = AddPlus(fp, AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
|
||
AddPower(beta, AddPossiblyNegativeConstant(i))),
|
||
AddVariable(mce_z1_symb_id, i)));
|
||
}
|
||
catch (SymbolTable::AlreadyDeclaredException& e)
|
||
{
|
||
cerr << "The variable/parameter '" << param_name
|
||
<< "' conflicts with a parameter that will be generated for the '" << name
|
||
<< "' PAC model. Please rename it." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
// Add diff nodes and eqs for pac_target_symb_id
|
||
const VariableNode* target_base_diff_node;
|
||
auto create_target_lag = [&](int lag) {
|
||
if (symbol_table.isAuxiliaryVariable(pac_target_symb_id))
|
||
// We know it is a log, see ExprNode::matchParamTimesTargetMinusVariable()
|
||
return AddLog(AddVariable(symbol_table.getOrigSymbIdForAuxVar(pac_target_symb_id), lag));
|
||
else
|
||
return dynamic_cast<ExprNode*>(AddVariable(pac_target_symb_id, lag));
|
||
};
|
||
|
||
expr_t diff_node_to_search = AddDiff(create_target_lag(0));
|
||
if (auto sit = diff_subst_table.find(diff_node_to_search); sit != diff_subst_table.end())
|
||
target_base_diff_node = sit->second;
|
||
else
|
||
{
|
||
int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
|
||
target_base_diff_node = AddVariable(symb_id);
|
||
auto neweq = AddEqual(const_cast<VariableNode*>(target_base_diff_node),
|
||
AddMinus(create_target_lag(0), create_target_lag(-1)));
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
neqs++;
|
||
}
|
||
|
||
map<int, VariableNode*> target_aux_var_to_add;
|
||
const VariableNode* last_aux_var = target_base_diff_node;
|
||
for (int i = 1; i <= pac_eq_max_lag; i++, neqs++)
|
||
{
|
||
expr_t this_diff_node = AddDiff(create_target_lag(i));
|
||
int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node,
|
||
last_aux_var->symb_id, 1);
|
||
VariableNode* current_aux_var = AddVariable(symb_id);
|
||
auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1));
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
last_aux_var = current_aux_var;
|
||
target_aux_var_to_add[i] = current_aux_var;
|
||
}
|
||
|
||
expr_t fs = Zero;
|
||
for (int k = 1; k <= pac_eq_max_lag; k++)
|
||
{
|
||
expr_t ssum = Zero;
|
||
for (int j = k + 1; j <= pac_eq_max_lag + 1; j++)
|
||
{
|
||
int alpha_j_symb_id = -1;
|
||
string param_name = "mce_alpha_" + name + "_" + to_string(j);
|
||
try
|
||
{
|
||
alpha_j_symb_id = symbol_table.getID(param_name);
|
||
}
|
||
catch (SymbolTable::UnknownSymbolNameException& e)
|
||
{
|
||
alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
|
||
}
|
||
ssum = AddPlus(ssum, AddTimes(AddVariable(alpha_j_symb_id),
|
||
AddPower(beta, AddPossiblyNegativeConstant(j))));
|
||
}
|
||
fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k]));
|
||
}
|
||
auto neweq = AddEqual(
|
||
AddVariable(mce_z1_symb_id),
|
||
AddMinus(AddTimes(A, AddMinus(const_cast<VariableNode*>(target_base_diff_node), fs)), fp));
|
||
addEquation(neweq, nullopt);
|
||
/* This equation is not added to the list of auxiliary equations, because it
|
||
is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
|
||
neqs++;
|
||
|
||
cout << "PAC Model Consistent Expectation: added " << neqs
|
||
<< " auxiliary variables and equations for model " << name << "." << endl;
|
||
|
||
/* The growth correction term is not added to the definition of Z₁
|
||
because the latter is recursive. Rather put it at the level of the
|
||
substition of pac_expectation operator. */
|
||
pac_expectation_substitution[name] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction_term);
|
||
}
|
||
|
||
void
|
||
DynamicModel::computePacBackwardExpectationSubstitution(
|
||
const string& name, const vector<int>& lhs, int max_lag, const string& aux_model_type,
|
||
expr_t growth_correction_term, string auxname, map<string, int>& pac_aux_var_symb_ids,
|
||
map<string, vector<int>>& pac_aux_param_symb_ids,
|
||
map<string, expr_t>& pac_expectation_substitution)
|
||
{
|
||
auto create_aux_param = [&](const string& param_name) {
|
||
try
|
||
{
|
||
return symbol_table.addSymbol(param_name, SymbolType::parameter);
|
||
}
|
||
catch (SymbolTable::AlreadyDeclaredException)
|
||
{
|
||
cerr << "ERROR: the variable/parameter '" << param_name
|
||
<< "' conflicts with some auxiliary parameter that will be generated for the '" << name
|
||
<< "' PAC model. Please rename that parameter." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
};
|
||
|
||
expr_t subExpr = Zero;
|
||
if (aux_model_type == "var")
|
||
{
|
||
/* If the auxiliary model is a VAR, add a parameter corresponding
|
||
to the constant. */
|
||
int new_param_symb_id = create_aux_param("h_" + name + "_constant");
|
||
pac_aux_param_symb_ids[name].push_back(new_param_symb_id);
|
||
subExpr = AddPlus(subExpr, AddVariable(new_param_symb_id));
|
||
}
|
||
for (int i = 1; i < max_lag + 1; i++)
|
||
for (auto lhsit : lhs)
|
||
{
|
||
int new_param_symb_id = create_aux_param("h_" + name + "_var_" + symbol_table.getName(lhsit)
|
||
+ "_lag_" + to_string(i));
|
||
pac_aux_param_symb_ids[name].push_back(new_param_symb_id);
|
||
subExpr
|
||
= AddPlus(subExpr, AddTimes(AddVariable(new_param_symb_id), AddVariable(lhsit, -i)));
|
||
}
|
||
|
||
subExpr = AddPlus(subExpr, growth_correction_term);
|
||
|
||
if (auxname.empty())
|
||
auxname = "pac_expectation_" + name;
|
||
int expect_var_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, subExpr);
|
||
expr_t neweq = AddEqual(AddVariable(expect_var_id), subExpr);
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
pac_aux_var_symb_ids[name] = expect_var_id;
|
||
pac_expectation_substitution[name] = AddVariable(expect_var_id);
|
||
}
|
||
|
||
void
|
||
DynamicModel::computePacBackwardExpectationSubstitutionWithComponents(
|
||
const string& name, const vector<int>& lhs, int max_lag, const string& aux_model_type,
|
||
vector<PacModelTable::target_component_t>& pac_target_components,
|
||
map<string, expr_t>& pac_expectation_substitution)
|
||
{
|
||
auto create_aux_param = [&](const string& param_name) {
|
||
try
|
||
{
|
||
return symbol_table.addSymbol(param_name, SymbolType::parameter);
|
||
}
|
||
catch (SymbolTable::AlreadyDeclaredException)
|
||
{
|
||
cerr << "ERROR: the variable/parameter '" << param_name
|
||
<< "' conflicts with some auxiliary parameter that will be generated for the '" << name
|
||
<< "' PAC model. Please rename that parameter." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
};
|
||
|
||
expr_t substexpr = Zero;
|
||
|
||
for (int component_idx {1};
|
||
auto& [component, growth, auxname, kind, coeff, growth_neutrality_param, h_indices,
|
||
original_growth, growth_info] : pac_target_components)
|
||
{
|
||
string name_component = name + "_component" + to_string(component_idx);
|
||
|
||
// Create the linear combination of the variables from the auxiliary model
|
||
expr_t auxdef = Zero;
|
||
if (aux_model_type == "var")
|
||
{
|
||
/* If the auxiliary model is a VAR, add a parameter corresponding
|
||
to the constant. */
|
||
int new_param_symb_id = create_aux_param("h_" + name_component + "_constant");
|
||
h_indices.push_back(new_param_symb_id);
|
||
auxdef = AddPlus(auxdef, AddVariable(new_param_symb_id));
|
||
}
|
||
for (int i = 1; i < max_lag + 1; i++)
|
||
for (auto lhsit : lhs)
|
||
{
|
||
int new_param_symb_id
|
||
= create_aux_param("h_" + name_component + "_var_" + symbol_table.getName(lhsit)
|
||
+ "_lag_" + to_string(i));
|
||
h_indices.push_back(new_param_symb_id);
|
||
auxdef
|
||
= AddPlus(auxdef, AddTimes(AddVariable(new_param_symb_id), AddVariable(lhsit, -i)));
|
||
}
|
||
|
||
// If needed, add the growth neutrality correction for this component
|
||
if (growth)
|
||
{
|
||
growth_neutrality_param
|
||
= create_aux_param(name_component + "_pac_growth_neutrality_correction");
|
||
auxdef = AddPlus(auxdef, AddTimes(growth, AddVariable(growth_neutrality_param)));
|
||
}
|
||
else
|
||
growth_neutrality_param = -1;
|
||
|
||
// Create the auxiliary variable for this component
|
||
int aux_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, auxdef);
|
||
expr_t auxvar = AddVariable(aux_id);
|
||
|
||
// Add the equation defining the auxiliary variable for this component
|
||
expr_t neweq = AddEqual(auxvar, auxdef);
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
|
||
// Update the expression to be substituted for the pac_expectation operator
|
||
substexpr = AddPlus(substexpr, AddTimes(coeff, auxvar));
|
||
|
||
component_idx++;
|
||
}
|
||
|
||
pac_expectation_substitution[name] = substexpr;
|
||
}
|
||
|
||
void
|
||
DynamicModel::substitutePacExpectation(const map<string, expr_t>& pac_expectation_substitution,
|
||
const map<string, string>& pac_eq_name)
|
||
{
|
||
for (auto& [model_name, substexpr] : pac_expectation_substitution)
|
||
{
|
||
optional<int> eq {equation_tags.getEqnByTag("name", pac_eq_name.at(model_name))};
|
||
auto substeq = dynamic_cast<BinaryOpNode*>(
|
||
equations[eq.value()]->substitutePacExpectation(model_name, substexpr));
|
||
assert(substeq);
|
||
equations[eq.value()] = substeq;
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::substitutePacTargetNonstationary(const string& pac_model_name, expr_t substexpr)
|
||
{
|
||
for (auto& eq : equations)
|
||
eq = dynamic_cast<BinaryOpNode*>(
|
||
eq->substitutePacTargetNonstationary(pac_model_name, substexpr));
|
||
}
|
||
|
||
void
|
||
DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder,
|
||
const eval_context_t& eval_context, bool no_tmp_terms, bool block,
|
||
bool use_dll)
|
||
{
|
||
initializeVariablesAndEquations();
|
||
|
||
// Prepare for derivation
|
||
computeDerivIDs();
|
||
|
||
// Computes dynamic jacobian columns, must be done after computeDerivIDs()
|
||
computeDynJacobianCols();
|
||
|
||
/* In both MATLAB and Julia, tensors for higher-order derivatives are stored
|
||
in matrices whose columns correspond to variable multi-indices. Since we
|
||
currently are limited to 32-bit signed integers (hence 31 bits) for matrix
|
||
indices, check that we will not overflow (see #89). Note that such a check
|
||
is not needed for parameter derivatives, since tensors for those are not
|
||
stored as matrices. This check cannot be done before since
|
||
getJacobianColsNbr() is not yet set.
|
||
We only do the check for the legacy representation, since the sparse
|
||
representation is not affected by this problem (TODO: thus the check can be
|
||
removed once the legacy representation is dropped). */
|
||
if (log2(getJacobianColsNbr(false)) * derivsOrder >= numeric_limits<int>::digits)
|
||
{
|
||
cerr << "ERROR: The derivatives matrix of the " << modelClassName()
|
||
<< " is too large. Please decrease the approximation order." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
// Compute derivatives w.r. to all endogenous, exogenous and exogenous deterministic
|
||
set<int> vars;
|
||
for (auto& [symb_lag, deriv_id] : deriv_id_table)
|
||
if (SymbolType type = symbol_table.getType(symb_lag.first); type == SymbolType::endogenous
|
||
|| type == SymbolType::exogenous
|
||
|| type == SymbolType::exogenousDet)
|
||
vars.insert(deriv_id);
|
||
|
||
// Launch computations
|
||
cout << "Computing " << modelClassName() << " derivatives (order " << derivsOrder << ")." << endl;
|
||
|
||
computeDerivatives(derivsOrder, vars);
|
||
|
||
if (derivsOrder > 1)
|
||
for (const auto& [indices, d2] : derivatives[2])
|
||
nonzero_hessian_eqs.insert(indices[0]);
|
||
|
||
if (paramsDerivsOrder > 0)
|
||
{
|
||
cout << "Computing " << modelClassName() << " derivatives w.r.t. parameters (order "
|
||
<< paramsDerivsOrder << ")." << endl;
|
||
computeParamsDerivatives(paramsDerivsOrder);
|
||
}
|
||
|
||
computeTemporaryTerms(!use_dll, no_tmp_terms);
|
||
|
||
if (paramsDerivsOrder > 0 && !no_tmp_terms)
|
||
computeParamsDerivativesTemporaryTerms();
|
||
|
||
/* Change block decomposition settings when the “block” keyword has not been
|
||
used and the model is purely backward/forward or static. This is meant to
|
||
be used with solve_algo={12,14}. */
|
||
if (!block && (max_endo_lag == 0 || max_endo_lead == 0))
|
||
{
|
||
time_recursive_block_decomposition = true;
|
||
mfs = 3; // Typically useful for equations of the form log(x)=RHS
|
||
}
|
||
computingPassBlock(eval_context, no_tmp_terms);
|
||
if (block_decomposed)
|
||
computeBlockDynJacobianCols();
|
||
if (!block_decomposed && block)
|
||
{
|
||
cerr << "ERROR: Block decomposition requested but failed." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeXrefs()
|
||
{
|
||
for (int i {0}; auto& equation : equations)
|
||
{
|
||
ExprNode::EquationInfo ei;
|
||
equation->computeXrefs(ei);
|
||
xrefs[i++] = ei;
|
||
}
|
||
|
||
for (int i {0}; const auto& [eq, eqinfo] : xrefs)
|
||
{
|
||
computeRevXref(xref_param, eqinfo.param, i);
|
||
computeRevXref(xref_endo, eqinfo.endo, i);
|
||
computeRevXref(xref_exo, eqinfo.exo, i);
|
||
computeRevXref(xref_exo_det, eqinfo.exo_det, i);
|
||
i++;
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeRevXref(map<pair<int, int>, set<int>>& xrefset,
|
||
const set<pair<int, int>>& eiref, int eqn)
|
||
{
|
||
for (const auto& it : eiref)
|
||
{
|
||
set<int> eq;
|
||
if (xrefset.contains(it))
|
||
eq = xrefset[it];
|
||
eq.insert(eqn);
|
||
xrefset[it] = eq;
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeXrefs(ostream& output) const
|
||
{
|
||
output << "M_.xref1.param = cell(1, M_.eq_nbr);" << endl
|
||
<< "M_.xref1.endo = cell(1, M_.eq_nbr);" << endl
|
||
<< "M_.xref1.exo = cell(1, M_.eq_nbr);" << endl
|
||
<< "M_.xref1.exo_det = cell(1, M_.eq_nbr);" << endl;
|
||
for (int i {1}; const auto& [eq, eqinfo] : xrefs)
|
||
{
|
||
output << "M_.xref1.param{" << i << "} = [ ";
|
||
for (const auto& [id, lag] : eqinfo.param)
|
||
output << symbol_table.getTypeSpecificID(id) + 1 << " ";
|
||
output << "];" << endl;
|
||
|
||
output << "M_.xref1.endo{" << i << "} = [ ";
|
||
for (const auto& [id, lag] : eqinfo.endo)
|
||
output << "struct('id', " << symbol_table.getTypeSpecificID(id) + 1 << ", 'shift', " << lag
|
||
<< ");";
|
||
output << "];" << endl;
|
||
|
||
output << "M_.xref1.exo{" << i << "} = [ ";
|
||
for (const auto& [id, lag] : eqinfo.exo)
|
||
output << "struct('id', " << symbol_table.getTypeSpecificID(id) + 1 << ", 'shift', " << lag
|
||
<< ");";
|
||
output << "];" << endl;
|
||
|
||
output << "M_.xref1.exo_det{" << i << "} = [ ";
|
||
for (const auto& [id, lag] : eqinfo.exo_det)
|
||
output << "struct('id', " << symbol_table.getTypeSpecificID(id) + 1 << ", 'shift', " << lag
|
||
<< ");";
|
||
output << "];" << endl;
|
||
|
||
i++;
|
||
}
|
||
|
||
output << "M_.xref2.param = cell(1, M_.param_nbr);" << endl
|
||
<< "M_.xref2.endo = cell(1, M_.endo_nbr);" << endl
|
||
<< "M_.xref2.exo = cell(1, M_.exo_nbr);" << endl
|
||
<< "M_.xref2.exo_det = cell(1, M_.exo_det_nbr);" << endl;
|
||
writeRevXrefs(output, xref_param, "param");
|
||
writeRevXrefs(output, xref_endo, "endo");
|
||
writeRevXrefs(output, xref_exo, "exo");
|
||
writeRevXrefs(output, xref_exo_det, "exo_det");
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeRevXrefs(ostream& output, const map<pair<int, int>, set<int>>& xrefmap,
|
||
const string& type) const
|
||
{
|
||
for (int last_tsid {-1}; const auto& [key, eqs] : xrefmap)
|
||
{
|
||
auto& [id, lag] = key;
|
||
int tsid = symbol_table.getTypeSpecificID(id) + 1;
|
||
output << "M_.xref2." << type << "{" << tsid << "} = [ ";
|
||
if (last_tsid == tsid)
|
||
output << "M_.xref2." << type << "{" << tsid << "}; ";
|
||
else
|
||
last_tsid = tsid;
|
||
|
||
for (int eq : eqs)
|
||
if (type == "param")
|
||
output << eq + 1 << " ";
|
||
else
|
||
output << "struct('shift', " << lag << ", 'eq', " << eq + 1 << ");";
|
||
output << "];" << endl;
|
||
}
|
||
}
|
||
|
||
map<tuple<int, int, int>, DynamicModel::BlockDerivativeType>
|
||
DynamicModel::determineBlockDerivativesType(int blk)
|
||
{
|
||
map<tuple<int, int, int>, BlockDerivativeType> derivType;
|
||
int size = blocks[blk].size;
|
||
int nb_recursive = blocks[blk].getRecursiveSize();
|
||
for (int lag {time_recursive_block_decomposition ? 0 : -blocks[blk].max_endo_lag};
|
||
lag <= (time_recursive_block_decomposition ? 0 : blocks[blk].max_endo_lead); lag++)
|
||
for (int eq = 0; eq < size; eq++)
|
||
{
|
||
set<pair<int, int>> endos_and_lags;
|
||
int eq_orig = getBlockEquationID(blk, eq);
|
||
equations[eq_orig]->collectEndogenous(endos_and_lags);
|
||
for (int var = 0; var < size; var++)
|
||
if (int var_orig = getBlockVariableID(blk, var); endos_and_lags.contains({var_orig, lag}))
|
||
{
|
||
if (getBlockEquationType(blk, eq) == EquationType::evaluateRenormalized
|
||
&& eq < nb_recursive)
|
||
/* It’s a normalized recursive equation, we have to recompute
|
||
the derivative using the chain rule */
|
||
derivType[{lag, eq, var}] = BlockDerivativeType::normalizedChainRule;
|
||
else if (!derivType.contains({lag, eq, var}))
|
||
derivType[{lag, eq, var}] = BlockDerivativeType::standard;
|
||
|
||
if (var < nb_recursive)
|
||
for (int feedback_var = nb_recursive; feedback_var < size; feedback_var++)
|
||
if (derivType.contains({lag, var, feedback_var}))
|
||
/* A new derivative needs to be computed using the chain rule
|
||
(a feedback variable appears in the recursive equation
|
||
defining the current variable) */
|
||
derivType[{lag, eq, feedback_var}] = BlockDerivativeType::chainRule;
|
||
}
|
||
}
|
||
return derivType;
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeChainRuleJacobian()
|
||
{
|
||
size_t nb_blocks {blocks.size()};
|
||
|
||
blocks_derivatives.resize(nb_blocks);
|
||
blocks_jacobian_sparse_column_major_order.resize(nb_blocks);
|
||
blocks_jacobian_sparse_colptr.resize(nb_blocks);
|
||
|
||
for (int blk {0}; blk < static_cast<int>(nb_blocks); blk++)
|
||
{
|
||
int nb_recursives = blocks[blk].getRecursiveSize();
|
||
int mfs_size {blocks[blk].mfs_size};
|
||
BlockSimulationType simulation_type {blocks[blk].simulation_type};
|
||
|
||
// Create a map from recursive vars to their defining (normalized) equation
|
||
map<int, BinaryOpNode*> recursive_vars;
|
||
for (int i = 0; i < nb_recursives; i++)
|
||
{
|
||
int deriv_id = getDerivID(
|
||
symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0);
|
||
if (getBlockEquationType(blk, i) == EquationType::evaluateRenormalized)
|
||
recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i);
|
||
else
|
||
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
|
||
}
|
||
|
||
// Compute the block derivatives
|
||
unordered_map<expr_t, set<int>> non_null_chain_rule_derivatives;
|
||
unordered_map<expr_t, map<int, expr_t>> chain_rule_deriv_cache;
|
||
for (const auto& [indices, derivType] : determineBlockDerivativesType(blk))
|
||
{
|
||
auto [lag, eq, var] = indices;
|
||
int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var);
|
||
int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), lag);
|
||
expr_t d {nullptr};
|
||
switch (derivType)
|
||
{
|
||
case BlockDerivativeType::standard:
|
||
if (auto it = derivatives[1].find({eq_orig, deriv_id}); it != derivatives[1].end())
|
||
d = it->second;
|
||
else
|
||
d = Zero;
|
||
break;
|
||
case BlockDerivativeType::chainRule:
|
||
d = equations[eq_orig]->getChainRuleDerivative(deriv_id, recursive_vars,
|
||
non_null_chain_rule_derivatives,
|
||
chain_rule_deriv_cache);
|
||
break;
|
||
case BlockDerivativeType::normalizedChainRule:
|
||
d = equation_type_and_normalized_equation[eq_orig].second->getChainRuleDerivative(
|
||
deriv_id, recursive_vars, non_null_chain_rule_derivatives,
|
||
chain_rule_deriv_cache);
|
||
break;
|
||
}
|
||
|
||
if (d != Zero && eq >= nb_recursives)
|
||
blocks_derivatives[blk][{eq, var, lag}] = d;
|
||
}
|
||
|
||
// Compute the sparse representation of the Jacobian
|
||
if (simulation_type != BlockSimulationType::evaluateForward
|
||
&& simulation_type != BlockSimulationType::evaluateBackward)
|
||
{
|
||
const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
|
||
|| simulation_type == BlockSimulationType::solveForwardSimple
|
||
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
||
|| simulation_type == BlockSimulationType::solveForwardComplete};
|
||
for (const auto& [indices, d1] : blocks_derivatives[blk])
|
||
{
|
||
auto& [eq, var, lag] {indices};
|
||
assert(lag >= -1 && lag <= 1 && eq >= nb_recursives);
|
||
if (var >= nb_recursives && !(one_boundary && lag != 0))
|
||
blocks_jacobian_sparse_column_major_order[blk].try_emplace(
|
||
{eq - nb_recursives,
|
||
var - nb_recursives + static_cast<int>(!one_boundary) * (lag + 1) * mfs_size},
|
||
d1);
|
||
}
|
||
blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(
|
||
blocks_jacobian_sparse_column_major_order[blk], (one_boundary ? 1 : 3) * mfs_size);
|
||
}
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeBlockDynJacobianCols()
|
||
{
|
||
size_t nb_blocks {blocks.size()};
|
||
// Structure used for lexicographic ordering over (lag, var ID)
|
||
vector<set<pair<int, int>>> dynamic_endo(nb_blocks);
|
||
|
||
for (auto& [indices, d1] : derivatives[1])
|
||
{
|
||
auto [eq_orig, deriv_id] {vectorToTuple<2>(indices)};
|
||
int block_eq {eq2block[eq_orig]};
|
||
int var {getTypeSpecificIDByDerivID(deriv_id)};
|
||
int lag {getLagByDerivID(deriv_id)};
|
||
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous && block_eq == endo2block[var])
|
||
dynamic_endo[block_eq].emplace(lag, getBlockInitialVariableID(block_eq, var));
|
||
}
|
||
|
||
// Compute Jacobian column indices
|
||
blocks_jacob_cols_endo.resize(nb_blocks);
|
||
for (size_t blk {0}; blk < nb_blocks; blk++)
|
||
for (int index {0}; auto [lag, var] : dynamic_endo[blk])
|
||
blocks_jacob_cols_endo[blk][{var, lag}] = index++;
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const string& mexext,
|
||
const filesystem::path& matlabroot, bool julia) const
|
||
{
|
||
filesystem::path model_dir {basename};
|
||
model_dir /= "model";
|
||
if (use_dll)
|
||
{
|
||
create_directories(model_dir / "src" / "sparse");
|
||
if (block_decomposed)
|
||
create_directories(model_dir / "src" / "sparse" / "block");
|
||
}
|
||
if (julia)
|
||
create_directories(model_dir / "julia");
|
||
else
|
||
{
|
||
auto plusfolder {packageDir(basename)};
|
||
/* The following is not a duplicate of the same call from
|
||
ModFile::writeMOutput(), because of planner_objective which needs its
|
||
+objective subdirectory */
|
||
create_directories(plusfolder);
|
||
|
||
auto sparsefolder {plusfolder / "+sparse"};
|
||
create_directories(sparsefolder);
|
||
if (block_decomposed)
|
||
create_directories(sparsefolder / "+block");
|
||
|
||
create_directories(plusfolder / "+debug");
|
||
}
|
||
create_directories(model_dir / "bytecode" / "block");
|
||
|
||
// Legacy representation
|
||
if (use_dll)
|
||
writeModelCFile<true>(basename, mexext, matlabroot);
|
||
else if (!julia) // M-files
|
||
writeDynamicMFile(basename);
|
||
// The legacy representation is no longer produced for Julia
|
||
|
||
/* Discretionary optimal policy models don’t have as many variables as
|
||
equations; bytecode does not support that case */
|
||
if (static_cast<int>(equations.size()) == symbol_table.endo_nbr())
|
||
writeDynamicBytecode(basename);
|
||
if (block_decomposed)
|
||
writeDynamicBlockBytecode(basename);
|
||
|
||
// Sparse representation
|
||
if (use_dll)
|
||
writeSparseModelCFiles<true>(basename, mexext, matlabroot);
|
||
else if (julia)
|
||
writeSparseModelJuliaFiles<true>(basename);
|
||
else // MATLAB/Octave
|
||
writeSparseModelMFiles<true>(basename);
|
||
|
||
writeSetAuxiliaryVariablesFile<true>(basename, julia);
|
||
|
||
// Support for model debugging
|
||
if (!julia)
|
||
writeDebugModelMFiles<true>(basename);
|
||
}
|
||
|
||
void
|
||
DynamicModel::clearEquations()
|
||
{
|
||
equations.clear();
|
||
equations_lineno.clear();
|
||
equation_tags.clear();
|
||
}
|
||
|
||
void
|
||
DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const
|
||
{
|
||
dynamic_model.clearEquations();
|
||
|
||
for (size_t i = 0; i < equations.size(); i++)
|
||
dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i]);
|
||
|
||
dynamic_model.equation_tags = equation_tags;
|
||
}
|
||
|
||
int
|
||
DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
|
||
{
|
||
cout << "Ramsey Problem: added " << equations.size() << " multipliers." << endl;
|
||
|
||
// Add aux LM to constraints in equations
|
||
// equation[i]->lhs = rhs becomes equation[i]->MULT_(i+1)*(lhs-rhs) = 0
|
||
for (int i {0}; i < static_cast<int>(equations.size()); i++)
|
||
{
|
||
auto substeq = dynamic_cast<BinaryOpNode*>(equations[i]->addMultipliersToConstraints(i));
|
||
assert(substeq);
|
||
equations[i] = substeq;
|
||
}
|
||
|
||
// Add Planner Objective to equations so that it appears in Lagrangian
|
||
assert(static_model.equations.size() == 1);
|
||
addEquation(static_model.equations[0]->clone(*this), nullopt);
|
||
|
||
// Get max endo lead and max endo lag
|
||
set<pair<int, int>> dynvars;
|
||
int max_eq_lead = 0;
|
||
int max_eq_lag = 0;
|
||
for (auto& equation : equations)
|
||
equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
|
||
|
||
for (const auto& [symb_id, lag] : dynvars)
|
||
{
|
||
max_eq_lead = max(lag, max_eq_lead);
|
||
max_eq_lag = max(-lag, max_eq_lag);
|
||
}
|
||
|
||
// Get Discount Factor
|
||
assert(symbol_table.exists("optimal_policy_discount_factor"));
|
||
int symb_id = symbol_table.getID("optimal_policy_discount_factor");
|
||
assert(symbol_table.getType(symb_id) == SymbolType::parameter);
|
||
expr_t discount_factor_node = AddVariable(symb_id, 0);
|
||
|
||
// Create (modified) Lagrangian (so that we can take the derivative once at time t)
|
||
expr_t lagrangian = Zero;
|
||
for (int i {0}; i < static_cast<int>(equations.size()); i++)
|
||
for (int lag = -max_eq_lag; lag <= max_eq_lead; lag++)
|
||
{
|
||
expr_t dfpower = nullptr;
|
||
stringstream lagstream;
|
||
lagstream << abs(lag);
|
||
if (lag < 0)
|
||
dfpower = AddNonNegativeConstant(lagstream.str());
|
||
else if (lag == 0)
|
||
dfpower = Zero;
|
||
else
|
||
dfpower = AddMinus(Zero, AddNonNegativeConstant(lagstream.str()));
|
||
|
||
lagrangian
|
||
= AddPlus(AddTimes(AddPower(discount_factor_node, dfpower),
|
||
equations[i]->getNonZeroPartofEquation()->decreaseLeadsLags(lag)),
|
||
lagrangian);
|
||
}
|
||
|
||
// Save line numbers and tags, see below
|
||
auto old_equations_lineno = equations_lineno;
|
||
auto old_equation_tags = equation_tags;
|
||
|
||
// Prepare derivation of the Lagrangian
|
||
clearEquations();
|
||
addEquation(AddEqual(lagrangian, Zero), nullopt);
|
||
computeDerivIDs();
|
||
|
||
/* Compute Lagrangian derivatives.
|
||
Also restore line numbers and tags for FOCs w.r.t. a Lagrange multiplier
|
||
(i.e. a FOC identical to an equation of the original model).
|
||
The guarantee given by the SymbolTable class that symbol IDs are
|
||
increasing, plus the fact that derivation IDs are increasing with symbol
|
||
IDs for a given lag, gives the expected ordering of the equations
|
||
(optimality FOCs first, then original equations a.k.a. constraints). */
|
||
vector<expr_t> neweqs;
|
||
vector<optional<int>> neweqs_lineno;
|
||
map<int, map<string, string>> neweqs_tags;
|
||
int orig_endo_nbr {0};
|
||
for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table)
|
||
{
|
||
auto& [symb_id, lag] = symb_id_and_lag;
|
||
if (symbol_table.getType(symb_id) == SymbolType::endogenous && lag == 0)
|
||
{
|
||
neweqs.push_back(
|
||
AddEqual(equations[0]->getNonZeroPartofEquation()->getDerivative(deriv_id), Zero));
|
||
if (optional<int> i = symbol_table.getEquationNumberForMultiplier(symb_id); i)
|
||
{
|
||
// This is a derivative w.r.t. a Lagrange multiplier
|
||
neweqs_lineno.push_back(old_equations_lineno[*i]);
|
||
neweqs_tags[neweqs.size() - 1] = old_equation_tags.getTagsByEqn(*i);
|
||
}
|
||
else
|
||
{
|
||
orig_endo_nbr++;
|
||
neweqs_lineno.emplace_back(nullopt);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Overwrite equations with the Lagrangian derivatives
|
||
clearEquations();
|
||
for (size_t i = 0; i < neweqs.size(); i++)
|
||
addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]);
|
||
|
||
return orig_endo_nbr;
|
||
}
|
||
|
||
bool
|
||
DynamicModel::ParamUsedWithLeadLag() const
|
||
{
|
||
return ParamUsedWithLeadLagInternal();
|
||
}
|
||
|
||
void
|
||
DynamicModel::createVariableMapping()
|
||
{
|
||
for (size_t ii = 0; ii < equations.size(); ii++)
|
||
{
|
||
set<int> eqvars;
|
||
equations[ii]->collectVariables(SymbolType::endogenous, eqvars);
|
||
equations[ii]->collectVariables(SymbolType::exogenous, eqvars);
|
||
for (auto eqvar : eqvars)
|
||
variableMapping[symbol_table.getUltimateOrigSymbID(eqvar)].emplace(ii);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::expandEqTags()
|
||
{
|
||
set<int> existing_tags = equation_tags.getEqnsByKey("name");
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
if (!existing_tags.contains(eq))
|
||
{
|
||
if (auto lhs_expr = dynamic_cast<VariableNode*>(equations[eq]->arg1);
|
||
lhs_expr && !equation_tags.exists("name", symbol_table.getName(lhs_expr->symb_id)))
|
||
equation_tags.add(eq, "name", symbol_table.getName(lhs_expr->symb_id));
|
||
else if (!equation_tags.exists("name", to_string(eq + 1)))
|
||
equation_tags.add(eq, "name", to_string(eq + 1));
|
||
else
|
||
{
|
||
cerr << "Error creating default equation tag: cannot assign default tag to equation "
|
||
"number "
|
||
<< eq + 1 << " because it is already in use" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
}
|
||
|
||
set<int>
|
||
DynamicModel::findUnusedEndogenous()
|
||
{
|
||
set<int> usedEndo, unusedEndo;
|
||
for (auto& equation : equations)
|
||
equation->collectVariables(SymbolType::endogenous, usedEndo);
|
||
for (auto& equation : static_only_equations)
|
||
equation->collectVariables(SymbolType::endogenous, usedEndo);
|
||
set<int> allEndo = symbol_table.getEndogenous();
|
||
set_difference(allEndo.begin(), allEndo.end(), usedEndo.begin(), usedEndo.end(),
|
||
inserter(unusedEndo, unusedEndo.begin()));
|
||
return unusedEndo;
|
||
}
|
||
|
||
set<int>
|
||
DynamicModel::findUnusedExogenous()
|
||
{
|
||
set<int> usedExo, unusedExo, unobservedExo;
|
||
for (auto& equation : equations)
|
||
equation->collectVariables(SymbolType::exogenous, usedExo);
|
||
for (auto& equation : static_only_equations)
|
||
equation->collectVariables(SymbolType::exogenous, usedExo);
|
||
set<int> observedExo = symbol_table.getObservedExogenous();
|
||
set<int> allExo = symbol_table.getExogenous();
|
||
set_difference(allExo.begin(), allExo.end(), observedExo.begin(), observedExo.end(),
|
||
inserter(unobservedExo, unobservedExo.begin()));
|
||
set_difference(unobservedExo.begin(), unobservedExo.end(), usedExo.begin(), usedExo.end(),
|
||
inserter(unusedExo, unusedExo.begin()));
|
||
return unusedExo;
|
||
}
|
||
|
||
void
|
||
DynamicModel::setLeadsLagsOrig()
|
||
{
|
||
set<pair<int, int>> dynvars;
|
||
|
||
for (auto& equation : equations)
|
||
{
|
||
equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::exogenousDet, dynvars);
|
||
|
||
max_lag_with_diffs_expanded_orig
|
||
= max(equation->maxLagWithDiffsExpanded(), max_lag_with_diffs_expanded_orig);
|
||
}
|
||
|
||
for (const auto& [symb_id, lag] : dynvars)
|
||
{
|
||
SymbolType type = symbol_table.getType(symb_id);
|
||
|
||
max_lead_orig = max(lag, max_lead_orig);
|
||
max_lag_orig = max(-lag, max_lag_orig);
|
||
|
||
switch (type)
|
||
{
|
||
case SymbolType::endogenous:
|
||
max_endo_lead_orig = max(lag, max_endo_lead_orig);
|
||
max_endo_lag_orig = max(-lag, max_endo_lag_orig);
|
||
break;
|
||
case SymbolType::exogenous:
|
||
max_exo_lead_orig = max(lag, max_exo_lead_orig);
|
||
max_exo_lag_orig = max(-lag, max_exo_lag_orig);
|
||
break;
|
||
case SymbolType::exogenousDet:
|
||
max_exo_det_lead_orig = max(lag, max_exo_det_lead_orig);
|
||
max_exo_det_lag_orig = max(-lag, max_exo_det_lag_orig);
|
||
break;
|
||
default:
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeDerivIDs()
|
||
{
|
||
set<pair<int, int>> dynvars;
|
||
|
||
for (auto& equation : equations)
|
||
{
|
||
equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::exogenousDet, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::parameter, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::trend, dynvars);
|
||
equation->collectDynamicVariables(SymbolType::logTrend, dynvars);
|
||
}
|
||
|
||
for (const auto& [symb_id, lag] : dynvars)
|
||
{
|
||
SymbolType type = symbol_table.getType(symb_id);
|
||
|
||
/* Setting maximum and minimum lags.
|
||
|
||
We don't want these to be affected by lead/lags on parameters: they
|
||
are accepted for facilitating variable flipping, but are simply
|
||
ignored. */
|
||
if (type != SymbolType::parameter)
|
||
{
|
||
max_lead = max(lag, max_lead);
|
||
max_lag = max(-lag, max_lag);
|
||
}
|
||
|
||
switch (type)
|
||
{
|
||
case SymbolType::endogenous:
|
||
max_endo_lead = max(lag, max_endo_lead);
|
||
max_endo_lag = max(-lag, max_endo_lag);
|
||
break;
|
||
case SymbolType::exogenous:
|
||
max_exo_lead = max(lag, max_exo_lead);
|
||
max_exo_lag = max(-lag, max_exo_lag);
|
||
break;
|
||
case SymbolType::exogenousDet:
|
||
max_exo_det_lead = max(lag, max_exo_det_lead);
|
||
max_exo_det_lag = max(-lag, max_exo_det_lag);
|
||
break;
|
||
default:
|
||
break;
|
||
}
|
||
|
||
// Create a new deriv_id
|
||
int deriv_id = deriv_id_table.size();
|
||
|
||
deriv_id_table[{symb_id, lag}] = deriv_id;
|
||
inv_deriv_id_table.emplace_back(symb_id, lag);
|
||
}
|
||
}
|
||
|
||
SymbolType
|
||
DynamicModel::getTypeByDerivID(int deriv_id) const noexcept(false)
|
||
{
|
||
return symbol_table.getType(getSymbIDByDerivID(deriv_id));
|
||
}
|
||
|
||
int
|
||
DynamicModel::getLagByDerivID(int deriv_id) const noexcept(false)
|
||
{
|
||
if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
|
||
throw UnknownDerivIDException();
|
||
|
||
return inv_deriv_id_table[deriv_id].second;
|
||
}
|
||
|
||
int
|
||
DynamicModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
|
||
{
|
||
if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
|
||
throw UnknownDerivIDException();
|
||
|
||
return inv_deriv_id_table[deriv_id].first;
|
||
}
|
||
|
||
int
|
||
DynamicModel::getTypeSpecificIDByDerivID(int deriv_id) const
|
||
{
|
||
return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
|
||
}
|
||
|
||
int
|
||
DynamicModel::getDerivID(int symb_id, int lag) const noexcept(false)
|
||
{
|
||
if (auto it = deriv_id_table.find({symb_id, lag}); it == deriv_id_table.end())
|
||
throw UnknownDerivIDException();
|
||
else
|
||
return it->second;
|
||
}
|
||
|
||
void
|
||
DynamicModel::addAllParamDerivId(set<int>& deriv_id_set)
|
||
{
|
||
for (size_t i = 0; i < inv_deriv_id_table.size(); i++)
|
||
if (symbol_table.getType(inv_deriv_id_table[i].first) == SymbolType::parameter)
|
||
deriv_id_set.insert(i);
|
||
}
|
||
|
||
void
|
||
DynamicModel::computeDynJacobianCols()
|
||
{
|
||
// Sort the dynamic endogenous variables by lexicographic order over (lag,
|
||
// type_specific_symbol_id)
|
||
map<pair<int, int>, int> ordered_dyn_endo;
|
||
for (const auto& [symb_lag, deriv_id] : deriv_id_table)
|
||
if (const auto& [symb_id, lag] = symb_lag;
|
||
symbol_table.getType(symb_id) == SymbolType::endogenous)
|
||
ordered_dyn_endo[{lag, symbol_table.getTypeSpecificID(symb_id)}] = deriv_id;
|
||
|
||
// Fill the dynamic jacobian columns for endogenous (legacy representation)
|
||
for (int sorted_id {0}; const auto& [ignore, deriv_id] : ordered_dyn_endo)
|
||
dyn_jacobian_cols_table[deriv_id] = sorted_id++;
|
||
|
||
/* Fill the dynamic columns for exogenous and exogenous deterministic (legacy
|
||
representation) */
|
||
for (const auto& [symb_lag, deriv_id] : deriv_id_table)
|
||
{
|
||
int symb_id {symb_lag.first};
|
||
int tsid {symbol_table.getTypeSpecificID(symb_id)}; // At this point, there is no trend_var
|
||
if (SymbolType type {symbol_table.getType(symb_id)}; type == SymbolType::exogenous)
|
||
dyn_jacobian_cols_table[deriv_id] = ordered_dyn_endo.size() + tsid;
|
||
else if (type == SymbolType::exogenousDet)
|
||
dyn_jacobian_cols_table[deriv_id] = ordered_dyn_endo.size() + symbol_table.exo_nbr() + tsid;
|
||
}
|
||
|
||
/* NB: the following could differ from dyn_jacobian_cols_table.size() if
|
||
there are unused exogenous (and “nostrict” option is given) */
|
||
dyn_jacobian_ncols
|
||
= ordered_dyn_endo.size() + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
|
||
}
|
||
|
||
void
|
||
DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t& eval_context)
|
||
{
|
||
for (auto& [symb_lag1, deriv_id1] : deriv_id_table)
|
||
if (auto& [symb_id1, lag1] = symb_lag1;
|
||
symbol_table.getType(symb_id1) == SymbolType::trend
|
||
|| symbol_table.getType(symb_id1) == SymbolType::logTrend)
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
expr_t homogeneq = AddMinus(equations[eq]->arg1, equations[eq]->arg2);
|
||
|
||
// Do not run the test if the term inside the log is zero
|
||
if (fabs(homogeneq->eval(eval_context)) > zero_band)
|
||
{
|
||
expr_t testeq = AddLog(homogeneq); // F = log(lhs-rhs)
|
||
testeq = testeq->getDerivative(deriv_id1); // d F / d Trend
|
||
for (auto& [symb_lag2, deriv_id2] : deriv_id_table)
|
||
if (auto& [symb_id2, lag2] = symb_lag2;
|
||
symbol_table.getType(symb_id2) == SymbolType::endogenous)
|
||
{
|
||
double nearZero = testeq->getDerivative(deriv_id2)->eval(
|
||
eval_context); // eval d F / d Trend d Endog
|
||
if (fabs(nearZero) > balanced_growth_test_tol)
|
||
{
|
||
cerr << "ERROR: trends not compatible with balanced growth path; the "
|
||
"second-order cross partial of equation "
|
||
<< eq + 1;
|
||
if (equations_lineno[eq])
|
||
cerr << " (line " << *equations_lineno[eq] << ") ";
|
||
cerr << "w.r.t. trend variable " << symbol_table.getName(symb_id1)
|
||
<< " and endogenous variable " << symbol_table.getName(symb_id2)
|
||
<< " is not null (abs. value = " << fabs(nearZero)
|
||
<< "). If you are confident that your trends are correctly specified, "
|
||
"you can raise the value of option 'balanced_growth_test_tol' in "
|
||
"the 'model' block."
|
||
<< endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeLatexFile(const string& basename, bool write_equation_tags) const
|
||
{
|
||
writeLatexModelFile(basename, "dynamic", ExprNodeOutputType::latexDynamicModel,
|
||
write_equation_tags);
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeLatexOriginalFile(const string& basename, bool write_equation_tags) const
|
||
{
|
||
writeLatexModelFile(basename, "original", ExprNodeOutputType::latexDynamicModel,
|
||
write_equation_tags);
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteEndoLeadGreaterThanTwo(bool deterministic_model)
|
||
{
|
||
substituteLeadLagInternal(AuxVarType::endoLead, deterministic_model, {});
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteEndoLagGreaterThanTwo(bool deterministic_model)
|
||
{
|
||
substituteLeadLagInternal(AuxVarType::endoLag, deterministic_model, {});
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteExoLead(bool deterministic_model)
|
||
{
|
||
substituteLeadLagInternal(AuxVarType::exoLead, deterministic_model, {});
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteExoLag(bool deterministic_model)
|
||
{
|
||
substituteLeadLagInternal(AuxVarType::exoLag, deterministic_model, {});
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteLeadLagInternal(AuxVarType type, bool deterministic_model,
|
||
const vector<string>& subset)
|
||
{
|
||
ExprNode::subst_table_t subst_table;
|
||
vector<BinaryOpNode*> neweqs;
|
||
|
||
// Substitute in used model local variables
|
||
set<int> used_local_vars;
|
||
for (auto& equation : equations)
|
||
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
||
|
||
for (int used_local_var : used_local_vars)
|
||
{
|
||
const expr_t value = local_variables_table.at(used_local_var);
|
||
expr_t subst;
|
||
switch (type)
|
||
{
|
||
case AuxVarType::endoLead:
|
||
subst = value->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
|
||
break;
|
||
case AuxVarType::endoLag:
|
||
subst = value->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
|
||
break;
|
||
case AuxVarType::exoLead:
|
||
subst = value->substituteExoLead(subst_table, neweqs, deterministic_model);
|
||
break;
|
||
case AuxVarType::exoLag:
|
||
subst = value->substituteExoLag(subst_table, neweqs);
|
||
break;
|
||
case AuxVarType::diffForward:
|
||
subst = value->differentiateForwardVars(subset, subst_table, neweqs);
|
||
break;
|
||
default:
|
||
cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
local_variables_table[used_local_var] = subst;
|
||
}
|
||
|
||
// Substitute in equations
|
||
for (auto& equation : equations)
|
||
{
|
||
expr_t subst;
|
||
switch (type)
|
||
{
|
||
case AuxVarType::endoLead:
|
||
subst = equation->substituteEndoLeadGreaterThanTwo(subst_table, neweqs,
|
||
deterministic_model);
|
||
break;
|
||
case AuxVarType::endoLag:
|
||
subst = equation->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
|
||
break;
|
||
case AuxVarType::exoLead:
|
||
subst = equation->substituteExoLead(subst_table, neweqs, deterministic_model);
|
||
break;
|
||
case AuxVarType::exoLag:
|
||
subst = equation->substituteExoLag(subst_table, neweqs);
|
||
break;
|
||
case AuxVarType::diffForward:
|
||
subst = equation->differentiateForwardVars(subset, subst_table, neweqs);
|
||
break;
|
||
default:
|
||
cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
auto substeq = dynamic_cast<BinaryOpNode*>(subst);
|
||
assert(substeq);
|
||
equation = substeq;
|
||
}
|
||
|
||
// Add new equations
|
||
for (auto& neweq : neweqs)
|
||
{
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
}
|
||
|
||
if (neweqs.size() > 0)
|
||
{
|
||
cout << "Substitution of ";
|
||
switch (type)
|
||
{
|
||
case AuxVarType::endoLead:
|
||
cout << "endo leads >= 2";
|
||
break;
|
||
case AuxVarType::endoLag:
|
||
cout << "endo lags >= 2";
|
||
break;
|
||
case AuxVarType::exoLead:
|
||
cout << "exo leads";
|
||
break;
|
||
case AuxVarType::exoLag:
|
||
cout << "exo lags";
|
||
break;
|
||
case AuxVarType::expectation:
|
||
cout << "expectation";
|
||
break;
|
||
case AuxVarType::diffForward:
|
||
cout << "forward vars";
|
||
break;
|
||
default:
|
||
cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
cout << ": added " << neweqs.size() << " auxiliary variables and equations." << endl;
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteAdl()
|
||
{
|
||
/* Contrary to other substitution methods, we do the substitution in MLV
|
||
definitions here, instead of doing it at the ExprNode method level,
|
||
because otherwise this would substitute MLV in the original model (see
|
||
#65). */
|
||
for (auto& [id, definition] : local_variables_table)
|
||
definition = definition->substituteAdl();
|
||
|
||
for (auto& equation : equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->substituteAdl());
|
||
|
||
for (auto& equation : static_only_equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->substituteAdl());
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteModelLocalVariables()
|
||
{
|
||
for (auto& equation : equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->substituteModelLocalVariables());
|
||
|
||
for (auto& equation : static_only_equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->substituteModelLocalVariables());
|
||
|
||
/* We can’t clear local_variables_table at this point, because in case of
|
||
ramsey_policy, the original model is saved via DynamicModel::operator=()
|
||
before computing the FOC. But since DataTree::operator=() clones all
|
||
nodes, it will try to clone nodes containing model-local variables, and
|
||
this will fail at the point where DataTree methods try to evaluate those
|
||
nodes to a numerical value. */
|
||
}
|
||
|
||
set<int>
|
||
DynamicModel::getEquationNumbersFromTags(const set<string>& eqtags) const
|
||
{
|
||
set<int> eqnumbers;
|
||
for (auto& eqtag : eqtags)
|
||
{
|
||
set<int> tmp = equation_tags.getEqnsByTag("name", eqtag);
|
||
if (tmp.empty())
|
||
{
|
||
cerr << "ERROR: looking for equation tag " << eqtag << " failed." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
eqnumbers.insert(tmp.begin(), tmp.end());
|
||
}
|
||
return eqnumbers;
|
||
}
|
||
|
||
set<int>
|
||
DynamicModel::findPacExpectationEquationNumbers() const
|
||
{
|
||
set<int> eqnumbers;
|
||
for (int i {0}; auto& equation : equations)
|
||
{
|
||
if (equation->containsPacExpectation())
|
||
eqnumbers.insert(i);
|
||
i++;
|
||
}
|
||
return eqnumbers;
|
||
}
|
||
|
||
pair<lag_equivalence_table_t, ExprNode::subst_table_t>
|
||
DynamicModel::substituteUnaryOps(VarExpectationModelTable& var_expectation_model_table,
|
||
PacModelTable& pac_model_table)
|
||
{
|
||
vector<int> eqnumbers(equations.size());
|
||
iota(eqnumbers.begin(), eqnumbers.end(), 0);
|
||
return substituteUnaryOps({eqnumbers.begin(), eqnumbers.end()}, var_expectation_model_table,
|
||
pac_model_table);
|
||
}
|
||
|
||
pair<lag_equivalence_table_t, ExprNode::subst_table_t>
|
||
DynamicModel::substituteUnaryOps(const set<int>& eqnumbers,
|
||
VarExpectationModelTable& var_expectation_model_table,
|
||
PacModelTable& pac_model_table)
|
||
{
|
||
lag_equivalence_table_t nodes;
|
||
ExprNode::subst_table_t subst_table;
|
||
|
||
// Mark unary ops to be substituted in model local variables that appear in selected equations
|
||
set<int> used_local_vars;
|
||
for (int eqnumber : eqnumbers)
|
||
equations[eqnumber]->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
||
for (int mlv : used_local_vars)
|
||
local_variables_table[mlv]->findUnaryOpNodesForAuxVarCreation(nodes);
|
||
|
||
// Mark unary ops to be substituted in selected equations
|
||
for (int eqnumber : eqnumbers)
|
||
equations[eqnumber]->findUnaryOpNodesForAuxVarCreation(nodes);
|
||
|
||
// Substitute in model local variables
|
||
vector<BinaryOpNode*> neweqs;
|
||
for (int mlv : used_local_vars)
|
||
local_variables_table[mlv]
|
||
= local_variables_table[mlv]->substituteUnaryOpNodes(nodes, subst_table, neweqs);
|
||
|
||
// Substitute in equations
|
||
for (int eq : eqnumbers)
|
||
{
|
||
auto substeq = dynamic_cast<BinaryOpNode*>(
|
||
equations[eq]->substituteUnaryOpNodes(nodes, subst_table, neweqs));
|
||
assert(substeq);
|
||
equations[eq] = substeq;
|
||
}
|
||
|
||
// Substitute in expressions of var_expectation_model
|
||
var_expectation_model_table.substituteUnaryOpsInExpression(nodes, subst_table, neweqs);
|
||
// Substitute in growth terms in pac_model and pac_target_info
|
||
pac_model_table.substituteUnaryOpsInGrowth(nodes, subst_table, neweqs);
|
||
|
||
// Add new equations
|
||
for (auto& neweq : neweqs)
|
||
{
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
}
|
||
|
||
if (subst_table.size() > 0)
|
||
cout << "Substitution of Unary Ops: added " << neweqs.size()
|
||
<< " auxiliary variables and equations." << endl;
|
||
|
||
return {nodes, subst_table};
|
||
}
|
||
|
||
pair<lag_equivalence_table_t, ExprNode::subst_table_t>
|
||
DynamicModel::substituteDiff(VarExpectationModelTable& var_expectation_model_table,
|
||
PacModelTable& pac_model_table)
|
||
{
|
||
/* Note: at this point, we know that there is no diff operator with a lead,
|
||
because they have been expanded by DataTree::AddDiff().
|
||
Hence we can go forward with the substitution without worrying about the
|
||
expectation operator. */
|
||
|
||
lag_equivalence_table_t diff_nodes;
|
||
ExprNode::subst_table_t diff_subst_table;
|
||
|
||
// Mark diff operators to be substituted in model local variables
|
||
set<int> used_local_vars;
|
||
for (auto equation : equations)
|
||
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
||
for (auto& [symb_id, expr] : local_variables_table)
|
||
if (used_local_vars.contains(symb_id))
|
||
expr->findDiffNodes(diff_nodes);
|
||
|
||
// Mark diff operators to be substituted in equations
|
||
for (auto equation : equations)
|
||
equation->findDiffNodes(diff_nodes);
|
||
|
||
pac_model_table.findDiffNodesInGrowth(diff_nodes);
|
||
|
||
// Substitute in model local variables
|
||
vector<BinaryOpNode*> neweqs;
|
||
for (auto& [symb_id, expr] : local_variables_table)
|
||
expr = expr->substituteDiff(diff_nodes, diff_subst_table, neweqs);
|
||
|
||
// Substitute in equations
|
||
for (auto& equation : equations)
|
||
{
|
||
auto substeq = dynamic_cast<BinaryOpNode*>(
|
||
equation->substituteDiff(diff_nodes, diff_subst_table, neweqs));
|
||
assert(substeq);
|
||
equation = substeq;
|
||
}
|
||
|
||
var_expectation_model_table.substituteDiffNodesInExpression(diff_nodes, diff_subst_table, neweqs);
|
||
pac_model_table.substituteDiffNodesInGrowth(diff_nodes, diff_subst_table, neweqs);
|
||
|
||
// Add new equations
|
||
for (auto neweq : neweqs)
|
||
{
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
}
|
||
|
||
if (diff_subst_table.size() > 0)
|
||
cout << "Substitution of Diff operator: added " << neweqs.size()
|
||
<< " auxiliary variables and equations." << endl;
|
||
|
||
return {diff_nodes, diff_subst_table};
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteExpectation(bool partial_information_model)
|
||
{
|
||
ExprNode::subst_table_t subst_table;
|
||
vector<BinaryOpNode*> neweqs;
|
||
|
||
// Substitute in model local variables
|
||
for (auto& [symb_id, expr] : local_variables_table)
|
||
expr = expr->substituteExpectation(subst_table, neweqs, partial_information_model);
|
||
|
||
// Substitute in equations
|
||
for (auto& equation : equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(
|
||
equation->substituteExpectation(subst_table, neweqs, partial_information_model));
|
||
assert(equation);
|
||
}
|
||
|
||
/* No need to substitute in static_only_equations, since expectation()
|
||
operators in [static] equations are forbidden at the parsing level. */
|
||
|
||
// Add new equations
|
||
for (auto neweq : neweqs)
|
||
{
|
||
addEquation(neweq, nullopt);
|
||
addAuxEquation(neweq);
|
||
}
|
||
|
||
if (subst_table.size() > 0)
|
||
{
|
||
if (partial_information_model)
|
||
cout << "Substitution of Expectation operator: added " << subst_table.size()
|
||
<< " auxiliary variables and " << neweqs.size() << " auxiliary equations." << endl;
|
||
else
|
||
cout << "Substitution of Expectation operator: added " << neweqs.size()
|
||
<< " auxiliary variables and equations." << endl;
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::transformPredeterminedVariables()
|
||
{
|
||
for (auto& [id, definition] : local_variables_table)
|
||
definition = definition->decreaseLeadsLagsPredeterminedVariables();
|
||
|
||
for (auto& equation : equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->decreaseLeadsLagsPredeterminedVariables());
|
||
assert(equation);
|
||
}
|
||
|
||
// No need to handle static_only_equations, since there are no leads/lags there
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteLogTransform()
|
||
{
|
||
for (int symb_id : symbol_table.getVariablesWithLogTransform())
|
||
{
|
||
expr_t aux_def = AddLog(AddVariable(symb_id));
|
||
int aux_symb_id = symbol_table.addLogTransformAuxiliaryVar(symb_id, 0, aux_def);
|
||
|
||
for (auto& [id, definition] : local_variables_table)
|
||
definition = definition->substituteLogTransform(symb_id, aux_symb_id);
|
||
|
||
for (auto& equation : equations)
|
||
equation
|
||
= dynamic_cast<BinaryOpNode*>(equation->substituteLogTransform(symb_id, aux_symb_id));
|
||
|
||
for (auto& equation : static_only_equations)
|
||
equation
|
||
= dynamic_cast<BinaryOpNode*>(equation->substituteLogTransform(symb_id, aux_symb_id));
|
||
|
||
/*
|
||
We add the following new equations:
|
||
+ X=exp(log_X) to the model
|
||
+ log_X=log(X) to the list of auxiliary equations
|
||
|
||
In this way:
|
||
+ statements like X=1 in initval/endval blocks will be correctly
|
||
handled (i.e. log_X will be initialized to 0 in this case), through
|
||
the set_auxiliary_variables.m and dynamic_set_auxiliary_series.m files
|
||
+ computation of X in perfect foresight simulations will be done by
|
||
simple evaluation when using block decomposition (X will belong to an
|
||
block of type “evaluate”, or maybe even the epilogue)
|
||
*/
|
||
addAuxEquation(AddEqual(AddVariable(aux_symb_id), aux_def));
|
||
addEquation(AddEqual(AddVariable(symb_id), AddExp(AddVariable(aux_symb_id))), nullopt, {});
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
|
||
{
|
||
set<int> endos;
|
||
for (int eq : eqnumbers)
|
||
equations[eq]->collectVariables(SymbolType::endogenous, endos);
|
||
|
||
const set<int>& with_log_transform = symbol_table.getVariablesWithLogTransform();
|
||
|
||
vector<int> intersect;
|
||
set_intersection(endos.begin(), endos.end(), with_log_transform.begin(), with_log_transform.end(),
|
||
back_inserter(intersect));
|
||
if (!intersect.empty())
|
||
{
|
||
cerr << "ERROR: the following variables are declared with var(log) and therefore cannot "
|
||
"appear in a VAR/TCM/PAC equation: ";
|
||
for (int symb_id : intersect)
|
||
cerr << symbol_table.getName(symb_id) << " ";
|
||
cerr << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::detrendEquations()
|
||
{
|
||
// We go backwards in the list of trend_vars, to deal correctly with I(2) processes
|
||
for (const auto& it : std::ranges::reverse_view(nonstationary_symbols_map))
|
||
{
|
||
for (auto& equation : equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(
|
||
equation->detrend(it.first, it.second.first, it.second.second));
|
||
assert(equation);
|
||
}
|
||
for (auto& equation : static_only_equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(
|
||
equation->detrend(it.first, it.second.first, it.second.second));
|
||
assert(equation);
|
||
}
|
||
}
|
||
|
||
for (auto& equation : equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->removeTrendLeadLag(trend_symbols_map));
|
||
assert(equation);
|
||
}
|
||
for (auto& equation : static_only_equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->removeTrendLeadLag(trend_symbols_map));
|
||
assert(equation);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::removeTrendVariableFromEquations()
|
||
{
|
||
for (auto& equation : equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->replaceTrendVar());
|
||
assert(equation);
|
||
}
|
||
for (auto& equation : static_only_equations)
|
||
{
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->replaceTrendVar());
|
||
assert(equation);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::differentiateForwardVars(const vector<string>& subset)
|
||
{
|
||
substituteLeadLagInternal(AuxVarType::diffForward, true, subset);
|
||
}
|
||
|
||
void
|
||
DynamicModel::fillEvalContext(eval_context_t& eval_context) const
|
||
{
|
||
// First, auxiliary variables
|
||
for (auto aux_equation : aux_equations)
|
||
{
|
||
assert(aux_equation->op_code == BinaryOpcode::equal);
|
||
auto auxvar = dynamic_cast<VariableNode*>(aux_equation->arg1);
|
||
assert(auxvar);
|
||
try
|
||
{
|
||
double val = aux_equation->arg2->eval(eval_context);
|
||
eval_context[auxvar->symb_id] = val;
|
||
}
|
||
catch (ExprNode::EvalException& e)
|
||
{
|
||
// Do nothing
|
||
}
|
||
}
|
||
|
||
// Second, model local variables
|
||
for (auto& [symb_id, expression] : local_variables_table)
|
||
{
|
||
try
|
||
{
|
||
double val = expression->eval(eval_context);
|
||
eval_context[symb_id] = val;
|
||
}
|
||
catch (ExprNode::EvalException& e)
|
||
{
|
||
// Do nothing
|
||
}
|
||
}
|
||
|
||
// Third, trend variables
|
||
for (int trendVar : symbol_table.getTrendVarIds())
|
||
eval_context[trendVar] = 2; // not <= 0 bc of log, not 1 bc of powers
|
||
}
|
||
|
||
void
|
||
DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
|
||
map<string, string> eq_tags)
|
||
{
|
||
auto beq = dynamic_cast<BinaryOpNode*>(eq);
|
||
assert(beq && beq->op_code == BinaryOpcode::equal);
|
||
|
||
static_only_equations_equation_tags.add(static_only_equations.size(), move(eq_tags));
|
||
static_only_equations.push_back(beq);
|
||
static_only_equations_lineno.push_back(lineno);
|
||
}
|
||
|
||
size_t
|
||
DynamicModel::staticOnlyEquationsNbr() const
|
||
{
|
||
return static_only_equations.size();
|
||
}
|
||
|
||
size_t
|
||
DynamicModel::dynamicOnlyEquationsNbr() const
|
||
{
|
||
return equation_tags.getDynamicEqns().size();
|
||
}
|
||
|
||
void
|
||
DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags,
|
||
const vector<string>& regimes_bind,
|
||
const vector<string>& regimes_relax)
|
||
{
|
||
auto beq = dynamic_cast<BinaryOpNode*>(eq);
|
||
assert(beq && beq->op_code == BinaryOpcode::equal);
|
||
|
||
// Construct the term to be added to the corresponding equation
|
||
expr_t basic_term = AddMinus(beq->arg1, beq->arg2);
|
||
expr_t term = basic_term;
|
||
for (auto& regime : regimes_bind)
|
||
{
|
||
int param_id = symbol_table.getID(ParsingDriver::buildOccbinBindParamName(regime));
|
||
term = AddTimes(term, AddVariable(param_id));
|
||
}
|
||
for (auto& regime : regimes_relax)
|
||
{
|
||
int param_id = symbol_table.getID(ParsingDriver::buildOccbinBindParamName(regime));
|
||
term = AddTimes(term, AddMinus(One, AddVariable(param_id)));
|
||
}
|
||
|
||
// Create or update the dynamic equation
|
||
optional<int> eqn {equation_tags.getEqnByTag("name", eq_tags.at("name"))};
|
||
if (eqn)
|
||
{
|
||
BinaryOpNode* orig_eq {equations[*eqn]};
|
||
/* In the following, we could have kept only orig_eq->arg1, but the
|
||
following adds a (somewhat bizarre) support for equation snippets
|
||
without “bind” nor “relax” */
|
||
equations[*eqn] = AddEqual(AddPlus(AddMinus(orig_eq->arg1, orig_eq->arg2), term), Zero);
|
||
// It’s unclear how to update lineno and tags, so don’t do it
|
||
}
|
||
else
|
||
{
|
||
auto eq_tags_dynamic = eq_tags;
|
||
eq_tags_dynamic["dynamic"] = "";
|
||
addEquation(AddEqual(term, Zero), lineno, eq_tags_dynamic);
|
||
}
|
||
|
||
// Create or update the static equation (corresponding to the pure relax regime)
|
||
if (regimes_bind.empty())
|
||
{
|
||
/* Similar remark as above. We could have entirely skipped this
|
||
equation updating, since normally there is only one such clause,
|
||
but the following adds a (somewhat bizarre) support for equation
|
||
snippets without “bind” nor “relax” */
|
||
optional<int> eqn_static {
|
||
static_only_equations_equation_tags.getEqnByTag("name", eq_tags.at("name"))};
|
||
if (eqn_static)
|
||
{
|
||
BinaryOpNode* orig_eq {static_only_equations[*eqn_static]};
|
||
static_only_equations[*eqn_static]
|
||
= AddEqual(AddPlus(AddMinus(orig_eq->arg1, orig_eq->arg2), basic_term), Zero);
|
||
// It’s unclear how to update lineno and tags, so don’t do it
|
||
}
|
||
else
|
||
{
|
||
eq_tags["static"] = "";
|
||
addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, move(eq_tags));
|
||
}
|
||
}
|
||
}
|
||
|
||
bool
|
||
DynamicModel::isChecksumMatching(const string& basename) const
|
||
{
|
||
stringstream buffer;
|
||
|
||
// Write equation tags
|
||
equation_tags.writeCheckSumInfo(buffer);
|
||
|
||
constexpr ExprNodeOutputType buffer_type {ExprNodeOutputType::CDynamicModel};
|
||
|
||
deriv_node_temp_terms_t tef_terms;
|
||
temporary_terms_t temp_term_union;
|
||
|
||
writeTemporaryTerms<buffer_type>(temporary_terms_derivatives[0], temp_term_union,
|
||
temporary_terms_idxs, buffer, tef_terms);
|
||
|
||
writeModelEquations<buffer_type>(buffer, temp_term_union);
|
||
|
||
size_t result = hash<string> {}(buffer.str());
|
||
|
||
// check whether basename directory exist. If not, create it.
|
||
// If it does, read old checksum if it exists, return if equal to result
|
||
fstream checksum_file;
|
||
auto filename = filesystem::path {basename} / "checksum";
|
||
if (!filesystem::create_directory(basename))
|
||
{
|
||
checksum_file.open(filename, ios::in | ios::binary);
|
||
if (checksum_file.is_open())
|
||
{
|
||
size_t old_checksum;
|
||
checksum_file >> old_checksum;
|
||
checksum_file.close();
|
||
if (old_checksum == result)
|
||
return true;
|
||
}
|
||
}
|
||
|
||
// write new checksum file if none or different from old checksum
|
||
checksum_file.open(filename, ios::out | ios::binary);
|
||
if (!checksum_file.is_open())
|
||
{
|
||
cerr << "ERROR: Can't open file " << filename.string() << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
checksum_file << result;
|
||
checksum_file.close();
|
||
return false;
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonOutput(ostream& output) const
|
||
{
|
||
deriv_node_temp_terms_t tef_terms;
|
||
writeJsonModelLocalVariables(output, false, tef_terms);
|
||
output << ", ";
|
||
writeJsonModelEquations(output, false);
|
||
output << ", ";
|
||
writeJsonXrefs(output);
|
||
output << ", ";
|
||
writeJsonAST(output);
|
||
output << ", ";
|
||
writeJsonVariableMapping(output);
|
||
output << R"(, "dynamic_tmp_nbr": [)";
|
||
for (bool printed_something {false}; const auto& tts : temporary_terms_derivatives)
|
||
{
|
||
if (exchange(printed_something, true))
|
||
output << ", ";
|
||
output << tts.size();
|
||
}
|
||
output << "], ";
|
||
writeJsonSparseIndicesHelper<true>(output);
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonAST(ostream& output) const
|
||
{
|
||
vector<pair<string, string>> eqtags;
|
||
output << R"("abstract_syntax_tree":[)" << endl;
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
if (eq != 0)
|
||
output << ", ";
|
||
|
||
output << R"({ "number":)" << eq;
|
||
if (equations_lineno[eq])
|
||
output << R"(, "line":)" << *equations_lineno[eq];
|
||
|
||
equation_tags.writeJsonAST(output, eq);
|
||
|
||
output << R"(, "AST": )";
|
||
equations[eq]->writeJsonAST(output);
|
||
output << "}";
|
||
}
|
||
output << "]";
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonVariableMapping(ostream& output) const
|
||
{
|
||
output << R"("variable_mapping":[)" << endl;
|
||
for (bool printed_something {false}; const auto& [var, eqs] : variableMapping)
|
||
{
|
||
if (exchange(printed_something, true))
|
||
output << ", ";
|
||
output << R"({"name": ")" << symbol_table.getName(var) << R"(", "equations":[)";
|
||
for (bool printed_something2 {false}; int it2 : eqs)
|
||
if (auto tmp = equation_tags.getTagValueByEqnAndKey(it2, "name"); tmp)
|
||
{
|
||
if (exchange(printed_something2, true))
|
||
output << ", ";
|
||
output << '"' << *tmp << '"';
|
||
}
|
||
output << "]}" << endl;
|
||
}
|
||
output << "]";
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonXrefsHelper(ostream& output,
|
||
const map<pair<int, int>, set<int>>& xrefmap) const
|
||
{
|
||
for (bool printed_something {false}; const auto& [symb_lag, eqs] : xrefmap)
|
||
{
|
||
if (exchange(printed_something, true))
|
||
output << ", ";
|
||
output << R"({"name": ")" << symbol_table.getName(symb_lag.first) << R"(")"
|
||
<< R"(, "shift": )" << symb_lag.second << R"(, "equations": [)";
|
||
for (bool printed_something2 {false}; int eq : eqs)
|
||
{
|
||
if (exchange(printed_something2, true))
|
||
output << ", ";
|
||
output << eq + 1;
|
||
}
|
||
output << "]}";
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonXrefs(ostream& output) const
|
||
{
|
||
output << R"("xrefs": {)"
|
||
<< R"("parameters": [)";
|
||
writeJsonXrefsHelper(output, xref_param);
|
||
output << "]"
|
||
<< R"(, "endogenous": [)";
|
||
writeJsonXrefsHelper(output, xref_endo);
|
||
output << "]"
|
||
<< R"(, "exogenous": [)";
|
||
writeJsonXrefsHelper(output, xref_exo);
|
||
output << "]"
|
||
<< R"(, "exogenous_deterministic": [)";
|
||
writeJsonXrefsHelper(output, xref_exo_det);
|
||
output << "]}" << endl;
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonOriginalModelOutput(ostream& output) const
|
||
{
|
||
writeJsonModelEquations(output, false);
|
||
output << ", ";
|
||
writeJsonAST(output);
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonDynamicModelInfo(ostream& output) const
|
||
{
|
||
output << R"("model_info": {)"
|
||
<< R"("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++)
|
||
{
|
||
if (endoID != 0)
|
||
output << ",";
|
||
output << "[";
|
||
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
|
||
{
|
||
if (lag != -max_endo_lag)
|
||
output << ",";
|
||
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 << "], "
|
||
<< R"("nstatic": )" << nstatic << ", " << endl
|
||
<< R"("nfwrd": )" << nfwrd << ", " << endl
|
||
<< R"("npred": )" << npred << ", " << endl
|
||
<< R"("nboth": )" << nboth << ", " << endl
|
||
<< R"("nsfwrd": )" << nfwrd + nboth << ", " << endl
|
||
<< R"("nspred": )" << npred + nboth << ", " << endl
|
||
<< R"("ndynamic": )" << npred + nboth + nfwrd << ", " << endl
|
||
<< R"("maximum_endo_lag": )" << max_endo_lag << ", " << endl
|
||
<< R"("maximum_endo_lead": )" << max_endo_lead << ", " << endl
|
||
<< R"("maximum_exo_lag": )" << max_exo_lag << ", " << endl
|
||
<< R"("maximum_exo_lead": )" << max_exo_lead << ", " << endl
|
||
<< R"("maximum_exo_det_lag": )" << max_exo_det_lag << ", " << endl
|
||
<< R"("maximum_exo_det_lead": )" << max_exo_det_lead << ", " << endl
|
||
<< R"("maximum_lag": )" << max_lag << ", " << endl
|
||
<< R"("maximum_lead": )" << max_lead << ", " << endl
|
||
<< R"("orig_maximum_endo_lag": )" << max_endo_lag_orig << "," << endl
|
||
<< R"("orig_maximum_endo_lead": )" << max_endo_lead_orig << "," << endl
|
||
<< R"("orig_maximum_exo_lag": )" << max_exo_lag_orig << "," << endl
|
||
<< R"("orig_maximum_exo_lead": )" << max_exo_lead_orig << "," << endl
|
||
<< R"("orig_maximum_exo_det_lag": )" << max_exo_det_lag_orig << "," << endl
|
||
<< R"("orig_maximum_exo_det_lead": )" << max_exo_det_lead_orig << "," << endl
|
||
<< R"("orig_maximum_lag": )" << max_lag_orig << "," << endl
|
||
<< R"("orig_maximum_lead": )" << max_lead_orig << "," << endl
|
||
<< R"("orig_maximum_lag_with_diffs_expanded": )" << max_lag_with_diffs_expanded_orig << ","
|
||
<< endl
|
||
<< R"("NNZDerivatives": [)";
|
||
for (int i = 1; i < static_cast<int>(NNZDerivatives.size()); i++)
|
||
{
|
||
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]);
|
||
if (i < static_cast<int>(NNZDerivatives.size()) - 1)
|
||
output << ", ";
|
||
}
|
||
output << "]}" << endl;
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonComputingPassOutput(ostream& output, bool writeDetails) const
|
||
{
|
||
auto [mlv_output, d_output] {writeJsonComputingPassOutputHelper<true>(writeDetails)};
|
||
|
||
if (writeDetails)
|
||
output << R"("dynamic_model": {)";
|
||
else
|
||
output << R"("dynamic_model_simple": {)";
|
||
output << mlv_output.str();
|
||
for (const auto& it : d_output)
|
||
output << ", " << it.str();
|
||
output << "}";
|
||
}
|
||
|
||
void
|
||
DynamicModel::writeJsonParamsDerivatives(ostream& output, bool writeDetails) const
|
||
{
|
||
if (!params_derivatives.size())
|
||
return;
|
||
|
||
auto [mlv_output, tt_output, rp_output, gp_output, rpp_output, gpp_output, hp_output,
|
||
g3p_output] {writeJsonParamsDerivativesHelper<true>(writeDetails)};
|
||
|
||
if (writeDetails)
|
||
output << R"("dynamic_model_params_derivative": {)";
|
||
else
|
||
output << R"("dynamic_model_params_derivatives_simple": {)";
|
||
output << mlv_output.str() << ", " << tt_output.str() << ", " << rp_output.str() << ", "
|
||
<< gp_output.str() << ", " << rpp_output.str() << ", " << gpp_output.str() << ", "
|
||
<< hp_output.str() << ", " << g3p_output.str() << "}";
|
||
}
|
||
|
||
void
|
||
DynamicModel::substituteVarExpectation(const map<string, expr_t>& subst_table)
|
||
{
|
||
for (auto& equation : equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->substituteVarExpectation(subst_table));
|
||
}
|
||
|
||
void
|
||
DynamicModel::checkNoRemainingPacExpectation() const
|
||
{
|
||
for (size_t eq = 0; eq < equations.size(); eq++)
|
||
if (equations[eq]->containsPacExpectation())
|
||
{
|
||
cerr << "ERROR: in equation "
|
||
<< equation_tags.getTagValueByEqnAndKey(eq, "name").value_or(to_string(eq + 1))
|
||
<< ", the pac_expectation operator references an unknown pac_model" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::simplifyEquations()
|
||
{
|
||
size_t last_subst_table_size = 0;
|
||
map<VariableNode*, NumConstNode*> subst_table;
|
||
// Equations with “mcp” tag are excluded, see dynare#1697
|
||
findConstantEquationsWithoutMcpTag(subst_table);
|
||
while (subst_table.size() != last_subst_table_size)
|
||
{
|
||
last_subst_table_size = subst_table.size();
|
||
for (auto& [id, definition] : local_variables_table)
|
||
definition = definition->replaceVarsInEquation(subst_table);
|
||
for (auto& equation : equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
|
||
for (auto& equation : static_only_equations)
|
||
equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
|
||
subst_table.clear();
|
||
findConstantEquationsWithoutMcpTag(subst_table);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::checkNoRemainingPacTargetNonstationary() const
|
||
{
|
||
for (size_t eq = 0; eq < equations.size(); eq++)
|
||
if (equations[eq]->containsPacTargetNonstationary())
|
||
{
|
||
cerr << "ERROR: in equation "
|
||
<< equation_tags.getTagValueByEqnAndKey(eq, "name").value_or(to_string(eq + 1))
|
||
<< ", the pac_target_nonstationary operator does not match a corresponding "
|
||
"'pac_target_info' block"
|
||
<< endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
void
|
||
DynamicModel::checkIsLinear() const
|
||
{
|
||
if (!nonzero_hessian_eqs.empty())
|
||
{
|
||
cerr << "ERROR: If the model is declared linear the second derivatives must be equal to zero."
|
||
<< endl
|
||
<< " The following equations have non-zero second derivatives:" << endl;
|
||
for (auto it : nonzero_hessian_eqs)
|
||
{
|
||
cerr << " * Eq # " << it + 1;
|
||
if (optional<string> eqname {equation_tags.getTagValueByEqnAndKey(it, "name")}; eqname)
|
||
cerr << " [" << *eqname << "]";
|
||
cerr << endl;
|
||
}
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|