Python时间序列分析：使用ARIMA模型预测销量

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller as ADF
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox

# 读取数据，指定日期为索引列
data = pd.read_excel('D:\\M_hua\\text.xlsx', index_col=0, header=None)  # 修改路径为你的文件路径，指定第一列为索引列

# 将数据转换为数值类型
data[1] = pd.to_numeric(data[1], errors='coerce')

# 绘图过程中
# 用来正常显示中文标签
plt.rcParams['font.sans-serif'] = ['SimHei']
# 用来正常显示负号
plt.rcParams['axes.unicode_minus'] = False

# 查看趋势图
data[1].plot()  # 有增长趋势，不平稳
plt.title('销量趋势图')
plt.xlabel('日期')
plt.ylabel('销量')
plt.show()

# 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
# 自相关图
plot_acf(data[1])
plt.show()

# 1 平稳性检测

def tagADF(t):
    result = pd.DataFrame(index=[
        'Test Statistic Value', 'p-value', 'Lags Used',
        'Number of Observations Used',
        'Critical Value(1%)', 'Critical Value(5%)', 'Critical Value(10%)'
    ], columns=['销量'])
    result['销量']['Test Statistic Value'] = t[0]
    result['销量']['p-value'] = t[1]
    result['销量']['Lags Used'] = t[2]
    result['销量']['Number of Observations Used'] = t[3]
    result['销量']['Critical Value(1%)'] = t[4]['1%']
    result['销量']['Critical Value(5%)'] = t[4]['5%']
    result['销量']['Critical Value(10%)'] = t[4]['10%']
    return result


print('原始序列的ADF检验结果为:', tagADF(ADF(data[1])))  # 添加标签后展现

# 2 进行数据差分，一般一阶差分就可以
D_data = data[1].diff(1).dropna()
D_data = D_data.to_frame(name='销量差分')

# 差分图趋势查看
D_data.plot()
plt.title('销量差分趋势图')
plt.xlabel('日期')
plt.ylabel('销量差分')
plt.show()

# 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
# 自相关图
plot_acf(D_data)
plt.show()

# 偏自相关图
plot_pacf(D_data)
plt.show()

# 3 平稳性检测
print('差分序列的ADF检验结果为：', tagADF(ADF(D_data['销量差分'])))

# 解释：Test Statistic Value值小于两个水平值，p值显著小于0.05，一阶差分后序列为平稳序列。

# 4 白噪声检验
# 返回统计量和p值
print('差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1))  # 分别为stat值（统计量）和P值
# P值小于0.05，所以一阶差分后的序列为平稳非白噪声序列。

# 5 p，q定阶
# 一般阶数不超过length/10
pmax = int(len(D_data) / 10)
# 一般阶数不超过length/10
qmax = int(len(D_data) / 10)

# bic矩阵
bic_matrix = []
for p in range(pmax + 1):
    tmp = []
    for q in range(qmax + 1):
        # 存在部分报错，所以用try来跳过报错。
        try:
            model = ARIMA(data[1], order=(p, 1, q)).fit()
            tmp.append(model.bic)
        except Exception:
            tmp.append(None)
    bic_matrix.append(tmp)

# 从中可以找出最小值
bic_matrix = pd.DataFrame(bic_matrix)

# 先用stack展平，然后用idxmin找出最小值位置。
p, q = bic_matrix.stack().idxmin()

print('BIC最小的p值和q值为：%s、%s' % (p, q))
# 取BIC信息量达到最小的模型阶数，结果p为0，q为1，定阶完成。

# 6 建立模型和预测
model = ARIMA(data[1], order=(p, 1, q)).fit()

# 给出一份模型报告
print(model.summary())

# 作为期5天的预测，返回预测结果、标准误差、置信区间。
print(model.forecast(steps=5))