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.