【机器学习笔记38】时序分析（ARMA模型）

【参考资料】
【1】《应用时序分析》
【2】https://pandas.pydata.org/pandas-docs/stable/api.html
【3】测试数据来源:智能运维挑战赛 http://iops.ai/
【4】https://blog.****.net/matrix_laboratory/article/details/53912312

当序列经过预处理被认为是一个平稳非白噪声的时序，那么我们可以通过建立一个线性模型来拟合这个时序。

一 数学方法

1.1 差分运算

1. p阶差分

对原始时序序列$x_t$xt​，记$\Delta x_t = x_t - x_{t-1}$Δxt​=xt​−xt−1​为一阶差分序列
对一阶差分序列$x_t^1$xt1​，记$\Delta x_t^2 = x_t^1 - x_{t-1}^1$Δxt2​=xt1​−xt−11​为二阶差分序列
同理类推，对p-1阶差分序列$x_t^{p-1}$xtp−1​，记$\Delta x_t^p = x_t^{p-1} - x_{t-1}^{p-1}$Δxtp​=xtp−1​−xt−1p−1​为p阶差分序列

2. k步差分

对原始时序序列$x_t$xt​，记$\Delta_k x_t = x_t - x_{t - k}$Δk​xt​=xt​−xt−k​为k步差分，其中k为采样点的间隔步数

3. 采用延迟算子表示

定义B为延迟算子，有$x_{t-1} = B x_t \quad x_{t-2} = B^2 x_t$xt−1​=Bxt​xt−2​=B2xt​依次类推

带入上述p阶差分公式，我们有:
一阶差分: $\Delta x_t = (1 - B) x_t $
二阶差分: $\Delta x_t^2 = (1-B)^2 x_t$Δxt2​=(1−B)2xt​
p阶差分 : $\Delta x_t^p = (1-B)^p x_t$Δxtp​=(1−B)pxt​
这里$(1 - B)^p = \sum\limits_{i=0}^{p}(-1)^iC_n^ix_{t-i}$(1−B)p=i=0∑p​(−1)iCni​xt−i​，其中$C_n^i = \dfrac{n!}{i!(n-i)!}$Cni​=i!(n−i)!n!​

同理我们表述k步差分 $\Delta_k x_t = (1 - B^k)x_t$Δk​xt​=(1−Bk)xt​

1.2 线性差分方程

详见《差分方程基础知识备注》笔记

二 ARMA 模型

2.1 AR 模型

具有如下结构的模型被称为p阶自回归模型，记作AR§
$\begin{cases} x_t = \phi_0 + \phi_1 x_{t-1} + \phi_2 x_{t-2} + ... + \phi_p x_{t-p} + \epsilon_t \\ \phi_p \ne 0 \\ E(\epsilon_t)=0, Var(\epsilon_t) = \sigma_{\epsilon_t}^2, E(\epsilon_t \epsilon_s) = 0 \\ E(x_s \epsilon_t) = 0 \end{cases}$⎩⎪⎪⎪⎨⎪⎪⎪⎧​xt​=ϕ0​+ϕ1​xt−1​+ϕ2​xt−2​+...+ϕp​xt−p​+ϵt​ϕp​̸​=0E(ϵt​)=0,Var(ϵt​)=σϵt​2​,E(ϵt​ϵs​)=0E(xs​ϵt​)=0​

理解: t时刻时序的值等于过去p个时间下时序值的线性和，再加上一个随机值，上面式子给出的三个条件，分表表征：
条件1: 保证是p阶
条件2: 保证$\epsilon$ϵ是白噪声随机变量
条件3: 保证$\epsilon$ϵ与需要拟合的时序$x_s$xs​无关

2.2 MA 模型

具有如下结构的模型被称为q阶移动平均模型，记作MA(q):

$\begin{cases} x_t = u_t + \epsilon_t - \theta_1 \epsilon_{t-1} - ... - \theta_q \epsilon_{t - q} \\ \theta_q \ne 0 \\ E(\epsilon_t)=0, Var(\epsilon_t) = \sigma_{\epsilon_t}^2, E(\epsilon_t \epsilon_s) = 0 \end{cases}$⎩⎪⎨⎪⎧​xt​=ut​+ϵt​−θ1​ϵt−1​−...−θq​ϵt−q​θq​̸​=0E(ϵt​)=0,Var(ϵt​)=σϵt​2​,E(ϵt​ϵs​)=0​
理解:t时刻时序的值等于过去值得均值加上q个时间下随机误差的线性组合。

2.3 ARMA 模型

$\begin{cases} x_t = \phi_0 + \phi_1 x_{t-1} + \phi_2 x_{t-2} + ... + \phi_p x_{t-p} + \epsilon_t - \theta_1 \epsilon_{t-1} - ... - \theta_q \epsilon_{t - q} \\ \phi_p \ne 0 \quad \theta_q \ne 0\\ E(\epsilon_t)=0, Var(\epsilon_t) = \sigma_{\epsilon_t}^2, E(\epsilon_t \epsilon_s) = 0 \\ E(x_s \epsilon_t) = 0 \end{cases}$⎩⎪⎪⎪⎨⎪⎪⎪⎧​xt​=ϕ0​+ϕ1​xt−1​+ϕ2​xt−2​+...+ϕp​xt−p​+ϵt​−θ1​ϵt−1​−...−θq​ϵt−q​ϕp​̸​=0θq​̸​=0E(ϵt​)=0,Var(ϵt​)=σϵt​2​,E(ϵt​ϵs​)=0E(xs​ϵt​)=0​
理解: ARMA模型是AR和MA模型的组合，把原来MA模型中的均值替换为前p个时间点时序值得线性组合。

三 平稳时序建模

建模依据

模型	自相关系数	偏自相关系数
AR§	拖尾	p阶截尾
MA(q)	q阶截尾	拖尾
ARMA(p, q)	拖尾	拖尾

代码实现

网上大部分找到的资料都是直接用statsmodels的库做的，这里根据自己的理解在sklearn上的多元线性回归实现。

ARMA模型本质上可以理解成一个多元线性回归的模型，模型的输出就是某一时刻的时序输出，即ARMA模型等式的左侧；而模型的输入就是其他时序的输出和噪声，即ARMA模型等式的右侧。

以构造一个ARMA(2,2)的多元线性回归模型为例，我们有:
$w_{11}*x_{11}+ w_{12}*x_{12} + w_{13}*x_{13} + w_{14}*x_{14} + w_{15}*x_{15} + w_{16}*x_{16}= y_1$w11​∗x11​+w12​∗x12​+w13​∗x13​+w14​∗x14​+w15​∗x15​+w16​∗x16​=y1​，其中：
$w_{11}$对应$\phi_0$ϕ0​，而$x_{11}$x11​恒为1
$w_{12}$为对应$x_{12}$x12​的参数，$x_{12}$x12​对应$\epsilon_t$ϵt​代码中采用一个0~1直接的随机变量
$x_{13}$x13​对应时刻n[t-2]
$x_{14}$x14​对应时刻n[t-1]
$x_{15}$x15​对应时刻t-2的噪声$\epsilon_{t-2}$ϵt−2​
$x_{16}$x16​对应时刻t-1]的噪声$\epsilon_{t-1}$ϵt−1​
$t_{1}$t1​对应时刻n[t]
根据当前时序序列依次构造完成输入矩阵X和输出向量Y。代码如下:

# -*- coding: utf-8 -*-
import numpy  as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from sklearn import linear_model

df = pd.read_csv('./data/train.csv', sep=',')

df_kpi = df[df['KPI ID'] == '02e99bd4f6cfb33f']

max_value = np.max(df_kpi['value'])
min_value = np.min(df_kpi['value'])

df_kpi['std_value'] = df_kpi.value.apply(lambda x : (x - min_value)/(max_value - min_value))

df_kpi.set_index(df_kpi['timestamp'])


#取100个点的时序
data = df_kpi.loc[0:200, ['value']].as_matrix().reshape(-1)


"""
以ARMA(2,2)为例，相当于做一个多元线性回归,如下:

w10*1 + w11*x11 + w12*x12 + w13*x13 + w14*x14 + w15*x15 = y1

y1: t = 2 时刻输出
w10: 对应参数0
w11: 恒为1
x11: 为随机白噪声
w12: 为x12的参数, 代表t = 0时刻输出
w13: 为x13的参数，代表t = 1时刻输出
w14: 为x14的参数，代表t = 0时刻噪声
w15: 为x15的参数，代表t = 1时刻噪声

"""

#一、构筑X和Y的输入、输出 
y = data[2:98]

x = np.ones((96, 6), dtype = float)

#1. 第一列为常数1不变
#2. 第二列为高斯分布随机变量
for i in range(0, 96):
    x[i, 1] = random.random()
    x[i, 4] = random.random()
    x[i, 5] = random.random()

#3. 第三列为n - 2 时刻输出
x[:, 2] = data[0:96]
#4. 第四列为n - 1 时刻输出
x[:, 3] = data[1:97]


#二、采用多元线性回归
regr = linear_model.LinearRegression()
regr.fit(x ,np.transpose(y))


y_pred = regr.predict(x)

#三、显示结果
plt.figure(figsize=(5, 6), facecolor='w')
plt.subplot(211)
plt.title('origin time series')
plt.plot(y)

plt.subplot(212)
plt.title('regession series')
plt.plot(y_pred, color='r')

plt.show()

尝试采用当前构建的模型来对后100个节点的输出进行预测

y_next = data[102:198]

x_next = np.ones((96, 6), dtype = float)

#1. 第一列为常数1不变
#2. 第二列为高斯分布随机变量
for i in range(0, 96):
    x_next[i, 1] = random.random()
    x_next[i, 4] = random.random()
    x_next[i, 5] = random.random()

#3. 第三列为n - 2 时刻输出
x_next[:, 2] = data[100:196]
#4. 第四列为n - 1 时刻输出
x_next[:, 3] = data[101:197]

y_pred = regr.predict(x_next)

#三、显示结果
plt.figure(figsize=(5, 6), facecolor='w')
plt.subplot(211)
plt.title('origin time series')
plt.plot(y_next)

plt.subplot(212)
plt.title('regession series')
plt.plot(y_pred, color='r')

plt.show()