【譯】時間序列預測初學者指南
這篇文章是《 ofollow,noindex" target="_blank">基於R語言的時間序列建模完整教程 》的後續文章，不同的是本文采用Python來進行講解。本文在原文基礎上刪除和修改了部分內容，如遇到不明白的，請檢視 原文 。
在 pandas 載入時間序列資料
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
data = pd.read_csv('AirPassengers.csv', parse_dates=['TravelDate'], index_col='TravelDate')
print(data.head())
ts = data['Passengers']
ts.head(10)
plt.plot(ts)
關於Pandas如何載入資料的，請檢視： 如何使用Pandas讀取Excel和CSV檔案資料
檢驗時間序列穩定性
前面的文章中已經講到穩定性的評判標準主要是均值、方差和協商差。具體可以通過下述兩種方法進行測試：
1、繪製滾動統計：繪製移動平均數和移動方差，觀察它是否隨著時間變化。
2、Dickey-Fuller檢驗：測試結果由測試統計量和一些置信區間的臨界值組成。如果“測試統計量”少於“臨界值”，並認為序列是穩定的。
以下方法定義了以上兩種方式(注意，為了保持單位和平均數相似，這裡採用了標準差來代替方差)
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(window=12).std()
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
test_stationarity(ts)
這裡解釋一下，DF測試的值怎麼看：
- Test Statistic的值如果比Critical Value (5%)小則滿足穩定性需求
- p-value越低（理論上需要低於05）證明序列越穩定。
這裡這個結果表明這些序列很不穩定，所以接下來考慮如何處理資料，使得序列相對穩定。
使時間序列平穩
有兩個主要的因素導致序列不穩定：
- 趨勢 Trend
- 季節性 Seasonality
消除趨勢
消除趨勢的第一個方法是轉換。在本例中我們可以清楚地看到有一個顯著的趨勢。所以我們可以通過變換，懲罰較高值而不是較小值。這可以採用取對數log,平方根,立方跟等。讓我們簡單在這兒轉換一個對數。
ts_log = np.log(ts) plt.plot(ts_log)
這裡我們可以明顯看到一個上升的趨勢，但是混雜在噪音當中，所以需要去除雜音。這裡簡單平滑一下資料。平滑的視窗取值12，因為一年有12個月。
在這個簡單的例子中,很容易看到一個向前的資料趨勢。但是它表現的不是很直觀。所以我們可以使用一些技術來估計或對這個趨勢建模,然後將它從序列中刪除。這裡有很多方法,最常用的有:
- 聚合-取一段時間的平均值（月/周平均值）
- 平滑-取滾動平均數
- 多項式迴歸分析-適合的迴歸模型
我在這兒討論將平滑,你也應該嘗試其他可以解決的問題的技術。平滑是指採取滾動估計,即考慮過去的幾個例項。有各種方法可以解決這些問題，但我將主要討論以下兩個。
移動平均數
在這個方法中，根據時間序列的頻率採用“K”連續值的平均數。我們可以採用過去一年的平均數，即過去12個月的平均數。
moving_avg = ts_log.rolling(12).mean() plt.plot(ts_log) plt.plot(moving_avg, color='red')
紅色表示了滾動平均數。讓我們從原始序列中減去這個平均數。注意，從我們採用過去12個月的值開始，滾動平均法還沒有對前11個月的值定義：
ts_log_moving_avg_diff = ts_log - moving_avg ts_log_moving_avg_diff.head(12)
TravelDate 1949-01-01 NaN 1949-02-01 NaN 1949-03-01 NaN 1949-04-01 NaN 1949-05-01 NaN 1949-06-01 NaN 1949-07-01 NaN 1949-08-01 NaN 1949-09-01 NaN 1949-10-01 NaN 1949-11-01 NaN 1949-12-01 -0.065494 Name: Passengers, dtype: float64
注意前11個月是NaN值，現在讓我們對這11個月排除後測試穩定性。
ts_log_moving_avg_diff.dropna(inplace=True) test_stationarity(ts_log_moving_avg_diff)
從上面的測試結果看，已經得到了一個穩定的序列。但是，這個方法有一個缺陷：需要先設定平滑的視窗期，實際在應用的時候像股票這樣的資料很難去設定視窗期，所以，我們採取“加權移動平均法”可以對最近的值賦予更高的權重。
加權移動平均法
指數加權移動平均法是很受歡迎的方法，所有的權重被指定給先前的值連同衰減係數。這可以通過pandas實現：
expwighted_avg = ts_log.ewm(halflife=12).mean() plt.plot(ts_log) plt.plot(expwighted_avg, color='red')
注意,這裡使用了引數“halflife”來定義指數衰減量。這只是一個假設，很大程度上取決於業務領域。其他引數，如跨度和質心也可以用來定義衰減。讓我們再來檢測下新得到的序列的穩定性：
ts_log_ewma_diff = ts_log - expwighted_avg test_stationarity(ts_log_ewma_diff)
檢測結果比移動平均的效果更好（Test Statistic的值比1%的臨界值還小），另外不會出現前面11個月資料遺漏問題。
消除季節性
之前討論來了簡單的趨勢減少技術不能在所有情況下使用，特別是在高季節性情況下。讓我們談論一下兩種消除趨勢和季節性的方法。
- 差分：採用一個特定時間差的差值
- 分解：建立有關趨勢和季節性的模型和從模型中刪除它們。
差分
處理趨勢和季節性的最常見的方法之一就是差分法。在這種方法中,我們採用特定瞬間和它前一個瞬間的不同的觀察結果。這主要是在提高平穩性。pandas可以實現一階差分：
ts_log_diff = ts_log.diff() //ts_log_diff = ts_log - ts_log.shift() plt.plot(ts_log_diff)
圖中可以看出很大程度上減少了趨勢。讓我們再來檢測下：
ts_log_diff.dropna(inplace=True) test_stationarity(ts_log_diff)
我們可以看到平均數和標準差隨著時間有小的變化。同時，DF檢驗統計量小於10% 的臨界值，因此該時間序列在90%的置信區間上是穩定的。我們同樣可以採取二階或三階差分在具體應用中獲得更好的結果。這些方法你可以自己嘗試。
分解
在這種方法中，趨勢和季節性都分別建模，並返回序列的其餘部分。
from statsmodels.tsa.seasonal import seasonal_decompose decomposition = seasonal_decompose(ts_log) trend = decomposition.trend seasonal = decomposition.seasonal residual = decomposition.resid plt.subplot(411) plt.plot(ts_log, label='Original') plt.legend(loc='best') plt.subplot(412) plt.plot(trend, label='Trend') plt.legend(loc='best') plt.subplot(413) plt.plot(seasonal,label='Seasonality') plt.legend(loc='best') plt.subplot(414) plt.plot(residual, label='Residuals') plt.legend(loc='best') plt.tight_layout()
趨勢和季節性，還有殘差值都被分解出來，然後我們計算殘差值的穩定性。
ts_log_decompose = residual ts_log_decompose.dropna(inplace=True) test_stationarity(ts_log_decompose)
DF測試統計量明顯低於1%的臨界值，這樣時間序列是非常接近穩定。
預測時間序列
我們看到，使用不同的技術都可以是的序列變得穩定，接下里我們以差分處理後的序列搭建模型，因為其相對來說更容易新增噪音及季節性，讓其回到預測值。。在執行趨勢和季節性預測上，有兩種情況：
- 不含依賴值的嚴格穩定系列。簡單的情況下,我們可以建立殘差模型作為白噪音（指功率譜密度在整個頻域內均勻分佈的噪聲）。但這是非常罕見的。
- 序列含有明顯的依賴值。在這種情況下,我們需要使用一些統計模型像ARIMA（差分自迴歸移動平均模型）來預測資料。
ARIMA（Auto-Regressive Integrated Moving Averages）模型。平穩時間序列的ARIMA預測的只不過是一個線性方程(如線性迴歸)。模型有三個主要引數：
- Number of AR (Auto-Regressive) terms (p)： 現在點使用多少個過往資料計算。AR條件僅僅是因變數的滯後。如：如果P等於5，那麼預測x（t）將是x（t-1）…x（t-5）。
- Number of MA (Moving Average) terms (q)：使用多少個過往的殘餘錯誤值。MA條件是預測方程的滯後預測錯誤。如：如果q等於5，預測x（t）將是e（t-1）…e(t-5)，e(i)是移動平均叔在第i個瞬間和實際值的差值。
- Number of Differences (d)：為時間序列成為平穩時所做的差分次數。有非季節性的差值，這種情況下我們採用一階差分。
在這裡一個重要的問題是如何確定“p”和“q”的值。我們使用兩張圖示來確定這些數字。
- 自相關函式（ACF）：這是時間序列和它自身滯後版本之間的相關性的測試。比如在自相關函式可以比較時間的瞬間‘t1’…’t2’以及序列的瞬間‘t1-5’…’t2-5’ (t1-5和t2 是結束點)。
- 部分自相關函式(PACF):這是時間序列和它自身滯後版本之間的相關性測試，但是在預測（已經通過比較干預得到解釋）的變數後。如：滯後值為5，它將檢查相關性，但是會刪除從滯後值1到4得到的結果。
時間序列的自迴歸函式和部分自迴歸函式可以在差分後繪製為：
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
在這個例子中，0刻度線上線的2條虛線為置信區間，用來確定“p”和“q”的值。
- p-部分自相關函式表第一次截斷的上層置信區間是滯後值。如果你仔細看，該值是p=0
- q-自相關函式表第一次截斷的上層置信區間是滯後值。如果你仔細看，該值是q=1
現在開始搭建ARIMA的三個模型，及計算各自的RSS（這裡的RSS是指殘差值,而不是實際序列）
自迴歸（AR）模型
model = ARIMA(ts_log, order=(0, 1, 0))
results_AR = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))
移動平均數（MA）模型
model = ARIMA(ts_log, order=(0, 1, 1))
results_MA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))
ARMA組合
在本例中，由於AR的p值為0，所以組合的結果與MA一致：
model = ARIMA(ts_log, order=(0, 1, 1))
results_ARIMA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))
轉化為原資料空間
前面模型採用的資料都是轉化後的資料，為了得到最終的結果，需要將資料轉化為原資料空間。第一步將預測結果儲存為序列。
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True) print(predictions_ARIMA_diff.head())
TravelDate 1949-02-01 0.009726 1949-03-01 0.020485 1949-04-01 0.034537 1949-05-01 -0.005924 1949-06-01 -0.006085 dtype: float64
注意這些是從‘1949-02-01’開始而不是第一個月。為什麼?這是因為我們將第一個月份取為滯後值，一月前面沒有可以減去的元素。一個簡單的解決方法是先確定索引的累積和，然後將它新增到基數。累積和的計算方式如下：
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum() print(predictions_ARIMA_diff_cumsum.head()) TravelDate 1949-02-01 0.009726 1949-03-01 0.030211 1949-04-01 0.064748 1949-05-01 0.058824 1949-06-01 0.052739 dtype: float64
將差分轉換為對數尺度的方法是這些差值連續地新增到基數。（第一個元素是基數本身）
predictions_ARIMA_log = pd.Series(ts_log.iloc[0], index=ts_log.index) predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0) predictions_ARIMA_log.head()
TravelDate 1949-01-01 4.718499 1949-02-01 4.728225 1949-03-01 4.748710 1949-04-01 4.783247 1949-05-01 4.777323 dtype: float64
最後一步是將指數與原序列比較。
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))
總結
本週，主要是對一些概念和流程有了更加深入的一些講解，但是對於如何預測未來一年等方法沒有講到位，後續繼續需要學習。