Horizontal density layers are commonly observed in stratified turbulence. Recent work (e.g. Taylor & Zhou, J. Fluid Mech., vol. 823, 2017, R5) has reinvigorated interest in the Phillips instability (PI), by which density layers form via negative diffusion if the turbulent buoyancy flux weakens as stratification increases. Theoretical understanding of PI is incomplete, in part because it remains unclear whether and by what mechanism the flux-gradient relationship for a given example of turbulence has the required negative-diffusion property. Furthermore, the difficulty of analysing the flux-gradient relation in evolving turbulence obscures the operating mechanism when layering is observed. These considerations motivate the search for an example of PI that can be analysed clearly. Here PI is shown to occur in two-dimensional Boussinesq sheared stratified turbulence maintained by stochastic excitation. PI is analysed using the second-order S3T closure of statistical state dynamics, in which the dynamics is written directly for statistical variables of the turbulence. The predictions of S3T are verified using nonlinear simulations. This analysis provides theoretical underpinning of PI based on the fundamental equations of motion that complements previous analyses based on phenomenological models of turbulence.