stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Block decomposition: ensure that the cutoff option only impacts normalization 

Previously, the cutoff option would also impact the block decomposition itself,
since it would had an influence on the incidence matrix used for computing the
blocks and their derivatives.

The problem is that, in the general case, it’s quite possible that an element
of the numerical Jacobian be zero at the evaluation point, while being quite
different from zero along the simulation path. A typical example is an
expression of the type x*y, where y is an endogenous and x is an exogenous not
present in the initval block (and hence initialized to zero).
issue#70
Sébastien Villemot 2020-04-28 18:00:26 +02:00
parent 470cb5fcb0
commit 39407083be
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
4 changed files with 72 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
src/DynamicModel.cc
View File

@ -3228,7 +3228,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
tmp_s.str(""); tmp_s.str("");
count_lead_lag_incidence = 0; count_lead_lag_incidence = 0;
dynamic_jacob_map_t reordered_dynamic_jacobian; map<tuple<int, int, int>, expr_t> reordered_dynamic_jacobian;
for (const auto &[idx, d] : blocks_derivatives[block]) for (const auto &[idx, d] : blocks_derivatives[block])
reordered_dynamic_jacobian[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d; reordered_dynamic_jacobian[{ get<2>(idx), get<1>(idx), get<0>(idx) }] = d;
output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];" << endl; output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];" << endl;
@ -4812,7 +4812,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian); computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
computePrologueAndEpilogue(static_jacobian); computePrologueAndEpilogue();
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous(); auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
@ -4820,7 +4820,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
cout << "Finding the optimal block decomposition of the model ..." << endl; cout << "Finding the optimal block decomposition of the model ..." << endl;
auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false); auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock();
reduceBlocksAndTypeDetermination(linear_decomposition); reduceBlocksAndTypeDetermination(linear_decomposition);
@ -4971,26 +4971,31 @@ DynamicModel::determineBlockDerivativesType(int blk)
int nb_recursive = blocks[blk].getRecursiveSize(); int nb_recursive = blocks[blk].getRecursiveSize();
for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++) for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++)
for (int eq = 0; eq < size; eq++) for (int eq = 0; eq < size; eq++)
for (int var = 0; var < size; var++) {
if (int eq_orig = getBlockEquationID(blk, eq), var_orig = getBlockVariableID(blk, var); set<pair<int, int>> endos_and_lags;
dynamic_jacobian.find({ lag, eq_orig, var_orig }) != dynamic_jacobian.end()) int eq_orig = getBlockEquationID(blk, eq);
{ equations[eq_orig]->collectEndogenous(endos_and_lags);
if (getBlockEquationType(blk, eq) == EquationType::evaluate_s for (int var = 0; var < size; var++)
&& eq < nb_recursive) if (int var_orig = getBlockVariableID(blk, var);
/* It’s a normalized recursive equation, we have to recompute endos_and_lags.find({ var_orig, lag }) != endos_and_lags.end())
the derivative using the chain rule */ {
derivType[{ lag, eq, var }] = BlockDerivativeType::normalizedChainRule; if (getBlockEquationType(blk, eq) == EquationType::evaluate_s
else if (derivType.find({ lag, eq, var }) == derivType.end()) && eq < nb_recursive)
derivType[{ lag, eq, var }] = BlockDerivativeType::standard; /* 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.find({ lag, eq, var }) == derivType.end())
derivType[{ lag, eq, var }] = BlockDerivativeType::standard;
if (var < nb_recursive) if (var < nb_recursive)
for (int feedback_var = nb_recursive; feedback_var < size; feedback_var++) for (int feedback_var = nb_recursive; feedback_var < size; feedback_var++)
if (derivType.find({ lag, var, feedback_var }) != derivType.end()) if (derivType.find({ lag, var, feedback_var }) != derivType.end())
/* A new derivative needs to be computed using the chain rule /* A new derivative needs to be computed using the chain rule
(a feedback variable appears in the recursive equation (a feedback variable appears in the recursive equation
defining the current variable) */ defining the current variable) */
derivType[{ lag, eq, feedback_var }] = BlockDerivativeType::chainRule; derivType[{ lag, eq, feedback_var }] = BlockDerivativeType::chainRule;
} }
}
return derivType; return derivType;
} }

129
src/ModelTree.cc
View File

@ -101,8 +101,6 @@ ModelTree::copyHelper(const ModelTree &m)
trend_symbols_map[it.first] = f(it.second); trend_symbols_map[it.first] = f(it.second);
for (const auto &it : m.nonstationary_symbols_map) for (const auto &it : m.nonstationary_symbols_map)
nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)}; nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)};
for (const auto &it : m.dynamic_jacobian)
dynamic_jacobian[it.first] = f(it.second);
for (const auto &it : m.equation_type_and_normalized_equation) for (const auto &it : m.equation_type_and_normalized_equation)
equation_type_and_normalized_equation.emplace_back(it.first, f(it.second)); equation_type_and_normalized_equation.emplace_back(it.first, f(it.second));
@ -200,7 +198,6 @@ ModelTree::operator=(const ModelTree &m)
trend_symbols_map.clear(); trend_symbols_map.clear();
nonstationary_symbols_map.clear(); nonstationary_symbols_map.clear();
dynamic_jacobian.clear();
eq_idx_block2orig = m.eq_idx_block2orig; eq_idx_block2orig = m.eq_idx_block2orig;
endo_idx_block2orig = m.endo_idx_block2orig; endo_idx_block2orig = m.endo_idx_block2orig;
eq_idx_orig2block = m.eq_idx_orig2block; eq_idx_orig2block = m.eq_idx_orig2block;
@ -285,7 +282,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
} }
void void
ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian) ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, double cutoff, const jacob_map_t &static_jacobian)
{ {
cout << "Normalizing the model..." << endl; cout << "Normalizing the model..." << endl;
@ -338,33 +335,6 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
TODO: Explain why symbolic_jacobian is not contemporaneous. */ TODO: Explain why symbolic_jacobian is not contemporaneous. */
auto symbolic_jacobian = computeSymbolicJacobian(); auto symbolic_jacobian = computeSymbolicJacobian();
found_normalization = computeNormalization(symbolic_jacobian, true); found_normalization = computeNormalization(symbolic_jacobian, true);
if (found_normalization)
{
/* Update the Jacobian matrices by ensuring that they have a
numerical value associated to each symbolic occurrence.
TODO: Explain why this is needed. Maybe in order to have the
right incidence matrix in computePrologueAndEpilogue? */
for (const auto &[eq_and_endo, ignore] : symbolic_jacobian)
{
if (static_jacobian.find(eq_and_endo) == static_jacobian.end())
static_jacobian[eq_and_endo] = 0;
if (dynamic_jacobian.find({ 0, eq_and_endo.first, eq_and_endo.second }) == dynamic_jacobian.end())
dynamic_jacobian[{ 0, eq_and_endo.first, eq_and_endo.second }] = Zero;
if (contemporaneous_jacobian.find(eq_and_endo) == contemporaneous_jacobian.end())
contemporaneous_jacobian[eq_and_endo] = 0;
try
{
if (derivatives[1].find({ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }) == derivatives[1].end())
derivatives[1][{ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }] = Zero;
}
catch (DataTree::UnknownDerivIDException &e)
{
cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second))
<< " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
exit(EXIT_FAILURE);
}
}
}
} }
if (!found_normalization) if (!found_normalization)
@ -375,7 +345,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
} }
pair<ModelTree::jacob_map_t, ModelTree::jacob_map_t> pair<ModelTree::jacob_map_t, ModelTree::jacob_map_t>
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) const
{ {
jacob_map_t contemporaneous_jacobian, static_jacobian; jacob_map_t contemporaneous_jacobian, static_jacobian;
int nb_elements_contemporaneous_jacobian = 0; int nb_elements_contemporaneous_jacobian = 0;
@ -423,17 +393,10 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
static_jacobian[{ eq, var }] += val; static_jacobian[{ eq, var }] += val;
else else
static_jacobian[{ eq, var }] = val; static_jacobian[{ eq, var }] = val;
dynamic_jacobian[{ lag, eq, var }] = d1;
} }
} }
} }
// Get rid of the elements of the Jacobian matrix below the cutoff
// TODO: Why?
for (const auto &it : jacobian_elements_to_delete)
derivatives[1].erase(it);
if (jacobian_elements_to_delete.size() > 0) if (jacobian_elements_to_delete.size() > 0)
{ {
cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded." << endl cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded." << endl
@ -514,7 +477,7 @@ ModelTree::computeNaturalNormalization()
} }
void void
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian) ModelTree::computePrologueAndEpilogue()
{ {
const int n = equations.size(); const int n = equations.size();
@ -533,17 +496,13 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian)
“eq_idx_block2orig” (resp. “endo_idx_block2orig”). Stored in row-major “eq_idx_block2orig” (resp. “endo_idx_block2orig”). Stored in row-major
order. */ order. */
vector<bool> IM(n*n, false); vector<bool> IM(n*n, false);
if (cutoff == 0) for (int i = 0; i < n; i++)
for (int i = 0; i < n; i++) {
{ set<pair<int, int>> endos_and_lags;
set<pair<int, int>> endos_and_lags; equations[i]->collectEndogenous(endos_and_lags);
equations[i]->collectEndogenous(endos_and_lags); for (auto [endo, lag] : endos_and_lags)
for (const auto &[endo, lag] : endos_and_lags) IM[i * n + endo2eq[endo]] = true;
IM[i * n + endo2eq[endo]] = true; }
}
else
for (const auto &[eq_and_endo, val] : static_jacobian)
IM[eq_and_endo.first * n + endo2eq[eq_and_endo.second]] = true;
bool something_has_been_done; bool something_has_been_done;
// Find the prologue equations and place first the AR(1) shock equations first // Find the prologue equations and place first the AR(1) shock equations first
@ -677,22 +636,24 @@ ModelTree::getVariableLeadLagByBlock(const vector<int> &endo2simblock) const
}; };
lag_lead_vector_t variable_lag_lead(nb_endo, { 0, 0 }), equation_lag_lead(nb_endo, { 0, 0 }); lag_lead_vector_t variable_lag_lead(nb_endo, { 0, 0 }), equation_lag_lead(nb_endo, { 0, 0 });
for (const auto &[key, value] : dynamic_jacobian) for (int eq = 0; eq < nb_endo; eq++)
{ {
auto [lag, eq, endo] = key; set<pair<int, int>> endos_and_lags;
if (belong_to_same_block(endo, eq)) equations[eq]->collectEndogenous(endos_and_lags);
{ for (auto [endo, lag] : endos_and_lags)
variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag); if (belong_to_same_block(endo, eq))
variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag); {
equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag); variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag);
equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag); variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag);
} equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag);
equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag);
}
} }
return { equation_lag_lead, variable_lag_lead }; return { equation_lag_lead, variable_lag_lead };
} }
lag_lead_vector_t lag_lead_vector_t
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_) ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock()
{ {
int nb_var = symbol_table.endo_nbr(); int nb_var = symbol_table.endo_nbr();
int nb_simvars = nb_var - prologue - epilogue; int nb_simvars = nb_var - prologue - epilogue;
@ -700,11 +661,10 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
/* Construct the graph representing the dependencies between all /* Construct the graph representing the dependencies between all
variables that do not belong to the prologue or the epilogue. variables that do not belong to the prologue or the epilogue.
For detecting dependencies between variables, use the static jacobian, For detecting dependencies between variables, use the symbolic adjacency
except when the cutoff is zero, in which case use the symbolic adjacency
matrix */ matrix */
VariableDependencyGraph G(nb_simvars); VariableDependencyGraph G(nb_simvars);
for (const auto &[key, value] : cutoff == 0 ? computeSymbolicJacobian() : static_jacobian) for (const auto &[key, value] : computeSymbolicJacobian())
{ {
auto [eq, endo] = key; auto [eq, endo] = key;
if (eq_idx_orig2block[eq] >= prologue if (eq_idx_orig2block[eq] >= prologue
@ -869,29 +829,14 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition)
set<pair<int, int>> endos_and_lags; set<pair<int, int>> endos_and_lags;
equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags); equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags);
for (const auto &[endo, lag] : endos_and_lags) for (const auto &[endo, lag] : endos_and_lags)
{ if (linear_decomposition ||
if (linear_decomposition) find(endo_idx_block2orig.begin()+first_eq,
{ endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo)
if (dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo }) != endo_idx_block2orig.begin()+first_eq+blocks[i].size)
!= dynamic_jacobian.end()) {
{ max_lead = max(lag, max_lead);
max_lead = max(lag, max_lead); max_lag = max(-lag, max_lag);
max_lag = max(-lag, max_lag); }
}
}
else
{
if (find(endo_idx_block2orig.begin()+first_eq,
endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo)
!= endo_idx_block2orig.begin()+first_eq+blocks[i].size
&& dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo })
!= dynamic_jacobian.end())
{
max_lead = max(lag, max_lead);
max_lag = max(-lag, max_lag);
}
}
}
} }
// Determine the block type // Determine the block type
@ -940,12 +885,10 @@ ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition)
for (int j = blocks[i-1].first_equation; for (int j = blocks[i-1].first_equation;
j < blocks[i-1].first_equation+blocks[i-1].size; j++) j < blocks[i-1].first_equation+blocks[i-1].size; j++)
{ {
if (dynamic_jacobian.find({ -1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] }) set<pair<int, int>> endos_and_lags;
!= dynamic_jacobian.end()) equations[eq_idx_block2orig[first_eq]]->collectEndogenous(endos_and_lags);
is_lag = true; is_lag = endos_and_lags.find({ endo_idx_block2orig[j], -1 }) != endos_and_lags.end();
if (dynamic_jacobian.find({ +1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] }) is_lead = endos_and_lags.find({ endo_idx_block2orig[j], 1 }) != endos_and_lags.end();
!= dynamic_jacobian.end())
is_lead = true;
} }
BlockSimulationType prev_Type = blocks[i-1].simulation_type; BlockSimulationType prev_Type = blocks[i-1].simulation_type;

21
src/ModelTree.hh
View File

@ -155,13 +155,6 @@ protected:
//! Nonstationary variables and their deflators //! Nonstationary variables and their deflators
nonstationary_symbols_map_t nonstationary_symbols_map; nonstationary_symbols_map_t nonstationary_symbols_map;
//! Sparse matrix of double to store the values of the Jacobian
/*! First index is lag, second index is equation number, third index is endogenous type specific ID */
using dynamic_jacob_map_t = map<tuple<int, int, int>, expr_t>;
//! The jacobian without the elements below the cutoff
dynamic_jacob_map_t dynamic_jacobian;
/* Maps indices of equations in the block-decomposition order into original /* Maps indices of equations in the block-decomposition order into original
equation IDs */ equation IDs */
vector<int> eq_idx_block2orig; vector<int> eq_idx_block2orig;
@ -255,7 +248,7 @@ protected:
//! Writes LaTeX model file //! Writes LaTeX model file
void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const; void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const;
//! Sparse matrix of double to store the values of the Jacobian //! Sparse matrix of double to store the values of the static Jacobian
/*! First index is equation number, second index is endogenous type specific ID */ /*! First index is equation number, second index is endogenous type specific ID */
using jacob_map_t = map<pair<int, int>, double>; using jacob_map_t = map<pair<int, int>, double>;
@ -287,24 +280,24 @@ protected:
If no matching is found with a zero cutoff, an error message is printed. If no matching is found with a zero cutoff, an error message is printed.
The resulting normalization is stored in endo2eq. The resulting normalization is stored in endo2eq.
*/ */
void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian); void computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, double cutoff, const jacob_map_t &static_jacobian);
//! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS //! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS
bool computeNaturalNormalization(); bool computeNaturalNormalization();
//! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff //! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff
/*! Returns a pair (contemporaneous_jacobian, static_jacobian). Also fills /*! Returns a pair (contemporaneous_jacobian, static_jacobian).
dynamic_jacobian. Elements below the cutoff are discarded. External functions are evaluated to 1. */ Elements below the cutoff are discarded. External functions are evaluated to 1. */
pair<jacob_map_t, jacob_map_t> evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose); pair<jacob_map_t, jacob_map_t> evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose) const;
//! Select and reorder the non linear equations of the model //! Select and reorder the non linear equations of the model
void select_non_linear_equations_and_variables(); void select_non_linear_equations_and_variables();
//! Search the equations and variables belonging to the prologue and the epilogue of the model //! Search the equations and variables belonging to the prologue and the epilogue of the model
void computePrologueAndEpilogue(const jacob_map_t &static_jacobian); void computePrologueAndEpilogue();
//! Determine the type of each equation of model and try to normalize the unnormalized equation //! Determine the type of each equation of model and try to normalize the unnormalized equation
void equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs); void equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs);
/* Compute the block decomposition and for a non-recusive block find the minimum feedback set /* Compute the block decomposition and for a non-recusive block find the minimum feedback set
Initializes the “blocks” structure, and fills the following fields: size, Initializes the “blocks” structure, and fills the following fields: size,
mfs_size, n_static, n_forward, n_backward, n_mixed. */ mfs_size, n_static, n_forward, n_backward, n_mixed. */
lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_); lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock();
/* Reduce the number of block by merging the same type of equations in the /* Reduce the number of block by merging the same type of equations in the
prologue and the epilogue, and determine the type of each block. prologue and the epilogue, and determine the type of each block.

4
src/StaticModel.cc
View File

@ -1125,7 +1125,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian); computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
computePrologueAndEpilogue(static_jacobian); computePrologueAndEpilogue();
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous(); auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
@ -1133,7 +1133,7 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
cout << "Finding the optimal block decomposition of the model ..." << endl; cout << "Finding the optimal block decomposition of the model ..." << endl;
auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false); auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock();
reduceBlocksAndTypeDetermination(false); reduceBlocksAndTypeDetermination(false);