Seasonal ARIMA (SARIMA) Model

A Seasonal ARIMA (SARIMA) Model extends ARIMA to capture both short-term autocorrelation and repeating seasonal patterns using seasonal differencing and seasonal AR/MA terms. It is important for forecasting periodic time series in ML pipelines, improving accuracy over non-seasonal models. Example: predicting weekly retail demand by modeling trend plus a 52-week seasonal cycle and recent shocks.

Imagine you run an ice cream shop: sales go up every summer and down every winter, and yesterday’s sales also give clues about today’s. A Seasonal ARIMA (SARIMA) model is a way to forecast data like this when it has both “recent-history” patterns and repeating seasonal cycles.

It learns from past values (so it can follow trends and short-term ups and downs) and also from a regular calendar rhythm—like weekly, monthly, or yearly repeats. In AI and machine learning, SARIMA is often used to predict things like demand, website traffic, or energy use when seasonality matters.

import numpy as np  # numerical utilities for feature construction
import pandas as pd  # time-indexed data handling
from statsmodels.tsa.statespace.sarimax import SARIMAX  # SARIMA implementation (SARIMAX class)
from sklearn.metrics import mean_absolute_error  # forecast accuracy metric

# Create a realistic monthly series (e.g., online retail orders) with trend + yearly seasonality + noise
rng = np.random.default_rng(7)  # reproducible randomness
dates = pd.date_range("2018-01-01", periods=72, freq="MS")  # 6 years of month-start timestamps
t = np.arange(len(dates))  # time index for trend construction
y = 200 + 1.8 * t + 25 * np.sin(2 * np.pi * t / 12) + rng.normal(0, 8, size=len(t))  # seasonal signal
ts = pd.Series(y, index=dates, name="orders")  # put values on a monthly DateTimeIndex

# Hold out the last 12 months to evaluate 1-year-ahead forecasting
train, test = ts.iloc[:-12], ts.iloc[-12:]  # train/test split that respects time order

# Fit a SARIMA model: non-seasonal (p,d,q) and seasonal (P,D,Q,s) with s=12 for monthly seasonality
model = SARIMAX(
    train,
    order=(1, 1, 1),              # non-seasonal ARIMA terms
    seasonal_order=(1, 1, 1, 12),  # *** KEY LINE *** seasonal ARIMA terms with yearly seasonality
    enforce_stationarity=False,    # allow the optimizer to explore non-stationary regions
    enforce_invertibility=False    # allow non-invertible MA terms during fitting
)
res = model.fit(disp=False)  # estimate parameters via maximum likelihood

# Forecast the next 12 months and compare against the held-out data
pred = res.get_forecast(steps=12)  # multi-step forecast distribution
mean_fcst = pred.predicted_mean  # point forecasts aligned to future dates
mae = mean_absolute_error(test, mean_fcst)  # evaluate forecast error on the holdout window

print(res.summary().tables[1])  # show estimated coefficients and their uncertainty
print(f"Holdout MAE (12 months): {mae:.2f}")  # quick accuracy check
print(pd.DataFrame({"actual": test, "forecast": mean_fcst}).head())  # inspect first few predictions

• The SARIMAX fit returns estimated coefficients for both seasonal and non-seasonal terms; the forecast uses those parameters to extrapolate 12 months ahead.
• The seasonal_order=(P, D, Q, s) parameter controls seasonal dynamics; s=12 is the seasonal period for monthly data with yearly seasonality.
• A common mistake is using the wrong seasonal period (e.g., s=12 for weekly data); this typically produces poor forecasts even if the model converges.
• In ML workflows, the holdout MAE is used to compare SARIMA against alternatives (e.g., ETS/Prophet/GBMs with lag features) under a time-based split.

library(forecast)  # ARIMA/SARIMA modeling and forecasting utilities

set.seed(7)  # reproducible randomness for the simulated series

# Create a realistic monthly series (e.g., online retail orders) with trend + yearly seasonality + noise
n <- 72  # 6 years of monthly observations
t <- 0:(n - 1)  # time index for trend construction
y <- 200 + 1.8 * t + 25 * sin(2 * pi * t / 12) + rnorm(n, mean = 0, sd = 8)  # seasonal signal
ts_y <- ts(y, start = c(2018, 1), frequency = 12)  # monthly time series with annual seasonality

# Hold out the last 12 months to evaluate 1-year-ahead forecasting
train <- window(ts_y, end = c(2022, 12))  # first 60 months (2018-01 .. 2022-12)
test  <- window(ts_y, start = c(2023, 1))  # last 12 months (2023-01 .. 2023-12)

# Fit a SARIMA model by specifying both non-seasonal and seasonal orders
fit <- Arima(
  train,
  order = c(1, 1, 1),                 # non-seasonal ARIMA terms
  seasonal = list(order = c(1, 1, 1), period = 12)  # *** KEY LINE *** seasonal ARIMA terms with yearly seasonality
)

# Forecast the next 12 months and compare against the held-out data
fc <- forecast(fit, h = 12)  # produce 12-step-ahead forecasts (+ intervals)
mae <- mean(abs(as.numeric(test) - as.numeric(fc$mean)))  # compute holdout MAE

print(fit)  # show fitted coefficients and diagnostics
cat(sprintf("Holdout MAE (12 months): %.2f\n", mae))  # quick accuracy check
print(head(cbind(actual = test, forecast = fc$mean)))  # inspect first few predictions

• Arima fits the specified seasonal and non-seasonal components; forecast generates multi-step predictions (and intervals) for the next 12 months.
• The seasonal=list(order=c(P,D,Q), period=s) argument is where the seasonal structure is applied; period=12 matches monthly data with yearly seasonality.
• A common mistake is evaluating with a random split; time series models should be validated with time-respecting splits (e.g., last 12 months as a holdout).
• In AI/ML pipelines, SARIMA forecasts can be used as a baseline or as features (e.g., residuals/forecast errors) for downstream anomaly detection or hybrid models.