In this talk, I present a closure theory of quasigeostrophic turbulence for arbitrary shear flows based on stochastic models. A stochastic model represents the eddy-mean flow interaction through a nonnormal dynamical operator, and parameterizes the eddy-eddy interactions by an effective dissipation and random excitation. In the context of a stochastic model, the main facts which a closure theory of turbulence must explain are (1) the form of the eddy dissipation operator, (2) the magnitude of the eddy dissipation, (3) the space-lag correlation of the random forcing at every point in space, (4) the rate of enstrophy transfer to subgrid-scales. We present a theory, called stochastic similarity theory, which accounts for all of these. The theory gives good quantitative agreement for the structure of eddy fluxes and variances obtained from fully nonlinear simulations, and predicts a -3 slope for the equilibrium energy spectrum. The theory is based on the hypothesis that the eddy-eddy nonlinear terms depend only on the local properties of the eddy statistics. Under a linearity assumption, this theory can be solved for all of the unknown parameters in the stochastic model except one, namely (4) listed above. However, the stochastic model predicts a -3 energy spectrum at asymptotically large wavenumbers, consistent with inertial range theory. By invoking inertial range theory for the functional dependence of the energy spectrum with enstrophy cascade rate, the stochastic model can be closed completely. This theory, which avoids phenomenological mixing-length arguments, can be applied to arbitrary shear-flows and gives reasonable statistical predictions of midlatitude eddy fluxes and variances.