[Python][pmdarima] ARIMA時間序列模型學習
背景
前段時間參與了一個快消行業需求預測的項目。其中，用到了移動平均法、ARIMA、Xgboost等方法進行預測，如今打算總結一下ARIMA。
由於項目的銷售數據屬於私密數據，這裏用網上找的一份案例數據用於展現。
構建ARIMA模型能夠用到statsmodels庫和pmdarima庫。我這裏用的是pmdarima庫，這個庫有一個優勢是可以自動地對ARIMA模型的參數進行自動肯定。
html
一、建模步驟
二、季節性ARIMA模型
A R I M A ( p , d , q ) ( P , D , Q ) m ARIMA (p,d,q) (P,D,Q)_m ARIMA(p,d,q)(P,D,Q)m 該模型須要提供的參數有以下兩類：
（1）趨勢參數
p：趨勢自迴歸階數。
d：趨勢差分階數。
q：趨勢移動平均階數。
（2）季節性參數
P：季節性自迴歸階數。
D：季節性差分階數。
Q：季節性移動平均階數。
m：單個季節期間的時間步數。
python
1、模塊導入及數據讀取
一、模塊及數據
導入相關模塊函數
from model.arimaModel import * # 導入自定義的方法 import pandas as pd import matplotlib.pyplot as plt import pmdarima as pm from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # 畫圖定階 from statsmodels.tsa.stattools import adfuller # ADF檢驗 from statsmodels.stats.diagnostic import acorr_ljungbox # 白噪聲檢驗 import warnings warnings.filterwarnings("ignore") # 選擇過濾警告
讀取本地數據，須要數據能夠留言。這裏我用的數據是月份數據，若是有的數據集是日的數據，當發生太大波動性性很難發現趨勢時，可使用resample()進行每個月重採樣。學習
io = 'D:/PythonProject/ARIMA/AirPassengers.csv' time_series = pd.DataFrame(pd.read_csv(io, header=0, index_col='Month')) # 將字符串索引轉換成時間格式 time_series.index = pd.to_datetime(time_series.index) # 輸出前5個數據，此時Date爲時間索引 print("樣本量爲{}個".format(len(time_series))) print(time_series.head())
二、季節性分析
# 時間序列圖生成函數 def sequencePlot(data, sequencePlot_name): data.plot(marker='o') # 設置座標軸標籤 plt.xlabel("時間", fontsize=15) plt.ylabel("銷量", fontsize=15) # 設置圖表標題 plt.title("{}時間序列圖".format(sequencePlot_name), fontsize=15) # 圖例位置 plt.legend(loc='best') # 生成網格 plt.grid(linestyle='--') plt.show()
orign_seqname = "原始" sequencePlot(time_series, orign_seqname)
用matplotlib輸出原始時間序列圖以下。明顯地，咱們能夠看到數據有上升的趨勢，能夠定性地判斷該序列非平穩。另外，該時間序列有比較明顯的季節性趨勢，即基本上每年數據都呈現先上升再降低。
同時，可使用seasonal_decompose()進行分析，將時間序列分解成長期趨勢項(Trend)、季節性週期項(Seansonal)和殘差項(Resid)這三部分。該方法的原理見這篇文章 。
測試
from statsmodels.tsa.seasonal import seasonal_decompose
# 分解數據查看季節性 freq爲週期 ts_decomposition = seasonal_decompose(time_series['Passengers'], freq=12) ts_decomposition.plot() plt.show()
由下子圖，確實每一年的數據呈現先升後降的特徵。
ui
2、平穩性檢驗
一、Augmented Dickey-Fuller Test(ADF檢驗)
ARIMA模型要求時間序列是平穩的。所謂平穩性，其基本思想是：決定過程特性的統計規律不隨着時間的變化而變化。由平穩性的定義：對於一切 t , s t,s t,s， Y t Y_t Yt和 Y s Y_s Ys的協方差對於時間的依賴之和時間間隔 ∣ t − s ∣ |t-s| ∣t−s∣有關，而與實際的時刻t和s無關。所以，平穩過程能夠簡化符號，其中 C o v Cov Cov爲自協方差函數， C o r r Corr Corr爲自相關函數，記爲： γ t = C o v ( Y t , Y t − k ) ， ρ k = C o r r ( Y t , Y t − k ) \gamma_t=Cov(Y_t,Y_{t-k}) ，\rho_k=Corr(Y_t,Y_{t-k}) γt=Cov(Yt,Yt−k)，ρk=Corr(Yt,Yt−k)另外，平穩分爲嚴平穩和寬平穩。定義以下：
（1）嚴平穩：特別地，若是對一切時滯 k k k和時點點 t 1 ， t 2 ， . . . ， t k t_1，t_2，...，t_k t1，t2，...，tk，都有 Y t 1 , Y t 2 , . . . , Y t n Y_{t_1},Y_{t_2},...,Y_{t_n} Yt1,Yt2,...,Ytn與 Y t 1 − k , Y t 2 − k , . . . , Y t n − k Y_{t_1-k},Y_{t_2-k},...,Y_{t_n-k} Yt1−k,Yt2−k,...,Ytn−k的聯合分佈相同，則稱過程{ Y t Y_t Yt}爲嚴平穩的。
（2）寬平穩：知足1）均值函數在全部時間上恆爲常數；2） γ t , t − k = γ 0 , k \gamma_{t,t-k}=\gamma_{0,k} γt,t−k=γ0,k，對全部的時間 t t t和滯後 k k k。
若是經過時間序列圖來用肉眼觀看的話可能會存在一些主觀性。ADF檢驗（又稱單位根檢驗）是一種比較經常使用的嚴格的統計檢驗方法。
ADF檢驗主要是經過判斷時間序列中是否含有單位根：若是序列平穩，就不存在單位根；不然，就會存在單位根。
spa
詳細的介紹能夠看這篇博客
https://blog.csdn.net/FrankieHello/article/details/86766625/
.net
二、Python實現
from statsmodels.tsa.stattools import adfuller # ADF檢驗
# 該函數參考了其餘文章 def stableCheck(timeseries): # 移動12期的均值和方差 rol_mean = timeseries.rolling(window=12).mean() rol_std = timeseries.rolling(window=12).std() # 繪圖 fig = plt.figure(figsize=(12, 8)) orig = plt.plot(timeseries, color='blue', label='Original') mean = plt.plot(rol_mean, color='red', label='Rolling Mean') std = plt.plot(rol_std, color='black', label='Rolling Std') plt.legend(loc='best') plt.title('Rolling Mean & Standard Deviation') plt.show() # 進行ADF檢驗 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('ADF檢驗結果:') print(dfoutput)
stableCheck_result1 = stableCheck(time_series)
檢驗結果以下。ADF檢驗的 H 0 H_0 H0假設就是存在單位根，若是獲得的顯著性檢驗統計量小於三個置信度（10%，5%，1%），則對應有（90%，95，99%）的把握來拒絕原假設。有如下兩種方式來判斷：
（1）這裏將Test Statistic與1%、%五、%10不一樣程度拒絕原假設的統計值進行比較。當ADF Test result 同時小於 1%、5%、10%即說明很是好地拒絕該假設。本數據中，0.815369並無都小於三個level的統計值，因此判斷爲該時間序列time series非平穩。
（2）觀察p-value，是否很是接近於0。這裏p-value約爲0.99>0.05，所以不能拒絕原假設。
3d
3、序列平穩化
差分法
對於非平穩時間序列，能夠經過d階差分運算使變成平穩序列。從統計學角度來說，只有平穩性的時間序列才能避免「僞迴歸」的存在，纔有經濟意義。
在這裏，我先進行了一階差分，再進行了一次季節性差分，才經過ADF檢驗。
code
# 差分處理非平穩序列，先進行一階差分 time_series_diff1 = time_series.diff(1).dropna() # 在一階差分基礎上進行季節性差分差分 time_series_diff2 = time_series_diff1.diff(12).dropna() stableCheck_result2 = stableCheck(time_series_diff2)
結果以下，p-value小於0.05，且Test Statistic遠小於Critical Value (5%)的值，能夠有95%的機率認爲時間序列平穩。
所以，咱們能夠肯定季節性ARIMA模型的 d = 1 ， D = 1 ， m = 12 d=1，D=1，m=12 d=1，D=1，m=12。
4、白噪聲檢驗
接下來，對d階差分後的的平穩序列進行白噪聲檢驗，若該平穩序列爲非白噪聲序列，則能夠進行第五步。
一、白噪聲定義
對於一序列 e 1 , e 2 , e 3 , . . . , e t e_1,e_2,e_3,...,e_t e1,e2,e3,...,et，若知足： E ( e t ) = 0 E(e_t)=0 E(et)=0 V a r ( e t ) = σ 2 Var(e_t)=σ^2 Var(et)=σ2 C o v ( e t , e t + k ) = 0 （ k ≠ 0 ） Cov(e_t,e_{t+k})=0 （k≠0） Cov(et,et+k)=0（k=0） 那麼該序列是一個白噪聲序列。
詳細的介紹能夠看這篇博客
https://www.cnblogs.com/travelcat/p/11400307.html
二、Ljung-Box檢驗
Ljung-Box檢驗即LB檢驗，是時間序列分析中檢驗序列自相關性的方法。LB檢驗的Q統計量爲：
Q ( m ) = T ( T + 2 ) ∑ l = 1 m ρ ^ l 2 T − l Q(m)=T(T+2) \sum_{l=1}^m \frac{\hat{ρ}_l^2}{T-l} \quad Q(m)=T(T+2)l=1∑mT−lρ^l2 用來檢驗 m m m階滯後範圍內序列的自相關性是否顯著，或序列是否爲白噪聲，Q統計量服從自由度爲 m m m的卡方分佈。
三、Python實現
def whiteNoiseCheck(data): result = acorr_ljungbox(data, lags=1) temp = result[1] print('白噪聲檢驗結果：', result) # 若是temp小於0.05，則能夠以95%的機率拒絕原假設，認爲該序列爲非白噪聲序列；不然，爲白噪聲序列，認爲沒有分析意義 print(temp) return result
ifwhiteNoise = whiteNoiseCheck(time_series_diff2)
檢驗結果以下：返回的是統計量和p-value，其中，延遲1階的p-value約爲0.00033<0.05，能夠以95%的機率拒絕原假設，認爲該序列爲非白噪聲序列。
5、時間序列定階
def draw_acf(data): # 利用ACF判斷模型階數 plot_acf(data) plt.title("序列自相關圖(ACF)") plt.show() def draw_pacf(data): # 利用PACF判斷模型階數 plot_pacf(data) plt.title("序列偏自相關圖(PACF)") plt.show() def draw_acf_pacf(data): f = plt.figure(facecolor='white') # 構建第一個圖 ax1 = f.add_subplot(211) # 把x軸的刻度間隔設置爲1，並存在變量裏 x_major_locator = MultipleLocator(1) plot_acf(data, ax=ax1) # 構建第二個圖 ax2 = f.add_subplot(212) plot_pacf(data, ax=ax2) plt.subplots_adjust(hspace=0.5) # 把x軸的主刻度設置爲1的倍數 ax1.xaxis.set_major_locator(x_major_locator) ax2.xaxis.set_major_locator(x_major_locator) plt.show()
draw_acf_pacf(time_series_diff2)
具體怎麼對季節性ARIMA模型進行定階數能夠看這裏。
（1）肯定 d d d和 D D D的階數：當對原序列進行了 d d d階差分和 l a g lag lag爲 m m m的 D D D階差分後序列爲平穩序列，則能夠肯定 d d d， D D D， m m m的值。
（2）肯定 p p p， q q q和 P P P， Q Q Q的階數：1）首先對平穩化後的時間序列繪製ACF和PACF圖；2）而後，能夠經過觀察季節性 l a g lag lag處的拖尾/截尾狀況來肯定 P , Q P,Q P,Q的值；3）觀察短時間非季節性 l a g lag lag處的拖尾/截尾狀況來肯定 p , q p,q p,q的值。
以前已經在步驟三處肯定了 d d d， D D D， m m m。對於剩下的參數，將依據以下的ACF和PACF圖肯定。下圖的橫座標 l a g lag lag是以月份爲單位。
非季節性部分：對於 p p p，在 l a g = 1 lag=1 lag=1後，ACF圖拖尾，PACF圖截尾。一樣地，對於 q q q，在 l a g = 1 lag=1 lag=1後，ACF圖截尾，PACF圖拖尾。
季節性部分： P ， Q P，Q P，Q的肯定和非季節性同樣，不過須要記得滯後的間隔爲12。
【補充參考】
AR模型：自相關係數拖尾，偏自相關係數截尾；
MA模型：自相關係數截尾，偏自相關函數拖尾；
ARMA模型：自相關函數和偏自相關函數均拖尾。
6、構建ARIMA模型及預測
這裏使用了pmdarima.auto_arima()方法。這個方法能夠幫助咱們自動肯定 A R I M A ( p , d , q ) ( P , D , Q ) m ARIMA(p,d,q)(P,D,Q)_m ARIMA(p,d,q)(P,D,Q)m的參數，也就是能夠直接跳過上述的步驟2~5，直接輸入數據，設置auto_arima()中的參數則可。
以前咱們是經過觀察ACF、PACF圖的拖尾截尾現象來定階，可是這樣可能不許確。實際上，每每須要結合圖像擬合多個模型，經過模型的AIC、BIC值以及殘差分析結果來選擇合適的模型。
一、構建模型
import pmdarima as pm
將數據分爲訓練集data_train和測試集data_test 。
split_point = int(len(time_series) * 0.85) # 肯定訓練集/測試集 data_train, data_test = time_series[0:split_point], time_series[split_point:len(time_series)] # 使用訓練集的數據來擬合模型 built_arimamodel = pm.auto_arima(data_train, start_p=0, # p最小值 start_q=0, # q最小值 test='adf', # ADF檢驗確認差分階數d max_p=5, # p最大值 max_q=5, # q最大值 m=12, # 季節性週期長度，當m=1時則不考慮季節性 d=None, # 經過函數來計算d seasonal=True, start_P=0, D=1, trace=True, error_action='ignore', suppress_warnings=True, stepwise=False # stepwise爲False則不進行徹底組合遍歷 ) print(built_arimamodel.summary())
相關結果以下，從Model能夠看出肯定的參數與以前肯定的參數大體相同。
二、模型識別
built_arimamodel.plot_diagnostics() plt.show()
三、需求預測
# 進行多步預測,再與測試集的數據進行比較 pred_list = [] for index, row in data_test.iterrows(): # 輸出索引，值 # print(index, row) pred_list += [built_arimamodel.predict(n_periods=1)] # 更新模型，model.update()函數，不斷用新觀測到的 value 更新模型 built_arimamodel.update(row) # 預測時間序列之外將來的一次 predict_f1 = built_arimamodel.predict(n_periods=1) print('將來一期的預測需求爲：', predict_f1[0])
三、結果評估
對預測結果進行評價，指標解釋看這。
# ARIMA模型評價 def forecast_accuracy(forecast, actual): mape = np.mean(np.abs(forecast - actual)/np.abs(actual)) me = np.mean(forecast - actual) mae = np.mean(np.abs(forecast - actual)) mpe = np.mean((forecast - actual)/actual) # rmse = np.mean((forecast - actual)**2)**0.5 # RMSE rmse_1 = np.sqrt(sum((forecast - actual) ** 2) / actual.size) corr = np.corrcoef(forecast, actual)[0, 1] mins = np.amin(np.hstack([forecast[:, None], actual[:, None]]), axis=1) maxs = np.amax(np.hstack([forecast[:, None], actual[:, None]]), axis=1) minmax = 1 - np.mean(mins/maxs) return ({'mape': mape, 'me': me, 'mae': mae, 'mpe': mpe, 'rmse': rmse_1, 'corr': corr, 'minmax': minmax })
# 畫圖觀測實際與測試的對比 test_predict = data_test.copy() for x in range(len(test_predict)): test_predict.iloc[x] = pred_list[x] # 模型評價 eval_result = forecast_accuracy(test_predict.values, data_test.values) print('模型評價結果\n', eval_result) # 畫圖顯示 plt.plot(origin_timeseries, 'b-', lw=1, label='True Data', marker='*') plt.plot(test_predict, 'r-', lw=2, label='Prediction', marker='o') plt.title('RMSE:{}'.format(eval_result['rmse'])) plt.legend(loc='best') plt.grid() # 生成網格 plt.show()
紅色的點爲預測結果，藍色的點爲原來的數據。
RMSE約爲16.4，均方根偏差（Root Mean Square Error，RMSE）是預測值與真實值誤差的平方與觀測次數n比值的平方根指標，衡量觀測值與真實值之間的誤差。通常是越小越好，可是實際上須要進行結果對比才能判斷。
7、小結
在作項目時候，發現ARIMA模型對數據的要求仍是挺高的。其中，數據量太少不行，當時作的時候發現一些門店的數據量太少，並且只有一年左右的跨度，很難發現一些週期規律。另外就是對平穩性的要求。 目前對時間序列分析仍是初學階段，不少地方的解釋不夠深刻，後續還須要繼續學習、完善。