ARIMA Model

An ARIMA Model (AutoRegressive Integrated Moving Average) is a parametric time-series model that explains a value using past values (AR), differencing to remove trends and achieve stationarity (I), and past forecast errors (MA). It is important for statistical forecasting and as a baseline in ML pipelines. Example: forecasting weekly demand by fitting ARIMA(p,d,q) to differenced sales data.

Imagine you’re trying to guess tomorrow’s temperature by looking at recent days: you notice a general direction (warming or cooling), and you also notice that today’s “surprise” compared to your guess often affects what happens next. An ARIMA model is a simple forecasting tool that works like this for any data measured over time (like sales, website visits, or energy use). It combines three ideas: autoregressive (using past values), integrated (removing trends by looking at changes), and moving average (using past prediction errors). It’s widely used to predict what comes next in a time series.

import numpy as np  # numerical arrays for synthetic time-series generation
import pandas as pd  # time index + tabular handling
from statsmodels.tsa.arima.model import ARIMA  # ARIMA estimation + forecasting
from statsmodels.tsa.statespace.sarimax import SARIMAX  # convenient for ARIMA with exogenous regressors

# --- Create a realistic monthly demand series with trend + seasonality + noise ---
rng = np.random.default_rng(7)  # reproducible randomness for consistent results
dates = pd.date_range("2018-01-01", periods=72, freq="MS")  # 6 years of monthly starts
t = np.arange(len(dates))  # time index for trend/seasonal construction
promo = (rng.random(len(dates)) < 0.18).astype(int)  # occasional promotion indicator (exogenous feature)
demand = 200 + 1.2 * t + 12 * np.sin(2 * np.pi * t / 12) + 18 * promo + rng.normal(0, 6, len(dates))  # observed series
y = pd.Series(demand, index=dates, name="demand")  # time-indexed series for modeling

# --- Train/test split to evaluate out-of-sample forecasting ---
y_train, y_test = y.iloc[:-12], y.iloc[-12:]  # hold out the last year
promo_train, promo_test = promo[:-12], promo[-12:]  # align exogenous feature with the split

# --- Fit ARIMA (via SARIMAX) with an exogenous regressor and forecast the holdout window ---
model = SARIMAX(y_train, order=(1, 1, 1), exog=promo_train, trend="c")  # ARIMA(p,d,q) + constant + exog
res = model.fit(disp=False)  # estimate parameters via maximum likelihood
pred = res.get_forecast(steps=len(y_test), exog=promo_test)  # *** KEY LINE *** generate out-of-sample forecasts
y_hat = pred.predicted_mean  # point forecasts aligned to forecast index
ci = pred.conf_int(alpha=0.05)  # 95% forecast interval for uncertainty quantification

# --- Quick forecast accuracy check (useful in ML workflows for model selection) ---
mae = float(np.mean(np.abs(y_test.values - y_hat.values)))  # mean absolute error on the holdout set
print(res.summary().tables[1])  # coefficients table (AR/MA terms + exog effect)
print(f"Holdout MAE: {mae:.2f}")  # scalar accuracy metric
print(ci.tail(3))  # show last few interval rows to verify shape/indexing

• The forecast is produced by get_forecast; predicted_mean gives point predictions and conf_int provides uncertainty bounds for each horizon.
• order=(1,1,1) controls the ARIMA structure; mis-setting d (differencing) is a common cause of poor forecasts or unstable parameters.
• When using exog, the exogenous array must be aligned in length and time with the target series for both training and forecasting, or you’ll get shape/index errors.
• In an AI/ML workflow, the holdout MAE is typically used to compare against alternatives (e.g., different orders, additional features, or other forecasting models) under a consistent backtest split.

set.seed(7)  # reproducible randomness for consistent results

# --- Create a realistic monthly demand series with trend + seasonality + noise ---
dates <- seq.Date(from = as.Date("2018-01-01"), by = "month", length.out = 72)  # 6 years monthly
t <- 0:(length(dates) - 1)  # time index for trend/seasonality
promo <- as.integer(runif(length(dates)) < 0.18)  # occasional promotion indicator (exogenous feature)
demand <- 200 + 1.2 * t + 12 * sin(2 * pi * t / 12) + 18 * promo + rnorm(length(dates), 0, 6)  # observed series
y <- ts(demand, start = c(2018, 1), frequency = 12)  # monthly time series object for ARIMA

# --- Train/test split to evaluate out-of-sample forecasting ---
y_train <- window(y, end = c(2022, 12))  # first 60 months (2018-01 .. 2022-12)
y_test  <- window(y, start = c(2023, 1))  # last 12 months (2023-01 .. 2023-12)
x_train <- promo[1:60]  # exogenous regressor aligned to training window
x_test  <- promo[61:72]  # exogenous regressor aligned to test window

# --- Fit ARIMA with an exogenous regressor and forecast the holdout window ---
fit <- arima(y_train, order = c(1, 1, 1), xreg = x_train, include.mean = TRUE)  # estimate ARIMA + exog
fc <- predict(fit, n.ahead = length(y_test), newxreg = x_test)  # *** KEY LINE *** out-of-sample forecasts
y_hat <- fc$pred  # point forecasts
se <- fc$se  # forecast standard errors for interval construction

# --- Quick forecast accuracy check (useful in ML workflows for model selection) ---
mae <- mean(abs(as.numeric(y_test) - as.numeric(y_hat)))  # mean absolute error on holdout
ci_low <- y_hat - 1.96 * se  # approximate 95% lower bound
ci_high <- y_hat + 1.96 * se  # approximate 95% upper bound
print(fit$coef)  # fitted coefficients (AR/MA + exogenous effect)
cat(sprintf("Holdout MAE: %.2f\n", mae))  # scalar accuracy metric
tail(cbind(point = y_hat, lo95 = ci_low, hi95 = ci_high), 3)  # sanity-check last intervals

• The forecast is produced by predict; it returns both $pred (point forecasts) and $se (standard errors), which are commonly used to form intervals.
• order=c(1,1,1) sets the ARIMA structure; if the differencing order is inappropriate, residuals often remain autocorrelated and forecast errors increase.
• With xreg/newxreg, the regressor must be provided for both training and each forecast step; mismatched lengths are a frequent source of runtime errors.
• In ML-style evaluation, the MAE on the held-out year is used for model comparison and for selecting orders/features before final refitting on all available data.