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.