  A general method is developed which enables the exact treatment of the
non-Markovian quantum dynamics of open systems through a Monte Carlo simulation
technique. The method is based on a stochastic formulation of the von Neumann
equation of the composite system and employs a pair of product states following
a Markovian random jump process. The performance of the method is illustrated
by means of stochastic simulations of the dynamics of open systems interacting
with a Bosonic reservoir at zero temperature and with a spin bath in the strong
coupling regime.
