Temporal clustering of large earthquakes in the Wellington region, New Zealand, has been investigated with a computer model that generates long synthetic seismicity catalogues. The model includes the elastic interactions between faults. Faults included in the model, besides the subduction thrust between the Australian and Pacific plates, are segments of the four major strike-slip faults that overlie the plate interface (Wairarapa, Wellington, Ohariu, and Wairau faults). Parameters of the model are adjusted to reproduce the geologically ohserved slip rates of the strike-slip faults. The seismic slip rate of the subduction thrust, which is unknown, is taken as 25% of the maximum predicted by the plate tectonic convergence rate, and its position fixed according to recent geodetic results. For comparison, the model was rerun with the elastic interactions suppressed, corresponding to the usual approach in the calculation of seismic hazard where each fault is considered in isolation. Considering earthquakes of magnitude 7.2 or more ("characteristic" events in the sense that they rupture most of a fault plane). the number of short (0-3 years) inter-event times is much higher with interactions than for the corresponding case without interactions (46% vs. 2% or all inter-event times). This reduces to 9% vs. 2% if the subduction thrust is removed from the models. Paleoseismic studies of the past seismic behaviour of the subduction thrust are clearly needed if the degree of clustering is to be tightly constrained. Although some other aspects of our model can he improved in future, we think that the probability of significant short-term clustering of large events, normally neglected in hazard studies, is very high. This has important implications for the engineering, insurance and emergency response communities.