时间序列分析之ARIMA模型预测餐厅销量

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

from matplotlib.pylab import style                                   # 自定义图表风格
style.use('ggplot')

# 解决中文的显示问题
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False

from statsmodels.graphics.tsaplots import plot_acf,plot_pacf         # 自相关图、偏自相关图
from statsmodels.tsa.stattools import adfuller as ADF                # 平稳性检验
from statsmodels.stats.diagnostic import acorr_ljungbox              # 白噪声检验
import statsmodels.api as sm                                         # D-W检验,一阶自相关检验
from statsmodels.graphics.api import qqplot                          # 画QQ图,检验一组数据是否服从正态分布
from statsmodels.tsa.arima_model import ARIMA

sale = pd.read_excel('./arima_data.xls', index_col='日期')
sale.head()

sale.info()
print('-----')
sale.销量 = sale.销量.astype('float')
sale.info()

plt.figure(figsize=(10,5))
sale.plot()
plt.show()
 
#解读：具有单调递增趋势，则是非平稳序列。

plot_acf(sale, lags=35).show()
 
#解读：自相关系数长期大于零，没有趋向于零，说明序列间具有很强的长期相关性。

#方法：单位根检验
 
print('原始序列的ADF检验结果为：',ADF(sale.销量))
 
#解读：P值大于显著性水平α（0.05），接受原假设（非平稳序列），说明原始序列是非平稳序列。

d1_sale = sale.diff(periods=1, axis=0).dropna()
 
#时序图
plt.figure(figsize=(10,5))
d1_sale.plot()
plt.show()
#解读：在均值附件比较平稳波动
 
#自相关图
plot_acf(d1_sale, lags=34).show()
#解读：有短期相关性，但趋向于零。
 
#平稳性检验
print('原始序列的ADF检验结果为：',ADF(d1_sale.销量))
 
#解读：P值小于显著性水平α（0.05），拒绝原假设（非平稳序列），说明一阶差分序列是平稳序列。

print('一阶差分序列的白噪声检验结果为：',acorr_ljungbox(d1_sale, lags=1))#返回统计量、P值
 
#解读：p值小于0.05，拒绝原假设（纯随机序列），说明一阶差分序列是非白噪声序列。

d1_sale = sale.diff(periods=1, axis=0).dropna()
 
#自相关图
plot_acf(d1_sale, lags=34).show()
 
#解读：有短期相关性，但趋向于零。
 
#偏自相关图
plot_pacf(d1_sale, lags=10).show()
 
 
#偏自相关图
plot_pacf(d1_sale, lags=17).show()
 
#解读：自相关图，1阶截尾；偏自相关图，拖尾。则ARIMA(p,d,q)=ARIMA(0,1,1)

pmax = int(len(d1_sale) / 10) #一般阶数不超过length/10
qmax = int(len(d1_sale) / 10) #一般阶数不超过length/10
pmax
qmax

bic_matrix = []
for p in range(pmax + 1):
    tmp = []
    for q in range(qmax + 1):
        try:
            tmp.append(ARIMA(tuple(sale), (p, 1, q)).fit().bic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)

bic_matrix = pd.DataFrame(bic_matrix)
bic_matrix

bic_matrix.stack()

p,q=bic_matrix.stack().idxmin() #最小值的索引
print('用BIC方法得到最优的p值是%d,q值是%d'%(p,q))

pmax = int(len(d1_sale )/ 10) #一般阶数不超过length/10
qmax = int(len(d1_sale) / 10) #一般阶数不超过length/10
 
aic_matrix = []
for p in range(pmax + 1):
    tmp = []
    for q in range(qmax + 1):
        try:
            tmp.append(ARIMA(tuple(sale), (p, 1, q)).fit().aic)
        except:
            tmp.append(None)
    aic_matrix.append(tmp)
aic_matrix = pd.DataFrame(aic_matrix)
p,q = aic_matrix.stack().idxmin() #最小值的索引

print('用AIC方法得到最优的p值是%d,q值是%d'%(p,q))

#创建模型
model = ARIMA(tuple(sale), (0, 1, 1)).fit()
#查看模型报告
model.summary2()

resid = model.resid
 
#自相关图
plot_acf(resid, lags=35).show()
 
#解读：有短期相关性，但趋向于零。
 
#偏自相关图
plot_pacf(resid, lags=10).show()
 
#偏自相关图
plot_pacf(resid, lags=17).show()

qqplot(resid, line='q', fit=True).show() 
 
#解读：残差服从正态分布，均值为零，方差为常数

print('D-W检验的结果为：',sm.stats.durbin_watson(resid.values))  
 
#解读：不存在一阶自相关

# 方法一
print('残差序列的白噪声检验结果为：',acorr_ljungbox(resid,lags=1))#返回统计量、P值
 
#解读：残差是白噪声

# 方法二
confint,qstat,pvalues = sm.tsa.acf(resid.values, qstat=True) #qstat is Ljung-Box Q-Statistic. confint is  Confidence intervals for the ACF
data = np.c_[range(1,36), confint[1:], qstat, pvalues]
table = pd.DataFrame(data, columns=['lag', "confint", "qstat", "pvalues(>Q)"])
print(table.set_index('lag'))

#预测
print('未来7天的销量预测：')
model.forecast(7) #预测、标准差、置信区间

forecast = pd.Series(model.forecast(7)[0], index=pd.date_range('2015-2-7', periods=7, freq='D'))
forecast

data = pd.concat((sale, forecast), axis=0)
data.columns = ['销量', '未来7天销量']
plt.figure(figsize = (10,5))
data.plot()
plt.show()