We study the impact of setting initial conditions in numerical simulations using the standard procedure based on the Zel'dovich approximation (ZA). As it is well known from the perturbation theory, ZA initial conditions have incorrect second- and higher-order growth and therefore excite long-lived transients in the evolution of the statistical properties of density and velocity fields. We also study the improvement brought by using more accurate initial conditions based on second-order Lagrangian perturbation theory (2LPT). We show that 2LPT initial conditions reduce transients significantly and thus are much more appropriate for numerical simulations devoted to precision cosmology. Using controlled numerical experiments with ZA and 2LPT initial conditions, we show that simulations started at redshift zi= 49 using the ZA underestimate the power spectrum in the non-linear regime by about 2, 4 and 8 per cent at z= 0, 1, and 3, respectively, whereas the mass function of dark matter haloes is underestimated by 5 per cent at m= 1015 M⊙ h−1 (z= 0) and 10 per cent at m= 2 × 1014 M⊙ h−1 (z= 1). The clustering of haloes is also affected to the few per cent level at z= 0. These systematics effects are typically larger than statistical uncertainties in recent mass function and power spectrum fitting formulae extracted from numerical simulations. At large scales, the measured transients in higher-order correlations can be understood from first principle calculations based on perturbation theory.