设观测时间序列 { Z 1 , Z 2 , Z 3 , . . . , Z n } \{Z_1, Z_2, Z_3, ..., Z_n\} {Z1,Z2,Z3,...,Zn}
平稳过程
序列满足三个条件:
- 均值恒定。 E ( Z t ) = μ E(Z_t)=\mu E(Zt)=μ
- 方差恒定。 V a r ( Z t ) = σ 2 Var(Z_t)=\sigma^2 Var(Zt)=σ2
- 协方差仅与时滞k相关。即 C o v ( Z t , Z t + k ) Cov(Z_t, Z_t+k) Cov(Zt,Zt+k)仅是k的函数
白噪声过程
白噪声是平稳过程的一个特例。即除了满足平稳过程的条件外,还满足如下条件:
弱白噪声:
E
(
ϵ
t
ϵ
s
)
=
0
,
t
̸
=
s
E(\epsilon_t \epsilon_s)=0, t\not= s
E(ϵtϵs)=0,t̸=s
独立同分布白噪声:
ϵ
t
与
ϵ
s
\epsilon_t 与\epsilon_s
ϵt与ϵs相互独立
高斯白噪声:
ϵ
t
\epsilon_t
ϵt服从高斯分布
之所以命名为"白"噪声,是因为噪声谱与白光类似。
差分
差分通常用来获取平稳性,可用来对非平稳序列做一个预处理。
p阶差分
相距一期序列元素值的差值为1阶差分,记为 ∇ Z t = Z t − Z t − 1 \nabla Z_t=Z_t-Z_{t-1} ∇Zt=Zt−Zt−1,对差分后得序列,再进行一次差分则为2阶差分,一次类推,p阶差分 ∇ p Z t = ∇ p − 1 Z t − ∇ p − 1 Z t − 1 = ( 1 − B ) p Z t \nabla^pZ_t=\nabla^{p-1}Z_t-\nabla^{p-1}Z_{t-1}=(1-B)^pZ_t ∇pZt=∇p−1Zt−∇p−1Zt−1=(1−B)pZt其中,B为滞后算子,即 B Z t = Z t − 1 BZ_t=Z_{t-1} BZt=Zt−1。
k步差分(季节差分)
相距k期序列元素值得差值,记为 ∇ k Z t = Z t − Z t − k = ( 1 − B k ) Z t \nabla_kZ_t=Z_t-Z_{t-k}=(1-B^k)Z_t ∇kZt=Zt−Zt−k=(1−Bk)Zt
自相关函数(ACF)
变量 Z t 和 Z t + k Z_t和Z_{t+k} Zt和Zt+k的皮尔逊相关系数。即 ρ k = C o v ( Z t , Z t + k ) V a r ( Z t ) V a r ( Z t + k ) \rho_k=\frac{Cov(Z_t, Z_{t+k})}{\sqrt{Var(Z_t)Var(Z_{t+k})}} ρk=Var(Zt)Var(Zt+k)Cov(Zt,Zt+k)
偏自相关函数(PACF)
变量 Z t 和 Z t + k Z_t和Z_{t+k} Zt和Zt+k除去共同线性干预变量 { Z t + 1 , . . . , Z t + k − 1 } \{Z_{t+1},...,Z_{t+k-1}\} {Zt+1,...,Zt+k−1}依赖后的自相关函数。即 P k = C o v ( Z t − Z ^ t , Z t + k − Z ^ t + k ) V a r ( Z t − Z ^ t ) V a r ( Z t + k − Z ^ t + k ) Z ^ t = ∑ i = 1 k − 1 α i ∗ Z t + i Z ^ t + k = ∑ i = 1 k − 1 β i ∗ Z t + k − i P_k = \frac{Cov(Z_t - \hat Z_t, Z_{t+k} - \hat Z_{t+k})}{\sqrt{Var(Z_t - \hat Z_t)Var(Z_{t+k} - \hat Z_{t+k})}} \\ \hat Z_t = \sum_{i=1}^{k-1}{\alpha_i*Z_{t+i}} \\ \hat Z_{t+k} = \sum_{i=1}^{k-1}{\beta_i*Z_{t+k-i}} Pk=Var(Zt−Z^t)Var(Zt+k−Z^t+k)Cov(Zt−Z^t,Zt+k−Z^t+k)Z^t=i=1∑k−1αi∗Zt+iZ^t+k=i=1∑k−1βi∗Zt+k−i
时间序列模型分解
时间序列数据可以分解为可预测部分以及随机波动部分。可预测部分可继续分为长期趋势T、季节变动S及循环变动C三个部分。
自回归模型(AR)
自回归模型对可预测部分进行建模,即对过去p个项进行线性回归:
Z
t
=
C
+
a
1
Z
t
−
1
+
a
2
Z
t
−
2
+
.
.
.
+
a
p
Z
t
−
p
Z_t = C+ a_1Z_{t-1}+ a_2Z_{t-2}+...+a_pZ_{t-p}
Zt=C+a1Zt−1+a2Zt−2+...+apZt−p或表示为
(
1
−
a
1
B
−
a
2
B
2
−
.
.
.
,
−
a
p
B
p
)
Z
t
=
C
(1- a_1B-a_2B^2-...,-a_pB^p)Z_t=C
(1−a1B−a2B2−...,−apBp)Zt=C
移动平均模型(MR)
移动平均模型对随机波动部分进行建模,即对过去q个预测误差项
ϵ
t
\epsilon_t
ϵt(已知)进行线性回归:
Z
t
=
μ
+
ϵ
t
+
θ
1
ϵ
t
−
1
+
θ
2
ϵ
t
−
2
+
.
.
.
+
θ
q
ϵ
t
−
q
Z_t = \mu+\epsilon_t+ \theta_1\epsilon_{t-1}+ \theta_2\epsilon_{t-2}+...+ \theta_q\epsilon_{t-q}
Zt=μ+ϵt+θ1ϵt−1+θ2ϵt−2+...+θqϵt−q或表示为:
Z
t
=
μ
+
(
1
+
θ
1
B
+
θ
2
B
2
+
.
.
.
+
θ
q
B
q
)
ϵ
t
Z_t = \mu + (1+ \theta_1B+ \theta_2B^2+...+ \theta_qB^q)\epsilon_t
Zt=μ+(1+θ1B+θ2B2+...+θqBq)ϵt
实际上,
ϵ
t
\epsilon_t
ϵt是历史实际值与预测值的差值,是白噪声,即服从
W
N
(
0
,
σ
2
)
WN(0,\sigma^2)
WN(0,σ2)分布。这里取名"移动平均",个人感觉不利于对模型的理解,其实是过去q个预测误差项的加权平均值。
自回归移动平均(ARMR)
自回归移动平均模型就是AR和MR的混合模型,非中心化的ARMR模型表示如下:
Z
t
=
μ
+
a
1
Z
t
−
1
+
a
2
Z
t
−
2
+
.
.
.
+
a
p
Z
t
−
p
+
ϵ
t
+
θ
1
ϵ
t
−
1
+
.
.
.
+
θ
q
ϵ
t
−
q
Z_t =\mu + a_1Z_{t-1}+a_2Z_{t-2}+...+a_pZ_{t-p}+\epsilon_t+\theta_1\epsilon_{t-1}+...+\theta_q\epsilon_{t-q}
Zt=μ+a1Zt−1+a2Zt−2+...+apZt−p+ϵt+θ1ϵt−1+...+θqϵt−q或表示为:
a
p
(
B
)
Z
t
=
θ
q
(
B
)
μ
t
a_p(B)Z_t=\theta_q(B)\mu_t
ap(B)Zt=θq(B)μt
该模型适合平稳过程建模,相当于在AR模型基础上引入q个误差项。模型含有p+q+2个参数,可用矩估计以及极大似然估计法估计参数。
自回归求和移动平均(ARIMR)
该模型适合于非平稳过程建模,即序列的均值或方差随时间变动。
齐次非平稳过程 :序列的不同部分的特性非常相似,只不过局部均值水平不同而已。
齐次非平稳过程可通过差分退化为一个平稳过程,即若
{
Z
t
}
\{Z_t\}
{Zt}是其次非平稳过程,那么
{
(
1
−
B
)
d
Z
t
}
\{(1-B)^dZ_t\}
{(1−B)dZt}是平稳过程,其中d是差分的阶。则可表示如下:
a
p
(
B
)
(
1
−
B
)
d
Z
t
=
μ
0
+
θ
q
(
B
)
ϵ
t
a_p(B)(1-B)^dZ_t=\mu_0 + \theta_q(B)\epsilon_t
ap(B)(1−B)dZt=μ0+θq(B)ϵt当d=0时,模型退化为ARMA模型,
μ
0
\mu_0
μ0与过程均值相关;当d>0时,
μ
0
\mu_0
μ0常常可以忽略。
季节性ARIMA模型(SARIMA)
考虑季节周期的影响,需要在ARIMA模型中引入季节差分,例如下文示例中,当前11月的客流量不仅与10月相关,还与去年的12月相关。因此我们需要考虑两个相关性,即周期内相关性(
Z
t
,
Z
t
−
1
,
.
.
.
,
Z
t
−
p
Z_t, Z_{t-1},...,Z_{t-p}
Zt,Zt−1,...,Zt−p)和周期间相关性(
Z
t
,
Z
t
−
s
,
.
.
.
,
Z
t
−
p
s
Z_t, Z_{t-s},...,Z_{t-ps}
Zt,Zt−s,...,Zt−ps,s为季节周期),则模型可表示如下:
A
P
(
B
s
)
a
p
(
B
)
(
1
−
B
s
)
D
(
1
−
B
)
d
Z
t
=
μ
0
+
Θ
Q
(
B
s
)
θ
q
(
B
)
ϵ
t
A_P(B^s)a_p(B)(1-B^s)^D(1-B)^dZ_t=\mu_0 +\Theta_Q(B^s) \theta_q(B)\epsilon_t
AP(Bs)ap(B)(1−Bs)D(1−B)dZt=μ0+ΘQ(Bs)θq(B)ϵt其中,
P
、
Q
、
D
P、Q、D
P、Q、D分别指周期间元素(相距s期)间的自回归项数、预测误差项数及差分阶数。
采用ARIMR建模步骤
- 可视化数据
- 使用ADF/KPSS判定平稳性,决定差分阶数d
- 使用ACF/PACF决定回归项数p及移动平均阶数q
- 尝试减小p、q,优化模型
一般原则
通常,差分消除正相关,但过度差分会引入负相关。
- 观察ACF,正的滞后项多,增加差分阶数d
- 若滞后第一项ACF为负,则可能差分过多
- 差分后,序列标准差最小,则最优
- 若原序列平稳,一阶差分则意味着原序列有一个为常数的平均趋势。 二阶差分则意味着原序列有一个依时间变化的趋势。
- 不差分建模,通常需要常数项;一阶差分建模,仅在平均趋势非零时加常数项;二阶差分无需常数项
- 观察PACF,滞后一项为正,且较大则表明差分不足,增加p观察PACF,滞后一项为正,且较大则表明差分不足,增加p
- 自回归项和预测误差项有可能会互相抵消,p、q同时过大可能会造成过拟合
- 如果自回归项的系数和接近1,那么这时候应该将自回归项减少一个, 同时增加一次差分操作
- 如果预测误差项的系数和接近1, 么这时候应该将预测误差项减少一个, 同时减少一次差分操作。
- 自回归或者移动平均部分有单位根通常表现为长期预测不稳定
示例
以international-airline-passengers.csv数据集为例:
可视化
显然,数据具有长期趋势以及季节变动。
差分
进行一阶差分,使用statsmodels.tsa.stattools.adfuller进行ADF检验,得p值0.05,仅刚好通过检验。使用statsmodels.graphics.tsaplots import plot_acf, plot_pacf 画ACF、PACF图观察ACF图以及PACF图:
依据一般原则判定,差分不足。另外具有明显的季节周期性,需要进行12期的季节差分。
进行二阶差分以及12期的季节差分后:
显然,过度差分了,后续可考虑增加1个预测误差项。
另外,抽取周期间数据,同样画ACF图以及PACF图,决定季节性自回归项数P、预测误差项数Q以及差分阶数D。
代码示例
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import matplotlib.pyplot as plt
import pandas as pd
df = pd.read_csv(r'international-airline-passengers.csv', skipfooter=0)
df.columns = ['time', 'data']
df.dropna(inplace=True)
df.index = pd.date_range(start='1949-01-31', end='1960-12-31', freq='M', )
y = df.data
model = SARIMAX(y[:120], order=(2, 2, 1), seasonal_order=(2, 1, 1, 12), freq='M').fit()
# print(model.summary())
y_pred = model.predict(start='1949-03-31', end='1960-12-31')
plt.plot(y)
plt.plot(y_pred)
plt.legend(['y_true', 'y_pred'])
plt.show()
参考资料
- 《时间序列分析》魏武雄
- 《Keras快速上手:基于Python的深度学习实战》谢梁,鲁颖,劳虹岚
注:如有不当之处,请指正。