  A minimal requirement for simulating multi-scale systems is to reproduce the
statistical behavior of the slow variables. In particular, a good numerical
method should accurately aproximate the probability density function of the
continuous-time slow variables. In this note we use results from homogenization
and from backward error analysis to quantify how errors of time integrators
affect the mean behavior of trajectories. We show that numerical simulations
converge, not to the exact probability density function (pdf) of the
homogenized multi-scale system, but rather to that of the homogenized modified
equations following from backward error analysis. Using homogenization theory
we find that the observed statistical bias is exacerbated for multi-scale
systems driven by fast chaotic dynamics that decorrelate insufficiently
rapidly. This suggests that to resolve the statistical behavior of trajectories
in certain multi-scale systems solvers of sufficiently high order are
necessary. Alternatively, backward error analysis suggests the form of an
amended vector field that corrects the lowest order bias in Euler's method. The
resulting scheme, a second order Taylor method, avoids any statistical drift
bias. We corroborate our analysis with a numerical example.
