  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.
