`
时间序列预测常用的方法有移动平均法,指数平均法,AR,MA,ARMA等,这类方法比较适用于小规模,单变量的预测,比如某门店的销量预测等,但是这些都要求序列是平稳的。
弱平稳序列
如果一个时间序列 x t x_t xt 满足以下两个条件,则它是弱平稳的:
对于所有的时刻
t
t
t ,有
E
[
x
t
]
=
μ
E[x_t]=\mu
E[xt]=μ ,其中
μ
\mu
μ是一个常数。
对于所有的时刻
t
t
t和任意的间隔
k
k
k,
x
t
x_t
xt和
x
t
−
k
x_{t-k}
xt−k 的协方差
σ
(
x
t
,
x
t
−
k
)
=
γ
k
\sigma(x_t,x_{t-k})=\gamma_k
σ(xt,xt−k)=γk ,其中
γ
k
\gamma_k
γk与时间
t
t
t无关,它仅仅依赖于间隔
k
k
k 。特别的,当
k
=
0
k=0
k=0时,这个特性意味着
σ
(
x
t
,
x
t
)
=
γ
k
\sigma(x_t,x_t)=\gamma_k
σ(xt,xt)=γk不随时间变化,等于一个与时间
t
t
t无关的常数
γ
0
\gamma_0
γ0,这称为方差平稳性(stationary in variance)。
平稳序列检验
单位根检验是指检验序列中是否存在单位根,因为存在单位根就是非平稳时间序列了。而迪基-福勒检验(Dickey-Fuller test)和扩展迪基-福勒检验(Augmented Dickey-Fuller test可以测试一个自回归模型是否存在单位根(unit root)。迪基-福勒检验模式是D. A迪基和W. A福勒建立的。
ADF检验的原假设是存在单位根,只要这个统计值是小于1%水平下的数字就可以极显著的拒绝原假设,认为数据平稳。注意,ADF值一般是负的,也有正的,但是它只有小于1%水平下的才能认为是及其显著的拒绝原假设。
ADF检验只适合AR( P)模型。且对方差齐性效果好,对异方差性效果不佳。异方差可用PP检验。
假设有这样一个数据,横坐标是时间,纵坐标是店铺的营销额度,
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
import warnings
from datetime import timedelta,datetime
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import itertools
import json
df = pd.read_csv('dj_ord.csv',index_col='pt')
df = df.sort_values(['pt'])
df.index = pd.to_datetime(df.index)
ord_cnt = df['cnt'].values
plt.plot(ord_cnt)
ord_cnt = df['cnt'].shift(-1) - df['cnt']
ord_cnt = ord_cnt.dropna().values
plt.plot(ord_cnt)
t = adfuller(ord_cnt)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
value | |
Test Statistic Value | -1.85795 |
p-value | 0.352095 |
Lags Used | 21 |
Number of Observations Used | 790 |
Critical Value(1%) | -3.43865 |
Critical Value(5%) | -2.86521 |
Critical Value(10%) | -2.56872 |
可以看到 t-statistic 的值 -1.85795 要大于1% 的值,所以无法拒绝原假设,另外,p-value的值也很大。接着我们做一个差分试试
ord_cnt = df['cnt'].shift(-1) - df['cnt']
ord_cnt = ord_cnt.dropna().values
value | |
Test Statistic Value | -5.35685 |
p-value | 4.1707e-06 |
Lags Used | 21 |
Number of Observations Used | 789 |
Critical Value(1%) | -3.43867 |
Critical Value(5%) | -2.86521 |
Critical Value(10%) | -2.56872 |
看到 t-statistic 的值 -5.35685 要小于1%的值,所以拒绝原假设,另外,p-value的值也很小。这说明差分之后的数据已经是平稳序列了。 |
AR模型
AR模型的思想非常的朴素,认为当前的结果
x
t
x_t
xt与之前的结果相关,简单的来说当前的结果等于历史数据的加权平均,然后再加上当前时刻
t
t
t的一个扰动
ε
t
\varepsilon_t
εt,我们只需要算出各个参数
ϕ
\phi
ϕ就可以了。
x
t
=
ϕ
0
+
ϕ
1
x
t
−
1
+
ϕ
2
x
t
−
2
+
.
.
.
.
.
.
ϕ
p
x
t
−
p
+
ε
t
x_t=\phi_0+\phi_1x_{t-1}+\phi_2x_{t-2}+...... \phi_px_{t-p}+\varepsilon_t
xt=ϕ0+ϕ1xt−1+ϕ2xt−2+......ϕpxt−p+εt
MA模型
MA模型的思想也非常的朴素,认为当前的结果
x
t
x_t
xt与之前的扰动相关,简单的来说就是认为
x
t
x_t
xt主要受过去q期的误差项影响,即当前的结果等于历史扰动的加权平均,然后再加上历史结果的均值
μ
\mu
μ,我们也是只需要算出各个参数
ϕ
\phi
ϕ就可以了。这个也非常好理解,如果咱们的数据是在均值附近不断波动的,是什么让数据波动的呢?就是扰动项,扰动使得数据不断的变化,最终变成当前的数据。
x
t
=
μ
+
ε
t
−
θ
1
ε
t
−
1
−
θ
2
ε
t
−
2
−
.
.
.
.
.
.
θ
p
ε
t
−
q
x_t=\mu+\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-...... \theta_p\varepsilon_{t-q}
xt=μ+εt−θ1εt−1−θ2εt−2−......θpεt−q
ARMA模型
ARMA模型是最常用的平稳序列拟合模型,随机变量Xt的取值xt不仅与以前p期的序列值有 关还与前q期的随机扰动有关。
x
t
=
ϕ
0
+
ϕ
1
x
t
−
1
+
ϕ
2
x
t
−
2
+
.
.
.
.
.
.
ϕ
p
x
t
−
p
+
ε
t
−
θ
1
ε
t
−
1
−
θ
2
ε
t
−
2
−
.
.
.
.
.
.
θ
p
ε
t
−
q
x_t=\phi_0+\phi_1x_{t-1}+\phi_2x_{t-2}+...... \phi_px_{t-p}+\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-...... \theta_p\varepsilon_{t-q}
xt=ϕ0+ϕ1xt−1+ϕ2xt−2+......ϕpxt−p+εt−θ1εt−1−θ2εt−2−......θpεt−q
ARIMA模型定阶
好了,现在讨论一下怎么算出这几个参数,这几个参数必然使得我们预测的误差最小,此时得到的即是最优的模型。通常来讲,akaike information criterion(AIC)和bayesian information criterion(BIC)是评价模型优良的两个指标。这两种评价指标不仅适用于事件序列模型,还广泛广泛运用于其他数学模型中。
A
I
C
=
−
2
l
n
(
L
)
+
2
k
AIC = -2 ln(L) + 2 k
AIC=−2ln(L)+2k
B
I
C
=
−
2
l
n
(
L
)
+
l
n
(
n
)
∗
k
BIC = -2 ln(L) + ln(n)*k
BIC=−2ln(L)+ln(n)∗k
让我们来理解AIC的含义,AIC由两部分组成,一部分是对数极大似然函数,另一部分则是参数的个数。极大似然函数是评价模型拟合优劣性的指标,值越大说明拟合的效果越好。然而使用过多的参数可以拟合的很好却会出现过度拟合的情况,这样的模型泛化能力很差,因此加上参数的个数实际上是对极大似然函数进行”惩罚“。选取AIC值最小的模型作为最优模型,实质上是平衡了欠拟合和过拟合。
# 创建训练数据
train_start = datetime(2019,10,1)
train_end = datetime(2020,10,1)
train_data = df[train_start:train_end]
#定阶
warnings.filterwarnings("ignore") # specify to ignore warning messages
pmax,qmax,dmax,smax = 8 ,8 ,2 ,7
# p是AR阶数,q是MA阶数, d是差分次数, s是周期
pdq = itertools.product(range(1,pmax+1),range(1,dmax+1),range(1,qmax+1),[smax])
aic_matrix = {} #aic矩阵
for p,d,q,s in pdq:
try:
model = SARIMAX(endog=train_data['cnt'],order=(p,d,q,s))
results = model.fit(disp=True)
aic_matrix['{},{},{}'.format(p,d,q)] = results.aic
print('ARIMA p:{} d:{} q:{} - AIC:{}'.format(p, d, q, results.aic))
except:
continue
p,d,q = sorted(aic_matrix.items(), key=lambda x: x[1])[0][0].split(',')
print('AIC最小的p : {}, d : {}, q : {}'.format(p,d,q),end='\r')
- ARIMA p:8 d:1 q:8 - AIC:9138.284284512645
- ARIMA p:8 d:2 q:1 - AIC:9165.776270275972
- ARIMA p:8 d:2 q:2 - AIC:9166.238801899044
- ARIMA p:8 d:2 q:3 - AIC:9166.731346671577
- ARIMA p:8 d:2 q:4 - AIC:9169.846039297874
- ARIMA p:8 d:2 q:5 - AIC:9167.21095209448
- ARIMA p:8 d:2 q:6 - AIC:9145.26033882199
- ARIMA p:8 d:2 q:7 - AIC:9129.24423753544
- ARIMA p:8 d:2 q:8 - AIC:9132.197927380188
- AIC最小的p : 6, d : 2, q : 7
模型训练
最优的参数求出来之后我们就可以训练一个模型进行预测了
model = SARIMAX(endog=train_data['cnt'], order=(6,2,7,7)) #建立ARIMA(7, 1,7)模型
result_ARIMA = model.fit(disp=-1)
predict=result_ARIMA.predict()
# 这里直接包含了差分,所以不需要咱们再自己手动相加了
# 预测的值(差分) + 上一时刻的值
plt.figure(figsize=(18,5),facecolor='white')
predict.plot(color='blue', label='Predict')
train_data['cnt'].plot(color='red', label='Original')
result_ARIMA.summary()
由SARIMAX的输出产生的SARIMAX返回大量的信息,但是我们将注意力集中在系数表上。 coef列显示每个特征的重量(即重要性)以及每个特征如何影响时间序列。 P>|z| 列通知我们每个特征重量的意义。 这里,系数的p值都低于0.05时是显著的 ,所以在我们的模型中保留小于0.05的权重特征。
在适合季节性ARIMA模型(以及任何其他模型)的情况下,运行模型诊断是非常重要的,以确保没有违反模型的假设。 plot_diagnostics对象允许我们快速生成模型诊断并调查任何异常行为。
result_ARIMA.plot_diagnostics(figsize=(16, 12));
左上图是训练数据的残差,右上图可以看到绿色的线是正态分布,橘色的线是残差的分布,和正态分布类似的钟形曲线,说明参数是随机序列(平稳),左下角是QQ图。
模型预测
pred_days = 60
test_start = train_end
test_end = train_end + timedelta(pred_days)
test_data = df[test_start:test_end]
y_forecasted =result_ARIMA.forecast(steps=pred_days, alpha=0.01) #作为期7天的预测
plt.figure(figsize=(18,5),facecolor='white')
y_forecasted.plot(color='blue', label='Predict')
test_data['cnt'].plot(color='red', label='Original')
做过时间序列的朋友可能常常会有这样的感受,用了某种算法做出来的测试集的平均绝对误差率或者r2系数都很好,但是把测试集的真实值及预测值画出来对比一下,就会发现t时刻的预测值往往是t-1时刻的真实值,也就是模型倾向于把上一时刻的真实值作为下一时刻的预测值,导致两条曲线存在滞后性,也就是真实值曲线滞后于预测值曲线,就像下图右边所显示的那样。之所以会这样,是因为序列存在自相关性,如一阶自相关指的是当前时刻的值与其自身前一时刻值之间的相关性。因此,如果一个序列存在一阶自相关,模型学到的就是一阶相关性。而消除自相关性的办法就是进行差分运算,也就是我们可以将当前时刻与前一时刻的差值作为我们的回归目标。但是,在其他任务进行特征选择的时候,我们是会把目标变量相关性低的特征去掉,留下相关性强的特征。
还有一点需要注意的是,单纯使用平均绝对误差率或者r2系数容易误导,因为即使指标效果很好,但是很有可能这个模型也是没有用的。一种做法是可以计算一个基准值,即如果全部预测值都采用上一时刻的真实值,这时候的平均绝对误差率或者r2系数是多少,如果你以后加了其他特征,依然没办法超过这个基准值或者提升不大,那就放弃吧,这个时间序列可能已经没办法预测了。
https://www.zhihu.com/collection/282158900
写给你的金融时间序列分析:基础篇
探索数据之自相关图分析(一)
时间序列分析导论书摘:自相关图意义分析
携程如何基于ARIMA时序分析做业务量的预测 - 携程技术的文章 - 知乎