Many applications across sciences and technologies require a careful quantification of non-deterministic effects to a system output, for example when evaluating the system's reliability or when gearing it towards more robust operation conditions. At the heart of these considerations lies an accurate characterization of uncertain system outputs. In this work we introduce and analyze novel multilevel Monte Carlo techniques for an efficient characterization of an uncertain system output's distribution. These techniques rely on accurately approximating general parametric expectations, i.e. expectations that depend on a parameter, uniformly on an interval. Applications of interest include, for example, the approximation of the characteristic function and of the cumulative distribution function of an uncertain system output. A further important consequence of the introduced approximation techniques for parametric expectations (i.e. for functions) is that they allow to construct multilevel Monte Carlo estimators for various robustness indicators, such as for a quantile (also known as value-at-risk) and for the conditional value-at-risk. These robustness indicators cannot be expressed as moments and are thus not easily accessible usually. In fact, here we provide a framework that allows to simultaneously estimate a cumulative distribution function, a quantile, and the associated conditional value-at-risk of an uncertain system output at the cost of a single multilevel Monte Carlo simulation, while each estimated quantity satisfies a prescribed tolerance goal.