Implementable coupling of L\'evy process and Brownian motion

We provide a simple algorithm for construction of Brownian paths approximating those of a L\'evy process on a finite time interval. It requires knowledge of the L\'evy process trajectory on a chosen regular grid and the law of its endpoint, or the ability to simulate from that. This algorithm is based on reordering of Brownian increments, and it can be applied in a recursive manner. We establish an upper bound on the mean squared maximal distance between the paths and determine a suitable mesh size in various asymptotic regimes. The analysis proceeds by reduction to the comonotonic coupling of increments. Applications to model risk and multilevel Monte Carlo are discussed in detail, and numerical examples are provided.