時(shí)間序列是最流行的數(shù)據(jù)類型之一。視頻，圖像，像素，信號(hào)，任何有時(shí)間成分的東西都可以轉(zhuǎn)化為時(shí)間序列。

在本文中將在分析時(shí)間序列時(shí)使用的常見(jiàn)的處理方法。這些方法可以幫助你獲得有關(guān)數(shù)據(jù)本身的見(jiàn)解，為建模做好準(zhǔn)備并且可以得出一些初步結(jié)論。

我們將分析一個(gè)氣象時(shí)間序列。利用逐時(shí)ERA5 Land[1]研究2023年西伯利亞?wèn)|南部點(diǎn)的2 m氣溫、總降水量、地表凈太陽(yáng)輻射和地表壓力。

首先我們導(dǎo)入相關(guān)的庫(kù)：

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 stats

matplotlib是可以設(shè)置不同的風(fēng)格的，這里我們使用 opinionated和 ambivalent來(lái)進(jìn)行風(fēng)格的設(shè)置

from ambivalent import STYLES
 import opinionated
 plt.style.use(STYLES['ambivalent'])
 plt.style.use("dark_background")

折線圖

要觀察一個(gè)時(shí)間序列，最簡(jiǎn)單的方法就是折線圖。為了處理地理空間多維數(shù)組，我們將使用xarray庫(kù)。

data = xr.open_dataset('Medium_data.nc')
 data

現(xiàn)在我們需要針對(duì)所選位置對(duì)數(shù)據(jù)進(jìn)行切片，并將其轉(zhuǎn)換為pandas DF，并創(chuàng)建一個(gè)線形圖:

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

從線形圖中可以清楚地看出，所有四個(gè)時(shí)間序列都有不同的特征，下面讓我們使用數(shù)學(xué)工具來(lái)研究它們。

分解與平穩(wěn)性

任何時(shí)間序列都有三個(gè)重要屬性需要考慮:

1、趨勢(shì)是時(shí)間序列中平穩(wěn)的長(zhǎng)期變化;

2、季節(jié)性指的是一個(gè)時(shí)間序列的平均值有規(guī)律的周期性變化;

3、噪聲(殘差)，它是均值為零的信號(hào)的隨機(jī)成分。

為了分別得到這些成分，可以使用經(jīng)典分解(加性或乘法)。該操作是通過(guò)應(yīng)用卷積濾波器產(chǎn)生的，因此每個(gè)時(shí)間序列分量被定義為

或者

這里的y為時(shí)間序列的值，S為季節(jié)分量，T為趨勢(shì)分量，n為噪聲。

為了進(jìn)行分解，除了選擇分解類之外，還需要設(shè)置一個(gè)季節(jié)周期(例如，p=1表示年度數(shù)據(jù)，p=4表示季度數(shù)據(jù)，p=12表示月度數(shù)據(jù)等)。

前面提到的經(jīng)典分解是一種非常幼稚和簡(jiǎn)單的方法。它具有明顯的局限性，如線性，無(wú)法捕捉動(dòng)態(tài)季節(jié)性和難以處理時(shí)間序列中的非平穩(wěn)性，但是就本文作為演示，這種方法是可以的。

為了進(jìn)行經(jīng)典的分解，我們將使用statmodels庫(kù)中的seasonal_decomposition函數(shù)，周期等于24，因?yàn)槲覀兲幚淼氖敲啃r(shí)的數(shù)據(jù):

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

你可以看到，對(duì)于所有的變量，季節(jié)性因素看起來(lái)都很混亂。這是因?yàn)槲覀兎治龅氖敲啃r(shí)的數(shù)據(jù)，這些季節(jié)變化是在一天內(nèi)觀察到的，并沒(méi)有直接的關(guān)聯(lián)。所以我們可以嘗試將數(shù)據(jù)重新采樣到每日間隔，并在一天的時(shí)間段內(nèi)進(jìn)行分解。

df_d = df.resample('1d').mean()

請(qǐng)注意到圖表右上角的Dickey-Fuller(ADF) 。這是一個(gè)平穩(wěn)性測(cè)試，使用的是adfuller函數(shù)。對(duì)于時(shí)間序列，平穩(wěn)性意味著時(shí)間序列的屬性不隨時(shí)間變化。我們這里說(shuō)的屬性是指：方差、季節(jié)性、趨勢(shì)和自相關(guān)性。

Dickey-Fuller (ADF)檢驗(yàn)的流程是：提出時(shí)間序列是非平穩(wěn)的零假設(shè)。然后我們選擇顯著性水平α，通常為5%。α是錯(cuò)誤地拒絕零假設(shè)的概率，而零假設(shè)實(shí)際上是正確的。所以在我們的例子中，α=5%有5%的風(fēng)險(xiǎn)得出時(shí)間序列是平穩(wěn)的，而實(shí)際上不是。

測(cè)試結(jié)果會(huì)給出一個(gè)p值。如果小于0.05，我們可以拒絕零假設(shè)。可以看到，根據(jù)ADF檢驗(yàn)所有4個(gè)變量都是平穩(wěn)的。

一般情況下要應(yīng)用時(shí)間序列預(yù)測(cè)模型，如ARIMA等，平穩(wěn)性是必須的。這也是我們選擇氣象數(shù)據(jù)的原因，因?yàn)樗鼈冊(cè)诖蠖鄶?shù)情況下是平穩(wěn)的，所以才會(huì)出現(xiàn)在不同的時(shí)間序列相關(guān)的學(xué)習(xí)材料中進(jìn)行分析。

分布

在得出所有時(shí)間序列都是平穩(wěn)的結(jié)論之后，讓我們來(lái)看看它們是如何分布的。我們將使用著名的seaborn庫(kù)及其函數(shù)pairplot，該函數(shù)允許使用歷史和kde創(chuàng)建信息豐富的圖。

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列)的示例。在分析核密度估計(jì)(kde)圖時(shí)，很明顯這個(gè)變量的分布是多模態(tài)的，這意味著它由2個(gè)或更多的“鐘形”組成。在本文的后續(xù)階段中，我們將嘗試將變量轉(zhuǎn)換為類似于正態(tài)分布的形式。

第一列和第一行中的其他圖是相同的，但它們的可視化方式不同。這些是散點(diǎn)圖，可以確定兩個(gè)變量是如何相關(guān)的。所以一個(gè)點(diǎn)的顏色越深，或者離中心圓越近，這個(gè)區(qū)域內(nèi)點(diǎn)的密度就越高。

Box-Cox轉(zhuǎn)換

由于我們已經(jīng)發(fā)現(xiàn)氣溫時(shí)間序列是平穩(wěn)的，但不是正態(tài)分布，所以可以嘗試使用Box-Cox變換來(lái)修復(fù)它。這里使用scipy包及其函數(shù)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])

圖的左邊部分是經(jīng)過(guò)BoxCox變換后的時(shí)間序列分布，可以看到，它還遠(yuǎn)遠(yuǎn)不能被稱為“正態(tài)”分布。但是如果我們把它和右邊的比較，我們可以說(shuō)的確更接近于“正態(tài)”。

我們還可以做的另一件事是確保執(zhí)行的轉(zhuǎn)換是有用的，可以創(chuàng)建一個(gè)概率圖：繪制理論分布的分位數(shù)(在我們的情況下是正態(tài))與經(jīng)驗(yàn)數(shù)據(jù)的樣本(即我們考慮的時(shí)間序列)。越靠近白線的點(diǎn)越好。

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

這個(gè)概率圖還有一個(gè)更常見(jiàn)的名字QQ圖

另外需要說(shuō)明的是，如果打算使用轉(zhuǎn)換后的時(shí)間序列進(jìn)行ML建模，不要忘記應(yīng)用反向BoxCox轉(zhuǎn)換，這樣才能的到最終的正確結(jié)果。

自相關(guān)

時(shí)間序列分析的最后一步是自相關(guān)。自相關(guān)函數(shù)(ACF)估計(jì)時(shí)間序列和滯后版本之間的相關(guān)性?；蛘邠Q句話說(shuō)，時(shí)間序列的特定值如何與不同時(shí)間間隔內(nèi)的其他先驗(yàn)值相關(guān)聯(lián)。繪制部分自相關(guān)函數(shù)(PACF)也可能有所幫助，它與自相關(guān)相同，但刪除了較短滯后的相關(guān)性。它估計(jì)某個(gè)時(shí)間戳內(nèi)值之間的相關(guān)性，但控制其他值的影響。

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

可以看到在地表壓力時(shí)間序列中有一個(gè)非常強(qiáng)的部分自相關(guān)，有1天的滯后。然后明顯減弱，3天后幾乎消失。這樣的分析可以幫助我們更好地理解正在處理的數(shù)據(jù)的性質(zhì)，從而得出更有意義的結(jié)論。

總結(jié)

以上就是在處理時(shí)間序列時(shí)進(jìn)行探索性數(shù)據(jù)分析時(shí)常用的方法，通過(guò)上面這些方法可以很好的了解到時(shí)間序列的信息，為我們后面的建模提供數(shù)據(jù)的支持。