We present a two-stage method to simulate the ground motions produced by an earthquake by using stochastic summation of small earthquakes. In this method, identical small earthquakes are multiplied by a scaling factor and summed together with time delays randomly distributed, during the two stages, over the source duration. The summation scheme is characterized by four fundamental parameters: the number of summed small earthquakes, the scaling factor, and both probability densities of time delays used in the first and second stages. By a proper choice of these parameters, this method generates a large number of synthetic time histories that, on average, agree exactly with the ω−2 model in the whole frequency band. The produced time histories are sufficiently realistic and different from each other to be associated with a multitude of rupture processes that could happen during an earthquake. However, because the extended target fault is approximated by a point source, this method does not take into account possible directivity effects and is not appropriate to simulate ground motions for near-source sites. We test this method on the Oaxaca earthquake (1999, Mw 7.5, Mexico) at regional distances and on the two mainshocks of the Umbria Marche crisis (1997, Mw 5.7 and Mw 6.0, Italy) at local distances. We found that the simulated ground motions fit the observed data well, both in time and in frequency domains. Within simulation context, only specification of seismic moment and stress drop is required for the target event. Because the magnitude and then the seismic moment are necessarily specified, the stress drop plays a major role in ground-motion simulation.