**Title: **Geometrically convergent simulation of the extrema of Lévy processes

**Abstract: **We develop a novel Monte Carlo algorithm for the simulation from the joint law of the position, the running supremum and the time of the supremum of a general Lévy process at an arbitrary finite time. We prove that the bias decays geometrically, in contrast to the power law for the random walk approximation (RWA). We identify the law of the error and, inspired by the recent work of Ivanovs (Zooming in on a Lévy process at its supremum, Ann. Appl. Probab. 28 (2018)) on RWA, characterise its asymptotic behaviour. We prove that the multilevel Monte Carlo (MLMC) estimator has optimal computational complexity (i.e. of order 1/ε^2 if the L2-norm of the error is at most ε) for locally Lipschitz and barrier-type functionals of the triplet. If the increments of the Lévy process cannot be sampled directly, we combine our algorithm with the Asmussen-Rosiński approximation by choosing the rate of decay of the cutoff level for small jumps so that the corresponding MC and MLMC estimators have minimal computational complexity. **arXiv:1810.11039**