We consider a model problem of diffusion or conductivity in a random medium (for instance, an heterogeneous domain obtained by mixing randomly two different phases, one being the matrix and the other the inclusions). The model that we use takes the form of an uniformly elliptic equation with high oscillatory (at scale S<<1) and random coefficient. Since the 70’ with the work of Kozlov, Papanicolaou and Varadhan, it is well known that the solution, provided the law of the coefficients is stationary and ergodic, can be approximated by its two-scale expansion: that is a first-order expansion in S taking account of the oscillation at scale S. The zero-order term corresponds to the limit as S goes to zero and solves a nice elliptic PDE with deterministic constant coefficients (called the homogenized coefficients) for which we have an explicit formula. However, the homogenized coefficients depend on the so-called first-order corrector which solves an elliptic PDE posed in the whole space and for which numerical computations are out of reach. The goal of this work is to analyse the approximation of the homogenized coefficients by the representative volume element method, that is when we replace the whole space by a large torus of size L. This operation makes the corrector equation easy to solve numerically and provides a natural approximation of the homogenized coefficients by the ones coming from the periodic homogenization theory. In the analysis of the error that we make in this approximation, the variance suffers from two types of error: a random one (that is the fluctuation around its expectation) and a systematic one (that is the difference between the homogenized coefficients and the expectation of the approximation). We focus in this work on the systematic error and we characterize its asymptotic behaviour as L goes to infinity. We show that, in the particular case where the law is generated by a stationary Gaussian field, the asymptotic of the systematic error is characterized by a deterministic matrix depending on the first and second-order correctors, the gradient of the covariance function and a fourth-order tensor involving the whole space Green function of the homogenized elliptic operator.