Accurate modeling of irradiation damage processes is important to predict materials performance in nuclear environments. The mean-field rate theory (RT) approach has been and remains the standard method for calculations of the kinetics of damage accumulation under irradiation. Despite its many advantages, when using RT, very large numbers of ordinary differential equations (ODEs) need to be solved if one is interested in the kinetics of complex defect populations containing defect clusters made up of constituent defect species of different types, e.g. He, H, O, alloying impurities, etc. Here we present a stochastic variant of RT, which we term stochastic cluster dynamics (SCD), intended as an alternative to the standard ODE-based implementation. SCD obviates the need to solve the exceedingly large sets of ODEs and relies instead on sparse stochastic sampling from the underlying kinetic Master Equation. The new method evolves an integer-valued defect population in a finite material volume. In this paper we describe the general ideas behind the SCD method, give essential details of our numerical implementation and present SCD simulations that verify the new method against standard RT implementations on two well characterized material models. We then apply the method to simulate triple beam irradiation experiments of heavy Fe ions, He, and H on model ferritic alloys up to a dose of 1dpa.