In this work we introduce and analyze a novel multilevel Monte Carlo (MLMC) estimator for the accurate approximation of central moments of system outputs affected by uncertainties. Central moments play a central role in many disciplines to characterize a random system output's distribution and are of primary importance in many prediction, optimization, and decision making processes under uncertainties. We detail how to effectively tune the MLMC algorithm for central moments of any order and present a complete practical algorithm that is implemented as part of a Python library [1]. In fact, we validate the methodology on selected reference problems and apply it to an aerodynamic relevant test case, namely the transonic RAE 2822 airfoil affected by operating and geometric uncertainties.