不會時間序列預測?不要緊,大神來教你

點擊上方「Python數據之道」，選擇「星標公眾號」本文內容較長，較為詳細的闡述了進行時間序列預測的步驟，有些內容可能暫時用不到或者看不懂，但不要緊，知道有這麼一個概念，後續碰到的時候，繼續深入學習以及使用就可以。一文弄懂時間序列預測的基本原理Photo by Adrian Schwarz on Unsplash我們被隨處可見的模式所包圍，人們可以注意到四季與天氣的關係模式，以交通量計算的交通高峰期的模式，你的心跳或者是股票市場和某些產品的銷售周期。分析時間序列數據對於發現這些模式和預測未來非常有用。有幾種方法可以創建這類預測，在本文中，我將介紹最基本且最傳統的方法概念。所有代碼都是用 Python 編寫的，並且在 GitHub 上可以看到所有的信息。https://nbviewer.jupyter.org/github/leandrovrabelo/tsmodels/blob/master/notebooks/english/Basic Principles for Time Series Forecasting.ipynb1平穩序列平穩時間序列是指統計特性，如均值、方差和自相關係數，隨時間相對恆定的序列。因此，非平穩序列是統計特性隨時間變化的序列。在開始任何預測建模之前，都有必要驗證這些統計屬性是否是常量，我將一一解釋下面的每個點：常數均值一個平穩序列在時間上具有一個相對穩定的均值，這個值沒有減少或者增加的趨勢。圍繞常數均值的小的變化，使我們更容易推測未來。在某些情況下，相對於平均值的變量比較小，使用它可以很好地預測未來。下圖顯示了變量與該常數平均值相對於時間變化的關係:在這種情況下，如果序列不是平穩的，對未來的預測將是無效的，因為平均值周圍的變量會顯著偏離，如下圖所示:在上圖中，我們可以明顯看到上升的趨勢，均值正在逐漸上升。在這種情況下，如果使用均值進行未來值的預測，誤差將非常大，因為預測價格會總是低於實際價格。常數方差當序列的方差為常數時，我們知道均值和標準差之間存在一種關係。當方差不為常數時（如下圖所示），預測在某些時期可能會有較大的誤差，而這些時期是不可預測的。可以預測到，隨著時間的推移直到未來，方差會保持不穩定。為了減小方差效應，可以採用對數變換。在本例中，也可以使用指數變換，如 Box-Cox 方法，或者使用膨脹率調整。自相關序列當兩個變量在時間上的標準差有相似的變化時，你可以說這些變量是相關的。例如，體重會隨著心臟疾病而增加，體重越大，心臟問題的發生率就越大。在這種情況下，相關性是正的，圖形應該是這樣的:負相關的情況類似於這樣：對工作安全措施的投入越多，工作相關的事故數量就越少。當談到自相關時，意思是某些先前時期與當前時期存在相關性，這種相關性是滯後的。例如，在以小時為單位的測量值序列中，今天 12:00 的溫度與 24 小時前的 12:00 的溫度非常相似。如果你比較 24 小時內的溫度變化，就會存在自相關，在本例中，我們將與第 24 小時前的時間存在自相關關係。自相關是使用單個變量創建預測的一種情況，因為如果沒有相關性，就不能使用過去的值來預測未來；當有多個變量時，則可以驗證因變量和獨立變量的滯後之間是否存在相關性。如果一個序列不存在自相關關係，那麼它就是隨機且不可預測的，做預測的最佳方法通常是使用前一天的值。我將在下面使用更詳細的圖表來解釋。從這裡開始我將分析 Esalq 上的每周含水乙醇價格（這是巴西談判含水乙醇的價格參考），數據可以在這裡(https://www.cepea.esalq.usp.br/br/indicador/etanol.aspx)下載。在開始任何分析之前，要將數據劃分為訓練集和測試集。劃分訓練集和測試集數據當我們要創建時序預測模型時，將數據劃分為兩部分至關重要：訓練集：這些數據將是定義模型係數/參數的主要依據；測試集：這些數據將被分離且對模型不可見，用於測試模型是否有效（通常將這些值與模型結果進行比較，最後測量平均誤差）。測試集的大小通常約為總樣本的20%，儘管這個百分比取決於你擁有的樣本大小以及你希望提前多少時間進行預測。理想情況下，測試集應至少與所需預測的最大範圍相同。與其他如分類和回歸等不受時間影響的預測方法不同，在時間序列中，不可以將訓練和測試數據從數據中隨機抽樣取出，我們必須遵循序列的時間標準，訓練數據應該始終是在測試數據之前。在本例中，我們有Esalq 含水乙醇的 856 周的價格數據，使用前 700 周的數據作為訓練集，後 156 周（3年，18%）的數據用作測試集：從現在開始，我們只使用訓練集來做研究，測試集僅用於驗證我們的預測。每一個時間序列可以分為三個部分：趨勢、季節性和殘差，殘差是將前兩部分從序列中去除後剩下的部分，使用這種分割方法之後：顯然，該序列具有上升趨勢，在每一年的年底到年初之間達到峰值，在4月和9月之間達到最低值（此時在巴西中南部開始甘蔗的壓榨）。我們仍然建議使用統計測試來確認序列是否是平穩的，這裡將使用兩個測試：Dickey-Fuller 測試和 KPSS 測試。首先，我們將使用 Dickey-Fuller 檢驗，我將使用 5% 的基礎 P 值，也就是說，如果 P 值低於 5% 這意味著這個序列在統計上是平穩的。此外，還有模型的統計檢驗，可以將檢驗值與 1%、5%、10% 的臨界值進行比較，如果統計檢驗低於選定的某個臨界值，就認為序列是平穩的:在本例中，Dickey-Fuller 檢驗結果表明序列不是平穩的（ P 值 36%，臨界值 5% 小於統計檢驗）。現在我們要用 KPSS 檢驗分析序列，與 Dickey-Fuller 檢驗不同，KPSS 檢驗已經假設序列是平穩的，只有當 P 值小於 5% 或統計檢驗小於某個臨界值時，序列才不是平穩的：KPSS 檢驗證實了 Dickey-Fuller 檢驗的正確性，同時也表明該序列不是平穩的，因為 P 值為 1%，統計檢驗高於任何臨界值。2
將序列轉換為平穩狀態差分差分法用來移除趨勢信號，也可以用來減少方差，它只是 T 周期的值與前一個 T-1 周期值的差值。為了更容易理解，下面我們只用一小部分的乙醇價格，以便更好地可視化，可以看到從 2005 年 5 月開始價格上漲，直到 2006 年 5 月中旬，價格每周都在上漲，這就累積了一個上升的趨勢，這種情況下，屬於非平穩序列。當進行一階微分時（如下圖），我們去除了序列的累積效應，並且僅顯示了整個系列中時段 T 相對於時段 T-1 的變化，因此如果 3 天前的價格為 800 BRL 且已漲到 850.00 BRL，差價將是 50.00 BRL，如果今天的價格是 860.00 BRL，那麼差價將是 - 10.00 BRL。通常只需要一階微分就足夠將序列轉換為平穩狀態，但如果需要，可以應用二階微分，在這種情況下，將對一階微分的值進行求導（幾乎沒有二階以上微分的情況）。同樣的例子，要進行二次微分，我們必須取 T 時刻減去 T-1 時刻的微分：2.9 BRL -5.5 BRL = - 2.6 BRL 等等。我們來做一下 Dickey-Fuller 測試，看看這個序列是否會在一階微分後是平穩的：在這種情況下，我們確定該序列是平穩的，P 值為零，並且當我們比較統計檢驗的值時，它遠遠低於臨界值。在下一個例子中，我們將嘗試調整通貨膨脹率將一個序列轉換到平穩狀態。膨脹率調整價格是相對於交易時間的，2002 年乙醇的價格是 680.00 BRL，如果現在產品的價格是這個價格，很多工廠肯定會倒閉，因為這個價格非常低。為了讓序列平穩，我將基於當前值使用巴西 IPCA 索引（巴西的 CPI 指數）調整整個序列，從訓練區間的結尾（2016年4月）到研究的開始，數據的來源是 IBGE 網站。如圖所示，上升趨勢已經消失，只剩下季節性振蕩， Dickey-Fuller 測試也證實了這個序列現在是穩定的。如果好奇，可以參閱下面的圖表，其中調整後的價格與原始系列的通貨膨脹率相對應。減小方差對數變換對數變換通常用於將指數增長的序列轉換為具有更趨於線性增長的序列，在本例中，我們將使用自然對數（Natural Logarithm，NL），其底數為 2.718 ，這種對數類型在經濟模型中被廣泛使用。轉換成 NL 值的差值近似等於原始序列值的百分比變化，作為降低不同價格序列的方差的基礎是很有效的，如下例:如果我們有一個產品在 2000 年價格上漲，從 50.00 BRL 到 52.50 BRL，幾年後（2019年），價格已經是 100.00 BRL，已經上漲到 105.00 BRL，價格之間的絕對差分別是 2.50 BRL 和 5.00 BRL，但兩者的百分比差為 5% 。當我們對這些價格中使用 NL 時，我們得到：NL (52.50) - NL(50.00) = 3.96 - 3.912 = 0.048 或 4.8%，同樣地，在第二個價格序列中使用 LN 時，我們得到：NL (105) - NL(100) = 4.654-4.605 = 0.049 或 4.9% 。在這個例子中，我們可以通過把幾乎所有的東西都放到相同的基上來減少差異值。

price1 = np.log(52.5) - np.log(50)

price2 = np.log(105) - np.log(100)

printf('The percentage variation of the first example is {round(price1*100,1)} and the second is {round(price2*100,1)}')

Box-Cox 變換（指數變換）Box-Cox 轉換也是一種轉換序列的方法，lambda（λ）的值是用於轉換序列的參數。簡而言之，這個函數是幾個指數變換函數的結合，我們需要找到轉換序列的 lambda 的最佳值，使其分布更接近正態高斯分布。使用此轉換的一個條件是序列只有正值，公式為：接下來我將繪製原始序列及其分布圖，然後用 lambda 最佳值繪製新的轉換序列及其分布圖，為了找到 lambda 的值，我們將使用庫 Scipy 的 boxcox 函數生成轉換的序列和理想 lambda 值：下面是一個交互式圖表，在圖中可以更改 lambda 值和檢查更改：此工具通常用於提高模型的性能，因為它使模型更趨於正態分布，記住在完成模型的預測後，必須根據以下公式反轉轉換到原始的基數：尋找相關時滯為了便於預測，具有單一變量的序列必須具有自相關性，即，當前時段必須是能夠基於較早的時段（滯後）而解釋的。由於這個序列每周為一周期，1 年大約 52 周，我將使用 60 的滯後期的自相關函數來驗證當前周期與這些滯後的相關性。通過對上述自相關圖的分析，似乎所有的滯後都可以用來為未來事件創建預測，因為它們的正相關接近 1 ，而且都在置信區間之外，但這一特徵屬於非平穩序列。另一個非常重要的函數是部分自相關函數，其中消除了先前的滯後對當前區間的影響，只保留了當前區間滯後的影響來分析，例如：第四個滯後的偏自相關將消除第一、第二和第三個滯後的影響。可以看到，幾乎沒有滯後對當前周期有影響，但是正如前面所演示的，沒有微分的序列不是平穩的，我們現在用一階微分的序列繪製這兩個函數來展示原理：自相關曲線變化顯著，表明該序列僅在第一個滯後期具有顯著相關，在第 26 個滯後(半年)左右具有負相關的季節效應。為了做出預測，我們必須注意一個找到相關的滯後現象的非常重要的細節，重要的是這種關聯背後的原因，因為如果沒有邏輯上的原因，就有可能是偶然的，當包含更多的數據時，這種關聯就會消失。另一個重點是自相關和部分自相關圖對異常值非常敏感，因此分析時間序列本身並與兩個自相關圖進行對比非常重要。在這個例子中，第一個滯後與當前周期具有高度相關性，因為前一周的價格歷史上沒有顯著變化，在相同的情況下，第 26 個滯後呈現負相關，表明與當前時期相反的趨勢，可能原因是一年內不同時期供需不同。隨著膨脹率調整後的序列已經趨於平穩，我們將使用它來創建我們的預測，下圖是調整後序列的自相關和部分自相關圖：想要了解更多信息的話，杜克大學教授 Robert Nau 的網站是與此主題相關的最佳網站之一。（http://people.duke.edu/~rnau/411home.htm)3 模型評價指標為了分析預測值是否接近當前值，必須對誤差進行測量，此種情況下的誤差（或殘差）基本上是 Yreal-YpredYreal-Ypred （這個暫時不知道怎麼翻譯， real 真實值， pred 預測值）。對訓練數據中的錯誤進行評估以驗證模型是否具有良好的確定性，然後通過檢查測試數據中的誤差（模型未「看到」的數據）來驗證模型。當將訓練數據與測試數據進行對比時，檢查誤差對於驗證你的模型是否過擬合或欠擬合非常重要。平均預測誤差——偏差（bias）它只是被評估序列的平均誤差，值可以是正的也可以是負的。該指標表明，模型傾向於預測實際值以上（負誤差）還是實際值以下（正誤差），因此也可以說平均預測誤差是模型的偏差。MAE——平均絕對誤差這個指標與上面提到的預測的平均誤差非常相似，唯一的區別是將誤差的負值轉化為正值，然後計算平均值。這個指標在時間序列中被廣泛使用，因為在一些情況下，負誤差可以抵消正誤差，使人誤以為模型是準確的，而在用 MAE 的情況下不會發生，因為這個指標顯示預測距離實際值有多遠，不管數值大還是小，示例如下：

a = np.array([1,2,3,4,5])

b = np.array([5,4,3,2,1])

error = a - b

MFE = error.mean()

MAE = np.abs(error).mean()

print(f'The error of each model value looks like this: {error}')

print(f'The MFE error was {MFE}, the MAE error was {MAE}')

與 MAE 和 MFE 不同，MSE 值是平方單位，而不是模型單位。RMSE——均方根誤差這個指標只是 MSE 的平方根，使誤差返回到模型的度量單位（BRL/m3），因為它對時間序列在平方過程中產生的較大誤差更為敏感而非常有用。MAPE——平均絕對百分誤差這是另一個可用的有趣的指標，它通常在管理報告中使用，因為誤差是以百分比度量的，所以產品 X 的錯誤可以與產品 Y 的誤差進行比較。該指標的計算取誤差的絕對值除以當前價格，然後計算平均值：我們來創建一個函數，用幾個評估指標來評估訓練數據和測試數據的誤差:

#Libraries to create the function:

from math import sqrt

from sklearn.metrics import mean_squared_error,mean_absolute_error, mean_absolute_error

def check_error(orig, pred, name_col='', index_name=''):

bias = np.mean(orig - pred)

mse = mean_squared_error(orig, pred)

rmse = sqrt(mean_squared_error(orig, pred))

mae = mean_absolute_error(orig, pred)

mape = np.mean(np.abs((orig - pred) / orig)) * 100

error_group = [bias, mse, rmse, mae, mape]

serie = pd.DataFrame(error_group, index=['BIAS','MSE','RMSE','MAE', 'MAPE'], columns=[name_col])

serie.index.name = index_name

return serie

分析這個圖是非常重要的，因為在這個圖中我們可以檢查模式，它可以告訴我們是否需要對模型進行一些修改，理想的情況是誤差沿著預測序列線性分布。https://en.wikipedia.org/wiki/Q–Q_plot總的來說這是一個顯示了殘差在理論上應該如何分布的圖形，遵循高斯分布，而不是實際情況。

def plot_error(data, figsize=(18,8)):

# Creating the column error

data['Error'] = data.iloc[:,0] -data.iloc[:,1]

plt.figure(figsize=figsize)

ax1 = plt.subplot2grid((2,2), (0,0))

ax2 = plt.subplot2grid((2,2), (0,1))

ax3 = plt.subplot2grid((2,2), (1,0))

ax4 = plt.subplot2grid((2,2), (1,1))

#Plotting actual and predicted values

ax1.plot(data.iloc[:,0:2])

ax1.legend(['Real','Pred'])

ax1.set_title('Real Value vs Prediction')

# Error vs Predicted value

ax2.scatter(data.iloc[:,1], data.iloc[:,2])

ax2.set_xlabel('Predicted Values')

ax2.set_ylabel('Residual')

ax2.set_title('Residual vs Predicted Values')

## Residual QQ Plot

sm.graphics.qqplot(data.iloc[:,2], line='r', ax=ax3)

# Autocorrelation Plot of residual

plot_acf(data.iloc[:,2], lags=60, zero=False, ax=ax4)

plt.tight_layout()

plt.show()

許多人還使用這種方法作為 基線（baseline），試圖用更複雜的模型來改進。QQ 圖顯示了有一些比理論上要大些（包括正負值）的殘差，這些是所謂的異常值，但在第一，第六和第七個滯後中仍然存在明顯的自相關，這可以用於改進模型。同樣地，我們現在將在測試數據中進行預測。預測序列的第一個值將是訓練數據的最後一個值，然後這些值將按照測試的當前值逐步更新，依此類推：RMSE 和 MAE 的誤差與訓練數據相似，QQ 圖與殘差更符合理論值，可能是由於與訓練數據相比樣本值較少。在對比殘差與預測值的圖表中，我們注意到當價格上漲時，誤差絕對值有增加的趨勢，可能用對數調整會減少誤差的擴大並完成殘差相關圖，表明由於第一個滯後有很強的相關性，因此仍有改進的空間，可能添加基於第一個滯後的回歸來改進預測。下一個模型是簡單平均值：簡單平均另一種預測方法是使用序列平均值，通常當數值在平均值附近振蕩時，具有常數的方差，沒有上升或下降趨勢時，這種預測形式是好的，但是也能使用更好的方法，其中可以使用季節模式進行預測。此模型使用數據首端直到分析的前一個時期的平均值，並且按天擴展到數據結束，最後，趨勢是一條直線，我們現在將此模型與第一個模型的誤差進行比較：在測試數據中，我將繼續使用訓練數據一開始的均值，並展開添加到測試數據上：簡單均值模型無法捕獲序列的相關信息，如真實值和預測值圖中所示，也可以在相關性和殘差和預測圖中看到。簡單滑動平均：滑動平均是針對給定周期（例如 5 天）計算的平均值，它是滑動的並始終使用此特定時段進行計算，在這種情況下，我們將始終使用過去 5 天的平均值來預測下一天的值。誤差低於簡單平均，但仍高於簡單模型，以下是測試模型:與訓練數據相似，滑動平均模型優於簡單平均模型，但尚未比簡單模型基礎有所增益。預測具有 2 個時滯的自相關性，並且相對於預測值有很大的方差誤差。指數滑動平均:上述簡單滑動平均模型具有同等地處理最後 X 個觀測值並完全忽略所有先前觀測值的特性。直觀地說，過去的數據應該逐漸打折，例如，理論上最近的觀測結果應該比第二近的更重要，而第二近的觀測應該比第三近的數據更重要，等等， 指數滑動平均（Exponential Moving Average，EMM）模型就是這樣做的。由於 α（alpha）是一個常數，其值介於 0 和 1 之間，因此我們將使用以下公式計算預測值：如果預測的第一個值是相應的當前值，其他值將更新為實際值與前一個時段的預測之差的 α 倍。當α為零時，我們根據第一個預測值得到一個常數，當 α 為 1 時，我們有一個簡單方法的模型，因為結果是前一個實際周期的值。EMM 預測中的平均數據周期為 1 /α。例如，當 α= 0.5 時，滯後相當於 2 個周期; 當 α= 0.2 時，滯後是 5 個周期; 當 α= 0.1 時，滯後是 10 個周期，依此類推。在這個模型中，我們將任意選用 α 值為 0.5 ，而你可以通過網格搜索算法查找在訓練集和驗證集中都中減少了錯誤的 α，數據大概應是這樣：這個模型的誤差與滑動平均的誤差相似，但是我們需要在測試集對模型進行驗證：在驗證數據中，目前為止的誤差在我們已經訓練過的模型中排名第二，而殘差圖的特徵與 5 天滑動平均模型的特徵非常相似。自回歸自回歸模型基本上是一個具有顯著相關滯後的線性回歸，首先要繪製自相關圖和部分自相關圖來驗證是否存在相關關係。下面是訓練序列的自相關圖和部分自相關圖，顯示了自回歸模型的特徵為具有 2 個時滯的顯著相關性:接下來我們將根據訓練數據創建模型，得到模型的係數後，將其乘以測試數據將要執行的值:這個模型與我們所訓練的其他模型相比，誤差最小，現在我們用它的係數對訓練數據進行逐步預測:注意，在測試數據中，誤差不會保持穩定，甚至會比簡單模型更差，可以看到圖中的預測值幾乎總是低於當前值，偏差測量顯示實際值比預測值高 50.19 BRL， 也許在訓練模型中調整一些參數，這種差異會減小。要改進這些模型，你可以應用多個轉換，例如本文中介紹的轉換，也可以添加外部變量作為預測源，但是，這已然超出本文內容了。4 結束語每個時間序列模型都有自己的特點，應該分別單獨分析，這樣我們就可以提取儘可能多的信息來做出好的預測，減少未來的不確定性。檢驗平穩度、轉換數據、在訓練數據中建立模型、驗證測試數據、檢驗殘差是建立良好時間序列預測的關鍵步驟。https://www.kaggle.com/leandrovrabelo/climate-change-forecast-sarima-model原文來源：https://towardsdatascience.com/basic-principles-to-create-a-time-series-forecast-6ae002d177a4本文來自公眾號讀者翻譯，歡迎各位童鞋向公號投稿，點擊下面圖片了解詳情！李潔，北京師範大學香港浸會大學聯合學院 數據科學系助教，香港科技大學電信學碩士。----End---- 「Python數據之道 」建立了讀者交流群，大家可以添加管理員微信進行加群