Universal Monte Carlo method for Lévy processes and their extrema

The article proposes a universal approach to constructing Monte Carlo methods for pricing options with payoffs depending on the joint distribution of the final position of the Lévy process $X_T$ and its infimum $\mathcal{I}_T$ (or supremum $\mathcal{S}_T$). We derive approximate formulas for the conditional cumulative distribution functions of the Lévy process ${\mathbf{P}(X_{T}<x|\mathcal{S}_{T}=y)}$ ($\mathbf{P}(X_{T}<x| \mathcal{I}_{T}=y)$), which are expressed through the partial derivative with respect to $y$ of the joint cumulative distribution function $\mathbf{P}(X_{T}\!<x, \mathcal{S}_{T}\!< y)$ ($\mathbf{P}(X_{T}\!< x, \mathcal{I}_{T}\!< y)\!$) and the density of the infimum (or supremum) at the final moment of time. By applying the Laplace transform to the joint cumulative distribution function of the Lévy process and its extremum, we use the approximate Wiener–Hopf factorization to represent the image of its partial derivative. By inverting the Laplace transform using the Gaver–Stehfest algorithm, we find the desired conditional cumulative distribution function. The developed algorithm for simulating the joint position of the Lévy process and its extremum at a given point in time consists of two key stages. At the first stage, we simulate the extremum value of the Lévy process based on the approximation of its cumulative distribution function $\mathbf{P}(\mathcal{S}_{T}<x)$ (or $\mathbf{P}(\mathcal{I}_{T}<x)$). In the second step, we simulate the final value of the Lévy process based on the approximation of the conditional cumulative distribution function of the final position of the Lévy process relative to its extremum. The universality of the Monte Carlo method we developed lies in the implementation of a uniform approach for a wide class of Lévy processes, in contrast to classical approaches, when simulations are essentially based on the features of the probability distribution associated with the simulated random process or its extrema. In our approach, it is enough to know the characteristic exponent of the Lévy process. The most time-consuming computational unit for simulating a random variable based on a known cumulative distribution function can be effectively implemented using neural networks and accelerated through parallel computing. Thus, on the one hand, the approach we propose is suitable for a wide class of Lévy models, on the other hand, it can be combined with machine learning methods.