Simulations of strongly stratified turbulence often exhibit coherent large-scale structures called vertically sheared horizontal flows (VSHFs). VSHFs emerge in both twodimensional (2D) and three-dimensional (3D) stratified turbulence with similar vertical structure. The mechanism responsible for VSHF formation is not fully understood. In this work, the formation and equilibration of VSHFs in a 2D Boussinesq model of stratified turbulence is studied using statistical state dynamics (SSD). In SSD, equations of motion are expressed directly in the statistical variables of the turbulent state. Restriction to 2D turbulence makes available an analytically and computationally attractive implementation of SSD referred to as S3T, in which the SSD is expressed by coupling the equation for the horizontal mean structure with the equation for the ensemble mean perturbation covariance. This second order SSD produces accurate statistics, through second order, when compared with fully nonlinear simulations. In particular, S3T captures the spontaneous emergence of the VSHF and associated density layers seen in simulations of turbulence maintained by homogeneous large-scale stochastic excitation. An advantage of the S3T system is that the VSHF formation mechanism, which is wavemean flow interaction between the emergent VSHF and the stochastically excited largescale gravity waves, is analytically understood in the S3T system. Comparison with fully nonlinear simulations verifies that S3T solutions accurately predict the scale selection, dependence on stochastic excitation strength, and nonlinear equilibrium structure of the VSHF. These results facilitate relating VSHF theory and geophysical examples of turbulent jets such as the ocean's equatorial deep jets.