  We present a first-principles scheme for incorporating many-body interactions
into the unified description of the quadratic optical response to light of
noncentrosymmetric crystals. The proposed method is based on time-dependent
current-density response theory and includes the electron-hole attraction
\textit{via} a tensorial long-range exchange-correlation kernel, which we
calculate self-consistently using the bootstrap method. By bridging with the
Wannier-interpolation of the independent-particle transition matrix elements,
the resulting numerical scheme is very general and allows resolving narrow
many-body spectral features at low computational cost. We showcase its
potential by inspecting the second-harmonic generation in the benchmark
zinc-blende semiconductor GaAs, the layered graphitic semiconductor BC$_{2}$N
and the Weyl semimetal TaAs. Our results show that excitonic effects can give
rise to large and sharply localized one- and two-photon resonances that are
absent in the independent-particle approximation. We find overall good
agreement with available experimental measurements, capturing the magnitude and
peak-structure of the spectrum as well as the angular dependence at fixed
photon energy. The implementation of the method in Wannier-based code packages
can serve as a basis for performing accurate theoretical predictions of
quadratic optical properties in a vast pool of materials.
