code: def Ikpw(dW, h, n=5):
      
      N = dW.shape[0]
      m = dW.shape[1]
      if dW.ndim < 3:
          dW = dW.reshape((N, -1, 1)) # change to array of shape (N, m, 1)
      if dW.shape[2] != 1 or dW.ndim > 3:
          raise(ValueError)
      A = _Aterm(N, h, m, 1, dW)
      for k in range(2, n+1):
          A += _Aterm(N, h, m, k, dW)
      A = (h/(2.0*np.pi))*A
      I = 0.5*(_dot(dW, _t(dW)) - np.diag(h*np.ones(m))) + A
      dW = dW.reshape((N, -1)) # change back to shape (N, m)
      return (A, I)
docs: matrix I approximating repeated Ito 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, I) where
          A: array of shape (N, m, m) giving the Levy areas that were used.
          I: array of shape (N, m, m) giving an m x m matrix of repeated Ito 
          integral values for each of the N time intervals.