  The \emph{ab initio} path integral Monte Carlo (PIMC) approach is one of the
most successful methods in quantum many-body theory. A particular strength of
this method is its straightforward access to imaginary-time correlation
functions (ITCF). For example, the well-known density-density ITCF
$F(\mathbf{q},\tau)$ allows one to estimate the linear response of a given
system for all wave vectors $\mathbf{q}$ from a single simulation of the
unperturbed system. Moreover, it constitutes the basis for the reconstruction
of the dynamic structure factor $S(\mathbf{q},\omega)$ -- a key quantity in
state-of-the-art scattering experiments. In this work, we present analogous
relations between the nonlinear density response in quadratic and cubic order
of the perturbation strength and generalized ITCFs measuring correlations
between up to four imaginary-time arguments. As a practical demonstration of
our new approach, we carry out simulations of the warm dense electron gas and
find excellent agreement with previous PIMC results that had been obtained with
substantially larger computational effort. In addition, we give a relation
between a cubic ITCF and the triple dynamic structure factor
$S(\mathbf{q}_1,\omega_1;\mathbf{q}_2,\omega_2)$, which evokes the enticing
possibility to study dynamic three-body effects on an \emph{ab initio} level.
