In this paper, we study high order correctors in stochastic homogenization. We consider elliptic equations in divergence form on $$\mathbb {Z}^d$$ , with the random coefficients constructed from i.i.d. random variables. We prove moment bounds on the high order correctors and their gradients under dimensional constraints. It implies the existence of stationary correctors and stationary gradients in high dimensions. As an application, we prove a two-scale expansion of the solutions to the random PDE, which identifies the first and higher order random fluctuations in a strong sense.