【数学建模笔记 24】数学建模的时间序列模型

news/2024/12/22 19:55:37/

24. 时间序列模型

定义

时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序列的方法构成数据分析的一个重要领域,即时间序列分析。

一个时间序列往往是以下几类变化形式的叠加。

  • 长期趋势变动 T t T_t Tt:朝一定方向的变化趋势;
  • 季节变动 S t S_t St
  • 循环变动 C t C_t Ct:周期一年以上,由非季节因素引起的涨落起伏波形相似的波动;
  • 不规则变动 R t R_t Rt

常见的确定性时间序列模型有

  • 加法模型 y t = T t + S t + C t + R t y_t=T_t+S_t+C_t+R_t yt=Tt+St+Ct+Rt​;
  • 乘法模型 y t = T t ⋅ S t ⋅ C t ⋅ R t y_t=T_t\cdot S_t\cdot C_t\cdot R_t yt=TtStCtRt​;
  • 混合模型。

其中 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+yt1++ytN+1)

= M t − 1 ( 1 ) + 1 N ( y t − t t − N ) . =M_{t-1}^{(1)}+\frac1N(y_t-t_{t-N}). =Mt1(1)+N1(ytttN).

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=TNt=N+1T(y^tyt)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++wNytn+1,tN,
其中 M t w M_{tw} Mtw t t t 期加权移动平均数, w i w_i wi y t − i + 1 y_{t-i+1} yti+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+yt1++ytN+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)++MtN+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)}). =Mt1(2)+N1(Mt(1)MtN(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,yt1=ytbt,yt2=yt2bt,

… \dots

y t − N + 1 = y t − ( N − 1 ) b t . y_{t-N+1}=y_t-(N-1)b_t. ytN+1=yt(N1)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++(N1)bt]

= y t − N − 1 2 b t . =y_t-\frac{N-1}{2}b_t. =yt2N1bt.


y t − M t ( 1 ) = N − 1 2 b t , y_t-M_t^{(1)}=\frac{N-1}{2}b_t, ytMt(1)=2N1bt,
同理有
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}, yt1Mt1(1)=2N1bt1,

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, ytyt1=Mt(1)Mt1(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)=2N1bt,
于是有
{ 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=N12(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α)St1(1)

= S t − 1 ( 1 ) + α ( y t − S t − 1 ( 1 ) ) . =S_{t-1}^{(1)}+\alpha(y_t-S_{t-1}^{(1)}). =St1(1)+α(ytSt1(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α)[αyt1+(1α)St2(1)]

= ⋯ = α ∑ j = 0 ∞ ( 1 − α ) j y t − j . =\dots=\alpha\sum_{j=0}^{\infty}(1-\alpha)^jy_{t-j}. ==αj=0(1α)jytj.

表明 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α)St1(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α)St1(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α)St1(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α)St1(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α)St1(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α[(65α)St(1)   2(54α)St(2)+(43α)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=ytyt1,

∇ 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=ytyt1,

∇ 2 y t = ∇ y t − ∇ y t − 1 \nabla^2 y_t=\nabla y_t-\nabla y_{t-1} 2yt=ytyt1

∇ 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+w2yt1++wNytN+1

= ∑ i = 1 N w i y t − i + 1 . =\sum_{i=1}^Nw_iy_{t-i+1}. =i=1Nwiyti+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+2kei+1yti+1.
其中 k k k 为学习常数, e i + 1 e_{i+1} ei+1 为第 t + 1 t+1 t+1 期预测误差。

步骤:

  1. 设定初始权数、学习常数和开始时期 t t t
  2. 根据预测公式计算 t + 1 t+1 t+1 期的预测值 y ^ t + 1 \hat{y}_{t+1} y^t+1
  3. 计算预测误差 e t + 1 = y t + 1 − y ^ t + 1 e_{t+1}=y_{t+1}-\hat{y}_{t+1} et+1=yt+1y^t+1
  4. 根据调整公式更新权数;
  5. t t t 进 1,回到步骤 2,直到 t t 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=1mj=1naij,k=mn,
    其中 i , j i,j i,j 分别表示年份、季度的序号;

  2. 计算同季度的算数平均值
    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, aj=m1i=1maij,j=1,2,,n,

  3. 计算季度系数 b j = a ‾ ⋅ j / a ‾ b_j=\overline{a}_{\cdot j}/\overline{a} bj=aj/a​;

  4. 预测下一年的年加权平均
    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=ni=1nwii=1mwiyi,

  5. 得到各季度的预测值
    y m + 1 , j = b j y ‾ m + 1 . y_{m+1,j}=b_j\overline{y}_{m+1}. ym+1,j=bjym+1.

ARMA 时间序列

设随机序列 { X t , t = 0 , ± 1 , ± 2 , … } \{X_t,t=0,\pm1,\pm2,\dots\} {Xt,t=0,±1,±2,} 满足

  • E ( X t ) = μ = 常数 E(X_t)=\mu=\text{常数} E(Xt)=μ=常数​,
  • γ t + k , t = γ k \gamma_{t+k,t}=\gamma_k γt+k,t=γk t t t​ 无关,

则称 X t X_t Xt 为平稳序列。

ARMA 序列,即自回归移动平均序列,由两部分组成:

  • AR 序列

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=ϕ1Xt1+ϕ2Xt2++ϕpXtp+ε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)

  • MA 序列

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εt1θ2εt2θqεtq,

其中 ε 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ϕ1Xt1+ϕ2Xt2++ϕpXtp

= ε 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εt1θ2εt2θqεtq,

其中 ε 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+ϕ2Xk1++ϕpXkp+1,X^k(2)=ϕ1X^k(1)+ϕ2Xk++ϕpXkp+2,X^k(m)=ϕ1X^k(m1)+ϕ2X^k(m2)++ϕpX^(mp),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θq1θq100001000010X^k(q)θ1θ2θqXk+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)=G1G2Gq1Gq+ϕq100ϕq1010ϕq2001ϕ1X^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}. +G1G2GqXk+1+00j=q+1pϕjXk+q+1j.

其中 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=0Gjεtj,可令初值 X ^ k 0 ( q ) = 0 \hat{X}_{k_0}^{(q)}=0 X^k0(q)=0

ARIMA 模型

对差分运算后得到的平稳序列用 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 的项数。

Python 代码

二次指数平滑法

对某厂钢产量作预测,数据如下

t钢产量一次平滑二次平滑预测值
1203120312031
222342091.92049.272031
325662234.132104.7282152.8
428202409.8912196.2772418.99
530062588.7242314.0112715.054
630932740.0072441.812981.171
732772901.1052579.5983166.002
835143084.9732731.2113360.4
937703290.4812898.9923590.348
1041073535.4373089.9253849.752
114171.882
124362.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 pltplt.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()

输出如下:

tproductions1s2predict
012031.02031.0000002031.000000NaN
122234.02091.9000002049.2700002031.000000
232566.02234.1300002104.7280002152.800000
342820.02409.8910002196.2769002418.990000
453006.02588.7237002314.0109402715.054000
563093.02740.0065902441.8096352981.170500
673277.02901.1046132579.5981283166.002240
783514.03084.9732292731.2106593360.399591
893770.03290.4812602898.9918393590.348330
9104107.03535.4368823089.9253523849.751862
011NaNNaNNaN3980.948412
012NaNNaNNaN4171.881925
013NaNNaNNaN4362.815438
014NaNNaNNaN4553.748951
015NaNNaNNaN4744.682464

ARMA 模型

根据 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:numNo. Observations:20
Model:ARMA(2, 2)Log Likelihood-88.462
Method:css-mleS.D. of innovations17.694
Date:Thu, 29 Jul 2021AIC188.925
Time:22:06:42BIC194.899
Sample:0HQIC190.091
coefstd errzP>|z|[0.0250.975]
const75.79773.13024.2170.00069.66381.932
ar.L1.num1.53990.10614.5710.0001.3331.747
ar.L2.num-0.85670.104-8.2450.000-1.060-0.653
ma.L1.num-0.56120.339-1.6580.097-1.2250.102
ma.L2.num-0.43860.291-1.5090.131-1.0080.131
RealImaginaryModulusFrequency
AR.10.8987-0.5996j1.0804-0.0936
AR.20.8987+0.5996j1.08040.0936
MA.11.0001+0.0000j1.00010.0000
MA.2-2.2796+0.0000j2.27960.5000

ARIMA 模型

对 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


http://www.ppmy.cn/news/190349.html

相关文章

Linux后台开发工具箱-葵花宝典

目录 目录 1 1. 前言 5 2. 脚本类工具 5 2.1. 双引号和单引号 5 2.2. 环境变量和变量 5 2.3. sed命令-字符串文本操作 5 2.4. sed和awk使用外部变量 6 2.5. awk给外部变量赋值 6 2.6. 设置shell模式 7 2.7. 设置shell提示 7 2.8. 浮点计算 7 2.9. 变量指定默认值 8…

时间序列的ARIMA预测分析

测得某一地区一口井7年的地下水埋深数据&#xff0c;试预报第8年全年的地下水埋深。 1,平稳性检验 画原始序列图和自相关图 x[9.40,8.81,8.65,10.01,11.07,11.54,12.73,12.43,11.64,11.39,11.10,10.85,...10.71,10.24,8.48,9.88,10.31,10.53,9.55,6.51,7.75,7.80,5.96,5.21,.…

利用JVM在线调试工具排查线上问题(超实用)

点击上方“后端技术精选”&#xff0c;选择“置顶公众号” 技术文章第一时间送达&#xff01; 来源&#xff1a;牛哥的博客 cnblogs.com/nxlhero/p/11660854.html 在生产上我们经常会碰到一些不好排查的问题&#xff0c;例如线程安全问题&#xff0c;用最简单的threaddump或者h…

360环视辅助驾驶硬件系统方案---OV2715+DS90UB913+DS90UB914+FPGA

1、前置摄像头 高级驾驶员辅助系统中的摄像头系统可以分析视频内容&#xff0c;以便提供车道偏离警告&#xff08;LDW&#xff09;、自动车道保持辅助&#xff08;LKA&#xff09;、远光灯/近光灯控制和交通标志识别&#xff08;TSR&#xff09;。在前视黑白摄像头中&#xff…

jsp java json解析,jsp中获取json字符串,并解析

JqueryDemo1 function showData() {var str={ "name": "John" };//json标准格式 var obj = eval(( + str + )); alert( obj.name); var str2="{ name: John }"; var obj2 = eval(( + str2 + )); alert( obj2.name); var str3={"GetUserPost…

HIT 2715 Matrix3(最大费用最大流)

Matrix3 My Tags (Edit) Source : bin3 Time limit : 5 sec Memory limit : 64 M Submitted : 302, Accepted : 74 Zhouguyue is a "驴友" and nowadays he likes traveling on an N * N matrix with a non-negative number in each grid, and each grid has a…

Ubuntu 18.04 cuda 9.0 双1080TI 只显示一张

追加&#xff1a;【已解决&#xff0c;有一张显卡硬件不稳定】 参考我的最终记录&#xff1a; https://blog.csdn.net/u012911347/article/details/82854018 这又是一篇关于cuda和nvidia的博客&#xff0c;暂时解决了显卡就只显示一张和cuda无法使用的问题。 如果你想了解更…

dell服务器接2k显示器,4K、2K已成主流DELL高分辨率显示器推荐

【IT168 资讯】高分辨率显示器的推出,让人们的视觉体验提升到了新的层次。而高分辨率显示器如今也正在从专业市场定位逐渐走进日常用户家中。目前,市面上可见到的桌面显示器的分辨率虽然已经达到5K,不过,如此之高的分辨率要普及还需时日,所以有着高分辨率需求的用户不如更…