Numerical simulation of high frequency waves in highly heterogeneous media is a challenging problem. Resolving the fine structure of the wave field typically requires extremely small time steps and spatial meshes. We show that capturing macroscopic quantities of the wave field, such as the wave energy density, is achievable with much coarser discretizations. We obtain such a result using a time splitting algorithm that solves separately and successively propagation and scattering in the simplified regime of the parabolic wave equation in a random medium. The mathematical theory of the convergence and statistical properties of the algorithm is based on the analysis of the Wigner transforms in random media. Our results provide a step toward understanding time and space discretizations that are needed in order for the numerical algorithm to capture the correct macroscopic statistics of the wave energy density in a random medium.