Geometrically convergent simulation of the extrema of Lévy processes

We develop a novel approximate simulation algorithm for 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. The law of the error is identifiable in simple terms and this error decays geometrically in L^p (for any p≥1) as a function of the computational cost, in contrast with the polynomial decay for the approximations available in the literature. This results in a multilevel Monte Carlo estimator that has optimal computational complexity (i.e. of order ϵ^{−2} if the mean squared error is at most ϵ^2) even for barrier-type functionals of the triplet. We illustrate the performance of the algorithm with numerical examples and discuss further results in multiple directions.