  We investigate a weighted Multilevel Richardson-Romberg extrapolation for the
ergodic approximation of invariant distributions of diffusions adapted from the
one introduced in~[Lemaire-Pag\`es, 2013] for regular Monte Carlo simulation.
In a first result, we prove under weak confluence assumptions on the diffusion,
that for any integer $R\ge2$, the procedure allows us to attain a rate
$n^{\frac{R}{2R+1}}$ whereas the original algorithm convergence is at a weak
rate $n^{1/3}$. Furthermore, this is achieved without any explosion of the
asymptotic variance. In a second part, under stronger confluence assumptions
and with the help of some second order expansions of the asymptotic error, we
go deeper in the study by optimizing the choice of the parameters involved by
the method. In particular, for a given $\varepsilon\textgreater{}0$, we exhibit
some semi-explicit parameters for which the number of iterations of the Euler
scheme required to attain a Mean-Squared Error lower than $\varepsilon^2$ is
about $\varepsilon^{-2}\log(\varepsilon^{-1})$. Finally, we numerically this
Multilevel Langevin estimator on several examples including the simple
one-dimensional Ornstein-Uhlenbeck process but also on a high dimensional
diffusion motivated by a statistical problem. These examples confirm the
theoretical efficiency of the method.
