'Generalized Poisson model Hessian matrix of the loglikelihood
Parameters
params : array-like
The parameters of the model
Returns
hess : ndarray, (k_vars, k_vars)
The Hessian, second derivative of loglikelihood function,
evaluated at `params`'

def hessian(self, params):

    if self._transparams:
        alpha = np.exp(params[(-1)])
    else:
        alpha = params[(-1)]
    params = params[:(-1)]
    p = self.parameterization
    exog = self.exog
    y = self.endog[:, None]
    mu = self.predict(params)[:, None]
    mu_p = np.power(mu, p)
    a1 = (1 + (alpha * mu_p))
    a2 = (mu + ((alpha * mu_p) * y))
    a3 = ((alpha * p) * (mu ** (p - 1)))
    a4 = (a3 * y)
    a5 = (p * (mu ** (p - 1)))
    dmudb = (mu * exog)
    dim = exog.shape[1]
    hess_arr = np.empty(((dim + 1), (dim + 1)))
    for i in range(dim):
        for j in range((i + 1)):
            hess_arr[(i, j)] = np.sum((((mu * exog[:, i, None]) * exog[:, j, None]) * ((mu * ((((((((((((a3 * a4) / (a1 ** 2)) - (((2 * (a3 ** 2)) * a2) / (a1 ** 3))) + (((2 * a3) * (a4 + 1)) / (a1 ** 2))) - ((a4 * p) / (mu * a1))) + (((a3 * p) * a2) / (mu * (a1 ** 2)))) + (a4 / (mu * a1))) - ((a3 * a2) / (mu * (a1 ** 2)))) + ((((y - 1) * a4) * (p - 1)) / (a2 * mu))) - (((y - 1) * ((1 + a4) ** 2)) / (a2 ** 2))) - ((a4 * (p - 1)) / (a1 * mu))) - (1 / (mu ** 2)))) + ((((((- a4) / a1) + ((a3 * a2) / (a1 ** 2))) + (((y - 1) * (1 + a4)) / a2)) - ((1 + a4) / a1)) + (1 / mu)))), axis=0)
    tri_idx = np.triu_indices(dim, k=1)
    hess_arr[tri_idx] = hess_arr.T[tri_idx]
    dldpda = np.sum(((((((((((2 * a4) * mu_p) / (a1 ** 2)) - ((((2 * a3) * mu_p) * a2) / (a1 ** 3))) - ((((mu_p * y) * (y - 1)) * (1 + a4)) / (a2 ** 2))) + ((mu_p * (1 + a4)) / (a1 ** 2))) + (((a5 * y) * (y - 1)) / a2)) - (((2 * a5) * y) / a1)) + ((a5 * a2) / (a1 ** 2))) * dmudb), axis=0)
    hess_arr[(-1), :(-1)] = dldpda
    hess_arr[:(-1), (-1)] = dldpda
    dldada = ((mu_p ** 2) * ((((3 * y) / (a1 ** 2)) - (((y ** 2) * (y - 1)) / (a2 ** 2))) - ((2 * a2) / (a1 ** 3))))
    hess_arr[((-1), (-1))] = dldada.sum()
    return hess_arr
