本文目录导读:

我来详细介绍如何使用Statsmodels进行时间序列分解。
基本安装和导入
# 安装必要的库 # pip install pandas numpy statsmodels matplotlib import pandas as pd import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.tsa.seasonal import seasonal_decompose from statsmodels.tsa.stattools import adfuller
加法模型分解(最常见)
# 创建示例数据:月度销售数据
np.random.seed(42)
dates = pd.date_range(start='2020-01-01', periods=36, freq='M')
# 生成合成时间序列:趋势 + 季节 + 噪声
trend = np.linspace(100, 200, 36) # 上升趋势
seasonal = 10 * np.sin(2 * np.pi * np.arange(36) / 12) # 12个月周期
noise = np.random.normal(0, 5, 36) # 随机噪声
data = trend + seasonal + noise
ts = pd.Series(data, index=dates, name='Sales')
# 执行时间序列分解
decomposition = seasonal_decompose(ts, model='additive', period=12)
# 获取分解结果
trend_component = decomposition.trend
seasonal_component = decomposition.seasonal
residual_component = decomposition.resid
print("时间序列分解结果:")
print(f"趋势成分: {trend_component.head()}")
print(f"季节成分: {seasonal_component.head()}")
print(f"残差成分: {residual_component.head()}")
# 可视化
fig, axes = plt.subplots(4, 1, figsize=(12, 10))
ts.plot(ax=axes[0], title='原始数据')
trend_component.plot(ax=axes[1], title='趋势成分')
seasonal_component.plot(ax=axes[2], title='季节成分')
residual_component.plot(ax=axes[3], title='残差成分')
plt.tight_layout()
plt.show()
乘法模型分解
# 创建乘法模型的数据 np.random.seed(42) trend_mul = np.linspace(100, 200, 36) seasonal_mul = 1 + 0.3 * np.sin(2 * np.pi * np.arange(36) / 12) noise_mul = np.random.normal(1, 0.05, 36) data_mul = trend_mul * seasonal_mul * noise_mul ts_mul = pd.Series(data_mul, index=dates, name='Sales_Multiplicative') # 乘法模型分解 decomposition_mul = seasonal_decompose(ts_mul, model='multiplicative', period=12) # 可视化 fig, axes = plt.subplots(4, 1, figsize=(12, 10)) ts_mul.plot(ax=axes[0], title='原始数据(乘法模型)') decomposition_mul.trend.plot(ax=axes[1], title='趋势成分') decomposition_mul.seasonal.plot(ax=axes[2], title='季节成分') decomposition_mul.resid.plot(ax=axes[3], title='残差成分') plt.tight_layout() plt.show()
实际数据案例:航空公司乘客数据
# 加载经典数据集
from statsmodels.datasets import get_rdataset
# 使用航空公司乘客数据集
air_passengers = get_rdataset('AirPassengers').data
air_passengers.index = pd.date_range(start='1949-01-01', periods=144, freq='M')
# 重命名列
ts_air = air_passengers['value']
# 执行分解
decomposition_air = seasonal_decompose(ts_air, model='multiplicative', period=12)
# 打印季节性因子
print("季节性因子(前12个月):")
print(decomposition_air.seasonal.head(12))
# 去季节化
detrended = ts_air / decomposition_air.trend
deseasonalized = ts_air / decomposition_air.seasonal
# 可视化
fig, axes = plt.subplots(3, 2, figsize=(15, 10))
# 原始数据
ts_air.plot(ax=axes[0, 0], title='原始航空公司乘客数据')
axes[0, 0].set_ylabel('乘客数')
# 分解成分
decomposition_air.trend.plot(ax=axes[0, 1], title='趋势成分')
axes[0, 1].set_ylabel('乘客数')
decomposition_air.seasonal.plot(ax=axes[1, 0], title='季节成分')
axes[1, 0].set_ylabel('季节因子')
decomposition_air.resid.plot(ax=axes[1, 1], title='残差成分')
axes[1, 1].set_ylabel('残差')
# 去趋势和去季节化后的数据
detrended.plot(ax=axes[2, 0], title='去趋势数据')
axes[2, 0].set_ylabel('乘客数')
deseasonalized.plot(ax=axes[2, 1], title='去季节化数据')
axes[2, 1].set_ylabel('乘客数')
plt.tight_layout()
plt.show()
使用STL分解(更稳健的方法)
from statsmodels.tsa.seasonal import STL
# STL分解(适用于更复杂的季节性模式)
stl = STL(ts_air, period=12, robust=True)
result = stl.fit()
# 获取结果
trend_stl = result.trend
seasonal_stl = result.seasonal
residual_stl = result.resid
# 可视化STL分解结果
fig, axes = plt.subplots(4, 1, figsize=(12, 10))
ts_air.plot(ax=axes[0], title='原始数据')
trend_stl.plot(ax=axes[1], title='STL趋势成分')
seasonal_stl.plot(ax=axes[2], title='STL季节成分')
residual_stl.plot(ax=axes[3], title='STL残差成分')
plt.tight_layout()
plt.show()
# 分析残差统计
print("残差统计信息:")
print(f"均值: {residual_stl.mean():.4f}")
print(f"标准差: {residual_stl.std():.4f}")
print(f"偏度: {residual_stl.skew():.4f}")
print(f"峰度: {residual_stl.kurtosis():.4f}")
异常检测应用
# 使用分解结果进行异常检测
def detect_anomalies(series, decomposition, threshold=3):
"""
基于时间序列分解的异常检测
"""
residuals = decomposition.resid
residual_std = residuals.std()
residual_mean = residuals.mean()
# 计算Z-score
z_scores = (residuals - residual_mean) / residual_std
# 标记异常点
anomalies = series[abs(z_scores) > threshold]
return anomalies, z_scores
# 检测异常
anomalies, z_scores = detect_anomalies(ts_air, decomposition_air, threshold=2.5)
print(f"检测到的异常点数量: {len(anomalies)}")
print("异常点详情:")
for date, value in anomalies.items():
print(f"{date.date()}: {value:.0f} (Z-score: {z_scores[date]:.2f})")
# 可视化异常点
fig, ax = plt.subplots(figsize=(12, 6))
ts_air.plot(ax=ax, label='原始数据')
ax.scatter(anomalies.index, anomalies, color='red', s=50, label='异常点')
ax.set_title('时间序列异常检测')
ax.legend()
plt.show()
季节性调整
# 季节性调整(去季节化)
def seasonal_adjustment(series, decomposition, model='additive'):
"""
季节性调整函数
"""
if model == 'additive':
adjusted = series - decomposition.seasonal
elif model == 'multiplicative':
adjusted = series / decomposition.seasonal
return adjusted
# 对乘法模型进行季节性调整
seasonally_adjusted = seasonal_adjustment(ts_air, decomposition_air, model='multiplicative')
# 对比调整前后的数据
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
ts_air.plot(ax=axes[0], title='原始数据')
axes[0].set_ylabel('乘客数')
seasonally_adjusted.plot(ax=axes[1], title='季节性调整后数据')
axes[1].set_ylabel('乘客数')
plt.tight_layout()
plt.show()
# 计算调整前后的统计数据
print("季节性调整效果对比:")
print(f"原始数据标准差: {ts_air.std():.2f}")
print(f"调整后标准差: {seasonally_adjusted.std():.2f}")
print(f"标准差异动减少: {(1 - seasonally_adjusted.std()/ts_air.std())*100:.1f}%")
实用的完整示例
class TimeSeriesDecomposer:
"""
时间序列分解封装类
"""
def __init__(self, series, period=None):
self.series = series
self.period = period or self._detect_period()
self.decomposition = None
self.stl_result = None
def _detect_period(self):
"""
自动检测周期(简单实现)
"""
if len(self.series) >= 24:
return 12 # 月度数据假设12个月周期
else:
return 4 # 季度数据
def decompose(self, model='additive'):
"""
执行分解
"""
self.decomposition = seasonal_decompose(
self.series,
model=model,
period=self.period
)
return self.decomposition
def stl_decompose(self, robust=True):
"""
STL分解
"""
stl = STL(self.series, period=self.period, robust=robust)
self.stl_result = stl.fit()
return self.stl_result
def get_components(self):
"""
获取所有成分
"""
if self.decomposition:
return {
'observed': self.series,
'trend': self.decomposition.trend,
'seasonal': self.decomposition.seasonal,
'residual': self.decomposition.resid
}
return None
def plot_decomposition(self, figsize=(12, 10)):
"""
可视化分解结果
"""
if not self.decomposition:
self.decompose()
fig, axes = plt.subplots(4, 1, figsize=figsize)
self.series.plot(ax=axes[0], title='原始数据')
self.decomposition.trend.plot(ax=axes[1], title='趋势')
self.decomposition.seasonal.plot(ax=axes[2], title='季节性')
self.decomposition.resid.plot(ax=axes[3], title='残差')
plt.tight_layout()
return fig
# 使用示例
decomposer = TimeSeriesDecomposer(ts_air, period=12)
components = decomposer.decompose(model='multiplicative')
decomposer.plot_decomposition()
plt.show()
关键要点
-
选择合适的模型:
- 加法模型(additive):季节波动幅度相对固定
- 乘法模型(multiplicative):季节波动随趋势变化
-
确定周期:
- 月度数据通常为12
- 季度数据通常为4
- 周数据通常为7
-
STL vs 经典分解:
- STL:更稳健,能处理非线性趋势
- 经典分解:简单快速,适合规则数据
这些方法可以帮助您理解时间序列的组成成分,并进行预测、异常检测等进一步分析。