  The multilevel Monte Carlo (MLMC) method is highly efficient for estimating
expectations of a functional of a solution to a stochastic differential
equation (SDE). However, MLMC estimators may be unstable and have a poor
(noncanonical) complexity in the case of low regularity of the functional. To
overcome this issue, we extend our previously introduced idea of numerical
smoothing in (Quantitative Finance, 23(2), 209-227, 2023), in the context of
deterministic quadrature methods to the MLMC setting. The numerical smoothing
technique is based on root-finding methods combined with one-dimensional
numerical integration with respect to a single well-chosen variable. This study
is motivated by the computation of probabilities of events, pricing options
with a discontinuous payoff, and density estimation problems for dynamics where
the discretization of the underlying stochastic processes is necessary. The
analysis and numerical experiments reveal that the numerical smoothing
significantly improves the strong convergence, and consequently, the complexity
and robustness (by making the kurtosis at deep levels bounded) of the MLMC
method. In particular, we show that numerical smoothing enables recovering the
MLMC complexities obtained for Lipschitz functionals due to the optimal
variance decay rate when using the Euler--Maruyama scheme. For the Milstein
scheme, numerical smoothing recovers the canonical MLMC complexity even for the
nonsmooth integrand mentioned above. Finally, our approach efficiently
estimates univariate and multivariate density functions.
