Seasonal Decomposition in Python — Deep Dive

How STL works under the hood

STL is an iterative algorithm that alternates between two loops:

Inner loop (repeated inner_iter times, default 2):

1. Detrend the series by subtracting the current trend estimate.
2. Apply seasonal smoothing (Loess on each subseries of same-period positions) to get the seasonal component.
3. Apply a low-pass filter to the seasonal component and subtract it to prevent leakage.
4. Deseasonalize the series by subtracting the seasonal estimate.
5. Apply trend smoothing (Loess on the deseasonalized series) to update the trend.

Outer loop (repeated outer_iter times, default 0 for non-robust, 15 for robust):

1. Compute residuals from the inner loop result.
2. Assign robustness weights that down-weight large residuals.
3. Re-run the inner loop using these weights.

The Loess (Locally Estimated Scatterplot Smoothing) regression at each step fits weighted polynomial regressions to local neighborhoods of the data, producing smooth but flexible estimates.

from statsmodels.tsa.seasonal import STL

stl = STL(
    series,
    period=7,
    seasonal=13,
    trend=None,      # auto-calculated
    low_pass=None,   # auto-calculated
    seasonal_deg=1,  # degree of Loess polynomial (0 or 1)
    trend_deg=1,
    low_pass_deg=1,
    robust=True,
    seasonal_jump=1,  # speed up: only compute every N-th point
    trend_jump=1,
    low_pass_jump=1,
)
result = stl.fit(inner_iter=5, outer_iter=15)

Parameter relationships

The seasonal window must be odd and at least 7. The trend window must be odd and at least:

trend ≥ ceil(1.5 × period / (1 − 1.5 / seasonal))

This ensures the trend smoother is wide enough to not absorb seasonal variation. When trend=None, statsmodels calculates this automatically.

MSTL: multiple seasonal decomposition

Real-world data often has overlapping seasonal patterns. Hourly electricity demand has daily (24h), weekly (168h), and yearly (8766h) cycles. MSTL (Multiple STL) handles this:

from statsmodels.tsa.seasonal import MSTL

mstl = MSTL(
    series,
    periods=[24, 168],           # daily and weekly
    windows=[25, 169],           # seasonal windows (odd, ≥ 7)
    stl_kwargs={"robust": True},
)
result = mstl.fit()

# result.trend — single trend
# result.seasonal — DataFrame with one column per seasonal period
# result.resid — residuals

MSTL works iteratively: it extracts the shortest-period seasonality first, then the next, and so on. The order matters — shorter periods are extracted before longer ones because they are easier to identify.

Custom MSTL with iterative STL

For more control, implement iterative decomposition manually:

import pandas as pd
from statsmodels.tsa.seasonal import STL

def multi_seasonal_decompose(series, periods, n_iterations=3, **stl_kwargs):
    """Iterative multi-seasonal STL decomposition."""
    seasonal_components = {p: pd.Series(0, index=series.index) for p in periods}
    remainder = series.copy()
    
    for iteration in range(n_iterations):
        for period in sorted(periods):
            # Remove other seasonal components
            adjusted = remainder + seasonal_components[period]
            
            stl = STL(adjusted, period=period, **stl_kwargs)
            result = stl.fit()
            
            seasonal_components[period] = result.seasonal
            remainder = result.resid + result.trend
    
    # Final trend extraction
    stl_final = STL(remainder, period=min(periods), **stl_kwargs)
    final_result = stl_final.fit()
    
    return {
        "trend": final_result.trend,
        "seasonal": seasonal_components,
        "resid": final_result.resid,
    }

Diagnostic workflows

Strength of trend and seasonality

Quantify how much of the variance each component explains:

import numpy as np

def decomposition_strength(result):
    """Compute strength of trend and seasonal components (0 to 1)."""
    resid_var = np.var(result.resid.dropna())
    deseason_var = np.var((result.trend + result.resid).dropna())
    detrend_var = np.var((result.seasonal + result.resid).dropna())
    
    trend_strength = max(0, 1 - resid_var / deseason_var)
    seasonal_strength = max(0, 1 - resid_var / detrend_var)
    
    return {
        "trend_strength": trend_strength,
        "seasonal_strength": seasonal_strength,
    }

# Values close to 1 = strong component
# Values close to 0 = weak component

This metric (from Hyndman & Athanasopoulos) is used by the tsfeatures package to classify time series characteristics at scale.

Residual diagnostics

from scipy import stats
from statsmodels.stats.diagnostic import acorr_ljungbox

def diagnose_residuals(residuals):
    """Full residual diagnostic suite."""
    resid = residuals.dropna()
    
    # Normality
    shapiro_stat, shapiro_p = stats.shapiro(resid[:5000])  # limit for large series
    
    # Autocorrelation
    lb = acorr_ljungbox(resid, lags=[7, 14, 21], return_df=True)
    
    # Heteroscedasticity (is variance changing?)
    half = len(resid) // 2
    first_var = np.var(resid[:half])
    second_var = np.var(resid[half:])
    var_ratio = max(first_var, second_var) / min(first_var, second_var)
    
    return {
        "normality_p": shapiro_p,
        "is_normal": shapiro_p > 0.05,
        "autocorrelation": lb.to_dict(),
        "variance_ratio": var_ratio,
        "is_homoscedastic": var_ratio < 2.0,
    }

If residuals show autocorrelation, the decomposition period may be wrong or the data has additional structure. If variance changes over time (heteroscedasticity), switch from additive to multiplicative decomposition.

Decomposition in production pipelines

Anomaly detection using decomposition

def decomposition_anomaly_detector(series, period, threshold_sigma=3.0):
    """Detect anomalies as points with extreme residuals."""
    stl = STL(series, period=period, robust=True)
    result = stl.fit()
    
    resid = result.resid
    median = resid.median()
    mad = np.median(np.abs(resid - median))  # robust std estimate
    robust_std = 1.4826 * mad
    
    z_scores = (resid - median) / robust_std
    anomalies = series[np.abs(z_scores) > threshold_sigma]
    
    return anomalies, z_scores

anomalies, scores = decomposition_anomaly_detector(series, period=7)
print(f"Found {len(anomalies)} anomalies")

Using MAD (Median Absolute Deviation) instead of standard deviation makes the detector robust to the very outliers it is trying to find.

Seasonal adjustment for reporting

def seasonally_adjust(series, period, multiplicative=False):
    """Remove seasonal component for cleaner trend analysis."""
    stl = STL(series, period=period, robust=True)
    result = stl.fit()
    
    if multiplicative:
        # For multiplicative: divide by seasonal factor
        seasonal_factor = result.seasonal / result.trend
        adjusted = series / (1 + seasonal_factor)
    else:
        adjusted = series - result.seasonal
    
    return adjusted, result

Government agencies like the US Census Bureau and Bureau of Labor Statistics use X-13ARIMA-SEATS for official seasonal adjustment, but STL-based adjustment is a practical alternative for business applications.

Streaming decomposition

For real-time data, full STL is impractical because it requires the entire series. Use an online approximation:

class OnlineSeasonalDecomposer:
    """Lightweight online decomposition using exponential smoothing."""
    
    def __init__(self, period, alpha=0.3):
        self.period = period
        self.alpha = alpha
        self.seasonal = [0.0] * period
        self.level = None
        self.step = 0
    
    def update(self, value):
        idx = self.step % self.period
        
        if self.level is None:
            self.level = value
            self.seasonal[idx] = 0.0
        else:
            deseasonalized = value - self.seasonal[idx]
            self.level = self.alpha * deseasonalized + (1 - self.alpha) * self.level
            self.seasonal[idx] = (
                self.alpha * (value - self.level)
                + (1 - self.alpha) * self.seasonal[idx]
            )
        
        residual = value - self.level - self.seasonal[idx]
        self.step += 1
        
        return {
            "trend": self.level,
            "seasonal": self.seasonal[idx],
            "residual": residual,
        }

This sacrifices accuracy for latency — suitable for monitoring dashboards and alerting where rough decomposition is better than none.

Period detection

When the seasonal period is unknown, use spectral analysis:

from scipy.signal import periodogram
import numpy as np

def detect_periods(series, max_periods=3):
    """Detect dominant seasonal periods using spectral analysis."""
    freqs, power = periodogram(series.dropna().values, detrend="linear")
    
    # Skip frequency 0 (DC component)
    freqs = freqs[1:]
    power = power[1:]
    
    # Find peaks
    from scipy.signal import find_peaks
    peaks, properties = find_peaks(power, height=np.median(power) * 5)
    
    if len(peaks) == 0:
        return []
    
    # Sort by power
    peak_order = np.argsort(power[peaks])[::-1][:max_periods]
    dominant_freqs = freqs[peaks[peak_order]]
    
    periods = [round(1.0 / f) for f in dominant_freqs if f > 0]
    return periods

Always validate detected periods against domain knowledge. A detected period of 11.7 for monthly data probably means yearly seasonality (12) with slight estimation noise.

The one thing to remember: Production-grade decomposition goes beyond calling seasonal_decompose — it requires choosing between STL and MSTL based on the number of seasonal patterns, tuning smoothing parameters for the right balance of flexibility and stability, and building diagnostic checks that catch when the decomposition fails to capture the data’s true structure.