Time Series Analysis

Python + Google Colab · 2 hours

Core concepts for working with economic time series: what makes them different, how to check stationarity, and how to fit and forecast with ARIMA.

What is time series data?

A sequence of observations indexed by time, usually at regular intervals.

Economics: GDP · CPI · unemployment · interest rates · exchange rates Elsewhere: daily temperature · weekly hospital visits · hourly heart rate

Two features define it:

• order matters
• observations are serially dependent — the past shapes the present

Why it's different from cross-section

The classical OLS model

$y_t = \beta_0 + \beta_1 x_t + u_t$

assumes independent observations, no serial correlation, constant variance.

Time series routinely violates all three.

What we care about instead:

• persistence — how strongly the past affects the present
• trends and cycles
• forecasting

What goes wrong with naive OLS

• spurious regression — high R², no real relationship
• biased standard errors from serial correlation
• invalid inference

Regress one trending series on another and you get significance from the trends alone.

What time series analysis does instead

Models temporal dependence directly: uses the variable's own past, accounts for trend and seasonality, focuses on dynamics.

Toolkit: AR · MA · ARIMA / SARIMA · forecasting

# =========================
# 0. Setup
# =========================
# !pip install pandas_datareader

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

import statsmodels.api as sm
from pandas_datareader.data import get_data_fred
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.statespace.sarimax import SARIMAX

from sklearn.metrics import mean_absolute_error, mean_squared_error

plt.rcParams["figure.figsize"] = (10, 4)

1. Load an economic series

U.S. macro data as the running example.

# ----------------------------
# 1) Load quarterly macro data
# ----------------------------
data = sm.datasets.macrodata.load_pandas().data.copy()

# Build a quarterly DatetimeIndex from year/quarter
# macrodata has columns: year, quarter, realgdp, infl, unem, ...
data["date"] = pd.PeriodIndex(year=data["year"].astype(int),
                             quarter=data["quarter"].astype(int),
                             freq="Q").to_timestamp(how="end")
data = data.set_index("date").sort_index()

# Target series: Real GDP
y = data["realgdp"].astype(float)

# Quick look
# print(data.head())
print(data[["realgdp", "infl", "unemp", "realint"]].head())

plt.figure(figsize=(10,4))
plt.plot(y.index, y, label="Real GDP")
plt.title("Quarterly Real GDP (level)")
plt.legend()
plt.show()

/tmp/ipython-input-174184877.py:8: FutureWarning: Constructing PeriodIndex from fields is deprecated. Use PeriodIndex.from_fields instead.
  data["date"] = pd.PeriodIndex(year=data["year"].astype(int),

                                realgdp  infl  unemp  realint
date                                                         
1959-03-31 23:59:59.999999999  2710.349  0.00    5.8     0.00
1959-06-30 23:59:59.999999999  2778.801  2.34    5.1     0.74
1959-09-30 23:59:59.999999999  2775.488  2.74    5.3     1.09
1959-12-31 23:59:59.999999999  2785.204  0.27    5.6     4.06
1960-03-31 23:59:59.999999999  2847.699  2.31    5.2     1.19

First questions for any series

• Trend — does the level drift up or down over time?
• Seasonality — a pattern that repeats on the calendar?
• Outliers — points far from the rest of the series?
• Long-run cycles — swings unrelated to the calendar?
• Constant variance — or does volatility change? Any abrupt breaks?

Reading quarterly real GDP (levels)

Trend — strong and upward. Non-stationary in levels; needs differencing.

Seasonality — none visible. The series is already seasonally adjusted.

Outliers — no isolated spikes, but large shocks: early 1980s (Volcker), 2008–09 (financial crisis). These are economic events, not statistical outliers.

Cycles — business cycles ride on the growth trend: irregular expansions and recessions, not fixed-length like seasonality.

Variance — grows with the level; recessions late in the sample move far more dollars than in the 1960s. Hence log or growth rates. There are also structural breaks — level and volatility both shift around 2008–09.

The plot alone tells us naive OLS is wrong here and points toward differencing and ARIMA.

Cycle vs shock vs structural break

• Shock — a one-time hit at a point in time (the innovation $\varepsilon_t$). Effects may last, but the event doesn't repeat.
• Cycle — repeated expansions and contractions, an irregular wave around the trend.
• Structural break — the data-generating process itself changes. Mean, volatility, or relationships shift for good.

# Shock vs cycle vs structural break with simple simulated series ---

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

np.random.seed(0)
T = 200
t = np.arange(T)

# Base trend
trend = 0.05 * t

# Cycle (smooth oscillation around trend)
cycle = 2.0 * np.sin(2 * np.pi * t / 40)

# Shock: one-time negative event at t=120
shock = np.zeros(T)
shock[120] = -8.0

# Structural break: higher trend slope after t=140 (regime change)
break_t = 140
trend2 = trend.copy()
trend2[break_t:] = trend2[break_t:] + 0.08 * (t[break_t:] - break_t)

# Noise
eps = np.random.normal(0, 0.8, size=T)

# Combine into one series
ys = trend2 + cycle + shock + eps

df = pd.DataFrame({"ys": ys}, index=pd.RangeIndex(T, name="t"))

plt.figure(figsize=(11,4))
plt.plot(df.index, df["ys"], label="Series")
plt.axvline(120, linestyle="--", label="Shock (t=120)")
plt.axvline(break_t, linestyle="--", label="Structural break (t=140)")
plt.title("Toy Example: Cycle vs Shock vs Structural Break")
plt.legend()
plt.show()

Reading the toy example

• Cycles — smooth recurring swings around the trend; persistent but temporary.
• Temporary shock (t ≈ 120) — sharp move, then the series returns to its old path.
• Structural break (t ≈ 140) — permanent shift to a new level or trend.

# --- Optional: quick diagnostics (rolling mean/volatility) to 'see' regime change ---
roll = 20
rolling_mean = df["ys"].rolling(roll).mean()
rolling_std  = df["ys"].rolling(roll).std()

plt.figure(figsize=(11,4))
plt.plot(rolling_mean, label=f"Rolling mean ({roll})")
plt.axvline(120, linestyle="--")
plt.axvline(break_t, linestyle="--")
plt.title("Rolling Mean (break shows up as a level/trend change)")
plt.legend()
plt.show()

plt.figure(figsize=(11,4))
plt.plot(rolling_std, label=f"Rolling std ({roll})")
plt.axvline(120, linestyle="--")
plt.axvline(break_t, linestyle="--")
plt.title("Rolling Volatility (break can also show as variance change)")
plt.legend()
plt.show()

Stationary vs non-stationary

Stationary — statistical properties don't change over time:

$$E[x_t] = \mu, \quad Var(x_t) = \sigma^2, \quad Cov(x_t, x_{t-k}) \text{ depends only on } k$$

Examples: GDP growth, inflation, temperature anomalies — they fluctuate around a stable level.

Non-stationary — a trend, structural breaks, or growing variance. Examples: GDP level, CPI level, stock prices. These drift rather than revert.

---

Why it matters. ARIMA, SARIMAX, VAR and state-space models all assume stationarity. Fit them to non-stationary data and you get unreliable forecasts, spurious correlations, and wrong standard errors.

The standard fix turns levels into growth rates:

$$\Delta \log(X_t)$$

GDP is non-stationary. GDP growth is roughly stationary.

Trend, and the rolling mean

A rolling mean averages within a fixed window that slides along the series:

$$\text{RollingMean}_t = \frac{1}{k} \sum_{i=0}^{k-1} x_{t-i}$$

It acts as a low-pass filter — strips high-frequency noise so slow-moving components stand out.

Use it to cut short-run noise, see trends, spot medium-term cycles and structural changes. Don't use it for inference; it does not make a series stationary.

Window size — small keeps short-run wiggles, large leaves only the long run. Match it to the data frequency: 12 months for monthly, 4 quarters for quarterly, 2–8 years for business cycles.

# ----------------------------
# 2) Rolling mean (trend proxy)
# ----------------------------
w = 4  # 4 quarters = 1 year window
current_gdp_series = data["realgdp"].astype(float)
roll_mean = current_gdp_series.rolling(window=w).mean()

plt.figure(figsize=(11,4))
plt.plot(current_gdp_series, label="Real GDP (level)")
plt.plot(roll_mean, label=f"{w}-quarter rolling mean", linestyle="--")
plt.title("Real GDP and Rolling Mean (Trend Evidence)")
plt.legend()
plt.show()

What the rolling mean shows

It climbs steadily and never returns to a fixed level — GDP is non-stationary in levels.

A flat rolling mean would be consistent with stationarity. Here the average itself is moving.

Trend or unit root? Detrend or difference?

Not all trends are the same, and the fix depends on where the trend comes from.

A · Deterministic trend

$$y_t = \alpha + \beta t + u_t$$

Smooth and predictable; deviations are temporary; stationary once the trend is removed. → Include a time trend, or detrend.

B · Stochastic trend (unit root)

$$y_t = y_{t-1} + \varepsilon_t$$

with $\varepsilon_t$ white noise. Shocks are permanent, the series never reverts, variance grows. → Difference the series, or log-difference for growth rates.

---

Deciding in practice — the ADF test

• $H_0$: the series has a unit root
• fail to reject → difference the data
• reject → detrend or add a time trend

from statsmodels.tsa.stattools import adfuller

adf_stat, p_value, *_ = adfuller(y)
print(f"ADF statistic (levels): {adf_stat:.3f}")
print(f"p-value: {p_value:.3f}")

ADF statistic (levels): 1.750
p-value: 0.998

Example: GDP

• rolling mean trends strongly upward
• ADF p ≈ 0.99 → fail to reject the unit root

GDP behaves like a random walk: the trend is stochastic, not deterministic. Difference or log-difference it — a time trend alone is not enough.

# # ----------------------------
# Trend: quantify average growth
# # ----------------------------
# Taking difference
d_gdp = y.diff()  # absolute change in GDP (billions)
d_gdp.plot(title="Taking the Difference of GDP")
plt.axhline(0, linewidth=1)
plt.show()

adf_stat, p_value, *_= adfuller(d_gdp.dropna())
print(f"ADF statistic (diff): {adf_stat:.3f}")
print(f"p-value: {p_value:.3f}")

# Use log to interpret changes as approx % growth and remove the trend
logy = np.log(y)
dlogy = logy.diff()  # quarterly growth (approx)
# annualized_growth_pct = (dlogy.mean() * 4) * 100
# print(f"Average annualized growth (approx, full sample): {annualized_growth_pct:.2f}%")

dlogy.plot(title="Taking the Log-difference of GDP")
plt.axhline(0, linewidth=1)
plt.show()

adf_stat, p_value, *_= adfuller(dlogy.dropna())
print(f"ADF statistic (diff): {adf_stat:.3f}")
print(f"p-value: {p_value:.3f}")

ADF statistic (diff): -6.306
p-value: 0.000

ADF statistic (diff): -6.973
p-value: 0.000

Why did Δlog(GDP) volatility fall?

The drop after the early 1980s is the Great Moderation.

Before 1980, stop–go monetary policy, high and volatile inflation, oil shocks and wage–price spirals meant output growth swung widely. In 1979–82 Volcker raised rates hard and broke inflation; the Fed then committed to price stability and systematic, Taylor-rule-type policy. Anchored expectations stabilized demand — and growth volatility.

Note the difference: Δlog(GDP) is a percentage growth rate, which policy can stabilize. ΔGDP in dollars keeps growing simply because the economy does.

ΔGDP vs Δlog(GDP)

Both can pass ADF and still behave very differently.

d_log_gdp = dlogy        # approx. growth rate
d_log_gdp_pct = dlogy * 100       # percent units


# ----------------------------
# Compare variability early vs late
# ----------------------------
# Split sample into early half and late half
mid = len(y) // 2
idx_mid = y.index[mid]

early_dgdp = d_gdp.loc[:idx_mid].dropna()
late_dgdp  = d_gdp.loc[idx_mid:].dropna()

early_dlog = d_log_gdp_pct.loc[:idx_mid].dropna()
late_dlog  = d_log_gdp_pct.loc[idx_mid:].dropna()

print("\n--- Variability comparison (early vs late) ---")
print(f"ΔGDP std early: {early_dgdp.std():.2f} | late: {late_dgdp.std():.2f} | ratio (late/early): {(late_dgdp.std()/early_dgdp.std()):.2f}")
print(f"Δlog(GDP)% std early: {early_dlog.std():.2f} | late: {late_dlog.std():.2f} | ratio (late/early): {(late_dlog.std()/early_dlog.std()):.2f}")

# ----------------------------
# Rolling volatility (visual proof)
# ----------------------------
roll = 20  # ~5 years of quarters
vol_dgdp = d_gdp.rolling(roll).std()
vol_dlog = d_log_gdp_pct.rolling(roll).std()

plt.figure(figsize=(11,4))
plt.plot(vol_dgdp.index, vol_dgdp, label=f"Rolling std of ΔGDP ({roll}q)")
plt.title("Rolling volatility: ΔGDP variance rises with GDP level")
plt.legend()
plt.show()

plt.figure(figsize=(11,4))
plt.plot(vol_dlog.index, vol_dlog, label=f"Rolling std of Δlog(GDP)% ({roll}q)")
plt.title("Rolling volatility: Δlog(GDP) is more variance-stable")
plt.legend()
plt.show()

# ----------------------------
# 7) (Optional) Scatter: level vs absolute change
# Shows scale dependence: higher GDP -> bigger ΔGDP swings
# ----------------------------
aligned = pd.DataFrame({"GDP": y, "dGDP": d_gdp, "dlogGDP_pct": d_log_gdp_pct}).dropna()

plt.figure(figsize=(6,5))
plt.scatter(aligned["GDP"], aligned["dGDP"], s=10)
plt.title("Scale effect: higher GDP → larger ΔGDP")
plt.xlabel("GDP level")
plt.ylabel("ΔGDP")
plt.show()

plt.figure(figsize=(6,5))
plt.scatter(aligned["GDP"], aligned["dlogGDP_pct"], s=10)
plt.title("Δlog(GDP)% is less tied to GDP level")
plt.xlabel("GDP level")
plt.ylabel("Δlog(GDP) × 100")
plt.show()

print("\nInterpretation:")
print("- ADF tells you whether the series has a unit root (trend). BOTH transforms can pass ADF.")
print("- But ΔGDP is in dollars, so fluctuations grow as the economy gets larger (non-constant variance).")
print("- Δlog(GDP) is a growth rate, so movements are comparable across decades (more stable variance).")

--- Variability comparison (early vs late) ---
ΔGDP std early: 52.77 | late: 64.10 | ratio (late/early): 1.21
Δlog(GDP)% std early: 1.08 | late: 0.61 | ratio (late/early): 0.56

Interpretation:
- ADF tells you whether the series has a unit root (trend). BOTH transforms can pass ADF.
- But ΔGDP is in dollars, so fluctuations grow as the economy gets larger (non-constant variance).
- Δlog(GDP) is a growth rate, so movements are comparable across decades (more stable variance).

2. Components and seasonality

$$y_t = \text{Trend}_t + \text{Seasonal}_t + \text{Residual}_t$$

Decomposition answers: is this series mostly long-run growth, a repeating calendar pattern, or what's left over?

Residual — shocks, measurement noise, misspecification. Strong remaining patterns mean the decomposition is incomplete.

Additive $y_t = T_t + S_t + R_t$ — seasonal swings stay the same size. Multiplicative $y_t = T_t \times S_t \times R_t$ — swings grow with the level, common in economic data.

Rule of thumb: if the series grows and the seasonal peaks grow with it, go multiplicative or take logs.

Decomposition is a diagnostic — it tells you whether you need a seasonal model (SARIMA) and which transformation to use.

decomp = seasonal_decompose(logy, model='multiplicative', period=4)
decomp.plot()
plt.show()

Reading the GDP decomposition

Observed — log GDP rises smoothly: long-run growth in percentage terms.

Trend — a smooth upward curve; GDP is strongly non-stationary even in logs.

Seasonal — a 4-quarter cycle of about ±0.1%, essentially nothing. The BEA already seasonally adjusts the series.

Residual — fluctuates around zero with no pattern: business-cycle shocks, recessions, policy disturbances, noise. The big negative spikes are the early-1980s recession and 2008–09.

Airline passengers

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url)

df["Month"] = pd.to_datetime(df["Month"])
df = df.set_index("Month")
df.head()

df.plot(figsize=(10,4), title="Monthly Airline Passengers")
plt.show()

Airline passengers at a glance

• Trend — passenger counts rise strongly
• Seasonality — clear 12-month pattern, summer peaks and winter troughs
• Growing swings — seasonal amplitude scales with the level → multiplicative
• Irregular — short-term deviations around the seasonal pattern

decomp = seasonal_decompose(df["Passengers"], model="multiplicative")
decomp.plot()
plt.show()

## If we use additive decomposition, the size of residuals grows over time, indicating heteroskedasticity.
## Because seasonal swings and residuals grow with the level of the series, the additive model is **misspecified**.

Reading the multiplicative decomposition

$$\text{Passengers}_t = \text{Trend}_t \times \text{Seasonal}_t \times \text{Residual}_t$$

Trend — climbs from ~120 to over 450 as incomes rise and flying becomes routine.

Seasonal — oscillates around 1, so these are percentage effects: 1.2 means 20% above normal, 0.8 means 20% below. Summers are consistently high, winters low — peak demand is proportional to the level.

Residual — stays near 1 with no systematic pattern and no growing variance. The multiplicative structure was the right call.

Under an additive model the seasonal swings would grow over time and leak into the residuals.

Takeaway: model this with a log-transformed SARIMA, not an additive seasonal model.

result = adfuller(df["Passengers"])
print("ADF statistic:", result[0])
print("p-value:", result[1])

ADF statistic: 0.8153688792060498
p-value: 0.991880243437641

df["diff"] = df["Passengers"].diff()
df["diff"].plot(figsize=(10,4), title="Differenced series")
plt.show()

result = adfuller(df["diff"].dropna())
print("ADF statistic:", result[0])
print("p-value:", result[1])

ADF statistic: -2.8292668241700047
p-value: 0.05421329028382478

p = df["Passengers"].astype(float)

# Log transform (multiplicative growth & seasonality)
p_log = np.log(p)

def adf_report(series, name):
    series = series.dropna()
    stat, pval, _, _, crit, _ = adfuller(series)
    print(f"ADF for {name}: stat={stat:.3f}, p-value={pval:.4f}")

d1 = p_log.diff(1)           # regular difference
D12 = p_log.diff(12)         # seasonal difference
d1D12 = p_log.diff(12).diff(1)
D12 = p_log.diff(12)


d1D12.plot(figsize=(10,4), title="Differenced series")
plt.show()

# adf_report(p, "p")
# adf_report(df["diff"].dropna(), "Δ p")
adf_report(p_log, "log(p)")
adf_report(d1, "Δ log(p)")
adf_report(D12, "Δ12 log(p)")
adf_report(d1D12, "Δ Δ12 log(p)")

ADF for log(p): stat=-1.717, p-value=0.4224
ADF for Δ log(p): stat=-2.717, p-value=0.0711
ADF for Δ12 log(p): stat=-2.710, p-value=0.0724
ADF for Δ Δ12 log(p): stat=-4.443, p-value=0.0002

The ARIMA model

$$\text{ARIMA}(p, d, q)$$

AR past values · I differencing · MA past shocks

---

AR(p) — today depends on the past

$$x_t = \phi_1 x_{t-1} + \cdots + \phi_p x_{t-p} + \varepsilon_t$$

Captures persistence. If airline demand was strong last month, it is probably still strong.

I(d) — differencing

$$y_t = \Delta^d x_t$$

Removes unit roots so the series fluctuates around a constant mean. $d=1$ is a first difference; quarterly data may need seasonal differencing.

MA(q) — today depends on past surprises

$$x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$$

$\varepsilon_t$ is a white-noise shock. A weather event or fuel-price spike moves demand for a few months, then fades.

---

Together: $\Delta^d x_t = \text{AR}(p) + \text{MA}(q) + \text{noise}$

Use ARIMA when the series is univariate, stationary after differencing, and driven mainly by its own past.

Why AR parameters must be below 1

$$x_t = \phi x_{t-1} + \varepsilon_t$$

$\phi$	behavior
0	no memory — pure noise
0 to 1	shocks fade — stable, mean-reverting
≈ 1	shocks persist — looks like a trend
1	random walk (unit root)

Why $|\phi| < 1$ gives stationarity. Iterating:

$$x_t = \phi^k x_{t-k} + \sum_{j=0}^{k-1} \phi^j \varepsilon_{t-j}$$

If $|\phi| < 1$ then $\phi^k \to 0$: old values die out, leaving a constant mean, finite variance, stable correlations.

At $\phi = 1$ shocks never decay, variance grows, the series drifts — a unit root. This is what GDP and CPI look like in levels, and why we difference:

$$\Delta x_t = x_t - x_{t-1}$$

That is the I(d) in ARIMA.

Why MA parameters also stay below 1

$$x_t = \varepsilon_t + \theta \varepsilon_{t-1}$$

An MA process is always stationary, so the restriction isn't about stationarity — it's invertibility: being able to recover the shocks $\varepsilon_t$ uniquely from the observed data.

With $|\theta| \ge 1$, different values of $\theta$ generate the same series, so estimation and forecasting are ambiguous.

AR: $|\phi| < 1$ for stability. MA: $|\theta| < 1$ for identification.

Invertibility, concretely

$$x_t = \varepsilon_t + \theta \varepsilon_{t-1}$$

The shocks are unobserved; we see only $x_t$. Solving backward:

$$\varepsilon_t = x_t - \theta x_{t-1} + \theta^2 x_{t-2} - \theta^3 x_{t-3} + \cdots$$

This converges only if $|\theta| < 1$ — and convergence is exactly what lets us recover the shocks from the data.

If $|\theta| \ge 1$ the series diverges and different shock sequences produce identical data: the model is not identified.

Choosing p and q with ACF and PACF

Once the series is stationary we still need the AR order p and the MA order q.

ACF

$$\text{ACF}(k) = \text{Corr}(x_t, x_{t-k})$$

This is total correlation, direct and indirect. If $x_t$ depends on $x_{t-1}$ and $x_{t-1}$ on $x_{t-2}$, the ACF at lag 2 is large even with no direct link.

Why the ACF gives q

In an MA(q), today depends only on the last $q$ shocks. For $k > q$ the shocks behind $x_t$ and $x_{t-k}$ don't overlap, and since shocks are independent:

$$\text{Corr}(x_t, x_{t-k}) = 0 \quad \text{for } k > q$$

So the ACF is nonzero through lag $q$ and cuts off sharply after it.

plot_acf(d1D12.dropna(), lags=24)
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipython-input-3514529897.py in <cell line: 0>()
----> 1 plot_acf(d1D12.dropna(), lags=24)
      2 plt.show()

NameError: name 'plot_acf' is not defined

Reading the ACF of ΔΔ₁₂ log(Passengers)

$$z_t = \Delta \Delta_{12}\log(\text{Passengers})$$

Differencing removed the trend ($\Delta$) and the seasonality ($\Delta_{12}$), so this shows what dependence is left.

• Lag 1, large negative spike → a positive surprise this month is partly reversed next month: $q = 1$
• Lag 12, large negative spike → observations a year apart are still correlated: seasonal $Q = 1$ with $s = 12$
• Everything else inside the band — no higher orders needed

With $d=1$ and $D=1$, this points to $\text{SARIMA}(p,1,1)\times(P,1,1)_{12}$. The PACF settles $p$ and $P$.

The PACF

The direct relationship between $x_t$ and $x_{t-k}$, with all intermediate lags removed:

$$\text{PACF}(k) = \text{Corr}\!\left(x_t, x_{t-k} \mid x_{t-1}, \dots, x_{t-k+1}\right)$$

How it's computed. Run

$$x_t = \phi_1 x_{t-1} + \cdots + \phi_k x_{t-k} + u_t$$

and take $\text{PACF}(k) = \hat{\phi}_k$, the coefficient on the last lag. It answers: once lags 1 to $k-1$ are accounted for, does lag $k$ still matter?

That is what the ACF cannot tell us — it shows chained correlation as if it were direct.

A PACF that cuts off after lag $p$ suggests AR($p$).

plot_pacf(d1D12.dropna(), lags=24)
plt.show()

Reading the PACF of ΔΔ₁₂ log(Passengers)

• Lag 1 significant → $z_t$ depends directly on $z_{t-1}$: $p = 1$
• Lag 12 significant → direct link one year apart: seasonal $P = 1$, $s = 12$
• Other lags inside the band → no AR(2), AR(3), or seasonal AR(2)

Combined with the ACF ($q = 1$, $Q = 1$), the natural candidate is

$$\boxed{\text{SARIMA}(1,1,1)\times(1,1,1)_{12}}$$

from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(
    p_log,
    order=(1,1,1),
    seasonal_order=(1,1,1,12),
    enforce_stationarity=False,
    enforce_invertibility=False
)

results = model.fit()
print(results.summary())

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipython-input-3123595262.py in <cell line: 0>()
      2 
      3 model = SARIMAX(
----> 4     p_log,
      5     order=(1,1,1),
      6     seasonal_order=(1,1,1,12),

NameError: name 'p_log' is not defined

SARIMAX results

$$SARIMAX(1,1,1)\times(1,1,1)_{12}$$

on monthly passenger data, 1949–1960. With $d=1$ and $D=1$, the model explains growth and de-seasonalized demand.

term	estimate	p	reading
MA(1)	−0.55	0.002	shocks are strongly corrected next period
Seasonal AR(12)	−0.33	0.043	same month last year matters — persistent seasonality
AR(1)	0.10	0.62	no short-run momentum
Seasonal MA(12)	−0.24	0.20	weak, not distinguishable from zero

Diagnostics — Ljung–Box p = 0.99 (no leftover autocorrelation) · Jarque–Bera p = 0.23 (roughly normal) · heteroskedasticity p ≈ 0.05 (mild). Residuals behave like white noise.

Implication: demand is driven by shock correction and seasonal persistence, not autoregressive momentum. The simpler $SARIMA(0,1,1)\times(1,1,0)_{12}$ should do about as well.

from statsmodels.stats.diagnostic import acorr_ljungbox


def fit_sarima(order, seasonal_order):
    model = SARIMAX(
        p_log,
        order=order,
        seasonal_order=seasonal_order,
        enforce_stationarity=False,
        enforce_invertibility=False
    )
    res = model.fit(disp=False)
    return res

# Small grid around your guessed structure
p_vals = [0,1,2]
q_vals = [0,1,2]
P_vals = [0,1]
Q_vals = [0,1]

d = 1
D = 1
s = 12

rows = []
for p in p_vals:
    for q in q_vals:
        for P in P_vals:
            for Q in Q_vals:
                order = (p,d,q)
                seasonal_order = (P,D,Q,s)
                try:
                    res = fit_sarima(order, seasonal_order)

                    # Ljung-Box on residuals (we want p-value > 0.05)
                    lb_p = acorr_ljungbox(res.resid.dropna(), lags=[12], return_df=True)["lb_pvalue"].iloc[0]

                    rows.append({
                        "order": order,
                        "seasonal_order": seasonal_order,
                        "AIC": res.aic,
                        "BIC": res.bic,
                        "LjungBox_p(12)": lb_p
                    })
                except Exception as e:
                    pass

results_df = pd.DataFrame(rows).sort_values("AIC").reset_index(drop=True)
results_df.head(10)

/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

	order	seasonal_order	AIC	BIC	LjungBox_p(12)
0	(1, 1, 0)	(0, 1, 0, 12)	-445.413971	-439.678902	0.000521
1	(0, 1, 1)	(0, 1, 0, 12)	-442.038312	-436.318687	0.000432
2	(1, 1, 1)	(0, 1, 0, 12)	-440.303889	-431.724451	0.000393
3	(2, 1, 2)	(0, 1, 0, 12)	-439.411104	-425.150953	0.000535
4	(2, 1, 0)	(0, 1, 0, 12)	-439.089395	-430.509958	0.000451
5	(1, 1, 0)	(0, 1, 1, 12)	-437.525820	-429.213766	0.000146
6	(0, 1, 1)	(1, 1, 0, 12)	-437.115916	-428.778545	0.028315
7	(0, 1, 2)	(0, 1, 0, 12)	-435.808965	-427.252875	0.000434
8	(2, 1, 0)	(0, 1, 1, 12)	-435.648119	-424.565380	0.000126
9	(0, 1, 1)	(0, 1, 1, 12)	-435.443521	-427.156999	0.000470

Selecting the final model

Two tools, both required:

1. AIC / BIC — fit
2. Ljung–Box — are the residuals white noise?

$H_0$: residuals are independent. p > 0.05 means the model captured the structure; p < 0.05 means leftover autocorrelation.

A model is invalid if Ljung–Box rejects it, however low the AIC.

In our runs, some low-AIC models had Ljung–Box p ≈ 0.000 — over-simplified, missing seasonal dynamics. Models with seasonal AR or MA terms had slightly higher AIC but passed.

Rule: lowest AIC among the models that pass Ljung–Box.

Final model

$$\boxed{\text{SARIMA}(0,1,1)\times(1,1,0)_{12}}$$

Lowest AIC among the models that pass Ljung–Box.

Forecasting with SARIMA

1. Split — train on the first 80%, hold out the last 20%
2. Fit on the training set, using the log level series; the model differences internally via $(d, D)$
3. Forecast the test window with confidence intervals
4. Back-transform — we modeled logs, so $\widehat{\text{Passengers}}_t = \exp(\widehat{\log \text{Passengers}}_t)$
5. Plot actual vs forecast with the interval band
6. Evaluate — RMSE, MAE, MAPE on the test period

A good forecast tracks both the trend and the seasonal swings, and most test points fall inside the band.

# ============================================================
# ✈️ Airline Passengers (LOG) Workshop Template
# Variable name: p
# Includes: train/test split, ACF/PACF, ARIMA, SARIMA, forecasting + plots
# ============================================================

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

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX

from sklearn.metrics import mean_absolute_error, mean_squared_error

# ----------------------------
# 1) Load Airline Passengers
# ----------------------------
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url)
df["Month"] = pd.to_datetime(df["Month"])
df = df.set_index("Month")

p = df["Passengers"].astype(float)

# Log transform (multiplicative growth & seasonality)
p_log = np.log(p)

# ----------------------------
# 2) Train/test split (80/20)
# ----------------------------
n = len(p_log)
train_size = int(0.8 * n)

train = p_log.iloc[:train_size]
test  = p_log.iloc[train_size:]

print("Train:", train.index.min().date(), "→", train.index.max().date(), f"(n={len(train)})")
print("Test: ", test.index.min().date(),  "→", test.index.max().date(),  f"(n={len(test)})")

plt.figure(figsize=(11,4))
plt.plot(train, label="Train (log passengers)")
plt.plot(test, label="Test (log passengers)", color="black")
plt.title("Airline Passengers (log) with 80/20 split")
plt.legend()
plt.show()

# ----------------------------
# 3) Stationarity checks + differencing helpers
# ----------------------------
def adf_report(series, name):
    series = series.dropna()
    stat, pval, _, _, crit, _ = adfuller(series)
    print(f"ADF for {name}: stat={stat:.3f}, p-value={pval:.4f}")

d1_train = train.diff(1)           # regular difference
D12_train = train.diff(12)         # seasonal difference
d1D12_train = train.diff(12).diff(1)

adf_report(train, "log(p)")
adf_report(d1_train, "Δ log(p)")
adf_report(D12_train, "Δ12 log(p)")
adf_report(d1D12_train, "Δ Δ12 log(p)")

# ----------------------------
# 4) ACF / PACF (on ΔΔ12 log(p))
# ----------------------------
series_for_acf = d1D12_train.dropna()

fig, axes = plt.subplots(1, 2, figsize=(12,4))
plot_acf(series_for_acf, lags=36, ax=axes[0])
axes[0].set_title("ACF of ΔΔ12 log(p) [train]")
plot_pacf(series_for_acf, lags=36, ax=axes[1])
axes[1].set_title("PACF of ΔΔ12 log(p) [train]")
plt.tight_layout()
plt.show()

# ============================================================
# 5) ARIMA on log(p) (non-seasonal baseline)
# ============================================================
arima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(0, 0, 0, 0),
    enforce_stationarity=False,
    enforce_invertibility=False
)
arima_fit = arima_model.fit(disp=False)
print(arima_fit.summary())

arima_fc = arima_fit.get_forecast(steps=len(test))
arima_pred = arima_fc.predicted_mean
arima_ci = arima_fc.conf_int()

arima_pred.index = test.index
arima_ci.index = test.index

plt.figure(figsize=(11,5))
plt.plot(train, label="Train (log)")
plt.plot(test, label="Test (log)", color="black")
plt.plot(arima_pred, label="ARIMA forecast (log)", linestyle="--")
plt.fill_between(arima_ci.index, arima_ci.iloc[:,0], arima_ci.iloc[:,1], alpha=0.3)
plt.title("ARIMA(1,1,1) on log(p): Forecast")
plt.legend()
plt.show()

mae_arima = mean_absolute_error(test, arima_pred)
rmse_arima = np.sqrt(mean_squared_error(test, arima_pred))
print(f"ARIMA on log scale: MAE={mae_arima:.4f}, RMSE={rmse_arima:.4f}")

# Back-transform
arima_pred_level = np.exp(arima_pred)
test_level = np.exp(test)
train_level = np.exp(train)

plt.figure(figsize=(11,5))
plt.plot(train_level, label="Train (level)")
plt.plot(test_level, label="Test (level)", color="black")
plt.plot(arima_pred_level, label="ARIMA forecast (level)", linestyle="--")
plt.title("ARIMA Forecast (levels)")
plt.legend()
plt.show()

# ============================================================
# 6) SARIMA on log(p)  (seasonal model)
# Classic: (0,1,1)(1,1,0,12)
# ============================================================
sarima_model = SARIMAX(
    train,
    order=(0, 1, 1),
    seasonal_order=(1, 1, 0, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
)
sarima_fit = sarima_model.fit(disp=False)
print(sarima_fit.summary())

sarima_fc = sarima_fit.get_forecast(steps=len(test))
sarima_pred = sarima_fc.predicted_mean
sarima_ci = sarima_fc.conf_int()

sarima_pred.index = test.index
sarima_ci.index = test.index

plt.figure(figsize=(11,5))
plt.plot(train, label="Train (log)")
plt.plot(test, label="Test (log)", color="black")
plt.plot(sarima_pred, label="SARIMA forecast (log)", linestyle="--")
plt.fill_between(sarima_ci.index, sarima_ci.iloc[:,0], sarima_ci.iloc[:,1], alpha=0.3)
plt.title("SARIMA(0,1,1)(0,1,1,12) on log(p): Forecast")
plt.legend()
plt.show()

mae_sarima = mean_absolute_error(test, sarima_pred)
rmse_sarima = np.sqrt(mean_squared_error(test, sarima_pred))
print(f"SARIMA on log scale: MAE={mae_sarima:.4f}, RMSE={rmse_sarima:.4f}")

# Back-transform
sarima_pred_level = np.exp(sarima_pred)

plt.figure(figsize=(11,5))
plt.plot(train_level, label="Train (level)")
plt.plot(test_level, label="Test (level)", color="black")
plt.plot(sarima_pred_level, label="SARIMA forecast (level)", linestyle="--")
plt.title("SARIMA Forecast (levels)")
plt.legend()
plt.show()

# ============================================================
# 7) Compare ARIMA vs SARIMA (levels)
# ============================================================
plt.figure(figsize=(11,5))
plt.plot(test_level, label="Test (level)", color="black")
plt.plot(arima_pred_level, label="ARIMA forecast", linestyle="--")
plt.plot(sarima_pred_level, label="SARIMA forecast", linestyle="--")
plt.title("Forecast Comparison (Levels)")
plt.legend()
plt.show()

Train: 1949-01-01 → 1958-07-01 (n=115)
Test:  1958-08-01 → 1960-12-01 (n=29)

ADF for log(p): stat=-1.574, p-value=0.4966
ADF for Δ log(p): stat=-2.636, p-value=0.0858
ADF for Δ12 log(p): stat=-1.991, p-value=0.2905
ADF for Δ Δ12 log(p): stat=-3.946, p-value=0.0017

/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

                               SARIMAX Results                                
==============================================================================
Dep. Variable:             Passengers   No. Observations:                  115
Model:               SARIMAX(1, 1, 1)   Log Likelihood                  99.320
Date:                Thu, 15 Jan 2026   AIC                           -192.640
Time:                        10:35:47   BIC                           -184.485
Sample:                    01-01-1949   HQIC                          -189.331
                         - 07-01-1958                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.5757      0.181     -3.173      0.002      -0.931      -0.220
ma.L1          0.8426      0.109      7.743      0.000       0.629       1.056
sigma2         0.0099      0.002      5.322      0.000       0.006       0.014
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                 5.62
Prob(Q):                              0.89   Prob(JB):                         0.06
Heteroskedasticity (H):               1.00   Skew:                             0.16
Prob(H) (two-sided):                  0.99   Kurtosis:                         1.95
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

ARIMA on log scale: MAE=0.1874, RMSE=0.2198

/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/usr/local/lib/python3.12/dist-packages/statsmodels/tsa/base/tsa_model.py:473: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

                                      SARIMAX Results                                      
===========================================================================================
Dep. Variable:                          Passengers   No. Observations:                  115
Model:             SARIMAX(0, 1, 1)x(1, 1, [], 12)   Log Likelihood                 165.656
Date:                             Thu, 15 Jan 2026   AIC                           -325.312
Time:                                     10:35:47   BIC                           -317.813
Sample:                                 01-01-1949   HQIC                          -322.288
                                      - 07-01-1958                                         
Covariance Type:                               opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.L1         -0.4206      0.094     -4.466      0.000      -0.605      -0.236
ar.S.L12      -0.4552      0.083     -5.495      0.000      -0.618      -0.293
sigma2         0.0015      0.000      7.249      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   0.01   Jarque-Bera (JB):                 2.04
Prob(Q):                              0.91   Prob(JB):                         0.36
Heteroskedasticity (H):               0.27   Skew:                             0.28
Prob(H) (two-sided):                  0.00   Kurtosis:                         3.48
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

SARIMA on log scale: MAE=0.0526, RMSE=0.0631

Beyond ARIMA

VAR — several macro variables moving together (GDP, inflation, unemployment); each depends on its own lags and the others'. For macro forecasting, impulse responses, policy analysis.

State-space and the Kalman filter — unobserved components (trend/cycle), missing data, dynamic factor models. For nowcasting, mixed-frequency data, time-varying dynamics.

Nonlinear and ML — gradient-boosted trees on lag features, LSTM/GRU/Transformers, or ARIMA plus ML on the residuals. For demand forecasting, web traffic, sensor data.

ARIMA is the baseline. Reach past it when the data has trend changes, several interacting series, time-varying volatility, structural breaks, or strong nonlinearity.