We study the statistical properties of random shock waves in stochastic Burgers equation subject to random space–time perturbations and random initial conditions. By using the response–excitation probability density function (PDF) method and the Mori–Zwanzig (MZ) formulation of irreversible statistical mechanics, we derive exact reduced-order equations for the one-point and two-point PDFs of the solution field. In particular, we compute the statistical properties of random shock waves in the inviscid limit by using an adaptive (shock-capturing) discontinuous Galerkin method in both physical and probability spaces. We consider stochastic flows generated by high-dimensional random initial conditions and random additive forcing terms, yielding multiple interacting shock waves collapsing into clusters and settling down to a similarity state. We also address the question of how random shock waves in space and time manifest themselves in probability space. The proposed new mathematical framework can be applied to different conservation laws, potentially leading to new insights into high-dimensional stochastic dynamical systems and more efficient computational algorithms.