  We develop a computational method for expected functionals of the drawdown
and its duration in exponential L\'evy models. It is based on a novel
simulation algorithm for the joint law of the state, supremum and time the
supremum is attained of the Gaussian approximation of a general L\'evy process.
We bound the bias for various locally Lipschitz and discontinuous payoffs
arising in applications and analyse the computational complexities of the
corresponding Monte Carlo and multilevel Monte Carlo estimators. Monte Carlo
methods for L\'evy processes (using Gaussian approximation) have been analysed
for Lipschitz payoffs, in which case the computational complexity of our
algorithm is up to two orders of magnitude smaller when the jump activity is
high. At the core of our approach are bounds on certain Wasserstein distances,
obtained via the novel SBG coupling between a L\'evy process and its Gaussian
approximation. Numerical performance, based on the implementation in the
dedicated GitHub repository, exhibits a good agreement with our theoretical
bounds.
