Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions

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 m2 integrals conditioned on the Brownian increments that achieves a mean square error of order 1/n2, 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.