score of the model (gradient of log-likelihood at theta=theta_0)

def psit(t, xp, x):
    """ score of the model (gradient of log-likelihood at theta=theta_0)
    """
    if t == 0:
        return -0.5 / sigma0**2 + \
            (0.5 * (1. - phi0**2) / sigma0**4) * (x - mu0)**2
    else:
        return -0.5 / sigma0**2 + (0.5 / sigma0**4) * \
            ((x - mu0) - phi0 * (xp - mu0))**2