We derive the evolution equations for two-time correlation functions of a generalized non-Markovian open quantum system based on a modified stochastic Schrödinger equation approach. We find that the two-time reduced propagator, an object that used to be characterized by two independent stochastic processes in the Hilbert space of the system, can be simplified and obtained by taking ensemble average over one single noise. This discovery can save the cost of computation, and accelerate the converging process when taking the average over noisy trajectories. As a result, our method can be widely applied to many open quantum models, especially large-scale systems and extend the quantum regression theory to the non-Markovian case. In the short-time simulations, it is observed a significant difference between Markovian and non-Markovian cases, which can be applied to realize the environmental spectrum detection and enhance the measurement sensitivity in varying open quantum systems.