---
title: Multilevel Seasonal SARMA Models
url: https://www.emergentmind.com/topics/multilevel-seasonal-sarma
type: topic
---

# Multilevel Seasonal SARMA Models

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” [2008.12340] and rule-based adaptations for special days in forecasting [1803.09590].

## 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:
  $$
  \phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p, \quad
  \theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q,
  $$
- Seasonal AR and MA polynomials:
  $$
  \Phi_{s_i}(B^{s_i}) = 1 - \Phi_{i,1} B^{s_i} - \cdots - \Phi_{i,P_i} B^{P_i s_i},
  $$
  $$
  \Theta_{s_i}(B^{s_i}) = 1 + \Theta_{i,1} B^{s_i} + \cdots + \Theta_{i,Q_i} B^{Q_i s_i}.
  $$

The multilevel seasonal SARMA (“MS”) model is:
$$
\phi(B) \cdot \prod_{i=1}^k \Phi_{s_i}(B^{s_i}) X_t =
\theta(B) \cdot \prod_{i=1}^k \Theta_{s_i}(B^{s_i}) \epsilon_t
$$
or, equivalently (explicit form),
$$
X_t - \sum_{m=1}^p \phi_m X_{t-m} - \sum_{i=1}^k \sum_{m=1}^{P_i} \Phi_{i,m} X_{t - m s_i}
= \epsilon_t + \sum_{n=1}^q \theta_n \epsilon_{t-n} + \sum_{i=1}^k \sum_{n=1}^{Q_i} \Theta_{i,n} \epsilon_{t-n s_i}
$$

This structure allows for interaction between multiple seasonalities and short-run dynamics, yielding a unified, interpretable, and flexible framework [2008.12340][1803.09590].

## 2. Parameter Estimation and Order Selection

All AR and MA coefficients $\{\phi_m, \theta_n, \Phi_{i,m}, \Theta_{i,n}\}$ 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:
$$
l(\Omega) = \sum_{t=1}^N \epsilon_t^2(\Omega),
$$
with $\epsilon_t$ computed from the model recursively. Optimization typically uses quasi-Newton techniques (e.g., BFGS) with analytical derivatives, ensuring joint estimation of all coefficients [2008.12340].

Lag order selection requires choosing values for $p, q, \{P_i, Q_i\}_{i=1}^k$. Common information criteria include:
- AIC: $-2\hat\ell + 2K$
- BIC: $-2\hat\ell + K \log N$
- Bridge Criterion (BC): $-2\hat\ell + 2\gamma(N) K$, with $\gamma(N) \to 1$ for small $N$ and $\gamma(N) \to \frac{1}{2}\log N$ for large $N$.

The total parameter count is $K = p + q + \sum_{i=1}^k (P_i + Q_i) + 1$ (variance). The model minimizing the selected information criterion is chosen [2008.12340].

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 [1803.09590].

## 3. Forecasting and Algorithmic Procedures

Given estimated parameters, $h$-step-ahead forecasts are generated by forward iteration of the SARMA equations, setting future innovations $\hat\epsilon_{t} = 0$ for $t > T$. The forecast recursion for $h \geq 1$:
$$
\hat X_{T+h} =
\sum_{m=1}^p \phi_m \hat X_{T+h-m}
+ \sum_{i=1}^k \sum_{m=1}^{P_i} \Phi_{i,m} \hat X_{T+h-m s_i}
+ \sum_{n=1}^q \theta_n \hat\epsilon_{T+h-n}
+ \sum_{i=1}^k \sum_{n=1}^{Q_i} \Theta_{i,n} \hat\epsilon_{T+h-n s_i}
$$
with the convention $\hat\epsilon_{t} = 0$ for $t > T$ [2008.12340].

For rule-based adaptations, such as forecasting French electricity load, special rules adapt the order and the annual seasonal lag $s_3(t)$ depending on whether $t$ is a normal or special day, and can propagate through the recursive forecasts [1803.09590]. 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 $IN_t$ flag normal vs. special day.
- The intrayear seasonal lag $s_3(t)$ 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 $s_3(t)$ switch depending on normal or special days.
- All parameters, including error variances for each regime, are estimated by maximizing the Gaussian log-likelihood
$$
\ell = -\frac{1}{2}\left[n_N \log(2\pi \sigma_N^2) + n_S \log(2\pi \sigma_S^2) + \sum_t IN_t \frac{\epsilon_t^2}{\sigma_N^2} + \sum_t (1-IN_t) \frac{\epsilon_t^2}{\sigma_S^2} \right]
$$
with $n_N, n_S$ the counts of normal/special-day observations. Optimization uses, e.g., Nelder–Mead simplex [1803.09590].

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 $z$ are outside the unit circle $|z| = 1$:
$$
\phi(z) \prod_{i=1}^k \Phi_{s_i}(z^{s_i}) = 0 \implies |z| > 1
$$
Similarly, invertibility requires all roots of the composite MA polynomial to satisfy $|z| > 1$:
$$
\theta(z) \prod_{i=1}^k \Theta_{s_i}(z^{s_i}) = 0 \implies |z| > 1
$$
When these conditions are satisfied for the estimated coefficients, the usual asymptotic theory and predictive properties of ARMA models apply [2008.12340].

## 6. Empirical Evidence and Comparative Performance

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

**Simulation results for 1-step forecasts** [2008.12340]:
- **Single-seasonal sine-wave plus ARMA noise (period 50):** $\mathrm{RMSE_{MS}} = 4.61$ vs. Prophet = 15.45, TBATS = 4.64
- **Double-seasonal sine-wave (periods 50 & 15):** $\mathrm{RMSE_{MS}} = 4.29$ vs. Prophet = 52.38, TBATS = 3.91
- **Nontrigonometric double seasonality (lags 50 & 15):** $\mathrm{RMSE_{MS}} = 31.80$ vs. Prophet = 68.04, TBATS = 45.32

**Real-world time series:**
- **PJM electricity load (hourly, daily + weekly):** Standardized CMSE at $n=24$: MS = 0.30 vs. Prophet = 0.55, TBATS = 0.38
- **Total sunspot number (monthly $\approx$11-year cycle):** 1-step standardized CMSE: MS = 8.26$\times 10^2$, Prophet = 47.7$\times 10^2$, TBATS = 9.96$\times 10^2$
- **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** [1803.09590]:
- **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., $s_1=48$, $s_2=336$, $s_3=17,520$ for half-hourly data).
- Orders: empirical search over $p, q, P_i, Q_i \leq 3$ 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 [2008.12340][1803.09590].

Source: https://www.emergentmind.com/topics/multilevel-seasonal-sarma