SARIMAX: Seasonal ARIMA with Exogenous Variables

• SARIMAX is a time series model that extends SARIMA by incorporating exogenous variables to capture external influences alongside seasonal patterns.
• It optimizes forecasts using maximum likelihood estimation, rolling re-estimation, and advanced techniques like MIQO for variable selection.
• Empirical case studies in network traffic, COVID-19 mortality, and energy load prediction demonstrate its enhanced accuracy over traditional models.

A Seasonal AutoRegressive Integrated Moving Average with eXogenous variables (SARIMAX) model is a flexible framework for modeling time series exhibiting both seasonal structure and dependence on external covariates. The SARIMAX generalizes SARIMA by incorporating arbitrary exogenous regressors linearly into the mean structure, allowing the model to capture variations driven by external factors beyond endogenous autocorrelation and seasonality. SARIMAX has found widespread application in domains such as network traffic forecasting, epidemiological surveillance, and energy load prediction, where both regular seasonal patterns and influential external drivers co-exist.

1. Formal Definition and Model Structure

The SARIMAX$\big((p,d,q)\times(P,D,Q)_s\big)$ model extends the classical SARIMA structure by adding an exogenous component: $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$ where:

• $L$ is the lag operator: $L\,y_t = y_{t-1}$,
• $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$ (non-seasonal AR),
• $\Phi_P(L^s) = 1-\Phi_1 L^s-\cdots-\Phi_P L^{sP}$ (seasonal AR, period $s$),
• $\Theta_q(L) = 1+\theta_1 L+\cdots+\theta_q L^q$ (non-seasonal MA),
• $$\Theta_Q(L^s) = 1+\Theta_1 L^s+\cdots+\Theta_Q L^{sQ}$$ (seasonal MA, period $s$),
• $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$0, $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$1 are non-seasonal and seasonal differencing of orders $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$2 and $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$3,
• $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$4 is the vector of exogenous regressors at time $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$5; $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$6 is their coefficient vector,
• $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$7 is a zero-mean white-noise error.

In many applications the SARIMAX can also be represented as: $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$8 where differencing is implicit when $$\Phi_p(L)\,\Phi_P(L^s)\,(1-L)^d\,(1-L^s)^D\,y_t \;=\; \Theta_q(L)\,\Theta_Q(L^s)\,\varepsilon_t \;+\; \beta^{\top} X_t$$9 or $L$0 (Saha et al., 2022, Toutiaee et al., 2021, Lowther et al., 2020, Bercu et al., 2012).

2. Role and Design of Exogenous Regressors

The exogenous covariate vector $L$1 operationalizes the dependence of $L$2 on known external signals. These regressors can be:

• Time calendar dummies: Binary indicators for time-of-day, weekday/weekend, or holiday effects (e.g., 8 time-of-day segments plus a weekend/weekday indicator as in (Saha et al., 2022)).
• Domain-driven signals: Physical measurements such as temperature for electricity demand (Bercu et al., 2012), or public health metrics (hospitalized and ICU counts) for epidemic surveillance (Toutiaee et al., 2021).
• Engineerable features: Lagged, smoothed, or transformed versions of raw covariates, possibly optimized for forecast accuracy via best-subset selection (Lowther et al., 2020).

Regression coefficients $L$3 quantify the linear effect of $L$4 on the level of $L$5 after accounting for all auto- and seasonally regressive terms.

The selection of exogenous features can be manual (driven by domain expertise), automated via grid search, or performed with explicit sparsity constraints using mixed integer quadratic optimization (MIQO) to ensure parsimony and interpretability (Lowther et al., 2020).

3. Estimation and Model Identification

Estimation of SARIMAX parameters involves determination of:

• Model orders: $L$6 and $L$7, and seasonal period $L$8.
• Regression coefficients: $L$9 governing exogenous effects.
• ARMA coefficients: for both non-seasonal and seasonal lags.

Typical workflows:

• Testing stationarity: Use the Augmented Dickey–Fuller (ADF) or KPSS tests, possibly on both $L\,y_t = y_{t-1}$0 and the residuals from regressing $L\,y_t = y_{t-1}$1 on $L\,y_t = y_{t-1}$2 (Saha et al., 2022, Bercu et al., 2012).
• Seasonality detection: Periodogram or Fourier analysis for dominant cycles (e.g., daily, weekly).
• ACF/PACF: Empirical autocorrelation/partial autocorrelation functions guide the choice of autoregressive and moving-average orders, though complex (e.g., sinusoidal) decay often requires information-criterion-based grid search.
• Automated tuning: Grid search over candidate orders with selection by Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) (Saha et al., 2022, Bercu et al., 2012, Toutiaee et al., 2021, Lowther et al., 2020).
• MIQO/Heterogeneous search: For high-dimensional or collinear $L\,y_t = y_{t-1}$3, MIQO enables simultaneous selection of regression subset and SARIMA orders by embedding variable selection directly in the likelihood optimization (Lowther et al., 2020).

Parameter estimation is typically via maximum likelihood, often implemented in state-space form using the Kalman filter (Bercu et al., 2012, Toutiaee et al., 2021), ensuring joint inference of regression and ARIMA terms.

4. Forecasting Protocols and Rolling Re-Estimation

Forecasting with SARIMAX can be performed in either static or dynamically re-estimated frameworks:

• Static fit: All parameters are fit on a training set and held fixed during forecasting on the test set (Toutiaee et al., 2021).
• Rolling or expanding window: After each new observation, the model is re-fit using all information up to the current time. This is effective for adapting to structural change, as demonstrated in traffic prediction and energy load forecasting (Saha et al., 2022, Bercu et al., 2012).

Rolling one-step-ahead forecasting involves, for each time $L\,y_t = y_{t-1}$4 in the test set:

1. Fit the SARIMAX model on all data up to $L\,y_t = y_{t-1}$5.
2. Forecast $L\,y_t = y_{t-1}$6 using the realized $L\,y_t = y_{t-1}$7.
3. Observe $L\,y_t = y_{t-1}$8, compute the forecast error, and update the training window.

This protocol is shown to significantly reduce mean absolute percentage error (MAPE) in nonstationary environments (Saha et al., 2022). For direct multi-step forecasting (e.g., $L\,y_t = y_{t-1}$9-steps ahead), exogenous regressors at times $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$0 must be known or forecasted ex ante (Toutiaee et al., 2021).

5. Empirical Performance and Case Studies

SARIMAX models have demonstrated competitive—and often superior—forecasting accuracy relative to SARIMA and other benchmarks across multiple application domains:

• Internet traffic: Incorporating time-of-day and weekend dummies as exogenous variables, SARIMAX(13,1,16)(1,0,1)[24] achieved a rolling-window MAPE of 6.83%, outperforming SARIMA (7.34%) and ARIMA (7.71%), and surpassing Holt–Winter by 0.76 percentage points (Saha et al., 2022).
• COVID-19 mortality: Adding “hospitalizedCurrently” and “inIcuCurrently” regressors, SARIMAX(4,1,4)(3,1,1)[7] achieved a 35.4% relative reduction in sMAPE over SARIMA (14.50% vs. 22.46%), and up to 64.6% reduction compared to deep graph neural network baselines (Toutiaee et al., 2021).
• Individual electric load: For thermosensitive customers, a dynamic coupled SARIMAX, separately regressing on lagged temperature and then modeling the residuals using SARIMA, minimized both absolute and relative forecast errors ($\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$1–0.24) with substantially reduced model order (Bercu et al., 2012).

A summary table of comparative performance (see (Saha et al., 2022)):

Model	Rolling MAPE	Notable exogenous variables
ARIMA	7.71%	none
SARIMA	7.34%	none
SARIMAX	6.83%	time-of-day, weekend dummies
Holt–Winter	7.59%	N/A

6. Algorithmic Advances: Subset Selection and Computational Strategies

Best-subset selection for exogenous regressors in SARIMAX models is computationally challenging due to the nonconvexity of $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$2 constraints and intricate dependence among AR, MA, and regression terms. Recent advances include:

• Mixed Integer Quadratic Optimization (MIQO): Embedding binary selection indicators $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$3 for $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$4 coefficients, enforcing cardinality constraints (e.g., $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$5) and sign/domain-specific restrictions on $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$6. This enables the discovery of sparse, interpretable SARIMAX models (Lowther et al., 2020).
• Warm starts and constraint programming: Utilizing previous solutions as initial points and structurally decomposing ARIMA order searches, greatly improving solver efficiency.
• State-space recursions and Kalman filtering: Efficient implementation for maximum likelihood estimation and real-time prediction intervals, leveraging innovations algorithms (Bercu et al., 2012, Toutiaee et al., 2021).

Practical recommendations include standardizing $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$7 and $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$8, blocking collinear variables, imposing domain-informed sign constraints, and utilizing high-performance solvers such as Gurobi or CPLEX for large $\Phi_p(L) = 1-\phi_1 L -\cdots-\phi_p L^p$9 and $\Phi_P(L^s) = 1-\Phi_1 L^s-\cdots-\Phi_P L^{sP}$0.

7. Limitations, Practical Guidelines, and Domain-Specific Insights

Limitations of SARIMAX include reliance on quality and forecastability of exogenous inputs (notably for weather-driven or health surveillance applications), and the inability to handle abrupt, non-linear shifts not anticipated by the model structure (Bercu et al., 2012, Saha et al., 2022). For some domains, domain-specific regressors may offer limited or no incremental predictive value (e.g., temperature for non-thermosensitive energy consumers) (Bercu et al., 2012).

Guidelines synthesized from empirical and methodological investigations:

• Always preprocess for stationarity and extricate seasonality before fitting SARIMAX (Saha et al., 2022, Toutiaee et al., 2021).
• Employ exogenous regressors when informed by domain knowledge or when forecast accuracy demonstrably improves.
• Select model orders and covariate subsets by penalized likelihood criteria rather than ACF/PACF alone, especially with multi-cyclic or complex autocorrelation.
• Favor rolling or dynamic re-estimation for adaptive forecasting in environments characterized by temporal drift or abrupt regime changes (Saha et al., 2022).
• Quantify predictive accuracy via scale-invariant metrics such as MAPE or sMAPE, facilitating robust comparison across models and domains.

In summary, SARIMAX provides a theoretically grounded, highly adaptable modeling framework that unifies autoregressive, seasonal, and exogenous information for univariate time series forecasting, with algorithmic innovations supporting scalable, interpretable model selection and strong empirical performance across diverse scientific and engineering applications (Saha et al., 2022, Bercu et al., 2012, Toutiaee et al., 2021, Lowther et al., 2020).