  In this paper, we present a minimal formalism for Stein operators which leads
to different probabilistic representations of solutions to Stein equations.
These in turn provide a wide family of Stein-Covariance identities which we put
to use for revisiting the very classical topic of bounding the variance of
functionals of random variables. Applying the Cauchy-Schwarz inequality yields
first order upper and lower Klaassen-type variance bounds. A probabilistic
representation of Lagrange's identity (i.e. Cauchy-Schwarz with remainder)
leads to Papathanasiou-type variance expansions of arbitrary order. A matrix
Cauchy-Schwarz inequality leads to Olkin-Shepp type covariance bounds. All
results hold for univariate target distribution under very weak assumptions (in
particular they hold for continuous and discrete distributions alike). Many
concrete illustrations are provided.
