The aim of this paper is to construct a nonequilibrium statistical- mechanics theory to study hysteresis in ferromagnetic systems. We study the hysteretic response of model spin systems to periodic magnetic fields H(t) as a function of the amplitude H0 and frequency Ω. At fixed H 0, we find conventional, squarelike hysteresis loops at low Ω, and rounded, roughly elliptical loops at high Ω, in agreement with experiments. For the O(N→∞), d=3, (Φ2)2 model with Langevin dynamics, we find a novel scaling behavior for the area A of the hysteresis loop, of the form A∝H0.660Ω 0.33. We carry out a Monte Carlo simulation of the hysteretic response of the two-dimensional, nearest-neighbor, ferromagnetic Ising model. These results agree qualitatively with the results obtained for the O(N) model.