【机器学习笔记39】时序分析（非平稳序列建模）_时序波形拟合算法

【参考资料】
【1】应用时序分析
【2】https://blog.csdn.net/u010414589/article/details/49622625
【3】http://www.statsmodels.org/devel /examples/notebooks/generated/tsa_arma_0.html
【4】https://www.jianshu.com/p/5c5a6293640a
【5】http://www.docin.com/p-241451370.html

1. 概述

---

1.1 非平稳序列的四大因素

长期趋势: 该因素的影响导致序列呈现出明显的长期趋势（递增、递减等）
循环波动: 该因素的影响导致序列从高到底再从底到高得反复循环波动
季节性变化: 该因素的影响导致序列呈现与季节相关的稳定周期变换
随机波动: 除了上述三种因素外，序列所呈现的随机性现象

2. 趋势分析

---

2.1 拟合法

拟合法是直接采用某已模型对一个明显的时序进行拟合，通常方法如下:

2.2 平滑法

平滑法利用修匀技术削弱短期随机波动序列的影响，类似于图像打马赛克。主要包括如下方法
a. 移动平滑法(n期中心平均):

x ^ t = { 1 n ( x t − ( n − 1 ) / 2 + x t − ( n − 1 ) / 2 + 1 + ⋯ + x t + ⋯ + x t − ( n + 1 ) / 2 ) , x i s o d d 1 n ( 1 2 x t − n / 2 + x t − n / 2 + 1 + ⋯ + x t + x t + n / 2 − 1 + 1 2 x t + n / 2 ) , x i s e v e n \hat{x}_t =

$$\begin{cases} \dfrac{1}{n}(x_{t-(n-1)/2} + x_{t-(n-1)/2+1} + \dots + x_t + \dots + x_{t-(n+1)/2}), & x \quad is \quad odd \\ \\ \dfrac{1}{n}(\dfrac{1}{2}x_{t-n/2} + x_{t-n/2+1}+\dots+x_t + x_{t+n/2-1} + \dfrac{1}{2}x_{t+n/2}), & x \quad is \quad even \end{cases}$$ x^t​=⎩⎪⎪⎪⎨⎪⎪⎪⎧​n1​(xt−(n−1)/2​+xt−(n−1)/2+1​+⋯+xt​+⋯+xt−(n+1)/2​),n1​(21​xt−n/2​+xt−n/2+1​+⋯+xt​+xt+n/2−1​+21​xt+n/2​),​xisoddxiseven​

b. 移动平滑法(n期移动平均):
${\hat{x}}_t=\frac{1}{n}(x_t+x_{t-1}+\cdots +x_{t-n+1})$

c. 简单指数平滑:
${\hat{x}}_t=a{x}_t+a(1-a)x_{t-1}+a(1-a{)}^2{x}_{t-2}+\dots$
其中a为平滑系数 0 < a < 1

d. Holt两参数指数平滑:
{ x ^ t = a x t + ( 1 − a ) ( x ^ t − 1 + r t − 1 ) r t = γ ( x ^ t − x ^ t − 1 ) + ( 1 − γ ) r t − 1

$$\begin{cases} \hat{x}_t=ax_t + (1-a)(\hat{x}_{t-1} + r_{t-1})\\ r_t = \gamma(\hat{x}_t - \hat{x}_{t-1}) + (1 - \gamma)r_{t-1} \end{cases}$$ {x^t​=axt​+(1−a)(x^t−1​+rt−1​)rt​=γ(x^t​−x^t−1​)+(1−γ)rt−1​​

3. 季节效应分析

---

季节效应分析针对存在周期性变化（季节是周期的一种）的时序。例如气温的波动由两个因素组成，一个是季节效应，一个是随机波动。

举例：

每个月的季节指数$S_1,S_2,...,S_{12}$
第i年的第j个月气温定义为: $x_{ij}=\bar{x}.S_j+I_{ij}$，其中$\bar{x}$是各月总的平均温度,$S_j$是第j个月的季节指数，$I_{ij}$是该月的随机波动。计算季节指数步骤如下:

假设m为一周期，共有n个周期

第一步: 计算周期内平均数（诸如所有6月的气温均值）
${\bar{x}}_k=\frac{\sum _{i=1}^n{x}_{ik}}{n}k=1,2,...,m$

第二步: 计算总平均数（所有月份气温均值）
$\bar{x}=\frac{\sum _{i=1}^n\sum _{k=1}^m{x}_{ik}}{nm}$

第二步: 计算各周期的季节指数（各个月气温相对于平均气温的偏差度）
$S_k=\frac{{\bar{x}}_k}{\bar{x}}k=1,2,...,m$

4. ARIMA模型

---

预备: 加载数据

df = pd.read_csv('./data/train.csv', sep=',')
df_kpi = df[df['KPI ID'] == '02e99bd4f6cfb33f']
data = df_kpi.loc[1000:1200, ['value']]
1
2
3

第一步: 利用差分将非平稳序列转换为平稳序列

一阶差分: $▽x_t=x_t-x_{t-1}$
二阶差分: 在++一阶差分的结果上++再次差分 ${▽}^2{x}_t=x_t-x_{t-1}$

def _cal_diff(data):
    output = [0 for col in range(len(data) - 1)] 
    for i in range(len(data) - 1):
        output[i] = data[i + 1] - data[i]
    return output


def _show_diff_data():

    plt.figure(figsize=(8, 6), facecolor='w')
    plt.subplot(121)
    plt.title('origin data')
    plt.plot(data)

    plt.subplot(122)
    plt.title('diff 1')
    plt.plot(_cal_diff(data.as_matrix().reshape(-1)), color='r')

    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

第二步: 采用ADF单位根检验判断差分后是否为平稳序列

单位根检验是将整个时序模型看成是一个差分方程，其平稳性意味着差分方程的稳定性条件。单位根检验就是检验该差分方程的特征方程的各个特征根是否小于1，若存在单位根则认为不是平稳序列。

我们检查原始数据、一阶差分和二阶差分的ADF结果如下:

data        = df_kpi.loc[1000:1200, ['value']].as_matrix().reshape(-1)
data_diff1  = _cal_diff(data)
data_diff2  = _cal_diff(data_diff1)

print(ts.adfuller(data))
print(ts.adfuller(data_diff1))
print(ts.adfuller(data_diff2))
1
2
3
4
5
6
7

具体结果做如下分析：

1 第一个值为ADF测试值，判断是否小于1%、5%、%10对应的统计量
2 第二个值pValue是否趋近于0
由上述两项可以判断一阶差分数据属于平稳序列

第三步: 确定ARMA的p、q参数

fig = plt.figure(figsize=(10, 6))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(data_diff1,lags=40,ax=ax1) 
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(data_diff1,lags=40,ax=ax2)
plt.show()
1
2
3
4
5
6

几个关键定义：
截尾: ACF或PACF在某一阶之后都为0
拖尾: ACF或PACF在某一阶之后不为0，但呈现趋于0的衰减
AR模型: ACF拖尾、PACF截尾
MA模型: ACF截尾、PACF拖尾
ARMA模型: ACF和PACF均为拖尾

参见笔记《时序分析（平稳序列建模）》

上图判断自相关图ACF存在2阶拖尾，PACF存在4阶拖尾，由于设定p=4、q=2。

第四步: 采用ARMA进行建模及预测

这里的代码重新开始全部用pandas下的API编写

# -*- coding: utf-8 -*-
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

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

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

data = df_kpi.loc[1001:1200, ['value']]
#根据pandas时序分析要求将index装换为date类型
data.index = pd.Index(pd.date_range(start='2018-12-12', periods=200, freq='S'))

dta  = data.diff(1) 

#这里从1开始时为了去掉NaN这个无效值
arma_mod42 = sm.tsa.ARMA(dta[1:],(4,2)).fit(disp=False)

predict_sunspots = arma_mod42.predict('2018-12-12 00:03:20', '2018-12-12 00:03:45', dynamic=True)

fig, ax = plt.subplots(figsize=(8, 6))
ax = dta.ix['2018-12-12 00:00:01':].plot(ax=ax)
predict_sunspots.plot(ax=ax)

plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

遗留一些概念尚比较模糊，待下一次笔记整理

---

1. ADF的原理
2. ACF和PACF与模型联系原理
3. ARMA模型的检验
4. 如何通过序列分解去除周期性和线性趋势