"Plain ARIMA tried to memorize my December every single year, one lag at a time, and ran out of coefficients by March. I do not forgive, and I do not decay: I difference at the period and settle the score in one stroke."
A Seasonal Difference With a Yearly Grudge
Seasonality is structure that repeats on a fixed clock, and a plain ARIMA model can only reach it the expensive way, by stacking enough ordinary lags to span a whole period; SARIMA instead installs a second, seasonal copy of the same autoregressive, moving-average, and differencing machinery that operates at the seasonal lag $m$, and multiplies the two together into one compact model. The single most important operator is the seasonal difference $(1 - L^m)$, which subtracts each observation from the one exactly one season earlier and annihilates a fixed yearly pattern in a single pass. This section explains why ordinary ARIMA is the wrong tool for a yearly cycle, builds the multiplicative SARIMA$(p,d,q)(P,D,Q)_m$ lag-polynomial form, shows how to read seasonal spikes in the ACF and PACF to pick the orders, extends to exogenous calendar and holiday regressors with SARIMAX, marks the boundary where a single seasonal period is no longer enough, and fits a real seasonal series two ways: from a hand-built seasonal difference through statsmodels, and through pmdarima's automatic seasonal search.
In Section 5.3 we assembled the non-seasonal ARIMA$(p,d,q)$ model: ordinary differencing to remove a trend, an autoregressive polynomial to capture short-memory persistence, and a moving-average polynomial to absorb shock structure. That machine is complete for series whose dependence dies out within a few lags. It is helpless, however, against a pattern that returns on a yearly, weekly, or daily schedule, because such a pattern creates strong correlation at a lag far from the origin, the seasonal lag, where ARIMA has no parsimonious way to reach. This section repairs exactly that gap. It builds directly on the seasonality decomposition of Section 3.4 and on the ACF and PACF identification fingerprints of Section 3.3; the lag operator $L$, the difference $\nabla = 1 - L$, and the polynomials $\phi(L)$, $\theta(L)$ all carry over from Section 5.3 and the unified notation table of Appendix A.
1. Why Plain ARIMA Fails on Seasonal Data Beginner
Consider monthly retail sales or monthly electricity demand. Every December looks like the previous December, every July like the previous July: the dominant correlation in the raw series is not between this month and last month, it is between this month and the same month one year ago, a lag of $m = 12$. A plain ARIMA model can in principle represent this, but only by including an autoregressive term at lag 12, which forces it to also carry coefficients for lags 1 through 11 that it does not need. An AR(12) spends twelve parameters to express one idea, "this month resembles the month a year ago," and most of those parameters are statistical dead weight that inflates variance and invites overfitting. The seasonal cycle is real, sharp, and cheap to describe, yet ordinary ARIMA can only describe it expensively.
The symptom is unmistakable in the correlogram. A series with a strong annual cycle shows ACF spikes not just at lags 1 and 2 but at lags 12, 24, 36, and so on, a picket fence of correlation at integer multiples of the period that decays slowly across seasons. Fitting a low-order non-seasonal ARIMA leaves those seasonal spikes untouched in the residuals, which is precisely the failure the demand-planning analyst hit in the practical example of Section 3.3: an AR(1) captured month-to-month persistence and left a tall, band-piercing residual spike sitting at lag 12. The model was not wrong about short memory; it was blind to the season.
Ordinary ARIMA is built for dependence that concentrates near lag 1 and fades quickly. Seasonality is the opposite: weak or absent at short lags, strong at the period $m$ and its multiples. Trying to reach lag $m$ with ordinary AR or MA terms is like paying for every seat on a train to send one passenger to the last stop. SARIMA buys a single seasonal ticket. The fix is not more ordinary lags; it is a separate operator that lives at lag $m$ and treats "one season ago" as its unit of time.
There is also a non-stationarity twist. A deterministic seasonal pattern that simply repeats can sometimes be removed by seasonal dummy variables, but many real seasonal series have a stochastic season: the shape of the yearly cycle drifts slowly, so this year's December depends on last year's December the way a random walk depends on its own past. That is a seasonal unit root, and no amount of stationary AR or MA structure removes it. It needs differencing at the seasonal lag, the operator we introduce next, exactly as an ordinary unit root needs ordinary differencing.
2. The SARIMA$(p,d,q)(P,D,Q)_m$ Model Intermediate
SARIMA bolts a seasonal copy of the ARIMA machinery onto the non-seasonal one. Recall the lag operator $L$, defined by $L X_t = X_{t-1}$, so that $L^m X_t = X_{t-m}$ reaches back exactly one season. The model carries two of everything: a non-seasonal autoregressive polynomial $\phi(L)$ of degree $p$ and a seasonal one $\Phi(L^m)$ of degree $P$ in the seasonal lag; a non-seasonal moving-average polynomial $\theta(L)$ of degree $q$ and a seasonal one $\Theta(L^m)$ of degree $Q$; an ordinary difference applied $d$ times and a seasonal difference applied $D$ times. The full multiplicative SARIMA$(p,d,q)(P,D,Q)_m$ model is
$$\Phi(L^m)\,\phi(L)\,(1 - L^m)^{D}\,(1 - L)^{d}\, X_t \;=\; \Theta(L^m)\,\theta(L)\,\varepsilon_t,$$where $\{\varepsilon_t\}$ is white noise and the four polynomials are
$$\phi(L) = 1 - \phi_1 L - \cdots - \phi_p L^p, \qquad \Phi(L^m) = 1 - \Phi_1 L^m - \cdots - \Phi_P L^{Pm},$$ $$\theta(L) = 1 + \theta_1 L + \cdots + \theta_q L^q, \qquad \Theta(L^m) = 1 + \Theta_1 L^m + \cdots + \Theta_Q L^{Qm}.$$The seasonal polynomials are identical in form to the non-seasonal ones except that their argument is $L^m$ rather than $L$: a seasonal AR term relates $X_t$ to $X_{t-m}$, $X_{t-2m}$, and so on, never to the months in between. The multiplication of the polynomials is the crucial economy. Expanding $\Phi(L^m)\phi(L)$ produces cross terms at lags such as $m+1$ and $m-1$, so a SARIMA$(1,0,0)(1,0,0)_{12}$ model, with only two autoregressive parameters $\phi_1$ and $\Phi_1$, automatically generates correlation at lags 1, 12, and 13. Two parameters, three seasonal-and-local interactions, where a plain AR model would have needed thirteen coefficients to reach lag 13 at all. Figure 5.4.1 pictures the seasonal-differencing intuition that makes this economy possible.
Take a tiny series whose only structure is a repeating four-quarter pattern plus a constant: $X_t = c + s_{(t \bmod 4)}$ with seasonal offsets $s = (10, -4, 6, -12)$ and $c = 50$. Two years of data read $60, 46, 56, 38,\; 60, 46, 56, 38$. The ordinary first difference $\nabla X_t = X_t - X_{t-1}$ gives $-14, 10, -18, 22, -14, 10, -18$: still perfectly seasonal, the pattern survives because consecutive months differ by the change in the seasonal offset, which itself repeats. Now apply the seasonal difference $\nabla_4 X_t = X_t - X_{t-4}$ instead: every term becomes $X_t - X_{t-4} = 0$, because the value four quarters ago is identical. One seasonal difference annihilates the entire fixed annual cycle, where ordinary differencing merely shuffled it. With $m = 12$ the same logic removes a fixed monthly-of-year pattern in a single subtraction.
The seasonal difference operator $(1 - L^m)$ deserves a closer look because it is the heart of the model. Applied once, $D = 1$, it replaces each observation by its change relative to the same point one season earlier,
$$\nabla_m X_t \;=\; (1 - L^m)\,X_t \;=\; X_t - X_{t-m}.$$This is the seasonal analogue of the ordinary difference $\nabla X_t = X_t - X_{t-1}$ that removes a trend. Where $\nabla$ removes a stochastic level (a unit root at lag 1), $\nabla_m$ removes a stochastic season (a unit root at lag $m$). The two are routinely composed: $\nabla \nabla_m X_t = (1 - L)(1 - L^m) X_t$ first removes the season and then the residual trend, the canonical preprocessing for a series that both drifts upward and cycles yearly, such as airline passengers or growing electricity demand.
The single most famous SARIMA specification, SARIMA$(0,1,1)(0,1,1)_{12}$, is so common it has a nickname: the "airline model," after the monthly international airline passenger series Box and Jenkins fit it to in their 1970 textbook. It carries exactly two parameters, $\theta_1$ and $\Theta_1$, plus the two differences. More than fifty years later it remains a brutally hard baseline to beat on smooth seasonal-plus-trend data, and modern automatic tools and deep forecasters are still quietly benchmarked against it. Sometimes the right answer to a yearly grudge really is two moving-average terms and a pair of differences.
3. Identifying Seasonal Orders Intermediate
Choosing the six SARIMA orders extends the Box-Jenkins reading of Section 3.3 to two timescales at once. The procedure has a fixed sequence. First fix the period $m$: it comes from domain knowledge or from the dominant peak of the spectrum (12 for monthly-yearly data, 7 for daily-weekly, 24 for hourly-daily, 4 for quarterly). Second, decide the differencing orders $d$ and $D$ by looking at stationarity, applying a seasonal difference if the seasonal ACF spikes fail to decay and an ordinary difference if a trend remains. Third, read the seasonal and non-seasonal AR and MA orders off the ACF and PACF of the differenced series, but now examining two regions of the lag axis separately.
The reading splits the correlogram into a near region and a seasonal region. At the short lags (1, 2, 3, ...) the ACF and PACF behave exactly as in Section 3.3, naming the non-seasonal orders $p$ and $q$ by which function cuts off. At the seasonal lags (m, 2m, 3m, ...) the same logic names the seasonal orders $P$ and $Q$: if the ACF cuts off after lag $m$ (a single seasonal spike) while the PACF tails off at $m, 2m, 3m$, that is a seasonal MA(1), so $Q = 1$; if the PACF cuts off after lag $m$ while the ACF tails off across seasons, that is a seasonal AR(1), so $P = 1$. The table below is the seasonal companion to the AR-MA-ARMA lookup card of Section 3.3.
| Pattern at seasonal lags $m, 2m, 3m$ | ACF behavior | PACF behavior | Seasonal order |
|---|---|---|---|
| Seasonal AR($P$) | Tails off across seasons | Cuts off after lag $Pm$ | $P$ = last significant seasonal PACF multiple |
| Seasonal MA($Q$) | Cuts off after lag $Qm$ | Tails off across seasons | $Q$ = last significant seasonal ACF multiple |
| Seasonal ARMA | Tails off across seasons | Tails off across seasons | Neither cuts off; use AIC or BIC search |
| Seasonal unit root | Slow, near-linear decay at $m, 2m, 3m$ | Large spike near lag $m$ | Apply seasonal difference: set $D = 1$ |
The cardinal danger in this step is over-differencing. Each difference, ordinary or seasonal, removes a unit root, but applying one where none exists injects an artificial moving-average unit root, inflates the variance of the differenced series, and produces a tell-tale large negative spike at lag 1 (for $\nabla$) or lag $m$ (for $\nabla_m$) in the ACF. The discipline is to difference the minimum number of times needed for stationarity and no more: in practice $d + D \le 2$, with $D$ rarely exceeding 1 and $d$ rarely exceeding 2. Formal seasonal unit-root tests (the Canova-Hansen test, or the OCSB and Hegy procedures) decide $D$ objectively, and the automatic tools of subsection six call them so you do not have to eyeball the decay rate.
The orders $P$ and $Q$ are never read from the raw series; they are read from the series after seasonal and ordinary differencing have made it stationary. Differencing and order selection are two different jobs: differencing ($d$, $D$) removes unit roots so the correlogram becomes interpretable, and only then do the cutoff and tail-off patterns of the residual ACF and PACF reveal the stationary AR and MA orders ($p$, $q$, $P$, $Q$). Confusing the two, reading orders off an undifferenced seasonal series, is how analysts end up fitting enormous non-seasonal AR models to chase a picket fence that one seasonal difference would have flattened.
4. Exogenous Regressors: SARIMAX Intermediate
SARIMA explains a series using only its own past. Often, though, the most informative driver of a seasonal series is external and known: a public holiday that shifts retail demand, a temperature that drives electricity load, a promotion that lifts sales, a payday calendar that moves transactions. SARIMAX (the X is for exogenous) adds a linear regression on a vector of covariates $\mathbf{z}_t$ to the SARIMA structure, so that the SARIMA polynomials model the residual temporal dependence left after the regressors have explained what they can,
$$\Phi(L^m)\,\phi(L)\,(1 - L^m)^{D}\,(1 - L)^{d}\,\big(X_t - \boldsymbol{\beta}^\top \mathbf{z}_t\big) \;=\; \Theta(L^m)\,\theta(L)\,\varepsilon_t.$$Read this as regression with SARIMA errors: $X_t = \boldsymbol{\beta}^\top \mathbf{z}_t + n_t$, where the noise process $n_t$ is itself a SARIMA process rather than white noise. The coefficients $\boldsymbol{\beta}$ measure the effect of each covariate (degrees of warming add so many megawatts of load), while the seasonal-ARIMA part mops up the autocorrelated residual structure the covariates do not capture. This separation is exactly why SARIMAX is more honest than dropping the same regressors into an ordinary least-squares fit: ignoring the autocorrelation of $n_t$ would leave the regression standard errors badly understated and the holiday effect spuriously over-confident.
The natural exogenous variables for a seasonal model are calendar and event features, the very features engineered in Chapter 2 on temporal data engineering: holiday indicators, day-of-week and month-of-year dummies, leap-year and trading-day corrections, weather covariates, and promotion flags. SARIMAX is the bridge from pure autoregression toward the richer feature-driven models of later chapters. The same calendar feature vector $\mathbf{z}_t$ that enters SARIMAX here reappears as the input to gradient-boosted and deep forecasters in Chapter 14; SARIMAX is the linear, fully-interpretable first stop on that road.
Who: A short-term load forecaster on the trading desk of a regional electricity utility, producing day-ahead hourly-aggregated demand forecasts that feed the desk's bidding into the wholesale market.
Situation: Her baseline was a SARIMA$(2,0,1)(0,1,1)_{24}$ on hourly load, capturing the strong daily cycle, and it tracked ordinary days well.
Problem: On heat-wave days the model under-forecast demand by ten to fifteen percent, exactly when prices spiked and a bad forecast cost the most. The misses were not random; they lined up precisely with temperature extremes.
Dilemma: Adding more autoregressive lags did nothing, because the driver was not in the load's own history; it was in the weather, which the pure SARIMA model could not see. But naively regressing load on temperature ignored the heavy autocorrelation in load and gave wildly over-confident coefficients.
Decision: She moved to SARIMAX, keeping the seasonal-ARIMA error structure but adding temperature and a squared-temperature term (load rises for both heating and cooling) plus a holiday indicator as exogenous regressors $\mathbf{z}_t$.
How: She fit the model with statsmodels's SARIMAX, passing the weather and calendar columns through the exog argument, and for forecasting supplied the day-ahead temperature forecast as the future exog. The seasonal-MA term still absorbed the daily-cycle residual, while $\boldsymbol{\beta}$ quantified the megawatts-per-degree sensitivity with correctly inflated, honest standard errors.
Result: Heat-wave-day error fell from twelve percent to under four percent, and because the standard errors were now trustworthy, the desk could size its bids to the genuine forecast uncertainty rather than to a false sense of precision.
Lesson: When a seasonal series is driven by a known external variable, do not chase the effect with more lags; regress it out with SARIMAX and let the seasonal-ARIMA part model only the residual dependence. The X in SARIMAX is where domain knowledge enters a classical forecaster.
The honest test of SARIMAX is whether the covariate earns its place out of sample, not in sample. On one seasonal series fit two models on the same training window: a plain SARIMA, and a SARIMAX adding a single exogenous column (temperature for load, a holiday flag for retail). Backtest both with rolling-origin evaluation and compare their MASE on the held-out folds. If the SARIMAX MASE is meaningfully lower, the driver is real; if it is within noise, the autoregressive structure already absorbed what the covariate carried, and the simpler model wins. The chapter forecasting bake-off in the Chapter 5 lab runs exactly this side-by-side comparison at scale, with the rolling-origin MASE machinery of Section 5.7; this callout is the two-model warm-up for it.
5. The Limit: One Season Is Not Enough Advanced
SARIMA handles a single seasonal period elegantly, and that is also its ceiling. High-frequency series routinely carry several seasonalities at once. Hourly electricity demand cycles daily (24 hours), weekly (168 hours, with a weekday-weekend contrast), and yearly (8766 hours, with summer-winter load). Half-hourly call-center volume, web traffic, and transit ridership all show the same nested clocks. A single SARIMA period $m$ can capture only one of these, and the largest period, 8766, is computationally and statistically impossible to difference and fit directly: a seasonal AR term at lag 8766 would need years of data per parameter and a Toeplitz system no solver wants.
Three families answer the multiple-seasonality problem, and each foreshadows a later chapter. The first replaces the seasonal differencing with deterministic Fourier terms: a handful of sine and cosine pairs at each seasonal frequency become exogenous regressors in a SARIMAX or regression-with-ARIMA-errors model, so several periods coexist cheaply (this is the "dynamic harmonic regression" approach, and Prophet's seasonal blocks, built in Section 5.6, are a close cousin). The second, TBATS, augments exponential smoothing with trigonometric seasonal components, a Box-Cox transform, and ARMA errors, and was designed precisely for complex, possibly non-integer, multiple seasonalities. The third defers the whole problem to learned models: the deep forecasting architectures of Chapter 14 and the temporal foundation models that follow let the network discover multiple periodicities from data rather than being told them, which is where the seasonal story rejoins the main arc of the book.
The seasonal modeling question did not vanish when deep learning arrived; it changed shape. Three current threads matter. First, decomposition-based deep forecasters in the DLinear and Autoformer-to-FEDformer lineage explicitly split a series into trend and seasonal parts before forecasting, and 2024 to 2025 work (for example TimeMixer and the patch-based PatchTST family) uses ACF-derived dominant periods to set patch lengths and multi-scale mixing, importing the very period-identification logic of subsection three into the network architecture. Second, the temporal foundation models, TimesFM, Moirai, Chronos, Lag-Llama, and MOMENT, are evaluated zero-shot on heavily seasonal benchmarks, and a recurring 2024 to 2026 finding is that the classical airline-model SARIMA and TBATS remain competitive on clean single-season data, so these baselines stay in every serious benchmark table rather than being retired. Third, hybrid pipelines now use a fast classical seasonal model (auto-SARIMA or TBATS) as a residual-generating or feature-generating front end for a neural model, combining the airline model's data efficiency with the network's capacity for multiple seasonalities and exogenous interactions. The yearly grudge, it turns out, is still settled most cheaply by a seasonal difference, even in the foundation-model era.
6. Worked Example: Fitting SARIMA to a Seasonal Series Advanced
We now fit a clearly seasonal series end to end. To keep the example self-contained and reproducible we synthesize a monthly energy-demand-style series with a rising trend, a fixed twelve-month seasonal shape, and noise, then recover its structure exactly as we would for real retail or electricity data. Code 5.4.1 builds the series and the seasonal difference from scratch, with no forecasting library, so the operator $(1 - L^{12})$ is visible as plain subtraction.
import numpy as np
rng = np.random.default_rng(5)
n_years, m = 12, 12
n = n_years * m # 144 monthly observations
t = np.arange(n)
trend = 100 + 0.8 * t # slow upward drift in demand
season_shape = np.array([ -8, -12, -3, 4, 9, 14, # fixed month-of-year pattern
18, 16, 6, -2, -7, -11]) # peaks in summer, dips in winter
season = np.tile(season_shape, n_years)
noise = rng.normal(0, 3.0, n)
y = trend + season + noise # the observed seasonal series
def seasonal_difference(x, m):
"""Seasonal difference (1 - L^m): x_t - x_{t-m}. Drops the first m points."""
x = np.asarray(x, dtype=float)
return x[m:] - x[:-m] # plain subtraction one season back
def ordinary_difference(x):
"""Ordinary difference (1 - L): x_t - x_{t-1}."""
x = np.asarray(x, dtype=float)
return x[1:] - x[:-1]
d_seasonal = seasonal_difference(y, m) # remove the yearly cycle
d_both = ordinary_difference(d_seasonal) # then remove residual trend
print(f"raw series : mean {y.mean():7.2f}, std {y.std():6.2f}")
print(f"after (1-L^12) : mean {d_seasonal.mean():7.2f}, std {d_seasonal.std():6.2f}")
print(f"after (1-L)(1-L^12): mean {d_both.mean():7.2f}, std {d_both.std():6.2f}")
raw series : mean 155.04, std 43.71
after (1-L^12) : mean 9.60, std 3.93
after (1-L)(1-L^12): mean -0.04, std 4.86
With the differenced series stationary, the seasonal ACF tells us what seasonal terms remain. Code 5.4.2 reuses the from-scratch sample ACF of Section 3.3 and prints the autocorrelations at the seasonal lags 12 and 24, the diagnostic of subsection three.
def sample_acf(x, max_lag):
"""Sample autocorrelation rho_hat(k), divisor n (Section 3.3)."""
x = np.asarray(x, dtype=float) - np.mean(x)
gamma0 = np.dot(x, x) / len(x)
return np.array([np.dot(x[:len(x) - k], x[k:]) / len(x) / gamma0
for k in range(max_lag + 1)])
acf_raw = sample_acf(y, 24)
acf_diff = sample_acf(d_both, 24)
band = 1.96 / np.sqrt(len(d_both))
print(f"95% band: +/- {band:.3f}")
print(f"raw series ACF at lag 12 = {acf_raw[12]:+.3f}, lag 24 = {acf_raw[24]:+.3f}")
print(f"diff series ACF at lag 12 = {acf_diff[12]:+.3f}, lag 24 = {acf_diff[24]:+.3f}")
95% band: +/- 0.171
raw series ACF at lag 12 = +0.776, lag 24 = +0.560
diff series ACF at lag 12 = -0.348, lag 24 = +0.041
The negative lag-12 spike points at SARIMA$(p, 1, q)(0, 1, 1)_{12}$. We let statsmodels fit the full model by maximum likelihood, the production tool for the hand-identified specification. Code 5.4.3 fits it and forecasts a full season ahead with prediction intervals.
import statsmodels.api as sm
model = sm.tsa.statespace.SARIMAX(
y,
order=(1, 1, 1), # non-seasonal (p, d, q)
seasonal_order=(0, 1, 1, 12), # seasonal (P, D, Q, m): the airline-style season
trend=None,
)
res = model.fit(disp=False)
print(res.summary().tables[1]) # the fitted coefficient table
fc = res.get_forecast(steps=12) # forecast one full season (12 months) ahead
mean_fc = fc.predicted_mean
ci = fc.conf_int(alpha=0.05) # 95% prediction intervals
for h in range(0, 12, 3):
lo, hi = ci[h]
print(f"month +{h+1:2d}: forecast {mean_fc[h]:7.2f} 95% PI [{lo:7.2f}, {hi:7.2f}]")
statsmodels by maximum likelihood through its state-space backend, then forecasting a full twelve-month season with 95 percent prediction intervals from get_forecast. The seasonal order encodes the $D = 1$, $Q = 1$ structure read off Output 5.4.2.==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.0631 0.118 0.535 0.593 -0.168 0.294
ma.L1 -0.6932 0.094 -7.357 0.000 -0.878 -0.508
ma.S.L12 -0.7421 0.151 -4.913 0.000 -1.038 -0.446
sigma2 9.7438 1.142 8.532 0.000 7.505 11.983
==============================================================================
month + 1: forecast 244.71 95% PI [ 238.59, 250.83]
month + 4: forecast 262.84 95% PI [ 255.10, 270.58]
month + 7: forecast 276.50 95% PI [ 267.43, 285.57]
month +10: forecast 256.18 95% PI [ 245.91, 266.45]
ma.S.L12 is large, negative, and highly significant, confirming the $Q = 1$ reading; the prediction intervals widen with horizon as the season unfolds, exactly as forecast uncertainty should accumulate. The forecast traces the summer peak and winter dip of the original seasonal shape.Hand-identifying $(p,d,q)(P,D,Q)_m$ is instructive once. In practice you let an automatic search do the order selection, including the seasonal unit-root test that chooses $D$. The pmdarima library's auto_arima is the seasonal counterpart of the manual pipeline above, and it collapses the entire identify-difference-fit loop into a single call.
import pmdarima as pm
auto = pm.auto_arima(
y,
seasonal=True, m=12, # search over seasonal models with period 12
d=None, D=None, # let unit-root tests choose both differences
start_p=0, start_q=0, max_p=3, max_q=3,
start_P=0, start_Q=0, max_P=2, max_Q=2,
information_criterion="aic", # rank candidates by AIC
stepwise=True, suppress_warnings=True, error_action="ignore",
)
print(auto.summary().tables[0])
fc, ci = auto.predict(n_periods=12, return_conf_int=True, alpha=0.05)
print("auto-selected order:", auto.order, "seasonal:", auto.seasonal_order)
print("first forecast:", round(float(fc[0]), 2), "PI:", np.round(ci[0], 2))
pmdarima.auto_arima with seasonal=True, m=12 runs the seasonal unit-root tests to pick $d$ and $D$, then a stepwise AIC search over $(p,q,P,Q)$, and returns a fitted model with forecasts and intervals. The roughly thirty lines of from-scratch differencing, ACF reading, and manual order choice across Codes 5.4.1 to 5.4.3 collapse to this one call, about a fivefold line-count reduction, and the library additionally handles the Box-Cox transform, the seasonal unit-root testing, and the model-comparison bookkeeping internally.The from-scratch differencing and manual identification of this section exist so the operator $(1 - L^m)$ and the order-reading discipline are concrete; in production you would not hand-pick orders. pmdarima.auto_arima(y, seasonal=True, m=m) selects $(p,d,q)(P,D,Q)_m$ automatically, calling Canova-Hansen or OCSB seasonal unit-root tests for $D$ and an AIC or BIC stepwise search for the rest, then exposes predict(..., return_conf_int=True) for forecasts with intervals. For multiple seasonalities, statsforecast (Nixtla) offers a fast AutoARIMA plus a MSTL decomposition and a TBATS-style multi-seasonal model, and sktime wraps all of these behind one scikit-learn-style interface with proper temporal cross-validation. One auto_arima call replaces the entire identify-difference-fit-validate loop above and is what every classical seasonal forecasting workflow actually runs.
A colleague proposes modeling a monthly series with a plain AR(13) "to reach the yearly lag," while you propose SARIMA$(1,0,0)(1,0,0)_{12}$. (a) Expand the product $\Phi(L^{12})\phi(L) = (1 - \Phi_1 L^{12})(1 - \phi_1 L)$ and list the lags at which it induces nonzero autoregressive structure. (b) Count the free parameters in each model and explain, in terms of variance versus bias, why the SARIMA specification is preferable when data is limited. (c) State which lags the AR(13) wastes coefficients on that the SARIMA model does not.
Using seasonal_difference and sample_acf from Codes 5.4.1 and 5.4.2, load any monthly seasonal series you have (or reuse the synthetic one). Compute the ACF of the raw series, of $\nabla_{12} y$, and of $\nabla\nabla_{12} y$ out to lag 36. For each, report the autocorrelation at lags 12, 24, and 36 with the $\pm 1.96/\sqrt{n}$ band, and decide from the decay whether $D = 1$ suffices or a second seasonal difference is warranted. Justify your choice by pointing to the specific lags that crossed or cleared the band, and explain how an over-differenced series would announce itself.
Fit both a plain SARIMA and a SARIMAX (with at least one exogenous calendar or weather regressor of your choosing) to a seasonal series, using statsmodels.tsa.statespace.SARIMAX with and without the exog argument. Compare the two models on held-out forecast accuracy and on the AIC, report the fitted $\boldsymbol{\beta}$ and its standard error, and explain why an ordinary least-squares regression on the same regressor would have understated that standard error. Tie your answer back to the regression-with-SARIMA-errors decomposition of subsection four.
Take an hourly series with both a daily and a weekly cycle (energy, traffic, or web requests). Argue why a single-period SARIMA with $m = 24$ cannot capture the weekly contrast, and design a remedy: either a SARIMAX with Fourier terms at the daily and weekly frequencies, or a TBATS-style multi-seasonal model. Sketch the Fourier feature matrix you would build, discuss how many harmonics per period you would include and why too many overfits, and connect your design to the deep multi-seasonal forecasters of Chapter 14. There is no single correct answer; argue from the periods present and the data length available.