  Using a diagrammatic approach to Eulerian perturbation theory, we
analytically calculate the variance and skewness of the density and velocity
divergence induced by gravitational evolution from Gaussian initial conditions,
including corrections *beyond* leading order. Except for the power spectrum,
previous calculations in cosmological perturbation theory have been confined to
leading order (tree level)-we extend these to include loop corrections. For
scale-free initial power spectra, the one-loop variance \sigma^2 = \sigma^2_l +
1.82 \sigma^4_l and the skewness S_3 = 34/7 + 9.8 \sigma^2_l, where \sigma_l is
the rms fluctuation of the linear density field. We also compute loop
corrections to the variance, skewness, and kurtosis for several non-linear
approximation schemes, where the calculation can be easily generalized to
1-point cumulants of higher order and arbitrary number of loops. We find that
the Zel'dovich approximation gives the best approximation to the loop
corrections of exact perturbation theory, followed by the Linear Potential
approximation (LPA) and the Frozen Flow approximation (FFA), in qualitative
agreement with the relative behavior of tree-level results. In LPA and FFA,
loop corrections are infrared divergent for spectral indices n < 0; this is
related to the breaking of Galilean invariance in these schemes.
