  We extend the Nakajima-Zwanzig projection operator technique to the
determination of multitime correlation functions of open quantum systems. The
correlation functions are expressed in terms of certain multitime homogeneous
and inhomogeneous memory kernels for which suitable equations of motion are
derived. We show that under the condition of finite memory times these
equations can be used to determine the memory kernels by employing an exact
stochastic unraveling of the full system-environment dynamics. The approach
thus allows to combine exact stochastic methods, feasible for short times, with
long-time master equation simulations. The applicability of the method is
demonstrated by numerical simulations of 2D-spectra for a donor-acceptor model,
and by comparison of the results with those obtained from the reduced hierarchy
equations of motion. We further show that the formalism is also applicable to
the time evolution of a periodically driven two-level system initially in
equilibrium with its environment.
