Using a real-time path integral approach we develop an algorithm to calculate multi-time correlation functions of open few-level quantum systems that is applicable to highly nonequilibrium dynamics. The calculational scheme fully keeps the non-Markovian memory introduced by the pure-dephasing type coupling to a continuum of oscillators. Furthermore, we discuss how to deal consistently with the simultaneous presence of non-Markovian and Markovian system reservoir interactions. We apply the method to a crucial test case, namely the evaluation of emission spectra of a laser-driven two-level quantum dot coupled to a continuum of longitudinal acoustic phonons, which give rise to non-Markovian dynamics. Here, we also account for the coupling to a photonic environment, which models radiative decay and can be treated as a Markovian bath. The phonon side bands are found on the correct side of the zero phonon line in our calculation, in contrast to known results where the quantum regression theorem is applied naively to non-Markovian dynamics. Combining our algorithm with a recently improved iteration scheme for performing the required sum over paths we demonstrate the numerical feasibility of our approach to systems with more than two levels. Results are shown for the second-order photonic two-time correlation function of a quantum dot-cavity system with seven states on the Jaynes-Cummings ladder taken into account.