def approx_hess_cs(x, f, epsilon=None, args=(), kwargs={}):
    # TODO: might want to consider lowering the step for pure derivatives
    n = len(x)
    h = _get_epsilon(x, 3, epsilon, n)
    ee = np.diag(h)
    hess = np.outer(h, h)

    n = len(x)

    for i in range(n):
        for j in range(i, n):
            hess[i, j] = (f(*((x + 1j*ee[i, :] + ee[j, :],) + args), **kwargs)
                          - f(*((x + 1j*ee[i, :] - ee[j, :],)+args),
                              **kwargs)).imag/2./hess[i, j]
            hess[j, i] = hess[i, j]

    return hess
Calculate Hessian with complex-step derivative approximation
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/tools/numdiff.py/approx_hess_cs