We have developed a Monte Carlo methodology for the estimation of seismic hazard at a site or across an area. This method uses a multitudinous resampling of an earthquake catalog, perhaps supplemented by parametric models, to construct synthetic earthquake catalogs and then to find earthquake ground motions from which the hazard values are found. Large earthquakes extrapolated from a Gutenberg-Richter recurrence relation and characteristic earthquakes can be included in the analysis. For the ground motion attenuation with distance, the method can use either a set of observed ground motion observations from which estimates are randomly selected, a table of ground motion values as a function of epicentral distance and magnitude, or a parametric ground motion attenuation relation. The method has been tested for sites in New England using an earthquake catalog for the northeastern United States and southeastern Canada, and it yields reasonable ground motions at standard seismic hazard values. This is true both when published ground motion attenuation relations and when a dataset of observed peak acceleration observations are used to compute the ground motion attenuation with distance. The hazard values depend to some extent on the duration of the synthetic catalog and the specific ground motion attenuation used, and the uncertainty in the ground motions increases with decreasing hazard probability. The program gives peak accelerations that are comparable to those of the 1996 U.S. national seismic hazard maps. The method can be adapted to compute seismic hazard for cases where there are temporal or spatial variations in earthquake occurrence rates or source parameters.