时序数据预处理

时序数据预处理完全指南

时间序列分析是数据科学中一个重要且广泛应用的领域。无论是股票价格预测、天气预报，还是用户行为分析，时序数据都扮演着关键角色。本文将系统介绍时序数据预处理的各个环节，帮助初学者建立完整的知识体系。

1. 什么是时序数据

时序数据（Time Series Data）是按照时间顺序排列的数据序列，其中每个数据点都与特定的时间戳相关联。与普通数据集不同，时序数据的一个显著特征是观测值之间存在时间依赖性——当前的值往往受到过去值的影响。

时序数据具有以下几个关键特征：

• 时间顺序性：数据点按照时间先后顺序排列，这种顺序不能随意改变
• 时间依赖性：相邻时间点的数据往往存在相关性
• 趋势性：数据可能表现出长期上升或下降的趋势
• 季节性：数据可能呈现周期性的模式
• 随机性：数据中包含无法预测的随机波动

让我们通过航空乘客数据来直观理解时序数据：

import pandas as pd
import matplotlib.pyplot as plt# 读取航空乘客数据
passengers = pd.read_csv("./data/AirPassengers.csv")
passengers['date'] = pd.to_datetime(passengers['date'])# 绘制时序图
plt.figure(figsize=(10, 6))
plt.plot(passengers["date"], passengers["value"])
plt.title("International Air Passengers Over Time")
plt.xlabel("Time")
plt.ylabel("Number of Passengers (in Thousands)")
plt.show()

从这个例子可以看出，航空乘客数据展现了明显的上升趋势和季节性模式，这是典型的时序数据特征。

2. 时序数据的重要统计指标

2.1 均值（Mean）

均值是时序数据最基本的统计量，它反映了数据的中心位置。对于时序数据，我们常常关注的不仅是整体均值，还包括局部时间窗口内的均值。

# 计算整体均值
overall_mean = passengers['value'].mean()
print(f"整体均值: {overall_mean:.2f}")# 计算滚动均值（12个月窗口）
passengers['rolling_mean'] = passengers['value'].rolling(window=12).mean()
# 整体均值: 280.30

2.2 方差（Variance）

方差衡量数据的离散程度。在时序分析中，方差的稳定性是判断数据平稳性的重要指标之一。如果方差随时间变化较大，可能需要进行数据变换。

# 计算整体方差
overall_variance = passengers['value'].var()
print(f"整体方差: {overall_variance:.2f}")
# 整体方差: 14391.92
# 计算滚动方差
passengers['rolling_var'] = passengers['value'].rolling(window=12).var()

2.3 协方差（Covariance）

协方差衡量两个变量之间的线性相关程度。在时序分析中，我们经常计算时间序列与其滞后值之间的协方差，这就是自协方差函数的基础。

# 使用欧洲股市数据展示多变量时序的协方差
stock = pd.read_csv("./data/EuStockMarkets.csv")# 计算不同股市指数之间的协方差
cov_matrix = stock[['DAX', 'SMI', 'CAC', 'FTSE']].cov()
print("股市指数协方差矩阵:")
print(cov_matrix)# 股市指数协方差矩阵:
#                DAX           SMI            CAC          FTSE
# DAX   1.176775e+06  1.788080e+06  608260.100153  1.033234e+06
# SMI   1.788080e+06  2.765657e+06  913749.179754  1.608010e+06
# CAC   6.082601e+05  9.137492e+05  336764.568483  5.190355e+05
# FTSE  1.033234e+06  1.608010e+06  519035.498220  9.539732e+05

2.4 数据平稳性（Stationarity）

平稳性是时序分析中的核心概念，它决定了我们可以使用哪些模型和分析方法。

"""
理解平稳性
"""
# %%  读取数据
passengers = pd.read_csv("./data/AirPassengers.csv")
# 将日期列转换为 datetime 格式
passengers['date'] = pd.to_datetime(passengers['date'])
# %%  绘制折线图plt.figure(figsize=(8, 5))
plt.plot(passengers["date"], passengers["value"])# 设置 x 轴刻度，只显示每隔24个点的年份
xticks = passengers["date"][::24]
plt.xticks(xticks, [tick.year for tick in xticks])# 添加图表标题和标签
plt.title("International Air Passengers Over Time")
plt.xlabel("Time")
plt.ylabel("Number of Passengers (in Thousands)")
plt.show()

严格平稳（Strict Stationarity）

如果时间序列的所有统计特性都不随时间变化，我们称之为严格平稳。这意味着任意时间点的联合分布函数与时间起点无关，只与时间间隔有关。严格平稳的条件非常苛刻，现实中很少有数据能够满足。

宽平稳（Weak Stationarity）

宽平稳也称为二阶平稳或协方差平稳，要求满足以下三个条件：

1. 均值不随时间变化：$E[X(t)]=\mu$（常数）
2. 方差不随时间变化：$Var[X(t)]={\sigma }^2$（常数）
3. 协方差只与时间间隔有关：$Cov[X(t),X(t+h)]=\gamma (h)$

大多数时序分析方法只要求数据满足宽平稳条件即可。

模型适用性

平稳性是许多经典时序模型（如$ARMA$模型）的基本假设。非平稳数据可能导致：

• 模型参数估计不一致
• 预测效果不佳
• 统计检验失效

准确预测性

平稳时序的统计特性保持稳定，使得我们可以用历史数据的模式来预测未来。如果数据非平稳，过去的模式可能在未来不再适用。

简化分析

平稳性大大简化了时序分析的复杂度。对于平稳序列，我们只需要估计有限的参数（如均值、方差、自相关函数），就能够刻画整个序列的特征。

2.5 伪相关性（Spurious Correlation）

伪相关是指两个本质上没有因果关系的时间序列，由于共同的趋势或其他因素，表现出虚假的相关性。这在非平稳时序中特别常见。

例如，如果两个时间序列都有上升趋势，即使它们之间没有任何实际联系，相关系数也可能很高。这就是为什么在分析时序相关性之前，通常需要先去除趋势或进行差分处理。

2.6 协整（Cointegration）

协整是描述非平稳时间序列之间长期均衡关系的概念。如果两个或多个非平稳序列的某个线性组合是平稳的，我们说这些序列是协整的。

协整关系在经济金融领域特别重要。例如，股票价格和股息、利率和通货膨胀率等经济变量之间往往存在协整关系。

2.7 时间依赖（Temporal Dependence）

时间依赖是时序数据的本质特征，指当前值与过去值之间的相关性。理解时间依赖对于建模和预测至关重要。

我们可以通过自相关函数（ACF）来量化时间依赖：

import statsmodels.api as sm# 绘制自相关函数
fig, ax = plt.subplots(figsize=(10, 6))
sm.graphics.tsa.plot_acf(passengers['value'], lags=40, ax=ax)
plt.title("Autocorrelation Function of Air Passengers Data")
plt.show()

3. 滑动平均

滑动平均（Moving Average）是时序数据处理中最基本也是最重要的技术之一。它是对数据集的一个时间范围的点进行平均，并随着时间的推移，连续移动这个取样的窗口，每一个新的平均值都是通过排除时间窗口的第一个数据点，并包含下一个数据点来计算的。

参数

• 窗长：平均值计算中，包括历史数据的数量。窗长越大，平滑效果越显著，减少随机波动。但可能导致滞后效应，即平均值对真实数据变化的响应太慢。
• 步长：滑动时窗口一定的数据点数。步长越小，平滑效果越显著，但计算量较大。

3.1 简单滑动平均

简单滑动平均对窗口内的所有数据点赋予相同的权重：

# 使用失业率数据演示滑动平均
unemp = pd.read_csv("data/unemployment_rate.csv")
unemp['DATE'] = pd.to_datetime(unemp['DATE'])# 计算不同窗口的滑动平均
unemp['MA_3'] = unemp['UNRATE'].rolling(window=3).mean()
unemp['MA_12'] = unemp['UNRATE'].rolling(window=12).mean()plt.figure(figsize=(12, 6))
plt.plot(unemp['DATE'], unemp['UNRATE'], label='Original', alpha=0.5)
plt.plot(unemp['DATE'], unemp['MA_3'], label='MA(3)', linewidth=2)
plt.plot(unemp['DATE'], unemp['MA_12'], label='MA(12)', linewidth=2)
plt.legend()
plt.title('Unemployment Rate with Moving Averages')
plt.show()

3.2 指数加权滑动平均

指数加权滑动平均（EWMA）给予近期数据更高的权重，权重按指数递减：

# 计算指数加权滑动平均
unemp['EWMA'] = unemp['UNRATE'].ewm(span=12, adjust=False).mean()

3.3 窗口类型

滑动窗口有多种类型，每种都有其适用场景：

# 演示不同的窗口类型
np.random.seed(0)
x = np.random.normal(loc=0, scale=10, size=100) + np.arange(1, 101)
x_series = pd.Series(x)# 滑动窗口（向后看）
rolling_mean = x_series.rolling(window=20).mean()# 居中窗口
center_mean = x_series.rolling(window=20, center=True).mean()# 扩展窗口（累积）
expanding_mean = x_series.expanding().mean()

4. 缺失值处理

时序数据中的缺失值处理比普通数据更具挑战性，因为我们需要考虑数据的时间连续性和时间依赖性。

4.1 缺失值检测

首先需要识别数据中的缺失模式：

# 创建带缺失值的时序数据
dates = pd.date_range(start=unemp['DATE'].min(), end=unemp['DATE'].max(), freq='MS')
unemp_full = unemp.set_index('DATE').reindex(dates)# 检测缺失值
missing_mask = unemp_full['UNRATE'].isna()
print(f"缺失值数量: {missing_mask.sum()}")
print(f"缺失值比例: {missing_mask.sum() / len(unemp_full):.2%}")

4.2 删除法

删除法是最简单的处理方式，但可能破坏时间连续性：

# 删除含有缺失值的行
unemp_dropna = unemp_full.dropna()

这种方法只适用于缺失值很少且随机分布的情况。

4.3 填充法

均值填充

使用整体均值或局部均值填充：

# 整体均值填充
mean_value = unemp_full['UNRATE'].mean()
unemp_mean_filled = unemp_full.fillna(mean_value)# 滚动均值填充
unemp_full['rolling_mean_filled'] = unemp_full['UNRATE'].fillna(unemp_full['UNRATE'].rolling(window=6, min_periods=1).mean()
)

中位数填充

中位数对异常值更稳健：

median_value = unemp_full['UNRATE'].median()
unemp_median_filled = unemp_full.fillna(median_value)

众数填充

适用于离散型时序数据：

mode_value = unemp_full['UNRATE'].mode()[0]
unemp_mode_filled = unemp_full.fillna(mode_value)

4.4 前向填充法

前向填充使用最近的历史值填充缺失值，特别适合时序数据：

"""
失业率时间序列数据缺失处理
"""# %%  导入库
import pandas as pd
import numpy as np
# import utils# %%  读取数据
unemp = pd.read_csv("data/unemployment_rate.csv")# 将日期列转换为 datetime 格式
unemp['DATE'] = pd.to_datetime(unemp['DATE'])# 设置日期列为索引
unemp.set_index('DATE', inplace=True)# %%  缺失数据生成# 生成一个数据随机缺失的数据集
rand_unemp_idx = np.random.choice(unemp.index, int(0.2 * len(unemp)), replace=False)
rand_unemp = unemp.drop(rand_unemp_idx)# 生成一个数据集，失业率高的数据更有可能缺失
high_unemp_idx = unemp[unemp['UNRATE'] > 8].index
num_to_select = int(0.2 * len(high_unemp_idx))
high_unemp_idx = np.random.choice(high_unemp_idx, num_to_select, replace=False)
bias_unemp = unemp.drop(high_unemp_idx)# %%  检测并报告缺失值# 生成所有日期的序列，从第一个日期到最后一个日期，按月递增，且固定为每个月的1号
all_dates = pd.date_range(start=unemp.index[0], end=unemp.index[-1], freq='MS')# 将 rand_unemp 和 bias_unemp 重新索引到 all_dates
rand_unemp_aligned = rand_unemp.reindex(all_dates)
bias_unemp_aligned = bias_unemp.reindex(all_dates)# 检测 rand_unemp 中的缺失值
rand_unemp_missing = rand_unemp_aligned['UNRATE'].isna()# 检测 bias_unemp 中的缺失值
bias_unemp_missing = bias_unemp_aligned['UNRATE'].isna()# 打印结果
print("Rand Unemp Missing Values:")
print(rand_unemp_missing[rand_unemp_missing].index)print("\nBias Unemp Missing Values:")
print(bias_unemp_missing[bias_unemp_missing].index)# %%  前向填充法rand_unemp_filled = rand_unemp_aligned.fillna(method='ffill')
bias_unemp_filled = bias_unemp_aligned.fillna(method='ffill')plot_forward_fill_ts(unemp, rand_unemp_filled, rand_unemp_missing, bias_unemp_filled, bias_unemp_missing)
plot_forward_fill_comparison(unemp, rand_unemp_filled, bias_unemp_filled)

4.5 移动平均填充

使用移动平均值填充，能够保持局部趋势：

# %%  移动平均法
# 定义一个自定义函数，用于移动窗口的操作
def custom_rolling_mean(x):if not np.isnan(x).any():return x[-1]else:return np.nanmean(x)# 在 rand_unemp_aligned 中计算没有前瞻的移动平均
rand_unemp_ma = (rand_unemp_aligned['UNRATE'].rolling(window=3, min_periods=1).apply(custom_rolling_mean, raw=True)
)# 在 bias_unemp_aligned 中计算没有前瞻的移动平均
bias_unemp_ma = (bias_unemp_aligned['UNRATE'].rolling(window=5, min_periods=1).apply(custom_rolling_mean, raw=True)
)
rand_unemp_ma = rand_unemp_ma.to_frame(name='UNRATE')
bias_unemp_ma = bias_unemp_ma.to_frame(name='UNRATE')
rand_unemp_ma.reset_index(inplace=True)
bias_unemp_ma.reset_index(inplace=True)plot_ma_fill_ts(unemp, rand_unemp_ma, bias_unemp_ma)
plot_ma_fill_comparison(unemp, rand_unemp_ma, bias_unemp_ma)

4.6 插值法

线性插值

线性插值假设缺失值在已知值之间呈线性变化：

# %%  插值法——线性插值法rand_unemp_impute_linear = rand_unemp_aligned['UNRATE'].interpolate(method='linear')
bias_unemp_impute_linear = bias_unemp_aligned['UNRATE'].interpolate(method='linear')
rand_unemp_impute_linear = rand_unemp_impute_linear.to_frame(name='UNRATE')
bias_unemp_impute_linear = bias_unemp_impute_linear.to_frame(name='UNRATE')
rand_unemp_impute_linear.reset_index(inplace=True)
bias_unemp_impute_linear.reset_index(inplace=True)plot_impute_linear_ts(unemp, rand_unemp_impute_linear, bias_unemp_impute_linear)
plot_impute_linear_comparison(unemp, rand_unemp_impute_linear, bias_unemp_impute_linear)

样条插值

样条插值能够产生更平滑的曲线：

# 三次样条插值
# %%  插值法——样条插值法rand_unemp_impute_spline = rand_unemp_aligned['UNRATE'].interpolate(method='spline', order=3)
bias_unemp_impute_spline = bias_unemp_aligned['UNRATE'].interpolate(method='spline', order=3)
rand_unemp_impute_spline = rand_unemp_impute_spline.to_frame(name='UNRATE')
bias_unemp_impute_spline = bias_unemp_impute_spline.to_frame(name='UNRATE')
rand_unemp_impute_spline.reset_index(inplace=True)
bias_unemp_impute_spline.reset_index(inplace=True)plot_impute_spline_ts(unemp, rand_unemp_impute_spline, bias_unemp_impute_spline)
plot_impute_spline_comparison(unemp, rand_unemp_impute_spline, bias_unemp_impute_spline)

4.7 机器学习方法

KNN填充

使用K最近邻算法，基于相似时间点的值进行填充：

from sklearn.impute import KNNImputer# 准备特征：使用时间特征
unemp_full['month'] = unemp_full.index.month
unemp_full['year'] = unemp_full.index.year# KNN填充
imputer = KNNImputer(n_neighbors=5)
features = unemp_full[['year', 'month', 'UNRATE']].values
filled_data = imputer.fit_transform(features)
unemp_full['UNRATE_knn'] = filled_data[:, 2]

5. 数据平滑

数据平滑技术用于减少时序数据中的噪声，帮助识别潜在的趋势和模式。

5.1 移动平均平滑

移动平均是最常用的平滑方法，取相邻观测值的平均值来平滑数据。

# 使用失业率数据演示不同的平滑方法
unemp = pd.read_csv("data/unemployment_rate.csv")
unemp['DATE'] = pd.to_datetime(unemp['DATE'])# 简单移动平均
unemp['MA'] = unemp['UNRATE'].rolling(window=12, min_periods=1).mean()

移动平均的窗口大小选择很重要：

• 窗口太小：平滑效果不明显
• 窗口太大：可能丢失重要的短期变化

5.2 指数平滑

指数平滑给予近期数据更高的权重，能够更快地响应数据变化：

# 指数平滑
unemp['ES'] = unemp['UNRATE'].ewm(span=12, adjust=False).mean()# 不同平滑参数的比较
unemp['ES_fast'] = unemp['UNRATE'].ewm(span=6, adjust=False).mean()
unemp['ES_slow'] = unemp['UNRATE'].ewm(span=24, adjust=False).mean()

5.3 Loess平滑

Loess（Locally Weighted Scatterplot Smoothing）是一种局部多项式回归平滑方法，能够自适应地调整平滑程度：

下面是几种平滑方式对unemployment_rate.csv失业率数据集处理后的效果情况。

"""
失业率时间序列数据平滑处理
"""# %%  读取数据
unemp = pd.read_csv("data/unemployment_rate.csv")# 将日期列转换为 datetime 格式
unemp['DATE'] = pd.to_datetime(unemp['DATE'])# %%  平滑
# 移动平均平滑
unemp['MA'] = unemp['UNRATE'].rolling(window=12, min_periods=1).mean()# 指数平滑
unemp['ES'] = unemp['UNRATE'].ewm(span=12, adjust=False).mean()# Loess 平滑
loess_smoothed = lowess(unemp['UNRATE'], unemp['DATE'], frac=0.1)
unemp['Loess'] = np.interp(unemp['DATE'], pd.to_datetime(loess_smoothed[:, 0]), loess_smoothed[:, 1])# 绘制结果
plt.figure(figsize=(14, 7))# 原始数据
plt.plot(unemp['DATE'], unemp['UNRATE'], label='Original Data', color='blue', alpha=0.5)# 移动平均平滑
plt.plot(unemp['DATE'], unemp['MA'], label='Moving Average (12 months)', color='green', linewidth=2)# 指数平滑
plt.plot(unemp['DATE'], unemp['ES'], label='Exponential Smoothing (12 months)', color='red', linewidth=2)# Loess 平滑
plt.plot(unemp['DATE'], unemp['Loess'], label='Loess Smoothing', color='orange', linewidth=2)# 图例和标题
plt.legend()
plt.xlabel('Date')
plt.ylabel('Unemployment Rate')
plt.title('Comparison of Smoothing Methods on Unemployment Rate')
plt.grid(True)# 显示图表
plt.show()

5.4 平滑参数选择

选择合适的平滑参数需要权衡偏差和方差：

• 过度平滑：可能丢失重要信息
• 平滑不足：噪声仍然明显

可以使用交叉验证或信息准则来选择最优参数。

6. 移动平均模型（MA）

移动平均模型是时序分析中的基础模型之一，它假设当前值是过去若干期随机扰动的线性组合。

6.1 MA模型的数学表达

$MA(q)$模型的数学形式为：
$$X_t=\mu +{\epsilon }_t+{\theta }_1\ast {\epsilon }_{t-1}+{\theta }_2\ast {\epsilon }_{t-2}+...+{\theta }_q\ast {\epsilon }_{t-q}$$

其中：

• $X_t$ 是时间 $t$ 的观测值
• $\mu$ 是均值
• ${\epsilon }_t$ 是时间 $t$ 的白噪声
• ${\theta }_i$ 是移动平均参数

6.2 MA模型的特征

$MA$模型具有以下重要特征：

• 自相关函数（$ACF$）在 $q$ 阶后截断
• 偏自相关函数（$PACF$）呈指数衰减
• 模型总是平稳的

# 模拟MA(2)过程
np.random.seed(42)
n = 1000
theta = [0.8, -0.3]  # MA参数# 生成白噪声
epsilon = np.random.normal(0, 1, n + 2)# 构造MA(2)序列
ma_series = np.zeros(n)
for t in range(n):ma_series[t] = epsilon[t+2] + theta[0]*epsilon[t+1] + theta[1]*epsilon[t]# 绘制ACF和PACF
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
sm.graphics.tsa.plot_acf(ma_series, lags=20, ax=ax1)
sm.graphics.tsa.plot_pacf(ma_series, lags=20, ax=ax2)
plt.tight_layout()
plt.show()

6.3 MA模型的应用

MA模型特别适合建模具有短期相关性的时间序列，如：

• 测量误差的建模
• 短期市场波动
• 随机冲击的影响

7. 数据平稳性检测

平稳性检验是时序分析的重要步骤，决定了后续使用的模型和方法。

7.1 ADF检验（Augmented Dickey-Fuller Test）

自相关图：自相关系数的快速衰减通常代表平稳性。

• 统计检测 ：
• 单位根检验：检验序列是否有单位根（非平稳特征）：ADF检验，PP检验
• KPSS检验：与ADF检验相反，假设序列是平稳的，检验序列中是否有趋势。

from statsmodels.tsa.stattools import adfuller# 对失业率数据进行ADF检验
result = adfuller(unemp['UNRATE'].dropna())
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:')
for key, value in result[4].items():print(f'\t{key}: {value}')# 解释结果
if result[1] <= 0.05:print("\n拒绝原假设：序列是平稳的")
else:print("\n不能拒绝原假设：序列是非平稳的")

7.2 PP检验（Phillips-Perron Test）

PP检验是ADF检验的改进版本，对序列相关和异方差更稳健：

from arch.unitroot import PhillipsPerron# PP检验
pp_test = PhillipsPerron(unemp['UNRATE'].dropna())
pp_result = pp_test.summary()
print(pp_result)

7.3 PP检验之Newey-West方法

Newey-West方法用于计算异方差和自相关一致的标准误：

# 使用Newey-West方法的PP检验
pp_test_nw = PhillipsPerron(unemp['UNRATE'].dropna(), lags=4)
print(pp_test_nw.summary())

7.4 KPSS检验

KPSS检验的原假设与ADF相反，原假设是序列平稳：

from statsmodels.tsa.stattools import kpss# KPSS检验
kpss_result = kpss(unemp['UNRATE'].dropna(), regression='c')
print('KPSS Statistic:', kpss_result[0])
print('p-value:', kpss_result[1])
print('Critical Values:')
for key, value in kpss_result[3].items():print(f'\t{key}: {value}')

7.5 综合判断

实践中常常综合多种检验结果：

def stationarity_test(timeseries, name="Series"):# ADF检验adf_result = adfuller(timeseries.dropna())# KPSS检验kpss_result = kpss(timeseries.dropna(), regression='c')print(f"=== {name} 平稳性检验结果 ===")print(f"ADF检验 p值: {adf_result[1]:.4f}")print(f"KPSS检验 p值: {kpss_result[1]:.4f}")if adf_result[1] < 0.05 and kpss_result[1] > 0.05:print("结论: 序列是平稳的")elif adf_result[1] > 0.05 and kpss_result[1] < 0.05:print("结论: 序列是非平稳的")else:print("结论: 检验结果不一致，需要进一步分析")

8. 归一化与标准化

数据归一化和标准化是预处理的重要步骤，特别是在使用机器学习方法时。

8.1 归一化（Normalization）

归一化将数据缩放到特定范围（通常是[0,1]）：

# Min-Max归一化
def min_max_normalize(series):return (series - series.min()) / (series.max() - series.min())unemp['UNRATE_normalized'] = min_max_normalize(unemp['UNRATE'])# 自定义范围归一化
def normalize_to_range(series, min_val=0, max_val=1):normalized = min_max_normalize(series)return normalized * (max_val - min_val) + min_val

8.2 标准化（Standardization）

标准化将数据转换为均值为0、标准差为1的分布：

# Z-score标准化
def z_score_standardize(series):return (series - series.mean()) / series.std()unemp['UNRATE_standardized'] = z_score_standardize(unemp['UNRATE'])

8.3 滑动窗口标准化

对于时序数据，使用滑动窗口进行标准化能够适应数据的动态变化：

# 滑动窗口标准化
window_size = 12
unemp['rolling_mean'] = unemp['UNRATE'].rolling(window=window_size).mean()
unemp['rolling_std'] = unemp['UNRATE'].rolling(window=window_size).std()
unemp['rolling_standardized'] = ((unemp['UNRATE'] - unemp['rolling_mean']) / unemp['rolling_std']
)# 绘制对比图
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(unemp['DATE'], unemp['UNRATE'], label='Original')
plt.plot(unemp['DATE'], unemp['rolling_mean'], label='Rolling Mean')
plt.fill_between(unemp['DATE'].dropna(), (unemp['rolling_mean'] - unemp['rolling_std']).dropna(),(unemp['rolling_mean'] + unemp['rolling_std']).dropna(), alpha=0.2, label='Rolling Std'
)
plt.legend()
plt.title('Original Data with Rolling Statistics')plt.subplot(2, 1, 2)
plt.plot(unemp['DATE'], unemp['rolling_standardized'], label='Rolling Standardized')
plt.axhline(0, color='gray', linestyle='--')
plt.axhline(1, color='gray', linestyle='--', alpha=0.5)
plt.axhline(-1, color='gray', linestyle='--', alpha=0.5)
plt.legend()
plt.title('Rolling Standardized Data')
plt.tight_layout()
plt.show()

9. 数据变换

数据变换用于稳定方差、使数据更接近正态分布，或者满足模型假设。

9.1 对数变换

对数变换适用于具有指数增长趋势或右偏分布的数据：

# 对航空乘客数据进行对数变换
passengers['log_value'] = np.log(passengers['value'])# 对比原始数据和对数变换后的数据
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))ax1.hist(passengers['value'], bins=30, alpha=0.7)
ax1.set_title('Original Distribution')ax2.hist(passengers['log_value'], bins=30, alpha=0.7, color='orange')
ax2.set_title('Log-Transformed Distribution')
plt.tight_layout()
plt.show()

9.2 平方根变换

平方根变换比对数变换温和，适用于轻度偏态的数据：

passengers['sqrt_value'] = np.sqrt(passengers['value'])# 对比原始数据和开跟变换后的数据
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))ax1.hist(passengers['value'], bins=30, alpha=0.7)
ax1.set_title('Original Distribution')ax2.hist(passengers['sqrt_value'], bins=30, alpha=0.7, color='orange')
ax2.set_title('Sqrt-Transformed Distribution')
plt.tight_layout()
plt.show()

9.3 Box-Cox变换

Box-Cox变换是一种更灵活的幂变换，能够自动选择最优的变换参数：

from scipy.stats import boxcox# Box-Cox变换
passengers['boxcox_value'], lambda_param = boxcox(passengers['value'])
print(f"最优λ参数: {lambda_param:.4f}")
# 最优λ参数: 0.1480
# 绘制变换效果
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.plot(passengers['date'], passengers['value'])
plt.title('Original Series')plt.subplot(2, 2, 2)
plt.plot(passengers['date'], passengers['log_value'])
plt.title('Log-Transformed')plt.subplot(2, 2, 3)
plt.plot(passengers['date'], passengers['sqrt_value'])
plt.title('Square Root Transformed')plt.subplot(2, 2, 4)
plt.plot(passengers['date'], passengers['boxcox_value'])
plt.title(f'Box-Cox Transformed (λ={lambda_param:.2f})')
plt.tight_layout()
plt.show()

"""
特征转换——对数变换和幂变换
"""
# 创建时间序列数据
dates = pd.date_range('2024-01-01', periods=100, freq='D')
data = {'value': np.random.exponential(scale=2, size=100)}  # 使用指数分布生成偏态数据
df = pd.DataFrame(data, index=dates)# 对数变换
df['log_value'] = np.log(df['value'] + 1)# 平方根变换
df['sqrt_value'] = np.sqrt(df['value'])# Box-Cox 变换
df['boxcox_value'], lambda_optimal = boxcox(df['value'] + 1)  # 注意: 数据必须为正值# 绘制原始数据和对数变换、平方根变换、Box-Cox 变换后的数据
plt.figure(figsize=(9, 9))
plt.subplot(2, 2, 1)
plt.hist(df['value'], bins=30, color='blue', alpha=0.5)
plt.title('Original Data Distribution')
plt.subplot(2, 2, 2)
plt.hist(df['log_value'], bins=30, color='orange', alpha=0.5)
plt.title('Log-Transformed Data Distribution')
plt.subplot(2, 2, 3)
plt.hist(df['sqrt_value'], bins=30, color='green', alpha=0.5)
plt.title('Square Root Transformed Data Distribution')
plt.subplot(2, 2, 4)
plt.hist(df['boxcox_value'], bins=30, color='purple', alpha=0.5)
plt.title(f'Box-Cox Transformed Data Distribution (λ={lambda_optimal:.2f})')plt.tight_layout()
plt.show()

9.4 差分变换

差分变换用于去除趋势，使非平稳序列变平稳：

unemp = pd.read_csv("data/unemployment_rate.csv")
# 一阶差分
unemp['diff_1'] = unemp['UNRATE'].diff(1)
# 二阶差分
unemp['diff_2'] = unemp['diff_1'].diff(1)
# 季节差分（12个月）
unemp['seasonal_diff'] = unemp['UNRATE'].diff(12)
# 绘制差分结果
fig, axes = plt.subplots(4, 1, figsize=(12, 10))
unemp['UNRATE'].plot(ax=axes[0], title='Original Series')
unemp['diff_1'].plot(ax=axes[1], title='First Difference')
unemp['diff_2'].plot(ax=axes[2], title='Second Difference')
unemp['seasonal_diff'].plot(ax=axes[3], title='Seasonal Difference (12)')
plt.tight_layout()
plt.show()

10. 数据降维

当处理多变量时序数据时，降维技术能够提取主要特征，简化分析。

10.1 主成分分析（PCA）

PCA通过线性变换将数据投影到低维空间：

from sklearn.decomposition import PCA# 使用欧洲股市数据进行PCA
stock = pd.read_csv("./data/EuStockMarkets.csv")
stock_data = stock[['DAX', 'SMI', 'CAC', 'FTSE']].values# 标准化数据
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
stock_scaled = scaler.fit_transform(stock_data)# 执行PCA
pca = PCA(n_components=2)
principal_components = pca.fit_transform(stock_scaled)# 查看解释方差比例
print("解释方差比例:", pca.explained_variance_ratio_)
print("累计解释方差:", pca.explained_variance_ratio_.cumsum())# 解释方差比例: [0.97324812 0.02302963]
# 累计解释方差: [0.97324812 0.99627775]# 绘制主成分
plt.figure(figsize=(10, 6))
plt.scatter(principal_components[:, 0], principal_components[:, 1], alpha=0.5)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
plt.title('PCA of European Stock Markets')
plt.show()

10.2 因子分析（Factor Analysis）

因子分析假设观测变量由少数潜在因子驱动：

from sklearn.decomposition import FactorAnalysis# 执行因子分析
fa = FactorAnalysis(n_components=2)
factors = fa.fit_transform(stock_scaled)# 查看因子载荷
loadings = pd.DataFrame(fa.components_.T,columns=['Factor 1', 'Factor 2'],index=['DAX', 'SMI', 'CAC', 'FTSE']
)
print("因子载荷矩阵:")
print(loadings)# 绘制因子得分
plt.figure(figsize=(10, 6))
plt.scatter(factors[:, 0], factors[:, 1], alpha=0.5, color='orange')
plt.xlabel('Factor 1')
plt.ylabel('Factor 2')
plt.title('Factor Analysis of European Stock Markets')
plt.show()# 因子载荷矩阵:
#      Factor 1  Factor 2
# DAX   0.994394  0.075102
# SMI   0.998292 -0.020644
# CAC   0.953440  0.241479
# FTSE  0.989301 -0.114054

10.3 时序特征提取

除了传统降维方法，还可以提取时序特有的特征：

def extract_time_features(series, window=12):"""提取时序统计特征"""features = pd.DataFrame(index=series.index)# 滚动统计量features['rolling_mean'] = series.rolling(window).mean()features['rolling_std'] = series.rolling(window).std()features['rolling_max'] = series.rolling(window).max()features['rolling_min'] = series.rolling(window).min()# 变化率features['pct_change'] = series.pct_change()features['diff'] = series.diff()# 技术指标features['ema'] = series.ewm(span=window).mean()return features# 应用特征提取
unemp_features = extract_time_features(unemp['UNRATE'])

11. 数据可视化

有效的可视化是理解时序数据的关键。

11.1折线图、直方图、散点图。

"""
绘制常用的三种可视化图表——折线图、直方图、散点图
"""
# %%  读取数据
stock = pd.read_csv("./data/EuStockMarkets.csv")# 查看数据前几行
stock.head()# 汇总统计
stock.describe()# %%  绘制折线图plt.figure(figsize=(10, 6))
plt.plot(stock["rownames"], stock["DAX"], label="DAX")
plt.plot(stock["rownames"], stock["SMI"], label="SMI")
plt.plot(stock["rownames"], stock["CAC"], label="CAC")
plt.plot(stock["rownames"], stock["FTSE"], label="FTSE")# 添加图表标题和标签
plt.title("European Stock Markets Over Time")
plt.xlabel("Time (days)")
plt.ylabel("Index Value")
plt.legend()
plt.grid(True)
plt.show()

# %%  绘制直方图
![在这里插入图片描述](https://i-blog.csdnimg.cn/direct/717f9499527548c5b87ac03a3d3134d6.png)
# 绘制4个子图的直方图
fig, axs = plt.subplots(2, 2, figsize=(12, 10))# DAX 子图
axs[0, 0].hist(stock["DAX"], bins=10, alpha=0.35, color='blue')
axs[0, 0].set_title('DAX')
axs[0, 0].set_xlabel('Index Value')
axs[0, 0].set_ylabel('Frequency')# SMI 子图
axs[0, 1].hist(stock['SMI'], bins=10, alpha=0.35, color='green')
axs[0, 1].set_title('SMI')
axs[0, 1].set_xlabel('Index Value')
axs[0, 1].set_ylabel('Frequency')# CAC 子图
axs[1, 0].hist(stock['CAC'], bins=10, alpha=0.35, color='red')
axs[1, 0].set_title('CAC')
axs[1, 0].set_xlabel('Index Value')
axs[1, 0].set_ylabel('Frequency')# FTSE 子图
axs[1, 1].hist(stock['FTSE'], bins=10, alpha=0.35, color='purple')
axs[1, 1].set_title('FTSE')
axs[1, 1].set_xlabel('Index Value')
axs[1, 1].set_ylabel('Frequency')# 调整布局
plt.tight_layout()
plt.show()

# %%  绘制差分处理直方图
# 绘制2个子图的直方图
fig, axs = plt.subplots(2, 1, figsize=(6, 5))
stock_smi_diff = stock['SMI'].diff()
# SMI 原始直方图
axs[0].hist(stock['SMI'], bins=10, alpha=0.35, color='green')
axs[0].set_title('SMI')
axs[0].set_xlabel('Index Value')
axs[0].set_ylabel('Frequency')# SMI 差分直方图
axs[1].hist(stock_smi_diff.dropna(), bins=10, alpha=0.35, color='green')
axs[1].set_title('SMI of First Difference')
axs[1].set_xlabel('Index Value')
axs[1].set_ylabel('Frequency of First Difference')# 调整布局
plt.tight_layout()
plt.show()

# %%  绘制散点图、差分散点图
stock_diff = stock.diff()# 绘制2个子图
fig, axs = plt.subplots(2, 1, figsize=(6, 5))axs[0].scatter(stock['DAX'], stock['SMI'], marker='o', facecolors='none', edgecolors='#9ecb9f')
axs[0].set_title('Scatter Plot of DAX vs SMI')
axs[0].set_xlabel('DAX Index Value')
axs[0].set_ylabel('SMI Index Value')axs[1].scatter(stock_diff['DAX'], stock_diff['SMI'], marker='o', facecolors='none', edgecolors='#9ecb9f')
axs[1].set_title('First Difference Scatter Plot of DAX vs SMI')
axs[1].set_xlabel('DAX Index Value of First Difference')
axs[1].set_ylabel('SMI Index Value of First Difference')plt.tight_layout()
plt.show()

# %%  绘制差分滞后1散点图plt.figure(figsize=(8, 6))
plt.scatter(stock_diff['SMI'].shift(1), stock_diff['DAX'], marker='o', facecolors='none', edgecolors='#9ecb9f')# 添加标题和标签
plt.title('Scatter Plot of Lagged SMI Difference vs DAX Difference')
plt.xlabel('Lagged SMI Difference')
plt.ylabel('DAX Difference')# 显示图表
plt.show()

11.2 时序图

时序图是最基本的可视化方法：

# 多变量时序图
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)for i, col in enumerate(['DAX', 'SMI', 'CAC', 'FTSE']):axes[i].plot(stock.index, stock[col])axes[i].set_ylabel(col)axes[i].grid(True, alpha=0.3)axes[-1].set_xlabel('Time Index')
plt.suptitle('European Stock Markets Time Series')
plt.tight_layout()
plt.show()

11.2 ACF和PACF图

自相关函数（ACF）和偏自相关函数（PACF）图用于识别时序模式：

# 使用中国银行数据
chinabank = pd.read_csv('./data/ChinaBank.csv', index_col='Date', parse_dates=['Date'])
train_data = chinabank['Close']['2014-01':'2014-03']fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))# ACF图
sm.graphics.tsa.plot_acf(train_data, lags=20, ax=ax1)
ax1.set_title('Autocorrelation Function')# PACF图
sm.graphics.tsa.plot_pacf(train_data, lags=20, ax=ax2)
ax2.set_title('Partial Autocorrelation Function')plt.tight_layout()
plt.show()

11.3 季节性分解图

from statsmodels.tsa.seasonal import seasonal_decompose# 对航空乘客数据进行季节性分解
passengers.set_index('date', inplace=True)
decomposition = seasonal_decompose(passengers['value'], model='multiplicative', period=12)# 绘制分解结果
fig, axes = plt.subplots(4, 1, figsize=(12, 10))
passengers['value'].plot(ax=axes[0], title='Original Series')
decomposition.trend.plot(ax=axes[1], title='Trend Component')
decomposition.seasonal.plot(ax=axes[2], title='Seasonal Component')
decomposition.resid.plot(ax=axes[3], title='Residual Component')
plt.tight_layout()
plt.show()

11.4 热力图

热力图适合展示时序数据的周期性模式：

# 将航空乘客数据重塑为年-月矩阵
passengers_matrix = passengers['value'].values.reshape(12, 12)
months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
years = [str(year) for year in range(1949, 1961)]# 绘制热力图
plt.figure(figsize=(10, 8))
sns.heatmap(passengers_matrix.T, xticklabels=months, yticklabels=years,cmap='YlOrRd',annot=True,fmt='d')
plt.title('Air Passengers Heatmap')
plt.xlabel('Month')
plt.ylabel('Year')
plt.show()

11.5 三维可视化

对于多维时序数据，三维可视化能够展示更多信息：

"""
三维可视化
"""
# %%  加载 AirPassengers 数据集
air_passengers = pd.read_csv("data/AirPassengers.csv")['value'].values# 将数据重塑为 12x12 的矩阵，并转换为 DataFrame
months = np.arange(1, 13)
ap = pd.DataFrame(air_passengers.reshape(12, 12), columns=[str(year) for year in range(1949, 1961)])# 添加 'month' 列
ap['month'] = months# 使用 melt 函数将数据转换为长格式
ap_long = pd.melt(ap, id_vars='month', var_name='year', value_name='count')# 将 year 列转换为数值型
ap_long['year'] = ap_long['year'].astype(int)# %%  绘制三维图fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111, projection='3d')# 绘制 3D 散点图
scatter = ax.scatter(ap_long['month'],ap_long['year'],ap_long['count'],c=ap_long['month'],cmap='viridis',s=50,alpha=0.8
)ax.set_xlabel('Month')
ax.set_ylabel('Year')
ax.set_zlabel('Passenger Count')
ax.set_title('AirPassengers 3D Visualization')
ax.yaxis.label.set_rotation(45)# 添加颜色条
cbar = fig.colorbar(scatter, ax=ax, pad=0.1)
cbar.set_label('Month')
cbar.ax.yaxis.label.set_rotation(0)plt.tight_layout()plt.show()# %%  加载 二维随机游走过程 数据集file_location = 'data/_3d-line-plot.csv'
data = pd.read_csv(file_location)# %%  绘制 二维随机游走过程 三维图fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111, projection='3d')# 绘制 3D 线图
ax.plot(data['x1'], data['y1'], data['z1'], color='#1f77b4', linewidth=1)# 设置轴标签
ax.set_xlabel('X Axis')
ax.set_ylabel('Y Axis')
ax.set_zlabel('Z Axis')
ax.set_title('Time-Varying Two-Dimensional Random Walk Process')plt.tight_layout()plt.show()

12. 数据分析与应用

12.1 自回归模型（AR）

AR模型假设当前值是过去值的线性组合：

# 使用statsmodels拟合AR模型
from statsmodels.tsa.ar_model import AutoReg# 拟合AR(2)模型
ar_model = AutoReg(train_data, lags=2)
ar_result = ar_model.fit()
print(ar_result.summary())# 预测
ar_forecast = ar_result.predict(start=len(train_data), end=len(train_data)+10)

12.2 移动平均模型（MA）

MA模型使用过去的预测误差：

# MA(q)模型已在前面详细介绍
# 这里展示如何在ARIMA框架下使用
from statsmodels.tsa.arima.model import ARIMA# 拟合MA(2)模型
ma_model = ARIMA(train_data, order=(0, 0, 2))
ma_result = ma_model.fit()
print(ma_result.summary())

12.3ARIMA模型详解

1. 什么是ARIMA模型？

ARIMA是"自回归综合移动平均"（AutoRegressive Integrated Moving Average）模型的缩写。它是一种专门用来分析和预测时间序列数据的统计方法。

简单理解： 想象你要预测明天的股价、下个月的销售额，或者未来几天的气温。ARIMA模型就是通过分析历史数据的规律，来帮你做出这些预测的工具。

2. ARIMA的三个组成部分

ARIMA模型由三个字母组成，每个字母代表一个重要概念：

A - 自回归（AutoRegressive，AR）

通俗解释： “今天的我受昨天的我影响”

• 自回归意思是用过去的值来预测现在的值
• 比如：今天的股价可能与昨天、前天的股价有关
• 数学表达：今天的值 = 昨天的值 × 系数1 + 前天的值 × 系数2 + …

生活例子： 你今天的心情可能受昨天心情的影响。如果昨天很开心，今天也更可能开心。

R - 差分（Integrated，I）

通俗解释： “让不稳定的数据变稳定”

• 有些数据有明显的上升或下降趋势，这叫"不平稳"
• 差分就是用"今天的值 - 昨天的值"来消除趋势
• 目的是让数据围绕一个固定水平波动

生活例子： 房价总是在上涨，但房价的涨幅可能是相对稳定的。差分就是关注涨幅而不是房价本身。

M - 移动平均（Moving Average，MA）

通俗解释： “今天的意外影响明天”

• 移动平均考虑的是随机冲击（意外事件）的影响
• 比如：昨天的一个突发新闻可能影响今天的股价
• 这些"意外"的影响会延续一段时间

生活例子： 你昨天看了一部很棒的电影，这个"意外的惊喜"可能让你今天也心情很好。

3. ARIMA模型的表示方法

ARIMA模型通常写作：ARIMA(p, d, q)

• p：自回归项的阶数（AR）- 用几天前的数据来预测
• d：差分次数（I）- 需要做几次差分才能让数据稳定
• q：移动平均项的阶数（MA）- 考虑几天前的随机冲击

常见例子解释

ARIMA(1,1,1)：

• p=1：用1天前的数据预测今天
• d=1：对原数据做1次差分
• q=1：考虑1天前的随机冲击影响

ARIMA(2,0,1)：

• p=2：用前2天的数据预测今天
• d=0：数据已经稳定，不需要差分
• q=1：考虑1天前的随机冲击影响

ARIMA模型是时间序列分析的经典工具，它通过分析数据的自相关性、处理数据的不平稳性、考虑随机冲击的影响，来实现对未来值的预测。虽然有一定的学习门槛，但掌握后在商业分析、经济预测等领域都有广泛应用。记住，任何模型都不是万能的，关键是要理解模型的假设和适用条件，结合实际情况灵活运用。

12.5 模型诊断

模型诊断是确保模型有效性的重要步骤：

# 残差分析
residuals = final_result.residfig, axes = plt.subplots(2, 2, figsize=(12, 8))# 残差时序图
axes[0, 0].plot(residuals)
axes[0, 0].set_title('Residuals')# 残差直方图
axes[0, 1].hist(residuals, bins=20, density=True, alpha=0.7)
axes[0, 1].set_title('Residual Distribution')# 残差ACF
sm.graphics.tsa.plot_acf(residuals, lags=20, ax=axes[1, 0])
axes[1, 0].set_title('Residual ACF')# Q-Q图
sm.qqplot(residuals, line='s', ax=axes[1, 1])
axes[1, 1].set_title('Normal Q-Q Plot')plt.tight_layout()
plt.show()# Ljung-Box检验
from statsmodels.stats.diagnostic import acorr_ljungbox
lb_test = acorr_ljungbox(residuals, lags=10, return_df=True)
print("Ljung-Box检验结果:")
print(lb_test)

---

总结

时序数据预处理是一个系统性的过程，涉及多个步骤和技术。本文详细介绍了从基础概念到高级应用的完整流程：

1. 理解数据特性：掌握时序数据的基本概念和统计特征
2. 处理缺失值：选择合适的方法填补或处理缺失数据
3. 数据平滑：减少噪声，识别潜在模式
4. 平稳性检验：确保数据满足模型假设
5. 数据变换：通过各种变换技术改善数据特性
6. 特征工程：提取和构造有意义的特征
7. 可视化分析：通过图表深入理解数据
8. 模型应用：选择和应用合适的时序模型

时序数据分析是一个不断发展的领域，新的方法和技术不断涌现。建议在掌握基础知识后，继续关注最新的研究进展，如深度学习在时序预测中的应用、高维时序数据的处理方法等。

没有一种方法适用于所有情况。成功的关键在于深入理解数据特性，选择合适的工具，并通过实践不断积累经验。