code: def psit(t, x, xf):
      
      if t == 0:
          return psi0(x) + psit(1, x, xf)
      else:
          return -0.5 / sigma0**2 + (0.5 / sigma0**4) * ((xf - mu0) - 
                                                         phi0 * (x - mu0))**2
docs: A function of t, X_t and X_{t+1} (f=future)