常用的時間序列分析方法總結和代碼示例
時間序列是最流行的數據類型之一。視頻，圖像，像素，信號，任何有時間成分的東西都可以轉化為時間序列。
在本文中將在分析時間序列時使用的常見的處理方法。這些方法可以幫助你獲得有關數據本身的見解，為建模做好準備并且可以得出一些初步結論。
我們將分析一個氣象時間序列。利用逐時ERA5 Land[1]研究2023年西伯利亞東南部點的2 m氣溫、總降水量、地表凈太陽輻射和地表壓力。
首先我們導入相關的庫：
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from scipy import statsmatplotlib是可以設置不同的風格的，這里我們使用 opinionated和 ambivalent來進行風格的設置
from ambivalent import STYLES
import opinionated
plt.style.use(STYLES['ambivalent'])
plt.style.use("dark_background")折線圖
要觀察一個時間序列，最簡單的方法就是折線圖。為了處理地理空間多維數組，我們將使用xarray庫。
data = xr.open_dataset('Medium_data.nc')
data
現在我們需要針對所選位置對數據進行切片，并將其轉換為pandas DF，并創建一個線形圖:
df = data.sel(latitude=52.53, lnotallow=101.63, method='pad').to_pandas().drop(['latitude', 'longitude'], axis=1)
fig, ax = plt.subplots(ncols = 2, nrows = 2, figsize=(16,9))
df['t2m'].plot(ax=ax[0,0])
ax[0,0].set_title('Air Temperature')
df['ssr'].plot(ax=ax[0,1])
ax[0,1].set_title('Surface Net Solar Radiation')
df['sp'].plot(ax=ax[1,0])
ax[1,0].set_title('Surface Pressure')
df['tp'].plot(ax=ax[1,1])
ax[1,1].set_title('Total Precipitation')
plt.tight_layout()
plt.show()
從線形圖中可以清楚地看出，所有四個時間序列都有不同的特征，下面讓我們使用數學工具來研究它們。
分解與平穩性
任何時間序列都有三個重要屬性需要考慮:
1、趨勢是時間序列中平穩的長期變化;
2、季節性指的是一個時間序列的平均值有規律的周期性變化;
3、噪聲(殘差)，它是均值為零的信號的隨機成分。
為了分別得到這些成分，可以使用經典分解(加性或乘法)。該操作是通過應用卷積濾波器產生的，因此每個時間序列分量被定義為
或者
這里的y為時間序列的值，S為季節分量，T為趨勢分量，n為噪聲。
為了進行分解，除了選擇分解類之外，還需要設置一個季節周期(例如，p=1表示年度數據，p=4表示季度數據，p=12表示月度數據等)。
前面提到的經典分解是一種非常幼稚和簡單的方法。它具有明顯的局限性，如線性，無法捕捉動態季節性和難以處理時間序列中的非平穩性，但是就本文作為演示，這種方法是可以的。
為了進行經典的分解，我們將使用statmodels庫中的seasonal_decomposition函數，周期等于24，因為我們處理的是每小時的數據:
vars = {'t2m': 'Air Temperature', 'tp': 'Total Precipitation', 'sp': 'Surface Pressure', 'ssr': 'Surface Net Solar Radiation'}
for var in df.columns:
result = sm.tsa.seasonal_decompose(df[var], model='additive', period = 24)
results_df = pd.DataFrame({'trend': result.trend, 'seasonal': result.seasonal, 'resid': result.resid, 'observed': result.observed})
fig, ax = plt.subplots(ncols = 2, nrows = 2,figsize=(16,9))
ax[0,0].plot(df.index, results_df.trend)
ax[0,0].set_title('Trend')
ax[0,0].set_ylabel('Value')
ax[0,1].plot(df.index, results_df.seasonal)
ax[0,1].set_title('Seasonal')
ax[1,0].plot(df.index, results_df.resid)
ax[1,0].set_title('Residual')
ax[1,0].set_ylabel('Value')
ax[1,0].set_xlabel('time')
ax[1,1].plot(df.index, results_df.observed)
ax[1,1].set_title('Observed')
ax[1,1].set_xlabel('time')
opinionated.set_title_and_suptitle(vars[var], f"Dickey-Fuller test: {round(sm.tsa.stattools.adfuller(df[var])[1],5)}", position_title=[0.45,1],
position_sub_title=[0.95, 1])
plt.tight_layout()
plt.savefig(f'Seasonal_{var}.png')
plt.show()
你可以看到，對于所有的變量，季節性因素看起來都很混亂。這是因為我們分析的是每小時的數據，這些季節變化是在一天內觀察到的，并沒有直接的關聯。所以我們可以嘗試將數據重新采樣到每日間隔，并在一天的時間段內進行分解。
df_d = df.resample('1d').mean()
請注意到圖表右上角的Dickey-Fuller(ADF) 。這是一個平穩性測試，使用的是adfuller函數。對于時間序列，平穩性意味著時間序列的屬性不隨時間變化。我們這里說的屬性是指：方差、季節性、趨勢和自相關性。
Dickey-Fuller (ADF)檢驗的流程是：提出時間序列是非平穩的零假設。然后我們選擇顯著性水平α，通常為5%。α是錯誤地拒絕零假設的概率，而零假設實際上是正確的。所以在我們的例子中，α=5%有5%的風險得出時間序列是平穩的，而實際上不是。
測試結果會給出一個p值。如果小于0.05，我們可以拒絕零假設。可以看到，根據ADF檢驗所有4個變量都是平穩的。
一般情況下要應用時間序列預測模型，如ARIMA等，平穩性是必須的。這也是我們選擇氣象數據的原因，因為它們在大多數情況下是平穩的，所以才會出現在不同的時間序列相關的學習材料中進行分析。
分布
在得出所有時間序列都是平穩的結論之后，讓我們來看看它們是如何分布的。我們將使用著名的seaborn庫及其函數pairplot，該函數允許使用歷史和kde創建信息豐富的圖。
ax = sns.pairplot(df, diag_kind='kde')
ax.map_upper(sns.histplot, bins=20)
ax.map_lower(sns.kdeplot, levels=5, color='.1')
plt.show()
讓我們考慮t2m(1行1列)的示例。在分析核密度估計(kde)圖時，很明顯這個變量的分布是多模態的，這意味著它由2個或更多的“鐘形”組成。在本文的后續階段中，我們將嘗試將變量轉換為類似于正態分布的形式。
第一列和第一行中的其他圖是相同的，但它們的可視化方式不同。這些是散點圖，可以確定兩個變量是如何相關的。所以一個點的顏色越深，或者離中心圓越近，這個區域內點的密度就越高。
Box-Cox轉換
由于我們已經發現氣溫時間序列是平穩的，但不是正態分布，所以可以嘗試使用Box-Cox變換來修復它。這里使用scipy包及其函數boxcox。
df_d['t2m_box'], _ = stats.boxcox(df_d.t2m)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(15,7))
sns.histplot(df_d.t2m_box, kde=True, ax=ax[0])
sns.histplot(df_d.t2m, kde=True, ax=ax[1])
圖的左邊部分是經過BoxCox變換后的時間序列分布，可以看到，它還遠遠不能被稱為“正態”分布。但是如果我們把它和右邊的比較，我們可以說的確更接近于“正態”。
我們還可以做的另一件事是確保執行的轉換是有用的，可以創建一個概率圖：繪制理論分布的分位數(在我們的情況下是正態)與經驗數據的樣本(即我們考慮的時間序列)。越靠近白線的點越好。
fig = plt.figure()
ax1 = fig.add_subplot(211)
prob = stats.probplot(df_d.t2m, dist=stats.norm, plot=ax1)
ax1.get_lines()[1].set_color('w')
ax1.get_lines()[0].set_color('#8dd3c7')
ax1.set_title('Probplot against normal distribution')
ax2 = fig.add_subplot(212)
prob = stats.probplot(df_d.t2m_box, dist=stats.norm, plot=ax2)
ax2.get_lines()[1].set_color('w')
ax2.get_lines()[0].set_color('#8dd3c7')
ax2.set_title('Probplot after Box-Cox transformation')
plt.tight_layout()fig = plt.figure()
ax1 = fig.add_subplot(211)
prob = stats.probplot(df_d.t2m, dist=stats.norm, plot=ax1)
ax1.set_title('Probplot against normal distribution')
ax2 = fig.add_subplot(212)
prob = stats.probplot(df_d.t2m_box, dist=stats.norm, plot=ax2)
ax2.set_title('Probplot after Box-Cox transformation')
plt.tight_layout()
這個概率圖還有一個更常見的名字QQ圖
另外需要說明的是，如果打算使用轉換后的時間序列進行ML建模，不要忘記應用反向BoxCox轉換，這樣才能的到最終的正確結果。
自相關
時間序列分析的最后一步是自相關。自相關函數(ACF)估計時間序列和滯后版本之間的相關性?；蛘邠Q句話說，時間序列的特定值如何與不同時間間隔內的其他先驗值相關聯。繪制部分自相關函數(PACF)也可能有所幫助，它與自相關相同，但刪除了較短滯后的相關性。它估計某個時間戳內值之間的相關性，但控制其他值的影響。
for var in df.columns[:-1]:
fig, (ax1, ax2) = plt.subplots(2,1,figsize=(10,8))
plot_acf(df_d.t2m, ax = ax1)
plot_pacf(df_d.t2m, ax = ax2)
opinionated.set_title_and_suptitle(vars[var], '',position_title=[0.38,1],
position_sub_title=[0.95, 1])
plt.tight_layout()
plt.show()
可以看到在地表壓力時間序列中有一個非常強的部分自相關，有1天的滯后。然后明顯減弱，3天后幾乎消失。這樣的分析可以幫助我們更好地理解正在處理的數據的性質，從而得出更有意義的結論。
總結
以上就是在處理時間序列時進行探索性數據分析時常用的方法，通過上面這些方法可以很好的了解到時間序列的信息，為我們后面的建模提供數據的支持。
