We develop a large-scale regularity theory of higher order for divergence-form elliptic equations with heterogeneous coefficient fields in the context of stochastic homogenization. Under the assumptions of stationarity and slightly quantified ergodicity of the ensemble, we derive a -"excess decay" estimate on large scales and a -Liouville principle: For a given -harmonic function on a ball , we show that its energy distance to the space of -harmonic "corrected quadratic polynomials" on some ball has the natural decay in the radius above some minimal (random) radius . Our Liouville principle states that the space of -harmonic functions growing at most quadratically has (almost surely) the same dimension as in the constant-coefficient case. The existence of -harmonic "corrected quadratic polynomials" - and therefore our regularity theory - relies on the existence of second-order correctors for the homogenization problem. By an iterative construction, we are able to establish existence of subquadratically growing second-order correctors.