diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 7db188b9570bf54ae72a5dabc8c40694b6b2e829..61295f9065d9fc37f3992d033482d18b4a05d203 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5386,333 +5386,33 @@ DynamicModel::writeJsonDynamicModelInfo(ostream &output) const
 void
 DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
 {
-  ostringstream model_local_vars_output; // Used for storing model local vars
-  vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
-
-  deriv_node_temp_terms_t tef_terms;
-  temporary_terms_t temp_term_union;
-
-  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
-
-  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
-  d_output[0] << ", ";
-  writeJsonModelEquations(d_output[0], true);
-
-  int ncols{getJacobianColsNbr()};
-  for (size_t i = 1; i < derivatives.size(); i++)
-    {
-      string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
-      writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
-      temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
-      d_output[i] << R"(, ")" << matrix_name  << R"(": {)"
-                  << R"(  "nrows": )" << equations.size()
-                  << R"(, "ncols": )" << ncols
-                  << R"(, "entries": [)";
-
-      for (bool printed_something{false};
-           const auto &[vidx, d] : derivatives[i])
-        {
-          if (exchange(printed_something, true))
-            d_output[i] << ", ";
-
-          int eq = vidx[0];
-
-          int col_idx = 0;
-          for (size_t j = 1; j < vidx.size(); j++)
-            {
-              col_idx *= getJacobianColsNbr();
-              col_idx += getJacobianCol(vidx[j]);
-            }
-
-          if (writeDetails)
-            d_output[i] << R"({"eq": )" << eq + 1;
-          else
-            d_output[i] << R"({"row": )" << eq + 1;
-
-          d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
-
-          if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
-            {
-              int col_idx_sym = getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1]);
-              d_output[i] << ", " << col_idx_sym + 1;
-            }
-          if (i > 1)
-            d_output[i] << "]";
-
-          if (writeDetails)
-            for (size_t j = 1; j < vidx.size(); j++)
-              d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << getNameByDerivID(vidx[j]) << R"(")"
-                          << R"(, "shift)" << (i > 1 ? to_string(j) : "") << R"(": )" << getLagByDerivID(vidx[j]);
-
-          d_output[i] << R"(, "val": ")";
-          d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
-          d_output[i] << R"("})" << endl;
-        }
-      d_output[i] << "]}";
-
-      ncols *= getJacobianColsNbr();
-    }
+  auto [mlv_output, d_output] { writeJsonComputingPassOutputHelper<true>(writeDetails) };
 
   if (writeDetails)
     output << R"("dynamic_model": {)";
   else
     output << R"("dynamic_model_simple": {)";
-  output << model_local_vars_output.str();
+  output << mlv_output.str();
   for (const auto &it : d_output)
     output << ", " << it.str();
   output << "}";
 }
 
 void
-DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
+DynamicModel::writeJsonParamsDerivatives(ostream &output, bool writeDetails) const
 {
   if (!params_derivatives.size())
     return;
 
-  ostringstream model_local_vars_output; // Used for storing model local vars
-  ostringstream model_output; // Used for storing model temp vars and equations
-  ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters
-  ostringstream gp_output; // 1st deriv. of Jacobian w.r.t. parameters
-  ostringstream rpp_output; // 2nd deriv of residuals w.r.t. parameters
-  ostringstream gpp_output; // 2nd deriv of Jacobian w.r.t. parameters
-  ostringstream hp_output; // 1st deriv. of Hessian w.r.t. parameters
-  ostringstream g3p_output; // 1st deriv. of 3rd deriv. matrix w.r.t. parameters
-
-  deriv_node_temp_terms_t tef_terms;
-  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
-
-  temporary_terms_t temp_term_union;
-  for (const auto &it : params_derivs_temporary_terms)
-    writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, "all");
-
-  rp_output << R"("deriv_wrt_params": {)"
-            << R"(  "neqs": )" << equations.size()
-            << R"(, "nparamcols": )" << symbol_table.param_nbr()
-            << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 0, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        rp_output << ", ";
-
-      auto [eq, param] = vectorToTuple<2>(vidx);
-
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        rp_output << R"({"eq": )" << eq + 1;
-      else
-        rp_output << R"({"row": )" << eq + 1;
-
-      rp_output << R"(, "param_col": )" << param_col;
-
-      if (writeDetails)
-        rp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      rp_output << R"(, "val": ")";
-      d->writeJsonOutput(rp_output, temp_term_union, tef_terms);
-      rp_output << R"("})" << endl;
-    }
-  rp_output << "]}";
-
-  gp_output << R"("deriv_jacobian_wrt_params": {)"
-            << R"(  "neqs": )" << equations.size()
-            << R"(, "nvarcols": )" << getJacobianColsNbr()
-            << R"(, "nparamcols": )" << symbol_table.param_nbr()
-            << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 1, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        gp_output << ", ";
-
-      auto [eq, var, param] = vectorToTuple<3>(vidx);
-
-      int var_col = getJacobianCol(var) + 1;
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        gp_output << R"({"eq": )" << eq + 1;
-      else
-        gp_output << R"({"row": )" << eq + 1;
-
-      gp_output << R"(, "var_col": )" << var_col
-                << R"(, "param_col": )" << param_col;
-
-      if (writeDetails)
-        gp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
-                  << R"(, "lag": )" << getLagByDerivID(var)
-                  << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      gp_output << R"(, "val": ")";
-      d->writeJsonOutput(gp_output, temp_term_union, tef_terms);
-      gp_output << R"("})" << endl;
-    }
-  gp_output << "]}";
-
-  rpp_output << R"("second_deriv_residuals_wrt_params": {)"
-             << R"(  "nrows": )" << equations.size()
-             << R"(, "nparam1cols": )" << symbol_table.param_nbr()
-             << R"(, "nparam2cols": )" << symbol_table.param_nbr()
-             << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 0, 2 })->second)
-    {
-      if (exchange(printed_something, true))
-        rpp_output << ", ";
-
-      auto [eq, param1, param2] = vectorToTuple<3>(vidx);
-
-      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
-      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
-
-      if (writeDetails)
-        rpp_output << R"({"eq": )" << eq + 1;
-      else
-        rpp_output << R"({"row": )" << eq + 1;
-      rpp_output << R"(, "param1_col": )" << param1_col
-                 << R"(, "param2_col": )" << param2_col;
-
-      if (writeDetails)
-        rpp_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
-                   << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
-
-      rpp_output << R"(, "val": ")";
-      d->writeJsonOutput(rpp_output, temp_term_union, tef_terms);
-      rpp_output << R"("})" << endl;
-    }
-  rpp_output << "]}";
-
-  gpp_output << R"("second_deriv_jacobian_wrt_params": {)"
-             << R"(  "neqs": )" << equations.size()
-             << R"(, "nvarcols": )" << getJacobianColsNbr()
-             << R"(, "nparam1cols": )" << symbol_table.param_nbr()
-             << R"(, "nparam2cols": )" << symbol_table.param_nbr()
-             << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 1, 2 })->second)
-    {
-      if (exchange(printed_something, true))
-        gpp_output << ", ";
-
-      auto [eq, var, param1, param2] = vectorToTuple<4>(vidx);
-
-      int var_col = getJacobianCol(var) + 1;
-      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
-      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
-
-      if (writeDetails)
-        gpp_output << R"({"eq": )" << eq + 1;
-      else
-        gpp_output << R"({"row": )" << eq + 1;
-
-      gpp_output << R"(, "var_col": )" << var_col
-                 << R"(, "param1_col": )" << param1_col
-                 << R"(, "param2_col": )" << param2_col;
-
-      if (writeDetails)
-        gpp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
-                   << R"(, "lag": )" << getLagByDerivID(var)
-                   << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
-                   << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
-
-      gpp_output << R"(, "val": ")";
-      d->writeJsonOutput(gpp_output, temp_term_union, tef_terms);
-      gpp_output << R"("})" << endl;
-    }
-  gpp_output << "]}" << endl;
-
-  hp_output << R"("derivative_hessian_wrt_params": {)"
-            << R"(  "neqs": )" << equations.size()
-            << R"(, "nvar1cols": )" << getJacobianColsNbr()
-            << R"(, "nvar2cols": )" << getJacobianColsNbr()
-            << R"(, "nparamcols": )" << symbol_table.param_nbr()
-            << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 2, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        hp_output << ", ";
-
-      auto [eq, var1, var2, param] = vectorToTuple<4>(vidx);
-
-      int var1_col = getJacobianCol(var1) + 1;
-      int var2_col = getJacobianCol(var2) + 1;
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        hp_output << R"({"eq": )" << eq + 1;
-      else
-        hp_output << R"({"row": )" << eq + 1;
-
-      hp_output << R"(, "var1_col": )" << var1_col
-                << R"(, "var2_col": )" << var2_col
-                << R"(, "param_col": )" << param_col;
-
-      if (writeDetails)
-        hp_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
-                  << R"(, "lag1": )" << getLagByDerivID(var1)
-                  << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")"
-                  << R"(, "lag2": )" << getLagByDerivID(var2)
-                  << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      hp_output << R"(, "val": ")";
-      d->writeJsonOutput(hp_output, temp_term_union, tef_terms);
-      hp_output << R"("})" << endl;
-    }
-  hp_output << "]}" << endl;
-
-  g3p_output << R"("derivative_g3_wrt_params": {)"
-             << R"(  "neqs": )" << equations.size()
-             << R"(, "nvar1cols": )" << getJacobianColsNbr()
-             << R"(, "nvar2cols": )" << getJacobianColsNbr()
-             << R"(, "nvar3cols": )" << getJacobianColsNbr()
-             << R"(, "nparamcols": )" << symbol_table.param_nbr()
-             << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 3, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        g3p_output << ", ";
-
-      auto [eq, var1, var2, var3, param] = vectorToTuple<5>(vidx);
-
-      int var1_col = getJacobianCol(var1) + 1;
-      int var2_col = getJacobianCol(var2) + 1;
-      int var3_col = getJacobianCol(var3) + 1;
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        g3p_output << R"({"eq": )" << eq + 1;
-      else
-        g3p_output << R"({"row": )" << eq + 1;
-
-      g3p_output << R"(, "var1_col": )" << var1_col
-                 << R"(, "var2_col": )" << var2_col
-                 << R"(, "var3_col": )" << var3_col
-                 << R"(, "param_col": )" << param_col;
-
-      if (writeDetails)
-        g3p_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
-                   << R"(, "lag1": )" << getLagByDerivID(var1)
-                   << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")"
-                   << R"(, "lag2": )" << getLagByDerivID(var2)
-                   << R"(, "var3": ")" << getNameByDerivID(var3) << R"(")"
-                   << R"(, "lag3": )" << getLagByDerivID(var3)
-                   << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      g3p_output << R"(, "val": ")";
-      d->writeJsonOutput(g3p_output, temp_term_union, tef_terms);
-      g3p_output << R"("})" << endl;
-    }
-  g3p_output << "]}" << endl;
+  auto [mlv_output, tt_output, rp_output, gp_output, rpp_output, gpp_output, hp_output, g3p_output]
+    { writeJsonParamsDerivativesHelper<true>(writeDetails) };
 
   if (writeDetails)
     output << R"("dynamic_model_params_derivative": {)";
   else
     output << R"("dynamic_model_params_derivatives_simple": {)";
-  output << model_local_vars_output.str()
-         << ", " << model_output.str()
+  output << mlv_output.str()
+         << ", " << tt_output.str()
          << ", " << rp_output.str()
          << ", " << gp_output.str()
          << ", " << rpp_output.str()
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index e800a7467e47818402bddc91df33b3e076ad9f87..972128f9e92d7927badae3332f4c46f9dbd8c05f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -342,8 +342,8 @@ public:
   //! Write JSON Output representation of dynamic model after computing pass
   void writeJsonComputingPassOutput(ostream &output, bool writeDetails) const;
 
-  //! Write JSON prams derivatives file
-  void writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const;
+  //! Write JSON params derivatives
+  void writeJsonParamsDerivatives(ostream &output, bool writeDetails) const;
 
   //! Write cross reference output if the xref maps have been filed
   void writeJsonXrefs(ostream &output) const;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index d21baa39c6c15358ffb95b6c905aee757a312931..08cb38d377141b25633c7f58ba11c21645ddcbdf 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1351,12 +1351,12 @@ ModFile::writeJsonComputingPassOutput(const string &basename, JsonFileOutputType
   dynamic_model.writeJsonComputingPassOutput(dynamic_output, !jsonderivsimple);
   dynamic_output << "}";
 
-  static_model.writeJsonParamsDerivativesFile(tmp_out, !jsonderivsimple);
+  static_model.writeJsonParamsDerivatives(tmp_out, !jsonderivsimple);
   if (!tmp_out.str().empty())
     static_paramsd_output << "{" << tmp_out.str() << "}" << endl;
 
   tmp_out.str("");
-  dynamic_model.writeJsonParamsDerivativesFile(tmp_out, !jsonderivsimple);
+  dynamic_model.writeJsonParamsDerivatives(tmp_out, !jsonderivsimple);
   if (!tmp_out.str().empty())
     dynamic_paramsd_output << "{" << tmp_out.str() << "}" << endl;
 
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index e7c92264bee7635d9fcc8629d56f6e7ed79a91e2..cdcdbd0183abcdbbb659812891285698391cf7a2 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -269,6 +269,18 @@ protected:
   tuple<ostringstream, ostringstream, ostringstream, ostringstream,
         ostringstream, ostringstream, ostringstream> writeParamsDerivativesFileHelper() const;
 
+  /* Helper for writing JSON output for residuals and derivatives.
+     Returns mlv and derivatives output at each derivation order. */
+  template<bool dynamic>
+  pair<ostringstream, vector<ostringstream>> writeJsonComputingPassOutputHelper(bool writeDetails) const;
+
+  /* Helper for writing JSON derivatives w.r.t. parameters.
+     Returns { mlv, tt, rp, gp, rpp, gpp, hp, g3p }.
+     g3p is empty if requesting a static output type. */
+  template<bool dynamic>
+  tuple<ostringstream, ostringstream, ostringstream, ostringstream, ostringstream,
+        ostringstream, ostringstream, ostringstream> writeJsonParamsDerivativesHelper(bool writeDetails) const;
+
   //! Writes JSON model equations
   //! if residuals = true, we are writing the dynamic/static model.
   //! Otherwise, just the model equations (with line numbers, no tmp terms)
@@ -992,4 +1004,341 @@ ModelTree::writeParamsDerivativesFileHelper() const
     move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output) };
 }
 
+template<bool dynamic>
+pair<ostringstream, vector<ostringstream>>
+ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
+{
+  ostringstream mlv_output; // Used for storing model local vars
+  vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
+
+  deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temp_term_union;
+
+  writeJsonModelLocalVariables(mlv_output, true, tef_terms);
+
+  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
+  d_output[0] << ", ";
+  writeJsonModelEquations(d_output[0], true);
+
+  int ncols { getJacobianColsNbr() };
+  for (size_t i {1}; i < derivatives.size(); i++)
+    {
+      string matrix_name { i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative"};
+      writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
+      temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
+      d_output[i] << R"(, ")" << matrix_name  << R"(": {)"
+                  << R"(  "nrows": )" << equations.size()
+                  << R"(, "ncols": )" << ncols
+                  << R"(, "entries": [)";
+
+      for (bool printed_something {false};
+           const auto &[vidx, d] : derivatives[i])
+        {
+          if (exchange(printed_something, true))
+            d_output[i] << ", ";
+
+          int eq { vidx[0] };
+
+          int col_idx {0};
+          for (size_t j {1}; j < vidx.size(); j++)
+            {
+              col_idx *= getJacobianColsNbr();
+              col_idx += getJacobianCol(vidx[j]);
+            }
+
+          if (writeDetails)
+            d_output[i] << R"({"eq": )" << eq + 1;
+          else
+            d_output[i] << R"({"row": )" << eq + 1;
+
+          d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
+
+          if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
+            {
+              int col_idx_sym { getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])};
+              d_output[i] << ", " << col_idx_sym + 1;
+            }
+          if (i > 1)
+            d_output[i] << "]";
+
+          if (writeDetails)
+            for (size_t j = 1; j < vidx.size(); j++)
+              {
+                d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << getNameByDerivID(vidx[j]) << R"(")";
+                if constexpr(dynamic)
+                   d_output[i] << R"(, "shift)" << (i > 1 ? to_string(j) : "") << R"(": )" << getLagByDerivID(vidx[j]);
+              }
+
+          d_output[i] << R"(, "val": ")";
+          d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
+          d_output[i] << R"("})" << endl;
+        }
+      d_output[i] << "]}";
+
+      ncols *= getJacobianColsNbr();
+    }
+
+  return { move(mlv_output), move(d_output) };
+}
+
+template<bool dynamic>
+tuple<ostringstream, ostringstream, ostringstream, ostringstream,
+      ostringstream, ostringstream, ostringstream, ostringstream>
+ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
+{
+  ostringstream mlv_output; // Used for storing model local vars
+  ostringstream tt_output; // Used for storing model temp vars and equations
+  ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters
+  ostringstream gp_output; // 1st deriv. of Jacobian w.r.t. parameters
+  ostringstream rpp_output; // 2nd deriv of residuals w.r.t. parameters
+  ostringstream gpp_output; // 2nd deriv of Jacobian w.r.t. parameters
+  ostringstream hp_output; // 1st deriv. of Hessian w.r.t. parameters
+  ostringstream g3p_output; // 1st deriv. of 3rd deriv. matrix w.r.t. parameters
+
+  deriv_node_temp_terms_t tef_terms;
+  writeJsonModelLocalVariables(mlv_output, true, tef_terms);
+
+  temporary_terms_t temp_term_union;
+  for (const auto &[order, tts] : params_derivs_temporary_terms)
+    writeJsonTemporaryTerms(tts, temp_term_union, tt_output, tef_terms, "all");
+
+  rp_output << R"("deriv_wrt_params": {)"
+            << R"(  "neqs": )" << equations.size()
+            << R"(, "nparamcols": )" << symbol_table.param_nbr()
+            << R"(, "entries": [)";
+  for (bool printed_something {false};
+       const auto &[vidx, d] : params_derivatives.find({ 0, 1 })->second)
+    {
+      if (exchange(printed_something, true))
+        rp_output << ", ";
+
+      auto [eq, param] { vectorToTuple<2>(vidx) };
+
+      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
+
+      if (writeDetails)
+        rp_output << R"({"eq": )" << eq + 1;
+      else
+        rp_output << R"({"row": )" << eq + 1;
+
+      rp_output << R"(, "param_col": )" << param_col;
+
+      if (writeDetails)
+        rp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
+
+      rp_output << R"(, "val": ")";
+      d->writeJsonOutput(rp_output, temp_term_union, tef_terms);
+      rp_output << R"("})" << endl;
+    }
+  rp_output << "]}";
+
+  gp_output << R"("deriv_jacobian_wrt_params": {)"
+            << R"(  "neqs": )" << equations.size()
+            << R"(, "nvarcols": )" << getJacobianColsNbr()
+            << R"(, "nparamcols": )" << symbol_table.param_nbr()
+            << R"(, "entries": [)";
+  for (bool printed_something {false};
+       const auto &[vidx, d] : params_derivatives.find({ 1, 1 })->second)
+    {
+      if (exchange(printed_something, true))
+        gp_output << ", ";
+
+      auto [eq, var, param] { vectorToTuple<3>(vidx) };
+
+      int var_col { getJacobianCol(var) + 1 };
+      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
+
+      if (writeDetails)
+        gp_output << R"({"eq": )" << eq + 1;
+      else
+        gp_output << R"({"row": )" << eq + 1;
+
+      gp_output << R"(, "var_col": )" << var_col
+                << R"(, "param_col": )" << param_col;
+
+      if (writeDetails)
+        {
+          gp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")";
+          if constexpr(dynamic)
+            gp_output << R"(, "lag": )" << getLagByDerivID(var);
+          gp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
+        }
+
+      gp_output << R"(, "val": ")";
+      d->writeJsonOutput(gp_output, temp_term_union, tef_terms);
+      gp_output << R"("})" << endl;
+    }
+  gp_output << "]}";
+
+  rpp_output << R"("second_deriv_residuals_wrt_params": {)"
+             << R"(  "nrows": )" << equations.size()
+             << R"(, "nparam1cols": )" << symbol_table.param_nbr()
+             << R"(, "nparam2cols": )" << symbol_table.param_nbr()
+             << R"(, "entries": [)";
+  for (bool printed_something {false};
+       const auto &[vidx, d] : params_derivatives.find({ 0, 2 })->second)
+    {
+      if (exchange(printed_something, true))
+        rpp_output << ", ";
+
+      auto [eq, param1, param2] { vectorToTuple<3>(vidx) };
+
+      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
+      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
+
+      if (writeDetails)
+        rpp_output << R"({"eq": )" << eq + 1;
+      else
+        rpp_output << R"({"row": )" << eq + 1;
+      rpp_output << R"(, "param1_col": )" << param1_col
+                 << R"(, "param2_col": )" << param2_col;
+
+      if (writeDetails)
+        rpp_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
+                   << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
+
+      rpp_output << R"(, "val": ")";
+      d->writeJsonOutput(rpp_output, temp_term_union, tef_terms);
+      rpp_output << R"("})" << endl;
+    }
+  rpp_output << "]}";
+
+  gpp_output << R"("second_deriv_jacobian_wrt_params": {)"
+             << R"(  "neqs": )" << equations.size()
+             << R"(, "nvarcols": )" << getJacobianColsNbr()
+             << R"(, "nparam1cols": )" << symbol_table.param_nbr()
+             << R"(, "nparam2cols": )" << symbol_table.param_nbr()
+             << R"(, "entries": [)";
+  for (bool printed_something {false};
+       const auto &[vidx, d] : params_derivatives.find({ 1, 2 })->second)
+    {
+      if (exchange(printed_something, true))
+        gpp_output << ", ";
+
+      auto [eq, var, param1, param2] { vectorToTuple<4>(vidx) };
+
+      int var_col { getJacobianCol(var) + 1 };
+      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
+      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
+
+      if (writeDetails)
+        gpp_output << R"({"eq": )" << eq + 1;
+      else
+        gpp_output << R"({"row": )" << eq + 1;
+
+      gpp_output << R"(, "var_col": )" << var_col
+                 << R"(, "param1_col": )" << param1_col
+                 << R"(, "param2_col": )" << param2_col;
+
+      if (writeDetails)
+        {
+          gpp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")";
+          if constexpr(dynamic)
+            gpp_output << R"(, "lag": )" << getLagByDerivID(var);
+          gpp_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
+                     << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
+        }
+
+      gpp_output << R"(, "val": ")";
+      d->writeJsonOutput(gpp_output, temp_term_union, tef_terms);
+      gpp_output << R"("})" << endl;
+    }
+  gpp_output << "]}" << endl;
+
+  hp_output << R"("derivative_hessian_wrt_params": {)"
+            << R"(  "neqs": )" << equations.size()
+            << R"(, "nvar1cols": )" << getJacobianColsNbr()
+            << R"(, "nvar2cols": )" << getJacobianColsNbr()
+            << R"(, "nparamcols": )" << symbol_table.param_nbr()
+            << R"(, "entries": [)";
+  for (bool printed_something {false};
+       const auto &[vidx, d] : params_derivatives.find({ 2, 1 })->second)
+    {
+      if (exchange(printed_something, true))
+        hp_output << ", ";
+
+      auto [eq, var1, var2, param] { vectorToTuple<4>(vidx) };
+
+      int var1_col { getJacobianCol(var1) + 1 };
+      int var2_col { getJacobianCol(var2) + 1 };
+      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
+
+      if (writeDetails)
+        hp_output << R"({"eq": )" << eq + 1;
+      else
+        hp_output << R"({"row": )" << eq + 1;
+
+      hp_output << R"(, "var1_col": )" << var1_col
+                << R"(, "var2_col": )" << var2_col
+                << R"(, "param_col": )" << param_col;
+
+      if (writeDetails)
+        {
+          hp_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")";
+          if constexpr(dynamic)
+            hp_output << R"(, "lag1": )" << getLagByDerivID(var1);
+          hp_output << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")";
+          if constexpr(dynamic)
+            hp_output << R"(, "lag2": )" << getLagByDerivID(var2);
+          hp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
+        }
+
+      hp_output << R"(, "val": ")";
+      d->writeJsonOutput(hp_output, temp_term_union, tef_terms);
+      hp_output << R"("})" << endl;
+    }
+  hp_output << "]}" << endl;
+
+  if constexpr(dynamic)
+    {
+      g3p_output << R"("derivative_g3_wrt_params": {)"
+                 << R"(  "neqs": )" << equations.size()
+                 << R"(, "nvar1cols": )" << getJacobianColsNbr()
+                 << R"(, "nvar2cols": )" << getJacobianColsNbr()
+                 << R"(, "nvar3cols": )" << getJacobianColsNbr()
+                 << R"(, "nparamcols": )" << symbol_table.param_nbr()
+                 << R"(, "entries": [)";
+      for (bool printed_something {false};
+           const auto &[vidx, d] : params_derivatives.find({ 3, 1 })->second)
+        {
+          if (exchange(printed_something, true))
+            g3p_output << ", ";
+
+          auto [eq, var1, var2, var3, param] { vectorToTuple<5>(vidx) };
+
+          int var1_col { getJacobianCol(var1) + 1 };
+          int var2_col { getJacobianCol(var2) + 1 };
+          int var3_col { getJacobianCol(var3) + 1 };
+          int param_col { getTypeSpecificIDByDerivID(param) + 1 };
+
+          if (writeDetails)
+            g3p_output << R"({"eq": )" << eq + 1;
+          else
+            g3p_output << R"({"row": )" << eq + 1;
+
+          g3p_output << R"(, "var1_col": )" << var1_col + 1
+                     << R"(, "var2_col": )" << var2_col + 1
+                     << R"(, "var3_col": )" << var3_col + 1
+                     << R"(, "param_col": )" << param_col + 1;
+
+          if (writeDetails)
+            g3p_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
+                       << R"(, "lag1": )" << getLagByDerivID(var1)
+                       << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")"
+                       << R"(, "lag2": )" << getLagByDerivID(var2)
+                       << R"(, "var3": ")" << getNameByDerivID(var3) << R"(")"
+                       << R"(, "lag3": )" << getLagByDerivID(var3)
+                       << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
+
+          g3p_output << R"(, "val": ")";
+          d->writeJsonOutput(g3p_output, temp_term_union, tef_terms);
+          g3p_output << R"("})" << endl;
+        }
+      g3p_output << "]}" << endl;
+    }
+
+  return { move(mlv_output), move(tt_output), move(rp_output), move(gp_output),
+    move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output) };
+}
+
 #endif
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index b8eb3913b390f25f880cf5dcbeccf2fe26bebebf..0d7af888a143643425330991619750df142c2be8 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1849,285 +1849,38 @@ StaticModel::writeJsonOutput(ostream &output) const
 void
 StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
 {
-  ostringstream model_local_vars_output; // Used for storing model local vars
-  vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
-
-  deriv_node_temp_terms_t tef_terms;
-  temporary_terms_t temp_term_union;
-
-  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
-
-  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
-  d_output[0] << ", ";
-  writeJsonModelEquations(d_output[0], true);
-
-  int ncols = symbol_table.endo_nbr();
-  for (size_t i = 1; i < derivatives.size(); i++)
-    {
-      string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
-      writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
-      temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
-      d_output[i] << R"(, ")" << matrix_name  << R"(": {)"
-                  << R"(  "nrows": )" << equations.size()
-                  << R"(, "ncols": )" << ncols
-                  << R"(, "entries": [)";
-
-      for (bool printed_something{false};
-           const auto &[vidx, d] : derivatives[i])
-        {
-          if (exchange(printed_something, true))
-            d_output[i] << ", ";
-
-          int eq = vidx[0];
-
-          int col_idx = 0;
-          for (size_t j = 1; j < vidx.size(); j++)
-            {
-              col_idx *= symbol_table.endo_nbr();
-              col_idx += getJacobianCol(vidx[j]);
-            }
-
-          if (writeDetails)
-            d_output[i] << R"({"eq": )" << eq + 1;
-          else
-            d_output[i] << R"({"row": )" << eq + 1;
-
-          d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
-
-          if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
-            {
-              int col_idx_sym = getJacobianCol(vidx[2]) * symbol_table.endo_nbr() + getJacobianCol(vidx[1]);
-              d_output[i] << ", " << col_idx_sym + 1;
-            }
-          if (i > 1)
-            d_output[i] << "]";
-
-          if (writeDetails)
-            for (size_t j = 1; j < vidx.size(); j++)
-              d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << getNameByDerivID(vidx[j]) << R"(")";
-
-          d_output[i] << R"(, "val": ")";
-          d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
-          d_output[i] << R"("})" << endl;
-        }
-      d_output[i] << "]}";
-
-      ncols *= symbol_table.endo_nbr();
-    }
+  auto [mlv_output, d_output] { writeJsonComputingPassOutputHelper<false>(writeDetails) };
 
   if (writeDetails)
     output << R"("static_model": {)";
   else
     output << R"("static_model_simple": {)";
-  output << model_local_vars_output.str();
+  output << mlv_output.str();
   for (const auto &it : d_output)
     output << ", " << it.str();
   output << "}";
 }
 
 void
-StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
+StaticModel::writeJsonParamsDerivatives(ostream &output, bool writeDetails) const
 {
   if (!params_derivatives.size())
     return;
 
-  ostringstream model_local_vars_output; // Used for storing model local vars
-  ostringstream model_output; // Used for storing model temp vars and equations
-  ostringstream jacobian_output; // Used for storing jacobian equations
-  ostringstream hessian_output; // Used for storing Hessian equations
-  ostringstream hessian1_output; // Used for storing Hessian equations
-  ostringstream third_derivs_output; // Used for storing third order derivatives equations
-  ostringstream third_derivs1_output; // Used for storing third order derivatives equations
-
-  deriv_node_temp_terms_t tef_terms;
-  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
-
-  temporary_terms_t temp_term_union;
-  for (const auto &it : params_derivs_temporary_terms)
-    writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, "all");
-
-  jacobian_output << R"("deriv_wrt_params": {)"
-                  << R"(  "neqs": )" << equations.size()
-                  << R"(, "nparamcols": )" << symbol_table.param_nbr()
-                  << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 0, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        jacobian_output << ", ";
-
-      auto [eq, param] = vectorToTuple<2>(vidx);
-
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        jacobian_output << R"({"eq": )" << eq + 1;
-      else
-        jacobian_output << R"({"row": )" << eq + 1;
-
-      if (writeDetails)
-        jacobian_output << R"(, "param_col": )" << param_col;
-
-      jacobian_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      jacobian_output << R"(, "val": ")";
-      d->writeJsonOutput(jacobian_output, temp_term_union, tef_terms);
-      jacobian_output << R"("})" << endl;
-    }
-  jacobian_output << "]}";
-
-  hessian_output << R"("deriv_jacobian_wrt_params": {)"
-                 << R"(  "neqs": )" << equations.size()
-                 << R"(, "nvarcols": )" << symbol_table.endo_nbr()
-                 << R"(, "nparamcols": )" << symbol_table.param_nbr()
-                 << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 1, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        hessian_output << ", ";
-
-      auto [eq, var, param] = vectorToTuple<3>(vidx);
-
-      int var_col { getTypeSpecificIDByDerivID(var) + 1 };
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        hessian_output << R"({"eq": )" << eq + 1;
-      else
-        hessian_output << R"({"row": )" << eq + 1;
-
-      if (writeDetails)
-        hessian_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
-                       << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
-
-      hessian_output << R"(, "var_col": )" << var_col
-                     << R"(, "param_col": )" << param_col
-                     << R"(, "val": ")";
-      d->writeJsonOutput(hessian_output, temp_term_union, tef_terms);
-      hessian_output << R"("})" << endl;
-    }
-  hessian_output << "]}";
-
-  hessian1_output << R"("second_deriv_residuals_wrt_params": {)"
-                  << R"(  "nrows": )" << equations.size()
-                  << R"(, "nparam1cols": )" << symbol_table.param_nbr()
-                  << R"(, "nparam2cols": )" << symbol_table.param_nbr()
-                  << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 0, 2 })->second)
-    {
-      if (exchange(printed_something, true))
-        hessian1_output << ", ";
-
-      auto [eq, param1, param2] = vectorToTuple<3>(vidx);
-
-      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
-      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
-
-      if (writeDetails)
-        hessian1_output << R"({"eq": )" << eq + 1;
-      else
-        hessian1_output << R"({"row": )" << eq + 1;
-
-      hessian1_output << R"(, "param1_col": )" << param1_col
-                      << R"(, "param2_col": )" << param2_col;
-
-      if (writeDetails)
-        hessian1_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
-                        << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
-
-      hessian1_output << R"(, "val": ")";
-      d->writeJsonOutput(hessian1_output, temp_term_union, tef_terms);
-      hessian1_output << R"("})" << endl;
-    }
-  hessian1_output << "]}";
-
-  third_derivs_output << R"("second_deriv_jacobian_wrt_params": {)"
-                      << R"(  "neqs": )" << equations.size()
-                      << R"(, "nvarcols": )" << symbol_table.endo_nbr()
-                      << R"(, "nparam1cols": )" << symbol_table.param_nbr()
-                      << R"(, "nparam2cols": )" << symbol_table.param_nbr()
-                      << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 1, 2 })->second)
-    {
-      if (exchange(printed_something, true))
-        third_derivs_output << ", ";
-
-      auto [eq, var, param1, param2] = vectorToTuple<4>(vidx);
-
-      int var_col { getTypeSpecificIDByDerivID(var) + 1 };
-      int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
-      int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
-
-      if (writeDetails)
-        third_derivs_output << R"({"eq": )" << eq + 1;
-      else
-        third_derivs_output << R"({"row": )" << eq + 1;
-      third_derivs_output << R"(, "var_col": )" << var_col
-                          << R"(, "param1_col": )" << param1_col
-                          << R"(, "param2_col": )" << param2_col;
-
-      if (writeDetails)
-        third_derivs_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
-                            << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
-                            << R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
-
-      third_derivs_output << R"(, "val": ")";
-      d->writeJsonOutput(third_derivs_output, temp_term_union, tef_terms);
-      third_derivs_output << R"("})" << endl;
-    }
-  third_derivs_output << "]}" << endl;
-
-  third_derivs1_output << R"("derivative_hessian_wrt_params": {)"
-                       << R"(  "neqs": )" << equations.size()
-                       << R"(, "nvar1cols": )" << symbol_table.endo_nbr()
-                       << R"(, "nvar2cols": )" << symbol_table.endo_nbr()
-                       << R"(, "nparamcols": )" << symbol_table.param_nbr()
-                       << R"(, "entries": [)";
-  for (bool printed_something{false};
-       const auto &[vidx, d] : params_derivatives.find({ 2, 1 })->second)
-    {
-      if (exchange(printed_something, true))
-        third_derivs1_output << ", ";
-
-      auto [eq, var1, var2, param] = vectorToTuple<4>(vidx);
-
-      int var1_col { getTypeSpecificIDByDerivID(var1) + 1 };
-      int var2_col { getTypeSpecificIDByDerivID(var2) + 1 };
-      int param_col { getTypeSpecificIDByDerivID(param) + 1 };
-
-      if (writeDetails)
-        third_derivs1_output << R"({"eq": )" << eq + 1;
-      else
-        third_derivs1_output << R"({"row": )" << eq + 1;
-
-      third_derivs1_output << R"(, "var1_col": )" << var1_col
-                           << R"(, "var2_col": )" << var2_col
-                           << R"(, "param_col": )" << param_col;
-
-      if (writeDetails)
-        third_derivs1_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
-                             << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")"
-                             << R"(, "param1": ")" << getNameByDerivID(param) << R"(")";
-
-      third_derivs1_output << R"(, "val": ")";
-      d->writeJsonOutput(third_derivs1_output, temp_term_union, tef_terms);
-      third_derivs1_output << R"("})" << endl;
-    }
-  third_derivs1_output << "]}" << endl;
+  auto [mlv_output, tt_output, rp_output, gp_output, rpp_output, gpp_output, hp_output, g3p_output]
+    { writeJsonParamsDerivativesHelper<false>(writeDetails) };
+  // g3p_output is ignored
 
   if (writeDetails)
     output << R"("static_model_params_derivative": {)";
   else
     output << R"("static_model_params_derivatives_simple": {)";
-  output << model_local_vars_output.str()
-         << ", " << model_output.str()
-         << ", " << jacobian_output.str()
-         << ", " << hessian_output.str()
-         << ", " << hessian1_output.str()
-         << ", " << third_derivs_output.str()
-         << ", " << third_derivs1_output.str()
+  output << mlv_output.str()
+         << ", " << tt_output.str()
+         << ", " << rp_output.str()
+         << ", " << gp_output.str()
+         << ", " << rpp_output.str()
+         << ", " << gpp_output.str()
+         << ", " << hp_output.str()
          << "}";
 }
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index dcd0118376c103488539fda6f1325f18abd11338..8948d846d0d2756528901479a7bd410ff8900d24 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -145,8 +145,8 @@ public:
   //! Write JSON representation of static model
   void writeJsonComputingPassOutput(ostream &output, bool writeDetails) const;
 
-  //! Writes file containing static parameters derivatives
-  void writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const;
+  //! Write JSON params derivatives
+  void writeJsonParamsDerivatives(ostream &output, bool writeDetails) const;
 
   //! Writes file containing static parameters derivatives
   template<bool julia>