Python案例如何用Statsmodels做时间序列分解

wen python案例 1

本文目录导读:

Python案例如何用Statsmodels做时间序列分解

  1. 基本安装和导入
  2. 加法模型分解(最常见)
  3. 乘法模型分解
  4. 实际数据案例:航空公司乘客数据
  5. 使用STL分解(更稳健的方法)
  6. 异常检测应用
  7. 季节性调整
  8. 实用的完整示例
  9. 关键要点

我来详细介绍如何使用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()

关键要点

  1. 选择合适的模型

    • 加法模型(additive):季节波动幅度相对固定
    • 乘法模型(multiplicative):季节波动随趋势变化
  2. 确定周期

    • 月度数据通常为12
    • 季度数据通常为4
    • 周数据通常为7
  3. STL vs 经典分解

    • STL:更稳健,能处理非线性趋势
    • 经典分解:简单快速,适合规则数据

这些方法可以帮助您理解时间序列的组成成分,并进行预测、异常检测等进一步分析。

抱歉,评论功能暂时关闭!