'score test for restrictions or for omitted variables
The covariance matrix for the score is based on the Hessian, i.e.
observed information matrix or optionally on the expected information
matrix..
Parameters
params_constrained : array_like
estimated parameter of the restricted model. This can be the
parameter estimate for the current when testing for omitted
variables.
k_constraints : int or None
Number of constraints that were used in the estimation of params
restricted relative to the number of exog in the model.
This must be provided if no exog_extra are given. If exog_extra is
not None, then k_constraints is assumed to be zero if it is None.
exog_extra : None or array_like
Explanatory variables that are jointly tested for inclusion in the
model, i.e. omitted variables.
observed : bool
If True, then the observed Hessian is used in calculating the
covariance matrix of the score. If false then the expected
information matrix is used.
Returns
chi2_stat : float
chisquare statistic for the score test
p-value : float
P-value of the score test based on the chisquare distribution.
df : int
Degrees of freedom used in the p-value calculation. This is equal
to the number of constraints.
Notes
not yet verified for case with scale not equal to 1.'

def score_test(self, params_constrained, k_constraints=None, exog_extra=None, observed=True):

    if (exog_extra is None):
        if (k_constraints is None):
            raise ValueError('if    exog_extra    is    None,    then    k_constraintsneeds    to    be    given')
        score = self.score(params_constrained)
        hessian = self.hessian(params_constrained, observed=observed)
    else:
        if (k_constraints is None):
            k_constraints = 0
        ex = np.column_stack((self.exog, exog_extra))
        k_constraints += (ex.shape[1] - self.exog.shape[1])
        score_factor = self.score_factor(params_constrained)
        score = (score_factor[:, None] * ex).sum(0)
        hessian_factor = self.hessian_factor(params_constrained, observed=observed)
        hessian = (- np.dot((ex.T * hessian_factor), ex))
    from scipy import stats
    chi2stat = (- score.dot(np.linalg.solve(hessian, score[:, None])))
    pval = stats.chi2.sf(chi2stat, k_constraints)
    return (chi2stat, pval, k_constraints)
