Geometrically Convergent Simulation of the Extrema of L\'{e}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\'evy process at an arbitrary finite time. We identify the law of the error in simple terms. We prove that the error decays geometrically in $L^p$ (for any $p\geq 1$) as a function of the computational cost, in contrast with the polynomial decay for the approximations available in the literature. We establish a central limit theorem and construct non-asymptotic and asymptotic confidence intervals for the corresponding Monte Carlo estimator. We prove that the multilevel Monte Carlo estimator has optimal computational complexity (i.e. of order $\epsilon^{-2}$ if the mean squared error is at most $\epsilon^2$) for locally Lipschitz and barrier-type functionals of the triplet and develop an unbiased version of the estimator. We illustrate the performance of the algorithm with numerical examples.