Reference: Time Series Anatomy

A reference for the structural pieces every time-series model is built out of: trend, seasonality, residual, autocorrelation, and the stationarity machinery that ties them together. Use this as a lookup table when you’re trying to decide whether to difference, model the trend explicitly, or just feed the raw series to a tree.

Roster

Signal / Quantity	One-line definition	Range / type	What it tells you	Where it’s used in the corpus
Trend $T_t$	Smooth, long-term drift in the level of the series	$\mathbb{R}$; often modeled as a polynomial or piecewise-linear path	Which direction the series is moving once you ignore the wiggles	ARIMA and SARIMA, Prophet vs. statistical models)
Seasonality $S_t$	Calendar-driven pattern that repeats at a fixed period (12 months, 52 weeks, 24 hours, etc.)	$\mathbb{R}$; sums to ~0 over one period in the additive case	Day-of-week or month-of-year regularity you can plan around	Why your forecast fails on Christmas)
Residual $R_t$	What’s left after trend and seasonality are removed: irregular noise, events, regime changes	$\mathbb{R}$; ideally behaves like white noise	Whether you still have unmodeled structure or outliers	How to detect and handle data drift)
Autocorrelation $\rho_k$	Correlation of the series with itself lagged by $k$ steps	$[-1, 1]$	How long the “memory” of the series is	Feature engineering for time series)
Partial autocorrelation $\phi_{kk}$	Direct correlation at lag $k$, after removing the effect of intermediate lags	$[-1, 1]$	Which lag to cut off at when picking an AR order	ARIMA and SARIMA
Stationarity (weak)	Mean and variance don’t change over time; autocovariance depends only on the lag, not the time	pass/fail test	Whether you can fit ARMA-style models directly on the raw series	Classical forecasting vs ML)
Differencing $\nabla y_t$	$y_t - y_{t-1}$ (first difference); sometimes $\nabla^2$	$\mathbb{R}$	The standard trick to remove a unit root / linear trend	Time series cross-validation)
ADF statistic	Augmented Dickey-Fuller test statistic for the null “there is a unit root”	$\mathbb{R}$, more negative = stronger rejection	Whether to reject non-stationarity at all	ARIMA and SARIMA

The three-signal mental model

Almost every classical time-series decomposition writes the observed series $y_t$ as

$y_t = T_t + S_t + R_t \quad \text{(additive)}$ $y_t = T_t \cdot S_t \cdot R_t \quad \text{(multiplicative)}$

The trend $T_t$ is the slow, smooth drift. The seasonality $S_t$ is a fixed-period oscillation — the up/down that comes back every 12 months because of the calendar, not because of the trend. The residual $R_t$ is everything else: noise, holiday spikes, news events, regime changes. The whole game of forecasting is to extract $T_t$ and $S_t$ well enough that $R_t$ is small and approximately white.

You usually want the additive form unless the size of the seasonal swing grows with the level of the series. Retail sales with seasonal multiplicative spikes (small in the early years, large once revenue grows) are the classic case where multiplicative decomposition — or equivalently a log transform followed by additive decomposition — wins.

The additive decomposition can be written as a regression onto a set of basis functions:

$y_t = \beta_0 + \beta_1 t + \sum_{j=1}^{K} \left[ \alpha_j \cos\!\left(\tfrac{2\pi j t}{m}\right) + \gamma_j \sin\!\left(\tfrac{2\pi j t}{m}\right) \right] + \varepsilon_t$

where $m$ is the seasonal period (e.g. 12 for monthly data) and $K \approx m/2$ harmonics capture the seasonal shape. STL replaces the polynomial trend with a locally weighted scatterplot smoother (LOESS) and uses a moving-average based seasonal smoother; Prophet replaces the trend with a piecewise linear (or saturating logistic) curve and the seasonal term with a Fourier series. They all produce the same three pieces — they just regularize them differently.

Plain English	Statistical symbol	Python equivalent
Observed value at time $t$	$y_t$	series.iloc[t]
Trend component	$T_t$	decomp.trend.iloc[t]
Seasonal component	$S_t$	decomp.seasonal.iloc[t]
Residual / remainder	$R_t$ or $\varepsilon_t$	decomp.resid.iloc[t]
Seasonal period	$m$	period=12
Additive identity	$y_t = T_t + S_t + R_t$	model="additive"
Multiplicative identity	$y_t = T_t \cdot S_t \cdot R_t$	model="multiplicative" or log-transform + additive

Autocorrelation, ACF, PACF, stationarity, differencing

Autocorrelation at lag $k$ is just the Pearson correlation between $y_t$ and $y_{t-k}$, computed over the overlapping part of the series. The ACF plots $\rho_k$ for $k = 0, 1, 2, \ldots$ and answers “how much does today still tell me about $k$ days from now?” The PACF removes the indirect effect that flows through intermediate lags: if $y_t$ is correlated with $y_{t-1}$ which is itself correlated with $y_{t-2}$, the PACF at lag 2 reports only the direct correlation between $y_t$ and $y_{t-2}$.

A series is stationary (in the weak sense) when its mean, variance, and autocovariance structure don’t change over time. Linear ARMA-style models assume this; if you feed them a series with a drifting mean, the model will happily fit and then forecast badly. The standard fix is differencing: replace $y_t$ with $\nabla y_t = y_t - y_{t-1}$, which kills a linear trend and any unit root, and re-test. If one difference isn’t enough, difference again (this is rare; twice usually means you’ve mis-specified something).

The autocovariance at lag $k$ is

$\gamma(k) = \mathbb{E}\!\left[(y_t - \mu)(y_{t+k} - \mu)\right]$

and the autocorrelation is the normalized version

$\rho_k = \frac{\gamma(k)}{\gamma(0)}.$

The sample ACF uses

$\hat{\rho}_k = \frac{\sum_{t=1}^{n-k}(y_t - \bar y)(y_{t+k} - \bar y)}{\sum_{t=1}^{n}(y_t - \bar y)^2}.$

The PACF at lag $k$ is the last coefficient $\phi_{kk}$ in the AR($k$) regression

$y_t = \phi_{k1} y_{t-1} + \phi_{k2} y_{t-2} + \cdots + \phi_{kk} y_{t-k} + \varepsilon_t.$

The first difference is $\nabla y_t = (1 - B) y_t = y_t - y_{t-1}$ where $B$ is the backshift operator; the $d$-th difference is $(1 - B)^d y_t$. A seasonal difference of period $m$ is $\nabla_m y_t = y_t - y_{t-m}$.

Plain English	Statistical symbol	Python equivalent
Lag-$k$ autocorrelation	$\rho_k$	acf(series, nlags=k)[k]
Lag-$k$ partial autocorrelation	$\phi_{kk}$	pacf(series, nlags=k, method="ywm")[k]
First difference	$\nabla y_t$	series.diff().dropna()
Seasonal difference (period $m$)	$\nabla_m y_t$	series.diff(m).dropna()
ADF null hypothesis	$H_0$: unit root present (non-stationary)	adfuller(series)[1] > 0.05 → don’t reject
Backshift	$B y_t = y_{t-1}$	series.shift(1)

Reading the plots

• ACF tails off slowly, PACF cuts off after lag $p$ → an AR($p$) signature; consider $p$ AR terms.
• PACF tails off slowly, ACF cuts off after lag $q$ → an MA($q$) signature; consider $q$ MA terms.
• Both tail off slowly → mixed ARMA; consider both, or you may have under-differenced.
• ACF has significant spikes at lags $m, 2m, 3m, \ldots$ → seasonal ARMA territory; you need a seasonal difference.

Worked Python example: three decompositions of the same series

Build a synthetic monthly series with a linear trend, a yearly seasonal cycle, and Gaussian noise. Then decompose it three ways: classical seasonal_decompose, STL, and Prophet. Compare what each method hands back.

import numpy as np
import pandas as pd

np.random.seed(7)
n = 60  # 5 years of monthly data
t = np.arange(n)
trend = 0.15 * t + 10
seasonality = 4 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 1.5, n)
y = trend + seasonality + noise

idx = pd.date_range("2019-01-01", periods=n, freq="MS")
series = pd.Series(y, index=idx, name="y")
print(series.head().round(3))
# 2019-01-01    12.536
# 2019-02-01    11.451
# 2019-03-01    13.813
# 2019-04-01    15.061
# 2019-05-01    12.881

Classical additive decomposition (moving-average based)

from statsmodels.tsa.seasonal import seasonal_decompose

cl = seasonal_decompose(series, model="additive", period=12)
classical_df = pd.DataFrame({
    "trend": cl.trend,
    "seasonal": cl.seasonal,
    "resid": cl.resid,
})
print(classical_df.head(14).round(3))
#               trend  seasonal  resid
# 2019-01-01      NaN    -0.108    NaN
# 2019-02-01      NaN     1.759    NaN
# 2019-03-01      NaN     4.666    NaN
# 2019-04-01      NaN     2.668    NaN
# 2019-05-01      NaN     4.518    NaN
# 2019-06-01      NaN     1.987    NaN
# 2019-07-01   10.819     0.171 -0.091
# 2019-08-01   10.908    -2.954 -1.536
# 2019-09-01   11.053    -3.595  1.804
# 2019-10-01   11.070    -4.091  1.272
# 2019-11-01   11.187    -4.384  0.295
# 2019-12-01   11.429    -0.637 -1.400
# 2020-01-01   11.604    -0.108  1.062
# 2020-02-01   11.785     1.759  0.013

The first and last period/2 rows of the trend are NaN because the classical method uses a centered moving average of width 12 (so it needs 6 observations on each side). The seasonal column is forced to sum to zero over one period and is exactly periodic. The residual is $y - T - S$.

STL (LOESS-based, robust)

from statsmodels.tsa.seasonal import STL

stl = STL(series, period=12, robust=True).fit()
stl_df = pd.DataFrame({
    "trend": stl.trend,
    "seasonal": stl.seasonal,
    "resid": stl.resid,
    "robust_weight": stl.weights,
})
print(stl_df.head(14).round(3))
#               trend  seasonal  resid  robust_weight
# 2019-01-01  10.194     2.079  0.263          0.978
# 2019-02-01  10.305     1.354 -0.208          0.984
# 2019-03-01  10.418     2.938  0.457          0.934
# 2019-04-01  10.533     4.059  0.470          0.894
# 2019-05-01  10.648     3.008 -0.775          0.802
# 2019-06-01  10.766     2.214 -0.227          0.981
# 2019-07-01  10.883     0.233 -0.218          0.982
# 2019-08-01  11.002    -4.479 -0.105          0.996
# 2019-09-01  11.121    -2.013  0.154          0.992
# 2019-10-01  11.242    -3.023  0.031          1.000
# 2019-11-01  11.366    -4.017 -0.251          0.977
# 2019-12-01  11.498    -2.252  0.146          0.993
# 2020-01-01  11.646     1.313 -0.401          0.948
# 2020-02-01  11.790     1.212  0.556          0.893

STL gives you a trend value at every observation (no NaN edges), a smooth and slowly-varying seasonal pattern that is not constrained to be perfectly periodic, and a robustness weight per observation that you can use to flag outliers — in this run, weights range from 0.0 to 1.0 and 31 of the 60 observations fall below 0.95, so the reweighting is doing real work, not sitting at a constant 1.0. The price you pay: STL has more knobs (trend and seasonal, the lengths of the trend and seasonal LOESS smoothers) and the components aren’t as cleanly separable as in the classical case.

Prophet (piecewise-linear trend + Fourier seasonality)

from prophet import Prophet

df = series.reset_index()
df.columns = ["ds", "y"]
m = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    seasonality_mode="additive",
)
m.fit(df)
forecast = m.predict(df)

prophet_df = forecast[["ds", "trend", "yearly", "additive_terms"]]
prophet_df = prophet_df.rename(columns={"yearly": "seasonal"})
print(prophet_df.head(14).round(3))
#            ds   trend  seasonal  additive_terms
# 0  2019-01-01   9.706     0.673           0.673
# 1  2019-02-01   9.862     1.645           1.645
# 2  2019-03-01  10.003     4.837           4.837
# 3  2019-04-01  10.159     4.061           4.061
# 4  2019-05-01  10.311     3.275           3.275
# 5  2019-06-01  10.467     2.219           2.219
# 6  2019-07-01  10.618     0.384           0.384
# 7  2019-08-01  10.774    -2.713          -2.713
# 8  2019-09-01  10.930    -4.282          -4.282
# 9  2019-10-01  11.081    -3.536          -3.536
# 10 2019-11-01  11.237    -3.690          -3.690
# 11 2019-12-01  11.388    -1.023          -1.023
# 12 2020-01-01  11.545     0.904           0.904
# 13 2020-02-01  11.701     1.167           1.167

Prophet fits a piecewise-linear trend with automatic changepoint detection, plus a Fourier-series seasonal term. The “residual” is whatever is in y - trend - additive_terms. The big win with Prophet is that it natively handles multiple seasonalities (yearly + weekly + daily), holidays, and saturating logistic trends — but the components it produces are optimized for forecasting, not for being clean structural descriptors of your series.

Synthetic series construction. The signal is 0.15 * t + 10 (a slow linear ramp from 10 to about 19 over 5 years) plus a 4-unit-amplitude sine of period 12 (one cycle per year in monthly data). The noise is Normal(0, 1.5) — about 1.5 units of jitter per observation, which is small relative to the 8-unit swing of the seasonal term. freq="MS" is month-start, so the index is clean calendar dates that Prophet’s ds column requires.

Classical decomposition. seasonal_decompose(series, model="additive", period=12) does two things: it estimates the seasonal pattern by averaging each calendar position across all cycles (detrended by a centered moving average first), then subtracts it and smoothes the rest to get the trend. The result has NaN at the edges because the moving average needs period/2 past and future observations. You cannot extrapolate the trend out of sample without writing your own code; this is a decomposition tool, not a forecaster.

STL. STL(series, period=12, robust=True).fit() uses LOESS for both the trend and the seasonal smoother. Two windows matter: trend (how smooth the trend is) and seasonal (how much the seasonal pattern is allowed to evolve over time) — both are lengths of the underlying LOESS smoothers, not separate _window-suffixed arguments. robust=True reweights outliers in the residual so they don’t contaminate the trend. .weights returns the per-observation robustness weight; values much less than 1 are outliers STL downweighted — on this series that’s 31 of 60 points, several down near 0.0.

Prophet. df = series.reset_index() and renaming to ds / y is non-negotiable — Prophet expects exactly those column names. yearly_seasonality=True enables a 10-term Fourier series for the annual cycle. The forecast DataFrame’s real columns are ds, trend, yhat_lower, yhat_upper, trend_lower, trend_upper, additive_terms (and its bounds), yearly (and its bounds; the Fourier seasonal contribution), multiplicative_terms (and its bounds; zero in additive mode), and yhat — there is no components column, so selecting one raises KeyError. additive_terms equals yearly here because yearly is the only additive term in play. The implicit residual is y - trend - additive_terms. Prophet is the only one of the three that can produce a forecast past the end of the training data without you writing extra code — that’s the main reason it’s tempting to reach for it.

What’s the same across all three, and what isn’t. All three trends grow from around 10 toward the high teens, and all three seasonal columns oscillate with roughly the article’s 4-unit amplitude — the shapes agree. The exact levels do not: the largest pairwise gap between any two trend columns on the overlapping dates is about 0.76 (classical vs. Prophet), roughly four times the size that “within ~0.2” would suggest, and each method places its own trend on a slightly different baseline (classical’s moving average has no intercept to fit; STL and Prophet each fit one differently). The residual is small for all three (mostly within about ±1 to ±2 of zero, consistent with the noise’s true standard deviation of 1.5). The three methods agree on the signals — there is a linear-ish trend and a roughly 12-month cycle — without agreeing closely enough on the levels that you should mix trend values from different methods.

What’s different. Classical gives you a clean, perfectly periodic seasonal pattern and a moving-average trend that’s NaN at the edges. STL gives you a smooth trend everywhere and a seasonal pattern that can drift. Prophet gives you a piecewise-linear trend with changepoints and a Fourier seasonal pattern, plus the ability to forecast, add holidays, and stack multiple seasonalities.

ACF / PACF / ADF on the same series

from statsmodels.tsa.stattools import acf, pacf, adfuller

acf_vals = acf(series, nlags=24)
pacf_vals = pacf(series, nlags=24, method="ywm")
print("lag  acf   pacf")
for k in range(6):
    print(f"{k:3d}  {acf_vals[k]:+.3f}  {pacf_vals[k]:+.3f}")
# lag  acf   pacf
#   0  +1.000  +1.000
#   1  +0.753  +0.753   <- ACF and PACF are always identical at lag 1: there
#                           are no intermediate lags left to partial out.
#   2  +0.573  +0.016
#   3  +0.322  -0.263
#   4  +0.068  -0.235
#   5  -0.087  +0.025

adf_stat, adf_p, *_ = adfuller(series)
print(f"ADF stat = {adf_stat:.3f}, p = {adf_p:.4f}")
# ADF stat = 0.134, p = 0.9683  -> not stationary (don't reject unit root)

diffed = series.diff().dropna()
_, p_diff, *_ = adfuller(diffed)
print(f"1st difference: p = {p_diff:.4f}")
# 1st difference: p = 0.0000  -> now stationary

The ACF here is not a smooth, monotone decay — it falls to -0.19 by lag 6, then climbs back to +0.518 by lag 12, a bump at the seasonal period. Read literally, the “Reading the plots” rule above (“ACF has significant spikes at lags $m, 2m, \ldots$ → seasonal ARMA territory; you need a seasonal difference”) points at seasonal differencing, not a plain first difference. What’s actually going on: this series’ seasonal component is deterministic (a fixed sine riding on the trend), not a stochastic seasonal unit root, so it shows up as an oscillation in the ACF without being the kind of non-stationarity a first difference exists to fix. Sure enough, one ordinary difference is enough — the ADF p-value drops from 0.9683 to effectively 0, because what needed removing was the linear trend. The sine doesn’t vanish (the differenced series’ own ACF still has a bump of 0.345 at lag 12, down from 0.518), it just isn’t non-stationary, so ADF stops objecting. The takeaway: an oscillating ACF alone doesn’t tell you whether to reach for a seasonal difference — check whether the oscillation survives an ordinary difference and the ADF test along with it before concluding you need $D=1$ as well as $d=1$.

Decision tree: is your series stationary, and what do you do about it?

Run the ADF (or KPSS) test, then look at the ACF plot, then decide.

• ADF p < 0.05 and ACF cuts off quickly (or tails off fast)
  • Series is stationary. Don’t difference.
  • Fit ARMA / ARIMA with $d=0$, or fit a model on the raw series.
• ADF p > 0.05 and ACF decays very slowly (linear-trend signature)
  • Series has a unit root or a deterministic linear trend.
  • Two options:
    • Difference once. Best when the trend is a random walk with drift, or when you’re going to fit ARIMA. The differenced series is stationary, fit ARMA on it, forecast, then “integrate” back. Cost: you lose the long-run level in the model and you can’t recover deterministic drift without drift terms.
    • Model the trend explicitly. Best when the trend is a clean deterministic function of time (linear, polynomial, piecewise-linear). Use Prophet, or a regression with time as a feature, or trend="additive" in a decomposition. Cost: extrapolating polynomials is dangerous; you have to pick a functional form.
• ADF p > 0.05 and ACF has spikes at multiples of $m$ (seasonal non-stationarity)
  • Apply a seasonal difference $\nabla_m y_t$ in addition to (or instead of) the regular difference.
  • This is what SARIMA’s $D=1$ term does. See ARIMA and SARIMA.
• Variance grows with the level (multiplicative)
  • Log-transform first, then treat as additive. ADF on np.log(series) instead of series.
• ADF is borderline (p ≈ 0.05) and you’re not sure
  • Plot rolling mean and rolling variance over a window of size period. If either drifts visibly, difference. If both look flat, don’t.
  • Run KPSS too; KPSS has the null of stationarity, so a low KPSS p-value means non-stationary. Cross-checking ADF and KPSS is more reliable than either alone.

When to difference vs. when to model the trend explicitly

• Difference when:
  • You’re going to fit ARIMA / SARIMA — differencing is literally what the $d$ and $D$ parameters do.
  • The trend is a stochastic random walk (e.g. stock prices, some sensor signals). The future level is genuinely unknown; differencing honestly represents that.
  • You want to forecast the change and then integrate, not the level directly.
• Model the trend explicitly when:
  • The trend is a clean deterministic shape (linear, piecewise linear, saturating logistic) and you trust it to continue.
  • You want interpretable components (“the trend is growing 0.15/month”) rather than a black-box differenced model.
  • You’re using Prophet, which fits the trend this way by design.
  • You’re using ML models (XGBoost, LightGBM) and want to hand them time and time^2 as features rather than differencing away the structure they need to learn from.
• Combine them when:
  • You have a deterministic trend and the residuals still show autocorrelation after detrending — neither pure differencing nor pure trend-modeling is removing all the structure on its own.
  • Detrend with a regression on time (or fit Prophet), then difference the residuals if they still show autocorrelation, or fit ARIMA on the residual directly. This is what Prophet’s additive mode + a residual ARIMA does well.
  • It costs a second fitted stage and compounds extrapolation error if the trend model is wrong — worth it specifically when a single approach leaves autocorrelated residuals behind, not as a default over the two options above. A pure random walk, for instance, is usually better served by differencing alone.

Edge cases and common mistakes

• Using multiplicative decomposition without checking the variance. If your series has roughly constant variance but a growing level, additive is correct and multiplicative will quietly inflate the seasonal swings in the later years. Plot the rolling std; if it’s flat, use additive.
• Trusting the ADF on short series. ADF has low power; on series shorter than ~50 observations it routinely fails to reject the unit-root null even when the series is stationary. Always look at the ACF plot too.
• Differencing twice when once was enough. Over-differencing introduces an MA unit root and inflates the forecast variance. If adfuller(series.diff().dropna()) says stationary, stop.
• Forgetting that seasonal_decompose can’t extrapolate the trend. It produces a trend only at the interior observations; you have to fit a regression on the trend yourself if you want to project it forward.
• Comparing STL and classical components as if they were the same object. Classical seasonal is exactly periodic and forced to sum to zero; STL seasonal is allowed to evolve. They’re not interchangeable; don’t be surprised when they disagree at the edges.
• Setting period wrong. period must be the number of observations in one cycle, not the cycle in human units. For monthly data with an annual cycle it’s 12, for daily data with a weekly cycle it’s 7, for hourly data with a daily cycle it’s 24. Getting this wrong makes the seasonal component look like noise.
• Prophet without changepoint_prior_scale tuning. The default often under-fits the trend (too few changepoints) or over-fits (too many). Plot the fitted trend against the data before trusting the components.
• Conflating stationarity with “no trend”. Stationarity is about the joint distribution, not just the mean. A series with constant mean but variance that doubles halfway through is non-stationary, and differencing won’t fix that — you need a variance-stabilizing transform (log, Box-Cox).
• Using method="ols" for PACF. statsmodels’s default is method="ywadjusted", not OLS — but explicitly passing method="ols" can give values outside $[-1, 1]$ on small samples; use method="ywm" (Yule-Walker, MLE-adjusted) for stable estimates.
• Treating the residual as “noise” without checking it. Plot the residual’s ACF; if it’s still autocorrelated, you’ve under-fit. The residual should look like white noise. If it has structure, go back and add an AR term or a different seasonal period.

Cross-references

• ARIMA and SARIMA — uses the stationarity machinery and the ACF/PACF decision rules above to pick $p$, $d$, $q$, $P$, $D$, $Q$.
• Prophet vs. statistical models) — the practical comparison that motivates the “model the trend explicitly” branch of the decision tree.
• Classical forecasting vs. machine learning) — when the structural decomposition approach wins vs. when feature-engineered tree models win.
• Feature engineering for time series) — uses lags, rolling windows, and the autocorrelation structure described here to build ML features.
• Why your forecast fails on Christmas) — the seasonal and holiday pieces of the decomposition, and what breaks when you treat them as ordinary seasonality.
• Time series cross-validation) — how the trend/seasonality/residual split interacts with proper backtesting; K-fold breaks because the trend leaks forward.
• How to detect and handle data drift in production) — uses the residual component and rolling-stat monitoring to spot drift.

Further reading

• Box, G. E. P., & Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control. Holden-Day. — the original ARIMA reference; the source of the differencing / ARMA identification machinery.
• Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed.), free online at https://otexts.com/fpp3/ — the modern practical textbook; chapters 1-3 cover trend, seasonality, and decomposition; chapter 9 covers ARIMA.
• Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A Seasonal-Trend Decomposition Procedure Based on LOESS. Journal of Official Statistics, 6(1), 3-73. — the original STL paper.
• Taylor, S. J., & Letham, B. (2018). Forecasting at Scale. The American Statistician, 72(1), 37-45. — the Prophet paper.
• statsmodels documentation: seasonal_decompose, STL, adfuller, acf/pacf.
• Prophet documentation — especially the sections on saturating forecasts and changepoints.
• Kaggle competition: M5 Forecasting - Uncertainty — Walmart retail sales with strong trend + multiple seasonalities + holiday effects; a good test bed for the decomposition choices in this reference.