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)