code: def psit(t, xp, x):
      
      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
docs: score of the model (gradient of log-likelihood at theta=theta_0)