Meet Nora, demand planner at a regional bakery chain. Her boss wants to know how many sourdough loaves the stores will sell next week. She has two years of daily sales data. Naturally, she reaches for a Long Short-Term Memory (LSTM) neural network. It’s complex. It sounds great in a meeting.

But here’s the catch: when you run the model, it predicts a flat line. Or worse, a wild spike that makes no sense. Meanwhile, the bakery owner’s simple rule—“we usually sell what we sold last week, plus a little more for the holiday”—beats the high-tech model by 20%.

What happened? More parameters do not always mean better predictions. Deep learning models need a lot of data—thousands, sometimes millions of rows—to grasp the shape of a time series. With only a few hundred rows, they overfit to noise or miss the obvious trend. This is where classical forecasting—ARIMA and SARIMA—comes in. They aren’t obsolete; they’re specialized tools for when you need to work with less.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Simulating a small retail dataset with a clear trend
np.random.seed(42)
time = np.arange(100)
sales = 10 + 0.5 * time + np.random.normal(0, 2, 100)

# A 'complex' model might struggle with such little data
# while a simple trend line captures the essence
plt.figure(figsize=(10, 5))
plt.plot(time, sales, label='Actual Sales (Small Data)')
plt.plot(time, 10 + 0.5 * time, '--', label='Simple Trend Line', color='red')
plt.title("Small Data: Complexity vs. Intuition")
plt.legend()
plt.show()

• np.random.seed(42) — locks the random number generator so the synthetic noise is reproducible; change the seed and the wiggles move.
• np.arange(100) — generates integers 0 through 99, serving as the time index.
• sales = 10 + 0.5 * time + np.random.normal(0, 2, 100) — builds the series from three pieces: an intercept of 10, a slope of 0.5 units per time step, and Gaussian noise with standard deviation 2. This is the classic “trend + noise” recipe.
• plt.plot(..., '--') — the dashed red line plots the deterministic trend (intercept + slope) without the noise, so you can visually compare what the data actually did versus the underlying signal a simple model would capture.

In the plot above, the red line represents the logic of a classical model. It sees the trend right away. A neural network might spend its entire training time trying to figure out if the random wiggles in the blue line matter. They don’t.

The Three Ingredients of a Time Series

You don’t need a math degree for ARIMA. Three ingredients do the work: Momentum, Cleaning, and Shocks. The name itself breaks down as Auto-Regressive (AR), Integrated (I), and Moving Average (MA).

An ARIMA$(p, d, q)$ model combines three components into one equation:

$\phi(B)\,(1 - B)^{d}\, y_t \;=\; \theta(B)\,\varepsilon_t$

where $B$ is the backshift operator ($By_t = y_{t-1}$), $\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$ is the AR polynomial, $\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$ is the MA polynomial, and $\varepsilon_t$ is white noise.

Plain English	Statistical symbol	Python equivalent
Momentum: “how many past values to remember”	$p$ (AR order), $\phi_i$	order=(p, …, …) in SARIMAX
Cleaning: “how many times to difference”	$d$ (differencing order)	order=(…, d, …) in SARIMAX
Shocks: “how many past errors to correct”	$q$ (MA order), $\theta_j$	order=(…, …, q) in SARIMAX
Seasonal cycle (SARIMA only)	$(P, D, Q, s)$	seasonal_order=(P, D, Q, s)

1. AR (Auto-Regressive): The ‘Yesterday Matters’ Factor

Treat this as momentum. Drive 60 mph yesterday, and you’ll likely be near 60 mph today. The model uses its own past values to predict what comes next. AR puts it plainly: “Tell me what happened yesterday, and I’ll tell you what happens today.”

2. I (Integrated): Making the Data Sit Still

Most time series drift up or down over time — that’s the trend. Classical models need the opposite: stationary data that fluctuates around a constant average. This trips people up, but the logic is straightforward once you see it.

The fix is differencing. Rather than predicting a house’s price directly, we predict how much the price changed from yesterday to today.

# Let's see what differencing actually does
data = pd.Series([10, 12, 15, 18, 22]) # A clear upward trend
diffed_data = data.diff().dropna()    # Subtracting yesterday from today

print(f"Original Data: {data.values}")
print(f"Differenced Data: {diffed_data.values}")
# The output [2, 3, 3, 4] is much 'flatter' than [10, 12, 15, 18, 22]

• pd.Series([...]) — wraps the raw list in a pandas Series so we get access to time-series methods like .diff().
• .diff() — computes $y_t - y_{t-1}$ for every position; the first element has no predecessor, so it becomes NaN.
• .dropna() — removes that leading NaN so downstream code doesn’t choke on a missing value.
• The result [2, 3, 3, 4] is the sequence of changes—a much flatter, more stationary series that a model can work with instead of the raw trending values.

3. MA (Moving Average): Dealing with Shocks

Sometimes, something weird happens—a sudden snowstorm closes the shop. That’s a “shock.” The MA part of the model doesn’t look at the past values themselves, but at the past errors. It asks: “How much was I wrong by yesterday? Let me use that mistake to adjust my guess for today.”

SARIMA: Adding the Calendar to the Mix

Standard ARIMA has a blind spot: it doesn’t understand calendars. Sell ice cream and you know July is busy, January is slow. Run regular ARIMA on that data, though, and it sees the July spike, assumes the trend is climbing, and predicts an even higher August.

SARIMA adds an ‘S’ for Seasonality. The model can now remember what happened exactly one year ago — or one week ago. It treats the seasonal pattern as a separate, repeating cycle.

from statsmodels.tsa.seasonal import seasonal_decompose

# Simulating seasonal data: Trend + Sin Wave (Seasonality) + Noise
t = np.arange(120)
seasonal_sales = 10 + 0.2 * t + 10 * np.sin(2 * np.pi * t / 12) + np.random.normal(0, 1, 120)
df = pd.Series(seasonal_sales, index=pd.date_range('2010-01-01', periods=120, freq='M'))

# Let's break it apart
analysis = seasonal_decompose(df, model='additive')
analysis.plot()
plt.show()

• np.sin(2 * np.pi * t / 12) — generates a smooth sine wave that completes one full cycle every 12 time steps, simulating yearly seasonality on monthly data.
• pd.date_range(..., freq='M') — creates month-end timestamps for 120 periods (10 years), giving the Series a proper datetime index that seasonal_decompose needs.
• seasonal_decompose(df, model='additive') — splits the series into three additive pieces: trend, seasonal, and residual. With model='additive' it assumes the components sum up; use model='multiplicative' when the seasonal swing grows with the trend.
• analysis.plot() — renders a four-panel chart (observed, trend, seasonal, residual) so you can eyeball whether the decomposition captured the structure.

In practice, we can tell the model: look at the trend, but also remember that every 12 months, the pattern repeats.

Let’s Build It: Forecasting Ice Cream Sales

Let’s put this into practice. We will use the statsmodels library to fit a SARIMA model. The summary table is the most intimidating part of the output. Here is how to read it.

import statsmodels.api as sm

# Hold out the last 12 months so we can evaluate on data the model never saw
train, holdout = df[:-12], df[-12:]

# Fit a SARIMA model
# (1,1,1) are the ARIMA parts, (1,1,1,12) are the Seasonal parts
model = sm.tsa.statespace.SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results = model.fit(disp=False)

print(results.summary())

• train, holdout = df[:-12], df[-12:] — splits the 120-point series into the first 108 months for fitting and the last 12 for evaluation, so the model never sees the months it’s later asked to forecast. Without this split, “evaluating” the forecast against the tail of the training data isn’t a real test — see the MAPE section below for what goes wrong if you skip it.
• sm.tsa.statespace.SARIMAX — the main SARIMA entry point in statsmodels; “X” means it also accepts exogenous regressors (like holidays), though none are passed here.
• order=(1, 1, 1) — sets $p=1$ (one past value), $d=1$ (difference once), $q=1$ (one past error). These are the non-seasonal ARIMA components.
• seasonal_order=(1, 1, 1, 12) — sets seasonal $P=1$, $D=1$, $Q=1$ with period $s=12$ (yearly cycle on monthly data).
• model.fit(disp=False) — runs maximum-likelihood estimation on the 108-month training set; disp=False suppresses the convergence iteration log so the output stays clean.

When you scan that big table of numbers, focus on these values:

1. P>|z| (P-values): A number under 0.05 means that part of the model is doing useful work. Anything higher suggests unnecessary noise.
2. AIC: Think of this as an efficiency score. Between two models, the one with the lower AIC is usually better. It rewards accuracy but penalizes complexity.
3. Ljung-Box (Prob): Check the bottom of the table. You want this number above 0.05. That means your residuals look like random static. If not, the model missed a pattern in the data.

Run the fit above and look at your own table: none of the five coefficients clear the 0.05 bar here (the closest, ar.L1, sits at roughly 0.38 — still well above it). That’s not a bug, and it’s not unusual. Once you’ve differenced away a deterministic trend and a clean yearly cycle, there’s often little autocorrelation left over for the AR/MA terms to explain — they’re mopping up whatever small amount of structure survives the differencing, and on this particular synthetic series there isn’t much. The Ljung-Box row is the one that actually tells you whether that’s a problem: at Prob(Q) ≈ 0.87 here, well above 0.05, the residuals still look like clean noise, which is the real signal that the model captured the structure — even though none of the individual coefficients are statistically significant on their own.

The Verdict: When to Ditch the Neural Network

So, when should you use SARIMA instead of a deep learning model? Here’s a simple framework:

• The 1000-Row Rule: If you have fewer than 1,000 time periods (days, months), SARIMA will almost always win. Neural networks need more examples to learn.
• Interpretability: If your boss asks, “Why did the forecast drop?”, with SARIMA you can say, “Because the seasonal dip from last year started.” With an LSTM, you have to say, “The weights in layer 4 changed.”
• Speed: A seasonal SARIMA fit on about 100 data points takes well under a second on a laptop CPU — measured at roughly 0.6-0.8 seconds for the model above, across repeated runs. You can still test dozens of model configurations in the time it takes to spin up a GPU for a neural network.

When deep learning (LSTM) might still win instead:

• You have tens of thousands of rows or more — enough data for the network to genuinely learn temporal patterns rather than memorize noise.
• You have rich exogenous features (weather, promotions, competitor prices) whose interactions are non-linear and hard to hand-specify.
• Multivariate forecasting across dozens or hundreds of correlated series — LSTMs can share representations across related time series.

Bottom line: With only a few hundred rows of bakery sales, SARIMA is the right default. Reach for deep learning when the data volume grows or the feature interactions become too tangled for a statistical model to capture.

Now let’s compare the SARIMA model’s error (MAPE — Mean Absolute Percentage Error) against a naive guess — just repeating the last observed value.

forecast = results.get_forecast(steps=12).predicted_mean

mape = np.mean(np.abs((holdout - forecast) / holdout)) * 100
print(f"SARIMA Forecast Error: {mape:.2f}%")
# A result under 10% is usually considered excellent in retail forecasting.

# Compare against a naive baseline: just repeat the last training value
naive_forecast = pd.Series(train.iloc[-1], index=holdout.index)
naive_mape = np.mean(np.abs((holdout - naive_forecast) / holdout)) * 100
print(f"Naive (last-value) Forecast Error: {naive_mape:.2f}%")

Running this prints SARIMA Forecast Error: 2.00% against Naive (last-value) Forecast Error: 21.51% — the seasonal model is about 10x more accurate than just repeating December’s number for the whole next year, which is exactly the payoff SARIMA is supposed to deliver on data with a real seasonal cycle.

• results.get_forecast(steps=12) — produces a 12-step-ahead out-of-sample prediction from the SARIMA model, which was fit only on train (the first 108 months). Because the model never saw the holdout months, this forecast and that holdout now land on the same 12-month date range — a genuine test, not a comparison against training data.
• .predicted_mean — extracts just the point forecast (the expected value at each step); the full forecast object also contains confidence intervals.
• np.abs((holdout - forecast) / holdout) — computes the absolute percentage error at each step; np.mean(... ) * 100 averages and scales to a percentage. This is the MAPE formula: $\text{MAPE} = \frac{100}{n}\sum_{t} \left|\frac{y_t - \hat{y}_t}{y_t}\right|$.
• naive_forecast = pd.Series(train.iloc[-1], index=holdout.index) — the simplest possible baseline: just repeat the last training value for all 12 holdout months. Any model worth deploying needs to beat this by a wide margin.

What we learned today:

• Classical models like ARIMA are often better for small-to-medium datasets.
• ARIMA combines three pieces: Momentum (AR), Cleaning (I), and Shocks (MA).
• SARIMA adds a memory for seasonal patterns, like holidays.
• Use the AIC score to pick the best model, and P-values to check what’s working.

ARIMA won this round for Nora, but she wants to know if that’s a fluke or a general pattern. Next, we’ll test classical forecasting against machine learning on a real benchmark and find out.

Check Your Understanding

Questions below move from simple recall up to open-ended design, roughly following Bloom’s Taxonomy.

Remember What do the three letters in ARIMA stand for, and what does each one do ?

Understand In your own words, explain why regular ARIMA would misinterpret a July ice cream sales spike, and how adding the “S” (Seasonality) in SARIMA fixes that.

Apply Using the article’s differencing example (data.diff()), calculate the differenced series for the data [100, 105, 115, 130].

Analyze The article says a lower AIC score is “usually better” because it “rewards accuracy but punishes the model for being too complex.” Walk through why a model that fits the training data perfectly (very low error) could still have a worse AIC than a slightly less accurate, simpler model.

Evaluate The article’s “1000-Row Rule” recommends SARIMA over neural networks below 1,000 time periods. Critique this as a hard threshold: what data characteristics (beyond just row count) might justify choosing a neural network even with, say, only 500 rows, or justify sticking with SARIMA even with 5,000 rows?

Create Design a forecasting choice for a new scenario: a startup has 18 months of weekly signups (about 78 data points) with a clear weekly cadence and a suspected but unconfirmed yearly cycle they don’t have enough history to verify. Would you reach for SARIMA or a deep learning model, and what specific SARIMA seasonal order would you try first, given the data you have?

Related articles

• Classical Forecasting vs. Machine Learning: When Does Each Win?) — the head-to-head benchmark follow-up: does ARIMA’s small-data victory generalize?
• Why Your Forecast Fails on Christmas: Handling Multiple Seasonalities) — Nora’s flagship problem: overlapping daily, weekly, and holiday cycles that break single-season models.

References & Further reading

• Box, G. E. P., & Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control. Holden-Day. — the foundational text that formalized the ARIMA framework.
• statsmodels ARIMA documentation: https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima.model.ARIMA.html
• statsmodels SARIMAX documentation: https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html

---

Apply What You Learned

Topic: How ARIMA/SARIMA decomposes a time series into momentum (AR), stationarity-fixing differencing (I), shock-correction (MA), and seasonal memory — and how model-selection diagnostics (AIC, Ljung-Box) and accuracy metrics (MAPE) govern the forecasting pipeline on small data where deep learning overfits.

Draw a mindmap (paper, Excalidraw, Miro — anything) with at least these nodes:

• ⟨AR (Auto-Regressive)⟩ — the “yesterday matters” / momentum component, order $p$
• ⟨I (Integrated / Differencing)⟩ — the “make the data sit still” component, order $d$
• ⟨MA (Moving Average)⟩ — the “past shocks” correction component, order $q$
• ⟨SARIMA Seasonality (P, D, Q, s)⟩ — the seasonal cycle extension, period $s=12$ for yearly on monthly
• ⟨Stationarity⟩ — the constant-mean property that differencing produces
• ⟨AIC⟩ — the efficiency score that punishes complexity; lower is better
• ⟨Ljung-Box test⟩ — residual whiteness check; want Prob $> 0.05$
• ⟨MAPE⟩ — Mean Absolute Percentage Error; under 10% is excellent in retail
• ⟨hypothesis-testing⟩ — the broader framework for the $P>|z|$ coefficients and the Ljung-Box significance check

and at least these edges:

• ⟨I (Integrated / Differencing)⟩ → ⟨Stationarity⟩
• ⟨AIC⟩ → ⟨SARIMA Seasonality (P, D, Q, s)⟩
• ⟨SARIMA Seasonality (P, D, Q, s)⟩ → ⟨MAPE⟩
• ⟨Ljung-Box test⟩ → ⟨hypothesis-testing⟩

Rubric: all 9 named nodes present; the 4 required edges drawn; one extra edge of your own with a one-sentence justification of why you added it.

You’re at the stakeholder review. The VP of Operations read a blog post about LSTMs and is demanding to know why Nora’s team shipped a SARIMA model instead. Write a 200–400 word memo defending the SARIMA choice. You must ground every claim in the article’s specific evidence:

• Cite the 1000-Row Rule by its correct threshold (fewer than 1,000 time periods) and explain why Nora’s two years of daily bakery data — roughly 730 points — falls in SARIMA’s wheelhouse, not the LSTM’s.
• Name the interpretability advantage concretely: when the boss asks “why did the forecast drop?”, SARIMA can point to the seasonal component (“the dip from last year started”), whereas the LSTM answer is “the weights in layer 4 changed.”
• Reference the MAPE benchmark — the article says a result under 10% is “excellent in retail forecasting” — as the accuracy standard you’re holding the model to.
• Use at least one diagnostic from the model summary table with the correct interpretation direction: AIC (lower is better), P>|z| (under 0.05 means the component is doing useful work), or Ljung-Box Prob (above 0.05 means residuals are clean).

Deliverable: A 200–400 word memo addressed to the VP of Operations.

Rubric:

• States the 1000-Row Rule with the correct threshold (< 1,000 periods) and correctly situates ~730 daily points below it
• Names the interpretability advantage concretely, quoting or paraphrasing the “seasonal dip” vs “weights in layer 4” contrast
• Cites the under-10% MAPE benchmark as the accuracy standard
• References at least one diagnostic (AIC, P>|z|, or Ljung-Box) with the correct direction (lower / under 0.05 / above 0.05)
• Word count between 200 and 400

Ship the article’s pipeline as a reusable forecasting module. Your service wraps the three steps from the article — fit (SARIMAX), forecast (get_forecast), evaluate (MAPE on holdout) — and adds a production safety net the article doesn’t cover: a rollback hook that detects dirty residuals and falls back to a naive last-value forecast.

Spec:

1. fit_forecast(series, order, seasonal_order, steps) returns a dict with:
   • predicted_mean: the steps-ahead point forecast
   • mape: MAPE computed by fitting on series with the last steps observations held out, then forecasting forward onto that holdout (matching the article’s train, holdout = df[:-12], df[-12:] pattern — never fit on the observations you’re scoring against)
   • aic: the model’s AIC (lower is better)
   • ljung_box_prob: the Ljung-Box Prob(Q) from the diagnostic summary
   • max_coef_pvalue: the largest P>|z| among the fitted coefficients
   • used_rollback: True if ljung_box_prob < 0.05 (residuals not clean — the article says you want this above 0.05), in which case populate naive_forecast with a last-value constant forecast
2. Validate on the article’s 120-point seasonal synthetic data (order=(1,1,1), seasonal_order=(1,1,1,12), steps=12).

Starter: projects/arima-and-sarima-intuitively-when-classical-foreca/sarima_service.py

Deliverable: A runnable sarima_service.py module plus the console output from run().

Rubric:

• fit_forecast returns a dict with all six keys (predicted_mean, mape, aic, ljung_box_prob, max_coef_pvalue, used_rollback)
• The model is fit only on the data before the held-out steps observations — mape is never computed against points the model was trained on
• ljung_box_prob is extracted from the actual summary table, not hardcoded
• When ljung_box_prob < 0.05, used_rollback is True and naive_forecast is populated with a last-value constant series
• run() executes without error on the article’s 120-point seasonal data with the specified orders
• The printed MAPE is a finite float (not NaN or inf)

The article presents the ARIMA equation from Box & Jenkins (1976) as:

$\phi(B)\,(1 - B)^{d}\, y_t \;=\; \theta(B)\,\varepsilon_t$

Implement this from scratch in numpy — no statsmodels.tsa.statespace.SARIMAX for the core fitting logic. You must build each operator the article describes in plain English:

• Backshift operator $B$: $By_t = y_{t-1}$ (the “yesterday” shift)
• Differencing $(1-B)^d$: the “I” component that makes data stationary (the article’s [10, 12, 15, 18, 22] → [2, 3, 3, 4] example)
• AR polynomial $\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$: the “momentum” component
• MA polynomial $\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$: the “shocks” component

Then break it: fit AR(1) on the article’s trending data (sales = 10 + 0.5 * time + np.random.normal(0, 2, 100), seed 42) with $d=1$ (differenced) vs $d=0$ (no differencing). The article claims differencing is what makes the model see the signal instead of chasing the trend. Prove it: show that $d=0$ residuals retain the linear trend while $d=1$ residuals look like noise.

Starter: projects/arima-and-sarima-intuitively-when-classical-foreca/arima_from_scratch.py

Deliverable: A runnable arima_from_scratch.py module with the backshift, differencing, AR, and MA operators, a least-squares AR fitter, and the $d=1$ vs $d=0$ ablation comparison.

Rubric:

• backshift(series, k) returns $y_{t-k}$ aligned to the same index; first k entries are NaN
• difference(series, d) applies $(1-B)^d$ iteratively; difference([10,12,15,18,22], d=1) yields [2,3,3,4] (matching the article)
• ar_polynomial(series, phi) computes $\phi(B)\,y_t$ for given AR coefficients
• fit_ar(series, p, d) returns AR coefficients via np.linalg.lstsq on the $d$-differenced series
• Ablation: on the article’s seed-42 trending data, residual_std_d0 > residual_std_d1 — the $d=0$ residuals are larger because they retain the trend
• d0_retains_trend is True (visually or by residual-vs-time correlation)
• No statsmodels.tsa.statespace.SARIMAX used for the core fitting (you may import statsmodels only for an optional validation comparison)