We design and implement a novel algorithm for computing a multilevel Monte Carlo (MLMC) estimator of the cumulative distribution function of a quantity of interest in problems with random input parameters or initial conditions. Our approach combines a standard MLMC method with stratified sampling by replacing standard Monte Carlo at each level with stratified Monte Carlo with proportional allocation. We show that the resulting stratified MLMC algorithm is more efficient than its standard MLMC counterpart, due to the reduction in variance at each level provided by the stratification of the random parameter's domain. A smoothing approximation for the indicator function based on kernel density estimation yields a more efficient algorithm compared to the typically used polynomial smoothing. The difference in computational cost between the smoothing methods depends on the required error tolerance.