We consider a discrete elliptic equation with random coefficients $A$, which (to fix ideas) are identically distributed and independent from grid point to grid point $x\in\mathbb{Z}^d$. On scales large w.\ r.\ t.\ the grid size (i.\ e.\ unity), the solution operator is known to behave like the solution operator of a (continuous) elliptic equation with constant deterministic coefficients. These symmetric ''homogenized'' coefficients $A_{hom}$ are characterized by % $$ \xi\cdot A_{hom}\xi\;=\;\langle\left((\xi+\nabla\phi)\cdot A(\xi+\nabla\phi)\right)(0)\rangle, \quad\xi\in\mathbb{R}^d, $$ % where the random field $\phi$ is the unique stationary solution of the ''corrector problem'' % $$ -\nabla\cdot A(\xi+\nabla\phi)\;=\;0 $$ % and $\langle\cdot\rangle$ denotes the ensemble average. \medskip It is known (''by ergodicity'') that the above ensemble average of the energy density $e=(\xi+\nabla\phi)\cdot A(\xi+\nabla\phi)$, which is a stationary random field, can be recovered by a system average. We quantify this by proving that the variance of a spatial average of $e$ on length scales $L$ is estimated as follows: % $$ {\rm var}\left[\sum_{x\in\mathbb{Z}^d}\eta_L(x)\,e(x)\right] \;\lesssim\;L^{-d}, $$ % where the averaging function (i.\ e.\ $\sum_{x\in\mathbb{Z}^d}\eta_L(x)=1$, ${\rm supp}\eta_L\subset[-L,L]^d$) has to be smooth in the sense that $|\nabla\eta_L|\lesssim L^{-1}$. In two space dimensions (i.\ e.\ $d=2$), there is a logarithmic correction. \medskip In other words, smooth averages of the energy density $e$ behave like as if $e$ would be independent from grid point to grid point (which it is not for $d>1$). This result is of practical significance, since it allows to estimate the error when numerically computing $A_{hom}$.