STL decomposition (seasonal-trend with loess)

Deep dive

How it works

In: series y, period np, spans, pass counts · Out: trend, seasonal, remainder, robustness weights · State: none (batch) · Cost: O(n) per pass with jumps and interpolation

Parameters: np period, ns seasonal span, nt trend span, nl low-pass span, ni inner passes, no outer passes. Start with trend T = 0 and robustness weights 1.

Inner loop (ni times):

1. Detrend: y - T.
2. For each phase, loess (span ns) over that phase's values, evaluated at every point plus one extra at each end: C, length n + 2*np.
3. Low-pass of C: moving averages of length np, np, 3, then loess (span nl): L, length n.
4. Seasonal S = C[np .. n+np] - L (L removes low-frequency content that belongs to the trend).
5. Trend T = loess (span nt) of y - S.

Outer loop: r = y - T - S, h = 6 * median(|r|), rho = (1 - (|r|/h)^2)^2 for |r| < h, else 0; rerun the inner loop with rho multiplying the loess weights in steps 2 and 5.

R defaults: nt = nextodd(ceil(1.5*np / (1 - 1.5/ns))), nl = smallest odd >= np (statsmodels: > np); seasonal loess degree 0 in R, 1 in statsmodels and the paper.

Code

Plain loess at every point, degree 1, tricube weights; loess_at is the only primitive.

/// Local linear fit at position `x` using the `q` nearest of the points 0..n-1 (tricube weights
/// times optional robustness weights `rw`). Returns None if all weights are zero.
fn loess_at(y: &[f64], rw: Option<&[f64]>, q: usize, x: f64) -> Option<f64> {
    let n = y.len();
    let q = q.min(n); // R's stl also widens h when q > n; clamping is enough for ns, nt <= n
    // window of the q consecutive points nearest to x
    let lo = ((x - (q as f64 - 1.0) / 2.0).round().max(0.0) as usize).min(n - q);
    let h = (x - lo as f64).max((lo + q - 1) as f64 - x).max(1e-9);
    let (mut sw, mut sx, mut sxx, mut sy, mut sxy) = (0.0, 0.0, 0.0, 0.0, 0.0);
    for i in lo..lo + q {
        let d = (i as f64 - x).abs() / h;
        let mut w = if d < 1.0 { (1.0 - d * d * d).powi(3) } else { 0.0 };
        if let Some(r) = rw { w *= r[i]; }
        let xi = i as f64 - x; // centre at x so the fit value is the intercept
        sw += w; sx += w * xi; sxx += w * xi * xi; sy += w * y[i]; sxy += w * xi * y[i];
    }
    if sw <= 1e-12 { return None; }
    let det = sw * sxx - sx * sx;
    // if the window is too narrow for a slope (det ~ 0), fall back to the weighted mean
    Some(if det.abs() > 1e-9 * sw * sw * h * h { (sxx * sy - sx * sxy) / det } else { sy / sw })
}

fn loess(y: &[f64], rw: Option<&[f64]>, q: usize, from: isize, to: isize) -> Vec<f64> {
    (from..=to).map(|x| loess_at(y, rw, q, x as f64).unwrap_or(y[(x.max(0) as usize).min(y.len() - 1)])).collect()
}

fn moving_avg(x: &[f64], k: usize) -> Vec<f64> {
    x.windows(k).map(|w| w.iter().sum::<f64>() / k as f64).collect()
}

pub struct Stl { pub trend: Vec<f64>, pub seasonal: Vec<f64>, pub remainder: Vec<f64> }

/// STL (Cleveland et al. 1990). np = period, ns = seasonal span (odd, >= 7), nt = trend span,
/// nl = low-pass span, ni inner passes, no outer robustness passes.
pub fn stl(y: &[f64], np: usize, ns: usize, nt: usize, nl: usize, ni: usize, no: usize) -> Stl {
    let n = y.len();
    let (mut trend, mut seasonal) = (vec![0.0; n], vec![0.0; n]);
    let mut rw: Option<Vec<f64>> = None;
    for outer in 0..=no {
        for _ in 0..ni {
            let detr: Vec<f64> = y.iter().zip(&trend).map(|(a, b)| a - b).collect();
            // 1-2: smooth each cycle-subseries, extended by one point at both ends -> n + 2*np values
            let mut c = vec![0.0; n + 2 * np];
            for s in 0..np {
                let idx: Vec<usize> = (s..n).step_by(np).collect();
                let sub: Vec<f64> = idx.iter().map(|&i| detr[i]).collect();
                let sw: Option<Vec<f64>> = rw.as_ref().map(|r| idx.iter().map(|&i| r[i]).collect());
                let sm = loess(&sub, sw.as_deref(), ns, -1, sub.len() as isize);
                for (j, v) in sm.iter().enumerate() { // j-th value sits at time s + (j-1)*np, shifted by np
                    let t = s + j * np;
                    if t < c.len() { c[t] = *v; }
                }
            }
            // 3: low-pass filter of c: MA(np), MA(np), MA(3), loess(nl) -> length n
            let lp = moving_avg(&moving_avg(&moving_avg(&c, np), np), 3);
            let low = loess(&lp, None, nl, 0, n as isize - 1);
            // 4: seasonal = smoothed subseries minus its low-frequency part
            for i in 0..n { seasonal[i] = c[i + np] - low[i]; }
            // 5-6: trend = loess of the deseasonalized series
            let ds: Vec<f64> = y.iter().zip(&seasonal).map(|(a, b)| a - b).collect();
            trend = loess(&ds, rw.as_deref(), nt, 0, n as isize - 1);
        }
        if outer < no { // bisquare robustness weights from the remainder
            let r: Vec<f64> = (0..n).map(|i| (y[i] - trend[i] - seasonal[i]).abs()).collect();
            let mut s = r.clone(); s.sort_by(|a, b| a.partial_cmp(b).unwrap());
            let h = 6.0 * 0.5 * (s[(n - 1) / 2] + s[n / 2]); // 6 * median |r|
            rw = Some(r.iter().map(|&v| {
                let u = v / h;
                if u < 0.001 { 1.0 } else if u < 0.999 { (1.0 - u * u).powi(2) } else { 0.0 }
            }).collect());
        }
    }
    let remainder = (0..n).map(|i| y[i] - trend[i] - seasonal[i]).collect();
    Stl { trend, seasonal, remainder }
}

/// Defaults of R's stl() for a seasonal window of ns (odd >= 7); robust => 1 inner, 15 outer.
pub fn stl_default(y: &[f64], np: usize, ns: usize, robust: bool) -> Stl {
    let nl = if np % 2 == 0 { np + 1 } else { np + 2 };           // smallest odd > np (statsmodels; R: >= np)
    let mut nt = (1.5 * np as f64 / (1.0 - 1.5 / ns as f64)).ceil() as usize;
    if nt % 2 == 0 { nt += 1; }
    if robust { stl(y, np, ns, nt, nl, 1, 15) } else { stl(y, np, ns, nt, nl, 2, 0) }
}

fn main() {
    let np = 12;
    let mut x: u64 = 7;
    let mut noise = || { x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
                         ((x >> 33) as f64 / (1u64 << 31) as f64 - 0.5) * 2.0 };
    let n = 120;
    let tr = |t: usize| 50.0 + 0.5 * t as f64;
    let se = |t: usize| 8.0 * (2.0 * std::f64::consts::PI * t as f64 / np as f64).sin();
    let mut y: Vec<f64> = (0..n).map(|t| tr(t) + se(t) + noise()).collect();
    for &t in &[30usize, 31, 77] { y[t] += 60.0; } // additive outliers
    let plain = stl_default(&y, np, 7, false);
    let robust = stl_default(&y, np, 7, true);
    for r in [&plain, &robust] {
        for i in 0..n { assert!((r.trend[i] + r.seasonal[i] + r.remainder[i] - y[i]).abs() < 1e-9); }
    }
    let mae = |s: &[f64]| (0..n).map(|t| (s[t] - se(t)).abs()).sum::<f64>() / n as f64;
    println!("seasonal MAE plain {:.3} robust {:.3}", mae(&plain.seasonal), mae(&robust.seasonal));
    println!("remainder at outliers (robust): {:.1} {:.1} {:.1}", robust.remainder[30], robust.remainder[31], robust.remainder[77]);
    println!("remainder at outliers (plain):  {:.1} {:.1} {:.1}", plain.remainder[30], plain.remainder[31], plain.remainder[77]);
    assert!(mae(&robust.seasonal) < 0.5 && mae(&robust.seasonal) < mae(&plain.seasonal));
    assert!(robust.remainder[30] > 45.0 && robust.remainder[77] > 45.0);
    let tmae = (6..n - 6).map(|t| (robust.trend[t] - tr(t)).abs()).sum::<f64>() / (n - 12) as f64;
    println!("trend MAE (robust, interior) {:.3}", tmae);
    assert!(tmae < 1.0);
}

The test checks trend + seasonal + remainder == y and that, with three +60 outliers, the robust fit recovers the seasonal shape (MAE 0.30 vs 2.41 plain) and leaves the spikes in the remainder (about 60 vs about 39). Compiled and run by me. Unlike R/statsmodels: no jump interpolation, window clamped when the span exceeds a subseries, no periodic mode.

Complexity and constants

• One loess evaluation visits q points (about 10 flops each); a pass costs roughly n * (nt + nl) plus n * ns / np per phase-smoothing.
• Passes: ni * (no + 1) = 2 non-robust, 16 robust (R defaults), so robust costs about 8x.
• Measured by me (Rust, Apple silicon, release, scalar, no jumps): n = 100,000, np = 24, ns = 7: about 27 ms non-robust, about 230 ms robust. One machine, order of magnitude.
• Memory: a few length-n f64 arrays. The trend span nt dominates cost: hundreds of points per window.

Efficiency

• Bound: compute over small windows; gather each phase into a contiguous array once per pass to avoid stride-np access.
• Jumps (R s.jump, t.jump, l.jump, default ceiling(window/10)): evaluate loess every jump-th point and interpolate linearly, cutting loess work by about the jump factor. statsmodels defaults to jumps of 1 (exact).
• Robustness is the cost multiplier: use it only when outliers exist.
• Interior tricube weights depend only on offsets, so precompute them per span (recomputed here for clarity).
• Many series: parallelize across series; for latest-remainder needs, decompose a recent window of a few periods plus nt.
• Cheaper alternative: classical decomposition (one centered moving average plus per-phase means, O(n), no robustness), or exponential smoothing for online use.

Implementation notes

• ns = 7 is a common start; fixed seasonality means a very large ns (R "periodic").
• nt must exceed the period or the trend absorbs the seasonal pattern; the default formula ensures it.
• Edge points use one-sided windows and are the noisiest; they change when new data arrives.
• Missing values: interpolate first, then give those points robustness weight 0.
• For anomalies, standardize the remainder with a robust scale (MAD or h/6): anomaly detection, robust statistics.

Pitfalls

• Wrong period (7 for daily data with weekly cycle, 24 for hourly with daily): the seasonal part looks like noise.
• Non-integer periods (365.25, 52.18) are unsupported; resample or use Fourier regressors.
• Multiplicative seasonality without a log transform leaves seasonal structure in the remainder.
• Too small ns or nt makes components follow noise and eats the remainder you wanted for anomalies.
• Level shifts are smoothed over, not located: change point detection.
• Decomposing including the test period leaks information; fit inside each fold (time-series cross-validation).
• R and statsmodels defaults differ (seasonal degree, jumps), so outputs differ slightly.

Variants and alternatives

In production

Seasonal adjustment of economic series such as US retail employment (FPP3 example), remainder-based alerting in monitoring, and "decompose, forecast the adjusted series, add back the seasonal naive" forecasting (R stlf), plus MSTL as a multi-seasonal baseline in forecasting toolkits.

Implementations

Sources

• Robert B. Cleveland, William S. Cleveland, Jean E. McRae, Irma Terpenning, STL: A Seasonal-Trend Decomposition Procedure Based on Loess, Journal of Official Statistics 6(1), 3-73, 1990. paper
• Rob J. Hyndman, George Athanasopoulos, Forecasting: Principles and Practice, 3rd ed., OTexts 2021, ch. 3.6 (STL decomposition). book
• Kasun Bandara, Rob J. Hyndman, Christoph Bergmeir, MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns, arXiv 2021. paper
• Qingsong Wen, Jingkun Gao, Xiaomin Song, Liang Sun, Huan Xu, Shenghuo Zhu, RobustSTL: A Robust Seasonal-Trend Decomposition Algorithm for Long Time Series, AAAI 2019. paper