季节-趋势分解（STL）方法详解

季节-趋势分解（STL）方法详解

在分析时间序列数据时，我们经常需要理解数据中隐藏的规律。比如零售商想知道销售额的增长是真实的业务增长还是仅仅是季节性因素，气候学家需要从温度数据中分离出长期变暖趋势和正常的季节变化，这些都需要一种强大的分解方法。STL（Seasonal and Trend decomposition using Loess）正是为此而生的统计方法，它能够将复杂的时间序列数据优雅地分解为三个易于理解的组成部分：趋势、季节性和余项。

数学原理与核心思想

STL的核心思想非常直观：任何时间序列都可以表示为三个加法组成部分的和。用数学公式表达就是：

Yν=Tν+Sν+RνY_\nu = T_\nu + S_\nu + R_\nuYν​=Tν​+Sν​+Rν​

其中YνY_\nuYν​代表在时间ν\nuν的观测值，TνT_\nuTν​是趋势分量，SνS_\nuSν​是季节分量，RνR_\nuRν​是余项分量。

STL方法的独特之处在于它使用了 LOESS （Locally Weighted Scatterplot Smoothing，局部加权散点图平滑）这一强大的非参数回归技术。LOESS的核心理念是"近朱者赤"——在估计某一点的值时，距离越近的点权重越大，距离越远的点权重越小。这种局部性使得STL能够灵活地适应数据的变化，而不像传统方法那样假设整个序列遵循同一个固定模式。

具体来说，LOESS使用一个优雅的权重函数，称为 三次立方权重函数 ：

W(u)={(1−∣u∣3)3当 ∣u∣≤10否则W(u) = \begin{cases} (1 - |u|^3)^3 & \text{当 } |u| \leq 1 \\ 0 & \text{否则} \end{cases}W(u)={(1−∣u∣3)30​当 ∣u∣≤1否则​

这个函数的妙处在于它在u=0u=0u=0时权重为1，随着距离增加权重平滑下降，到达边界时恰好为0，避免了突兀的截断。对于每个待估计的点xxx，我们首先找到qqq个最近的邻居点（qqq是平滑参数），然后计算每个邻居点的权重：

vi=W(∣xi−x∣λq(x))v_i = W\left(\frac{|x_i - x|}{\lambda_q(x)}\right)vi​=W(λq​(x)∣xi​−x∣​)

其中λq(x)\lambda_q(x)λq​(x)是到第qqq个最远邻居的距离。接下来，我们用这些权重进行局部多项式拟合，通常是线性或二次多项式，通过最小化加权误差平方和来估计该点的值。

算法的精妙设计：双循环结构

STL算法采用了一个巧妙的双循环结构，这种设计既保证了算法的稳健性，又提供了足够的灵活性。

内循环：交替估计的艺术

内循环是STL的核心，它通过反复迭代来交替估计季节分量和趋势分量。这个过程就像是两位艺术家轮流作画，每次都在对方的基础上进行改进，直到作品趋于完美。

第一步：去趋势 。我们从原始序列中减去当前的趋势估计Tν(k)T_ν^(k)Tν(​k)，得到去趋势序列。这一步的目的是让季节模式更加清晰，不受趋势干扰。

第二步：循环子序列平滑 。这是STL的独特创新。假设我们分析的是月度数据，那么我们把所有的一月份数据放在一起，所有的二月份数据放在一起，以此类推。对每个这样的子序列应用LOESS平滑，平滑参数为nsn_sns​（季节窗口）。这种方法允许季节模式随时间缓慢变化，比如夏季销售高峰可能逐年提前或推迟。

第三步：低通滤波 。对上一步得到的结果进行低通滤波，目的是去除季节分量中可能包含的趋势成分。这通过连续两次长度为3的移动平均，再加上一次LOESS平滑来实现。

第四步：计算季节分量。将循环子序列平滑的结果减去低通滤波的结果，得到纯净的季节分量：

Sν(k+1)=Cν(k+1)−Lν(k+1)S_\nu^{(k+1)} = C_\nu^{(k+1)} - L_\nu^{(k+1)}Sν(k+1)​=Cν(k+1)​−Lν(k+1)​

第五步：趋势平滑 。从原始序列中减去新的季节分量，对得到的序列进行LOESS平滑，平滑参数为ntn_tnt​（趋势窗口），得到新的趋势估计。

这个过程通常重复2-3次就能收敛，每次迭代都让趋势和季节的估计更加准确。

外循环：稳健性的保证

外循环是STL处理异常值的秘密武器。在现实数据中，异常值是不可避免的——可能是测量错误、特殊事件或者真实的极端情况。传统方法往往会被这些异常值严重影响，而STL通过外循环巧妙地解决了这个问题。

外循环使用双平方权重函数来计算稳健性权重：

ρν=B(∣Rν∣h)\rho_\nu = B\left(\frac{|R_\nu|}{h}\right)ρν​=B(h∣Rν​∣​)

其中B(u)=(1−u2)2B(u) = (1 - u^2)^2B(u)=(1−u2)2当0≤u≤10 \leq u \leq 10≤u≤1，否则为0，h=6×median(∣Rν∣)h = 6 \times \text{median}(|R_\nu|)h=6×median(∣Rν​∣)是基于余项中位数的尺度参数。

这个设计的精妙之处在于：正常的观测值（余项较小）获得接近1的权重，而异常值（余项较大）的权重会被大幅降低甚至降为0。这样，异常值主要影响余项分量，而不会扭曲趋势和季节的估计。

参数选择的艺术与科学

STL的灵活性很大程度上来自于其参数设置。理解这些参数如何影响分解结果，是成功应用STL的关键。

季节窗口nsn_sns​ 控制着季节分量的灵活性。较小的值（如7）允许季节模式快速变化，适合季节性不太稳定的数据；较大的值则假设季节模式相对稳定。选择时需要在灵活性和稳定性之间权衡。如果设置为"periodic"，则假设季节模式完全不变，这在某些情况下是合理的简化。

趋势窗口ntn_tnt​ 决定了趋势的平滑程度。STL提供了一个经验公式来计算默认值：最小的奇数大于1.5×np/(1−1.5/ns)1.5 \times n_p/(1-1.5/n_s)1.5×np​/(1−1.5/ns​)，其中npn_pnp​是季节周期。这个公式确保趋势变化的时间尺度大于季节变化，避免趋势捕获到季节性变化。

低通窗口nln_lnl​ 通常设置为最小的奇数不小于季节周期npn_pnp​。这确保低通滤波能够有效去除所有的季节性成分。

LOESS回归的深度数学推导

要理解STL的核心，我们必须深入探究LOESS（Locally Weighted Scatterplot Smoothing）的数学原理。LOESS的本质是在每个目标点附近构造一个局部加权最小二乘问题。

对于目标点x0x_0x0​，我们寻求最优的局部线性拟合：

f^(x0)=β^0+β^1(x0−x0)=β^0\hat{f}(x_0) = \hat{\beta}_0 + \hat{\beta}_1(x_0 - x_0) = \hat{\beta}_0f^​(x0​)=β^​0​+β^​1​(x0​−x0​)=β^​0​

其中系数通过最小化加权残差平方和获得：

SSE(β0,β1)=∑i=1nwi(x0)[yi−β0−β1(xi−x0)]2\text{SSE}(\beta_0, \beta_1) = \sum_{i=1}^n w_i(x_0)[y_i - \beta_0 - \beta_1(x_i - x_0)]^2SSE(β0​,β1​)=i=1∑n​wi​(x0​)[yi​−β0​−β1​(xi​−x0​)]2

权重函数wi(x0)w_i(x_0)wi​(x0​)基于距离计算。首先定义标准化距离：

ui=∣xi−x0∣h(x0)u_i = \frac{|x_i - x_0|}{h(x_0)}ui​=h(x0​)∣xi​−x0​∣​

其中h(x0)=∣x(k)−x0∣h(x_0) = |x_{(k)} - x_0|h(x0​)=∣x(k)​−x0​∣是到第kkk个最近邻居的距离（k=⌈αn⌉k = \lceil\alpha n\rceilk=⌈αn⌉，α\alphaα是平滑参数）。

三次立方权重函数的完整定义为：

W(u)={(1−u3)3if 0≤u≤10if u>1W(u) = \begin{cases} (1 - u^3)^3 & \text{if } 0 \leq u \leq 1 \\ 0 & \text{if } u > 1 \end{cases}W(u)={(1−u3)30​if 0≤u≤1if u>1​

这个函数具有优美的数学性质：

• W(0)=1W(0) = 1W(0)=1，在中心点权重最大
• W′(0)=0W'(0) = 0W′(0)=0，在原点处斜率为0，保证平滑性
• W(1)=0W(1) = 0W(1)=0，在边界处权重为0
• W′(1)=0W'(1) = 0W′(1)=0，边界处斜率也为0，避免突跳

对加权最小二乘问题求导并令其为0：

∂SSE∂β0=−2∑iwi(x0)[yi−β0−β1(xi−x0)]=0\frac{\partial \text{SSE}}{\partial \beta_0} = -2\sum_i w_i(x_0)[y_i - \beta_0 - \beta_1(x_i - x_0)] = 0∂β0​∂SSE​=−2i∑​wi​(x0​)[yi​−β0​−β1​(xi​−x0​)]=0

∂SSE∂β1=−2∑iwi(x0)(xi−x0)[yi−β0−β1(xi−x0)]=0\frac{\partial \text{SSE}}{\partial \beta_1} = -2\sum_i w_i(x_0)(x_i - x_0)[y_i - \beta_0 - \beta_1(x_i - x_0)] = 0∂β1​∂SSE​=−2i∑​wi​(x0​)(xi​−x0​)[yi​−β0​−β1​(xi​−x0​)]=0

记W=diag{w1(x0),w2(x0),…,wn(x0)}\mathbf{W} = \text{diag}\{w_1(x_0), w_2(x_0), \ldots, w_n(x_0)\}W=diag{w1​(x0​),w2​(x0​),…,wn​(x0​)}为权重对角矩阵，X\mathbf{X}X为设计矩阵：

X=[1(x1−x0)1(x2−x0)⋮⋮1(xn−x0)]\mathbf{X} = \begin{bmatrix} 1 & (x_1-x_0) \\ 1 & (x_2-x_0) \\ \vdots & \vdots \\ 1 & (x_n-x_0) \end{bmatrix}X=​11⋮1​(x1​−x0​)(x2​−x0​)⋮(xn​−x0​)​​

则正规方程组可写为矩阵形式：

(XTWX)β=XTWy(\mathbf{X}^T \mathbf{W} \mathbf{X})\boldsymbol{\beta} = \mathbf{X}^T \mathbf{W} \mathbf{y}(XTWX)β=XTWy

解得：β^=(XTWX)−1XTWy\hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{y}β^​=(XTWX)−1XTWy

局部拟合值为：f^(x0)=e1Tβ^=e1T(XTWX)−1XTWy\hat{f}(x_0) = \mathbf{e}_1^T \hat{\boldsymbol{\beta}} = \mathbf{e}_1^T (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{y}f^​(x0​)=e1T​β^​=e1T​(XTWX)−1XTWy

其中e1=[1,0]T\mathbf{e}_1 = [1, 0]^Te1​=[1,0]T。这表明LOESS实际上是原始观测值的线性组合，具有良好的统计性质。

稳健性权重的统计学理论

STL的稳健性来自于M-估计理论。定义稳健权重函数：

ρ(r)={r22if ∣r∣≤cc∣r∣−c22if ∣r∣>c\rho(r) = \begin{cases} \frac{r^2}{2} & \text{if } |r| \leq c \\ c|r| - \frac{c^2}{2} & \text{if } |r| > c \end{cases}ρ(r)={2r2​c∣r∣−2c2​​if ∣r∣≤cif ∣r∣>c​

其中ccc是调节参数。对应的影响函数为：

ψ(r)=ρ′(r)={rif ∣r∣≤cc⋅sign(r)if ∣r∣>c\psi(r) = \rho'(r) = \begin{cases} r & \text{if } |r| \leq c \\ c \cdot \text{sign}(r) & \text{if } |r| > c \end{cases}ψ(r)=ρ′(r)={rc⋅sign(r)​if ∣r∣≤cif ∣r∣>c​

权重函数定义为：w(r)=ψ(r)rw(r) = \frac{\psi(r)}{r}w(r)=rψ(r)​

STL使用双平方（bisquare）权重函数，它是Tukey的M-估计的特例：

w(r)={(1−(rc)2)2if ∣r∣≤c0if ∣r∣>cw(r) = \begin{cases} \left(1 - \left(\frac{r}{c}\right)^2\right)^2 & \text{if } |r| \leq c \\ 0 & \text{if } |r| > c \end{cases}w(r)=⎩⎨⎧​(1−(cr​)2)20​if ∣r∣≤cif ∣r∣>c​

其中c=6×median{∣ri∣}c = 6 \times \text{median}\{|r_i|\}c=6×median{∣ri​∣}是基于中位数绝对偏差的尺度估计。

这种权重的理论基础来自于最大化稳健似然函数：

L(μ)=∏i=1ng(riσ)L(\mu) = \prod_{i=1}^n g\left(\frac{r_i}{\sigma}\right)L(μ)=i=1∏n​g(σri​​)

其中g(u)∝exp⁡(−ρ(u))g(u) \propto \exp(-\rho(u))g(u)∝exp(−ρ(u))是稳健的密度函数。

低通滤波器的频域分析

STL中的低通滤波步骤涉及复杂的卷积运算。设季节长度为ppp，低通滤波器HHH由以下操作复合而成：

1. 移动平均MA§：HMA(Z)=1p∑j=0p−1ZjH_{MA}(Z) = \frac{1}{p} \sum_{j=0}^{p-1} Z^jHMA​(Z)=p1​∑j=0p−1​Zj
2. 移动平均MA§：再次应用相同滤波器
3. 移动平均MA(3)：H3(Z)=13(1+Z+Z2)H_3(Z) = \frac{1}{3}(1 + Z + Z^2)H3​(Z)=31​(1+Z+Z2)
4. LOESS平滑器：设为HL(Z)H_L(Z)HL​(Z)

复合滤波器的Z变换为：

H(Z)=HL(Z)⋅H3(Z)⋅HMA(Z)⋅HMA(Z)H(Z) = H_L(Z) \cdot H_3(Z) \cdot H_{MA}(Z) \cdot H_{MA}(Z)H(Z)=HL​(Z)⋅H3​(Z)⋅HMA​(Z)⋅HMA​(Z)

=HL(Z)⋅13(1+Z+Z2)⋅1p2(∑j=0p−1Zj)2= H_L(Z) \cdot \frac{1}{3}(1 + Z + Z^2) \cdot \frac{1}{p^2}\left(\sum_{j=0}^{p-1} Z^j\right)^2=HL​(Z)⋅31​(1+Z+Z2)⋅p21​(j=0∑p−1​Zj)2

其频率响应为：

H(eiω)=HL(eiω)⋅13∣1+eiω+ei2ω∣⋅1p2∣∑j=0p−1eijω∣2H(e^{i\omega}) = H_L(e^{i\omega}) \cdot \frac{1}{3}|1 + e^{i\omega} + e^{i2\omega}| \cdot \frac{1}{p^2}\left|\sum_{j=0}^{p-1} e^{ij\omega}\right|^2H(eiω)=HL​(eiω)⋅31​∣1+eiω+ei2ω∣⋅p21​​j=0∑p−1​eijω​2

对于季节频率ωk=2πkp\omega_k = \frac{2\pi k}{p}ωk​=p2πk​ (k=1,2,…,p−1k = 1, 2, \ldots, p-1k=1,2,…,p−1)，有：

∣∑j=0p−1eijωk∣2=0\left|\sum_{j=0}^{p-1} e^{ij\omega_k}\right|^2 = 0​j=0∑p−1​eijωk​​2=0

这确保了低通滤波器在所有季节频率处的增益为0，从而完全去除季节性成分。

收敛性分析与不动点理论

STL的内循环可以表述为不动点迭代。定义算子：

T:(S,T)↦(S′,T′)\mathcal{T}: (S, T) \mapsto (S', T')T:(S,T)↦(S′,T′)

其中S′S'S′和T′T'T′分别是通过STL步骤得到的新季节和趋势估计。

算子T\mathcal{T}T的构造如下：

1. 季节提取算子SopS_{op}Sop​：从去趋势序列中提取季节分量
2. 趋势提取算子TopT_{op}Top​：从去季节序列中提取趋势分量

则有：T(S,T)=(Sop(Y−T),Top(Y−Sop(Y−T)))\mathcal{T}(S, T) = (S_{op}(Y - T), T_{op}(Y - S_{op}(Y - T)))T(S,T)=(Sop​(Y−T),Top​(Y−Sop​(Y−T)))

STL的收敛性等价于证明存在不动点(S∗,T∗)(S^*, T^*)(S∗,T∗)使得：

T(S∗,T∗)=(S∗,T∗)\mathcal{T}(S^*, T^*) = (S^*, T^*)T(S∗,T∗)=(S∗,T∗)

在适当的函数空间中，可以证明T\mathcal{T}T是压缩映射。设∥⋅∥\|\cdot\|∥⋅∥为某个适当的范数，存在0<λ<10 < \lambda < 10<λ<1使得：

∥T(S1,T1)−T(S2,T2)∥≤λ∥(S1,T1)−(S2,T2)∥\|\mathcal{T}(S_1, T_1) - \mathcal{T}(S_2, T_2)\| \leq \lambda\|(S_1, T_1) - (S_2, T_2)\|∥T(S1​,T1​)−T(S2​,T2​)∥≤λ∥(S1​,T1​)−(S2​,T2​)∥

由Banach不动点定理，T\mathcal{T}T有唯一不动点，且迭代序列{(Sn,Tn)}\{(S_n, T_n)\}{(Sn​,Tn​)}收敛到该不动点。

最优性理论与变分方法

STL可以从变分的角度理解。考虑以下优化问题：

min⁡S,T∑iL(Yi−Si−Ti)+λ1R1(S)+λ2R2(T)\min_{S,T} \sum_i L(Y_i - S_i - T_i) + \lambda_1 R_1(S) + \lambda_2 R_2(T)S,Tmin​i∑​L(Yi​−Si​−Ti​)+λ1​R1​(S)+λ2​R2​(T)

其中：

• L(⋅)L(\cdot)L(⋅)是损失函数（通常是平方损失或稳健损失）
• R1(S)R_1(S)R1​(S)是季节分量的正则化项
• R2(T)R_2(T)R2​(T)是趋势分量的正则化项
• λ1,λ2\lambda_1, \lambda_2λ1​,λ2​是正则化参数

对于季节分量，自然的正则化是周期平滑性：

R1(S)=∑i∑j=1kωj(Si−Si+jp)2R_1(S) = \sum_i \sum_{j=1}^k \omega_j(S_i - S_{i+jp})^2R1​(S)=i∑​j=1∑k​ωj​(Si​−Si+jp​)2

其中kkk是考虑的滞后阶数，ωj\omega_jωj​是权重。

对于趋势分量，使用二阶差分的惩罚：

R2(T)=∑i(∇2Ti)2R_2(T) = \sum_i (\nabla^2 T_i)^2R2​(T)=i∑​(∇2Ti​)2

其中∇2Ti=Ti+1−2Ti+Ti−1\nabla^2 T_i = T_{i+1} - 2T_i + T_{i-1}∇2Ti​=Ti+1​−2Ti​+Ti−1​是二阶差分算子。

通过变分法，可以得到最优性条件：

∂J∂S=−L′(R)+λ1∂R1∂S=0\frac{\partial J}{\partial S} = -L'(R) + \lambda_1 \frac{\partial R_1}{\partial S} = 0∂S∂J​=−L′(R)+λ1​∂S∂R1​​=0

∂J∂T=−L′(R)+λ2∂R2∂T=0\frac{\partial J}{\partial T} = -L'(R) + \lambda_2 \frac{\partial R_2}{\partial T} = 0∂T∂J​=−L′(R)+λ2​∂T∂R2​​=0

这些方程的解正是STL算法近似求解的对象。

参数选择的理论基础

STL参数的选择涉及偏差-方差权衡的深层理论。对于LOESS平滑器，在点x0x_0x0​处的偏差和方差分别为：

Bias2[f^(x0)]≈h24(f′′(x0))2(∫−11u2W(u)du)2\text{Bias}^2[\hat{f}(x_0)] \approx \frac{h^2}{4}(f''(x_0))^2 \left(\int_{-1}^1 u^2 W(u) du\right)^2Bias2[f^​(x0​)]≈4h2​(f′′(x0​))2(∫−11​u2W(u)du)2

Var[f^(x0)]≈σ2nh∫−11W2(u)du\text{Var}[\hat{f}(x_0)] \approx \frac{\sigma^2}{nh} \int_{-1}^1 W^2(u) duVar[f^​(x0​)]≈nhσ2​∫−11​W2(u)du

其中hhh是带宽，nnn是样本量，σ2\sigma^2σ2是噪声方差。

最优带宽通过最小化均方误差MSE=Bias2+Variance\text{MSE} = \text{Bias}^2 + \text{Variance}MSE=Bias2+Variance得到：

hopt∝n−1/5h_{opt} \propto n^{-1/5}hopt​∝n−1/5

这是非参数回归中的经典结果。

对于STL，季节窗口nsn_sns​和趋势窗口ntn_tnt​的选择需要考虑季节和趋势的相对变化速度。如果季节模式变化较快，需要较小的nsn_sns​；如果趋势变化较平缓，可以使用较大的ntn_tnt​。

理论上，最优窗口大小满足：

ns∝(Var[S]∣∇2S∣2)1/5n_s \propto \left(\frac{\text{Var}[S]}{|\nabla^2 S|^2}\right)^{1/5}ns​∝(∣∇2S∣2Var[S]​)1/5

nt∝(Var[T]∣∇2T∣2)1/5n_t \propto \left(\frac{\text{Var}[T]}{|\nabla^2 T|^2}\right)^{1/5}nt​∝(∣∇2T∣2Var[T]​)1/5

其中∇2\nabla^2∇2表示二阶差分，反映了函数的局部曲率。

谱密度分析

从频域角度看，STL实现了自适应的谱估计。对于平稳时间序列，其谱密度为：

S(ω)=12π∑k=−∞∞γ(k)e−ikωS(\omega) = \frac{1}{2\pi} \sum_{k=-\infty}^{\infty} \gamma(k) e^{-ik\omega}S(ω)=2π1​k=−∞∑∞​γ(k)e−ikω

其中γ(k)\gamma(k)γ(k)是自协方差函数。

STL的分解等价于将谱密度分解为三部分：

S(ω)=ST(ω)+SS(ω)+SR(ω)S(\omega) = S_T(\omega) + S_S(\omega) + S_R(\omega)S(ω)=ST​(ω)+SS​(ω)+SR​(ω)

其中：

• ST(ω)S_T(\omega)ST​(ω)对应趋势分量的谱，主要集中在低频
• SS(ω)S_S(\omega)SS​(ω)对应季节分量的谱，表现为季节频率处的峰
• SR(ω)S_R(\omega)SR​(ω)对应余项的谱，理想情况下是白噪声谱

STL的滤波器设计确保了这种分解的有效性：低通滤波器抑制高频成分，季节滤波器增强季节频率，而稳健处理保证了余项的随机性。