  Multiple time correlation functions are found in the dynamical description of
different phenomena. They encode and describe the fluctuations of the dynamical
variables of a system. In this paper we formulate a theory of non-Markovian
multiple-time correlation functions (MTCF) for a wide class of systems. We
derive the dynamical equation of the {\it reduced propagator}, an object that
evolve state vectors of the system conditioned to the dynamics of its
environment, which is not necessarily at the vacuum state at the initial time.
Such reduced propagator is the essential piece to obtain multiple-time
correlation functions. An average over the different environmental histories of
the reduced propagator permits us to obtain the evolution equations of the
multiple-time correlation functions. We also study the evolution of MTCF within
the weak coupling limit and it is shown that the multiple-time correlation
function of some observables satisfy the Quantum Regression Theorem (QRT),
whereas other correlations do not. We set the conditions under which the
correlations satisfy the QRT. We illustrate the theory in two different cases;
first, solving an exact model for which the MTCF are explicitly given, and
second, presenting the results of a numerical integration for a system coupled
with a dissipative environment through a non-diagonal interaction.
