ARIMAX Model

An ARIMAX Model extends ARIMA by adding exogenous (external) predictors to explain a time series while still modeling autocorrelation and differencing for nonstationarity. It is important in AI/ML forecasting when outcomes depend on both past values and known drivers, improving accuracy and interpretability. Example: forecasting electricity demand using past demand plus temperature and calendar indicators.

Imagine you’re trying to predict tomorrow’s ice cream sales. You’d look at past sales (hot days usually mean more sales), but you might also want to include an extra clue like tomorrow’s temperature forecast. An ARIMAX model does exactly this for data that changes over time: it predicts future values using patterns in the past (trends, repeating cycles, and “momentum”) and also uses outside influences called exogenous variables (extra inputs), like weather, holidays, or ad spend. It’s a way to make time-based predictions smarter by adding relevant real-world factors.

import numpy as np  # numerical work (simulate a realistic time series)
import pandas as pd  # time-indexed data handling
from statsmodels.tsa.statespace.sarimax import SARIMAX  # ARIMAX via SARIMAX with exogenous regressors
from sklearn.metrics import mean_absolute_error  # simple forecast accuracy metric

# --- Create a realistic monthly demand series with a promo (exogenous) driver ---
rng = np.random.default_rng(7)  # reproducible randomness
dates = pd.date_range("2018-01-01", periods=72, freq="MS")  # 6 years of monthly data
promo = (rng.random(len(dates)) < 0.25).astype(int)  # 25% of months have promotions (0/1)
noise = rng.normal(0, 3, size=len(dates))  # idiosyncratic demand shocks

# Build demand with trend + seasonality + promo lift + autocorrelated structure
base = 120 + 0.6 * np.arange(len(dates))  # slow upward trend
season = 8 * np.sin(2 * np.pi * np.arange(len(dates)) / 12)  # annual seasonality
demand = base + season + 10 * promo + noise  # promo increases demand on average

df = pd.DataFrame({"demand": demand, "promo": promo}, index=dates)  # assemble time series table

# --- Train/test split to evaluate out-of-sample forecasting ---
train = df.iloc[:-12]  # keep last 12 months for testing
test = df.iloc[-12:]  # holdout period

# Fit ARIMAX: ARIMA(p,d,q) with exogenous regressor(s) in X (here: promo)
# *** KEY LINE ***
model = SARIMAX(train["demand"], exog=train[["promo"]], order=(1, 1, 1), enforce_stationarity=False)
res = model.fit(disp=False)  # estimate parameters via maximum likelihood

# Forecast the holdout period using known/assumed future promo values
# *** KEY LINE ***
fc = res.get_forecast(steps=len(test), exog=test[["promo"]])  # exog must align with forecast horizon
pred = fc.predicted_mean  # point forecasts
mae = mean_absolute_error(test["demand"], pred)  # evaluate forecast error

print(res.summary().tables[1])  # show coefficient table (incl. promo effect)
print(f"Holdout MAE: {mae:.2f}")  # report forecast accuracy

• The key output is the fitted parameter table from SARIMAX, including the exogenous coefficient (promo lift) and ARIMA terms; forecasts come from get_forecast.
• In SARIMAX, ARIMAX is implemented by passing exog=...; the forecast call must provide future exogenous values with the same column structure.
• Common mistake: forgetting to pass exog at forecast time (or passing misaligned indices/shapes), which causes errors or incorrect forecasts.
• In ML workflows, ARIMAX-style models are often used as strong baselines and for feature-attribution-like interpretation of known drivers (e.g., promos, holidays) alongside autocorrelation.

library(forecast)  # ARIMA with exogenous regressors via Arima()
library(ggplot2)   # quick visualization if needed

# --- Create a realistic monthly demand series with a promo (exogenous) driver ---
set.seed(7)  # reproducibility
n <- 72  # 6 years of monthly observations
dates <- seq.Date(as.Date("2018-01-01"), by = "month", length.out = n)  # monthly index
promo <- as.integer(runif(n) < 0.25)  # 25% promo months (0/1)
noise <- rnorm(n, mean = 0, sd = 3)  # demand shocks

# Build demand with trend + seasonality + promo lift
base <- 120 + 0.6 * (0:(n - 1))  # upward trend
season <- 8 * sin(2 * pi * (0:(n - 1)) / 12)  # annual seasonality
demand <- base + season + 10 * promo + noise  # promo increases demand

y <- ts(demand, frequency = 12, start = c(2018, 1))  # monthly time series object
x <- promo  # exogenous regressor (must align 1:1 with y)

# --- Train/test split (last 12 months as holdout) ---
y_train <- window(y, end = c(2022, 12))  # first 60 months
y_test  <- window(y, start = c(2023, 1))  # last 12 months
x_train <- x[1:length(y_train)]  # exog for training window
x_test  <- x[(length(y_train) + 1):length(y)]  # exog for forecast horizon

# Fit ARIMAX: ARIMA(p,d,q) with xreg (here: promo)
# *** KEY LINE ***
fit <- Arima(y_train, order = c(1, 1, 1), xreg = x_train, include.drift = TRUE)

# Forecast using future promo values (xreg must be provided for h steps)
# *** KEY LINE ***
fc <- forecast(fit, h = length(y_test), xreg = x_test)

mae <- mean(abs(as.numeric(y_test) - as.numeric(fc$mean)))  # simple accuracy metric
print(fit)  # shows coefficients incl. xreg effect
cat(sprintf("Holdout MAE: %.2f\n", mae))  # report error

• The fitted object from Arima includes the exogenous coefficient (xreg) and ARIMA terms; the forecasted values are in fc$mean from forecast.
• forecast::forecast requires xreg for the full horizon h; its length must match h exactly.
• Common mistake: using an xreg vector that is off by one due to windowing, which silently misaligns the driver with the target and degrades forecasts.
• In AI/ML pipelines, the MAE computed on a time-based holdout is a quick sanity check before moving to more complex models or adding richer exogenous features.