code: def approx_hess1(
      x, f, epsilon=None, args=(), kwargs={}, return_grad=False
  ):
      n = len(x)
      h = _get_epsilon(x, 3, epsilon, n)
      ee = np.diag(h)
  
      f0 = f(*((x,) + args), **kwargs)
      # Compute forward step
      g = np.zeros(n)
      for i in range(n):
          g[i] = f(*((x + ee[i, :],) + args), **kwargs)
  
      hess = np.outer(h, h)  # this is now epsilon**2
      # Compute "double" forward step
      for i in range(n):
          for j in range(i, n):
              hess[i, j] = (
                  f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs)
                  - g[i]
                  - g[j]
                  + f0
              ) / hess[i, j]
              hess[j, i] = hess[i, j]
      if return_grad:
          grad = (g - f0) / h
          return hess, grad
      else:
          return hess
  
file_path: statsmodels/statsmodels/tools/numdiff.py
imports: import numpy
  from statsmodels.compat.pandas import Appender, Substitution
label_desc: Use image node_id 7 for calling a global function with example usage: approx_hess1(x, f, epsilon, args, kwargs, return_grad) and returns: hess,  grad, hess