Hey guys, 第七篇就來實作一遍，「以傳統統計方法」預測多變量時間序列吧

雖然 VAR 的準確度和複雜度不像我們後面天數要介紹的神經網絡一樣，但這個實作流程的重點，我覺得在於學習「更了解時間序列的特性和預測的關係」，比方說，平穩性、從時間序列中找自迴歸係數等。
對基本的東西更熟悉可以避免說，只是會使用神經網絡，但不了解時間序列的特性。

像我一開始寫這個系列也在掙扎說，要不要直接從機器學習迴歸開始講，或是直接介紹各種神經網絡、Attention is all you need? 不是嗎XD
後來還是覺得應該從基礎開始，即便多半是統計和數學，也至少確認自己是理解基本理論的。

好了，不廢話了，開始我們的 Python 實作吧

今天的大綱如下：

• 資料集介紹
  • 預測目標
• 時間序列平穩性檢測、原因
• 訓練/測試集切分
• 自迴歸模型係數搜尋
• 建立模型
• 預測並判讀結果

資料集介紹

我們使用 R 語言的一份範例資料集 Raotbl6，來自於套件庫 urca (Unit Root and Cointegration Tests for Time Series Data，一個計量經濟學的常用R套件)。這個資料集也是 Yash P. Mehra 曾在其 1994 年論文中使用的資料。

載入資料集

dataset = 'https://raw.githubusercontent.com/selva86/datasets/master/Raotbl6.csv'
df = pd.read_csv(dataset, parse_dates=['date'], index_col='date')
print(df.shape)
df.head()

• rgnp: 實際國民生產總值(real GNP)
• pgnp: 潛在國民生產總值(potential GNP)
• ulc: 單位勞動成本(unit labor cost)
• gdfco: 核心通膨率(core PCE)（剔除季節性食品及能源價格後的個人消費支出物價指數）：一個國家在不同時期內，個人消費支出總水平變動程度的經濟指數
• gdf: 國民生產總值(GNP) 縮減指數
• gdfim: 進口縮減指數
• gdfcf: 食品價格通膨率
• gdfce: 能源價格通膨率

預測目標

瞭解這 8 種經濟指標之間的關係，以及納入指標間的相互影響，預測未來數值

時間序列平穩性檢測

何謂平穩性（Stationarity）？

• 資料分布不會隨著時間改變而改變，平均數與變異數維持固定，例如白噪音
• 大部分時間序列需要去掉趨勢性、做 differencing 等動作，平穩性才會顯現

為什麼需要檢測時間序列的平穩性？

• 自迴歸模型中，僅以歷史變量和當前變量預測未來；因此需要時間序列變量的基本特性能長時間維持不變，否則以過去預測未來的思路就不成立了。

檢測平穩性的方法？

EX: Unit Root Tests

• 單位根檢驗用於測試時間序列是否非平穩，並具有單位根
• 虛無假設為，假定時間序列存在 unit root（表示是 non-stationary），計算 p-value 驗證是否推翻虛無假設

以下我們使用 Augmented Dickey-Fuller Test 進行平穩性檢測：

def adf_test(series, title=''):
    print(f'Augmented Dickey-Fuller Test: {title}')
    result = adfuller(series.dropna(), autolag='AIC') # .dropna() handles differenced data
    labels = ['ADF test statistic', 'p-value', 'Number of lags used', 'Number of observations used']
    out = pd.Series(result[0:4], index=labels)
    for k, v in result[4].items():
        out[f'critical value ({k})'] = v

    print(out)
    if result[1] <= 0.05:  # 有顯著性，推翻虛無假設
        print("Data has no unit root and is stationary")
    else:
        print("Data has a unit root and is non-stationary")

for col in df.columns:
    adf_test(df[col], title=col)
    print()

幾乎所有變量都是不平穩序列，differencing!

# differencing(1)
df_differenced = df.diff().dropna()

# do ADF test on differencing(1) dataframe
for col in df_differenced.columns:
    adf_test(df_differenced[col], title=col)
    print()

還是有部分時間序列變量不平穩，再做一次 differencing
（自動一點的方法可以寫成當所有時間序列變量的 p-value 都 ≤ 0.05 時停止 differencing，保存做過差分的次數；這邊為了快速示範，就沒有特別寫～）

# Second Differencing
df_differenced = df_differenced.diff().dropna()

# do ADF test on differencing(2) dataframe
for col in df_differenced.columns:
    adf_test(df_differenced[col], title=col)
    print()

所有時間序列變量都符合平穩性了！

訓練/測試集切分

接著我們把要用來 forecasting 的部分切出來保存

steps = 4
train, test = df_differenced[0:-steps], df_differenced[-steps:]
print(train.shape) 
print(test.shape) 

自迴歸模型係數搜尋

接著搜尋自迴歸係數 p：

for i in range(1, 11):
    model = VAR(train)
    results = model.fit(i)
    print('Number of Lag  = ', i)
    print('AIC: ', results.aic)
    print('BIC: ', results.bic, '\n')

根據數據，選擇 AIC 及 BIC 最小值所在的 Lag = 1

P.S. 這邊不一定要用迴圈找，也可以使用 VAR model 的 model.select_order(maxlags) 搜尋最佳 order(p)：

model.select_order(maxlags=12).summary()

建立模型

建立 order(p=1) 的自迴歸模型：

model_fitted = model.fit(1)
model_fitted.summary()

這邊呈現了模型的擬合情況和各項係數

預測並判讀結果

預測

lag_num = model_fitted.k_ar
test_x = df_differenced.values[-lag_num:]

# forecasting
pred = model_fitted.forecast(y=test_x, steps=steps)
# 加 "_2d" 是為了提醒，預測出的結果是經過 2 次 differencing 的
df_pred = pd.DataFrame(pred, index=df.index[-steps:], columns=df.columns + '_2d')

# 將預測結果還原
for col in df.columns:
    df_pred[f'{col}_1d'] = (df[col].iloc[-steps-1]-df[col].iloc[-steps-2]) + df_pred[f'{col}_2d'].cumsum()
    df_pred[f'{col}_forecast'] = df[col].iloc[-steps-1] + df_pred[f'{col}_1d'].cumsum()
    
test_original = df[-steps:]
test_original.index = pd.to_datetime(test_original.index)

將預測結果和實際畫出比較：

fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(15,10))
for i, ax in enumerate(axes.flatten()):
  ax.plot(test_original[df.columns[i]], color='blue', linewidth=1)
  ax.plot(df_forecast[f'{df.columns[i]}_Forecast'], color='red', linewidth=1)
  ax.set_title(df.columns[i])
  ax.spines['top'].set_alpha(0)
  ax.tick_params(labelsize=10)
  ax.legend(['actual', 'forecast'])
fig.tight_layout();

根據圖表，我們發現在 rgnp 和 pgnp 的預測效果較好

你也可以引入 R2, RMSE, MSE 等指標來呈現成效好壞～

參考資料

• 資料序列stationary檢定方法
• 如何理解时间序列的平稳性？
• Unit root test

好的，總結本篇教學重點：

1. 瞭解如何使用 python 套件建立自迴歸模型
2. 瞭解如何利用資訊指標搜尋最佳自迴歸係數
3. 瞭解何謂時間序列平穩性，以及為什麼傳統計量經濟上需要先做平穩化後，才能預測

感謝你的收看！明天見！