时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序列的方法构成数据分析的一个重要领域，即时间序列分析。
一个时间序列往往是以下几类变化形式的叠加。
常见的确定性时间序列模型有
其中 y t y_t yt 是目标的观测记录， E ( R t ) = 0 , E ( R t 2 ) = σ 2 E(R_t)=0,E(R_t^2)=\sigma^2 E(Rt)=0,E(Rt2)=σ2。
当时间序列由于受周期变动和不规则变动的影响，起伏较大，不易显示出发展趋势时，可用移动平均法消除这些因素的影响。
设观测序列为
y
1
,
…
,
y
T
y_1,\dots,y_T
y1,…,yT，取移动平均项数
N
<
T
N<T
N<T，一次简单移动平均值计算公式为
M
t
(
1
)
=
1
N
(
y
t
+
y
t
−
1
+
⋯
+
y
t
−
N
+
1
)
M_t^{(1)}=\frac1N(y_t+y_{t-1}+\dots+y_{t-N+1})
Mt(1)=N1(yt+yt−1+⋯+yt−N+1)
= M t − 1 ( 1 ) + 1 N ( y t − t t − N ) . =M_{t-1}^{(1)}+\frac1N(y_t-t_{t-N}). =Mt−1(1)+N1(yt−tt−N).
令
y
^
t
+
1
=
M
t
(
1
)
\hat{y}_{t+1}=M_t^{(1)}
y^t+1=Mt(1)，得到预测模型，并有误差
S
=
∑
t
=
N
+
1
T
(
y
^
t
−
y
t
)
2
T
−
N
.
S=\sqrt{\frac{\sum_{t=N+1}^T(\hat{y}_t-y_t)^2}{T-N}}.
S=T−N∑t=N+1T(y^t−yt)2
.
设时间序列为
y
1
,
…
,
y
t
,
…
y_1,\dots,y_t,\dots
y1,…,yt,…，加权移动平均公式为
M
t
w
=
w
1
y
t
+
w
2
y
2
+
⋯
+
w
N
y
t
−
n
+
1
w
1
+
w
2
+
⋯
+
w
N
,
t
≥
N
,
M_{tw}=\frac{w_1y_t+w_2y_2+\dots+w_Ny_{t-n+1}}{w_1+w_2+\dots+w_N},t\ge N,
Mtw=w1+w2+⋯+wNw1yt+w2y2+⋯+wNyt−n+1,t≥N,
其中
M
t
w
M_{tw}
Mtw 为
t
t
t 期加权移动平均数，
w
i
w_i
wi 为
y
t
−
i
+
1
y_{t-i+1}
yt−i+1 的权数，得预测公式
y
^
t
+
1
=
M
t
w
.
\hat{y}_{t+1}=M_{tw}.
y^t+1=Mtw.
一次移动得平均数为
M
t
(
1
)
=
1
N
(
y
t
+
y
t
−
1
+
⋯
+
y
t
−
N
+
1
)
,
M_t^{(1)}=\frac1N(y_t+y_{t-1}+\dots+y_{t-N+1}),
Mt(1)=N1(yt+yt−1+⋯+yt−N+1),
在一次移动平均的基础上再进行一次移动平均，就是二次移动平均，计算公式为
M
t
(
2
)
=
1
N
(
M
t
(
1
)
+
⋯
+
M
t
−
N
+
1
(
1
)
)
M_t^{(2)}=\frac1N(M_t^{(1)}+\dots+M_{t-N+1}^{(1)})
Mt(2)=N1(Mt(1)+⋯+Mt−N+1(1))
= M t − 1 ( 2 ) + 1 N ( M t ( 1 ) − M t − N ( 1 ) ) . =M_{t-1}^{(2)}+\frac1N(M_t^{(1)}-M_{t-N}^{(1)}). =Mt−1(2)+N1(Mt(1)−Mt−N(1)).
设时间序列某时期开始具有直线趋势，则设
y
^
t
+
T
=
a
t
+
b
t
T
,
T
=
1
,
2
,
…
\hat{y}_{t+T}=a_t+b_tT,T=1,2,\dots
y^t+T=at+btT,T=1,2,…
则有
a
t
=
y
t
,
y
t
−
1
=
y
t
−
b
t
,
y
t
−
2
=
y
t
−
2
b
t
,
a_t=y_t,y_{t-1}=y_t-b_t,y_{t-2}=y_t-2b_t,
at=yt,yt−1=yt−bt,yt−2=yt−2bt,
… \dots …
y t − N + 1 = y t − ( N − 1 ) b t . y_{t-N+1}=y_t-(N-1)b_t. yt−N+1=yt−(N−1)bt.
故
M
t
(
1
)
=
N
y
t
−
[
1
+
⋯
+
(
N
−
1
)
b
t
]
N
M_t^{(1)}=\frac{Ny_t-[1+\dots+(N-1)b_t]}{N}
Mt(1)=NNyt−[1+⋯+(N−1)bt]
= y t − N − 1 2 b t . =y_t-\frac{N-1}{2}b_t. =yt−2N−1bt.
即
y
t
−
M
t
(
1
)
=
N
−
1
2
b
t
,
y_t-M_t^{(1)}=\frac{N-1}{2}b_t,
yt−Mt(1)=2N−1bt,
同理有
y
t
−
1
−
M
t
−
1
(
1
)
=
N
−
1
2
b
t
−
1
,
y_{t-1}-M_{t-1}^{(1)}=\frac{N-1}{2}b_{t-1},
yt−1−Mt−1(1)=2N−1bt−1,
故
y
t
−
y
t
−
1
=
M
t
(
1
)
−
M
t
−
1
(
1
)
=
b
t
,
y_t-y_{t-1}=M_t^{(1)}-M_{t-1}^{(1)}=b_t,
yt−yt−1=Mt(1)−Mt−1(1)=bt,
同理有
M
t
(
1
)
−
M
t
(
2
)
=
N
−
1
2
b
t
,
M_t^{(1)}-M_t^{(2)}=\frac{N-1}{2}b_t,
Mt(1)−Mt(2)=2N−1bt,
于是有
{
a
t
=
2
M
t
(
1
)
−
M
t
(
2
)
,
b
t
=
2
N
−
1
(
M
t
(
1
)
−
M
t
(
2
)
)
.
\left\{\begin{aligned} &a_t=2M_t^{(1)}-M_t^{(2)},\\ &b_t=\frac{2}{N-1}(M_t^{(1)}-M_t^{(2)}). \end{aligned}\right.
⎩⎪⎨⎪⎧at=2Mt(1)−Mt(2),bt=N−12(Mt(1)−Mt(2)).
一次移动平均实际上认为最近 N 期数据对未来影响相同，而 N 期以前的数据对未来没有影响。高阶移动平均则认为两端项权数小，中间项权数大，是对称的。
一般说来，历史数据对未来的影响是递减的，指数平滑法可满足这一要求。
设时间序列
y
1
,
y
2
,
…
,
y
t
y_1,y_2,\dots,y_t
y1,y2,…,yt，
α
\alpha
α 为加权系数，一次指数平滑公式为
S
t
(
1
)
=
α
y
t
+
(
1
−
α
)
S
t
−
1
(
1
)
S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}
St(1)=αyt+(1−α)St−1(1)
= S t − 1 ( 1 ) + α ( y t − S t − 1 ( 1 ) ) . =S_{t-1}^{(1)}+\alpha(y_t-S_{t-1}^{(1)}). =St−1(1)+α(yt−St−1(1)).
将式展开，有
S
t
(
1
)
=
α
y
t
+
(
1
−
α
)
[
α
y
t
−
1
+
(
1
−
α
)
S
t
−
2
(
1
)
]
S_t^{(1)}=\alpha y_t+(1-\alpha)[\alpha y_{t-1}+(1-\alpha)S_{t-2}^{(1)}]
St(1)=αyt+(1−α)[αyt−1+(1−α)St−2(1)]
= ⋯ = α ∑ j = 0 ∞ ( 1 − α ) j y t − j . =\dots=\alpha\sum_{j=0}^{\infty}(1-\alpha)^jy_{t-j}. =⋯=αj=0∑∞(1−α)jyt−j.
表明 S t ( 1 ) S_t^{(1)} St(1) 是全部历史数据的加权平均，系数分别为 α , α ( 1 − α ) , α ( 1 − α ) 2 , … \alpha,\alpha(1-\alpha),\alpha(1-\alpha)^2,\dots α,α(1−α),α(1−α)2,…，符合指数规律。
有预测模型
y
^
t
+
1
=
S
t
(
1
)
=
α
y
t
+
(
1
−
α
)
y
^
t
.
\hat{y}_{t+1}=S_t^{(1)}=\alpha y_t+(1-\alpha)\hat{y}_t.
y^t+1=St(1)=αyt+(1−α)y^t.
一次指数平滑法存在明显的滞后偏差，可以再作二次指数平滑，有
S
t
(
1
)
=
α
y
t
+
(
1
−
α
)
S
t
−
1
(
1
)
,
S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)},
St(1)=αyt+(1−α)St−1(1),
S t ( 2 ) = α S t ( 1 ) + ( 1 − α ) S t − 1 ( 2 ) . S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)}. St(2)=αSt(1)+(1−α)St−1(2).
设时间序列从某时期开始具有直线趋势，可得
y
^
t
+
T
=
a
t
+
b
t
T
,
T
=
1
,
2
,
…
\hat{y}_{t+T}=a_t+b_tT,T=1,2,\dots
y^t+T=at+btT,T=1,2,…
{ a t = 2 S t ( 1 ) − S t ( 2 ) , b t = α 1 − α ( S t ( 1 ) − S t ( 2 ) ) . \left\{\begin{aligned} &a_t=2S_t^{(1)}-S_t^{(2)},\\ &b_t=\frac{\alpha}{1-\alpha}(S_t^{(1)}-S_t^{(2)}). \end{aligned}\right. ⎩⎨⎧at=2St(1)−St(2),bt=1−αα(St(1)−St(2)).
时间序列变动表现为二次曲线趋势时，需要用三次指数平滑法，即
S
t
(
1
)
=
α
y
t
+
(
1
−
α
)
S
t
−
1
(
1
)
,
S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)},
St(1)=αyt+(1−α)St−1(1),
S t ( 2 ) = α S t ( 1 ) + ( 1 − α ) S t − 1 ( 2 ) . S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)}. St(2)=αSt(1)+(1−α)St−1(2).
S t ( 3 ) = α S t ( 2 ) + ( 1 − α ) S t − 1 ( 3 ) . S_t^{(3)}=\alpha S_t^{(2)}+(1-\alpha)S_{t-1}^{(3)}. St(3)=αSt(2)+(1−α)St−1(3).
预测模型为
y
^
t
+
T
=
a
t
+
b
t
T
+
C
t
T
2
,
T
=
1
,
2
,
…
\hat{y}_{t+T}=a_t+b_tT+C_tT^2,T=1,2,\dots
y^t+T=at+btT+CtT2,T=1,2,…
{ a t = 3 S t ( 1 ) − 3 S t ( 2 ) + S t ( 3 ) , b t = α 2 ( 1 − α ) 2 [ ( 6 − 5 α ) S t ( 1 ) − 2 ( 5 − 4 α ) S t ( 2 ) + ( 4 − 3 α ) S t ( 3 ) ] , c t = α 2 2 ( 1 − α ) 2 [ S t ( 1 ) − 2 S t ( 2 ) + S t ( 3 ) ] . \left\{\begin{aligned} &a_t=3S_t^{(1)}-3S_t^{(2)}+S_t^{(3)},\\ &b_t=\frac{\alpha}{2(1-\alpha)^2}[(6-5\alpha)S_t^{(1)}\\ &\ \ \ -2(5-4\alpha)S_t^{(2)}+(4-3\alpha)S_t^{(3)}],\\ &c_t=\frac{\alpha^2}{2(1-\alpha)^2}[S_t^{(1)}-2S_t^{(2)}+S_t^{(3)}]. \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧at=3St(1)−3St(2)+St(3),bt=2(1−α)2α[(6−5α)St(1) −2(5−4α)St(2)+(4−3α)St(3)],ct=2(1−α)2α2[St(1)−2St(2)+St(3)].
当时间序列变动具有直线趋势时，用一次指数平滑法会出现滞后偏差。可以先对源数据进行差分，构成一个平稳的新序列 (消除直线趋势)，再使用指数平滑法预测，再与变量当前实际值迭加作为预测值。
∇
y
t
=
y
t
−
y
t
−
1
,
\nabla y_t=y_t-y_{t-1},
∇yt=yt−yt−1,
∇ y ^ t + 1 = α ∇ y t + ( 1 − α ) ∇ y ^ t , \nabla\hat{y}_{t+1}=\alpha\nabla y_t+(1-\alpha)\nabla\hat{y}_t, ∇y^t+1=α∇yt+(1−α)∇y^t,
y ^ t + 1 = ∇ y ^ t + 1 + y t . \hat{y}_{t+1}=\nabla\hat{y}_{t+1}+y_t. y^t+1=∇y^t+1+yt.
同样的，当时间序列呈现二次曲线增长时，可用二阶差分指数平滑模型预测
∇
y
t
=
y
t
−
y
t
−
1
,
\nabla y_t=y_t-y_{t-1},
∇yt=yt−yt−1,
∇ 2 y t = ∇ y t − ∇ y t − 1 \nabla^2 y_t=\nabla y_t-\nabla y_{t-1} ∇2yt=∇yt−∇yt−1
∇ 2 y ^ t + 1 = α ∇ 2 y t + ( 1 − α ) ∇ 2 y ^ t , \nabla^2\hat{y}_{t+1}=\alpha\nabla^2y_t+(1-\alpha)\nabla^2\hat{y}_t, ∇2y^t+1=α∇2yt+(1−α)∇2y^t,
y ^ t + 1 = ∇ 2 y ^ t + 1 + ∇ y t + y t . \hat{y}_{t+1}=\nabla^2\hat{y}_{t+1}+\nabla y_{t}+y_t. y^t+1=∇2y^t+1+∇yt+yt.
基本预测公式为
y
^
t
+
1
=
w
1
y
t
+
w
2
y
t
−
1
+
⋯
+
w
N
y
t
−
N
+
1
\hat{y}_{t+1}=w_1y_t+w_2y_{t-1}+\dots+w_Ny_{t-N+1}
y^t+1=w1yt+w2yt−1+⋯+wNyt−N+1
= ∑ i = 1 N w i y t − i + 1 . =\sum_{i=1}^Nw_iy_{t-i+1}. =i=1∑Nwiyt−i+1.
调整权数的公式为
w
i
=
w
i
+
2
k
⋅
e
i
+
1
y
t
−
i
+
1
.
w_i=w_i+2k\cdot e_{i+1}y_{t-i+1}.
wi=wi+2k⋅ei+1yt−i+1.
其中
k
k
k 为学习常数，
e
i
+
1
e_{i+1}
ei+1 为第
t
+
1
t+1
t+1 期预测误差。
步骤：
使用季节系数法，步骤如下：
计算每年所有季度的算数平均值
a
‾
=
1
k
∑
i
=
1
m
∑
j
=
1
n
a
i
j
,
k
=
m
n
,
\overline{a}=\frac1k\sum_{i=1}^m\sum_{j=1}^na_{ij},k=mn,
a=k1i=1∑mj=1∑naij,k=mn,
其中
i
,
j
i,j
i,j 分别表示年份、季度的序号；
计算同季度的算数平均值
a
‾
⋅
j
=
1
m
∑
i
=
1
m
a
i
j
,
j
=
1
,
2
,
…
,
n
,
\overline{a}_{\cdot j}=\frac1m\sum_{i=1}^ma_{ij},j=1,2,\dots,n,
a⋅j=m1i=1∑maij,j=1,2,…,n,
计算季度系数 b j = a ‾ ⋅ j / a ‾ b_j=\overline{a}_{\cdot j}/\overline{a} bj=a⋅j/a；
预测下一年的年加权平均
y
‾
m
+
1
=
∑
i
=
1
m
w
i
y
i
n
∑
i
=
1
n
w
i
,
\overline{y}_{m+1}=\frac{\sum_{i=1}^mw_iy_i}{n\sum_{i=1}^nw_i},
ym+1=n∑i=1nwi∑i=1mwiyi,
得到各季度的预测值
y
m
+
1
,
j
=
b
j
y
‾
m
+
1
.
y_{m+1,j}=b_j\overline{y}_{m+1}.
ym+1,j=bjym+1.
设随机序列 { X t , t = 0 , ± 1 , ± 2 , … } \{X_t,t=0,\pm1,\pm2,\dots\} {Xt,t=0,±1,±2,…} 满足
则称 X t X_t Xt 为平稳序列。
ARMA 序列，即自回归移动平均序列，由两部分组成：
X t = ϕ 1 X t − 1 + ϕ 2 X t − 2 + ⋯ + ϕ p X t − p + ε t , X_t=\phi_1X_{t-1}+\phi_2X_{t-2}+\dots+\phi_pX_{t-p}+\varepsilon_t, Xt=ϕ1Xt−1+ϕ2Xt−2+⋯+ϕpXt−p+εt,
其中 ε t \varepsilon_t εt 是均值为 0，方差为 σ ε 2 \sigma_{\varepsilon}^2 σε2 的平稳白噪声，称 X t X_t Xt 是阶数为 p p p 的自回归序列，记 A R ( p ) AR(p) AR(p)。
X t = ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q , X_t=\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\dots-\theta_q\varepsilon_{t-q}, Xt=εt−θ1εt−1−θ2εt−2−⋯−θqεt−q,
其中 ε t \varepsilon_t εt 是均值为 0，方差为 σ ε 2 \sigma_{\varepsilon}^2 σε2 的平稳白噪声，称 X t X_t Xt 是阶数为 q q q 的移动平均序列，记 M A ( q ) MA(q) MA(q)。
于是有
X
t
−
ϕ
1
X
t
−
1
+
ϕ
2
X
t
−
2
+
⋯
+
ϕ
p
X
t
−
p
X_t-\phi_1X_{t-1}+\phi_2X_{t-2}+\dots+\phi_pX_{t-p}
Xt−ϕ1Xt−1+ϕ2Xt−2+⋯+ϕpXt−p
= ε t − θ 1 ε t − 1 − θ 2 ε t − 2 − ⋯ − θ q ε t − q , =\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\dots-\theta_q\varepsilon_{t-q}, =εt−θ1εt−1−θ2εt−2−⋯−θqεt−q,
其中 ε t \varepsilon_t εt 是均值为 0，方差为 σ ε 2 \sigma_{\varepsilon}^2 σε2 的平稳白噪声，称 X t X_t Xt 是阶数为 p , q p,q p,q 的自回归移动平均序列，记 A R M A ( p , q ) ARMA(p,q) ARMA(p,q)。
AIC 准则：选
p
,
q
p,q
p,q 使得
min
A
I
C
=
n
ln
σ
^
ε
2
+
2
(
p
+
q
+
1
)
,
\min AIC=n\ln\hat{\sigma}_{\varepsilon}^2+2(p+q+1),
minAIC=nlnσ^ε2+2(p+q+1),
其中
n
n
n 是样本容量，
σ
^
ε
2
\hat{\sigma}_{\varepsilon}^2
σ^ε2 是
σ
ε
2
\sigma_{\varepsilon}^2
σε2 的估计值。
可以使用矩估计、最小二乘估计、最大似然估计等方法。
记
X
k
+
m
X_{k+m}
Xk+m 的估计量为
X
^
k
(
m
)
\hat{X}_k(m)
X^k(m)，对于 AR 序列
{
X
^
k
(
1
)
=
ϕ
1
X
k
+
ϕ
2
X
k
−
1
+
⋯
+
ϕ
p
X
k
−
p
+
1
,
X
^
k
(
2
)
=
ϕ
1
X
^
k
(
1
)
+
ϕ
2
X
k
+
⋯
+
ϕ
p
X
k
−
p
+
2
,
…
X
^
k
(
m
)
=
ϕ
1
X
^
k
(
m
−
1
)
+
ϕ
2
X
^
k
(
m
−
2
)
+
⋯
+
ϕ
p
X
^
(
m
−
p
)
,
m
>
p
.
\left\{\begin{aligned} &\hat{X}_k(1)=\phi_1X_k+\phi_2X_{k-1}+\dots+\phi_pX_{k-p+1},\\ &\hat{X}_k(2)=\phi_1\hat{X}_k(1)+\phi_2X_{k}+\dots+\phi_pX_{k-p+2},\\ &\dots\\ &\hat{X}_k(m)=\phi_1\hat{X}_k(m-1)+\phi_2\hat{X}_k(m-2)+\dots+\phi_p\hat{X}(m-p),m>p.\\ \end{aligned}\right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧X^k(1)=ϕ1Xk+ϕ2Xk−1+⋯+ϕpXk−p+1,X^k(2)=ϕ1X^k(1)+ϕ2Xk+⋯+ϕpXk−p+2,…X^k(m)=ϕ1X^k(m−1)+ϕ2X^k(m−2)+⋯+ϕpX^(m−p),m>p.
对于 MA 序列，有
{
X
^
k
+
1
(
1
)
=
θ
1
X
^
k
(
1
)
+
X
^
k
(
2
)
−
θ
1
X
k
+
1
,
X
^
k
+
1
(
2
)
=
θ
2
X
^
k
(
1
)
+
X
^
k
(
3
)
−
θ
2
X
k
+
1
,
…
X
^
k
+
1
(
q
)
=
θ
q
X
^
k
(
1
)
+
X
^
k
(
q
)
−
θ
q
X
k
+
1
.
\left\{\begin{aligned} &\hat{X}_{k+1}(1)=\theta_1\hat{X}_k(1)+\hat{X}_k(2)-\theta_1X_{k+1},\\ &\hat{X}_{k+1}(2)=\theta_2\hat{X}_k(1)+\hat{X}_k(3)-\theta_2X_{k+1},\\ &\dots\\ &\hat{X}_{k+1}(q)=\theta_{q}\hat{X}_k(1)+\hat{X}_k(q)-\theta_{q}X_{k+1}.\\ \end{aligned}\right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧X^k+1(1)=θ1X^k(1)+X^k(2)−θ1Xk+1,X^k+1(2)=θ2X^k(1)+X^k(3)−θ2Xk+1,…X^k+1(q)=θqX^k(1)+X^k(q)−θqXk+1.
令
X
^
k
(
q
)
=
(
X
^
k
(
1
)
,
…
,
X
^
k
(
q
)
)
T
\hat{X}_k^{(q)}=(\hat{X}_k(1),\dots,\hat{X}_k(q))^T
X^k(q)=(X^k(1),…,X^k(q))T，得
X
^
k
+
1
(
q
)
=
(
θ
1
1
0
…
0
θ
2
0
1
…
0
⋮
⋮
⋱
⋮
⋮
θ
q
−
1
0
0
…
1
θ
q
0
0
…
0
)
X
^
k
(
q
)
−
(
θ
1
θ
2
⋮
θ
q
)
X
k
+
1.
\hat{X}_{k+1}^{(q)}=\begin{pmatrix} \theta_1&1&0&\dots&0\\ \theta_2&0&1&\dots&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\\ \theta_{q-1}&0&0&\dots&1\\ \theta_q&0&0&\dots&0\\ \end{pmatrix}\hat{X}_k^{(q)}-\begin{pmatrix} \theta_1\\ \theta_2\\ \vdots\\ \theta_q \end{pmatrix}X_{k+1.}
X^k+1(q)=⎝⎜⎜⎜⎜⎜⎛θ1θ2⋮θq−1θq10⋮0001⋱00……⋮……00⋮10⎠⎟⎟⎟⎟⎟⎞X^k(q)−⎝⎜⎜⎜⎛θ1θ2⋮θq⎠⎟⎟⎟⎞Xk+1.
可取初值
X
^
k
0
(
q
)
=
0
\hat{X}_{k_0}^{(q)}=0
X^k0(q)=0，当
n
n
n 充分大后，初始误差的影响逐渐消失。
因此对于 ARMA 序列，令
ϕ
j
∗
=
{
ϕ
j
,
j
=
1
,
2
,
…
,
p
,
0
,
j
>
p
.
\phi_j^*=\left\{\begin{aligned} &\phi_j,j=1,2,\dots,p,\\ &0,j>p. \end{aligned}\right.
ϕj∗={ϕj,j=1,2,…,p,0,j>p.
有递推式
X
^
k
+
1
(
q
)
=
(
−
G
1
1
0
…
0
−
G
2
0
1
…
0
⋮
⋮
⋱
⋮
⋮
−
G
q
−
1
0
0
…
1
−
G
q
+
ϕ
q
∗
ϕ
q
−
1
∗
ϕ
q
−
2
∗
…
ϕ
1
∗
)
X
^
k
(
q
)
\hat{X}_{k+1}^{(q)}=\begin{pmatrix} -G_1&1&0&\dots&0\\ -G_2&0&1&\dots&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\\ -G_{q-1}&0&0&\dots&1\\ -G_q+\phi_q^*&\phi_{q-1}^*&\phi_{q-2}^*&\dots&\phi_1^*\\ \end{pmatrix}\hat{X}_k^{(q)}
X^k+1(q)=⎝⎜⎜⎜⎜⎜⎛−G1−G2⋮−Gq−1−Gq+ϕq∗10⋮0ϕq−1∗01⋱0ϕq−2∗……⋮……00⋮1ϕ1∗⎠⎟⎟⎟⎟⎟⎞X^k(q)
+ ( G 1 G 2 ⋮ G q ) X k + 1 + ( 0 0 ⋮ ∑ j = q + 1 p ϕ j ∗ X k + q + 1 − j ) . +\begin{pmatrix} G_1\\ G_2\\ \vdots\\ G_q \end{pmatrix}X_{k+1}+\begin{pmatrix} 0\\ 0\\ \vdots\\ \sum_{j=q+1}^p\phi_j^*X_{k+q+1-j} \end{pmatrix}. +⎝⎜⎜⎜⎛G1G2⋮Gq⎠⎟⎟⎟⎞Xk+1+⎝⎜⎜⎜⎛00⋮∑j=q+1pϕj∗Xk+q+1−j⎠⎟⎟⎟⎞.
其中 G j G_j Gj 满足 X t = ∑ j = 0 ∞ G j ε t − j X_t=\sum_{j=0}^{\infty}G_j\varepsilon_{t-j} Xt=∑j=0∞Gjεt−j，可令初值 X ^ k 0 ( q ) = 0 \hat{X}_{k_0}^{(q)}=0 X^k0(q)=0。
对差分运算后得到的平稳序列用 ARMA 模型拟合，称 A R I M A ( p , d , q ) ARIMA(p,d,q) ARIMA(p,d,q) 模型，其中 p , d , q p,d,q p,d,q 分别表示 A R AR AR 的项数、差分阶数、 M A MA MA 的项数。
对某厂钢产量作预测，数据如下
| t | 钢产量 | 一次平滑 | 二次平滑 | 预测值 |
|---|---|---|---|---|
| 1 | 2031 | 2031 | 2031 | |
| 2 | 2234 | 2091.9 | 2049.27 | 2031 |
| 3 | 2566 | 2234.13 | 2104.728 | 2152.8 |
| 4 | 2820 | 2409.891 | 2196.277 | 2418.99 |
| 5 | 3006 | 2588.724 | 2314.011 | 2715.054 |
| 6 | 3093 | 2740.007 | 2441.81 | 2981.171 |
| 7 | 3277 | 2901.105 | 2579.598 | 3166.002 |
| 8 | 3514 | 3084.973 | 2731.211 | 3360.4 |
| 9 | 3770 | 3290.481 | 2898.992 | 3590.348 |
| 10 | 4107 | 3535.437 | 3089.925 | 3849.752 |
| 11 | 4171.882 | |||
| 12 | 4362.815 |
输出如下：
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date:
# @ function: ARIMA 模型
# %%
import numpy as np
import pandas as pd
# %%
# 源数据
df = pd.DataFrame({
't': [i for i in range(1, 11)],
'production': [2031, 2234, 2566, 2820, 3006,
3093, 3277, 3514, 3770, 4107],
})
# 设 alpha=.3，计算一次、二次指数平滑
alpha = .3
s1, s2 = [int(df['production'][0])],\
[int(df['production'][0])]
for i in range(1, len(df['t'])):
s1.append(alpha*df['production'][i] + (1-alpha)*s1[i-1])
s2.append(alpha*s1[i] + (1-alpha)*s2[i-1])
df['s1'] = s1
df['s2'] = s2
# %%
# 计算过去年的预测值，以及未来年的线性表达式
predict_list = [None]
for i in range(len(df['t'])-1):
a = 2*df['s1'][i] - df['s2'][i]
b = (alpha / (1-alpha)) * (df['s1'][i] - df['s2'][i])
predict_list.append(a + b)
t = 10
a = 2*df['s1'][t-1] - df['s2'][t-1]
b = (alpha / (1-alpha)) * (df['s1'][t-1] - df['s2'][t-1])
df['predict'] = predict_list
print('at =', a)
print('bt =', b)
# %%
# 计算未来年的预测值
pred_df = df.copy()
for i in range(5):
pred_df = pd.concat([pred_df, pd.DataFrame({
't': [10+i+1],
'predict': a + b*i,
})])
pred_df
# %%
# 画图
import matplotlib.pyplot as plt
plt.plot(pred_df['t'], pred_df['production'],label='production')
plt.scatter(pred_df['t'], pred_df['production'])
plt.plot(pred_df['t'], pred_df['s1'],label='s1')
plt.scatter(pred_df['t'], pred_df['s1'])
plt.plot(pred_df['t'], pred_df['s2'],label='s2')
plt.scatter(pred_df['t'], pred_df['s2'])
plt.plot(pred_df['t'], pred_df['predict'],label='predict')
plt.scatter(pred_df['t'], pred_df['predict'])
plt.legend()
输出如下：
| t | production | s1 | s2 | predict | |
|---|---|---|---|---|---|
| 0 | 1 | 2031.0 | 2031.000000 | 2031.000000 | NaN |
| 1 | 2 | 2234.0 | 2091.900000 | 2049.270000 | 2031.000000 |
| 2 | 3 | 2566.0 | 2234.130000 | 2104.728000 | 2152.800000 |
| 3 | 4 | 2820.0 | 2409.891000 | 2196.276900 | 2418.990000 |
| 4 | 5 | 3006.0 | 2588.723700 | 2314.010940 | 2715.054000 |
| 5 | 6 | 3093.0 | 2740.006590 | 2441.809635 | 2981.170500 |
| 6 | 7 | 3277.0 | 2901.104613 | 2579.598128 | 3166.002240 |
| 7 | 8 | 3514.0 | 3084.973229 | 2731.210659 | 3360.399591 |
| 8 | 9 | 3770.0 | 3290.481260 | 2898.991839 | 3590.348330 |
| 9 | 10 | 4107.0 | 3535.436882 | 3089.925352 | 3849.751862 |
| 0 | 11 | NaN | NaN | NaN | 3980.948412 |
| 0 | 12 | NaN | NaN | NaN | 4171.881925 |
| 0 | 13 | NaN | NaN | NaN | 4362.815438 |
| 0 | 14 | NaN | NaN | NaN | 4553.748951 |
| 0 | 15 | NaN | NaN | NaN | 4744.682464 |
根据 1971~1990 年的太阳黑子数量预测 1991~2000 年的太阳黑子数量，代码如下：
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-29
# @ function: ARMA 模型
# %%
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
# %%
# 源数据
df = pd.DataFrame({
'year': [i for i in range(1971, 1991)],
'num': [66.6, 68.9, 38, 34.5, 15.5,
12.6, 27.5, 92.5, 155.4, 154.6,
140.4, 115.9, 66.6, 45.9, 17.9,
3.4, 29.4, 100.2, 157.6, 142.6],
})
# %%
# 画 acf 图
plot_acf(df['num'])
# %%
# 画 pacf 图
plot_pacf(df['num'], lags=9)
# %%
# 建立模型，参考 acf、pacf 代入 p、q，观察 aic
str_list = []
for p in range(1, 6):
for q in range(1, 3):
model = sm.tsa.ARMA(df['num'], (p, q)).fit()
str_list.append('p = {}, q = {}, aic = {}'.format(p, q, model.aic))
for each in str_list:
print(each)
# %%
# 发现 p=2,q=2 时 aic 最小，取 p=2,q=2
model = sm.tsa.ARMA(df['num'], (2, 2)).fit()
model.summary()
# %%
# 预测和画图
plt.plot(df['year'], df['num'])
plt.scatter(df['year'], df['num'], label='actual')
year_list = [i for i in range(1971, 2001)]
plt.plot(year_list, model.predict(0, len(year_list)-1))
plt.scatter(year_list, model.predict(0, len(year_list)-1), label='predict')
plt.legend()
输出如下：
p = 1, q = 1, aic = 194.84701061998135
p = 1, q = 2, aic = 194.31427174057512
p = 2, q = 1, aic = 188.97865479513626
p = 2, q = 2, aic = 188.92457119015984
p = 3, q = 1, aic = 191.31192040477467
p = 3, q = 2, aic = 192.40349950348713
p = 4, q = 1, aic = 193.19411553305596
p = 4, q = 2, aic = 193.1505615705116
p = 5, q = 1, aic = 194.3858472180579
p = 5, q = 2, aic = 194.78028983920348
| Dep. Variable: | num | No. Observations: | 20 |
|---|---|---|---|
| Model: | ARMA(2, 2) | Log Likelihood | -88.462 |
| Method: | css-mle | S.D. of innovations | 17.694 |
| Date: | Thu, 29 Jul 2021 | AIC | 188.925 |
| Time: | 22:06:42 | BIC | 194.899 |
| Sample: | 0 | HQIC | 190.091 |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| const | 75.7977 | 3.130 | 24.217 | 0.000 | 69.663 | 81.932 |
| ar.L1.num | 1.5399 | 0.106 | 14.571 | 0.000 | 1.333 | 1.747 |
| ar.L2.num | -0.8567 | 0.104 | -8.245 | 0.000 | -1.060 | -0.653 |
| ma.L1.num | -0.5612 | 0.339 | -1.658 | 0.097 | -1.225 | 0.102 |
| ma.L2.num | -0.4386 | 0.291 | -1.509 | 0.131 | -1.008 | 0.131 |
| Real | Imaginary | Modulus | Frequency | |
|---|---|---|---|---|
| AR.1 | 0.8987 | -0.5996j | 1.0804 | -0.0936 |
| AR.2 | 0.8987 | +0.5996j | 1.0804 | 0.0936 |
| MA.1 | 1.0001 | +0.0000j | 1.0001 | 0.0000 |
| MA.2 | -2.2796 | +0.0000j | 2.2796 | 0.5000 |
对 austa 数据集使用 ARIMA 模型作预测，代码如下：
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date:
# @ function: ARIMA 模型
# %%
import numpy as np
import pandas as pd
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
# %%
# 源数据
df = pd.DataFrame({
'year': [i for i in range(1980, 2011)],
'val': [0.82989428, 0.85951092, 0.87668916, 0.86670716, 0.932052,
1.04826364, 1.3111932, 1.63756228, 2.0641074, 1.91268276,
2.03544572, 2.17721128, 2.38968344, 2.75059208, 3.0906664,
3.42664028, 3.83064908, 3.97190864, 3.83160036, 4.143101,
4.566551, 4.47541, 4.462796, 4.384829, 4.796861,
5.046211, 5.098759, 5.196519, 5.166843, 5.174744,
5.440894],
})
# %%
# 作一阶差分
df['val_diff1'] = df['val'].diff()
plt.plot(df['year'], df['val'])
plt.plot(df['year'], df['val_diff1'])
# %%
# ACF 图
plot_acf(df['val_diff1'][1:])
# %%
# PACF 图
plot_pacf(df['val_diff1'][1:], lags=14)
# %%
# 参考 ACF 和 PCAF 代入 p,q
str_list = []
for p in range(1, 4):
for q in range(0, 4):
model = ARIMA(df['val'], order=(p, 1, q))
res = model.fit(disp=0)
str_list.append('p = {}, q = {}, aic = {}'.format(p, q, res.aic))
for each in str_list:
print(each)
# %%
# p=2,q=0 使 aic 最小
model = ARIMA(df['val'], order=(2, 1, 0))
res = model.fit(disp=0)
res.summary()
# %%
# 作预测图
res.plot_predict(end=40)
输出如下：
p = 1, q = 0, aic = -13.718554829816128
p = 1, q = 1, aic = -12.644935232547795
p = 1, q = 2, aic = -13.023448862467816
p = 1, q = 3, aic = -11.034475643866699
p = 2, q = 0, aic = -13.534487801792835
p = 2, q = 1, aic = -11.848263794395535
p = 2, q = 2, aic = -11.032058544143538
p = 2, q = 3, aic = -9.036565460115241
p = 3, q = 0, aic = -11.784035132385696
p = 3, q = 1, aic = -9.84986571141387
p = 3, q = 2, aic = -8.617052624191913
p = 3, q = 3, aic = -5.945170779675806
