Multiple seasonalities often appear in high-frequency data. In this context multiple seasonal components are usually modelled in a deterministic way by trigonometric functions or dummy variables. This assumption may be too strict. Instead, a more flexible model is to allow the seasonality to slowly change as a seasonal Autoregressive Integrated Moving Average model, where the seasonality is modelled as a stochastic processes. In this study, we propose to model them iteratively, combining different seasonal Autoregressive Integrated Moving Average models. To this end, we test the proposed methodology with Madrid's NO2 hourly measurements of pollutants with daily, weekly and annual seasonalities, due to human activity and weather conditions. Here, we demonstrate the usefulness of our approach by comparing it with other methodological approaches proposed for this type of data. In an extensive exercise involving 15-year hourly forecasts, we show that the proposed procedure performs very well in predicting hourly pollution over a 24-h horizon and improves on alternative procedures. Additionally, the impact on the predictions of covariates such as wind speed, temperature and festivities were evaluated.