We use a novel suite of cosmological hydrodynamical simulations to study the formation of stellar spheroids of Milky Way-mass disc galaxies. The simulations contain accurate treatments of metal-dependent radiative cooling, star formation, supernova feedback, and chemodynamics, and the large volumes that have been simulated yield an unprecedentedly large sample of ∼400 simulated L* disc galaxies. The simulated galaxies are surrounded by low-mass, low-surface brightness stellar haloes that extend out to ∼100 kpc and beyond. The diffuse stellar distributions bear a remarkable resemblance to those observed around the Milky Way, M31 and other nearby galaxies, in terms of mass density, surface brightness, and metallicity profiles. We show that in situ star formation dominates the stellar spheroids (by mass) at radii of ∼30kpc, whereas accretion of stars dominates at larger radii and this change in origin induces a change in slope of the surface brightness and metallicity profiles, which is also present in the observational data. We also show that the in situ mass fraction of the spheroid is linked to the formation epoch of the system. Dynamically older systems have, on average, larger contributions from in situ star formation, although there is significant system-to-system scatter in this relationship. Thus, in situ star formation likely represents the solution to the longstanding failure of pure accretion-based models to reproduce the observed properties of the inner spheroid.