  It is well known that, when numerically simulating solutions to SDEs,
achieving a strong convergence rate better than O(\sqrt{h}) (where h is the
step size) requires the use of certain iterated integrals of Brownian motion,
commonly referred to as its "L\'{e}vy areas". However, these stochastic
integrals are difficult to simulate due to their non-Gaussian nature and for a
d-dimensional Brownian motion with d > 2, no fast almost-exact sampling
algorithm is known.
  In this paper, we propose L\'{e}vyGAN, a deep-learning-based model for
generating approximate samples of L\'{e}vy area conditional on a Brownian
increment. Due to our "Bridge-flipping" operation, the output samples match all
joint and conditional odd moments exactly. Our generator employs a tailored
GNN-inspired architecture, which enforces the correct dependency structure
between the output distribution and the conditioning variable. Furthermore, we
incorporate a mathematically principled characteristic-function based
discriminator. Lastly, we introduce a novel training mechanism termed
"Chen-training", which circumvents the need for expensive-to-generate training
data-sets. This new training procedure is underpinned by our two main
theoretical results.
  For 4-dimensional Brownian motion, we show that L\'{e}vyGAN exhibits
state-of-the-art performance across several metrics which measure both the
joint and marginal distributions. We conclude with a numerical experiment on
the log-Heston model, a popular SDE in mathematical finance, demonstrating that
high-quality synthetic L\'{e}vy area can lead to high order weak convergence
and variance reduction when using multilevel Monte Carlo (MLMC).
