今天云朵君和大家一起學習如何使用探索性數據分析從時間序列數據中獲取信息，并使用 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 模型或基于樹的模型）。主要的證據來自季節圖包括以下幾點：

1. 年度消耗量在不同年份之間的變化不大，這意味著可以利用年度季節性特征，例如滯后變量或外生變量。
2. 周消費量在各月份中的變化規律相似，這表明可以利用周特征，如滯后變量或外生變量。
3. 日常消費與平日和周末有所不同，因此應當使用分類特征來區分平日和非平日。

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 個模型的結果。也就是說，每個小時模型都有其特殊性，其中大部分都會考慮到重要的滯后因素。

在繼續討論之前，讓我們先定義一下在進行滯后分析時可以處理的兩種滯后類型：

1. 自回歸滯后期：接近滯后期 0 的滯后期，我們預期其值較高（最近的滯后期更有可能預測現值）。它們代表了序列的趨勢程度。
2. 季節滯后期：指季節性的滯后期。當按小時分割數據時，它們通常代表每周的季節性。

請注意，自動回歸滯后期 1 也可以作為序列的**日季節性滯后期。

現在我們來討論一下上面打印的 PACF 圖。

1. 夜間時間夜間時間（0，4）的消費更依賴于自回歸滯后期而不是周滯后期，因為最重要的都集中在前五個滯后期。7、14、21、28 等季節性時段似乎不太重要，這建議我們在進行特征工程時特別關注滯后期 1 至 5。
2. 日時日間時段（8、12、16、20）的消費量表現出自回歸和季節性滯后。這一點在 8 和 12 小時尤為明顯，因為這兩個小時的消費量特別高，而在臨近夜間時，季節滯后就變得不那么重要了。對于這些子集，我們還應該包括季節滯后和自回歸滯后。

最后，這里有一些工程滯后特征的提示：

• 不要考慮過多的滯后期，因為這可能會導致過度擬合。一般來說，自動回歸滯后期為 1 至 7，而周滯后期應為 7、14、21 和 28。但并不是一定要把每個滯后期都作為特征。
• 考慮非自動回歸或季節性滯后通常是個壞主意，因為它們也可能帶來過度擬合。相反，應盡量理解某個滯后期的重要性。
• 對滯后期進行轉換通?？梢垣@得更強大的特征。例如，可以使用加權平均值對季節性滯后進行聚合，以創建代表序列季節性的單一特征。

寫在最后

本文構建了一個全面的探索性數據分析框架、旨在為時間序列預測提供參考。

探索性數據分析是數據科學研究的基礎步驟、能夠揭示數據的本質特征、為后續特征工程奠定基礎、從而提高模型性能。

我們介紹了常用的時間序列EDA方法、包括統計/數學分析和可視化分析。該框架僅供參考、實際應用需要根據具體的時間序列類型和業務場景進行適當調整和擴展。