基本概念
一阶差分:时间序列在t 与t-1 时刻函数值的差值,提升时序数据的平稳性(ARIMA算法对数据平稳性有要求)
二阶差分:在一阶差分的基础上再做一次(一般时序数据最多做两阶,再多则预测意义不大)
自回归模型: f ( t ) , t ∈ [ 0 , m ] f(t),t\in[0,m] f(t),t∈[0,m]已知历史数据,预测 f ( t ′ ) , t ′ > m f(t'),t'>m f(t′),t′>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图后观察两者构成的点分布在一条直线上,就证明样本数据与正态分布存在线性相关性,即服从正态分布。
算法步骤
- ADF检验平稳性,分析t值,确定是否可以显著拒绝序列不平稳假设 ( p < 0.05 或 0.01 ) (p<0.05或0.01) (p<0.05或0.01)
- 偏自相关分析、自相关分析,根据截尾情况估算p、q值(使用AIC BIC选择简单的模型)
- 检验模型残差是否为白噪声(纯随机),Q统计量的p值检验 或 信息准则AIC和BIC 或 ACF/PACF图分析。
- 模型参数表→模型公式
- 向后预测输出指定阶数(向后预测几个时间点的值)
NOTE:前两步骤主要是研究严谨性所需。
AR模型
自回归模型,构建历史数据到未来数据的映射关系
差分得到平稳序列
p阶自回归过程公式定义:
y t = μ 常 数 + Σ i = 1 p 阶 数 γ i y t − i + ϵ t 误 差 y_t=\mu常数+\Sigma^{p阶数}_{i=1}\gamma_iy_{t-i}+\epsilon_t误差 yt=μ常数+Σi=1p阶数γiyt−i+ϵt误差
γ
\gamma
γ是自相关系数,表示i时刻和t时刻的相关性,相关性越大则系数越大,i时刻的y值对t时刻的y值影响越大。
p
p
p阶数表示t时刻与前面p个时刻的y值相关,可通过计算自相关系数确定p的大小,若i时刻自相关系数<0.5,则p就截止到i+1时刻。
MA模型
移动平均模型
自回归模型中误差项会累加,所以预测的模型公式中加入MA项来消除预测时的随机波动
Σ i = 1 q 阶 数 θ i ϵ t − i \Sigma^{q阶数}_{i=1}\theta_i\epsilon_{t-i} Σi=1q阶数θiϵt−i
ARMA模型
自回归移动平均模型
y t = μ 常 数 + Σ i = 1 p 阶 数 γ i y t − i + ϵ t 误 差 + Σ i = 1 q 阶 数 θ i ϵ t − i y_t=\mu常数+\Sigma^{p阶数}_{i=1}\gamma_iy_{t-i}+\epsilon_t误差+\Sigma^{q阶数}_{i=1}\theta_i\epsilon_{t-i} yt=μ常数+Σi=1p阶数γiyt−i+ϵt误差+Σi=1q阶数θiϵt−i
(p,d,q)都需要指定,pq是阶数,也称为滞后值,d=1或2(通常取1),对应I字母,所以:
ARIMA模型
ARIMA模型是差分自回归移动平均模型
(Autoregressive Integrated Moving Average Model)
比ARMA模型多了个开始的差分处理
ACF自相关函数
autocorrelation function
反映序列在不同时间点上的取值相关性
A
C
F
(
k
)
=
ρ
k
=
C
o
v
(
y
t
,
y
t
−
k
)
/
V
a
r
(
y
t
)
∈
[
−
1
,
1
]
ACF(k)=\rho_k=Cov(y_t,y_{t-k})/Var(y_t)\in[-1,1]
ACF(k)=ρk=Cov(yt,yt−k)/Var(yt)∈[−1,1]
-1负相关 +1 正相关 0 不相关
可以绘制
A
C
F
(
k
)
ACF(k)
ACF(k)的曲线图
NOTE:图来自网络博客 https://blog.csdn.net/y1163307648/article/details/107832178。
PACF偏/部分自相关函数
partial autocorrelation function
- ACF得到的并不是 y t , y t − k y_t,y_{t-k} yt,yt−k间的严格关系,因为 y t y_t yt还受前k-1个y的影响,即ACF包括直接和间接的相关性,所以PACF是剔除干扰后的相关性计算。
- PACF找的是残差 R t R_t Rt与前第k个滞后值 y t − k y_{t-k} yt−k的相关性。只描述观测值 y t y_t yt和其滞后项 y t − k y_{t-k} yt−k的直接关系
- AR过程:在ACF图中显示出逐渐减少的趋势,因为作为一个AR过程,其当前与过去的滞后项具有良好的相关性;PACF在滞后项阶数后会急剧下降,因为这些接近当前项的滞后项可以很好地捕获变化,因此不需要很多过去的滞后项来预测当前项
- MA过程:在ACF 图中在阶数q 之后迅速下降,能够画出相邻的滞后项之间的良好的相关关系,且和过去的滞后项没有相关关系;在PACF中,曲线逐渐下降,临近的滞后项不能预测当前项,更远的滞后项有良好的相关关系。
- 计算公式复杂,参考博客:点击查看。
pq确定表
AIC/BIC/HQ模型选择
AIC
赤化信息准则 Akaike Information Criterion 人名命名
A
I
C
=
2
k
−
2
l
n
L
AIC=2k-2lnL
AIC=2k−2lnL
BIC
贝叶斯信息准则 Bayesian Information Criterion
B
I
C
=
k
l
n
N
−
2
l
n
L
BIC=klnN-2lnL
BIC=klnN−2lnL
HQ
Hannan-Quinn Criterion
H
Q
=
k
l
n
(
l
n
N
)
−
2
l
n
L
HQ=kln(lnN)-2lnL
HQ=kln(lnN)−2lnL
k:模型参数的个数,模型参数的差别就在于p和q,所以可认为k=p+q
N:样本数量
L:模型函数的(最大)似然函数值
NOTE:比较时只采用其中的一种原则;均满足参数越多越不好的哲学原理;原则的值越小则说明模型越简单。
模型残差检验
检验残差是否是均值=0 方差=constant的正态分布
- 绘制残差图,定量计算均值和方差,并根据样本分布使用假设检验理论来计算是否符合正态分布
- 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。
备注
算法的代码还是以基础代码实现为好,使用库函数则参数难以调整。
- 只实现了AR模型的预测,因为MA模型的θ相关性系数不清楚如何计算
- AR模型预测效果在d=0或1的时候较稳定(指的是De_Diff后),在d=2的时候差分序列预测结果较好,但逆差分还原后到一阶差分结果时预测结果一般,还原到零阶差分结果时预测结果很差,所以d不能过大。实际应用中需要的结果是零阶的结果(即原时序数据)而不是二阶的结果,所以如果逆差分的结果偏差大,则即使未逆差分时的结果拟合程度较好(实测的结果,d=2时两阶差分拟合程度高,但逆差分到零阶后拟合程度却很差),仍无法使用该模型。
- Q-Q图也是基础代码实现,没有调用高级库,未展现在本博客中。思路:横坐标是分位点,纵坐标是分位点处的序列值(归一化后),也可以将多组序列数据绘制到一张Q-Q图中,观察是否相似,如果相似则说明统计上服从的分布接近。根据实际图可以发现和正态分布接近,所以残差检验通过。