Stochastic parameterizations of fast-evolving, subgrid-scale processes are increasingly being used in a range of models from conceptual models to general circulation models. However, stochastic terms are generally included in an ad hoc fashion. In this study, a systematic method-''Hasselmann's method''-of stochastic parameterization is developed through the direct application of rigorously justified limit theorems that predict the effective slow dynamics in systems with coupled slow and fast variables. The multiple Hasselmann models form a hierarchy of models ordered by the time scales over which they are expected to provide good approximations to the slowly evolving variables. Adaptable, efficient algorithms for integrating these reduced models are developed that require minimal changes to the unreduced model.Hasselmann's method is tested on an O(10 000)-dimensional (planetary and synoptic scale) quasigeostrophic model of atmospheric low-frequency variability. Low-dimensional deterministic and stochastic models in the planetary-scale modes alone are derived, which accurately generate the statistics of the corresponding modes of the unreduced model, including the statistical signatures of jet regime behavior. It is shown that deterministic nonlinearity through slow forcing averaged with respect to the fast modes distribution dominates over multiplicative noise in generating the regime behavior.