docs: def Iwik(dW, h, n=5):
      """matrix I approximating repeated Ito integrals for each of N time
      intervals, using the method of Wiktorsson (2001).
  
      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:
        (Atilde, I) where
          Atilde: array of shape (N,m(m-1)//2,1) giving the area integrals 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.
      """
      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)
      if m == 1:
          return (np.zeros((N, 1, 1)), (dW*dW - h)/2.0)
      Pm0 = _P(m)
      Km0 = _K(m)
      M = m*(m-1)//2
      Atilde_n = _AtildeTerm(N, h, m, 1, dW, Km0, Pm0)
      for k in range(2, n+1):
          Atilde_n += _AtildeTerm(N, h, m, k, dW, Km0, Pm0)
      Atilde_n = (h/(2.0*np.pi))*Atilde_n # approximation after n terms
      S = _sigmainf(N, h, m, dW, Km0, Pm0)
      normdW2 = np.sum(np.abs(dW)**2, axis=1)
      radical = np.sqrt(1.0 + normdW2/h).reshape((N, 1, 1))
      IM = broadcast_to(np.eye(M), (N, M, M))
      Im = broadcast_to(np.eye(m), (N, m, m))
      Ims0 = np.eye(m**2)
      sqrtS = (S + 2.0*radical*IM)/(np.sqrt(2.0)*(1.0 + radical))
      G = np.random.normal(0.0, 1.0, (N, M, 1))
      tailsum = h/(2.0*np.pi)*_a(n)**0.5*_dot(sqrtS, G)
      Atilde = Atilde_n + tailsum # our final approximation of the areas
      factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
      vecI = 0.5*(_kp(dW, dW) - _vec(h*Im)) + _dot(factor3, Atilde)
      I = _unvec(vecI)
      dW = dW.reshape((N, -1)) # change back to shape (N, m)
      return (Atilde, I)