简介

时间序列分析无疑是数据科学和机器学习领域最广泛的主题之一：无论是预测金融事件、能源消耗、产品销售还是股市趋势，这一领域一直是企业非常感兴趣的领域。

显然，数据可用性的大幅提高，加上机器学习模型的不断进步，使得这一话题如今更加引人关注。除了源自统计的传统预测方法（如回归模型、ARIMA 模型、指数平滑法）外，与机器学习（如基于树的模型）和深度学习（如 LSTM 网络、CNN、基于变换器的模型）相关的技术也已出现了一段时间。

尽管这些技术之间存在巨大差异，但无论采用何种模型，都必须完成一个初步步骤： 探索性数据分析。

在统计学中，探索性数据分析（EDA）是一门对数据进行分析和可视化的学科，目的是总结数据的主要特征并从中获取相关信息。这在数据科学领域相当重要，因为它可以为另一个重要步骤：特征工程奠定基础。也就是从数据集中创建、转换和提取特征，从而使模型发挥最大作用的实践。

因此，本文的目的是定义一个清晰的探索性数据分析模板，该模板侧重于时间序列，可以总结和突出数据集的最重要特征。为此，我们将使用一些常用的 Python 库，如 Pandas、Seaborn 和 Statsmodel。

数据

首先让我们定义一下数据集：在本文中，我们将使用 Kaggle 的每小时能源消耗数据。该数据集与 PJM 每小时能源消耗数据有关，PJM 是美国的一个区域输电组织，为特拉华州、伊利诺伊州、印第安纳州、肯塔基州、马里兰州、密歇根州、新泽西州、北卡罗来纳州、俄亥俄州、宾夕法尼亚州、田纳西州、弗吉尼亚州、西弗吉尼亚州和哥伦比亚特区提供电力服务。

每小时电力消耗数据来自 PJM 网站，单位为兆瓦 (MW)。

探索性数据分析

现在让我们来定义一下在处理时间序列时需要进行的最重要的分析。

当然，最重要的一点是绘制数据图：图表可以突出许多特征，如模式、异常观测、随时间的变化以及变量之间的关系。如前所述，从这些图表中得出的见解必须尽可能纳入预测模型。此外，一些数学工具，如描述性统计和时间序列分解，也会非常有用。

我在本文中提出的 EDA 包括六个步骤： 描述性统计、时间图、季节图、箱形图、时间序列分解、滞后分析。

1. 描述性统计

描述性统计是一种总结性统计，用于定量描述或总结结构化数据集合的特征。

常用于描述数据集的指标有：中心倾向度量（如平均值、中位数）、离散度量（如范围、标准差）和位置度量（如百分位数、四分位数）。所有这些都可以用所谓的五数汇总来概括，其中包括：分布的最小值、第一四分位数 (Q1)、中位数或第二四分位数 (Q2)、第三四分位数 (Q3) 和最大值。

在 Python 中，可以使用 Pandas 中众所周知的 describe 方法轻松获取这些信息：

import pandas as pd
# Loading and preprocessing steps
df = pd.read_csv('../input/hourly-energy-consumption/PJME_hourly.csv')
df = df.set_index('Datetime')
df.index = pd.to_datetime(df.index)
df.describe()

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')

这幅图已经提供了一些信息：

1. 正如我们所预料的那样，这种模式显示出每年的季节性。
2. 如果只看一年，似乎会发现更多的规律。由于耗电量较大，冬季和夏季的耗电量很可能各有一个高峰。
3. 多年来，该系列并没有表现出明显的增减趋势，平均消耗量保持稳定。
4. 在 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, fontsize=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, fontsize=12, color=colors[i])
# Setting labels
plt.gca().set(ylabel= 'PJME_MW', xlabel = 'Month')
plt.yticks(fontsize=12, alpha=.7)
plt.title("Seasonal Plot - Monthly Consumption", fontsize=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Month')
plt.show()

从图中可以看出，每年的消耗量实际上都有一个非常固定的模式：冬季消耗量明显增加，夏季消耗量达到峰值（由于供暖/制冷系统），而春季和秋季通常不需要供暖或制冷，因此消耗量为最小值。

此外，这幅图还告诉我们，不同年份的总体消耗量并没有明显的增减模式。

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.month==y, :].shape[0]-.9, df_plot.loc[df_plot.month==y, 'PJME_MW'][-1:].values[0], y, fontsize=12, color=colors[i])
        else:
            plt.text(df_plot.loc[df_plot.month==y, :].shape[0]-.9, df_plot.loc[df_plot.month==y, 'PJME_MW'][-1:].values[0], y, fontsize=12, color=colors[i])

# Setting Labels
plt.gca().set(ylabel= 'PJME_MW', xlabel = 'Month')
plt.yticks(fontsize=12, alpha=.7)
plt.title("Seasonal Plot - Weekly Consumption", fontsize=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Month')
plt.show()

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", fontsize=20)
plt.ylabel('Consumption [MW]')
plt.xlabel('Hour')
plt.legend()

通常情况下，这种曲线图会显示出一种非常典型的模式，有人称之为 "M 型曲线"，因为白天的消耗量似乎呈 "M "型。有时这种模式很明显，有时则不明显（比如本例）。

不过，这种曲线图通常会在一天的中间（上午 10 点到下午 2 点）出现一个相对的峰值，然后是一个相对的峰值（下午 2 点到下午 6 点）和另一个峰值（下午 6 点到晚上 8 点）。最后，它还显示了周末和其他日子的消耗量差异。

3.4 季节图--特征工程

现在我们来看看如何将这些信息用于特征工程。假设我们正在使用一些需要高质量特征的 ML 模型（如 ARIMA 模型或基于树的模型）。

这些是来自季节图的主要证据：

1. 年度消费在不同年份变化不大：这表明可以使用来自滞后变量或外生变量的年度季节性特征。
2. 各月的周消费遵循相同的模式：这表明可以使用来自滞后变量或外生变量的周特征。
3. 日常消费与平日和周末不同：这表明可以使用分类特征来识别平日和非平日。

4. 方框图

方框图是确定数据分布情况的有效方法。简而言之，方框图描述了百分位数（代表分布的第 1 个四分位数（Q1）、第 2 个四分位数（Q2/中位数）和第 3 个四分位数（Q3））和边线（代表数据的范围）。超出晶须的每一个值都可以被视为离群值，更深入地说，晶须的计算方法通常是：

4.1 箱形图 - 总消耗量

我们首先来计算总消耗量的箱形图，这可以通过 Seaborn 轻松完成：

plt.figure(figsize=(8,5))8,5))
sns.boxplot(data=df, x='PJME_MW')
plt.xlabel('Consumption [MW]')
plt.title(f'Boxplot - Consumption Distribution');

尽管这幅图看起来信息量不大，但它告诉我们，我们正在处理的是一个类似高斯分布的 分布，其尾部更偏向右侧。

4.2 箱形图--日月分布

日/月箱形图是一个非常有趣的图。通过创建一个 "日-月 "变量，并根据该变量对消耗量进行分组，可以得到该图。以下是代码，仅涉及 2017 年：

df['year'] = [x for x in df.index.year]'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(rotation=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')

可以看出，夏/冬两季（即气温达到峰值时）的消费量不确定性较小，而春/秋两季（即气温变化较大时）的消费量较为分散。最后，2018 年夏季的消费量高于 2017 年，这可能是由于夏季较为温暖的缘故。在进行特征工程设计时，请记住将温度曲线（如果有的话）包括在内，因为它可能可以用作外生变量。

4.3 方框图--日分布

另一种有用的曲线图是指一周内的消耗量分布，这与周消耗量季节曲线图类似。

df_plot = df[['day_str', 'day', 'PJME_MW']].sort_values(by='day')'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.4 箱形图--小时分布

最后让我们来看看小时分布箱形图。它与每日消费季节图类似，因为它提供了消费在一天中的分布情况。代码如下

plt.title(f'Boxplot Hour Distribution');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')

请注意，之前看到的 "M "形现在变得更粗了。此外，还有很多离群值，这说明数据不仅依赖于每日的季节性（例如，今天上午 12 点的消耗量与昨天上午 12 点的消耗量相似），还依赖于其他因素，可能是温度或湿度等外生气候特征。

5. 时间序列分解

如前所述，时间序列数据可以表现出多种模式。通常情况下，将一个时间序列拆分成几个部分是有帮助的，每个部分代表一个基本模式类别。

我们可以把时间序列看作由三个部分组成：趋势部分、季节部分和剩余部分（包含时间序列中的任何其他部分）。对于某些时间序列（如能源消耗序列），可能存在不止一个季节性成分，分别对应于不同的季节性时期（日、周、月、年）。

分解方法主要有两种：加法分解和乘法分解。

对于加法分解，我们将一个序列（?）表示为季节成分（?）、趋势（?）和余数（?）的总和：

同样，乘法分解可以写成

一般来说，加法分解最能代表方差恒定的序列，而乘法分解最适合方差非平稳的时间序列。

在 Python 中，使用 Statsmodel 库可以轻松实现时间序列分解：

df_plot = df[df['year'] == 2017].reset_index()'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('', fontsize=22)
plt.xticks(rotation=45)
result_mul.plot().suptitle('', fontsize=22)
plt.xticks(rotation=45)
plt.show()

上图是 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('', fontsize=22)
plt.xticks(rotation=45)
result_mul.plot().suptitle('', fontsize=22)
plt.xticks(rotation=45)
plt.show()

6. 滞后分析

在时间序列预测中，滞后期就是序列的过去值。例如，对于日序列，第一个滞后值指的是该序列前一天的值，第二个滞后值指的是前一天的值，以此类推。

滞后分析的基础是计算序列与序列本身的滞后版本之间的相关性，这也称为自相关性。对于一个 k 个滞后版本的序列，我们将自相关系数定义为

其中，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()

正如你所看到的，PACF 简单来说就是绘制不同滞后期的皮尔逊部分自相关系数。当然，非滞后序列本身显示出完美的自相关性，因此滞后 0 总是 1。蓝色区域代表置信区间：如果滞后期超过该区域，则在统计上具有显著性，我们可以断言其具有重要意义。

6.1 滞后分析--特征工程

滞后分析是对时间序列特征工程影响最大的研究之一。如前所述，相关性高的滞后期是序列的重要滞后期，因此应加以考虑。

一种广泛使用的特征工程技术是按小时划分数据集。也就是说，将数据分成 24 个子集，每个子集指一天中的一个小时。这样做的效果是使信号正则化和平滑化，从而使预测更加简单。

然后对每个子集进行特征设计、训练和微调。最终的预测将综合这 24 个模型的结果。也就是说，每个小时模型都有其特殊性，其中大部分都会考虑到重要的滞后因素。

在继续讨论之前，我们先来定义一下滞后分析中的两种滞后情况：

1. 自回归滞后期：接近滞后期 0 的滞后期，我们预期其值较高（最近的滞后期更有可能预测现值）。自回归滞后期代表了序列的趋势程度。
2. 季节滞后期：滞后期指的是季节性时期。当按小时分割数据时，它们通常代表每周的季节性。

现在我们来讨论上面打印的 PACF 图。

夜间时间

夜间时间（0，4）的消费更依赖于自回归滞后期而不是周滞后期，因为最重要的都集中在前五个滞后期。7、14、21、28 等季节性时段似乎不太重要，这建议我们在进行特征工程时特别关注滞后期 1 至 5。

日间时段

日间时段（8、12、16、20）的消耗量既表现出自回归滞后，也表现出季节性滞后。这一点在 8 和 12 小时尤为明显，因为这两个小时的消费量特别高，而在临近夜间时，季节性滞后就变得不那么重要了。对于这些子集，我们还应该包括季节滞后和自回归滞后。

结论

本文旨在为时间序列预测提供一个全面的探索性数据分析模板。

在任何类型的数据科学研究中，EDA 都是一个基本步骤，因为它可以了解数据的性质和特殊性，并为特征工程奠定基础，而特征工程反过来又可以显著提高模型性能。