Discovery of the KZ spectra was so powerful that it dominated the WT theory for decades thereafter. Such spectra were found for a large variety of physical situations, from quantum to astrophysical applications, and a great effort was put in their numerical and experimental verification. For a long time, studies of spectra dominated WT literature. A detailed account of these works was given in  which is the only book so far written on this subject. Work on these lines has continued till now and KZ spectra were found in new applications, particularly in astrophysics , ocean interior  and even cosmology . However, the spectra do not tell the whole story about the turbulence statistics. In particular, they do not tell us if the wavefield statistics is Gaussian or not and, if not, in what way. This question is of general importance in the field of Turbulence because it is related with the intermittency phenomenon, - an anomalously high probability of large fluctuations. Such ``bursts'' of turbulent wavefields were predicted based on a scaling analysis in  and they were linked to formation of coherent structures, such as whitecaps on sea  or collapses in optical turbulence . To study these problems qualitatively, the kinetic equation description is not sufficient and one has to deal directly with the probability density functions (PDF).
In fact, such a description in terms of the PDF appeared already in the the Peierls 1929 paper simultaneously with the kinetic equation for waves . This result was largely forgotten by the WT community because fine statistical details and intermittency had not interested turbulence researchers until relatively recently and also because, perhaps, this result got in the shade of the KZ spectrum discovery. However, this line of investigation was continued by Brout and Prigogine  who derived an evolution equation for the multi-mode PDF commonly known as the Brout-Prigogine equation. This approach was applied to the study of randomness underlying the WT closures by Zaslavski and Sagdeev . All of these authors, Peierls, Brout and Prigogine and Zaslavski and Sagdeev restricted their consideration to the nonlinear interaction arising from the potential energy only (i.e. the interaction Hamiltonian involves coordinates but not momenta). This restriction leaves out the capillary water waves, Alfven, internal and Rossby waves, as well as many other interesting WT systems. Recently, this restriction was removed by considering the most general three-wave Hamiltonian systems . It was shown that the multi-mode PDF still obeys the Peierls-Brout-Prigogine (PBP) equation in this general case. This work will be described in the present review. We will also present, for the first time, a derivation of the evolution equation for the multi-mode PDF for the general case of four wave-systems. This equation is applicable, for example, to WT of the deep water surface gravity waves and waves in Bose-Einstein condensates or optical media described by the nonlinear Schroedinger (NLS) equation.
We will also describe the analysis of papers  of the randomness assumptions underlying the statistical WT closures. Previous analyses in this field examined validity of the random phase assumption [20,21] without devoting much attention to the amplitude statistics. Such ``asymmetry'' arised from a common mis-conception that the phases evolve much faster than amplitudes in the system of nonlinear dispersive waves and, therefore, the averaging may be made over the phases only ``forgetting'' that the amplitudes are statistical quantities too (see e.g. ). This statement become less obvious if one takes into account that we are talking not about the linear phases but about the phases of the Fourier modes in the interaction representation. Thus, it has to be the nonlinear frequency correction that helps randomising the phases . On the other hand, for three-wave systems the period associated with the nonlinear frequency correction is of the same order in small nonlinearity as the nonlinear evolution time and, therefore, phase randomisation cannot occur faster that the nonlinear evolution of the amplitudes. One could hope that the situation is better for 4-wave systems because the nonlinear frequency correction is still but the nonlinear evolution appears only in the order. However, in order to make the asymptotic analysis consistent, such correction has to be removed from the interaction-representation amplitudes and the remaining phase and amplitude evolutions are, again, at the same time scale (now ). This picture is confirmed by the numerical simulations of the 4-wave systems [26,23] which indicate that the nonlinear phase evolves at the same timescale as the amplitude. Thus, to proceed theoretically one has to start with phases which are already random (or almost random) and hope that this randomness is preserved over the nonlinear evolution time. In most of the previous literature prior to  such preservation was assumed but not proven. Below, we will describe the analysis of the extent to which such an assumption is valid made in .
We will also describe the results of  who derived the time evolution equation for higher-order moments of the Fourier amplitude, and its application to description of statistical wavefields with long correlations and associated ``noisiness'' of the energy spectra characteristic to typical laboratory and numerical experiments. We will also describe the results of  about the time evolution of the one-mode PDF and their consequences for the intermittency of stochastic nonlinear wavefields. In particular, we will discuss the relation between intermittency and the probability fluxes in the amplitude space.