下圖顯示了關于不同類型葡萄酒銷量的月度多元時間序列。每種葡萄酒類型都是時間序列中的一個變量。

假設要預測其中一個變量。比如，sparkling wine。如何建立一個模型來進行預測呢？

一種常見的方法是將該變量其視為單變量時間序列。這樣就有很多方法可以用來模擬這些系列。比如 ARIMA、指數平滑或 Facebook 的 Prophet，還有自回歸的機器學習方法也可以使用。

但是其他變量可能包含sparkling wine未來銷售的重要線索。看看下面的相關矩陣。

可以看到sparkling wine的銷量(第二排)與其他葡萄酒的銷量有相當的相關性。所以在模型中包含這些變量可能是一個好主意。

本文將介紹可以通過一種稱為自回歸分布滯后(ARDL)的方法來做到這一點。

Auto-Regressive Distributed Lag

ARDL模型采用自回歸。自回歸是大多數單變量時間序列模型的基礎。它主要分為兩個步驟。

首先將(單變量)時間序列從一個值序列轉換為一個矩陣?？梢杂糜醚訒r嵌入法（time delay embedding）來做到這一點。盡管名字很花哨，但這種方法非常簡單。它基于之前的最近值對每個值進行建模。然后建立一個回歸模型。未來值表示目標變量。解釋變量是過去最近的值。

多元時間序列的思路與此類似，我們可以將其他變量的過去值添加到解釋變量中。這就是了被稱為自回歸分布式滯后方法。分布式滯后的意思指的是使用額外變量的滯后。

現在我們把他們進行整合，時間序列中一個變量的未來值取決于它自身的滯后值以及其他變量的滯后值。

代碼實現

多變量時間序列通常是指許多相關產品的銷售數據。我們這里以葡萄酒銷售時間序列為例。當然ARDL方法也適用于零售以外的其他領域。

轉換時間序列

首先使用下面的腳本轉換時間序列。

import pandas as pd
 
 # https://github.com/vcerqueira/blog/
 from src.tde import time_delay_embedding
 
 wine = pd.read_csv('data/wine_sales.csv', parse_dates=['date'])
 
 # setting date as index
 wine.set_index('date', inplace=True)
 
 # you can simulate some data with the following code
 # wine = pd.DataFrame(np.random.random((100, 6)),
 #           columns=['Fortified','Drywhite','Sweetwhite',
 #                     'Red','Rose','Sparkling'])
 
 # create data set with lagged features using time delay embedding
 wine_ds = []
 for col in wine:
  
  
 
 # concatenating all variables
 wine_df = pd.concat(wine_ds, axis=1).dropna()
 
 # defining target (Y) and explanatory variables (X)
 predictor_variables = wine_df.columns.str.contains('\(t\-')
 target_variables = wine_df.columns.str.contains('Sparkling\(t\+')
 
 X = wine_df.iloc[:, predictor_variables]
 Y = wine_df.iloc[:, target_variables]

將 time_delay_embedding 函數應用于時間序列中的每個變量（第 18-22 行）。第 23 行將結果與我們的數據集進行合并。

解釋變量 (X) 是每個變量在每個時間步長的最后 12 個已知值（第 29 行）。以下是它們如何查找滯后 t-1（為簡潔起見省略了其他滯后值）：

目標變量在第30行中定義。這指的是未來銷售的6個值:

建立模型

準備好數據之后，就可以構建模型了。使用隨機森林進行一個簡單的訓練和測試循環。

from sklearn.model_selection import train_test_split
 from sklearn.metrics import mean_absolute_error as mae
 from sklearn.ensemble import RandomForestRegressor
 
 
 # train/test split
 X_tr, X_ts, Y_tr, Y_ts = train_test_split(X, Y, test_size=0.3, shuffle=False)
 
 # fitting a RF model
 model = RandomForestRegressor()
 model.fit(X_tr, Y_tr)
 
 # getting forecasts for the test set
 preds = model.predict(X_ts)
 
 # computing MAE error
 print(mae(Y_ts, preds))
 # 288.13

擬合模型之后(第11行)，得到了測試集中的預測(第14行)。該模型的平均絕對誤差為288.13。

滯后參數的選擇

上面的基線使用每個變量的 12 個滯后作為解釋變量。這是在函數 time_delay_embedding 的參數 n_lags 中定義的。那么應該如何設置這個參數的值呢？

很難先驗地說應該包括多少值，因為 這取決于輸入數據和特定變量。

解決這個問題的一種簡單方法是使用特征選擇。從相當數量的值開始，然后根據重要性評分或預測性能來修改這個數字，或者直接使用GridSearch進行超參數的搜索。

我們這里將簡單的演示一個判斷的過程：根據隨機森林的重要性得分選擇前 10 個特征。

# getting importance scores from previous model
 importance_scores = pd.Series(dict(zip(X_tr.columns, model.feature_importances_)))
 
 # getting top 10 features
 top_10_features = importance_scores.sort_values(ascending=False)[:10]
 top_10_features_nm = top_10_features.index
 
 X_tr_top = X_tr[top_10_features_nm]
 X_ts_top = X_ts[top_10_features_nm]
 
 # re-fitting the model
 model_top_features = RandomForestRegressor()
 model_top_features.fit(X_tr_top, Y_tr)
 
 # getting forecasts for the test set
 preds_topf = model_top_features.predict(X_ts_top)
 
 # computing MAE error
 print(mae(Y_ts, preds_topf))
 # 274.36

前10個特征比原始預測顯示出更好的預測性能。以下是這些功能的重要性:

目標變量(Sparkling)的滯后是最重要的。但是其他變量的一些滯后也是相關的。

ARDL 的擴展

多個目標變量預測，目前為止，我們都在預測單個變量（sparkling wine）。如果我們想要同時預測幾個變量呢？

這種方法被稱為：向量自回歸 (VAR)

就像在 ARDL 中一樣，每個變量都是根據其滯后和其他變量的滯后建模的。當想要預測多個變量而不僅僅是一個變量時，將使用 VAR。

與全局預測模型的關系

值得注意的是，ARDL并不等同于全局預測模型（Global Forecasting Models）。

在ARDL的前提下，每個變量的信息被添加到解釋變量中。變量的數量通常很少，且大小相同。

全局預測模型匯集了許多時間序列的歷史觀測結果。模型通過這些所有觀察結果進行建模。每一個新的時間序列都是作為新的觀察結果加入到數據中。全局預測模型通常涉及多達數千個時間序列量級也很大。

總結

本文的主要內容如下：多變量時間序列包含兩個或多個變量；ARDL 方法可用于多變量時間序列的監督學習；使用特征選擇策略優化滯后數。如果要預測多個變量，可以使用 VAR 方法。

最后本文的數據集在這里：

https://rdrr.io/cran/Rssa/man/AustralianWine.html