diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 6116b85a..5168145b 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1691,7 +1691,7 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
 }
 
 void
-DynamicModel::writeDynamicCFile(const string &basename, int order) const
+DynamicModel::writeDynamicCFile(const string &basename) const
 {
   filesystem::create_directories(basename + "/model/src");
   string filename = basename + "/model/src/dynamic.c";
@@ -1780,7 +1780,7 @@ DynamicModel::writeDynamicCFile(const string &basename, int order) const
                   << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
                   << "{" << endl
                   << "  /* Check that no derivatives of higher order than computed are being requested */" << endl
-                  << "  if (nlhs > " << order + 1 << ")" << endl
+                  << "  if (nlhs > " << computed_derivs_order << ")" << endl
                   << R"(    mexErrMsgTxt("Derivatives of higher order than computed have been requested");)" << endl
                   << "  /* Create a pointer to the input matrix y. */" << endl
                   << "  double *y = mxGetPr(prhs[0]);" << endl
@@ -3080,7 +3080,7 @@ DynamicModel::includeExcludeEquations(const string &eqs, bool exclude_eqs)
 }
 
 void
-DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const
+DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const
 {
   /* Writing initialisation for M_.lead_lag_incidence matrix
      M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -3689,7 +3689,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   // Use -1 if the derivatives have not been computed
   output << modstruct << (julia ? "nnzderivatives" : "NNZDerivatives") << " = [";
   for (int i = 1; i < static_cast<int>(NNZDerivatives.size()); i++)
-    output << (i > order ? -1 : NNZDerivatives[i]) << "; ";
+    output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
   output << "];" << endl;
 
   // Write Pac Model Consistent Expectation parameter info
@@ -4959,11 +4959,8 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
   computeDerivatives(derivsOrder, vars);
 
   if (derivsOrder > 1)
-    {
-      hessian_computed = true;
-      for (const auto &[indices, d2] : derivatives[2])
-        nonzero_hessian_eqs.insert(indices[0]);
-    }
+    for (const auto &[indices, d2] : derivatives[2])
+      nonzero_hessian_eqs.insert(indices[0]);
 
   if (paramsDerivsOrder > 0)
     {
@@ -5442,7 +5439,7 @@ DynamicModel::collectBlockVariables()
 }
 
 void
-DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, int order, bool julia) const
+DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
 {
   if (block && bytecode)
     writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
@@ -5456,7 +5453,7 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_d
     writeSparseDynamicMFile(basename);
   else if (use_dll)
     {
-      writeDynamicCFile(basename, order);
+      writeDynamicCFile(basename);
       compileDll(basename, "dynamic", mexext, matlabroot, dynareroot);
     }
   else if (julia)
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 487017a6..8eea2c3c 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -99,8 +99,6 @@ private:
 
   //! Nonzero equations in the Hessian
   set<int> nonzero_hessian_eqs;
-  //! Whether the hessian has actually been computed
-  bool hessian_computed{false};
 
   //! Number of columns of dynamic jacobian
   /*! Set by computeDerivID()s and computeDynJacobianCols() */
@@ -123,7 +121,7 @@ private:
   void writeDynamicJuliaFile(const string &dynamic_basename) const;
   //! Writes dynamic model file (C version)
   /*! \todo add third derivatives handling */
-  void writeDynamicCFile(const string &basename, int order) const;
+  void writeDynamicCFile(const string &basename) const;
   //! Writes dynamic model file when SparseDLL option is on
   void writeSparseDynamicMFile(const string &basename) const;
   //! Writes the dynamic model equations and its derivatives
@@ -327,7 +325,7 @@ public:
   void computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
                      const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool linear_decomposition);
   //! Writes model initialization and lead/lag incidence matrix to output
-  void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
+  void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const;
 
   //! Write JSON AST
   void writeJsonAST(ostream &output) const;
@@ -363,7 +361,7 @@ public:
   inline bool
   isHessianComputed() const
   {
-    return hessian_computed;
+    return computed_derivs_order >= 2;
   }
   //! Returns equations that have non-zero second derivatives
   inline set<int>
@@ -415,7 +413,7 @@ public:
   void Write_Inf_To_Bin_File_Block(const string &basename,
                                    int num, int &u_count_int, bool &file_open, bool is_two_boundaries, bool linear_decomposition) const;
   //! Writes dynamic model file
-  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, int order, bool julia) const;
+  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const;
   //! Writes file containing parameters derivatives
   void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 47707d04..20c27c37 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2006-2019 Dynare Team
+ * Copyright © 2006-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1018,8 +1018,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   if (dynamic_model.equation_number() > 0)
     {
       if (linear_decomposition)
-        non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
-      dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
+        non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
+      dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
       if (!no_static)
         static_model.writeOutput(mOutputFile, block);
     }
@@ -1126,11 +1126,11 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
           if (linear_decomposition)
             {
-              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mexext, matlabroot, dynareroot, mod_file_struct.order_option, false);
+              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mexext, matlabroot, dynareroot, false);
               non_linear_equations_dynamic_model.writeParamsDerivativesFile(basename, false);
             }
 
-          dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mexext, matlabroot, dynareroot, mod_file_struct.order_option, false);
+          dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mexext, matlabroot, dynareroot, false);
 
           dynamic_model.writeParamsDerivativesFile(basename, false);
 
@@ -1245,7 +1245,6 @@ ModFile::writeExternalFilesJulia(const string &basename) const
   if (dynamic_model.equation_number() > 0)
     {
       dynamic_model.writeOutput(jlOutputFile, basename, false, false, false, false,
-                                mod_file_struct.order_option,
                                 mod_file_struct.estimation_present, false, true);
       if (!no_static)
         {
@@ -1253,7 +1252,7 @@ ModFile::writeExternalFilesJulia(const string &basename) const
           static_model.writeParamsDerivativesFile(basename, true);
         }
       dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll,
-                                     "", {}, {}, mod_file_struct.order_option, true);
+                                     "", {}, {}, true);
       dynamic_model.writeParamsDerivativesFile(basename, true);
     }
   steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, true);
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 3695cc86..f6ec39b3 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -112,6 +112,7 @@ ModelTree::ModelTree(const ModelTree &m) :
   equations_lineno{m.equations_lineno},
   equation_tags{m.equation_tags},
   equation_tags_xref{m.equation_tags_xref},
+  computed_derivs_order{m.computed_derivs_order},
   NNZDerivatives{m.NNZDerivatives},
   equation_reordered{m.equation_reordered},
   variable_reordered{m.variable_reordered},
@@ -138,6 +139,7 @@ ModelTree::operator=(const ModelTree &m)
   aux_equations.clear();
   equation_tags = m.equation_tags;
   equation_tags_xref = m.equation_tags_xref;
+  computed_derivs_order = m.computed_derivs_order;
   NNZDerivatives = m.NNZDerivatives;
 
   derivatives.clear();
@@ -1245,6 +1247,8 @@ ModelTree::computeDerivatives(int order, const set<int> &vars)
 {
   assert(order >= 1);
 
+  computed_derivs_order = order;
+
   // Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
   derivatives.resize(max(static_cast<size_t>(order+1), derivatives.size()));
   NNZDerivatives.resize(max(static_cast<size_t>(order+1), NNZDerivatives.size()), 0);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 040b2696..5d7ad1cb 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -106,6 +106,9 @@ protected:
       tests/optimal_policy/nk_ramsey_expectation.mod */
   vector<BinaryOpNode *> aux_equations;
 
+  //! Maximum order at which (endogenous) derivatives have been computed
+  int computed_derivs_order{0};
+
   //! Stores derivatives
   /*! Index 0 is not used, index 1 contains first derivatives, ...
      For each derivation order, stores a map whose key is a vector of integer: the