A combined method was developed to forecast the spatial and the temporal probability of occurrence of rainfall-induced shallow landslides over large areas. The method also allowed to estimate the dynamic change of this probability during a rainfall event. The model, developed through a data-driven approach basing on Multivariate Adaptive Regression Splines technique, was based on a joint probability between the spatial probability of occurrence (susceptibility) and the temporal one. The former was estimated on the basis of geological, geomorphological, and hydrological predictors. The latter was assessed considering short-term cumulative rainfall, antecedent rainfall, soil hydrological conditions, expressed as soil saturation degree, and bedrock geology. The predictive capability of the methodology was tested for past triggering events of shallow landslides occurred in representative catchments of Oltrepò Pavese, in northern Italian Apennines. The method provided excellently to outstanding performance for both the really unstable hillslopes (area under ROC curve until 0.92, true positives until 98.8%, true negatives higher than 80%) and the identification of the triggering time (area under ROC curve of 0.98, true positives of 96.2%, true negatives of 94.6%). The developed methodology allowed us to obtain feasible results using satellite-based rainfall products and data acquired by field rain gauges. Advantages and weak points of the method, in comparison also with traditional approaches for the forecast of shallow landslides, were also provided.