code: def deriv2(self, p):
          
          from statsmodels.tools.numdiff import approx_fprime
          p = np.atleast_1d(p)
          # Note: special function for norm.ppf does not support complex
          return np.diag(approx_fprime(p, self.deriv, centered=True))
docs: Second derivative of the link function g''(p)
  
          implemented through numerical differentiation