时间序列分解理论

概述

时间序列分解（Time Series Decomposition）是时间序列分析的基础，其核心思想是将一个复杂的时间序列分解为若干个具有明确物理意义的成分。这种分解不仅有助于理解数据的生成机制，还能为预测建模提供特征工程的基础。¹

时间序列分解在实际应用中有广泛价值：

• 趋势分析：识别长期变化方向
• 季节性识别：发现周期性模式
• 异常检测：残差中识别离群点
• 预测增强：分别预测各成分再合成

---

一、时间序列分解的基本框架

1.1 加法与乘法分解

设原始时间序列为 $Y_t$，最基本的分解模型有两种：

加法分解（Additive Decomposition）：

$$Y_t=T_t+S_t+R_t$$

其中 $T_t$ 是趋势成分，$S_t$ 是季节成分，$R_t$ 是残差成分。

乘法分解（Multiplicative Decomposition）：

$$Y_t=T_t\times S_t\times R_t$$

1.2 如何选择分解方式

特征	加法分解	乘法分解
季节波动	不随均值变化	随均值成比例变化
适用场景	经济数据、平稳序列	增长型序列、销售数据
判别方法	观察季节波动幅度	季节波动随均值扩大

实践判断：当季节波动的绝对幅度不随时间序列水平变化时，使用加法分解；当季节波动的相对幅度保持恒定时，使用乘法分解。

1.3 经典分解方法

移动平均法（Moving Average）

经典的加法分解使用 $k$ 阶移动平均来估计趋势：

$$T_t=\frac{1}{k}\sum _{j=-q}^q{Y}_{t+j}$$

其中 $k=2q+1$ 为移动平均窗口长度。

中心移动平均：对于偶数阶移动平均，使用加权平均：

$$T_t=\frac{1}{2}\cdot \frac{1}{k}\sum _{j=-k/2}^{k/2-1}{Y}_{t+j}+\frac{1}{2}\cdot \frac{1}{k}\sum _{j=-k/2+1}^{k/2}{Y}_{t+j}$$

季节指数法

对于乘法分解，通过计算各周期的季节指数（Seasonal Index）来提取季节成分：

1. 计算所有周期的移动平均值（趋势）
2. 计算同期平均比率：$\frac{Y_t}{T_t}$
3. 按季节期求平均，得到各期季节指数 $S_i$
4. 归一化季节指数使其乘积为1

---

二、STL分解算法

2.1 算法原理

STL（Seasonal and Trend decomposition using Loess）是由 Cleveland 等人于1990年提出的经典分解算法。STL 具有强大的鲁棒性和灵活性，能处理任何类型的季节性数据。²

核心特点：

• 使用 Loess（局部加权回归）进行平滑
• 可以控制趋势和季节性成分的变化速率
• 对异常值具有鲁棒性
• 可以处理非整数的季节周期

2.2 数学表述

给定时间序列 $Y_t$，STL 将其分解为：

$$Y_t=T_t+S_t+R_t$$

其中 $T_t$ 是趋势成分，$S_t$ 是季节成分，$R_t$ 是残差成分。

Loess 局部加权回归：

对于点 $x_0$ 的估计，Loess 给出：

$${\hat{y}}_0=\sum _{i=1}^n{w}_i(x_0)\cdot y_i$$

其中权重函数 $w_i(x_0)$ 由距离核函数 $K$ 决定：

$$w_i(x_0)=K\left(\frac{|x_0-x_i|}{h}\right)\cdot {\nu }_i$$

这里 $h$ 是带宽参数，${\nu }_i$ 是鲁棒性权重（用于处理异常值）。

2.3 算法流程

STL 使用嵌套迭代方式，交替估计趋势和季节成分：

外层循环：处理异常值的鲁棒性权重

内层循环（对趋势和季节成分交替优化）：

步骤1：去季节化
    去趋势序列：Y_t - T_t^{(old)}
    
步骤2：周期子序列平滑
    对每个季节相位，将去趋势序列堆叠
    用Loess平滑，得到季节成分 S_t^{(new)}
    
步骤3：去季节化
    去季节序列：Y_t - S_t^{(new)}
    
步骤4：趋势平滑
    对去季节序列用Loess平滑，得到趋势成分 T_t^{(new)}
    
步骤5：更新残差
    R_t = Y_t - T_t^{(new)} - S_t^{(new)}

2.4 参数选择

参数	说明	选择建议
period	季节周期长度	根据数据频率：日数据(7/365)、月数据(12)、季度数据(4)
season	季节成分的平滑参数	奇数，范围$[1.5n_{season},10n_{season}]$
trend	趋势成分的平滑参数	奇数，通常是season的1.5-15倍
robust	是否使用鲁棒拟合	默认为True，处理异常值

2.5 Python实现

import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import STL
 
# 示例数据
np.random.seed(42)
n = 200
t = np.arange(n)
trend = 0.05 * t
seasonal = 10 * np.sin(2 * np.pi * t / 50)
noise = np.random.randn(n)
y = trend + seasonal + noise
 
# STL分解
stl = STL(y, period=50, robust=True)
result = stl.fit()
 
# 分解结果
trend_component = result.trend      # 趋势成分
seasonal_component = result.seasonal # 季节成分
residual_component = result.resid   # 残差成分
 
# 可视化
import matplotlib.pyplot as plt
fig, axes = plt.subplots(4, 1, figsize=(12, 8))
axes[0].plot(y, label='Original')
axes[0].legend()
axes[1].plot(trend_component, label='Trend')
axes[1].legend()
axes[2].plot(seasonal_component, label='Seasonal')
axes[2].legend()
axes[3].plot(residual_component, label='Residual')
axes[3].legend()
plt.tight_layout()
plt.show()

---

三、季节性强度的数学定义

3.1 趋势强度

趋势强度（Strength of Trend）衡量趋势成分相对于残差的变异程度：

$$F_T=\max \left(0,1-\frac{\text{Var}(R_t)}{\text{Var}(T_t+R_t)}\right)$$

• $F_T\approx 1$：序列有强趋势
• $F_T\approx 0$：几乎没有趋势

3.2 季节强度

季节强度（Strength of Seasonality）衡量季节成分相对于残差的变异程度：

$$F_S=\max \left(0,1-\frac{\text{Var}(R_t)}{\text{Var}(S_t+R_t)}\right)$$

• $F_S\approx 1$：序列有强季节性
• $F_S\approx 0$：几乎没有季节性

3.3 计算示例

def strength_measures(decomposition):
    """计算趋势强度和季节强度"""
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    
    # 去除NaN
    valid_mask = ~(np.isnan(trend) | np.isnan(residual))
    r = residual[valid_mask]
    t = trend[valid_mask]
    s = seasonal[valid_mask]
    
    # 趋势强度
    var_residual = np.var(r)
    var_trend_residual = np.var(t + r)
    F_T = max(0, 1 - var_residual / var_trend_residual)
    
    # 季节强度
    var_season_residual = np.var(s + r)
    F_S = max(0, 1 - var_residual / var_season_residual)
    
    return F_T, F_S
 
F_T, F_S = strength_measures(result)
print(f"趋势强度: {F_T:.3f}")
print(f"季节强度: {F_S:.3f}")

---

四、多周期分解：MSTL

4.1 问题背景

实际数据往往存在多个季节周期。例如：

• 每小时数据可能有日周期（24小时）和周周期（168小时）
• 交通流量数据可能有小时周期和工作日/周末周期

4.2 MSTL算法

MSTL（Multiple STL）由 Bandura 等人于2021年提出，能够自动处理多个季节周期。

算法思想：递归应用 STL，从最短周期到最长周期依次分解。

from statsmodels.tsa.seasonal import MSTL
 
# 多周期分解
# period=(24, 24*7) 表示日周期和周周期
mstl = MSTL(data, periods=(24, 24*7), robust=True)
result = mstl.fit()
 
# 结果包含多个季节成分
trend = result.trend
seasonal_24 = result.seasonal['24']  # 日周期
seasonal_168 = result.seasonal['168']  # 周周期
residual = result.resid

4.3 分解层次

对于包含日周期和周周期的小时数据，MSTL 的分解结构为：

$$Y_t=T_t+S_t^{(日)}+S_t^{(周)}+R_t$$

其中：

• $T_t$：趋势成分（捕捉长期变化）
• $S_t^{(日)}$：日周期成分（每天的重复模式）
• $S_t^{(周)}$：周周期成分（工作日/周末差异）
• $R_t$：残差成分

---

五、频域视角的分解

5.1 傅里叶级数表示

对于周期为 $s$ 的季节成分，可用傅里叶级数表示：

$$S_t=\sum _{k=1}^K\left[a_k\cos \left(\frac{2\pi kt}{s}\right)+b_k\sin \left(\frac{2\pi kt}{s}\right)\right]$$

其中 $K$ 是傅里叶项的阶数，控制季节性的复杂程度。

5.2 与时域分解的关系

方法	分解视角	优点	缺点
STL	时域平滑	鲁棒、可处理非线性	计算较慢
傅里叶级数	频域分析	数学优雅、可解释	假设固定周期
小波分解	时频联合	多尺度分析	实现复杂

5.3 傅里叶季节性实现

from scipy import signal
 
def fourier_seasonal(y, period, n_harmonics=5):
    """
    使用傅里叶级数估计季节成分
    """
    n = len(y)
    t = np.arange(n)
    
    # 构建傅里叶特征矩阵
    X = np.zeros((n, 2 * n_harmonics))
    for k in range(1, n_harmonics + 1):
        X[:, 2*(k-1)] = np.cos(2 * np.pi * k * t / period)
        X[:, 2*(k-1) + 1] = np.sin(2 * np.pi * k * t / period)
    
    # 最小二乘估计
    from numpy.linalg import lstsq
    coeffs, _, _, _ = lstsq(X, y, rcond=None)
    
    # 季节成分预测
    seasonal = X @ coeffs
    
    return seasonal, coeffs
 
# 示例
y_deseasonal = y - trend  # 先去趋势
seasonal, coeffs = fourier_seasonal(y_deseasonal, period=50, n_harmonics=3)

---

六、实践指南

6.1 分解方法选择

graph TD
    A[选择分解方法] --> B{数据特点}
    B -->|单一周期| C{异常值?}
    C -->|是| D[STL robust=True]
    C -->|否| E[STL]
    B -->|多周期| F[MSMSTL]
    B -->|固定周期| G[傅里叶级数]

6.2 参数自动选择

def auto_period_detection(y, fs=1):
    """
    自动检测季节周期
    fs: 采样频率
    """
    from scipy import signal
    
    # 计算自相关
    n = min(len(y), 1000)
    autocorr = np.correlate(y[:n], y[:n], mode='full')
    autocorr = autocorr[n-1:]
    autocorr /= autocorr[0]
    
    # 找峰值（排除lag=0）
    peaks, _ = signal.find_peaks(autocorr[1:], height=0.1)
    
    if len(peaks) > 0:
        return peaks[0] + 1  # 返回第一个显著周期
    return None

6.3 常见问题

问题	原因	解决方案
趋势不连续	period选择不当	增大trend参数
季节性消失	季节成分过平滑	减小season参数
残差有模式	分解不完全	检查多周期或非线性趋势

---

七、相关算法与扩展

7.1 相关算法

算法	特点	适用场景
X-13-ARIMA-SEATS	官方统计方法	政府经济数据
SEATS	基于ARIMA的分解	季节调整
HP滤波器	Hodrick-Prescott滤波	经济趋势-周期分解
B-K滤波器	Baxter-King滤波	宏观经济分析

7.2 与预测的结合

分解后的成分可以分别预测再合成：

def decompose_forecast(y, forecast_horizon):
    """
    分解-预测-合成流程
    """
    # 分解
    stl = STL(y, period=50, robust=True)
    result = stl.fit()
    
    # 预测各成分
    # 趋势：线性或多项式外推
    trend_model = np.polyfit(range(len(y)), result.trend, deg=1)
    trend_forecast = np.polyval(trend_model, 
                                range(len(y), len(y) + forecast_horizon))
    
    # 季节性：周期延拓
    seasonal_forecast = np.tile(result.seasonal[:50], 
                               forecast_horizon // 50 + 1)[:forecast_horizon]
    
    # 残差：假设均值为0或用简单模型
    residual_forecast = np.zeros(forecast_horizon)
    
    # 合成
    forecast = trend_forecast + seasonal_forecast + residual_forecast
    
    return forecast

---

八、参考

---

相关阅读

• 高级平稳性理论
• 时序频域分析
• N-BEATS与N-HiTS架构
• 循环神经网络
• 长短期记忆网络

Footnotes

1. Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice, 3rd edition. OTexts. ↩

2. Cleveland, R.B., Cleveland, W.S., McRae, J.E., & Terpenning, I. (1990). STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Journal of Official Statistics, 6(1), 3-73. ↩