今天云朵君和大家一起学习如何使用探索性数据分析从时间序列数据中获取信息，并使用 Python 加强特征工程设计。
简介
时间序列预测是数据科学和机器学习领域中极其重要的应用场景，广泛运用于金融、能源、零售等众多行业，对于企业来说具有重大价值。随着数据获取能力的提升和机器学习模型的不断进化，时间序列预测技术也日趋丰富和成熟。
传统的统计预测方法，如回归模型、ARIMA模型和指数平滑等，一直是该领域的基础。近年来，机器学习算法如基于树的模型，以及深度学习技术如LSTM网络、卷积神经网络和基于Transformer的模型，也逐步应用于时间序列预测，都取得了不错的成绩。
尽管上述各种模型和技术存在显著差异，但无论采用何种方法，探索性数据分析(Exploratory Data Analysis，EDA)都是时间序列预测不可或缺的第一步。探索性数据分析是一门数据分析和可视化技巧，旨在总结数据的主要统计特征并从中提取有价值的信息。在数据科学中，EDA为后续的特征工程奠定了基础，有助于从原始数据集中创建、转换和提取最有效的特征，从而最大限度地发挥机器学习模型的潜力。
本文算是定义了一个针对时间序列数据的探索性数据分析模板，全面总结和突出时间序列数据集的关键特征。这里我们将使用流行的Python数据分析库，如Pandas、Seaborn和Statsmodels等，来实现这一目标。
数据
在本文中，我们将使用 Kaggle 的 数据。该数据集与 PJM 每小时能耗数据有关，PJM 是美国的一个区域输电组织，为Delaware、Illinois、Indiana、Kentucky、Maryland、Michigan、New Jersey、North Carolina、Ohio、Pennsylvania、Tennessee、Virginia、West Virginia 和 the District of Columbia（特拉华州、伊利诺伊州、印第安纳州、肯塔基州、马里兰州、密歇根州、新泽西州、北卡罗来纳州、俄亥俄州、宾夕法尼亚州、田纳西州、弗吉尼亚州、西弗吉尼亚州和哥伦比亚特区）提供电力服务。
每小时电力消耗数据来自 PJM 网站，单位为兆瓦 (MW)。
探索性数据分析
时间序列分析的关键步骤包括绘制数据图，利用图表突出特征、模式、不寻常的观察结果，以及变量之间的关系。这些图表的见解必须纳入预测模型中，同时还可以利用描述性统计和时间序列分解等数学工具来提高分析效果。
因此，我在本文中提出的 EDA 包括六个步骤：描述性统计、时间图、季节图、箱形图、时间序列分解、滞后分析。
1. 描述性统计
描述性统计是一种用于定量描述或总结结构化数据集合特征的汇总统计方法。常用的指标包括中心倾向度量（如平均值、中位数）、离散度量（如范围、标准偏差）和位置度量（如百分位数、四分位数）。这些指标包括分布的最小值、第一四分位数 (Q1)、中位数或第二四分位数 (Q2)、第三四分位数 (Q3) 和最大值。
在 Python 中，可以使用 Pandas 中广为人知的 describe 方法轻松获取这些信息：
import pandas as pd
# Loading and preprocessing steps
df = pd.read_csv('hourly-energy-consumption/PJME_hourly.csv')
df = df.set_index('Datetime')
df.index = pd.to_datetime(df.index)
df.describe()
1.PJME 统计摘要。
2. 时间图
首先要绘制的图形显然是时间图。也就是说，将观测值与观测时间相对应，用线条连接连续的观测值。
在 Python 中，我们可以使用 Pandas 和 Matplotlib：
import matplotlib.pyplot as plt
# Set pyplot style
plt.style.use("seaborn")
# Plot
df['PJME_MW'].plot(title='PJME - Time Plot', figsize=(10,6))
plt.ylabel('Consumption [MW]')
plt.xlabel('Date')
图片
2.1 PJME 消耗时间图
该曲线图提供了一些信息：
- 年度季节性模式，体现在冬夏两季耗电量高于其他季节。
- 多年来，整体耗电量未显现出明显的上升或下降趋势，平均消耗量保持稳定水平。
- 2023年前后存在一个异常值，在建模时需予以考虑。
- 除此之外，单个年份内还可能存在其他影响耗电量的因素。
3. 季节图
季节图从根本上说是一种时间图，其中的数据是根据其所属系列的各个 "季节" 绘制的。
在能源消耗方面，我们通常有每小时的数据，因此可能会有几种季节性： 年、周、日。在深入研究这些图表之前，先在 Pandas 数据框中设置一些变量：
# Defining required fields
df['year'] = [x for x in df.index.year]
df['month'] = [x for x in df.index.month]
df = df.reset_index()
df['week'] = df['Datetime'].apply(lambda x:x.week)
df = df.set_index('Datetime')
df['hour'] = [x for x in df.index.hour]
df['day'] = [x for x in df.index.day_of_week]
df['day_str'] = [x.strftime('%a') for x in df.index]
df['year_month'] = [str(x.year) + '_' + str(x.month) for x in df.index]
图片
3.1 季节性曲线图--年度消耗量
这个图表按照年份和月份对能源消耗进行了分组，展现了每年的季节性变化，并展示了多年来的上升或下降趋势。
下面是 Python 代码：
import numpy as np
# Defining colors palette
np.random.seed(42)
df_plot = df[['month', 'year', 'PJME_MW']].dropna().groupby(['month', 'year']).mean()[['PJME_MW']].reset_index()
years = df_plot['year'].unique()
colors = np.random.choice(list(mpl.colors.XKCD_COLORS.keys()), len(years), replace=False)
# Plot
plt.figure(figsize=(16,12))
for i, y in enumerate(years):
if i > 0:
plt.plot('month', 'PJME_MW', data=df_plot[df_plot['year'] == y], color=colors[i], label=y)
if y == 2018:
plt.text(df_plot.loc[df_plot.year==y, :].shape[0]+0.3, df_plot.loc[df_plot.year==y, 'PJME_MW'][-1:].values[0], y, fnotallow=12, color=colors[i])
else:
plt.text(df_plot.loc[df_plot.year==y, :].shape[0]+0.1, df_plot.loc[df_plot.year==y, 'PJME_MW'][-1:].values[0], y, fnotallow=12, color=colors[i])
# Setting labels
plt.gca().set(ylabel= 'PJME_MW', xlabel = 'Month')
plt.yticks(fnotallow=12, alpha=.7)
plt.title("Seasonal Plot - Monthly Consumption", fnotallow=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Month')
plt.show()
图片
3.1 PJME 年度季节图
可以观察到，每年都有一个非常明显的模式：冬季耗电量急剧增加，夏季耗电量达到峰值，而春秋两季通常不需要供暖或制冷，因此耗电量最小。
另外，图表也表明，不同年份的总体消耗量并没有明显的增减趋势。
3.2 季节图--每周消耗量
周曲线图是一种有用的曲线图类型，它展示了每周消耗量的变化情况，并能够揭示一年中每周的消耗量变化趋势。
用 Python 来计算：
# Defining colors palette
np.random.seed(42)
df_plot = df[['month', 'day_str', 'PJME_MW', 'day']].dropna().groupby(['day_str', 'month', 'day']).mean()[['PJME_MW']].reset_index()
df_plot = df_plot.sort_values(by='day', ascending=True)
months = df_plot['month'].unique()
colors = np.random.choice(list(mpl.colors.XKCD_COLORS.keys()), len(months), replace=False)
# Plot
plt.figure(figsize=(16,12))
for i, y in enumerate(months):
if i > 0:
plt.plot('day_str', 'PJME_MW', data=df_plot[df_plot['month'] == y], color=colors[i], label=y)
if y == 2018:
plt.text(df_plot.loc[df_plot.mnotallow==y, :].shape[0]-.9, df_plot.loc[df_plot.mnotallow==y, 'PJME_MW'][-1:].values[0], y, fnotallow=12, color=colors[i])
else:
plt.text(df_plot.loc[df_plot.mnotallow==y, :].shape[0]-.9, df_plot.loc[df_plot.mnotallow==y, 'PJME_MW'][-1:].values[0], y, fnotallow=12, color=colors[i])
# Setting Labels
plt.gca().set(ylabel= 'PJME_MW', xlabel = 'Month')
plt.yticks(fnotallow=12, alpha=.7)
plt.title("Seasonal Plot - Weekly Consumption", fnotallow=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Month')
plt.show()
3.2 PJME 每周季节图
3.3 季节性曲线图--日消耗量
最后一个季节性曲线图要展示的是日消耗量曲线图。如您所猜测的那样，它显示了一天中消耗量的变化。数据被按星期分组并取平均值进行汇总。
下面是代码：
import seaborn as sns
# Defining the dataframe
df_plot = df[['hour', 'day_str', 'PJME_MW']].dropna().groupby(['hour', 'day_str']).mean()[['PJME_MW']].reset_index()
# Plot using Seaborn
plt.figure(figsize=(10,8))
sns.lineplot(data = df_plot, x='hour', y='PJME_MW', hue='day_str', legend=True)
plt.locator_params(axis='x', nbins=24)
plt.title("Seasonal Plot - Daily Consumption", fnotallow=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Hour')
plt.legend()
图片
3.3 PJME 日季节图
通常情况下，这个图表展示了一种常见的模式，被称为"M型曲线"，因为它似乎在一天中描绘出了一个"M"的形状。有时这种模式非常明显，而有时则不那么明显（就像在这个例子中）。
然而，这种曲线通常表现为一天中间（上午10点到下午2点）出现一个相对高峰，接着是另一个高峰（下午2点到下午6点）和最后一个高峰（下午6点到晚上8点）。最后，它还展示了周末和其他日子的用电量差异。
3.4 季节图--特征工程
探讨如何将这些信息应用于特征工程。假设我们正在使用一些需要高质量特征的 ML 模型（如 ARIMA 模型或基于树的模型）。主要的证据来自季节图包括以下几点：
- 年度消耗量在不同年份之间的变化不大，这意味着可以利用年度季节性特征，例如滞后变量或外生变量。
- 周消费量在各月份中的变化规律相似，这表明可以利用周特征，如滞后变量或外生变量。
- 日常消费与平日和周末有所不同，因此应当使用分类特征来区分平日和非平日。
4. 箱形图
箱形图是一种有效的方法来确定数据分布情况。简而言之，它描述了百分位数，包括第一四分位数（Q1）、第二四分位数（Q2/中位数）和第三四分位数（Q3），以及箱图代表的数据范围。超出箱图的每一个值都可以被视为离群值。更详细地说，箱图通常是通过以下方式计算的：
箱图公式
4.1 箱形图 - 总消耗量
我们首先来计算总消耗量的箱形图，这可以通过 Seaborn 轻松完成：
plt.figure(figsize=(8,5))
sns.boxplot(data=df, x='PJME_MW')
plt.xlabel('Consumption [MW]')
plt.title(f'Boxplot - Consumption Distribution');
图片
4.1 PJME 方框图
尽管这幅图看起来信息量不大，但它告诉我们，我们面对的是一个类似高斯分布的 分布，其尾部更偏向右侧。
4.2 箱形图--日月分布
箱形图非常有趣，它利用 "日-月" 变量对消耗量进行分组来展现数据。以下是仅涉及 2017 年的相关代码：
df['year'] = [x for x in df.index.year]
df['month'] = [x for x in df.index.month]
df['year_month'] = [str(x.year) + '_' + str(x.month) for x in df.index]
df_plot = df[df['year'] >= 2017].reset_index().sort_values(by='Datetime').set_index('Datetime')
plt.title(f'Boxplot Year Month Distribution');
plt.xticks(rotatinotallow=90)
plt.xlabel('Year Month')
plt.ylabel('MW')
sns.boxplot(x='year_month', y='PJME_MW', data=df_plot)
plt.ylabel('Consumption [MW]')
plt.xlabel('Year Month')
4.2 PJME 年/月箱形图
可以观察到，在夏季和冬季月份（即出现高峰时）的消费量的不确定性较小，而在春季和秋季月份（即气温变化较大时）的消费量较为分散。值得注意的是，2018年夏季的消费量高于2017年，这可能是由于夏季较为温暖的原因。在进行特征工程设计时，请务必考虑将温度曲线（如果有的话）纳入考虑范围，或许它可以作为外生变量。
4.3 箱形图--日分布
另一种有用的曲线图是一周内的消耗量分布图，这与每周消耗量季节曲线图类似。
df_plot = df[['day_str', 'day', 'PJME_MW']].sort_values(by='day')
plt.title(f'Boxplot Day Distribution')
plt.xlabel('Day of week')
plt.ylabel('MW')
sns.boxplot(x='day_str', y='PJME_MW', data=df_plot)
plt.ylabel('Consumption [MW]')
plt.xlabel('Day of week')
4.3 PJME 日箱形图
如前所述，周末的消耗量明显较低。无论如何，有几个异常值表明，"星期" 等日历特征肯定是有用的，但不能完全解释这一系列数据。
4.4 箱形图--小时分布
最后让我们来看看小时分布箱形图。它与每日消费季节图相似，因为它提供了消费在一天中的分布情况。代码如下
plt.title(f'Boxplot Hour Distribution');
plt.xlabel('Hour')
plt.ylabel('MW')
sns.boxplot(x='hour', y='PJME_MW', data=df)
plt.ylabel('Consumption [MW]')
plt.xlabel('Hour')
4.4 PJME 小时箱形图
请注意，之前看到的 "M" 形现在变得更粗了。此外，还有很多离群值，这说明数据不仅依赖于每日的季节性（例如，今天上午 12 点的消耗量与昨天上午 12 点的消耗量相似），还依赖于其他因素，可能是温度或湿度等外生气候特征。
5. 时间序列分解
如之前所述，时间序列数据能够展示出多种模式。通常情况下，将时间序列分解成几个部分是非常有帮助的，每个部分代表一个基本模式类别。
时间序列可以被分解成三个部分：趋势部分、季节部分和残差部分（包含时间序列中的任何其他成分）。对于一些时间序列（例如能源消耗序列），可能会存在不止一个季节性成分，分别对应不同的季节性周期（日、周、月、年）。
分解的主要类型有两种：加法和乘法。
对于加法分解，我们将一个序列（𝑦）表示为季节成分（𝑆）、趋势（𝑇）和余数（𝑅）的总和：
同样，乘法分解可以写成
一般来说，加法分解最能代表方差恒定的序列，而乘法分解最适合方差非平稳的时间序列。
在 Python 中，Statsmodel 库可以轻松实现时间序列分解：
df_plot = df[df['year'] == 2017].reset_index()
df_plot = df_plot.drop_duplicates(subset=['Datetime']).sort_values(by='Datetime')
df_plot = df_plot.set_index('Datetime')
df_plot['PJME_MW - Multiplicative Decompose'] = df_plot['PJME_MW']
df_plot['PJME_MW - Additive Decompose'] = df_plot['PJME_MW']
# Additive Decomposition
result_add = seasonal_decompose(df_plot['PJME_MW - Additive Decompose'], model='additive', period=24*7)
# Multiplicative Decomposition
result_mul = seasonal_decompose(df_plot['PJME_MW - Multiplicative Decompose'], model='multiplicative', period=24*7)
# Plot
result_add.plot().suptitle('', fnotallow=22)
plt.xticks(rotatinotallow=45)
result_mul.plot().suptitle('', fnotallow=22)
plt.xticks(rotatinotallow=45)
plt.show()
5.1 PJME 系列分解--加法分解
5.2 PJME 系列分解--乘法分解
2017年的数据显示了两种情况下的趋势，发现有几个局部峰值，其中夏季值偏高。从季节性成分来看，该序列实际上呈现出几个周期性，图表更加突出了周周期性，但如果我们关注同一年的特定月份（例如1月），也会出现日周期性。
df_plot = df[(df['year'] == 2017)].reset_index()
df_plot = df_plot[df_plot['month'] == 1]
df_plot['PJME_MW - Multiplicative Decompose'] = df_plot['PJME_MW']
df_plot['PJME_MW - Additive Decompose'] = df_plot['PJME_MW']
df_plot = df_plot.drop_duplicates(subset=['Datetime']).sort_values(by='Datetime')
df_plot = df_plot.set_index('Datetime')
# Additive Decomposition
result_add = seasonal_decompose(df_plot['PJME_MW - Additive Decompose'], model='additive', period=24*7)
# Multiplicative Decomposition
result_mul = seasonal_decompose(df_plot['PJME_MW - Multiplicative Decompose'], model='multiplicative', period=24*7)
# Plot
result_add.plot().suptitle('', fnotallow=22)
plt.xticks(rotatinotallow=45)
result_mul.plot().suptitle('', fnotallow=22)
plt.xticks(rotatinotallow=45)
plt.show()
5.3 PJME 系列分解--加法分解，聚焦 2017 年 1 月
5.4 PJME 系列分解--乘法分解，聚焦 2017 年 1 月
6. 滞后分析
在时间序列预测中，滞后期就是序列的过去值。例如，对于日序列，第一个滞后期指的是序列前一天的值，第二个滞后期指的是前一天的值，以此类推。
滞后分析的基础是计算序列与序列本身的滞后版本之间的相关性，这也称为*自相关：
图片
其中y条代表序列的平均值，k代表滞后期。
自相关系数构成了序列的自相关函数（ACF），展现了自相关系数与所考虑的滞后期数的关系的曲线图。
当数据具有趋势性时，较小滞后期的自相关系数通常较大且为正，因为时间上接近的观测值在数值上也接近。当数据具有季节性时，与季节性滞后期（和季节性周期的倍数）相对应的自相关值会比其他滞后期大。同时，具有趋势和季节性的数据将显示这些效应的组合。
实际上，更有用的函数是部分自相关函数（PACF）。这与ACF相似，只是它仅显示两个滞后期之间的直接自相关。例如，滞后3的部分自相关指的是滞后1和滞后2无法解释的唯一相关性。换句话说，部分相关指的是某个滞后期对当前时间值的直接影响。
在开始Python代码之前，需要强调的是，如果序列是稳定的，自相关系数会更加明显。因此，最好先将序列区分开来，以识别稳定信号。
下面是绘制一天中不同时段 PACF 的代码：
from statsmodels.graphics.tsaplots import plot_pacf
actual = df['PJME_MW']
hours = range(0, 24, 4)
for hour in hours:
plot_pacf(actual[actual.index.hour == hour].diff().dropna(), lags=30, alpha=0.01)
plt.title(f'PACF - h = {hour}')
plt.ylabel('Correlation')
plt.xlabel('Lags')
plt.show()
6.1 PJME 滞后分析 - 部分自相关函数（h=0）
6.2 PJME 滞后分析 - 部分自相关函数（h=4）
6.3 PJME 滞后分析 - 部分自相关函数（h=8）
6.4 PJME 滞后分析 - 部分自相关函数（h=12）
6.5 PJME 滞后分析 - 部分自相关函数（h=16）
6.6 PJME 滞后分析 - 部分自相关函数（h=20）
PACF 只包括绘制不同滞后期的皮尔逊部分自相关系数。非滞后序列本身显示出完美的自相关性，因此滞后 0 始终为 1。蓝色区域代表置信区间： 统计上具有显著性的滞后期将超过该区域，从而可以断言其具有重要意义。
6.1 滞后分析--特征工程
滞后分析是对时间序列特征工程影响最大的研究之一。如前所述，相关性高的滞后期是序列的重要滞后期，因此应加以考虑。
广泛使用的特征工程技术包括对数据集进行小时分割。也就是说，将数据分成 24 个子集，每个子集指一天中的一个小时。这样做的效果是使信号正则化和平滑化，从而使预测更加简单。
然后对每个子集进行特征设计、训练和微调。最终的预测将综合这 24 个模型的结果。也就是说，每个小时模型都有其特殊性，其中大部分都会考虑到重要的滞后因素。
在继续讨论之前，让我们先定义一下在进行滞后分析时可以处理的两种滞后类型：
- 自回归滞后期：接近滞后期 0 的滞后期，我们预期其值较高（最近的滞后期更有可能预测现值）。它们代表了序列的趋势程度。
- 季节滞后期：指季节性的滞后期。当按小时分割数据时，它们通常代表每周的季节性。
请注意，自动回归滞后期 1 也可以作为序列的**日季节性滞后期。
现在我们来讨论一下上面打印的 PACF 图。
- 夜间时间夜间时间（0，4）的消费更依赖于自回归滞后期而不是周滞后期，因为最重要的都集中在前五个滞后期。7、14、21、28 等季节性时段似乎不太重要，这建议我们在进行特征工程时特别关注滞后期 1 至 5。
- 日时日间时段（8、12、16、20）的消费量表现出自回归和季节性滞后。这一点在 8 和 12 小时尤为明显，因为这两个小时的消费量特别高，而在临近夜间时，季节滞后就变得不那么重要了。对于这些子集，我们还应该包括季节滞后和自回归滞后。
最后，这里有一些工程滞后特征的提示：
- 不要考虑过多的滞后期，因为这可能会导致过度拟合。一般来说，自动回归滞后期为 1 至 7，而周滞后期应为 7、14、21 和 28。但并不是一定要把每个滞后期都作为特征。
- 考虑非自动回归或季节性滞后通常是个坏主意，因为它们也可能带来过度拟合。相反，应尽量理解某个滞后期的重要性。
- 对滞后期进行转换通常可以获得更强大的特征。例如，可以使用加权平均值对季节性滞后进行聚合，以创建代表序列季节性的单一特征。
写在最后
本文构建了一个全面的探索性数据分析框架、旨在为时间序列预测提供参考。
探索性数据分析是数据科学研究的基础步骤、能够揭示数据的本质特征、为后续特征工程奠定基础、从而提高模型性能。
我们介绍了常用的时间序列EDA方法、包括统计/数学分析和可视化分析。该框架仅供参考、实际应用需要根据具体的时间序列类型和业务场景进行适当调整和扩展。
