diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 074e9af2..52b62373 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5438,6 +5438,22 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
       hessian1_output << ";" << endl;
 
       i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                          << "=rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",4"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+          i++;
+        }
     }
 
   i = 1;
@@ -5465,6 +5481,24 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
       third_derivs_output << ";" << endl;
 
       i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                              << "=gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+          i++;
+        }
     }
 
   i = 1;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index db780b53..dcdb5059 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -2033,73 +2033,51 @@ ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, Expr
 void
 ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
 {
-  if (!(paramsDerivsOrder == 1 || paramsDerivsOrder == 2))
-    return;
+  assert(paramsDerivsOrder >= 1);
+
   set<int> deriv_id_set;
   addAllParamDerivId(deriv_id_set);
 
+  // First-order derivatives w.r.t. params
   for (int param : deriv_id_set)
     {
       for (int eq = 0; eq < (int) equations.size(); eq++)
         {
-          expr_t d1 = equations[eq]->getDerivative(param);
-          if (d1 == Zero)
+          expr_t d = equations[eq]->getDerivative(param);
+          if (d == Zero)
             continue;
-          params_derivatives[{ 0, 1 }][{ eq, param }] = d1;
+          params_derivatives[{ 0, 1 }][{ eq, param }] = d;
         }
 
-      if (paramsDerivsOrder == 2)
-        for (const auto &it : params_derivatives[{ 0, 1 }])
+      for (int endoOrd = 1; endoOrd < (int) derivatives.size(); endoOrd++)
+        for (const auto &it : derivatives[endoOrd])
           {
-            int eq, param1;
-            tie(eq, param1) = vectorToTuple<2>(it.first);
-            expr_t d1 = it.second;
-
-            expr_t d2 = d1->getDerivative(param);
-            if (d2 == Zero)
+            expr_t d = it.second->getDerivative(param);
+            if (d == Zero)
               continue;
-            params_derivatives[{ 0, 2 }][{ eq, param1, param }] = d2;
+            vector<int> indices{it.first};
+            indices.push_back(param);
+            params_derivatives[{ endoOrd, 1 }][indices] = d;
           }
-
-      for (const auto &it : derivatives[1])
-        {
-          int eq, var;
-          tie(eq, var) = vectorToTuple<2>(it.first);
-          expr_t d1 = it.second;
-
-          expr_t d2 = d1->getDerivative(param);
-          if (d2 == Zero)
-            continue;
-          params_derivatives[{ 1, 1 }][{ eq, var, param }] = d2;
-        }
-
-      if (paramsDerivsOrder == 2)
-        {
-          for (const auto &it : params_derivatives[{ 1, 1 }])
-            {
-              int eq, var, param1;
-              tie(eq, var, param1) = vectorToTuple<3>(it.first);
-              expr_t d1 = it.second;
-
-              expr_t d2 = d1->getDerivative(param);
-              if (d2 == Zero)
-                continue;
-              params_derivatives[{ 1, 2 }][{ eq, var, param1, param }] = d2;
-            }
-
-          for (const auto &it : derivatives[2])
-            {
-              int eq, var1, var2;
-              tie(eq, var1, var2) = vectorToTuple<3>(it.first);
-              expr_t d1 = it.second;
-
-              expr_t d2 = d1->getDerivative(param);
-              if (d2 == Zero)
-                continue;
-              params_derivatives[{ 2, 1 }][{ eq, var1, var2, param }] = d2;
-            }
-        }
     }
+
+  // Higher-order derivatives w.r.t. parameters
+  for (int endoOrd = 0; endoOrd < (int) derivatives.size(); endoOrd++)
+    for (int paramOrd = 2; paramOrd <= paramsDerivsOrder; paramOrd++)
+      for (const auto &it : params_derivatives[{ endoOrd, paramOrd-1 }])
+        for (int param : deriv_id_set)
+          {
+            if (it.first.back() > param)
+              continue;
+
+            expr_t d = it.second->getDerivative(param);
+            if (d == Zero)
+              continue;
+            vector<int> indices{it.first};
+            indices.push_back(param);
+            // At this point, indices of both endogenous and parameters are sorted in non-decreasing order
+            params_derivatives[{ endoOrd, paramOrd }][indices] = d;
+          }
 }
 
 void
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 8612ded8..aec7fe64 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -105,7 +105,8 @@ protected:
   derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian
   differentiated twice w.r.t. to parameters.
   In inner maps, the vector of integers consists of: the equation index, then
-  the derivation IDs of endogenous (in non-decreasing order), then the IDs of parameters */
+  the derivation IDs of endogenous (in non-decreasing order),
+  then the IDs of parameters (in non-decreasing order)*/
   map<pair<int,int>, map<vector<int>, expr_t>> params_derivatives;
 
   //! Storage for temporary terms in block/bytecode mode
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 64229a76..b2d9641a 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -2710,6 +2710,22 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
       hessian1_output << ";" << endl;
 
       i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                          << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                          << "=rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",4"
+                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+          i++;
+        }
     }
 
   i = 1;
@@ -2737,6 +2753,24 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
       third_derivs_output << ";" << endl;
 
       i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                              << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                              << "=gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
+                              << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+          i++;
+        }
     }
 
   i = 1;