diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 32ff9f2cc88fa1563a337ee33715a63eecc4f68e..4178c44738eefd91a9a85f44aac45e58e4bc0b3a 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 eb6c6f14a0a82840301e3b811d32ad24655cbd0c..4bce7c3d63d3a06f940274ffe48c960829c9a62c 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 1e3c38b9699f863c171cfdb3789eece539924201..63ef9a43dad20543d9a5d0dbfb78c5909290c6aa 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 7fc35150f96f84971abbbd38d7d5f5c61770ab41..7d5559268691b33597d506bd09267c7ae435e30a 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 038ea375e40a22ed0cb1eea55c1f6c9858151709..affe3c51bccb26bd0cad3124fada20cdf39a91fe 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 771013046e7e8912981681f469c1350cf0ab891d..96f10d3a6562944f30b2232b2ed3a2526c66ea92 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 37e3841db6f42beb255a9e5bc33f8ada464f7e31..d7cee2bf79c1ad4410d9dbd04af02c4ad84779d2 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 e8e787be7fa29c449ee5ee7e100a65176688849c..316882e9d84ef06dcc886879ef943bafe6d90588 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