We use direct numerical simulations (DNSs) of the forced Navier–Stokes equation for a three-dimensional incompressible fluid in order to test recent theoretical predictions. We study the two- and three-point spatiotemporal correlation functions of the velocity field in stationary, isotropic, and homogeneous turbulence. We compare our numerical results to the predictions from the Functional Renormalization Group (FRG) which were obtained in the large wavenumber limit. DNSs are performed at various Reynolds numbers and the correlations are analyzed in different time regimes focusing on the large wavenumbers. At small time delays, we find that the two-point correlation function decays as a Gaussian in the variable kt, where k is the wavenumber and t the time delay. We compute a triple correlation from the time-dependent advection-velocity correlations, and find that it also follows a Gaussian decay at small t with the same prefactor as the one of the two-point functions. These behaviors are in precise agreement with the FRG results, and can be simply understood as a consequence of sweeping. At large time delays, the FRG predicts a crossover to an exponential in k2t, which we were not able to resolve in our simulations. However, we analyze the two-point spatiotemporal correlations of the modulus of the velocity and show that they exhibit this crossover from a Gaussian to an exponential decay, although we lack of a theoretical understanding in this case. This intriguing phenomenon calls for further theoretical investigation.