diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index a3dd2d721d2929004c179bfac0dabb2dd7e9504d..e22bd2cc4c106e8f84fd4302b77e914cfa11388f 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4756,344 +4756,6 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context
         }
 }
 
-void
-DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) const
-{
-  if (!params_derivatives.size())
-    return;
-
-  ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaDynamicModel : ExprNodeOutputType::matlabDynamicModel);
-  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
-
-  temporary_terms_t temp_term_union;
-  deriv_node_temp_terms_t tef_terms;
-
-  writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
-  for (const auto &it : params_derivs_temporary_terms)
-    writeTemporaryTerms(it.second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
-
-  for (const auto & [indices, d1] : params_derivatives.find({ 0, 1 })->second)
-    {
-      auto [eq, param] = vectorToTuple<2>(indices);
-
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      rp_output << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << param_col
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
-      d1->writeOutput(rp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      rp_output << ";" << endl;
-    }
-
-  for (const auto & [indices, d2] : params_derivatives.find({ 1, 1 })->second)
-    {
-      auto [eq, var, param] = vectorToTuple<3>(indices);
-
-      int var_col = getJacobianCol(var) + 1;
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col
-                << ", " << param_col << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
-      d2->writeOutput(gp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      gp_output << ";" << endl;
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 0, 2 })->second)
-    {
-      auto [eq, param1, param2] = vectorToTuple<3>(indices);
-
-      int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
-      int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
-
-      rpp_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) << "=" << param1_col << ";" << endl
-                 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
-                 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(rpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      rpp_output << ";" << endl;
-
-      i++;
-
-      if (param1 != param2)
-        {
-          // Treat symmetric elements
-          rpp_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++;
-        }
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 1, 2 })->second)
-    {
-      auto [eq, var, param1, param2] = vectorToTuple<4>(indices);
-
-      int var_col = getJacobianCol(var) + 1;
-      int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
-      int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
-
-      gpp_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) << "=" << param1_col << ";" << endl
-                 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
-                 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(gpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      gpp_output << ";" << endl;
-
-      i++;
-
-      if (param1 != param2)
-        {
-          // Treat symmetric elements
-          gpp_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++;
-        }
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 2, 1 })->second)
-    {
-      auto [eq, var1, var2, param] = vectorToTuple<4>(indices);
-
-      int var1_col = getJacobianCol(var1) + 1;
-      int var2_col = getJacobianCol(var2) + 1;
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
-                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
-                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
-                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
-                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(hp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      hp_output << ";" << endl;
-
-      i++;
-
-      if (var1 != var2)
-        {
-          // Treat symmetric elements
-          hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
-                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
-                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
-                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
-                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                    << "=hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
-                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
-          i++;
-        }
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 3, 1 })->second)
-    {
-      auto [eq, var1, var2, var3, param] = vectorToTuple<5>(indices);
-
-      int var1_col = getJacobianCol(var1) + 1;
-      int var2_col = getJacobianCol(var2) + 1;
-      int var3_col = getJacobianCol(var3) + 1;
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
-                 << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
-                 << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
-                 << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var3_col << ";" << endl
-                 << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
-                 << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",6"
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(g3p_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      g3p_output << ";" << endl;
-
-      i++;
-    }
-
-  string filename = julia ? basename + "DynamicParamsDerivs.jl" : packageDir(basename) + "/dynamic_params_derivs.m";
-  ofstream paramsDerivsFile{filename, ios::out | ios::binary};
-  if (!paramsDerivsFile.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  if (!julia)
-    {
-      // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
-      map<string, string> tmp_paren_vars;
-      bool message_printed = false;
-      fixNestedParenthesis(tt_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(rp_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(gp_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(rpp_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(gpp_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(hp_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(g3p_output, tmp_paren_vars, message_printed);
-      paramsDerivsFile << "function [rp, gp, rpp, gpp, hp, g3p] = dynamic_params_derivs(y, x, params, steady_state, it_, ss_param_deriv, ss_param_2nd_deriv)" << endl
-                       << "%" << endl
-                       << "% Compute the derivatives of the dynamic model with respect to the parameters" << endl
-                       << "% Inputs :" << endl
-                       << "%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored" << endl
-                       << "%                                                 in M_.lead_lag_incidence; see the Manual" << endl
-                       << "%   x         [nperiods by M_.exo_nbr] double     matrix of exogenous variables (in declaration order)" << endl
-                       << "%                                                 for all simulation periods" << endl
-                       << "%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order" << endl
-                       << "%   steady_state  [M_.endo_nbr by 1] double       vector of steady state values" << endl
-                       << "%   it_       scalar double                       time period for exogenous variables for which to evaluate the model" << endl
-                       << "%   ss_param_deriv     [M_.eq_nbr by #params]     Jacobian matrix of the steady states values with respect to the parameters" << endl
-                       << "%   ss_param_2nd_deriv [M_.eq_nbr by #params by #params] Hessian matrix of the steady states values with respect to the parameters" << endl
-                       << "%" << endl
-                       << "% Outputs:" << endl
-                       << "%   rp        [M_.eq_nbr by #params] double    Jacobian matrix of dynamic model equations with respect to parameters " << endl
-                       << "%                                              Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
-                       << "%   gp        [M_.endo_nbr by #dynamic variables by #params] double    Derivative of the Jacobian matrix of the dynamic model equations with respect to the parameters" << endl
-                       << "%                                                           rows: equations in order of declaration" << endl
-                       << "%                                                           columns: variables in order stored in M_.lead_lag_incidence" << endl
-                       << "%   rpp       [#second_order_residual_terms by 4] double   Hessian matrix of second derivatives of residuals with respect to parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: number of the first parameter in derivative" << endl
-                       << "%                                                              3rd column: number of the second parameter in derivative" << endl
-                       << "%                                                              4th column: value of the Hessian term" << endl
-                       << "%   gpp      [#second_order_Jacobian_terms by 5] double   Hessian matrix of second derivatives of the Jacobian with respect to the parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: column number of variable in Jacobian of the dynamic model" << endl
-                       << "%                                                              3rd column: number of the first parameter in derivative" << endl
-                       << "%                                                              4th column: number of the second parameter in derivative" << endl
-                       << "%                                                              5th column: value of the Hessian term" << endl
-                       << "%   hp      [#first_order_Hessian_terms by 5] double   Jacobian matrix of derivatives of the dynamic Hessian with respect to the parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: column number of first variable in Hessian of the dynamic model" << endl
-                       << "%                                                              3rd column: column number of second variable in Hessian of the dynamic model" << endl
-                       << "%                                                              4th column: number of the parameter in derivative" << endl
-                       << "%                                                              5th column: value of the Hessian term" << endl
-                       << "%   g3p      [#first_order_g3_terms by 6] double   Jacobian matrix of derivatives of g3 (dynamic 3rd derivs) with respect to the parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: column number of first variable in g3 of the dynamic model" << endl
-                       << "%                                                              3rd column: column number of second variable in g3 of the dynamic model" << endl
-                       << "%                                                              4th column: column number of third variable in g3 of the dynamic model" << endl
-                       << "%                                                              5th column: number of the parameter in derivative" << endl
-                       << "%                                                              6th column: value of the Hessian term" << endl
-                       << "%" << endl
-                       << "%" << endl
-                       << "% Warning : this file is generated automatically by Dynare" << endl
-                       << "%           from model file (.mod)" << endl << endl
-                       << "T = NaN(" << params_derivs_temporary_terms_idxs.size() << ",1);" << endl
-                       << tt_output.str()
-                       << "rp = zeros(" << equations.size() << ", "
-                       << symbol_table.param_nbr() << ");" << endl
-                       << rp_output.str()
-                       << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl
-                       << gp_output.str()
-                       << "if nargout >= 3" << endl
-                       << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
-                       << rpp_output.str()
-                       << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
-                       << gpp_output.str()
-                       << "end" << endl
-                       << "if nargout >= 5" << endl
-                       << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
-                       << hp_output.str()
-                       << "end" << endl
-                       << "if nargout >= 6" << endl
-                       << "g3p = zeros(" << params_derivatives.find({ 3, 1 })->second.size() << ",6);" << endl
-                       << g3p_output.str()
-                       << "end" << endl
-                       << "end" << endl;
-    }
-  else
-    paramsDerivsFile << "module " << basename << "DynamicParamsDerivs" << endl
-                     << "#" << endl
-                     << "# NB: this file was automatically generated by Dynare" << endl
-                     << "#     from " << basename << ".mod" << endl
-                     << "#" << endl
-                     << "export params_derivs" << endl << endl
-                     << "function params_derivs(y, x, paramssteady_state, it_, "
-                     << "ss_param_deriv, ss_param_2nd_deriv)" << endl
-		     << "@inbounds begin" << endl
-                     << tt_output.str()
-		     << "end" << endl
-                     << "rp = zeros(" << equations.size() << ", "
-                     << symbol_table.param_nbr() << ");" << endl
-		     << "@inbounds begin" << endl
-                     << rp_output.str()
-		     << "end" << endl
-                     << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl
-		     << "@inbounds begin" << endl
-                     << gp_output.str()
-		     << "end" << endl
-                     << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
-		     << "@inbounds begin" << endl
-                     << rpp_output.str()
-		     << "end" << endl
-                     << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
-		     << "@inbounds begin" << endl
-                     << gpp_output.str()
-		     << "end" << endl
-                     << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
-		     << "@inbounds begin" << endl
-                     << hp_output.str()
-		     << "end" << endl
-                     << "g3p = zeros(" << params_derivatives.find({ 3, 1 })->second.size() << ",6);" << endl
-		     << "@inbounds begin" << endl
-                     << g3p_output.str()
-		     << "end" << endl
-                     << "(rp, gp, rpp, gpp, hp, g3p)" << endl
-                     << "end" << endl
-                     << "end" << endl;
-
-  paramsDerivsFile.close();
-}
-
 void
 DynamicModel::writeLatexFile(const string &basename, bool write_equation_tags) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 5cdaf360be07a63476bebed973a64e002cbb3e7c..15adc46bc8517d6e138406608bb299be5a7eee53 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -392,8 +392,10 @@ public:
 
   //! Writes dynamic model file (+ bytecode)
   void writeDynamicFile(const string &basename, bool block, 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;
+  template<bool julia>
+  void writeParamsDerivativesFile(const string &basename) const;
 
   //! Writes file containing coordinates of non-zero elements in the Jacobian
   /*! Used by the perfect_foresight_problem MEX */
@@ -656,4 +658,146 @@ public:
   // Returns the set of equations (as numbers) which have a pac_expectation operator
   set<int> findPacExpectationEquationNumbers() const;
 };
+
+template<bool julia>
+void
+DynamicModel::writeParamsDerivativesFile(const string &basename) const
+{
+  if (!params_derivatives.size())
+    return;
+
+  constexpr ExprNodeOutputType output_type { julia ? ExprNodeOutputType::juliaDynamicModel : ExprNodeOutputType::matlabDynamicModel };
+
+  auto [tt_output, rp_output, gp_output, rpp_output, gpp_output, hp_output, g3p_output]
+    { writeParamsDerivativesFileHelper<output_type>() };
+
+  string filename { julia ? basename + "DynamicParamsDerivs.jl" : packageDir(basename) + "/dynamic_params_derivs.m" };
+  ofstream paramsDerivsFile { filename, ios::out | ios::binary };
+  if (!paramsDerivsFile.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if constexpr(!julia)
+    {
+      paramsDerivsFile << "function [rp, gp, rpp, gpp, hp, g3p] = dynamic_params_derivs(y, x, params, steady_state, it_, ss_param_deriv, ss_param_2nd_deriv)" << endl
+                       << "%" << endl
+                       << "% Compute the derivatives of the dynamic model with respect to the parameters" << endl
+                       << "% Inputs :" << endl
+                       << "%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored" << endl
+                       << "%                                                 in M_.lead_lag_incidence; see the Manual" << endl
+                       << "%   x         [nperiods by M_.exo_nbr] double     matrix of exogenous variables (in declaration order)" << endl
+                       << "%                                                 for all simulation periods" << endl
+                       << "%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order" << endl
+                       << "%   steady_state  [M_.endo_nbr by 1] double       vector of steady state values" << endl
+                       << "%   it_       scalar double                       time period for exogenous variables for which to evaluate the model" << endl
+                       << "%   ss_param_deriv     [M_.eq_nbr by #params]     Jacobian matrix of the steady states values with respect to the parameters" << endl
+                       << "%   ss_param_2nd_deriv [M_.eq_nbr by #params by #params] Hessian matrix of the steady states values with respect to the parameters" << endl
+                       << "%" << endl
+                       << "% Outputs:" << endl
+                       << "%   rp        [M_.eq_nbr by #params] double    Jacobian matrix of dynamic model equations with respect to parameters " << endl
+                       << "%                                              Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
+                       << "%   gp        [M_.endo_nbr by #dynamic variables by #params] double    Derivative of the Jacobian matrix of the dynamic model equations with respect to the parameters" << endl
+                       << "%                                                           rows: equations in order of declaration" << endl
+                       << "%                                                           columns: variables in order stored in M_.lead_lag_incidence" << endl
+                       << "%   rpp       [#second_order_residual_terms by 4] double   Hessian matrix of second derivatives of residuals with respect to parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: number of the first parameter in derivative" << endl
+                       << "%                                                              3rd column: number of the second parameter in derivative" << endl
+                       << "%                                                              4th column: value of the Hessian term" << endl
+                       << "%   gpp      [#second_order_Jacobian_terms by 5] double   Hessian matrix of second derivatives of the Jacobian with respect to the parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: column number of variable in Jacobian of the dynamic model" << endl
+                       << "%                                                              3rd column: number of the first parameter in derivative" << endl
+                       << "%                                                              4th column: number of the second parameter in derivative" << endl
+                       << "%                                                              5th column: value of the Hessian term" << endl
+                       << "%   hp      [#first_order_Hessian_terms by 5] double   Jacobian matrix of derivatives of the dynamic Hessian with respect to the parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: column number of first variable in Hessian of the dynamic model" << endl
+                       << "%                                                              3rd column: column number of second variable in Hessian of the dynamic model" << endl
+                       << "%                                                              4th column: number of the parameter in derivative" << endl
+                       << "%                                                              5th column: value of the Hessian term" << endl
+                       << "%   g3p      [#first_order_g3_terms by 6] double   Jacobian matrix of derivatives of g3 (dynamic 3rd derivs) with respect to the parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: column number of first variable in g3 of the dynamic model" << endl
+                       << "%                                                              3rd column: column number of second variable in g3 of the dynamic model" << endl
+                       << "%                                                              4th column: column number of third variable in g3 of the dynamic model" << endl
+                       << "%                                                              5th column: number of the parameter in derivative" << endl
+                       << "%                                                              6th column: value of the Hessian term" << endl
+                       << "%" << endl
+                       << "%" << endl
+                       << "% Warning : this file is generated automatically by Dynare" << endl
+                       << "%           from model file (.mod)" << endl << endl
+                       << "T = NaN(" << params_derivs_temporary_terms_idxs.size() << ",1);" << endl
+                       << tt_output.str()
+                       << "rp = zeros(" << equations.size() << ", "
+                       << symbol_table.param_nbr() << ");" << endl
+                       << rp_output.str()
+                       << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl
+                       << gp_output.str()
+                       << "if nargout >= 3" << endl
+                       << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+                       << rpp_output.str()
+                       << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+                       << gpp_output.str()
+                       << "end" << endl
+                       << "if nargout >= 5" << endl
+                       << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+                       << hp_output.str()
+                       << "end" << endl
+                       << "if nargout >= 6" << endl
+                       << "g3p = zeros(" << params_derivatives.find({ 3, 1 })->second.size() << ",6);" << endl
+                       << g3p_output.str()
+                       << "end" << endl
+                       << "end" << endl;
+    }
+  else
+    paramsDerivsFile << "module " << basename << "DynamicParamsDerivs" << endl
+                     << "#" << endl
+                     << "# NB: this file was automatically generated by Dynare" << endl
+                     << "#     from " << basename << ".mod" << endl
+                     << "#" << endl
+                     << "export params_derivs" << endl << endl
+                     << "function params_derivs(y, x, paramssteady_state, it_, "
+                     << "ss_param_deriv, ss_param_2nd_deriv)" << endl
+		     << "@inbounds begin" << endl
+                     << tt_output.str()
+		     << "end" << endl
+                     << "rp = zeros(" << equations.size() << ", "
+                     << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
+                     << rp_output.str()
+		     << "end" << endl
+                     << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
+                     << gp_output.str()
+		     << "end" << endl
+                     << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+		     << "@inbounds begin" << endl
+                     << rpp_output.str()
+		     << "end" << endl
+                     << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
+                     << gpp_output.str()
+		     << "end" << endl
+                     << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
+                     << hp_output.str()
+		     << "end" << endl
+                     << "g3p = zeros(" << params_derivatives.find({ 3, 1 })->second.size() << ",6);" << endl
+		     << "@inbounds begin" << endl
+                     << g3p_output.str()
+		     << "end" << endl
+                     << "(rp, gp, rpp, gpp, hp, g3p)" << endl
+                     << "end" << endl
+                     << "end" << endl;
+
+  paramsDerivsFile.close();
+}
+
 #endif
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 58dfdbcdccae39f733d8045d2e224536265ef8ac..9bbc2ad044580582f7f9aad636149941b3750650 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1072,12 +1072,12 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
           if (!no_static)
             {
               static_model.writeStaticFile(basename, block, use_dll, mexext, matlabroot, dynareroot, false);
-              static_model.writeParamsDerivativesFile(basename, false);
+              static_model.writeParamsDerivativesFile<false>(basename);
             }
 
           dynamic_model.writeDynamicFile(basename, block, use_dll, mexext, matlabroot, dynareroot, false);
 
-          dynamic_model.writeParamsDerivativesFile(basename, false);
+          dynamic_model.writeParamsDerivativesFile<false>(basename);
 
           dynamic_model.writeDynamicJacobianNonZeroElts(basename);
         }
@@ -1102,10 +1102,10 @@ ModFile::writeJuliaOutput(const string &basename) const
       if (!no_static)
         {
           static_model.writeStaticFile(basename, false, false, "", {}, {}, true);
-          static_model.writeParamsDerivativesFile(basename, true);
+          static_model.writeParamsDerivativesFile<true>(basename);
         }
       dynamic_model.writeDynamicFile(basename, block, use_dll, "", {}, {}, true);
-      dynamic_model.writeParamsDerivativesFile(basename, true);
+      dynamic_model.writeParamsDerivativesFile<true>(basename);
     }
   steady_state_model.writeSteadyStateFile(basename, true);
 }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 8ff7d07c00685b6c2de2c1f27bc9ed31e49d6123..90513662550a2e0d2e6ef5fc9d3e24f5e6acebd9 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -262,6 +262,13 @@ protected:
   template<ExprNodeOutputType output_type>
   pair<vector<ostringstream>, vector<ostringstream>> writeModelFileHelper() const;
 
+  /* Helper for writing derivatives w.r.t. parameters.
+     Returns { tt, rp, gp, rpp, gpp, hp, g3p }.
+     g3p is empty if requesting a static output type. */
+  template<ExprNodeOutputType output_type>
+  tuple<ostringstream, ostringstream, ostringstream, ostringstream,
+        ostringstream, ostringstream, ostringstream> writeParamsDerivativesFileHelper() 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)
@@ -657,4 +664,223 @@ ModelTree::writeModelFileHelper() const
   return { move(d_output), move(tt_output) };
 }
 
+template<ExprNodeOutputType output_type>
+tuple<ostringstream, ostringstream, ostringstream, ostringstream,
+      ostringstream, ostringstream, ostringstream>
+ModelTree::writeParamsDerivativesFileHelper() const
+{
+  static_assert(!isCOutput(output_type), "C output is not implemented");
+
+  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 (only in dynamic case)
+
+  temporary_terms_t temp_term_union;
+  deriv_node_temp_terms_t tef_terms;
+
+  writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
+  for (const auto &[order, tts] : params_derivs_temporary_terms)
+    writeTemporaryTerms(tts, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
+
+  for (const auto &[indices, d1] : params_derivatives.find({ 0, 1 })->second)
+    {
+      auto [eq, param] { vectorToTuple<2>(indices) };
+
+      int param_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1 };
+
+      rp_output << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << param_col
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d1->writeOutput(rp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+      rp_output << ";" << endl;
+    }
+
+  for (const auto &[indices, d2] : params_derivatives.find({ 1, 1 })->second)
+    {
+      auto [eq, var, param] { vectorToTuple<3>(indices) };
+
+      int var_col { getJacobianCol(var) + 1 };
+      int param_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1 };
+
+      gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col
+                << ", " << param_col << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d2->writeOutput(gp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+      gp_output << ";" << endl;
+    }
+
+  for (int i {1};
+       const auto &[indices, d2] : params_derivatives.find({ 0, 2 })->second)
+    {
+      auto [eq, param1, param2] { vectorToTuple<3>(indices) };
+
+      int param1_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1 };
+      int param2_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1 };
+
+      rpp_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) << "=" << param1_col << ";" << endl
+                 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(rpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+      rpp_output << ";" << endl;
+
+      i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          rpp_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++;
+        }
+    }
+
+  for (int i {1};
+       const auto &[indices, d2] : params_derivatives.find({ 1, 2 })->second)
+    {
+      auto [eq, var, param1, param2] { vectorToTuple<4>(indices) };
+
+      int var_col { getJacobianCol(var) + 1 };
+      int param1_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1 };
+      int param2_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1 };
+
+      gpp_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) << "=" << param1_col << ";" << endl
+                 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(gpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+      gpp_output << ";" << endl;
+
+      i++;
+
+      if (param1 != param2)
+        {
+          // Treat symmetric elements
+          gpp_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++;
+        }
+    }
+
+  for (int i {1};
+       const auto &[indices, d2] : params_derivatives.find({ 2, 1 })->second)
+    {
+      auto [eq, var1, var2, param] { vectorToTuple<4>(indices) };
+
+      int var1_col { getJacobianCol(var1) + 1 };
+      int var2_col { getJacobianCol(var2) + 1 };
+      int param_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1 };
+
+      hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
+                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
+                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
+                << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(hp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+      hp_output << ";" << endl;
+
+      i++;
+
+      if (var1 != var2)
+        {
+          // Treat symmetric elements
+          hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
+                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
+                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
+                    << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                    << "=hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
+                    << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
+          i++;
+        }
+    }
+
+  if constexpr(output_type == ExprNodeOutputType::matlabDynamicModel
+               || output_type == ExprNodeOutputType::juliaDynamicModel)
+    for (int i {1};
+         const auto &[indices, d2] : params_derivatives.find({ 3, 1 })->second)
+      {
+        auto [eq, var1, var2, var3, param] { vectorToTuple<5>(indices) };
+
+        int var1_col { getJacobianCol(var1) + 1 };
+        int var2_col { getJacobianCol(var2) + 1 };
+        int var3_col { getJacobianCol(var3) + 1 };
+        int param_col { symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1 };
+
+        g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                   << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
+                   << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
+                   << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var3_col << ";" << endl
+                   << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
+                   << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",6"
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+        d2->writeOutput(g3p_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
+        g3p_output << ";" << endl;
+
+        i++;
+      }
+
+  if constexpr(isMatlabOutput(output_type))
+    {
+      // Check that we don't have more than 32 nested parenthesis because MATLAB does not support this. See Issue #1201
+      map<string, string> tmp_paren_vars;
+      bool message_printed {false};
+      fixNestedParenthesis(tt_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(rp_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(gp_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(rpp_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(gpp_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(hp_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(g3p_output, tmp_paren_vars, message_printed);
+    }
+
+  return { 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 9ed6106349b9c8b3bc098a56fc411ad8a887aa5f..0a276ad200a11cdf3cd254c13576b849d6f2e1fd 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1928,291 +1928,6 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
     }
 }
 
-void
-StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) const
-{
-  if (!params_derivatives.size())
-    return;
-
-  ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel);
-
-  ostringstream tt_output; // Used for storing temporary terms
-  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
-
-  temporary_terms_t temp_term_union;
-  deriv_node_temp_terms_t tef_terms;
-
-  writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
-  for (const auto &it : params_derivs_temporary_terms)
-    writeTemporaryTerms(it.second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
-
-  for (const auto & [indices, d1] : params_derivatives.find({ 0, 1 })->second)
-    {
-      auto [eq, param] = vectorToTuple<2>(indices);
-
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      jacobian_output << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                      <<  eq+1 << ", " << param_col
-                      << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
-      d1->writeOutput(jacobian_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      jacobian_output << ";" << endl;
-    }
-
-  for (const auto & [indices, d2] : params_derivatives.find({ 1, 1 })->second)
-    {
-      auto [eq, var, param] = vectorToTuple<3>(indices);
-
-      int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      hessian_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                     << eq+1 << ", " << var_col << ", " << param_col
-                     << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
-      d2->writeOutput(hessian_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      hessian_output << ";" << endl;
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 0, 2 })->second)
-    {
-      auto [eq, param1, param2] = vectorToTuple<3>(indices);
-
-      int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
-      int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
-
-      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) << "=" << param1_col << ";" << endl
-                      << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                      << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
-                      << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                      << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(hessian1_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      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++;
-        }
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 1, 2 })->second)
-    {
-      auto [eq, var, param1, param2] = vectorToTuple<4>(indices);
-
-      int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
-      int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
-      int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
-
-      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) << "=" << param1_col << ";" << endl
-                          << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
-                          << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                          << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(third_derivs_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      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++;
-        }
-    }
-
-  for (int i{1};
-       const auto &[indices, d2] : params_derivatives.find({ 2, 1 })->second)
-    {
-      auto [eq, var1, var2, param] = vectorToTuple<4>(indices);
-
-      int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1;
-      int var2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var2)) + 1;
-      int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
-
-      third_derivs1_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
-                           << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
-                           << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
-                           << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
-                           << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                           << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
-                           << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                           << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
-                           << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                           << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
-      d2->writeOutput(third_derivs1_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
-      third_derivs1_output << ";" << endl;
-
-      i++;
-
-      if (var1 != var2)
-        {
-          // Treat symmetric elements
-          third_derivs1_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
-                               << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
-                               << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
-                               << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
-                               << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                               << "=hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
-                               << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
-          i++;
-        }
-    }
-
-  string filename = julia ? basename + "StaticParamsDerivs.jl" : packageDir(basename) + "/static_params_derivs.m";
-  ofstream paramsDerivsFile{filename, ios::out | ios::binary};
-  if (!paramsDerivsFile.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  if (!julia)
-    {
-      // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
-      map<string, string> tmp_paren_vars;
-      bool message_printed = false;
-      fixNestedParenthesis(tt_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(jacobian_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(hessian_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(hessian1_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(third_derivs_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(third_derivs1_output, tmp_paren_vars, message_printed);
-
-      paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = static_params_derivs(y, x, params)" << endl
-                       << "%" << endl
-                       << "% Status : Computes derivatives of the static model with respect to the parameters" << endl
-                       << "%" << endl
-                       << "% Inputs : " << endl
-                       << "%   y         [M_.endo_nbr by 1] double    vector of endogenous variables in declaration order" << endl
-                       << "%   x         [M_.exo_nbr by 1] double     vector of exogenous variables in declaration order" << endl
-                       << "%   params    [M_.param_nbr by 1] double   vector of parameter values in declaration order" << endl
-                       << "%" << endl
-                       << "% Outputs:" << endl
-                       << "%   rp        [M_.eq_nbr by #params] double    Jacobian matrix of static model equations with respect to parameters " << endl
-                       << "%                                              Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
-                       << "%   gp        [M_.endo_nbr by M_.endo_nbr by #params] double    Derivative of the Jacobian matrix of the static model equations with respect to the parameters" << endl
-                       << "%                                                           rows: variables in declaration order" << endl
-                       << "%                                                           rows: equations in order of declaration" << endl
-                       << "%   rpp       [#second_order_residual_terms by 4] double   Hessian matrix of second derivatives of residuals with respect to parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: number of the first parameter in derivative" << endl
-                       << "%                                                              3rd column: number of the second parameter in derivative" << endl
-                       << "%                                                              4th column: value of the Hessian term" << endl
-                       << "%   gpp      [#second_order_Jacobian_terms by 5] double   Hessian matrix of second derivatives of the Jacobian with respect to the parameters;" << endl
-                       << "%                                                              rows: respective derivative term" << endl
-                       << "%                                                              1st column: equation number of the term appearing" << endl
-                       << "%                                                              2nd column: column number of variable in Jacobian of the static model" << endl
-                       << "%                                                              3rd column: number of the first parameter in derivative" << endl
-                       << "%                                                              4th column: number of the second parameter in derivative" << endl
-                       << "%                                                              5th column: value of the Hessian term" << endl
-                       << "%" << endl
-                       << "%" << endl
-                       << "% Warning : this file is generated automatically by Dynare" << endl
-                       << "%           from model file (.mod)" << endl << endl
-                       << "T = NaN(" << params_derivs_temporary_terms_idxs.size() << ",1);" << endl
-                       << tt_output.str()
-                       << "rp = zeros(" << equations.size() << ", "
-                       << symbol_table.param_nbr() << ");" << endl
-                       << jacobian_output.str()
-                       << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
-                       << symbol_table.param_nbr() << ");" << endl
-                       << hessian_output.str()
-                       << "if nargout >= 3" << endl
-                       << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
-                       << hessian1_output.str()
-                       << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
-                       << third_derivs_output.str()
-                       << "end" << endl
-                       << "if nargout >= 5" << endl
-                       << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
-                       << third_derivs1_output.str()
-                       << "end" << endl
-                       << "end" << endl;
-    }
-  else
-    paramsDerivsFile << "module " << basename << "StaticParamsDerivs" << endl
-                     << "#" << endl
-                     << "# NB: this file was automatically generated by Dynare" << endl
-                     << "#     from " << basename << ".mod" << endl
-                     << "#" << endl
-                     << "export params_derivs" << endl << endl
-                     << "function params_derivs(y, x, params)" << endl
-		     << "@inbounds begin" << endl
-                     << tt_output.str()
-		     << "end" << endl
-		     << "rp = zeros(" << equations.size() << ", "
-                     << symbol_table.param_nbr() << ");" << endl
-		     << "@inbounds begin" << endl
-                     << jacobian_output.str()
-		     << "end" << endl
-                     << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
-                     << symbol_table.param_nbr() << ");" << endl
-		     << "@inbounds begin" << endl
-                     << hessian_output.str()
-		     << "end" << endl
-                     << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
-		     << "@inbounds begin" << endl
-                     << hessian1_output.str()
-		     << "end" << endl
-                     << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
-		     << "@inbounds begin" << endl
-                     << third_derivs_output.str()
-		     << "end" << endl
-                     << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
-		     << "@inbounds begin" << endl
-                     << third_derivs1_output.str()
-		     << "end" << endl
-                     << "(rp, gp, rpp, gpp, hp)" << endl
-                     << "end" << endl
-                     << "end" << endl;
-
-  paramsDerivsFile.close();
-}
-
 void
 StaticModel::writeJsonOutput(ostream &output) const
 {
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 28ecc8cc38b720dd4b42763e37f1b3f208994102..d781d0a35d532e2987681ead6cc9edbc751dcdfe 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -149,7 +149,8 @@ public:
   void writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const;
 
   //! Writes file containing static parameters derivatives
-  void writeParamsDerivativesFile(const string &basename, bool julia) const;
+  template<bool julia>
+  void writeParamsDerivativesFile(const string &basename) const;
 
   //! Writes LaTeX file with the equations of the static model
   void writeLatexFile(const string &basename, bool write_equation_tags) const;
@@ -171,4 +172,119 @@ public:
   void addAllParamDerivId(set<int> &deriv_id_set) override;
 };
 
+template<bool julia>
+void
+StaticModel::writeParamsDerivativesFile(const string &basename) const
+{
+  if (!params_derivatives.size())
+    return;
+
+  constexpr ExprNodeOutputType output_type { julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel };
+
+  auto [tt_output, rp_output, gp_output, rpp_output, gpp_output, hp_output, g3p_output]
+    { writeParamsDerivativesFileHelper<output_type>() };
+  // g3p_output is ignored
+
+  string filename { julia ? basename + "StaticParamsDerivs.jl" : packageDir(basename) + "/static_params_derivs.m" };
+  ofstream paramsDerivsFile { filename, ios::out | ios::binary };
+  if (!paramsDerivsFile.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if constexpr(!julia)
+    {
+      paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = static_params_derivs(y, x, params)" << endl
+                       << "%" << endl
+                       << "% Status : Computes derivatives of the static model with respect to the parameters" << endl
+                       << "%" << endl
+                       << "% Inputs : " << endl
+                       << "%   y         [M_.endo_nbr by 1] double    vector of endogenous variables in declaration order" << endl
+                       << "%   x         [M_.exo_nbr by 1] double     vector of exogenous variables in declaration order" << endl
+                       << "%   params    [M_.param_nbr by 1] double   vector of parameter values in declaration order" << endl
+                       << "%" << endl
+                       << "% Outputs:" << endl
+                       << "%   rp        [M_.eq_nbr by #params] double    Jacobian matrix of static model equations with respect to parameters " << endl
+                       << "%                                              Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
+                       << "%   gp        [M_.endo_nbr by M_.endo_nbr by #params] double    Derivative of the Jacobian matrix of the static model equations with respect to the parameters" << endl
+                       << "%                                                           rows: variables in declaration order" << endl
+                       << "%                                                           rows: equations in order of declaration" << endl
+                       << "%   rpp       [#second_order_residual_terms by 4] double   Hessian matrix of second derivatives of residuals with respect to parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: number of the first parameter in derivative" << endl
+                       << "%                                                              3rd column: number of the second parameter in derivative" << endl
+                       << "%                                                              4th column: value of the Hessian term" << endl
+                       << "%   gpp      [#second_order_Jacobian_terms by 5] double   Hessian matrix of second derivatives of the Jacobian with respect to the parameters;" << endl
+                       << "%                                                              rows: respective derivative term" << endl
+                       << "%                                                              1st column: equation number of the term appearing" << endl
+                       << "%                                                              2nd column: column number of variable in Jacobian of the static model" << endl
+                       << "%                                                              3rd column: number of the first parameter in derivative" << endl
+                       << "%                                                              4th column: number of the second parameter in derivative" << endl
+                       << "%                                                              5th column: value of the Hessian term" << endl
+                       << "%" << endl
+                       << "%" << endl
+                       << "% Warning : this file is generated automatically by Dynare" << endl
+                       << "%           from model file (.mod)" << endl << endl
+                       << "T = NaN(" << params_derivs_temporary_terms_idxs.size() << ",1);" << endl
+                       << tt_output.str()
+                       << "rp = zeros(" << equations.size() << ", "
+                       << symbol_table.param_nbr() << ");" << endl
+                       << rp_output.str()
+                       << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
+                       << symbol_table.param_nbr() << ");" << endl
+                       << gp_output.str()
+                       << "if nargout >= 3" << endl
+                       << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+                       << rpp_output.str()
+                       << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+                       << gpp_output.str()
+                       << "end" << endl
+                       << "if nargout >= 5" << endl
+                       << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+                       << hp_output.str()
+                       << "end" << endl
+                       << "end" << endl;
+    }
+  else
+    paramsDerivsFile << "module " << basename << "StaticParamsDerivs" << endl
+                     << "#" << endl
+                     << "# NB: this file was automatically generated by Dynare" << endl
+                     << "#     from " << basename << ".mod" << endl
+                     << "#" << endl
+                     << "export params_derivs" << endl << endl
+                     << "function params_derivs(y, x, params)" << endl
+		     << "@inbounds begin" << endl
+                     << tt_output.str()
+		     << "end" << endl
+                     << "rp = zeros(" << equations.size() << ", "
+                     << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
+                     << rp_output.str()
+		     << "end" << endl
+                     << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
+                     << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
+                     << gp_output.str()
+		     << "end" << endl
+                     << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+		     << "@inbounds begin" << endl
+                     << rpp_output.str()
+		     << "end" << endl
+                     << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
+                     << gpp_output.str()
+		     << "end" << endl
+                     << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
+                     << hp_output.str()
+		     << "end" << endl
+                     << "(rp, gp, rpp, gpp, hp)" << endl
+                     << "end" << endl
+                     << "end" << endl;
+
+  paramsDerivsFile.close();
+}
+
 #endif