如果你在尋找時間序列是什么？如何實現時間序列？那么請看這篇博客，將以通俗易懂的語言，全面的闡述時間序列及其python實現。
時間序列算法理論詳見我的另一篇博客：時間序列算法理論及python實現 - 知-青 - 博客園
5 Python實現ARIMA模型
下面應用以上理論知識，對表6中2015/1/1~2015/2/6某餐廳的銷售數據進行建模。
就餐飲企業而言，經常會碰到如下問題。
由於餐飲行業是勝場和銷售同時進行的，因此銷售預測對於餐飲企業十分必要。如何基於菜品歷史銷售數據，做好餐銷售預測，以便減少菜品脫銷現象和避免因備料不足而造成的生產延誤，從而減少菜品生產等待時間，提供給客戶更優質的服務，同事可以減少安全庫存量，做到生產准時制，降低物流成本
餐飲銷售預測可以看作是基於時間序列的短期數據預測，預測對象為具體菜品銷售量
表6 原序列數據
5.1 環境配置
1 import pandas as pd 2 import matplotlib.pyplot as plt 3 from matplotlib.pylab import style 4 from statsmodels.tsa.stattools import adfuller as ADF 5 from statsmodels.stats.diagnostic import acorr_ljungbox # 白噪聲檢驗 6 from statsmodels.tsa.arima_model import ARIMA 7 import statsmodels.tsa.api as smt 8 import seaborn as sns 9 style.use('ggplot') 10 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用來正常顯示中文標簽 11 plt.rcParams['axes.unicode_minus'] = False # 用來正常顯示負號
要安裝的環境有點小多，需要提前安裝好。
5.2 導入數據
1 # 參數初始化 2 discfile = './data/arima_data.xls' 3 forecastnum = 5 4 5 # 讀取數據，指定日期列為指標，Pandas自動將“日期”列識別為Datetime格式 6 data = pd.read_excel(discfile, index_col=u'日期')
代碼和數據將會公布在Github，請到文末鏈接。
5.3 檢驗序列的平穩性
1 # 時序圖 2 import matplotlib.pyplot as plt 3 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用來正常顯示中文標簽 4 plt.rcParams['axes.unicode_minus'] = False # 用來正常顯示負號 5 data.plot() 6 plt.show() 7 8 # 自相關圖 9 from statsmodels.graphics.tsaplots import plot_acf 10 plot_acf(data).show() 11 12 # 平穩性檢測 13 from statsmodels.tsa.stattools import adfuller as ADF 14 print(u'原始序列的ADF檢驗結果為：', ADF(data[u'銷量'])) 15 # 返回值依次為adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
圖3 原始序列的時序圖
圖4 原始序列的自相關圖
原始時間序列的單位根檢驗
表7 原始序列的單位根檢驗
圖3時序圖顯示該序列具有明顯的單調遞增趨勢，可以判斷為是非平穩序列；圖4的自相關圖顯示自相關系數長期大於零，說明序列間具有很強的長期相關性；表7單位根檢驗統計量對應的P值顯著大於0.05，最終將該序列判斷為非平穩序列（非平穩序列一定不是白噪聲序列）。
5.4 對原始序列進行一階差分，並進行平穩性和白噪聲檢驗
5.4.1 對一階差分后的序列再次做平穩性判斷
1 # 差分后的結果 2 D_data = data.diff().dropna() 3 D_data.columns = [u'銷量差分'] 4 D_data.plot() # 時序圖 5 plt.show() 6 plot_acf(D_data).show() # 自相關圖 7 from statsmodels.graphics.tsaplots import plot_pacf 8 plot_pacf(D_data).show() # 偏自相關圖 9 print(u'差分序列的ADF檢驗結果為：', ADF(D_data[u'銷量差分'])) # 平穩性檢測
圖5 一階差分之后序列的時序圖
圖6 一階差分之后序列的自相關圖
一階差分之后序列的單位根檢驗
表8 一階差分之后序列的單位根檢驗
結果顯示，一階差分之后的序列的時序圖在均值附近比較平穩的波動、自相關圖有很強的短期相關性、單位根檢驗P值小於0.05，所以一階差分之后的序列是平穩序列。
5.4.2 對一階差分后的序列做白噪聲檢驗（結果見表5-28）
from statsmodels.stats.diagnostic import acorr_ljungbox print(u'差分序列的白噪聲檢驗結果為：', acorr_ljungbox(D_data, lags=1)) # 返回統計量和p值
表9 一階差分后的序列的白噪聲檢驗
輸出的P值遠遠小於0.05，所以一階差分之后的序列是平穩非白噪聲序列。
5.5 對一階差分之后的平穩非白噪聲序列擬合ARMA模型
下面進行模型定階，模型定階就是確定p和q。
5.5.1 人為識別實現模型定階
一階差分后自相關圖（見圖6）顯示出1階截尾，偏自相關圖顯示出拖尾性，所以可以考慮用MA(1)模型擬合1階差分后的序列，即對原始序列建立ARIMA(0,1,1)模型。
圖7 一階差分后序列的偏自相關圖
5.5.2 相對最優模型識別
計算ARMA(p,q)。當p和q均小於等於3的所有組合的BIC信息量，取其中BIC信息量達到最小的模型階數。
1 from statsmodels.tsa.arima_model import ARIMA 2 3 data[u'銷量'] = data[u'銷量'].astype(float) 4 # 定階 5 pmax = int(len(D_data) / 10) # 一般階數不超過length/10 6 qmax = int(len(D_data) / 10) # 一般階數不超過length/10 7 bic_matrix = [] # bic矩陣 8 for p in range(pmax + 1): 9 tmp = [] 10 for q in range(qmax + 1): 11 try: # 存在部分報錯，所以用try來跳過報錯。 12 tmp.append(ARIMA(data, (p, 1, q)).fit().bic) 13 except: 14 tmp.append(None) 15 bic_matrix.append(tmp) 16 17 bic_matrix = pd.DataFrame(bic_matrix) # 從中可以找出最小值 18 19 p, q = bic_matrix.stack().idxmin() # 先用stack展平，然后用idxmin找出最小值位置。 20 print(u'BIC最小的p值和q值為：%s、%s' % (p, q))
計算完成BIC矩陣如下（繪制程序在主程序，以上程序僅僅只有計算）
圖8 矩陣熱度圖
P值為0、q值為1時最小BIC值為：430.1374。p、q定階完成！
5.6 模型檢驗
用AR(1)模型擬合一階差分后的序列，即對原始序列建立ARIMA(0,1,1)模型。雖然兩種方法建立的模型是一樣，但模型是非唯一的，可以檢驗ARIMA(1,1,0)和ARIMA(1,1,1)，這兩個模型也能通過檢驗。
下面對一階差分后的序列擬合AR(1)模型進行分析。
（1）模型檢驗。殘差為白噪聲序列，p值為：0.627016
（2）參數檢驗和參數估計見表10。
表10 模型參數
5.7 模型預測
1 model = ARIMA(data, (p, 1, q)).fit() # 建立ARIMA(0, 1, 1)模型 2 model.summary2() # 給出一份模型報告 3 model.forecast(5) # 作為期5天的預測，返回預測結果、標准誤差、置信區間。
應用ARIMA(0,1,1)對表11中的2015/1/1~2015/2/6某餐廳的銷售數據做為期5天的預測，結果如下。
表11 預測結果
需要說明的是，利用模型向前預測的時期越長，預測誤差將會越大，這是時間預測的典型特點。
6 文獻
王黎明，王連等. 應用時間序列分析
張良均，王路，譚立雲，蘇劍林. Python數據分析與挖掘實戰
Complete guide to create a Time Series Forecast (with Codes in Python)
時間序列預測如何變成有監督學習問題？ - 雲+社區 - 騰訊雲
7 附錄：程序及數據
說明：為了方便調用，我把所有程序都封裝成函數，調用極其方便只用改動很小的參數。
1 # -*- coding:utf-8 -*- 2 # @Time : 2018/7/11 15:18 3 # @Author : yuanjing liu 4 # @Email : lauyuanjing@163.com 5 # @File : ts_arima.py 6 # @Software: PyCharm 7 # arima時序模型 8 9 import pandas as pd 10 import matplotlib.pyplot as plt 11 from matplotlib.pylab import style 12 from statsmodels.tsa.stattools import adfuller as ADF 13 from statsmodels.stats.diagnostic import acorr_ljungbox # 白噪聲檢驗 14 from statsmodels.tsa.arima_model import ARIMA 15 import statsmodels.tsa.api as smt 16 import seaborn as sns 17 style.use('ggplot') 18 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用來正常顯示中文標簽 19 plt.rcParams['axes.unicode_minus'] = False # 用來正常顯示負號 20 21 22 # 對原始數據進行ACF、PACF檢驗 23 def tsplot(y, lags=None, title='', figsize=(14, 8)): 24 fig = plt.figure(figsize=figsize) 25 layout = (2, 2) 26 ts_ax = plt.subplot2grid(layout, (0, 0)) 27 hist_ax = plt.subplot2grid(layout, (0, 1)) 28 acf_ax = plt.subplot2grid(layout, (1, 0)) 29 pacf_ax = plt.subplot2grid(layout, (1, 1)) 30 31 y.plot(ax=ts_ax) 32 ts_ax.set_title(title) 33 y.plot(ax=hist_ax, kind='hist', bins=25) 34 hist_ax.set_title('Histogram') 35 smt.graphics.plot_acf(y, lags=lags, ax=acf_ax) 36 smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax) 37 [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]] 38 sns.despine() 39 fig.tight_layout() 40 plt.show() 41 return ts_ax, acf_ax, pacf_ax 42 43 44 # 平穩性檢測(P值大於0.05，則存在單位根，是不平穩時間序列） 45 # adf_jy返回值依次為adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore 46 def steady(sdata): 47 adf_jy = ADF(sdata) # data[u'銷量'] 48 adf_p_value = adf_jy[1] 49 return adf_jy, adf_p_value 50 51 52 # 白噪聲檢驗 53 def w_noise(wdata): 54 w_noise = acorr_ljungbox(wdata, lags=1) # 返回統計量和p值 55 w_p_value = float(w_noise[1]) 56 return w_noise, w_p_value 57 58 59 # 差分后的結果（如果不平穩） 60 def ts_diff(ddata): 61 D_data = ddata.diff().dropna() # dropna是缺失值處理 62 D_data.columns = [u'1階差分'] 63 return D_data 64 65 66 def ts_arima(tsdata, forenum=5): 67 tsdata = tsdata.astype(float) 68 # 定階 69 D_data = ts_diff(tsdata) 70 pmax = int(len(D_data) / 10) # 一般階數不超過length/10 71 qmax = int(len(D_data) / 10) # 一般階數不超過length/10 72 bic_matrix = [] # bic矩陣 73 for p in range(pmax + 1): 74 tmp = [] 75 for q in range(qmax + 1): 76 try: # 存在部分報錯，所以用try來跳過報錯。 77 tmp.append(ARIMA(tsdata, (p, 1, q)).fit().bic) 78 except: 79 tmp.append(None) 80 bic_matrix.append(tmp) 81 82 bic_matrix = pd.DataFrame(bic_matrix) # 從中可以找出最小值 83 84 # 可視化BIC矩陣 85 fig, ax = plt.subplots(figsize=(10, 8)) 86 ax = sns.heatmap(bic_matrix, 87 mask=bic_matrix.isnull(), 88 ax=ax, 89 annot=True, 90 fmt='.2f', 91 ) 92 ax.set_title('BIC') 93 plt.show() 94 95 p, q = bic_matrix.stack().idxmin() # 先用stack展平，然后用idxmin找出最小值位置。 96 # print(u'BIC最小的p值和q值為：%s、%s' % (p, q)) 97 98 model = ARIMA(tsdata, (p, 1, q)).fit() # 建立ARIMA(0, 1, 1)模型 99 summary = model.summary2() # 給出一份模型報告 100 forecast = model.forecast(forenum) # 作為期forenum天的預測，返回預測結果、標准誤差、置信區間。 101 return bic_matrix, p, q, model, summary, forecast 102 103 104 # 測試 105 # 讀取數據 106 discfile = '../data/arima_data.xls' 107 forecastnum = 5 108 data = pd.read_excel(discfile, index_col=u'日期') 109 ddata = data[u'銷量'] 110 # 檢驗 111 ts_ap = tsplot(ddata, title='A Given Training Series', lags=20) # ACF 和 PACF 檢驗 112 s_total, s_p = steady(ddata) # 平穩性檢驗 113 w_total, w_p = w_noise(ddata) 114 # 差分 115 dif_data = ts_diff(ddata) 116 # arima模型 117 bic_matrix1, p1, q1, model1, summary, forecast = ts_arima(ddata)
轉載說明
1、本人博客純屬技術積累和分享，歡迎大家評論和交流以求共同進步。
2、在無明確說明下，博客可以轉載以供個人學習和交流，但是要附上出處。
3、如果原創博客使用涉及商業/公司行為請郵件(1547364995@qq.com)告知，一般情況均會及時回復同意。
4、如果個人博客中涉及他人文章我會盡力注明出處，但受限於能力並不能保證所有引用之處均能夠注明出處，如有冒犯，請您及時郵件告知以便修改，並於此提前向您道歉。
5、轉載過程中如有涉及他人作品請您與作者聯系。
6、所有文章（不限於原創）僅為個人見解，個人只能盡量保證正確，如有錯誤您需要自負責任，並請您留下評論提出錯誤之處以便及時更正，惠澤他人，謝謝