API Reference / algorithms / timeseries /

sarimax.py

Seasonal ARIMA with eXogenous regressors, estimated by Kalman-filter MLE.

Classes

SARIMAX

class algorithms.timeseries.sarimax.SARIMAX(Regressor)

Seasonal AutoRegressive Integrated Moving Average model with eXogenous regressors, estimated by exact Gaussian maximum likelihood through a Kalman filter.

SARIMAX is the full (p, d, q) \times (P, D, Q, m) specification plus a regression on external covariates. Unlike a hand-rolled conditional-least-squares ARIMA, the model is written in state-space form and the exact likelihood is evaluated by the Kalman filter, so the first observations are used rather than discarded, and the optimiser is kept inside the stationary/invertible region by reparameterising every AR and MA block through partial autocorrelations (Monahan/Jones transform).
Constructor 
__init__(
    self,
    order: Tuple[int, int, int] = (),
    seasonal_order: Tuple[int, int, int, int] = (),
    trend: Optional[str] = None,
    maxiter: int = 50,
    tol: float = 1e-06,
    enforce_stationarity: bool = True,
    enforce_invertibility: bool = True,
)

Overview

The estimation procedure runs as follows:

  1. Difference the series d times at lag 1 and D times at
lag m, applying the same operators to the exogenous design.
  1. Expand the multiplicative seasonal polynomials
\phi(L)\Phi(L^m) and \theta(L)\Theta(L^m) into flat lag polynomials.
  1. Place the resulting ARMA process in Harvey's companion state-space
form with a single innovation.
  1. Initialise the state covariance from the stationary solution of the
discrete Lyapunov equation P = T P T' + R R'.
  1. Run the Kalman filter, concentrating \sigma^2 out of the
likelihood, and maximise the resulting profile log-likelihood over the transformed ARMA parameters and the regression coefficients.
  1. Forecast by iterating the state transition forward, then invert the
differencing operators to return to the original scale.

Theory

The model is

\phi(L)\, \Phi(L^{m})\, (1 - L)^{d} (1 - L^{m})^{D} \left( y_t - \beta^{\top} x_t \right) = \theta(L)\, \Theta(L^{m})\, \varepsilon_t , \qquad \varepsilon_t \sim N(0, \sigma^2).

Writing the differenced, regression-adjusted series as w_t, the state-space representation is

\begin{aligned} \alpha_{t+1} &= T \alpha_t + R \varepsilon_t \\ w_t &= Z \alpha_t \end{aligned}

with T the companion matrix of the expanded AR polynomial, R = (1, \theta_1, \dots, \theta_{r-1})^{\top} and Z = (1, 0, \dots, 0). The Kalman recursions give the one-step-ahead innovations v_t and their variances \sigma^2 F_t, and the exact log-likelihood is

\log L = -\frac{n}{2}\log(2\pi\sigma^2) -\frac{1}{2}\sum_{t=1}^{n}\log F_t -\frac{1}{2\sigma^2}\sum_{t=1}^{n} \frac{v_t^2}{F_t}.

Concentrating out \sigma^2 yields \hat{\sigma}^2 = n^{-1}\sum_t v_t^2 / F_t, so only the ARMA and regression parameters remain to be optimised numerically.

Stationarity is enforced structurally. Each AR block is parameterised by partial autocorrelations r_k = \tanh(u_k) \in (-1, 1) which the Levinson-Durbin recursion maps to coefficients whose polynomial has all roots outside the unit circle, so the optimiser physically cannot reach an explosive region.

Parameters

order
tuple of (int, int, int), 0, 0) = (1
The non-seasonal :math:`(p, d, q)` order.
seasonal_order
tuple of (int, int, int, int), 0, 0, 0) = (0
The seasonal :math:`(P, D, Q, m)` order. A period m of 0 or 1 disables the seasonal component.
trend
{"c", "t", "ct", None} = None
Deterministic terms added to the differenced series: "c" a constant, "t" a linear time index, "ct" both.
maxiter
int = 50
Maximum number of optimiser iterations. Kept small by default so that a fit on a short series stays well under a second.
tol
float = 1e-6
Convergence tolerance passed to the optimiser.
enforce_stationarity
bool = True
If True, autoregressive blocks are reparameterised through partial autocorrelations. If False, the raw coefficients are optimised directly (faster, but the optimiser may wander).
enforce_invertibility
bool = True
Same as above for the moving-average blocks.

Attributes

ar_params_
np.ndarray of shape (p,)
Fitted non-seasonal autoregressive coefficients.
ma_params_
np.ndarray of shape (q,)
Fitted non-seasonal moving-average coefficients.
seasonal_ar_params_
np.ndarray of shape (P,)
Fitted seasonal autoregressive coefficients.
seasonal_ma_params_
np.ndarray of shape (Q,)
Fitted seasonal moving-average coefficients.
exog_params_
np.ndarray of shape (k_exog,)
Regression coefficients on the exogenous columns, in input order.
trend_params_
np.ndarray
Coefficients of the deterministic terms selected by trend.
sigma2_
float
Concentrated innovation variance :math:`\hat{\sigma}^2`.
loglik_
float
Maximised exact Gaussian log-likelihood.
aic_
float
Akaike information criterion.
bic_
float
Bayesian information criterion.
resid_
np.ndarray
Standardised one-step-ahead prediction errors :math:`v_t/\sqrt{F_t}`.
state_
np.ndarray of shape (r,)
Predicted state at the first out-of-sample time point.
state_cov_
np.ndarray of shape (r, r)
Covariance of state_, in units of :math:`\sigma^2`.
n_obs_
int
Number of observations supplied to fit.

Notes

Complexity:

  • Training: O(\text{maxiter} \cdot k \cdot n r^3) where
r = \max(p^*, q^* + 1) is the state dimension of the expanded polynomials, n the sample size and k the number of free parameters (numerical gradients cost one filter pass each).
  • Prediction: O(h r^2) for h steps.
When to use SARIMAX:
  • You have exogenous regressors (price, temperature, a promotion
dummy) that should drive the series alongside its own dynamics.
  • You need a genuine seasonal (P, D, Q, m) component rather
than a plain ARIMA on a de-seasonalised series.
  • You want exact maximum likelihood and calibrated forecast
intervals from the Kalman variance rather than point forecasts alone.

Relationship to ARIMA: the two are deliberately not interchangeable. ARIMA is the light, fast option: it estimates a non-seasonal (p, d, q) model by Yule-Walker plus conditional-sum-of-squares refinement, it ignores exogenous input (its fit signature spells the argument _X), and its seasonal_order argument is accepted but not acted upon. Reach for SARIMAX whenever you actually need exogenous regressors, a working seasonal specification, exact-likelihood estimates, or forecast intervals; reach for ARIMA when you want a cheap non-seasonal point forecast and nothing more.

Interval caveat: predict_interval builds the forecast variance from the MA(\infty) weights of the expanded model, which includes the differencing operators, so intervals widen correctly with the horizon under differencing. Parameter-estimation uncertainty is not included -- the intervals condition on the fitted parameters.

References

Box2015
Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control, 5th ed. Wiley. :doi:`10.1111/jtsa.12194`
Harvey1990
Harvey, A. C. (1990). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press. :doi:`10.1017/CBO9781107049994`
Monahan1984
Monahan, J. F. (1984). A note on enforcing stationarity in autoregressive-moving average models. Biometrika, 71(2), 403-404. :doi:`10.1093/biomet/71.2.403`
Jones1980
Jones, R. H. (1980). Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics, 22(3), 389-395. :doi:`10.2307/1268324`

See Also

ARIMA

Lighter non-seasonal model with no exogenous support; see the Notes above for how to choose.
TBATS

Trigonometric seasonal state-space model for high-frequency or non-integer periods.
ExponentialSmoothing

Holt-Winters smoothing with integer seasonality.
python 
>>> import numpy as np
>>> from tuiml.algorithms.timeseries.sarimax import SARIMAX
>>> rng = np.random.default_rng(0)
>>> eps = rng.normal(scale=0.5, size=400)
>>> y = np.zeros(400)
>>> for t in range(1, 400):
...     y[t] = 0.7 * y[t - 1] + eps[t]
>>> model = SARIMAX(order=(1, 0, 0)).fit(y)
>>> bool(abs(model.ar_params_[0] - 0.7) < 0.05)
True
>>> model.predict(steps=3).shape
(3,)

An exogenous regression recovers the true coefficient:

python 
>>> x = rng.normal(size=(300, 1))
>>> y2 = 3.0 * x[:, 0] + 0.01 * rng.normal(size=300)
>>> m2 = SARIMAX(order=(0, 0, 0)).fit(y2, x)
>>> bool(abs(m2.exog_params_[0] - 3.0) < 0.01)
True
>>> bool(np.all(np.abs(m2.predict(steps=2, X=np.zeros((2, 1)))) < 0.01))
True

Calling predict without the future exogenous values is an error:

python 
>>> m2.predict(steps=2)
Traceback (most recent call last):
    ...
ValueError: This SARIMAX was fitted with 1 exogenous regressor(s); predict() requires X with the future values for those regressors.

Methods

get_parameter_schema (cls) -> Dict[str, Dict[str, Any]]

Return JSON Schema for constructor parameters.

get_capabilities (cls) -> List[str]

Return supported capabilities.

get_complexity (cls) -> str

Return complexity analysis.

get_references (cls) -> List[str]

Return academic citations.

fit (self, y: np.ndarray, X: Optional[np.ndarray]=None) -> 'SARIMAX'

Fit the model by exact maximum likelihood.

Parameters
y
np.ndarray of shape (n_samples,)
The time series to model.
X
np.ndarray of shape (n_samples, k_exog) = None
Exogenous regressors aligned row-wise with y.
Returns
self
SARIMAX
Fitted estimator.
predict (self, steps: int=1, X: Optional[np.ndarray]=None) -> np.ndarray

Forecast future values.

Parameters
steps
int = 1
Number of future time steps to forecast.
X
np.ndarray of shape (steps, k_exog) = None
Future values of the exogenous regressors. Required if and only if the model was fitted with exogenous regressors.
Returns
forecast
np.ndarray of shape (steps,)
Point forecasts on the scale of the original series.
forecast_variance (self, steps: int=1) -> np.ndarray

Return the forecast error variance for each horizon.

Parameters
steps
int = 1
Forecast horizon.
Returns
var
np.ndarray of shape (steps,)
Forecast error variances.
predict_interval (self, steps: int=1, alpha: float=0.05, X: Optional[np.ndarray]=None)

Forecast with a Gaussian prediction interval.

Parameters
steps
int = 1
Number of future time steps to forecast.
alpha
float = 0.05
Significance level; 0.05 gives a 95% interval.
X
np.ndarray of shape (steps, k_exog) = None
Future exogenous regressors, required if the model uses them.
Returns
forecast
np.ndarray of shape (steps,)
Point forecasts.
lower
np.ndarray of shape (steps,)
Lower interval bounds.
upper
np.ndarray of shape (steps,)
Upper interval bounds.
fit_predict (self, y: np.ndarray, steps: int=1) -> np.ndarray

Fit the model and forecast in one call.

Parameters
y
np.ndarray of shape (n_samples,)
Time series to fit.
steps
int = 1
Number of future time steps to forecast.
Returns
forecast
np.ndarray of shape (steps,)
Forecasted values.