【时间序列预测算法】——ARIMA 算法介绍及代码实现

基本概念

一阶差分：时间序列在t 与t-1 时刻函数值的差值，提升时序数据的平稳性（ARIMA算法对数据平稳性有要求）
二阶差分：在一阶差分的基础上再做一次（一般时序数据最多做两阶，再多则预测意义不大）
自回归模型：$f(t)，t\in [0,m]$已知历史数据，预测$f(t^{\prime }),t^{\prime }>m$ 。用变量自身的历史数据对自身进行预测 自回归模型必须满足平稳性的要求
ARIMA：全称“自回归移动平均模型”。记作ARIMA(p,d,q)，用于时序预测的模型，通常适用单列时序数据分析，前提是时序数据平稳（围绕均值有限波动，方差有限，且均值和方差几乎不变），不能有明显上升/下降趋势（如果有，则要差分处理），可使用ADF检验稳定性（待研究）。
严平稳：时序的统计特性不随时间变化而时移（实际中难以满足）
弱平稳：期望值基本不变，相关性基本不变（未来的值与历史值的依赖关系不变，否则无法预测）（实际中通常情况）
截尾：自相关函数（ACF）或偏自相关函数（PACF）在某阶后瞬间衰减为0的性质，可通过计算相邻点的斜率超过预设阈值来自动判定，或通过计算折线图中积分值占据总积分（面积）的比例超过预设比例阈值如90%来自动判定。
拖尾：ACF或PACF在某阶后逐渐衰减为0的性质。
QQ图：quantile-quantile plot，用于检验一组数据是否服从某一分布；检验两个分布是否服从同一分布。原理是用图形的方式比较两个概率分布，把两组数据的分位数放在一起绘图比较——首先选好分位数间隔，图中的每个点（x,y）反映出其中一组数据的分布分位数（y坐标）和相同分位数下另一组数据的值（x坐标）。如果两个分布相似，则该Q-Q图趋近于落在y=x线上。如果两分布线性相关，则点在Q-Q图上趋近于落在一条直线y=kx上，如果要利用Q-Q图来对数据进行正态分布的检验，则可以令一组数据为正态分布的标准数据，另一组数据与之计算分位数，绘QQ图后观察两者构成的点分布在一条直线上，就证明样本数据与正态分布存在线性相关性，即服从正态分布。

算法步骤

1. ADF检验平稳性，分析t值，确定是否可以显著拒绝序列不平稳假设$(p<0.05或0.01)$
2. 偏自相关分析、自相关分析，根据截尾情况估算p、q值（使用AIC BIC选择简单的模型）
3. 检验模型残差是否为白噪声（纯随机），Q统计量的p值检验 或 信息准则AIC和BIC 或 ACF/PACF图分析。
4. 模型参数表→模型公式
5. 向后预测输出指定阶数（向后预测几个时间点的值）

NOTE：前两步骤主要是研究严谨性所需。

AR模型

自回归模型，构建历史数据到未来数据的映射关系
差分得到平稳序列
p阶自回归过程公式定义：

$y_t=\mu 常数+{\Sigma }_{i=1}^{p阶数}{\gamma }_i{y}_{t-i}+{ϵ}_t误差$

$\gamma$是自相关系数，表示i时刻和t时刻的相关性，相关性越大则系数越大，i时刻的y值对t时刻的y值影响越大。
$p$阶数表示t时刻与前面p个时刻的y值相关，可通过计算自相关系数确定p的大小，若i时刻自相关系数<0.5，则p就截止到i+1时刻。

MA模型

移动平均模型
自回归模型中误差项会累加，所以预测的模型公式中加入MA项来消除预测时的随机波动

${\Sigma }_{i=1}^{q阶数}{\theta }_i{ϵ}_{t-i}$

ARMA模型

自回归移动平均模型

$y_t=\mu 常数+{\Sigma }_{i=1}^{p阶数}{\gamma }_i{y}_{t-i}+{ϵ}_t误差+{\Sigma }_{i=1}^{q阶数}{\theta }_i{ϵ}_{t-i}$

(p,d,q)都需要指定，pq是阶数，也称为滞后值，d=1或2（通常取1），对应I字母，所以：

ARIMA模型

ARIMA模型是差分自回归移动平均模型
（Autoregressive Integrated Moving Average Model）
比ARMA模型多了个开始的差分处理

ACF自相关函数

autocorrelation function
反映序列在不同时间点上的取值相关性
$ACF(k)={\rho }_k=Cov(y_t,y_{t-k})/Var(y_t)\in [-1,1]$
-1负相关 +1 正相关 0 不相关
可以绘制$ACF(k)$的曲线图

NOTE:图来自网络博客 https://blog.csdn.net/y1163307648/article/details/107832178。

PACF偏/部分自相关函数

partial autocorrelation function

1. ACF得到的并不是$y_t,y_{t-k}$间的严格关系，因为$y_t$还受前k-1个y的影响，即ACF包括直接和间接的相关性，所以PACF是剔除干扰后的相关性计算。
2. PACF找的是残差$R_t$与前第k个滞后值$y_{t-k}$的相关性。只描述观测值$y_t$和​其滞后项$y_{t-k}$的直接关系
3. AR过程：在ACF图中显示出逐渐减少的趋势，因为作为一个AR过程，其当前与过去的滞后项具有良好的相关性；PACF在滞后项阶数后会急剧下降，因为这些接近当前项的滞后项可以很好地捕获变化，因此不需要很多过去的滞后项来预测当前项
4. MA过程：在ACF 图中在阶数q 之后迅速下降，能够画出相邻的滞后项之间的良好的相关关系，且和过去的滞后项没有相关关系；在PACF中，曲线逐渐下降，临近的滞后项不能预测当前项，更远的滞后项有良好的相关关系。
5. 计算公式复杂，参考博客：点击查看。

pq确定表

AIC/BIC/HQ模型选择

AIC

赤化信息准则 Akaike Information Criterion 人名命名
$AIC=2k-2lnL$

BIC

贝叶斯信息准则 Bayesian Information Criterion
$BIC=klnN-2lnL$

HQ

Hannan-Quinn Criterion
$HQ=kln(lnN)-2lnL$

k:模型参数的个数，模型参数的差别就在于p和q，所以可认为k=p+q
N:样本数量
L:模型函数的(最大)似然函数值

NOTE：比较时只采用其中的一种原则；均满足参数越多越不好的哲学原理；原则的值越小则说明模型越简单。

模型残差检验

检验残差是否是均值=0 方差=constant的正态分布

1. 绘制残差图，定量计算均值和方差，并根据样本分布使用假设检验理论来计算是否符合正态分布
2. QQ图绘制，如果呈现线性则符合正态分布。

代码实现

差分使用pandas：

df.diff(1)#按照间隔1个点计算差值
df.plot(subplots=True)#绘制曲线

使用statsmodels库提供的函数可以计算相关性系数等。

pip install statsmodels

使用基础库计算：

matplotlib.rc('xtick', labelsize=40)
matplotlib.rc('ytick', labelsize=40)
# 使用库进行处理更简单
from statsmodels.tsa.stattools import acf, pacf

# ------------------参数定义------------------
#d=1时 p,d,q=6,1,6 p+q=12
#d=2时 p,d,q=4,2,2 p+q=6<12 根据模型选择原则选择d=2
d = 2
K_MAX = 10
epochs = 10
Period = 100
epsilon = [np.random.normal() for _ in range(Period)]
# ------------------时序数据构造------------------
t = np.linspace(0, 10, 500)
# adding normally distributed series in exponential series
y = np.random.normal(0, 5, 500) + np.exp(t ** 0.5)
# 原始时序数据
plt.figure(figsize=(16, 7))
plt.plot(t, y)
plt.show()


# ------------------差分函数------------------
def diff(y, d):
    for d_index in range(d):
        y = [y[i] - y[i - 1] for i in range(1, len(y))]
    return y


diff_y = diff(y, d)
plt.plot(t[:len(diff_y)], diff_y)
plt.title("Diff of %d order" % d)
plt.show()


# ------------------ACF函数计算------------------
def ACF(y, k):
    y_t = np.array([y[i] for i in range(k, len(y))])
    y_t_k = np.array([y[i] for i in range(k - k, len(y) - k)])
    avg_y_t = np.mean(y_t)
    avg_y_t_k = np.mean(y_t_k)
    # * 相当于 np.multiply对应元素相乘 要求两array的shape一致
    cov = np.mean((y_t - avg_y_t) * (y_t_k - avg_y_t_k))  # 应该是÷(n-1)
    var = np.mean((y_t - avg_y_t) ** 2)  # 应该是÷(n-1)
    return cov / var


acf_list = []
for k in range(K_MAX):
    acf_list.append(ACF(diff_y, k))
plt.plot(acf_list)
plt.title("ACF-K Curve")
plt.show()
for k in range(K_MAX):
    plt.plot(np.ones(2) * k, [0, acf_list[k]])
# 绘制水平辅助线 ax h line → axhline
# 绘制竖直辅助线 ax v line → axvline
plt.axhline(y=0, linestyle='--')
plt.title("ACF-k Point")
plt.show()
# ------------------PACF函数计算------------------
# 参考博客：https://www.cnblogs.com/nkuhyx/p/12162637.html
from scipy.linalg import toeplitz

# 1.使用statsmodels计算：
import statsmodels.tsa.stattools as stattools


def default_pacf(ts, k):
    return statools.PACF(ts, nlags=k, unbiased=True)


# 2.使用原生库计算：
def PACF(ts, k):
    ''' Compute partial autocorrelation coefficients for given time series，unbiased
    '''

    def yule_walker(ts, order):
        ''' Solve yule walker equation
        '''
        x = np.array(ts) - np.mean(ts)
        n = x.shape[0]

        r = np.zeros(order + 1, np.float64)  # to store acf
        r[0] = x.dot(x) / n  # r(0)
        for k in range(1, order + 1):
            r[k] = x[:-k].dot(x[k:]) / (n - k)  # r(k)

        R = toeplitz(r[:-1])

        return np.linalg.solve(R, r[1:])  # solve `Rb = r` to get `b`

    res = [1.]
    for i in range(1, k + 1):
        res.append(yule_walker(ts, i)[-1])
    return res


# pacf_list = []
# for k in range(K_MAX):
#     pacf_list.append()
pacf_list = PACF(diff_y, K_MAX)
plt.plot(pacf_list)
plt.title("PACF-K Curve")
plt.show()
for k in range(K_MAX):
    plt.plot(np.ones(2) * k, [0, pacf_list[k]])
# 绘制水平辅助线 ax h line → axhline
# 绘制竖直辅助线 ax v line → axvline
plt.axhline(y=0, linestyle='--')
plt.title("PACF-k Point")
plt.show()
# --------------截尾自动判别-------------
# 计算绝对值，在>0的范围上进行计算斜率
# 斜率截尾时满足↗↘，转折点就是阶数
acf_list_abs = np.abs(acf_list)
pacf_list_abs = np.abs(pacf_list)
# 计算斜率k = 差值δy÷δx
acf_k = np.abs(diff(acf_list_abs, 1)) / 1
pacf_k = np.abs(diff(pacf_list_abs, 1)) / 1
q = 0
for i in range(1, len(acf_k) - 1):
    if (acf_k[i] > acf_k[i - 1]) and (acf_k[i] > acf_k[i + 1]):
        q = i + 1
        break
p = 0
for i in range(1, len(pacf_k) - 1):
    if (pacf_k[i] > pacf_k[i - 1]) and (pacf_k[i] > pacf_k[i + 1]):
        p = i + 1
        break
print("p:", p,end=" ")
print("d:",d,end=" ")
print("q:", q,end=" ")
# -----------AR模型----------------
def AR(miu,yt,p,epsilon,t):
    # t:yt输入的第0项对应的时刻点
    # yt:y0 y1 y2 -> y3
    # acf_list:k=0 1 2 3
    ar = miu
    for i in range(1,p+1):
        ar += acf_list[i]*yt[-i]+epsilon[(t+i-1)%len(epsilon)]
    return ar
data = []
label = []
predict = []
for i in range(len(diff_y)-p):
    data.append(diff_y[i:i+p])
    label.append(diff_y[i+p])
miu=0
print()
for epoch in range(epochs):
    predict = []
    for i in range(len(data)):
        predict.append(AR(miu,data[i],p,epsilon,i))
    # 计算偏差
    bias_loss = np.mean(-np.array(predict)+np.array(label))
    miu = bias_loss#自动学习 更新模型参数
    loss = np.mean((np.array(predict)-np.array(label))**2)
    print("epoch:",epoch+1,"loss:",loss,"bias_loss:",bias_loss,"miu:",miu)
print("学习后 miu:",miu)
plt.plot(label,label='True')
plt.plot(predict,label='Pred')
plt.legend(loc=1)
plt.title("AR predictions")
plt.show()

AR模型的预测效果：

10个分位点时绘制QQ图

20个分位点时绘制QQ图

50个分位点时绘制QQ图：

NOTE：必过(0,0) (1,1)两点，因为0时对应min值0（首先要归一化 否则关系比较难看出来），1时对应max值1。

备注

算法的代码还是以基础代码实现为好，使用库函数则参数难以调整。

1. 只实现了AR模型的预测，因为MA模型的θ相关性系数不清楚如何计算
2. AR模型预测效果在d=0或1的时候较稳定（指的是De_Diff后），在d=2的时候差分序列预测结果较好，但逆差分还原后到一阶差分结果时预测结果一般，还原到零阶差分结果时预测结果很差，所以d不能过大。实际应用中需要的结果是零阶的结果（即原时序数据）而不是二阶的结果，所以如果逆差分的结果偏差大，则即使未逆差分时的结果拟合程度较好（实测的结果，d=2时两阶差分拟合程度高，但逆差分到零阶后拟合程度却很差），仍无法使用该模型。
3. Q-Q图也是基础代码实现，没有调用高级库，未展现在本博客中。思路：横坐标是分位点，纵坐标是分位点处的序列值（归一化后），也可以将多组序列数据绘制到一张Q-Q图中，观察是否相似，如果相似则说明统计上服从的分布接近。根据实际图可以发现和正态分布接近，所以残差检验通过。