What if a hit song
was a time series?

Part I · The dataset

1.What we are actually modelling

Each track in the dataset carries a popularity score $P_i \in \{0, 1, \dots, 100\}$ that Spotify recomputes daily as a function of recent play counts, recency-weighted. The score is non-stationary at the individual-track level (most tracks decay, a few go viral), but if we aggregate to the population level $\bar{P}_t = \frac{1}{|T_t|} \sum_{i \in T_t} P_i$ where $T_t$ is the set of tracks released in period $t$, we get a smooth, low-dimensional series that summarises how the average track from each cohort scores on the platform's modern popularity algorithm.

Alongside $P_i$, the dataset carries Spotify's 12 audio features per track: danceability, energy, valence (positivity), tempo, loudness, acousticness, instrumentalness, liveness, speechiness, key, mode, and time signature. We treat these as exogenous covariates in the EDA and as the input vector for the audio synth in this section.

Part II · From raw tracks to a time series

2.Aggregation

We collapse the 586,672 tracks into a yearly mean-popularity series$\bar P_t$, plotted below. The line shows the platform's well-documented recency bias: the modern popularity algorithm scores recent releases far higher than 1970s catalogue. By 2019 the average track had a popularity score of $\bar P \approx 44.9$, compared to $\bar P \approx 3.4$in 1950. The 2021 dip is an artefact: the dataset was snapshotted mid-year, and 2021's tracks had less time to accumulate plays.

3.Stationarity

ARIMA needs the series (or some difference of it) to be stationary. We run an Augmented Dickey-Fuller test on the raw level series and on the first-difference $\Delta \bar P_t = \bar P_t - \bar P_{t-1}$. The level series fails ADF (visible trend); the first difference passes at the 5% level, so we set $d = 1$ for the integration order. The auto.arima search confirms the same choice independently.

Part III · ARIMA

4.The model class

An ARIMA($p, d, q$) model expresses the $d$-differenced series as a linear combination of its own lagged values and lagged shocks:

where $L$ is the lag operator and $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$is white noise. ACF and PACF plots of the differenced series suggest two candidate orders: a parsimonious (1, 1, 1) and a slightly richer (2, 1, 1). We let R's auto.arima do the AIC search and confirm the parsimonious model.

5.Forecast and uncertainty bands

The forecast point estimate is the conditional mean $\mathbb{E}[P_{T+h} \mid \mathcal{F}_T]$, and the standard 80% / 95% bands widen as $\sqrt{h}$: each additional period of forecast horizon adds independent noise variance, so the band scales like $\sigma \sqrt{h}$. The studio below lets you slide $(p, d, q)$ manually and watch how the forecast shape changes.

Part IV · GARCH for popularity volatility

6.Why a vol model

The ARIMA residuals show heteroskedastic clustering: large absolute residuals tend to be followed by more large absolute residuals. That violates the constant-variance assumption baked into the standard ARIMA error term, and it is the standard motivation for stacking a conditional-variance model on top. The GARCH(1,1) specification with student-t innovations is the workhorse:

We fit three competing specs: plain sGARCH(1,1) with t innovations, EGARCH(1,1) with normal innovations, and EGARCH(1,1) with t innovations. The EGARCH adds a leverage term that lets negative shocks generate more future variance than positive ones, which is exactly the pattern observed in popularity returns: tracks that suddenly lose popularity are followed by more volatile periods than tracks that gain it. The Student-t EGARCH wins on AIC.

7.Per-genre risk vs reward

Refitting ARIMA per genre, then computing the forecast mean and the return-series standard deviation, gives a clean two-dimensional signature for each genre: forecasted growth on the y-axis, volatility on the x-axis. The desirable quadrant is upper-left: high growth, low volatility. The studio below plots all 14 most-populated genres.

Part V · Conclusion

8.Findings

The same stack quants use for asset returns ports cleanly onto Spotify popularity. The level series is non-stationary by visual inspection and by ADF (p > 0.05); the first difference is stationary; auto.arima converges on ARIMA(1, 1, 1) as the parsimonious overall fit; per-genre models are dominated by simple (1, 0, 1) or (0, 1, 1) structures with a handful of outlier orders. Stacking a Student-t EGARCH(1,1) on the residual returns wins on AIC against both plain GARCH and normal-EGARCH, which says the popularity-return series has both fat tails and leverage — the same pair of features that motivate using these specs on equity returns.

9.What this is and what it is not

This is not a music-industry forecasting service. The dataset truncates at 2021 and the popularity score is a moving target redefined by Spotify's ranking engine. What it is, is a clean demonstration that ARIMA and GARCH are tools, not financial instruments: they apply to any unevenly-clustered, autocorrelated time series with a meaningful first difference. The audio synth in Part I extends the same point in the opposite direction — Spotify's 12 audio features are a feature vector you can do other things with, including procedural sound generation that has nothing to do with forecasting at all.

10.References

1. [1]Box, G. E. P. & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day.
2. [2]Hyndman, R. J. & Khandakar, Y. (2008). Automatic Time Series Forecasting: The forecast Package for R. Journal of Statistical Software, 27(3).
3. [3]Bollerslev, T. (1986). Generalized Autoregressive Conditional Heteroskedasticity. Journal of Econometrics, 31, 307–327.
4. [4]Nelson, D. B. (1991). Conditional Heteroskedasticity in Asset Returns: A New Approach. Econometrica, 59(2), 347–370.
5. [5]Pandya, M. (2021). Spotify Tracks Dataset. Kaggle public release.