From 0169240f76532df6d507b01647e8576d1094e95d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 24 Mar 2023 18:58:12 +0100
Subject: [PATCH] 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
---
 src/DynamicModel.cc |  24 +++++--
 src/DynamicModel.hh |   9 ++-
 src/ModFile.cc      |  18 ++++-
 src/Statement.hh    |   7 +-
 src/StaticModel.cc  | 169 ++++++++++++++++++++++++++++++++++++++++++++
 src/StaticModel.hh  |  61 ++++++++++++++++
 src/SymbolTable.cc  |  10 +++
 src/SymbolTable.hh  |   8 ++-
 8 files changed, 293 insertions(+), 13 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 32ff9f2c..4178c447 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -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
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index eb6c6f14..4bce7c3d 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -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();
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 1e3c38b9..63ef9a43 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -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);
diff --git a/src/Statement.hh b/src/Statement.hh
index 7fc35150..7d555926 100644
--- a/src/Statement.hh
+++ b/src/Statement.hh
@@ -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
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 038ea375..affe3c51 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -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);
+}
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 77101304..96f10d3a 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -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
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index 37e3841d..d7cee2bf 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -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;
+}
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index e8e787be..316882e9 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -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
-- 
GitLab