Multilevel Seasonal SARMA Models

• Multilevel Seasonal SARMA models are extensions of SARMA that incorporate multiple, interacting seasonal cycles in univariate time series.
• They use maximum likelihood estimation and criteria like AIC/BIC to jointly identify nonseasonal and seasonal AR/MA orders in complex settings.
• Empirical evidence shows these models enhance forecast accuracy over standard methods, especially when adapting to special day effects.

Multilevel Seasonal SARMA (Seasonal Autoregressive Moving Average) models extend conventional SARMA techniques to accommodate multiple, potentially interacting, seasonal cycles within univariate time series. This structure is especially suited for data exhibiting complex seasonalities, such as those arising in hourly electricity load, meteorological records, retail demand, and other domains featuring several periodic components with short- and long-term patterns. Multilevel seasonal SARMA models have received recent methodological elaboration, notably as in “Forecasting with Multiple Seasonality” (Xie et al., 2020) and rule-based adaptations for special days in forecasting (Arora et al., 2018).

1. Model Formulation

A multilevel seasonal SARMA model is defined for a univariate series $\{X_t\}$ with innovations $\epsilon_t \sim$ i.i.d. $N(0, \sigma^2)$, and is constructed as follows:

Let $B$ denote the backshift operator: $B X_t = X_{t-1}$. Suppose $k$ seasonal cycles with periods $s_1, \dots, s_k$, nonseasonal AR and MA orders $(p, q)$, and seasonal orders $(P_i, Q_i)$ for each $i = 1, \ldots, k$.

Define:

• Nonseasonal AR and MA polynomials:

$\epsilon_t \sim$0

• Seasonal AR and MA polynomials:

$\epsilon_t \sim$1

$\epsilon_t \sim$2

The multilevel seasonal SARMA (“MS”) model is:

$\epsilon_t \sim$3

or, equivalently (explicit form),

$\epsilon_t \sim$4

This structure allows for interaction between multiple seasonalities and short-run dynamics, yielding a unified, interpretable, and flexible framework (Xie et al., 2020, Arora et al., 2018).

2. Parameter Estimation and Order Selection

All AR and MA coefficients $\epsilon_t \sim$5 and innovation variance are jointly estimated by maximum likelihood under Gaussianity. This is equivalent to minimizing the sum of squared one-step-ahead forecast errors:

$\epsilon_t \sim$6

with $\epsilon_t \sim$7 computed from the model recursively. Optimization typically uses quasi-Newton techniques (e.g., BFGS) with analytical derivatives, ensuring joint estimation of all coefficients (Xie et al., 2020).

Lag order selection requires choosing values for $\epsilon_t \sim$8. Common information criteria include:

• AIC: $\epsilon_t \sim$9
• BIC: $N(0, \sigma^2)$0
• Bridge Criterion (BC): $N(0, \sigma^2)$1, with $N(0, \sigma^2)$2 for small $N(0, \sigma^2)$3 and $N(0, \sigma^2)$4 for large $N(0, \sigma^2)$5.

The total parameter count is $N(0, \sigma^2)$6 (variance). The model minimizing the selected information criterion is chosen (Xie et al., 2020).

In the context of special day forecasting, rule-based triple seasonal SARMA models may further select model orders via manual grid search or AIC, and encode additional rules for seasonal period selection or parameter switching (Arora et al., 2018).

3. Forecasting and Algorithmic Procedures

Given estimated parameters, $N(0, \sigma^2)$7-step-ahead forecasts are generated by forward iteration of the SARMA equations, setting future innovations $N(0, \sigma^2)$8 for $N(0, \sigma^2)$9. The forecast recursion for $B$0:

$B$1

with the convention $B$2 for $B$3 (Xie et al., 2020).

For rule-based adaptations, such as forecasting French electricity load, special rules adapt the order and the annual seasonal lag $B$4 depending on whether $B$5 is a normal or special day, and can propagate through the recursive forecasts (Arora et al., 2018). For probabilistic/density forecasts, Monte Carlo methods simulate future innovations from appropriate residual variance distributions conditioned on day type.

4. Special-Day Adaptation and Rule-Based Mechanisms

Rule-based multilevel SARMA implementations extend the model to accommodate regime changes such as holidays or special days:

• Binary indicators $B$6 flag normal vs. special day.
• The intrayear seasonal lag $B$7 varies by a rule linking each special day to the closest previous occurrence of the same type, partitioned into categories A–G (e.g., by weekday/weekend).
• AR and MA parameters at seasonal lag $B$8 switch depending on normal or special days.
• All parameters, including error variances for each regime, are estimated by maximizing the Gaussian log-likelihood

$B$9

with $B X_t = X_{t-1}$0 the counts of normal/special-day observations. Optimization uses, e.g., Nelder–Mead simplex (Arora et al., 2018).

This rule-based adaptation significantly improves forecast accuracy on special and proximity days by dynamically adjusting seasonal memory and model parameters.

5. Theoretical Properties

Stationarity and invertibility properties follow classical ARMA/SARMA theory. The multilevel seasonal SARMA model is stationary if all roots of the joint AR polynomial (including all seasonal AR) in $B X_t = X_{t-1}$1 are outside the unit circle $B X_t = X_{t-1}$2:

$B X_t = X_{t-1}$3

Similarly, invertibility requires all roots of the composite MA polynomial to satisfy $B X_t = X_{t-1}$4:

$B X_t = X_{t-1}$5

When these conditions are satisfied for the estimated coefficients, the usual asymptotic theory and predictive properties of ARMA models apply (Xie et al., 2020).

6. Empirical Evidence and Comparative Performance

Extensive empirical studies evaluate multilevel seasonal SARMA (MS) models against established approaches:

Simulation results for 1-step forecasts (Xie et al., 2020):

• Single-seasonal sine-wave plus ARMA noise (period 50): $B X_t = X_{t-1}$6 vs. Prophet = 15.45, TBATS = 4.64
• Double-seasonal sine-wave (periods 50 & 15): $B X_t = X_{t-1}$7 vs. Prophet = 52.38, TBATS = 3.91
• Nontrigonometric double seasonality (lags 50 & 15): $B X_t = X_{t-1}$8 vs. Prophet = 68.04, TBATS = 45.32

Real-world time series:

• PJM electricity load (hourly, daily + weekly): Standardized CMSE at $B X_t = X_{t-1}$9: MS = 0.30 vs. Prophet = 0.55, TBATS = 0.38
• Total sunspot number (monthly $k$011-year cycle): 1-step standardized CMSE: MS = 8.26$k$1, Prophet = 47.7$k$2, TBATS = 9.96$k$3
• Facebook Events (daily data, weekly + yearly): 1-step standardized CMSE: MS = 0.11 vs. Prophet = 0.45, TBATS = 0.13

Special day French electricity load (Arora et al., 2018):

• Rule-based triple seasonal SARMA (RB-SARMA): Outperformed benchmarks and neural networks at all horizons; MAPE from 0.53% (short horizon) to 3.22% (24h ahead).
• RB-SARMA models reduced CRPS density forecast error by half relative to un-adapted SARMA.

The MS framework consistently matches or outperforms additive models (Prophet), advanced state-space methods (TBATS), and neural architectures (LSTM, ANN), with particular strength in settings with unknown, multiple, and interacting seasonalities and in regimes requiring flexible regime switching.

7. Practical Implementation Considerations

Model specification for practical workloads:

• Common seasonal cycles: intraday, intraweek, intrayear (e.g., $k$4, $k$5, $k$6 for half-hourly data).
• Orders: empirical search over $k$7 suffices in many applications.
• Grid search or BC/AIC-based order selection.
• Efficient likelihood optimization using BFGS or simplex algorithms.
• Rule-based models require special-day calendars, indicator coding, and lag calculation for recurrences.

The flexibility of the multilevel seasonal SARMA model, particularly in its rule-based adaptations, enables unified handling of complex calendar effects and regime changes within a single univariate modeling framework (Xie et al., 2020, Arora et al., 2018).