'Notes
Need to generalize for AR(p) and for a constant.
Not correct yet.  Returns numerical gradient.  Depends on package
numdifftools.'

def score(self, params):

    y = self.endog
    ylag = self.exog
    nobs = self.nobs
    diffsumsq = sumofsq((y - np.dot(ylag, params)))
    dsdr = ((((1 / nobs) * (-2)) * np.sum((ylag * (y - np.dot(ylag, params))[:, None]))) + ((2 * params) * (ylag[0] ** 2)))
    sigma2 = ((1 / nobs) * (diffsumsq - ((ylag[0] ** 2) * (1 - (params ** 2)))))
    gradient = ((((((((- nobs) / (2 * sigma2)) * dsdr) + (params / (1 - (params ** 2)))) + ((1 / sigma2) * np.sum((ylag * (y - np.dot(ylag, params))[:, None])))) + (((0.5 * (sigma2 ** (-2))) * diffsumsq) * dsdr)) + (((ylag[0] ** 2) * params) / sigma2)) + ((((ylag[0] ** 2) * (1 - (params ** 2))) / (2 * (sigma2 ** 2))) * dsdr))
    if self.penalty:
        pass
    j = Jacobian(self.loglike)
    return j(params)
