We perform detailed comparison of the semianalytic halo model predictions with measurements in numerical simulations of the two- and three-point correlation functions, as well as power spectrum and bispectrum. We discuss the accuracy and self-consistency of the halo model description of gravitational clustering in the nonlinear regime and constrain halo model parameters. We exploit the recently proposed multipole expansion of three-point statistics that expresses rotation invariance in the most natural way. This not only offers technical advantages by reducing the integrals required for the halo model predictions, but also amounts to a convenient way of compressing the information contained in the three-point correlation functions (3PCFs). We find that, with an appropriate choice of the halo boundary and mass function cutoff, halo model predictions are in good agreement with the bispectrum measured in numerical simulations. However, the halo model predicts less than the observed configuration dependence of the 3PCF on ~megaparsec scales. This effect is mainly due to quadrupole moment deficit, possibly related to the assumption of spherical halo geometry. Our analysis shows that using its harmonic decomposition, the full configuration dependence of the 3PCF in the nonlinear regime can be compressed into just a few numbers, the lowest multipoles. Moreover, these multipoles are closely related to the highest signal-to-noise eigenmodes of the 3PCF. Therefore this estimator may simplify future analyses aimed at constraining cosmological and halo model parameters from observational data.