diff --git a/DynamicModel.cc b/DynamicModel.cc
index 5a62d00ddad409f0b3acc58873637624009c7dd5..76b5988f97cd46667cf65d3bbc566987dca7a37a 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -2055,7 +2055,6 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       int var = it->first.second;
       expr_t d1 = it->second;
 
-      jacobian_output << "g1";
       jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
       jacobian_output << "=";
       d1->writeOutput(jacobian_output, output_type, temporary_terms, tef_terms);
@@ -2178,18 +2177,18 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
                     << "if nargout >= 2," << endl
                     << "  g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl
                     << endl
-                    << "%" << endl
-                    << "% Jacobian matrix" << endl
-                    << "%" << endl
+                    << "  %" << endl
+                    << "  % Jacobian matrix" << endl
+                    << "  %" << endl
                     << endl
                     << jacobian_output.str()
                     << "end" << endl;
 
       // Initialize g2 matrix
       DynamicOutput << "if nargout >= 3," << endl
-                    << "%" << endl
-                    << "% Hessian matrix" << endl
-                    << "%" << endl
+                    << "  %" << endl
+                    << "  % Hessian matrix" << endl
+                    << "  %" << endl
                     << endl;
       if (second_derivatives.size())
         DynamicOutput << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl
@@ -2201,9 +2200,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
 
       // Initialize g3 matrix
       DynamicOutput << "if nargout >= 4," << endl
-                    << "%" << endl
-                    << "% Third order derivatives" << endl
-                    << "%" << endl
+                    << "  %" << endl
+                    << "  % Third order derivatives" << endl
+                    << "  %" << endl
                     << endl;
       int ncols = hessianColsNbr * dynJacobianColsNbr;
       if (third_derivatives.size())
@@ -3690,28 +3689,6 @@ DynamicModel::writeLatexFile(const string &basename) const
   writeLatexModelFile(basename + "_dynamic.tex", oLatexDynamicModel);
 }
 
-void
-DynamicModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
-{
-  output << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
-    output << eq_nb + 1 << "," << col_nb + 1;
-  else
-    output << eq_nb + col_nb *equations.size();
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
-void
-DynamicModel::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
-{
-  output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
-    output << row_nb + 1 << "," << col_nb + 1;
-  else
-    output << row_nb + col_nb * NNZDerivatives[order-1];
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
 void
 DynamicModel::substituteEndoLeadGreaterThanTwo(bool deterministic_model)
 {
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 8298ece8d5635e11671c7d92d11a291a27d06fc2..f26f93b98a5c6c1632ce1cf0b9f8ce77d6ab9ebc 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -168,15 +168,6 @@ private:
   /*! Also computes max_{endo,exo}_{lead_lag}, and initializes dynJacobianColsNbr to the number of dynamic endos */
   void computeDerivIDs();
 
-  //! Helper for writing the Jacobian elements in MATLAB and C
-  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
-  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
-
-  //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
-  /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]]
-    If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */
-  void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
-
   //! Write chain rule derivative of a recursive equation w.r. to a variable
   void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
 
diff --git a/ModelTree.cc b/ModelTree.cc
index 90b19941fece38bbe4893d4984042c6aa7f9f133..eb9b749b9a96dcf03bea1d7e1faa0942b8125c7c 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -1423,3 +1423,25 @@ ModelTree::set_cutoff_to_zero()
 {
   cutoff = 0;
 }
+
+void
+ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
+{
+  output << "  g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
+  if (IS_MATLAB(output_type))
+    output << eq_nb + 1 << "," << col_nb + 1;
+  else
+    output << eq_nb + col_nb *equations.size();
+  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
+}
+
+void
+ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
+{
+  output << "  v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
+  if (IS_MATLAB(output_type))
+    output << row_nb + 1 << "," << col_nb + 1;
+  else
+    output << row_nb + col_nb * NNZDerivatives[order-1];
+  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
+}
diff --git a/ModelTree.hh b/ModelTree.hh
index 72fddb753cd0af6bd7bfe6cf5735f770f551d6ae..37413ae6bb3713ef350ee12a4605f8c60ff42e15 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -262,7 +262,14 @@ public:
   //! Adds a nonstationary variable with its deflator
   void addNonstationaryVariables(vector<int> nonstationary_vars, expr_t deflator) throw (TrendException);
   void set_cutoff_to_zero();
-  
+  //! Helper for writing the Jacobian elements in MATLAB and C
+  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
+  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
+  //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
+  /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]]
+    If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */
+  void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
+
   inline static std::string
   c_Equation_Type(int type)
   {
diff --git a/StaticModel.cc b/StaticModel.cc
index b30fe30bf250cd6bae1f1801dfe9665d0d33492b..09bfeb573dca292ba3af85b9d4178047830f4a59 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -1168,7 +1168,6 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       int symb_id = it->first.second;
       expr_t d1 = it->second;
 
-      jacobian_output << "  g1";
       jacobianHelper(jacobian_output, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
       jacobian_output << "=";
       d1->writeOutput(jacobian_output, output_type, temporary_terms, tef_terms);
@@ -1192,13 +1191,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
       int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
 
-      hessianHelper(hessian_output, k, 0, output_type);
+      sparseHelper(2, hessian_output, k, 0, output_type);
       hessian_output << "=" << eq + 1 << ";" << endl;
 
-      hessianHelper(hessian_output, k, 1, output_type);
+      sparseHelper(2, hessian_output, k, 1, output_type);
       hessian_output << "=" << col_nb + 1 << ";" << endl;
 
-      hessianHelper(hessian_output, k, 2, output_type);
+      sparseHelper(2, hessian_output, k, 2, output_type);
       hessian_output << "=";
       d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
       hessian_output << ";" << endl;
@@ -1208,15 +1207,15 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       // Treating symetric elements
       if (symb_id1 != symb_id2)
         {
-          hessianHelper(hessian_output, k, 0, output_type);
+          sparseHelper(2, hessian_output, k, 0, output_type);
           hessian_output << "=" << eq + 1 << ";" << endl;
 
-          hessianHelper(hessian_output, k, 1, output_type);
+          sparseHelper(2, hessian_output, k, 1, output_type);
           hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
 
-          hessianHelper(hessian_output, k, 2, output_type);
+          sparseHelper(2, hessian_output, k, 2, output_type);
           hessian_output << "=";
-          hessianHelper(hessian_output, k-1, 2, output_type);
+          sparseHelper(2, hessian_output, k-1, 2, output_type);
           hessian_output << ";" << endl;
 
           k++;
@@ -1700,28 +1699,6 @@ StaticModel::writeLatexFile(const string &basename) const
   writeLatexModelFile(basename + "_static.tex", oLatexStaticModel);
 }
 
-void
-StaticModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
-{
-  output << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
-    output << eq_nb + 1 << "," << col_nb + 1;
-  else
-    output << eq_nb + col_nb *equations.size();
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
-void
-StaticModel::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
-{
-  output << "v2" << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
-    output << row_nb + 1 << ", " << col_nb + 1;
-  else
-    output << row_nb + col_nb * NNZDerivatives[1];
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
-}
-
 void
 StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const
 {
diff --git a/StaticModel.hh b/StaticModel.hh
index 81020d8aff8f63e36a00ff434646f87a772ee196..6da13047782288156eeec1ebd2af672f14e580c5 100644
--- a/StaticModel.hh
+++ b/StaticModel.hh
@@ -107,14 +107,6 @@ private:
   //! Collect only the first derivatives
   map<pair<int, pair<int, int> >, expr_t> collect_first_order_derivatives_endogenous();
 
-  //! Helper for writing the Jacobian elements in MATLAB and C
-  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
-  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
-
-  //! Helper for writing the sparse Hessian elements in MATLAB and C
-  /*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */
-  void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
-
   //! Write chain rule derivative of a recursive equation w.r. to a variable
   void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;