簡介
自回歸
多變量時間序列包含兩個或多個變量，研究這些數據集的目的是預測一個或多個變量，參見下面的示例。
上圖是包含9個變量的多變量時間序列。這些是智能浮標捕捉到的海洋狀況。
大多數預測模型都是基于自回歸的。這相當于解決了一個監督學習回歸任務。該序列的未來值是目標變量。輸入的解釋變量是每個變量最近的過去值。
自回歸在一個主要假設下工作。最近的過去值包含了關于未來的足夠信息。但這可能不一定是真的。我們可以嘗試從最近的數據中提取更多的信息。例如，滾動匯總統計信息有助于描述最近的動態。
自動化特征工程
特征工程包括提取和生成解釋變量，這是任何數據科學項目的關鍵。特征的質量是模型性能的一個核心方面，所以數據科學家在這個過程中花費了大量的時間。
特性工程通常是一個特別的過程：數據科學家基于他們的領域知識和專業知識創建特性，如果該過程的能夠自動化化處理將會為我們節省很多的時間。讓我們看看如何在多元時間序列中做到這一點。
基線模型
讀取數據
我們將使用從智能浮標收集的多元時間序列作為本文的數據集 [1]。 這個浮標位于愛爾蘭海岸。 它捕獲了 9 個與海洋條件相關的變量。 其中包括海水溫度、波浪高度和海水流速等。 上面的圖 1 顯示了 2022 年第一個月的情況。
以下是使用 pandas 讀取這些數據的方法：
import pandas as pd
# skipping second row, setting time column as a datetime column
# dataset available here: https://github.com/vcerqueira/blog/tree/main/data
buoy = pd.read_csv('data/smart_buoy.csv',
skiprows=[1],
parse_dates=['time'])
# setting time as index
buoy.set_index('time', inplace=True)
# resampling to hourly data
buoy = buoy.resample('H').mean()
# simplifying column names
buoy.columns = [
'PeakP', 'PeakD', 'Upcross',
'SWH', 'SeaTemp', 'Hmax', 'THmax',
'MCurDir', 'MCurSpd'
]
這個數據集研究的目標是預測SWH(顯著波高)變量的未來值。這個變量常被用來量化海浪的高度。這個問題的一個用例是估計海浪發電的大小，因為這種能源是一種越來越受歡迎的替代不可再生能源。
自回歸模型
時間序列是多元的，所以可以使用ARDL(Auto-regressive distributed lags)方法來解決這個任務。我們在之前也介紹過則個方法。下面是這個方法的實現：
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_percentage_error as mape
from sklearn.multioutput import MultiOutputRegressor
from lightgbm import LGBMRegressor
# https://github.com/vcerqueira/blog/blob/main/src/tde.py
from src.tde import time_delay_embedding
target_var = 'SWH'
colnames = buoy.columns.tolist()
# create data set with lagged features using time delay embedding
buoy_ds = []
for col in buoy:
col_df = time_delay_embedding(buoy[col], n_lags=24, horizon=12)
buoy_ds.append(col_df)
# concatenating all variables
buoy_df = pd.concat(buoy_ds, axis=1).dropna()
# defining target (Y) and explanatory variables (X)
predictor_variables = buoy_df.columns.str.contains('\(t\-')
target_variables = buoy_df.columns.str.contains(f'{target_var}\(t\+')
X = buoy_df.iloc[:, predictor_variables]
Y = buoy_df.iloc[:, target_variables]
# train/test split
X_tr, X_ts, Y_tr, Y_ts = train_test_split(X, Y, test_size=0.3, shuffle=False)
# fitting a lgbm model without feature engineering
model_wo_fe = MultiOutputRegressor(LGBMRegressor())
model_wo_fe.fit(X_tr, Y_tr)
# getting forecasts for the test set
preds_wo_fe = model_wo_fe.predict(X_ts)
# computing the MAPE error
mape(Y_ts, preds_wo_fe)
# 0.238
首先將時間序列轉化為一個自回歸問題。這是通過函數time_delay_embedding完成的。預測的目標是預測未來12個SWH值(horiznotallow=12)。解釋變量是序列中每個變量的過去的24個值(n_lag =24)。
我們這里直接使用LightGBM對每個預測層位進行訓練。這種方法法是一種常用的多步超前預測方法。它在scikit-learn中也有實現，名為MultiOutputRegressor。
上面的代碼構建和測試一個自回歸模型。解釋變量只包括每個變量最近的過去值。結果的平均絕對百分比誤差為0.238。
我們把這個結果作為基類對比，讓我們看看是否可以通過特性工程來提高。
多元時間序列的特征工程
本文本將介紹兩種從多元時間序列中提取特征的方法:
- 單變量特征提取。計算各變量的滾動統計。例如，滾動平均可以用來消除虛假的觀測;
- 二元特征提取。計算變量對的滾動統計，以總結它們的相互作用。例如，兩個變量之間的滾動協方差。
單變量特征提取
我們可以總結每個變量最近的過去值。例如，計算滾動平均來總結最近的情況?；蛘邼L動差量來了解最近的分散程度。
import numpy as np
SUMMARY_STATS = {
'mean': np.mean,
'sdev': np.std,
}
univariate_features = {}
# for each column in the data
for col in colnames:
# get lags for that column
X_col = X.iloc[:, X.columns.str.startswith(col)]
# for each summary stat
for feat, func in SUMMARY_STATS.items():
# compute that stat along the rows
univariate_features[f'{col}_{feat}'] = X_col.apply(func, axis=1)
# concatenate features into a pd.DF
univariate_features_df = pd.concat(univariate_features, axis=1)
如果能需要添加更多的統計數據?？梢韵騍UMMARY_STATS字典添加函數來實現這一點。將這些函數放在一個字典中可以保持代碼整潔。
二元特征提取
單變量統計漏掉了不同變量之間潛在的相互作用。所以我們可以使用二元特征提取過程捕獲這些信息。
這個想法是為不同的變量對計算特征?？梢允褂枚y計總結了這些對的聯合動態。
有兩種方法可以做到這一點:
- 滾動二元統計。計算以變量對作為輸入的統計信息。例如，滾動協方差或滾動相關性滾動二元統計的例子包括協方差、相關性或相對熵。
- 滾動二元變換，然后單變量統計。這將一對變量轉換為一個變量，并對該變量進行統計。例如，計算元素相互關系，然后取其平均值。有許多二元轉換的方法。例如，百分比差異、相互關聯或成對變量之間的線性卷積。通過第一步操作后，用平均值或標準偏差等統計數據對這些轉換進行匯總。
下面是用于性完成這兩個過程的代碼：
import itertools
import pandas as pd
from scipy.spatial.distance import jensenshannon
from scipy import signal
from scipy.special import rel_entr
from src.feature_extraction import covariance, co_integration
BIVARIATE_STATS = {
'covariance': covariance,
'co_integration': co_integration,
'js_div': jensenshannon,
}
BIVARIATE_TRANSFORMATIONS = {
'corr': signal.correlate,
'conv': signal.convolve,
'rel_entr': rel_entr,
}
# get all pairs of variables
col_combs = list(itertools.combinations(colnames, 2))
bivariate_features = []
# for each row
for i, _ in X.iterrows():
# feature set in the i-th time-step
feature_set_i = {}
for col1, col2 in col_combs:
# features for pair of columns col1, col2
# getting the i-th instance for each column
x1 = X.loc[i, X.columns.str.startswith(col1)]
x2 = X.loc[i, X.columns.str.startswith(col2)]
# compute each summary stat
for feat, func in BIVARIATE_SUMMARY_STATS.items():
feature_set_i[f'{col1}|{col2}_{feat}'] = func(x1, x2)
# for each transformation
for trans_f, t_func in BIVARIATE_TRANSFORMATIONS.items():
# apply transformation
xt = t_func(x1, x2)
# compute summary stat
for feat, s_func in SUMMARY_STATS.items():
feature_set_i[f'{col1}|{col2}_{trans_f}_{feat}'] = s_func(xt)
bivariate_features.append(feature_set_i)
bivariate_features_df = pd.DataFrame(bivariate_features, index=X.index)
字典bivariate_transforms或BIVARIATE_STATS中添加其他的函數，可以添加額外的轉換或統計信息。
在提取所有特征之后，我們將將它們連接到原始解釋變量。訓練和測試的過程和之前的是一樣的，只不過我們增加了一些人工生成的變量。
# concatenating all features with lags
X_with_features = pd.concat([X, univariate_features_df, bivariate_features_df], axis=1)
# train/test split
X_tr, X_ts, Y_tr, Y_ts = train_test_split(X_with_features, Y, test_size=0.3, shuffle=False)
# fitting a lgbm model with feature engineering
model_w_fe = MultiOutputRegressor(LGBMRegressor())
model_w_fe.fit(X_tr, Y_tr)
# getting forecasts for the test set
preds_w_fe = model_w_fe.predict(X_ts)
# computing MAPE error
print(mape(Y_ts, preds_w_fe))
# 0.227
得到了0.227的平均絕對百分比誤差，這是一個小小的提高，因為我們的基線是0.238。
特征選擇
以上提取過程共得到了558個解釋變量。根據變量和匯總統計信息的數量，這可能會產生高維問題。因此，從數據集中刪除糟糕或冗余的特征是很重要的。
我們將找到一些重要特征并重新訓練
# getting the importance of each feature in each horizon
avg_imp = pd.DataFrame([x.feature_importances_
# getting the top 100 features
n_top_features = 100
importance_scores = pd.Series(dict(zip(X_tr.columns, avg_imp)))
top_features = importance_scores.sort_values(ascending=False)[:n_top_features]
top_features_nm = top_features.index
# subsetting training and testing sets by those features
X_tr_top = X_tr[top_features_nm]
X_ts_top = X_ts[top_features_nm]
# re-fitting the lgbm model
model_top_features = MultiOutputRegressor(LGBMRegressor())
model_top_features.fit(X_tr_top, Y_tr)
# getting forecasts for the test set
preds_top_feats = model_top_features.predict(X_ts_top)
# computing MAE error
mape(Y_ts, preds_top_feats)
# 0.229
可以看到前100個特性與完整的558個特性的性能相似。以下是前15個特征的重要性(為了簡潔起見省略了其他特征):
可以看到最重要的特征是目標變量的第一個滯后值。一些提取的特征也出現在前15名中。例如第三個特征SWH|Hmax_js_div。這表示目標變量的滯后與Hmax的滯后之間的Jensen-Shannon散度。第五個特性是SeaTemp_sdev，表示海洋溫度的標準偏差滯后。
另一種去除冗余特征的方法是應用相關性過濾器。刪除高度相關的特征以減少數據的維數，這里我們就不進行演示了。
總結
本文側重于多變量時間序列的預測問題。特征提取過程應用于時間序列的多個子序列，在每個時間步驟中，都要用一組統計數據總結過去24小時的數據。
我們也可以用這些統計來一次性描述整個時間序列。如果我們目標是將一組時間序列聚類，那么這可能是很有用。用特征提取總結每個時間序列。然后對得到的特征應用聚類算法。
用幾句話總結本文的關鍵點：
- 多變量時間序列預測通常是一個自回歸過程
- 特征工程是數據科學項目中的一個關鍵步驟。
- 可以用特征工程改進多元時間序列數據。這包括計算單變量和雙變量轉換和匯總統計信息。
- 提取過多的特征會導致高維問題?？梢允褂锰卣鬟x擇方法來刪除不需要的特征。
