'Evaluate the Hessian function at a given point.
Parameters
params : array-like
The parameter vector at which the Hessian is computed.
scale : float or None
If None, return the profile (concentrated) log likelihood
(profiled over the scale parameter), else return the
log-likelihood using the given scale value.
Returns
The Hessian matrix.'

def hessian(self, params, scale=None):

    if (not hasattr(self, '_wexog_xprod')):
        self._setup_score_hess()
    xtxb = np.dot(self._wexog_xprod, params)
    if (scale is None):
        ssr = (self._wendog_xprod - (2 * np.dot(self._wexog_x_wendog.T, params)))
        ssr += np.dot(params, xtxb)
        ssrp = (((-2) * self._wexog_x_wendog) + (2 * xtxb))
        hm = ((self._wexog_xprod / ssr) - (np.outer(ssrp, ssrp) / (ssr ** 2)))
        return (((- self.nobs) * hm) / 2)
    else:
        return ((- self._wexog_xprod) / scale)
