The main goal of this work is to provide a probability map for the next moderate to large earthquakes (M ≥ 5.5) in Italy. For this purpose we apply a new nonparametric multivariate model to characterize the spatiotemporal distribution of earthquakes. The method is able to account for tectonics/physics parameters that can potentially influence the spatiotemporal variability and tests their relative importance; moreover, it allows straightforward testing of a variety of hypothesis, such as Seismic Gap, Cluster, and Poisson hypothesis. The method has been applied to Italian seismicity of the last four centuries for earthquakes with M ≥ 5.5. Italy has been divided into 61 irregular zones representing areas with homogeneous tectonic regime resulting from active stress data. Besides the magnitude and the time of the earthquakes, the model includes information on the tectonic stress regime, the homogeneity of its orientation, the number of active faults, the dimension of the area and the homogeneity of the topography. The time distribution of the M ≥ 5.5 earthquakes appears clustered in time for a few years after an event, and then the distribution becomes similar to a memoryless Poisson process, leading to a time-dependent probability map. This map shows that the most likely regions where the next moderate to large earthquakes may occur are Friuli, Umbria-Marche, and part of Southern Apennines and the Calabrian arc.