時間是決定企業(yè)興衰的最關(guān)鍵因素。這就是為什么我們看到商店和電子商務(wù)平臺的銷售與節(jié)日一致。這些企業(yè)分析多年的消費數(shù)據(jù)，以了解打開大門的最佳時間，并看到消費支出的增加。
但是，作為一個數(shù)據(jù)科學(xué)家，你怎么能進行這種分析呢？別擔(dān)心，你不需要建造一臺時間機器！時間序列建模是一種強大的技術(shù)，是理解與預(yù)測趨勢和模式的門戶。
但是即使是時間序列模型也有不同的方面。我們在Web上看到的大多數(shù)例子都是用單變量時間序列來處理的。不幸的是，現(xiàn)實世界的用例并不是這樣工作的。有多個變量在起作用，同時處理所有這些變量是數(shù)據(jù)科學(xué)家體現(xiàn)其價值的地方。
在這篇文章中，我們將了解什么是多元時間序列，以及如何處理它。我們還將在Python中進行案例研究并實現(xiàn)它，以使您對該主題有實際的理解。
目錄表
1.單變量與多變量時間序列
— —單變量時間序列
— —多元時間序列
2.多變量時間序列-向量自回歸（VAR）的處理
3.為什么我們需要VAR？
4.多變量時間序列的平穩(wěn)性
5.訓(xùn)練驗證拆分
6.Python實現(xiàn)
1.單變量與多變量時間序列
本文假定讀者熟悉單變量時間序列的性質(zhì)，以及用于預(yù)測的各種技術(shù)。由于本文將重點討論多變量時間序列，因此我建議您閱讀以下文章，這些文章可以作為對單變量時間序列的良好介紹：
但是，在討論多元時間序列的細節(jié)之前，我將快速地介紹一下什么是單變量時間序列。讓我們逐一觀察它們的差異。
1.1單變量時間序列
一個單變量時間序列，顧名思義，是一個具有單一時間依賴變量的序列。
例如，在過去的2年中，查看下面的樣本數(shù)據(jù)集，該數(shù)據(jù)集包括溫度值（每小時）。這里，溫度是因變量（取決于時間）。
如果要求我們預(yù)測未來幾天的溫度，我們將查看過去的值，并嘗試測量和提取一個模式。我們會注意到早晨和晚上的溫度較低，下午則達到峰值。此外，如果你有過去幾年的數(shù)據(jù)，你會發(fā)現(xiàn)，在11月到1月份，天氣比較冷，而在4月至6月比較熱。
這樣的觀察將有助于我們預(yù)測未來的價值觀。你注意到我們只使用了一個變量（過去2年的溫度）嗎？因此，這被稱為單變量時間序列分析/預(yù)測。
1.2多元時間序列（MTS）
多變量時間序列具有一個以上的時間依賴變量。每個變量不僅取決于其過去的值，而且還對其他變量有一定的依賴性。這種依賴性用于預(yù)測未來的價值。聽起來很復(fù)雜？讓我解釋一下。
考慮上面的例子。現(xiàn)在假設(shè)我們的數(shù)據(jù)集包括過去兩年的汗水百分比、露點、風(fēng)速、云層覆蓋率等，以及溫度值。在這種情況下，有多個變量被認為是最佳預(yù)測溫度。像這樣的系列將屬于多元時間序列的范疇。下面是一個例證：
現(xiàn)在我們了解了多元時間序列的樣子，讓我們了解如何利用它來建立預(yù)測。
2.多元時間序列的處理——VAR
在本節(jié)中，我將介紹多變量時間序列預(yù)測中最常用的方法之一——向量自回歸（VAR）。
在VAR模型中，每個變量是其自身過去值和所有其他變量的過去值的線性函數(shù)。為了更好地解釋這一點，我將使用一個簡單的視覺例子：
我們有兩個變量y1和y2。我們需要從過去N值的給定數(shù)據(jù)預(yù)測時間t上的這兩個變量的值。為了簡單起見，我已經(jīng)考慮滯后值為1。
為了計算y1（t），我們將使用y1和y2的過去值。類似地，為了計算y2（t），將使用y1和y2的過去值。下面是表示這種關(guān)系的簡單數(shù)學(xué)方法：
在這里，a1和a2是常數(shù)項，w11、w12、w21和w22是系數(shù)，e1和e2是誤差項。
這些方程類似于AR方程的過程。由于AR過程用于單變量時間序列數(shù)據(jù)，所以未來值只是它們自己的過去值的線性組合?？紤]AR（1）過程：
y(t) = a + w*y(t-1) +e
在這種情況下，我們只有一個變量—y、常數(shù)項—a、誤差項—e和系數(shù)—w。為了適應(yīng)VAR的每個方程中的多變量項，我們將使用向量。我們可以用下列形式寫出方程式（1）和（2）：
這兩個變量是Y1和Y2，其次是常數(shù)、系數(shù)度量、滯后值和誤差度量。這是VAR（1）過程的向量方程。對于VAR（2）過程，將另一個時間向量項（T-2）加入到方程中，以推廣P滯后：
上述方程表示具有變量y1、y2…yk的var（p）過程。同樣可以寫成：
方程中εt表示多元矢量白噪聲。對于一個多變量時間序列，εt應(yīng)該是滿足下列條件的連續(xù)隨機向量：
1.E（εt）＝0
誤差向量的期望值為0。
2.E（εt1，εt2′）＝12
εt和εt′的期望值是級數(shù)的標準偏差。
3.我們?yōu)槭裁葱枰猇AR？
回想一下我們早些時候看到的溫度預(yù)測例子，可以將其作為一個多變量單變量序列進行討論。我們可以使用簡單的單變量預(yù)測方法如AR來解決這個問題，因為其目的是預(yù)測溫度，所以我們可以簡單地去除其他變量（除了溫度），并在剩余的單變量序列上擬合一個模型。
另一個簡單的想法是使用我們已經(jīng)知道的技術(shù)來逐個預(yù)測每個系列的值。這將使工作非常直截了當！那你為什么還要學(xué)習(xí)另一種預(yù)測技術(shù)呢？這個話題已經(jīng)足夠復(fù)雜了嗎？
從上面的方程（1）和（2）可以清楚地看到，每個變量都使用每個變量的過去值來進行預(yù)測。與AR不同，VAR能夠理解和使用多個變量之間的關(guān)系。這有助于描述數(shù)據(jù)的動態(tài)行為，并提供更好的預(yù)測結(jié)果。此外，實現(xiàn)VAR與使用任何其他單變量技術(shù)一樣簡單（您將在最后一節(jié)中看到）。
4.多變量時間序列的平穩(wěn)性
通過研究單變量概念，我們知道，平穩(wěn)時間序列往往會給我們一組更好的預(yù)測。如果您不熟悉平穩(wěn)性的概念，請先閱讀本文：處理非平穩(wěn)時間序列的溫和介紹。
總之，對于給定的單變量時間序列：
y(t) = c*y(t-1) + ε t
如果c＜1，則該級數(shù)是固定的?，F(xiàn)在，回憶一下我們的VAR過程方程：
對于一系列不動點，π（L）- 1π的特征值應(yīng)小于1。考慮到派生變量的數(shù)量，這可能看起來很復(fù)雜。這個想法已經(jīng)在視頻中用一個簡單的數(shù)值例子來解釋。我高度鼓勵看它鞏固你的理解。
類似于單變量序列的增強Dickey-Fuller檢驗，我們有用來檢驗任何多變量時間序列數(shù)據(jù)的平穩(wěn)性的Johansen檢驗。我們將在本文的最后一節(jié)中看到如何進行測試。
5.訓(xùn)練驗證拆分
如果你以前使用過單變量時間序列數(shù)據(jù)，你會知道訓(xùn)練驗證集。創(chuàng)建驗證集的想法是在使用模型進行預(yù)測之前，分析模型的性能。
創(chuàng)建時間序列問題的驗證集是棘手的，因為我們必須考慮到時間成分。不能直接使用train_test_split 分割或K-折疊驗證，因為這會破壞序列中的模式。應(yīng)考慮日期和時間值創(chuàng)建驗證集。
假設(shè)我們必須使用過去兩年的數(shù)據(jù)來預(yù)測未來兩個月的溫帶氣候、露點、云百分比等。一種可能的方法是把過去兩個月的數(shù)據(jù)放在一邊，并在剩下的22個月內(nèi)訓(xùn)練模型。
一旦模型被訓(xùn)練，我們就可以使用它來對驗證集進行預(yù)測?；谶@些預(yù)測和實際值，我們可以檢查模型執(zhí)行得有多好，以及模型沒有很好地執(zhí)行的變量。為了進行最終的預(yù)測，使用完整的數(shù)據(jù)集（訓(xùn)練和驗證集相結(jié)合）。
6.Python實現(xiàn)
在本節(jié)中，我們將在玩具數(shù)據(jù)集上實現(xiàn)向量AR模型。我已經(jīng)使用了空氣質(zhì)量數(shù)據(jù)集，你可以從這里下載。
#import required packages
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
#read the data
df = pd.read_csv("AirQualityUCI.csv", parse_dates=[['Date', 'Time']])
#check the dtypes
df.dtypes
Date_Time? ? ? ? object
CO(GT)? ? ? ? ? ? int64
PT08.S1(CO)? ? ? int64
NMHC(GT)? ? ? ? ? int64
C6H6(GT)? ? ? ? ? int64
PT08.S2(NMHC)? ? int64
NOx(GT)? ? ? ? ? int64
PT08.S3(NOx)? ? ? int64
NO2(GT)? ? ? ? ? int64
PT08.S4(NO2)? ? ? int64
PT08.S5(O3)? ? ? int64
T? ? ? ? ? ? ? ? int64
RH? ? ? ? ? ? ? ? int64
AH? ? ? ? ? ? ? ? int64
dtype: object
Date_Time列的數(shù)據(jù)類型是Object，我們需要將其更改為datetime。此外，為了準備數(shù)據(jù)，我們需要索引具有datetime。遵循以下命令：
df['Date_Time'] = pd.to_datetime(df.Date_Time , format = '%d/%m/%Y %H.%M.%S')
data = df.drop(['Date_Time'], axis=1)
data.index = df.Date_Time
下一步是處理缺失的值。由于數(shù)據(jù)中缺失的值被替換為值-200，我們將不得不用更好的數(shù)字來計算缺失的值。考慮這個——如果當前露點值丟失，我們可以安全地假設(shè)它將接近前一個小時的值。有道理，對吧？在這里，我將用前一個值來填充-200。
您可以選擇使用前面幾個值的平均值來替換該值，或者同時使用前一天的值（您可以在下面的評論區(qū)中共享您關(guān)于輸入缺失值的想法）。
#missing value treatment
cols = data.columns
for j in cols:
? ? for i in range(0,len(data)):
? ? ? if data[j][i] == -200:
? ? ? ? ? data[j][i] = data[j][i-1]
#checking stationarity
from statsmodels.tsa.vector_ar.vecm import coint_johansen
#since the test works for only 12 variables, I have randomly dropped
#in the next iteration, I would drop another and check the eigenvalues
johan_test_temp = data.drop([ 'CO(GT)'], axis=1)
coint_johansen(johan_test_temp,-1,1).eig
下面是測試的結(jié)果：
array([ 0.17806667,? 0.1552133 ,? 0.1274826 ,? 0.12277888,? 0.09554265,
? ? ? ? 0.08383711,? 0.07246919,? 0.06337852,? 0.04051374,? 0.02652395,
? ? ? ? 0.01467492,? 0.00051835])
現(xiàn)在我們可以繼續(xù)創(chuàng)建驗證集來適應(yīng)模型，并測試模型的性能：
#creating the train and validation set
train = data[:int(0.8*(len(data)))]
valid = data[int(0.8*(len(data))):]
#fit the model
from statsmodels.tsa.vector_ar.var_model import VAR
model = VAR(endog=train)
model_fit = model.fit()
# make prediction on validation
prediction = model_fit.forecast(model_fit.y, steps=len(valid))
預(yù)測是數(shù)組的形式，其中每個列表代表行的預(yù)測。我們將把這個轉(zhuǎn)變成一個更具代表性的格式。
#converting predictions to dataframe
pred = pd.DataFrame(index=range(0,len(prediction)),columns=[cols])
for j in range(0,13):
? ? for i in range(0, len(prediction)):
? ? ? pred.iloc[i][j] = prediction[i][j]
#check rmse
for i in cols:
? ? print('rmse value for', i, 'is : ', sqrt(mean_squared_error(pred[i], valid[i])))
以上代碼的輸出：
rmse value for CO(GT) is :? 1.4200393103392812
rmse value for PT08.S1(CO) is :? 303.3909208229375
rmse value for NMHC(GT) is :? 204.0662895081472
rmse value for C6H6(GT) is :? 28.153391799471244
rmse value for PT08.S2(NMHC) is :? 6.538063846286176
rmse value for NOx(GT) is :? 265.04913993413805
rmse value for PT08.S3(NOx) is :? 250.7673347152554
rmse value for NO2(GT) is :? 238.92642219826683
rmse value for PT08.S4(NO2) is :? 247.50612831072633
rmse value for PT08.S5(O3) is :? 392.3129907890131
rmse value for T is :? 383.1344361254454
rmse value for RH is :? 506.5847387424092
rmse value for AH is :? 8.139735443605728
在驗證集的測試之后，讓模型在完整的數(shù)據(jù)集上擬合
#make final predictions
model = VAR(endog=data)
model_fit = model.fit()
yhat = model_fit.forecast(model_fit.y, steps=1)
print(yhat)
小結(jié)
在我開始這篇文章之前，使用多變量時間序列的想法似乎在其范圍內(nèi)令人望而生畏。這是一個復(fù)雜的話題，所以要花時間了解細節(jié)。最好的學(xué)習(xí)方法是練習(xí)，所以我希望上面的Python實現(xiàn)會對你有用。
我建議您在您選擇的數(shù)據(jù)集上使用此方法。這將進一步鞏固你對這個復(fù)雜而又非常有用的話題的理解。免費視頻教程：www.mlxs.top
