docs: def Jkpw(dW, h, n=5):
      """matrix J approximating repeated Stratonovich integrals for each of N
      time intervals, based on the method of Kloeden, Platen and Wright (1992).
  
      Args:
        dW (array of shape (N, m)): giving m independent Weiner increments for
          each time step N. (You can make this array using sdeint.deltaW())
        h (float): the time step size
        n (int, optional): how many terms to take in the series expansion
  
      Returns:
        (A, J) where
          A: array of shape (N, m, m) giving the Levy areas that were used.
          J: array of shape (N, m, m) giving an m x m matrix of repeated
          Stratonovich integral values for each of the N time intervals.
      """
      m = dW.shape[1]
      A, I = Ikpw(dW, h, n)
      J = I + 0.5*h*np.eye(m).reshape((1, m, m))
      return (A, J)