We consider all two-times iterated Itô integrals obtained by pairing m independent standard Brownian motions. First we calculate the conditional joint characteristic function of these integrals, given the Brownian increments over the integration interval, and show that it has a form entirely similar to what is obtained in the univariate case. Then we propose an algorithm for the simultaneous simulation of the $m^2$ integrals conditioned on the Brownian increments that achieves a mean square error of order $1/n^2$, where n is the number of terms in a truncated sum. The algorithm is based on approximation of the tail-sum distribution, which is a multivariate normal variance mixture, by a multivariate normal distribution.
"Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions." Ann. Appl. Probab. 11 (2) 470 - 487, May 2001. https://doi.org/10.1214/aoap/1015345301