We introduce an exact open system method to describe the dynamics of quantum systems that are strongly coupled to specific types of environments comprising of spins, such as central spin systems. Our theory is similar to the established non-Markovian quantum state diffusion theory, but for a spin bath instead of a Gaussian bath. The method allows us to represent the time-evolved reduced state of the system as an ensemble average of stochastically evolving pure states. We present a comprehensive theory for arbitrary linear spin environments at both zero and finite temperatures. Furthermore, we introduce a hierarchical expansion method that enables the numerical computation of the time evolution of the stochastic pure states, facilitating a numerical solution of the open system problem in relevant strong coupling regimes.