[Python][pmdarima] 季节性ARIMA模型学习

[Python][pmdarima] ARIMA时间序列模型学习

• 背景
• 一、建模步骤
• 二、季节性ARIMA模型
• 1、模块导入及数据读取
• 一、模块及数据
• 二、季节性分析
• 2、平稳性检验
• 一、Augmented Dickey-Fuller Test(ADF检验)
• 二、Python实现
• 3、序列平稳化
• 差分法
• 4、白噪声检验
• 一、白噪声定义
• 二、Ljung-Box检验
• 三、Python实现
• 5、时间序列定阶
• 6、构建ARIMA模型及预测
• 一、构建模型
• 二、模型识别
• 三、需求预测
• 三、结果评估
• 7、小结

背景

前段时间参与了一个快消行业需求预测的项目。其中，用到了移动平均法、ARIMA、Xgboost等方法进行预测，如今打算总结一下ARIMA。
由于项目的销售数据属于私密数据，这里用网上找的一份案例数据用于展现。
构建ARIMA模型能够用到statsmodels库和pmdarima库。我这里用的是pmdarima库，这个库有一个优势是可以自动地对ARIMA模型的参数进行自动肯定。

一、建模步骤

二、季节性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：单个季节期间的时间步数。

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()

由下子图，确实每一年的数据呈现先升后降的特征。

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检验主要是经过判断时间序列中是否含有单位根：若是序列平稳，就不存在单位根；不然，就会存在单位根。

详细的介绍能够看这篇博客
https://blog.csdn.net/FrankieHello/article/details/86766625/

二、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，所以不能拒绝原假设。

3、序列平稳化

差分法

      对于非平稳时间序列，能够经过d阶差分运算使变成平稳序列。从统计学角度来说，只有平稳性的时间序列才能避免“伪回归”的存在，才有经济意义。
      在这里，我先进行了一阶差分，再进行了一次季节性差分，才经过ADF检验。

# 差分处理非平稳序列，先进行一阶差分
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∑m​T−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模型对数据的要求仍是挺高的。其中，数据量太少不行，当时作的时候发现一些门店的数据量太少，并且只有一年左右的跨度，很难发现一些周期规律。另外就是对平稳性的要求。       目前对时间序列分析仍是初学阶段，不少地方的解释不够深刻，后续还须要继续学习、完善。