'Probit model Hessian matrix of the log-likelihood
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`
Notes
.. math:: \frac{\partial^{2}\ln L}{\partial\beta\partial\beta^{\prime}}=-\lambda_{i}\left(\lambda_{i}+x_{i}^{\prime}\beta\right)x_{i}x_{i}^{\prime}
where
.. math:: \lambda_{i}=\frac{q_{i}\phi\left(q_{i}x_{i}^{\prime}\beta\right)}{\Phi\left(q_{i}x_{i}^{\prime}\beta\right)}
and :math:`q=2y-1`'

def hessian(self, params):

    X = self.exog
    XB = np.dot(X, params)
    q = ((2 * self.endog) - 1)
    L = ((q * self.pdf((q * XB))) / self.cdf((q * XB)))
    return np.dot((((- L) * (L + XB)) * X.T), X)
