'Multinomial logit Hessian matrix of the log-likelihood
Parameters
params : array-like
The parameters of the model
Returns
hess : ndarray, (J*K, J*K)
The Hessian, second derivative of loglikelihood function with
respect to the flattened parameters, evaluated at `params`
Notes
.. math:: \frac{\partial^{2}\ln L}{\partial\beta_{j}\partial\beta_{l}}=-\sum_{i=1}^{n}\frac{\exp\left(\beta_{j}^{\prime}x_{i}\right)}{\sum_{k=0}^{J}\exp\left(\beta_{k}^{\prime}x_{i}\right)}\left[\boldsymbol{1}\left(j=l\right)-\frac{\exp\left(\beta_{l}^{\prime}x_{i}\right)}{\sum_{k=0}^{J}\exp\left(\beta_{k}^{\prime}x_{i}\right)}\right]x_{i}x_{l}^{\prime}
where
:math:`\boldsymbol{1}\left(j=l\right)` equals 1 if `j` = `l` and 0
otherwise.
The actual Hessian matrix has J**2 * K x K elements. Our Hessian
is reshaped to be square (J*K, J*K) so that the solvers can use it.
This implementation does not take advantage of the symmetry of
the Hessian and could probably be refactored for speed.'

def hessian(self, params):

    params = params.reshape(self.K, (-1), order='F')
    X = self.exog
    pr = self.cdf(np.dot(X, params))
    partials = []
    J = self.J
    K = self.K
    for i in range((J - 1)):
        for j in range((J - 1)):
            if (i == j):
                partials.append((- np.dot(((pr[:, (i + 1)] * (1 - pr[:, (j + 1)]))[:, None] * X).T, X)))
            else:
                partials.append((- np.dot(((pr[:, (i + 1)] * (- pr[:, (j + 1)]))[:, None] * X).T, X)))
    H = np.array(partials)
    H = np.transpose(H.reshape((J - 1), (J - 1), K, K), (0, 2, 1, 3)).reshape(((J - 1) * K), ((J - 1) * K))
    return H
