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 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.PACS numbers: 3.65 Yz, 42.50 Lc Introduction and motivation. Many research contexts are focused on the dynamics of a system (S) that is affected by an environment (B) from which it cannot be considered isolated. Examples of such situations are encountered in statistical physics, condensed matter and quantum optics. We found a concrete example in the description of the dynamics of an atom (S) immersed in an electromagnetic field (B) [1,2].In some circumstances, the analysis of the dynamics of the system is done using the expectation values of its observables over state vectors of the whole system, and then averaging over the environmental degrees of freedom. However, in some other situations, like when studying the response of a system to an external EM field, some additional information is needed. In particular, for the analysis of the spectroscopic properties of a system some multiple-time correlation function (MTCF) has to be computed, usually a two-time correlation function.The dynamics of the system S is usually described through its reduced density operator. Such operator verifies some master equation that in the Markovian case is of Lindblad type [1,3,4,5,6,7]. Complementary to the master equation approach, a series of Monte-Carlo type of approaches based on the so called stochastic Schrödinger equations [1,4,8,9,10] have been developed in the last decade. In such schemes, the dynamics of system state vectors is integrated, and after an average is made over many realizations of environment histories that eventually are understood as a noise and takes into account the environment influence on S. In the non-Markovian case, within the context of nuclear magnetic resonance, the ...