  Within the framework of probability distributions on projective Hilbert space
a scheme for the calculation of multitime correlation functions is developed.
The starting point is the Markovian stochastic wave function description of an
open quantum system coupled to an environment consisting of an ensemble of
harmonic oscillators in arbitrary pure or mixed states. It is shown that matrix
elements of reduced Heisenberg picture operators and general time-ordered
correlation functions can be expressed by time-symmetric expectation values of
extended operators in a doubled Hilbert space. This representation allows the
construction of a stochastic process in the doubled Hilbert space which enables
the determination of arbitrary matrix elements and correlation functions. The
numerical efficiency of the resulting stochastic simulation algorithm is
investigated and compared with an alternative Monte Carlo wave function method
proposed first by Dalibard et al. [Phys. Rev. Lett. {\bf 68}, 580 (1992)]. By
means of a standard example the suggested algorithm is shown to be more
efficient numerically and to converge faster. Finally, some specific examples
from quantum optics are presented in order to illustrate the proposed method,
such as the coupling of a system to a vacuum, a squeezed vacuum within a finite
solid angle, and a thermal mixture of coherent states.
