  This manuscript continues and extends in various directions the result in
arXiv:2104.11204, which gave a full derivation of the wave kinetic equation
(WKE) from the nonlinear Schr\"{o}dinger (NLS) equation in dimensions $d\geq
3$. The wave kinetic equation describes the effective dynamics of the second
moments of the Fourier modes of the NLS solution at the kinetic timescale, and
in the kinetic limit in which the size of the system diverges to infinity and
the strength of the nonlinearity vanishes asymptotically according to a
specified scaling law.
  Here, we investigate the behavior of the joint distribution of these Fourier
modes and derive their effective limit dynamics at the kinetic timescale. In
particular, we prove propagation of chaos in the wave setting: initially
independent Fourier modes retain this independence in the kinetic limit. Such
statements are central to the formal derivations of all kinetic theories,
dating back to the work of Boltzmann (Stosszahlansatz). We obtain this by
deriving the asymptotics of the higher Fourier moments, which are given by
solutions of the wave kinetic heirarchy (WKH) with factorized initial data. As
a byproduct, we also provide a rigorous justification of this hierarchy for
general (not necessarily factorized) initial data.
  We treat both Gaussian and non-Gaussian initial distributions. In the
Gaussian setting, we prove propagation of Gaussianity as we show that the
asymptotic distribution retains the Gaussianity of the initial data in the
limit. In the non-Gaussian setting, we derive the limiting equations for the
higher order moments, as well as for the density function (PDF) of the
solution. Some of the results we prove were conjectured in the physics
literature, others appear to be new. This gives a complete description of the
statistics of the solutions in the kinetic limit.
