The study focuses on the estimation and modeling of the temporal rates of major explosions and paroxysms at Stromboli volcano (also named small-scale and large-scale paroxysms respectively). The analysis was further motivated by the paroxysm of July 3rd 2019, which raised, once again, the attention of the scientific community and civil protection authorities on the volcanic hazards of Stromboli. In fact, at the present state of knowledge, major explosions and paroxysms cannot be forecasted based on monitoring data, and a full probabilistic assessment based on past eruption data would be quite useful for scientific and civil protection purposes. In the study we perform a time series analysis either considering the last ~150 years of reconstructed activity and the most recent 35 years. We included the estimation of event rates and rate changes in time. Results clearly highlight that the activity is non-homogeneus in time, with a significant low of activity between about 1960 and 1990. Maximum values of event rates were computed during the first half of last century, for both major explosions and paroxysms, whereas the rate of paroxysms is significantly lower in the last decades with respect to maximum rates. We also accomplish a statistical analysis of the inter-event times, enabling us to determine if the data can be modeled as a Poisson process or not, e.g. if it shows time dependent distributions, recurring cycles, or temporal clusters. The uncertainty quantification on the current and future rates is mainly related to the choice of the modeling assumptions. The study represents a crucial progress towards quantitative hazard and risk assessments at Stromboli, which is particularly relevant for the thousands of people (e.g. tourists, guides and volcanologists) that regularly climb the volcano every year.