ARIMA模型

ARIMA(p,d,q)模型全称为：差分自回归移动平均模型

AR是自回归，p是自回归项；MA是移动平均，q是移动平均项；d是时间序列成为平稳时所做的差分次数

原理：将非平稳时间序列转化为平稳时间序列，然后将因变量仅对它的滞后值（Lag，阶数）以及随机误差项的现值和滞后值进行回归所建立的模型

github示例

平稳性（I）

平稳性就是要求经由样本时间序列所得到的拟合曲线在未来的一段期间内仍能顺着现有的形态“惯性”地延续下去

平稳性要求序列的均值和方差不发生明显变化

严平稳：严平稳表示的分布不随时间的改变而改变，如：白噪音（正态），无论怎么取，都是期望为0，方差为1

弱平稳：期望与相关系数（依赖性）不变，未来某时刻的t的值$X_t$就要依赖它的过去信息，所以需要依赖性

差分法：时间序列在t与t-1时刻的差值

使用差分法可以是原数据变得更平稳，通常用一阶差分或二阶差分即可

# 参数1表示行间隔为1，即：下一行-当前行
# 一阶差分
df['diff_1'] = df['quantity'].diff(1)
# 二阶差分
df['diff_2'] = df['diff_1'].diff(1)
df.plot(subplots=True, figsize=(18,12))

自回归模型（AR）

描述当前值与历史值之间的关系，用变量自身的历史时间数据对自身进行预测

自回归模型必须满足平稳性的要求

必须具有自相关性，如果自相关性系数小于0.5，则不宜采用

自回归只适用于预测与自身前期相关的现象

p阶自回归过程的公式定义：$y_t = \mu + \displaystyle\sum_{i=1}^p \gamma_i y_{t-i} + \epsilon_t$

$y_t$是当前值、$\mu$是常数项、p是阶数、$\gamma_i$是自相关系数、$\epsilon_t$是误差

移动平均模型（MA）

移动平均模型关注的是自回归模型中的误差项的累加

q阶自回归过程的公式定义：$y_t = \mu + \epsilon_t + \displaystyle\sum_{i=1}^q \theta_i \epsilon_{t-i}$

移动平均法能有效地消除预测中的随机波动

自回归移动平均模型（ARMA）

自回归与移动平均的结合

公式定义：$y_t = \mu + \displaystyle\sum_{i=1}^p \gamma_i y_{t-i} + \epsilon_t + \displaystyle\sum_{i=1}^q \theta_i \epsilon_{t-i}$

自相关函数ACF

ACF：autocorrelation function

有序的随机变量序列与其自身相比较，自相关函数反映了同一序列在不同时序的取值之间的相关性

公式：$ACF(k)=p_k= \frac{Cov(y_t,y_{t-k})}{Var(y_t)}$

$p_k$的取值范围为：[-1,1]

ACF图的阴影区间就是置信区间，95%置信区间表示95%的数据在正态分布范围内

偏自相关函数（PACF）

PACF：partial autocorrelation function

对于一个平稳AR（p）模型，求出滞后k自相关系数p(k)时实际上得到的并不是x(t)与x(t-k)之间单纯的相关关系

x(t)同时还会受到中间k-1个随机变量x(t-1)、x(t-2)、...、x(t-k+1)的影响，而这k-1个随机变量又都和x(t-k)具有相关关系

所以自相关系数p(k)里面实际掺杂了其他变量对x(t)与x(t-k)的影响

PACF剔除了中间k-1个随机变量x(t-1)、x(t-2)、...、x(t-k+1)的干扰之后，x(t-k)对x(t)的影响相关程度

即：ACF还包含了其他变量的影响，PACF是严格的这两个变量之间的相关性

根据ACF、PACF确定p、q

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12, 8))
# 211表示2行1列，第一个子图
# ax是axes的缩写，表示一个坐标轴，即：子图所在的位置
ax1 = fig.add_subplot(211)
plot_acf(data['quantity'], lags=8, ax=ax1)
ax1.xaxis.set_ticks_position('bottom')
fig.tight_layout()

# 212表示2行1列，第二个子图
ax2 = fig.add_subplot(212)
plot_pacf(data['quantity'], lags=8, ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
fig.tight_layout()

plt.show()

根据AIC、BIC确定p、q

选择更简单的模型，值越小越好

AIC：赤池信息准则（Akaike Information Criterion, AIC）

AIC = 2k - 2ln(L)

BIC：贝叶斯信息准则（Bayesian Information Criterion, BIC）

AIC = kln(n) - 2ln(L)

k为模型参数个数、n为样本数量、L为似然函数

import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import itertools
from math import inf


def aic_bic(data, flag='BIC', p_min=0, p_max=3, d_min=0, d_max=2, q_min=0, q_max=3, seasonal_order=(0, 0, 0, 0)):
    min_res = {
        "min": inf,
        "pq": (0, 0)
    }
    df = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min, p_max + 1)],
                      columns=['MA{}'.format(i) for i in range(q_min, q_max + 1)]
                      )
    for p, d, q in itertools.product(range(p_min, p_max + 1), range(d_min, d_max + 1), range(q_min, q_max + 1)):
        if p == 0 and d == 0 and q == 0:
            df.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
            continue
        try:
            model = SARIMAX(data, order=(p, d, q), seasonal_order=seasonal_order)
            results = model.fit()
            if flag == 'BIC':
                df.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
                temp = results.bic
            else:
                df.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.aic
                temp = results.aic
            if temp < min_res['min']:
                min_res['min'] = temp
                min_res['pq'] = (p, q)
        except BaseException as e:
            continue
    df = df[df.columns].astype(float)
    fig, ax = plt.subplots(figsize=(10, 8))
    ax = sns.heatmap(df, mask=df.isnull(), ax=ax, annot=True, fmt=".2f")
    ax.set_title(flag)
    plt.show()
    return min_res

# df = xxx
res = aic_bic(df['quantity'], flag='BIC', seasonal_order=(1, 0, 1, 12))
print(res)

残差标准

ARIMA模型的残差是否是平均值为0且方差为常数的正态分布

QQ图：线性即为正态分布

from statsmodels.tsa.statespace.sarimax import SARIMAX
import matplotlib.pyplot as plt

model = SARIMAX(df['quantity'], order=(0, 2, 2), seasonal_order=(1, 0, 2, 12))
model_results = model.fit()
model_results.plot_diagnostics(figsize=(16, 12))
plt.show()

比较预测值误差

from statsmodels.tsa.statespace.sarimax import SARIMAX, SARIMAXResultsWrapper
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

model = SARIMAX(df['quantity'][0:-2], order=(1, 2, 2))
model_results = model.fit()  # type: SARIMAXResultsWrapper
res = model_results.forecast(steps=2)
# 计算均方根误差，越小越好
# 小于1：预测精度高，1到10：预测精度中等，大于10：预测精度可能较低
# 然而，这些阈值仅供参考，具体的阈值可能因数据的特性而异。比较 RMSE 的大小时，最好将其与你的实际数据范围相对比，以更好地理解误差的大小
# 此外，还可以通过其他指标（如平均绝对误差 MAE、平均绝对百分比误差 MAPE 等）来综合考虑模型的性能。在选择模型时，不仅仅要关注 RMSE，还要考虑其他指标，以便更全面地评估模型的表现
rmse = mean_squared_error(df['quantity'][-2:], res, squared=False)

# 绘制实际值和预测值的时间序列图
plt.plot(df['quantity'][0:-2], label='train_data')
plt.plot(df['quantity'][-2:], label='test_data')
plt.plot(res, label='predict_data', color='red')
plt.legend()
plt.show()