A general theoretical framework for describing the thermally induced sequential decay in atomic clusters is presented. The scheme relies on a full treatment of individual dissociation steps based on phase space theory (PST), built into a kinetic Monte Carlo (kMC) procedure. This combined PST/kMC approach allows one to follow the evolution of several statistical properties such as the size, the angular momentum, or the temperature of the cluster over arbitrarily long time scales. Quantitative accuracy is achieved by incorporating anharmonicities of the vibrational densities of states, the rigorous conservation of angular momentum via the effective dissociation potential, and a proper calibration of the rate constants. The approach is tested and validated on selected Lennard-Jones clusters in various situations. Several approximations, including a mean-field rate equation treatment, are critically discussed; possible extensions are presented.