​When approximating the expectation of a functional of a certain stochastic process, the efficiency and performance of deterministic quadrature methods, and hierarchical variance reduction methods such as multilevel Monte Carlo (MLMC), is highly deteriorated in different ways by the low regularity of the integrand with respect to the input parameters. To overcome this issue, a smoothing procedure is needed to uncover the available regularity and improve the performance of the aforementioned methods. In this work, we consider cases where we cannot perform an analytic smoothing. Thus, we introduce a novel numerical smoothing technique based on root-finding combined with a one dimensional integration with respect to a single well-chosen variable. We prove that under appropriate conditions, the resulting function of the remaining variables is highly smooth, potentially allowing a higher efficiency of adaptive sparse grids quadrature (ASGQ), in particular when combined with hierarchical representations to treat the high dimensionality effectively. Our study is motivated by option pricing problems and our main focus is on dynamics where a discretization of the asset price is needed. Our analysis and numerical experiments illustrate the advantage of combining numerical smoothing with ASGQ compared to the Monte Carlo method. Furthermore, we demonstrate how numerical smoothing significantly reduces the kurtosis at the deep levels of MLMC, and also improves the strong convergence rate, when using Euler scheme. Due to the complexity theorem of MLMC, and given a pre-selected tolerance, $\text{TOL}$, this results in an improvement of the complexity from $\Ordo{\text{TOL}^{-2.5}}$ in the standard case to $\Ordo{\text{TOL}^{-2} \log(\text{TOL})^2}$. Finally, we show how our numerical smoothing combined with MLMC enables us also to estimate density functions, which standard MLMC (without smoothing) fails to achieve.