时间序列分析
时间序列分析大致可分成 三大部分:描述过去、分析规律和预测未来。而根据时间和数值性质不同,可分为两类:
- 时期序列:数值要素反映现象在一定时期内发展的结果
- 时点序列:数值要素反映现象在一定时点上的瞬时水平
时间序列分解:
- 长期趋势T:表现在相当长的一段时间内,收到长期趋势影响因素的影响
- 季节趋势S:季节是广义的,可以按月,季,周为单位,指由于季节转变使指标发生周期性变动
- 循环变动C:通常以若干年为周期,变现为波浪式的周期变动
- 不规则变动I:由于某些随机因素导致的数值变化
数据:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
# 1. 生成时间序列数据
np.random.seed(42)
n_samples = 100
time_steps = np.arange(n_samples)
ts_data = np.random.randn(n_samples).cumsum()
plt.figure()
plt.plot(ts_data)
plt.xlabel('时间')
plt.ylabel('值')
plt.show()
-
自相关系数(Autocorrelation Coefficient): 自相关系数是衡量时间序列在不同滞后阶数下的相关性的统计指标。在时间序列中,可能存在当前时刻和之前某些时刻之间的相关性,自相关系数就是用来度量这种相关性的。自相关系数的范围在 -1 到 1 之间,其中 1 表示完全正相关,-1 表示完全负相关,0 表示无相关性。
自相关系数可以用来检测时间序列数据中的季节性、周期性和其他模式。常用的自相关系数包括一阶自相关系数(ACF(1))、二阶自相关系数(ACF(2))等,表示时间序列与自身在不同滞后阶数下的相关性。
-
滞后阶数(Lag Order): 滞后阶数是指在计算自相关系数时,选择的时间间隔或滞后步数。它表示当前时刻的观测值与之前第几个时刻的观测值进行比较,用于计算相关系数。选择适当的滞后阶数对于准确描述时间序列的相关性非常重要。
ADF平稳性检验
ADF( Augmented Dickey-Fuller test )检验是一种常用的平稳性检验方法。它基于单位根检验,通过检验时间序列的单位根是否存在来判断序列的平稳性。如果拒绝原假设(存在单位根),则序列不平稳。
# 2. ADF检验
def adf_test(timeseries):
result = adfuller(timeseries)
print("ADF 统计量:", result[0])
print("p-value:", result[1])
if result[1] <= 0.05:
print("平稳")
else:
print("不平稳")
print("ADF检验结果:")
adf_test(ts_data)
Ljung-Box纯随机性检验
Ljung-Box检验的原假设是序列中的自相关系数为零,即序列呈现纯随机性。检验的统计量计算了一系列滞后阶数下的样本自相关系数的平方和,然后与自由度为滞后阶数的卡方分布进行比较。
# 3. 纯随机数检验
def pure_random_test(timeseries):
n = len(timeseries)
acf_vals = np.correlate(timeseries, timeseries, mode='full') / np.sum(timeseries**2)
lag = int(np.floor(n/2))
if acf_vals[lag] <= 0.05:
print("不纯随机")
else:
print("纯随机")
print("\n纯随机数检验:")
pure_random_test(ts_data)
指数平滑模型
分为季节性模型和非季节性模型
一元时间序列分析模型
时间序列的平稳性
若时间序列{xt}满足以下三个条件:
-
E ( x t ) = E ( x t − s ) = u ( 均值固定为常数 ) \large E(x_t)=E(x_{t-s})=u(均值固定为常数) E(xt)=E(xt−s)=u(均值固定为常数)
-
V a r ( x i ) = V a r ( x t − s ) = σ 2 ( 方差存在且为常数 ) \large Var(x_{i})=Var(x_{t-s})=\sigma^2(方差存在且为常数) Var(xi)=Var(xt−s)=σ2(方差存在且为常数)
-
C o v ( x t , x t − s ) = γ s ( 协方差只和间隔 s 有关,与 t 无关 ) \large Cov(x_t,x_{t-s})=\gamma_s(协方差只和间隔s有关,与t无关) Cov(xt,xt−s)=γs(协方差只和间隔s有关,与t无关)
则称{xt}为协方差平稳,又称弱平稳
如果对任意的t和h,多为随机变量 x t 1 , x t 2 , . . . , x t k \large x_{t_1},x_{t_2},...,x_{t_k} xt1,xt2,...,xtk和 x t 1 + h , x t 2 + h , . . . x t k + h \large x_{t_1+h},x_{t_2+h},...x_{t_k+h} xt1+h,xt2+h,...xtk+h的联合分布相同,则称{xt}为严格平稳
若时间序列{xt}满足以下三个条件:
-
E ( x t ) = E ( t t − s ) = 0 \large E(x_t)=E(t_{t-s})=0 E(xt)=E(tt−s)=0
-
V a r ( x i ) = V a r ( x t − s ) = σ 2 ( 方差存在且为常数 ) \large Var(x_{i})=Var(x_{t-s})=\sigma^2(方差存在且为常数) Var(xi)=Var(xt−s)=σ2(方差存在且为常数)
-
C o v ( x t , x t − s ) = γ s \large Cov(x_t,x_{t-s})=\gamma_s Cov(xt,xt−s)=γs
则称{xt}为白噪声序列
差分方程
将某个时间序列变量表示为该变量的滞后项、时间和其他变量的函数,这个函数方程被称为差分方程
-
y t = α 0 + α 1 y t − 1 + α 2 y t 2 + . . . + α p y y − p + ε t ( 自回归 A R ( p ) 模型 ) \large y_t=\alpha_0+\alpha_1y_{t-1}+\alpha_2y_{t_2}+...+\alpha_py_{y-p}+\varepsilon_t(自回归AR(p)模型) yt=α0+α1yt−1+α2yt2+...+αpyy−p+εt(自回归AR(p)模型)
-
y t = ε t + β 1 ε t − 1 + β 2 ε t − 2 + . . . + β q ε t − q ( 移动平均 M A ( q ) 模型 ) \large y_t=\varepsilon_t+\beta_1\varepsilon_{t-1}+\beta_2\varepsilon_{t-2}+...+\beta_q\varepsilon_{t-q}(移动平均MA(q)模型) yt=εt+β1εt−1+β2εt−2+...+βqεt−q(移动平均MA(q)模型)
-
y t = α 0 + ∑ i = 1 p = 1 α i y t − i + ε t + ∑ i = 1 q β i ε t − i ( 自回归移动平均 A R M A ( p , q ) 模型 ) \large y_t=\alpha_0+\sum_{i=1}^{p=1}\alpha_iy_{t-i}+\varepsilon_{t}+\sum_{i=1}^{q}\beta_i\varepsilon_{t-i} (自回归移动平均ARMA(p,q)模型) yt=α0+i=1∑p=1αiyt−i+εt+i=1∑qβiεt−i(自回归移动平均ARMA(p,q)模型)
-
y t = a + b y t − 1 + c z t + d z t − 1 + ε t \large y_t=a+by_{t-1}+cz_t+dz_{t-1}+\varepsilon_t yt=a+byt−1+czt+dzt−1+εt
-
y t = a + b y t − 1 + c t + ε t \large y_t=a+by_{t-1}+ct+\varepsilon_t yt=a+byt−1+ct+εt
ε t \varepsilon_t εt不做说明一般指白噪声序列
差分方程的特征方程
差分方程的齐次部分:只包含该变量自身和它的滞后项的式子
y
t
=
α
0
+
∑
i
=
1
p
α
i
y
t
−
i
+
ε
t
+
∑
i
=
1
q
β
i
ε
t
−
i
(
自回归移动平均
A
R
M
A
(
p
,
q
)
模型
)
\large y_t=\alpha_0+\sum_{i=1}^p\alpha_iy_{t-i}+\varepsilon_{t}+\sum_{i=1}^{q}\beta_i\varepsilon_{t-i} (自回归移动平均ARMA(p,q)模型)
yt=α0+i=1∑pαiyt−i+εt+i=1∑qβiεt−i(自回归移动平均ARMA(p,q)模型)
齐次部分:
y
t
=
∑
i
=
1
p
α
i
y
t
−
i
\large y_t=\sum_{i=1}^p\alpha_iy_{t-i}
yt=∑i=1pαiyt−i
将其次部分转化为特征方程:令yt=xt带入:
x
t
=
∑
i
=
1
p
α
i
x
t
−
i
⇒
x
p
=
α
1
x
p
−
1
+
α
2
x
p
−
2
+
.
.
+
α
p
\large x^t=\sum_{i=1}^p\alpha_ix^{t-i}\Rightarrow x^p=\alpha_1x^{p-1}+\alpha_2x^{p-2}+..+\alpha_p
xt=i=1∑pαixt−i⇒xp=α1xp−1+α2xp−2+..+αp
特征方程是一个p阶多项式,对应可求出p个解(虚根和实根都有),p个解的模长大小决定了ARMA(p,q)模型的{yi}是否平稳
滞后算子
用符号L表示滞后算子lag operator:
L
i
y
t
=
y
t
−
i
\large L^iy_t=y_{t-i}
Liyt=yt−i,满足如下性质:
(
1
)
L
C
=
C
(
C
为常数
)
(
2
)
(
L
i
+
L
j
)
y
t
=
y
t
−
i
+
y
t
−
j
(
3
)
L
i
L
j
y
t
=
y
t
−
i
−
j
(1)LC=C(C为常数)\\ (2)(L^i+L^j)y_t=y_{t-i}+y_{t-j}\\ (3)L^iL^jy_t=y_{t-i-j}
(1)LC=C(C为常数)(2)(Li+Lj)yt=yt−i+yt−j(3)LiLjyt=yt−i−j
-
一阶差分
Δ y t = y t − y t − 1 = ( 1 − L ) y t \large \Delta y_t=y_t-y_{t-1}=(1-L)y_t Δyt=yt−yt−1=(1−L)yt -
二阶差分
Δ 2 y t = Δ ( Δ y t ) = Δ ( y t − y t − 1 ) = ( 1 − L ) 2 y t \large\Delta^2y_t=\Delta(\Delta y_t)=\Delta(y_t-y_{t-1})=(1-L)^2y_t Δ2yt=Δ(Δyt)=Δ(yt−yt−1)=(1−L)2yt -
n阶差分
Δ n y t = ( 1 − L ) n y t \large \Delta^ny_t=(1-L)^ny_t Δnyt=(1−L)nyt
p,d,q
- 自回归阶数 §: 自回归阶数表示模型中考虑的历史观测值的数量,也就是当前值与过去多少个时刻的值有关。例如,ARIMA(p, d, q)中的p就表示自回归阶数。当p为0时,模型不考虑历史观测值,只根据当前时刻的值进行预测;当p为1时,模型根据当前时刻和前一个时刻的值进行预测,以此类推。
自回归模型基于时间序列的过去观测值,将其作为预测当前时刻值的依据。增加p的值可以考虑更多过去时刻的观测值,但过多的自回归项可能会导致模型过拟合。
- 差分阶数 (d): 差分阶数表示在进行建模之前,需要对时间序列进行多少次差分以使其平稳。平稳性是ARIMA模型的前提,它指的是序列的均值、方差和自协方差都不随时间变化。如果时间序列不是平稳的,就需要通过差分来减少或消除趋势、季节性等。
当d为0时,说明原始时间序列已经平稳,无需差分;当d为1时,说明进行一次一阶差分;以此类推。通常,需要适当的差分次数才能将时间序列转换为平稳序列,以便构建ARIMA模型。
- 滑动平均阶数 (q): 滑动平均阶数表示模型中考虑的过去滑动平均项的数量,也就是当前值与过去多少个时刻的平均值有关。例如,ARIMA(p, d, q)中的q就表示滑动平均阶数。和自回归类似,滑动平均模型也用于考虑过去观测值对当前时刻值的影响。
增加q的值可以考虑更多过去时刻的平均值,但同样要注意过多的滑动平均项可能导致模型复杂性过高。
AR§
AR§是p阶的自回归模型,将自己的1到p阶滞后项视为自变量进行回归,只适用于平稳的时间序列。
y
t
=
α
0
+
α
1
y
t
−
1
+
α
2
y
t
2
+
.
.
.
+
α
p
y
y
−
p
+
ε
t
\large y_t=\alpha_0+\alpha_1y_{t-1}+\alpha_2y_{t_2}+...+\alpha_py_{y-p}+\varepsilon_t
yt=α0+α1yt−1+α2yt2+...+αpyy−p+εt
自回归只能适用于预测与自身前期相关的经济现象,即受自身历史因素影响较大的经济现象,如矿的开采量,各种自然资源产量等。对于受社会因素影响较大的经济现象, 不宜采用自回归,而应使用可纳入其他变量的向量自回归模型(多元时间序列)。
将齐次部分转化为特征方程后,有:
x
p
=
α
1
x
p
−
1
+
α
2
x
p
−
2
+
.
.
.
+
α
p
\large x^p=\alpha_1x^{p-1}+\alpha_2x^{p-2}+...+\alpha_p
xp=α1xp−1+α2xp−2+...+αp
这样可以求出p个解:
- 若p个解模长均小于1,则{yt}序列平稳,AR§模型平稳
- 若p个解有k个解模长等于1,则{yt}为k阶单位根过程,可经过k阶差分变为平稳序列
- 若任意一个解模长大于1,则{yt}为爆炸过程
MA(q)
q阶移动平均过程
y
t
=
ε
t
+
β
1
ε
t
−
1
+
β
2
ε
t
−
2
+
.
.
.
+
β
q
ε
t
−
q
\large y_t=\varepsilon_t+\beta_1\varepsilon_{t-1}+\beta_2\varepsilon_{t-2}+...+\beta_q\varepsilon_{t-q}
yt=εt+β1εt−1+β2εt−2+...+βqεt−q
可知,移动平均模型具有可逆性:一阶MA可以转换为无穷阶自回归AR。
类似的,可逆性条件下,MA(q)可以转化为无穷阶自回归过程
平稳性:只要q是常数,MA(q)一定平稳。
ARMA(p,q)
自回归移动平均模型(Autoregressive Moving Average),将MA和AR结合起来
y
t
=
α
0
+
∑
i
=
1
p
=
1
α
i
y
t
−
i
+
ε
t
+
∑
i
=
1
q
β
i
ε
t
−
i
\large y_t=\alpha_0+\sum_{i=1}^{p=1}\alpha_iy_{t-i}+\varepsilon_{t}+\sum_{i=1}^{q}\beta_i\varepsilon_{t-i}
yt=α0+i=1∑p=1αiyt−i+εt+i=1∑qβiεt−i
平稳性只和自回归AR§部分有关
ACF自相关系数/PACF偏自相关系数
ACF表示时间序列与其自身在不同滞后阶数下的相关性。它衡量了时间序列当前时刻与之前各个时刻之间的相关程度。ACF函数的图像可以帮助我们观察时间序列的周期性和季节性,以及确定最合适的滞后阶数。
-
相关系数:X和Y是两个随机变量, ρ = c o v ( X , Y ) V a r ( x ) V a r ( y ) \large\rho=\frac{cov(X,Y)}{\sqrt{Var(x)}\sqrt{Var(y)}} ρ=Var(x)Var(y)cov(X,Y)表示X和Y的相关系数(线性关系)
-
自相关系数(autocorrelation):假定时间序列平稳,定义自相关系数为:
ρ s = c o v ( x , x t − s ) V a r ( x t ) V a r ( x t − s ) = γ s γ 0 , γ 0 = C o v ( x t , x t ) = V a r ( x t ) \large\rho_s=\frac{cov(x,x_{t-s})}{\sqrt{Var(x_t)}\sqrt{Var(x_{t-s})}}=\frac{\gamma_s}{\gamma_0},\gamma_0=Cov(x_t,x_t)=Var(x_t) ρs=Var(xt)Var(xt−s)cov(x,xt−s)=γ0γs,γ0=Cov(xt,xt)=Var(xt)
将 ρ \rho ρ视为s的函数,则称该函数为自相关函数,函数图像为自相关图
由于真实的自相关系数未知,只能依赖样本数据进行估计
对于白噪声序列:
ρ s = { 1 , s = 0 0 , s ≠ 0 \large\rho_s= \begin{cases} 1, s=0\\ 0, s\ne0 \end{cases} ρs=⎩ ⎨ ⎧1,s=00,s=0
PACF表示两个时刻之间的相关性,消除了中间时刻的影响,即在考虑其他滞后阶数的影响时,直接比较两个时刻之间的关系。PACF可以帮助我们确定时间序列的阶数(order),即在AR(自回归)或MA(移动平均)模型中需要考虑的滞后阶数。
model | ACF | PACF |
---|---|---|
AR§ | 逐渐衰减,即拖尾 | p阶后截尾 |
MA(q) | q阶后截尾 | 逐渐衰减,即拖尾 |
ARMA(p,q) | 逐渐衰减,即拖尾 | 逐渐衰减,即拖尾 |
# 4. 绘制ACF图和PACF图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
plot_acf(ts_data, ax=ax1)
plot_pacf(ts_data, ax=ax2)
plt.show()
AIC/BIC准则
Akaike information criterion:AIC=2(模型中参数的个数)-2ln(模型的极大似然函数值)
Bayesian information criterion:BIC=ln(T)(模型参数的个数)-2ln(模型的极大似然函数值)
- 样本个数:T
- 模型中参数的个数n:反应模型复杂程度
- 模型的极大似然函数值:反应模型对于数据解释程度
使AIC或BIC最小,BIC对模型复杂程度的惩罚系数更大,因此BIC往往比AIC选择的模型更简洁
模型检验
估计完成时间序列模型后,我们需要对残差进行白噪声检验,如果残差是白噪声,则说明我们选取的模型能完全识别出时间序列数据的规律,即模型可接受;如果残差不是白噪声,则说明还有部分信息没有被模型所识别,我们需要修正模型来识别这一部分的信息
Q检验:
H
0
:
ρ
1
=
ρ
2
=
.
.
.
=
ρ
s
=
0
,
H
1
:
ρ
i
(
i
=
1
,
2
,
.
.
.
s
)
至少有一个不为
0
若
H
0
成立,统计量
Q
=
T
(
T
+
2
)
∑
k
=
1
s
r
k
2
T
−
k
∼
χ
s
−
n
2
H_0:\rho_1=\rho_2=...=\rho_s=0,H_1:\rho_i(i=1,2,...s)至少有一个不为0\\ 若H_0成立,统计量Q=T(T+2)\sum_{k=1}^s\frac{r_k^2}{T-k}\sim\chi^2_{s-n}
H0:ρ1=ρ2=...=ρs=0,H1:ρi(i=1,2,...s)至少有一个不为0若H0成立,统计量Q=T(T+2)k=1∑sT−krk2∼χs−n2
p<0/05:拒绝原假设,模型需要修正
- T:样本个数
- n:未知参数个数,例如ARMA(p,q):n=p+q+1
- s一般取8,16,24等,SPSS里取的是18
ARIMA(p,d,q)
y t ′ t = α 0 + ∑ i = 1 p = 1 α i y t − i ′ + ε t + ∑ i = 1 q β i ε t − i 且 y t ′ = Δ d y t = ( 1 − L ) d y t \large y_t't=\alpha_0+\sum_{i=1}^{p=1}\alpha_iy'_{t-i}+\varepsilon_{t}+\sum_{i=1}^{q}\beta_i\varepsilon_{t-i}\\且y'_t=\Delta^dy_t=(1-L)^dy_t yt′t=α0+i=1∑p=1αiyt−i′+εt+i=1∑qβiεt−i且yt′=Δdyt=(1−L)dyt
SARIMA(p,d,q)
# 5. 构建ARIMA模型并进行预测
# p 是自回归阶数,d 是差分阶数(在这里设置为0,即不进行差分),q 是滑动平均阶数。
p, d, q = 1, 0, 1
model = ARIMA(ts_data, order=(p, d, q))
model_fit = model.fit()
# 预测的步数
forecast_steps = 5
forecast = model_fit.forecast(steps=forecast_steps)
# 6. 预测结果可视化
plt.figure()
plt.plot(time_steps, ts_data, label='原始数据')
plt.plot(np.arange(n_samples, n_samples+forecast_steps), forecast, label='预测数据', color='red')
plt.xlabel('时间')
plt.ylabel('值')
plt.legend()
plt.show()