Ramsey: write derivatives of static model w.r.t. Lagrange multipliers in a new file

The computing of the Ramsey steady state relies on the fact that Lagrange
multipliers appear linearly in the system to be solved. Instead of directly
solving for the Lagrange multipliers along with the other variables,
dyn_ramsey_static.m reduces the size of the problem by always computing the
value of the multipliers that minimizes the residuals, given the other
variables (using a minimum norm solution, easy to compute because of the
linearity property). That function thus needs the derivatives of the optimality
FOCs with respect to the multipliers. The problem is that, when multipliers
appear in an auxiliary variable related to a lead/lag, then those derivatives
need to be retrieved by a chain rule derivation, which cannot be easily done
with the regular static file.

This commit implements the creation of a new file,
ramsey_multipliers_static_g1.{m,mex}, that provides exactly the needed
derivatives w.r.t. Lagrange multipliers through chain rule derivation.

Ref. dynare#633, dynare#1119, dynare#1133
master
Sébastien Villemot 2023-03-24 18:58:12 +01:00
parent fe3f18947e
commit 0169240f76
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
8 changed files with 293 additions and 13 deletions

View File

@ -2548,19 +2548,19 @@ DynamicModel::replaceMyEquations(DynamicModel &dynamic_model) const
dynamic_model.equation_tags = equation_tags;
}
void
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
int i;
for (i = 0; i < static_cast<int>(equations.size()); i++)
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;
}
cout << "Ramsey Problem: added " << i << " Multipliers." << endl;
// Add Planner Objective to equations so that it appears in Lagrangian
assert(static_model.equations.size() == 1);
@ -2587,7 +2587,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
// Create (modified) Lagrangian (so that we can take the derivative once at time t)
expr_t lagrangian = Zero;
for (i = 0; i < static_cast<int>(equations.size()); i++)
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;
@ -2615,10 +2615,15 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
/* 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) */
(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;
@ -2633,7 +2638,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
neweqs_tags[neweqs.size()-1] = old_equation_tags.getTagsByEqn(*i);
}
else
neweqs_lineno.push_back(nullopt);
{
orig_endo_nbr++;
neweqs_lineno.push_back(nullopt);
}
}
}
@ -2641,6 +2649,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
clearEquations();
for (size_t i = 0; i < neweqs.size(); i++)
addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]);
return orig_endo_nbr;
}
bool

View File

@ -396,8 +396,13 @@ public:
void removeEquations(const vector<pair<string, string>> &listed_eqs_by_tag, bool exclude_eqs,
bool excluded_vars_change_type);
//! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
void computeRamseyPolicyFOCs(const StaticModel &static_model);
/* Replaces model equations with derivatives of Lagrangian w.r.t. endogenous.
The optimality FOCs (derivatives w.r.t. ordinary endogenous) will appear
first, followed by the constraints (derivatives w.r.t. multipliers).
Returns the number of optimality FOCs, which is by construction equal to
the number of endogenous before adding the Lagrange multipliers
(internally called ramsey_endo_nbr). */
int computeRamseyPolicyFOCs(const StaticModel &static_model);
//! Clears all equations
void clearEquations();

View File

@ -497,7 +497,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
orig_ramsey_dynamic_model = dynamic_model;
DynamicModel ramsey_FOC_equations_dynamic_model {symbol_table, num_constants, external_functions_table, trend_component_model_table, var_model_table};
ramsey_FOC_equations_dynamic_model = dynamic_model;
ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective);
mod_file_struct.ramsey_orig_endo_nbr = ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective);
ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
}
@ -673,6 +673,9 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
paramsDerivsOrder = params_derivs_order;
static_model.computingPass(derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll);
if (mod_file_struct.ramsey_model_present)
static_model.computeRamseyMultipliersDerivatives(mod_file_struct.ramsey_orig_endo_nbr,
!use_dll, no_tmp_terms);
}
// Set things to compute for dynamic model
if (mod_file_struct.perfect_foresight_solver_present
@ -937,7 +940,11 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
{
dynamic_model.writeDriverOutput(mOutputFile, compute_xrefs);
if (!no_static)
static_model.writeDriverOutput(mOutputFile);
{
static_model.writeDriverOutput(mOutputFile);
if (mod_file_struct.ramsey_model_present)
static_model.writeDriverRamseyMultipliersDerivativesSparseIndices(mOutputFile);
}
}
if (onlymodel || gui)
@ -1044,6 +1051,13 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
{
static_model.writeStaticFile(basename, use_dll, mexext, matlabroot, false);
static_model.writeParamsDerivativesFile<false>(basename);
if (mod_file_struct.ramsey_model_present)
{
if (use_dll)
static_model.writeRamseyMultipliersDerivativesCFile(basename, mexext, matlabroot, mod_file_struct.ramsey_orig_endo_nbr);
else
static_model.writeRamseyMultipliersDerivativesMFile(basename, mod_file_struct.ramsey_orig_endo_nbr);
}
}
dynamic_model.writeDynamicFile(basename, use_dll, mexext, matlabroot, false);

View File

@ -1,5 +1,5 @@
/*
* Copyright © 2006-2022 Dynare Team
* Copyright © 2006-2023 Dynare Team
*
* This file is part of Dynare.
*
@ -135,6 +135,11 @@ public:
int orig_eq_nbr{0};
//! Stores the number of equations added to the Ramsey model
int ramsey_eq_nbr{0};
/* The number of endogenous variables in the model present just before adding
the Lagrange multipliers and computing the Ramsey FOC; it is by
construction equal to the number of equations that will be added by the
process of computing the FOCs */
int ramsey_orig_endo_nbr {0};
//! Whether there was a steady_state_model block
bool steady_state_model_present{false};
//! Whether there is a write_latex_steady_state_model statement present

View File

@ -862,3 +862,172 @@ StaticModel::writeJsonParamsDerivatives(ostream &output, bool writeDetails) cons
<< ", " << hp_output.str()
<< "}";
}
void
StaticModel::computeRamseyMultipliersDerivatives(int ramsey_orig_endo_nbr, bool is_matlab,
bool no_tmp_terms)
{
// Compute derivation IDs of Lagrange multipliers
set<int> mult_symb_ids { symbol_table.getLagrangeMultipliers() };
vector<int> mult_deriv_ids;
for (int symb_id : mult_symb_ids)
mult_deriv_ids.push_back(getDerivID(symb_id, 0));
// Compute the list of aux vars for which to apply the chain rule derivation
map<int, BinaryOpNode *> recursive_variables;
for (auto aux_eq : aux_equations)
{
auto varexpr { dynamic_cast<VariableNode *>(aux_eq->arg1) };
assert(varexpr && symbol_table.isAuxiliaryVariable(varexpr->symb_id));
/* Determine whether the auxiliary variable has been added after the last
Lagrange multiplier. We use the guarantee given by SymbolTable that
symbol IDs are increasing. */
if (varexpr->symb_id > *mult_symb_ids.crbegin())
recursive_variables.emplace(getDerivID(varexpr->symb_id, 0), aux_eq);
}
// Compute the chain rule derivatives w.r.t. multipliers
map<expr_t, set<int>> non_null_chain_rule_derivatives;
map<pair<expr_t, int>, expr_t> cache;
for (int eq {0}; eq < ramsey_orig_endo_nbr; eq++)
for (int mult {0}; mult < static_cast<int>(mult_deriv_ids.size()); mult++)
if (expr_t d { equations[eq]->getChainRuleDerivative(mult_deriv_ids[mult], recursive_variables,
non_null_chain_rule_derivatives, cache) };
d != Zero)
ramsey_multipliers_derivatives.try_emplace({ eq, mult }, d);
// Compute the temporary terms
map<pair<int, int>, temporary_terms_t> temp_terms_map;
map<expr_t, pair<int, pair<int, int>>> reference_count;
for (const auto &[row_col, d] : ramsey_multipliers_derivatives)
d->computeTemporaryTerms({ 1, 0 }, temp_terms_map, reference_count, is_matlab);
/* If the user has specified the notmpterms option, clear all temporary
terms, except those that correspond to external functions (since they are
not optional) */
if (no_tmp_terms)
for (auto &it : temp_terms_map)
erase_if(it.second,
[](expr_t e) { return !dynamic_cast<AbstractExternalFunctionNode *>(e); });
ramsey_multipliers_derivatives_temporary_terms = move(temp_terms_map[{ 1, 0 }]);
for (int idx {0};
auto it : ramsey_multipliers_derivatives_temporary_terms)
ramsey_multipliers_derivatives_temporary_terms_idxs[it] = idx++;
// Compute the CSC format
ramsey_multipliers_derivatives_sparse_colptr = computeCSCColPtr(ramsey_multipliers_derivatives,
mult_deriv_ids.size());
}
void
StaticModel::writeDriverRamseyMultipliersDerivativesSparseIndices(ostream &output) const
{
output << "M_.ramsey_multipliers_static_g1_sparse_rowval = int32([";
for (auto &[row_col, d] : ramsey_multipliers_derivatives)
output << row_col.first+1 << ' ';
output << "]);" << endl
<< "M_.ramsey_multipliers_static_g1_sparse_colval = int32([";
for (auto &[row_col, d] : ramsey_multipliers_derivatives)
output << row_col.second+1 << ' ';
output << "]);" << endl
<< "M_.ramsey_multipliers_static_g1_sparse_colptr = int32([";
for (int it : ramsey_multipliers_derivatives_sparse_colptr)
output << it+1 << ' ';
output << "]);" << endl;
}
void
StaticModel::writeRamseyMultipliersDerivativesMFile(const string &basename, int ramsey_orig_endo_nbr) const
{
constexpr auto output_type { ExprNodeOutputType::matlabStaticModel };
filesystem::path filename {packageDir(basename) / "ramsey_multipliers_static_g1.m"};
ofstream output_file{filename, ios::out | ios::binary};
if (!output_file.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output_file << "function g1m = ramsey_multipliers_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr)" << endl
<< "g1m_v=NaN(" << ramsey_multipliers_derivatives.size() << ",1);" << endl;
writeRamseyMultipliersDerivativesHelper<output_type>(output_file);
// On MATLAB < R2020a, sparse() does not accept int32 indices
output_file << "if ~isoctave && matlab_ver_less_than('9.8')" << endl
<< " sparse_rowval = double(sparse_rowval);" << endl
<< " sparse_colval = double(sparse_colval);" << endl
<< "end" << endl
<< "g1m = sparse(sparse_rowval, sparse_colval, g1m_v, " << ramsey_orig_endo_nbr << ", " << symbol_table.getLagrangeMultipliers().size() << ");" << endl
<< "end" << endl;
output_file.close();
}
void
StaticModel::writeRamseyMultipliersDerivativesCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, int ramsey_orig_endo_nbr) const
{
constexpr auto output_type { ExprNodeOutputType::CStaticModel };
const filesystem::path model_src_dir {filesystem::path{basename} / "model" / "src"};
const int xlen { symbol_table.exo_nbr()+symbol_table.exo_det_nbr() };
const int nzval { static_cast<int>(ramsey_multipliers_derivatives.size()) };
const int ncols { static_cast<int>(symbol_table.getLagrangeMultipliers().size()) };
const filesystem::path p {model_src_dir / "ramsey_multipliers_static_g1.c"};
ofstream output{p, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "#include <math.h>" << endl << endl
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
<< endl;
writePowerDeriv(output);
output << endl
<< "void ramsey_multipliers_static_g1(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict g1m_v)" << endl
<< "{" << endl;
writeRamseyMultipliersDerivativesHelper<output_type>(output);
output << "}" << endl
<< endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " if (nrhs != 6)" << endl
<< R"( mexErrMsgTxt("Accepts exactly 6 input arguments");)" << endl
<< " if (nlhs != 1)" << endl
<< R"( mexErrMsgTxt("Accepts exactly 1 output argument");)" << endl
<< " if (!(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0]) && !mxIsSparse(prhs[0]) && mxGetNumberOfElements(prhs[0]) == " << symbol_table.endo_nbr() << "))" << endl
<< R"( mexErrMsgTxt("y must be a real dense numeric array with )" << symbol_table.endo_nbr() << R"( elements");)" << endl
<< " const double *restrict y = mxGetPr(prhs[0]);" << endl
<< " if (!(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) && !mxIsSparse(prhs[1]) && mxGetNumberOfElements(prhs[1]) == " << xlen << "))" << endl
<< R"( mexErrMsgTxt("x must be a real dense numeric array with )" << xlen << R"( elements");)" << endl
<< " const double *restrict x = mxGetPr(prhs[1]);" << endl
<< " if (!(mxIsDouble(prhs[2]) && !mxIsComplex(prhs[2]) && !mxIsSparse(prhs[2]) && mxGetNumberOfElements(prhs[2]) == " << symbol_table.param_nbr() << "))" << endl
<< R"( mexErrMsgTxt("params must be a real dense numeric array with )" << symbol_table.param_nbr() << R"( elements");)" << endl
<< " const double *restrict params = mxGetPr(prhs[2]);" << endl
<< " if (!(mxIsInt32(prhs[3]) && mxGetNumberOfElements(prhs[3]) == " << nzval << "))" << endl
<< R"( mexErrMsgTxt("sparse_rowval must be an int32 array with )" << nzval << R"( elements");)" << endl
<< " if (!(mxIsInt32(prhs[5]) && mxGetNumberOfElements(prhs[5]) == " << ncols+1 << "))" << endl
<< R"( mexErrMsgTxt("sparse_colptr must be an int32 array with )" << ncols+1 << R"( elements");)" << endl
<< "#if MX_HAS_INTERLEAVED_COMPLEX" << endl
<< " const int32_T *restrict sparse_rowval = mxGetInt32s(prhs[3]);" << endl
<< " const int32_T *restrict sparse_colptr = mxGetInt32s(prhs[5]);" << endl
<< "#else" << endl
<< " const int32_T *restrict sparse_rowval = (int32_T *) mxGetData(prhs[3]);" << endl
<< " const int32_T *restrict sparse_colptr = (int32_T *) mxGetData(prhs[5]);" << endl
<< "#endif" << endl
<< " plhs[0] = mxCreateSparse(" << ramsey_orig_endo_nbr << ", " << ncols << ", " << nzval << ", mxREAL);" << endl
<< " mwIndex *restrict ir = mxGetIr(plhs[0]), *restrict jc = mxGetJc(plhs[0]);" << endl
<< " for (mwSize i = 0; i < " << nzval << "; i++)" << endl
<< " *ir++ = *sparse_rowval++ - 1;" << endl
<< " for (mwSize i = 0; i < " << ncols+1 << "; i++)" << endl
<< " *jc++ = *sparse_colptr++ - 1;" << endl
<< " mxArray *T_mx = mxCreateDoubleMatrix(" << ramsey_multipliers_derivatives_temporary_terms.size() << ", 1, mxREAL);" << endl
<< " ramsey_multipliers_static_g1(y, x, params, mxGetPr(T_mx), mxGetPr(plhs[0]));" << endl
<< " mxDestroyArray(T_mx);" << endl
<< "}" << endl;
output.close();
compileMEX(packageDir(basename), "ramsey_multipliers_static_g1", mexext, { p }, matlabroot);
}

View File

@ -34,6 +34,25 @@ class DynamicModel;
class StaticModel : public ModelTree
{
private:
/* First-order derivatives of equations w.r.t. Lagrange multipliers, using
chain rule derivation for auxiliary variables added after the multipliers
(so that derivatives of optimality FOCs w.r.t. multipliers with lead or
lag 2 are self-contained, which is required by dyn_ramsey_static.m).
Only used if 'ramsey_model' or 'ramsey_policy' is present.
The first index of the key is the equation number (NB: auxiliary equations
added after the multipliers do not appear).
The second index is the index of the Lagrange multiplier (ordered by
increasing symbol ID) */
SparseColumnMajorOrderMatrix ramsey_multipliers_derivatives;
/* Column indices for the derivatives w.r.t. Lagrange multipliers in
Compressed Sparse Column (CSC) storage (corresponds to the jc vector in
MATLAB terminology) */
vector<int> ramsey_multipliers_derivatives_sparse_colptr;
// Temporary terms for ramsey_multipliers_derivatives
temporary_terms_t ramsey_multipliers_derivatives_temporary_terms;
// Stores, for each temporary term, its index in the MATLAB/Octave vector
temporary_terms_idxs_t ramsey_multipliers_derivatives_temporary_terms_idxs;
// Writes static model file (MATLAB/Octave version, legacy representation)
void writeStaticMFile(const string &basename) const;
@ -94,6 +113,10 @@ private:
// Write the block structure of the model in the driver file
void writeBlockDriverOutput(ostream &output) const;
// Helper for writing ramsey_multipliers_derivatives
template<ExprNodeOutputType output_type>
void writeRamseyMultipliersDerivativesHelper(ostream &output) const;
protected:
string
modelClassName() const override
@ -158,6 +181,18 @@ public:
int getDerivID(int symb_id, int lag) const noexcept(false) override;
void addAllParamDerivId(set<int> &deriv_id_set) override;
// Fills the ramsey_multipliers_derivatives structure (see the comment there)
void computeRamseyMultipliersDerivatives(int ramsey_orig_endo_nbr, bool is_matlab, bool no_tmp_terms);
// Writes the sparse indices of ramsey_multipliers_derivatives to the driver file
void writeDriverRamseyMultipliersDerivativesSparseIndices(ostream &output) const;
// Writes ramsey_multipliers_derivatives (MATLAB/Octave version)
void writeRamseyMultipliersDerivativesMFile(const string &basename, int ramsey_orig_endo_nbr) const;
// Writes ramsey_multipliers_derivatives (C version)
void writeRamseyMultipliersDerivativesCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, int ramsey_orig_endo_nbr) const;
};
template<bool julia>
@ -264,4 +299,30 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
}
}
template<ExprNodeOutputType output_type>
void
StaticModel::writeRamseyMultipliersDerivativesHelper(ostream &output) const
{
// Write temporary terms (which includes external function stuff)
deriv_node_temp_terms_t tef_terms;
temporary_terms_t unused_tt_copy;
writeTemporaryTerms<output_type>(ramsey_multipliers_derivatives_temporary_terms,
unused_tt_copy,
ramsey_multipliers_derivatives_temporary_terms_idxs,
output, tef_terms);
// Write chain rule derivatives
for (int k {0};
auto &[row_col, d] : ramsey_multipliers_derivatives)
{
output << "g1m_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
d->writeOutput(output, output_type, ramsey_multipliers_derivatives_temporary_terms,
ramsey_multipliers_derivatives_temporary_terms_idxs, tef_terms);
output << ";" << endl;
k++;
}
}
#endif

View File

@ -1026,3 +1026,13 @@ SymbolTable::getVariablesWithLogTransform() const
{
return with_log_transform;
}
set<int>
SymbolTable::getLagrangeMultipliers() const
{
set<int> r;
for (const auto &aux_var : aux_vars)
if (aux_var.type == AuxVarType::multiplier)
r.insert(aux_var.symb_id);
return r;
}

View File

@ -91,7 +91,11 @@ struct AuxVarInfo
//! Stores the symbol table
/*!
A symbol is given by its name, and is internally represented by a unique integer.
A symbol is given by its name, and is internally represented by a unique
integer, called a symbol ID.
There is a guarantee that symbol IDs are increasing, i.e. if symbol A is
added after symbol B, then the ID of A is greater than the ID of B.
When method freeze() is called, computes a distinct sequence of IDs for some types
(endogenous, exogenous, parameters), which are used by the Matlab/Octave functions.
@ -407,6 +411,8 @@ public:
const AuxVarInfo &getAuxVarInfo(int symb_id) const;
// Returns the set of all endogenous declared with “var(log)”
const set<int> &getVariablesWithLogTransform() const;
// Returns all Lagrange multipliers
set<int> getLagrangeMultipliers() const;
};
inline void