Concerned with elliptic operator with stationary random coefficients of integrable correlations and bounded Lipschitz domains, arising from stochastic homogenization theory, this paper is mainly devoted to studying Calderón-Zygmund estimates. As an application, we obtain the homogenization error in the sense of oscillation and fluctuation, and these results are optimal up to a quantity O(ln(1/ε)), which is coming from the quantified sublinearity of correctors in dimension two and less smoothness of the boundary of the domain.The main scheme relies on (weighted) annealed Calderón-Zygmund estimates, which was recently developed by Josien and Otto [27] via a robust non-perturbative argument, independent of the quenched Calderón-Zygmund estimate that originally developed by Armstrong and Daniel [1] and Gloria, Neukamm and Otto [19].Enlightened by Duerinckx and Otto's job [9], the main innovation of the present work is to apply Shen's real arguments (of weighted version) to study elliptic systems with Dirichlet or Neumann boundary conditions. We start from a qualitative description of homogenization theorem in a local way (without using boundary correctors), and find a new form of minimal radius, which proved to be a suitable key to open quantitative homogenization on the boundary value problems, if adopting Gloria-Neukamm-Otto's strategy that originally inspired by Naddaf and Spencer's pioneering work.