We introduce an algorithm based on semidefinite programming that yields increasing (resp.
decreasing) sequences of lower (resp. upper) bounds on polynomial stationary averages of diffusions
with polynomial drift vector and diffusion coefficients. The bounds are obtained by
optimising an objective, determined by the stationary average of interest, over the set of real
vectors defined by certain linear equalities and semidefinite inequalities which are satisfied by
the moments of any stationary measure of the diffusion. We exemplify the use of the approach
through several applications: a Bayesian inference problem; the computation of Lyapunov exponents
of linear ordinary differential equations perturbed by multiplicative white noise; and a
reliability problem from structural mechanics. Additionally, we prove that the bounds converge
to the infimum and supremum of the set of stationary averages for certain SDEs associated with
the computation of the Lyapunov exponents, and we provide numerical evidence of convergence
in more general settings.