一阶差分：时间序列在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图后观察两者构成的点分布在一条直线上，就证明样本数据与正态分布存在线性相关性，即服从正态分布。
NOTE：前两步骤主要是研究严谨性所需。
自回归模型，构建历史数据到未来数据的映射关系
差分得到平稳序列
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项来消除预测时的随机波动
Σ i = 1 q 阶 数 θ i ϵ t − i \Sigma^{q阶数}_{i=1}\theta_i\epsilon_{t-i} Σi=1q阶数θiϵt−i
自回归移动平均模型
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模型是差分自回归移动平均模型
（Autoregressive Integrated Moving Average Model）
比ARMA模型多了个开始的差分处理
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。
partial autocorrelation function
赤化信息准则 Akaike Information Criterion 人名命名
A
I
C
=
2
k
−
2
l
n
L
AIC=2k-2lnL
AIC=2k−2lnL
贝叶斯信息准则 Bayesian Information Criterion
B
I
C
=
k
l
n
N
−
2
l
n
L
BIC=klnN-2lnL
BIC=klnN−2lnL
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的正态分布
差分使用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。
算法的代码还是以基础代码实现为好，使用库函数则参数难以调整。