時間序列算法理論及python實現（2-python實現）

時間 2019-11-17

原文原文鏈接

若是你在尋找時間序列是什麼？如何實現時間序列？那麼請看這篇博客，將以通俗易懂的語言，全面的闡述時間序列及其python實現。html

時間序列算法理論詳見個人另外一篇博客：時間序列算法理論及python實現 - 知-青 - 博客園python

5 Python實現ARIMA模型

下面應用以上理論知識，對錶6中2015/1/1~2015/2/6某餐廳的銷售數據進行建模。git

就餐飲企業而言，常常會碰到以下問題。github

因爲餐飲行業是勝場和銷售同時進行的，所以銷售預測對於餐飲企業十分必要。如何基於菜品歷史銷售數據，作好餐銷售預測，以便減小菜品脫銷現象和避免因備料不足而形成的生產延誤，從而減小菜品生產等待時間，提供給客戶更優質的服務，同事能夠減小安全庫存量，作到生產準時制，下降物流成本算法

餐飲銷售預測能夠看做是基於時間序列的短時間數據預測，預測對象爲具體菜品銷售量api

表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  # 用來正常顯示負號

要安裝的環境有點小多，須要提早安裝好。app

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 原始序列的時序圖ide

圖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數據分析與挖掘實戰

python時間序列分析 - 大熊貓淘沙 - 博客園

機器學習_時間序列預測分析算法 | FEI's Blog

ARIMA模型的拖尾截尾問題 - CSDN博客

時間序列初級理論篇 - CSDN博客

大數據時間序列分析、建模與預測系列第一部分: 數據準備

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)