Higher-order spring-coupled multilevel Monte Carlo method for invariant measures

A higher-order change-of-measure multilevel Monte Carlo (MLMC) method is developed for computing weak approximations of the invariant measures of SDE with drift coefficients that do not satisfy the contractivity condition. This is achieved by introducing a spring term in the pairwise coupling of the MLMC trajectories employing the order 1.5 strong It^o--Taylor method. Through this, we can recover the contractivity property of the drift coefficient while still retaining the telescoping sum property needed for implementing the MLMC method. We show that the variance of the change-of-measure MLMC method grows linearly in time $T$ for all $T > 0$, and for all sufficiently small timestep size $h > 0$. For a given error tolerance $\epsilon > 0$, we prove that the method achieves a mean-square-error accuracy of $O(\epsilon^2)$ with a computational cost of $O(\epsilon^{-2} \big\vert \log \epsilon \big\vert^{3/2} (\log \big\vert \log \epsilon \big\vert)^{1/2})$ for uniformly Lipschitz continuous payoff functions and $O \big( \epsilon^{-2} \big\vert \log \epsilon \big\vert^{5/3 + \xi} \big)$ for discontinuous payoffs, respectively, where $\xi > 0$. We also observe an improvement in the constant associated with the computational cost of the higher-order change-of-measure MLMC method, marking an improvement over the Milstein change-of-measure method in the aforementioned seminal work by M. Giles and W. Fang. Several numerical tests were performed to verify the theoretical results and assess the robustness of the method.