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

Koorye

2803人浏览 · 2021-07-29 23:33:12

Koorye · 2021-07-29 23:33:12 发布

24. 时间序列模型

定义

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

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

长期趋势变动 $T_t$ ：朝一定方向的变化趋势；
季节变动 $S_t$ ；
循环变动 $C_t$ ：周期一年以上，由非季节因素引起的涨落起伏波形相似的波动；
不规则变动 $R_t$ 。

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

加法模型 $y_t=T_t+S_t+C_t+R_t$ ；
乘法模型 $yt=Tt⋅St⋅Ct⋅Rty_t=T_t\cdot S_t\cdot C_t\cdot R_t$ ；
混合模型。

其中 $y_t$ 是目标的观测记录， $E(Rt)=0,E(Rt2)=σ2E(R_t)=0,E(R_t^2)=\sigma^2$ 。

移动平均法

当时间序列由于受周期变动和不规则变动的影响，起伏较大，不易显示出发展趋势时，可用移动平均法消除这些因素的影响。

简单移动平均法

设观测序列为 $y1,…,yTy_1,\dots,y_T$ ，取移动平均项数 $N < T$ ，一次简单移动平均值计算公式为
$M_t^{(1)}=\frac1N(y_t+y_{t-1}+\dots+y_{t-N+1})$

$=M_{t-1}^{(1)}+\frac1N(y_t-t_{t-N}).$

令 $y^t+1=Mt(1)\hat{y}_{t+1}=M_t^{(1)}$ ，得到预测模型，并有误差
$S=∑t=N+1T(y^t−yt)2T−N. S=\sqrt{\frac{\sum_{t=N+1}^T(\hat{y}_t-y_t)^2}{T-N}}.$

加权移动平均法

设时间序列为 $y1,…,yt,…y_1,\dots,y_t,\dots$ ，加权移动平均公式为
$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,$
其中 $M_{tw}$ 为 $t$ 期加权移动平均数， $w_i$ 为 $y_{t-i+1}$ 的权数，得预测公式
$y^t+1=Mtw. \hat{y}_{t+1}=M_{tw}.$

趋势移动平均法

一次移动得平均数为
$M_t^{(1)}=\frac1N(y_t+y_{t-1}+\dots+y_{t-N+1}),$
在一次移动平均的基础上再进行一次移动平均，就是二次移动平均，计算公式为
$M_t^{(2)}=\frac1N(M_t^{(1)}+\dots+M_{t-N+1}^{(1)})$

$=M_{t-1}^{(2)}+\frac1N(M_t^{(1)}-M_{t-N}^{(1)}).$

设时间序列某时期开始具有直线趋势，则设
$y^t+T=at+btT,T=1,2,… \hat{y}_{t+T}=a_t+b_tT,T=1,2,\dots$
则有
$a_t=y_t,y_{t-1}=y_t-b_t,y_{t-2}=y_t-2b_t,$

$\dots$

$y_{t-N+1}=y_t-(N-1)b_t.$

故
$M_t^{(1)}=\frac{Ny_t-[1+\dots+(N-1)b_t]}{N}$

$=y_t-\frac{N-1}{2}b_t.$

即
$y_t-M_t^{(1)}=\frac{N-1}{2}b_t,$
同理有
$y_{t-1}-M_{t-1}^{(1)}=\frac{N-1}{2}b_{t-1},$
故
$y_t-y_{t-1}=M_t^{(1)}-M_{t-1}^{(1)}=b_t,$
同理有
$M_t^{(1)}-M_t^{(2)}=\frac{N-1}{2}b_t,$
于是有
$\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.$
一次移动平均实际上认为最近 N 期数据对未来影响相同，而 N 期以前的数据对未来没有影响。高阶移动平均则认为两端项权数小，中间项权数大，是对称的。

指数平滑法

一般说来，历史数据对未来的影响是递减的，指数平滑法可满足这一要求。

一次指数平滑法

设时间序列 $y1,y2,…,yty_1,y_2,\dots,y_t$ ， $α\alpha$ 为加权系数，一次指数平滑公式为
$S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}$

$=S_{t-1}^{(1)}+\alpha(y_t-S_{t-1}^{(1)}).$

将式展开，有
$S_t^{(1)}=\alpha y_t+(1-\alpha)[\alpha y_{t-1}+(1-\alpha)S_{t-2}^{(1)}]$

$=\dots=\alpha\sum_{j=0}^{\infty}(1-\alpha)^jy_{t-j}.$

表明 $S_t^{(1)}$ 是全部历史数据的加权平均，系数分别为 $α,α(1−α),α(1−α)2,…\alpha,\alpha(1-\alpha),\alpha(1-\alpha)^2,\dots$ ，符合指数规律。

有预测模型
$y^t+1=St(1)=αyt+(1−α)y^t. \hat{y}_{t+1}=S_t^{(1)}=\alpha y_t+(1-\alpha)\hat{y}_t.$

二次指数平滑法

一次指数平滑法存在明显的滞后偏差，可以再作二次指数平滑，有
$S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)},$

$S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)}.$

设时间序列从某时期开始具有直线趋势，可得
$y^t+T=at+btT,T=1,2,… \hat{y}_{t+T}=a_t+b_tT,T=1,2,\dots$

$\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.$

三次指数平滑法

时间序列变动表现为二次曲线趋势时，需要用三次指数平滑法，即
$S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)},$

$S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_{t-1}^{(2)}.$

$S_t^{(3)}=\alpha S_t^{(2)}+(1-\alpha)S_{t-1}^{(3)}.$

预测模型为
$y^t+T=at+btT+CtT2,T=1,2,… \hat{y}_{t+T}=a_t+b_tT+C_tT^2,T=1,2,\dots$

$\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.$

差分指数平滑法

当时间序列变动具有直线趋势时，用一次指数平滑法会出现滞后偏差。可以先对源数据进行差分，构成一个平稳的新序列 (消除直线趋势)，再使用指数平滑法预测，再与变量当前实际值迭加作为预测值。
$\nabla y_t=y_t-y_{t-1},$

$∇y^t+1=α∇yt+(1−α)∇y^t, \nabla\hat{y}_{t+1}=\alpha\nabla y_t+(1-\alpha)\nabla\hat{y}_t,$

$y^t+1=∇y^t+1+yt. \hat{y}_{t+1}=\nabla\hat{y}_{t+1}+y_t.$

同样的，当时间序列呈现二次曲线增长时，可用二阶差分指数平滑模型预测
$\nabla y_t=y_t-y_{t-1},$

$\nabla^2 y_t=\nabla y_t-\nabla y_{t-1}$

$∇2y^t+1=α∇2yt+(1−α)∇2y^t, \nabla^2\hat{y}_{t+1}=\alpha\nabla^2y_t+(1-\alpha)\nabla^2\hat{y}_t,$

$y^t+1=∇2y^t+1+∇yt+yt. \hat{y}_{t+1}=\nabla^2\hat{y}_{t+1}+\nabla y_{t}+y_t.$

自适应滤波法

基本预测公式为
$y^t+1=w1yt+w2yt−1+⋯+wNyt−N+1 \hat{y}_{t+1}=w_1y_t+w_2y_{t-1}+\dots+w_Ny_{t-N+1}$

$=\sum_{i=1}^Nw_iy_{t-i+1}.$

调整权数的公式为
$w_i=w_i+2k\cdot e_{i+1}y_{t-i+1}.$
其中 $k$ 为学习常数， $e_{i+1}$ 为第 $t + 1$ 期预测误差。

步骤：

设定初始权数、学习常数和开始时期 $t$ ；
根据预测公式计算 $t + 1$ 期的预测值 $y^t+1\hat{y}_{t+1}$ ；
计算预测误差 $et+1=yt+1−y^t+1e_{t+1}=y_{t+1}-\hat{y}_{t+1}$ ；
根据调整公式更新权数；
$t$ 进 1，回到步骤 2，直到 $t$ 超出序列的最后一个时期，迭代结果。

季节性时间序列预测

使用季节系数法，步骤如下：

计算每年所有季度的算数平均值
$\overline{a}=\frac1k\sum_{i=1}^m\sum_{j=1}^na_{ij},k=mn,$
其中 $i, j$ 分别表示年份、季度的序号；
计算同季度的算数平均值
$\overline{a}_{\cdot j}=\frac1m\sum_{i=1}^ma_{ij},j=1,2,\dots,n,$
计算季度系数 $bj=a‾⋅j/a‾b_j=\overline{a}_{\cdot j}/\overline{a}$ ；
预测下一年的年加权平均
$\overline{y}_{m+1}=\frac{\sum_{i=1}^mw_iy_i}{n\sum_{i=1}^nw_i},$
得到各季度的预测值
$y_{m+1,j}=b_j\overline{y}_{m+1}.$

ARMA 时间序列

设随机序列 $}\{X_t,t=0,\pm1,\pm2,\dots\}$ 满足

$E(Xt)=μ=常数E(X_t)=\mu=\text{常数}$ ，
$γt+k,t=γk\gamma_{t+k,t}=\gamma_k$ 与 $t$ 无关，

则称 $X_t$ 为平稳序列。

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

AR 序列

$X_t=\phi_1X_{t-1}+\phi_2X_{t-2}+\dots+\phi_pX_{t-p}+\varepsilon_t,$

其中 $εt\varepsilon_t$ 是均值为 0，方差为 $σε2\sigma_{\varepsilon}^2$ 的平稳白噪声，称 $X_t$ 是阶数为 $p$ 的自回归序列，记 $A R (p)$ 。

MA 序列

$X_t=\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\dots-\theta_q\varepsilon_{t-q},$

其中 $εt\varepsilon_t$ 是均值为 0，方差为 $σε2\sigma_{\varepsilon}^2$ 的平稳白噪声，称 $X_t$ 是阶数为 $q$ 的移动平均序列，记 $M A (q)$ 。

于是有
$X_t-\phi_1X_{t-1}+\phi_2X_{t-2}+\dots+\phi_pX_{t-p}$

$=\varepsilon_t-\theta_1\varepsilon_{t-1}-\theta_2\varepsilon_{t-2}-\dots-\theta_q\varepsilon_{t-q},$

其中 $εt\varepsilon_t$ 是均值为 0，方差为 $σε2\sigma_{\varepsilon}^2$ 的平稳白噪声，称 $X_t$ 是阶数为 $p, q$ 的自回归移动平均序列，记 $A R M A (p, q)$ 。

模型构建

AIC 准则：选 $p, q$ 使得
$min⁡AIC=nln⁡σ^ε2+2(p+q+1), \min AIC=n\ln\hat{\sigma}_{\varepsilon}^2+2(p+q+1),$
其中 $n$ 是样本容量， $σ^ε2\hat{\sigma}_{\varepsilon}^2$ 是 $σε2\sigma_{\varepsilon}^2$ 的估计值。

参数估计

可以使用矩估计、最小二乘估计、最大似然估计等方法。

模型预测

记 $X_{k+m}$ 的估计量为 $X^k(m)\hat{X}_k(m)$ ，对于 AR 序列
${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. \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.$
对于 MA 序列，有
${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. \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(q)=(X^k(1),…,X^k(q))T\hat{X}_k^{(q)}=(\hat{X}_k(1),\dots,\hat{X}_k(q))^T$ ，得
$X^k+1(q)=(θ110…0θ201…0⋮⋮⋱⋮⋮θq−100…1θq00…0)X^k(q)−(θ1θ2⋮θq)Xk+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^k0(q)=0\hat{X}_{k_0}^{(q)}=0$ ，当 $n$ 充分大后，初始误差的影响逐渐消失。

因此对于 ARMA 序列，令
$\phi_j^*=\left\{\begin{aligned} &\phi_j,j=1,2,\dots,p,\\ &0,j>p. \end{aligned}\right.$
有递推式
$X^k+1(q)=(−G110…0−G201…0⋮⋮⋱⋮⋮−Gq−100…1−Gq+ϕ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)}$

$+\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}.$

其中 $G_j$ 满足 $Xt=∑j=0∞Gjεt−jX_t=\sum_{j=0}^{\infty}G_j\varepsilon_{t-j}$ ，可令初值 $X^k0(q)=0\hat{X}_{k_0}^{(q)}=0$ 。

ARIMA 模型

对差分运算后得到的平稳序列用 ARMA 模型拟合，称 $A R I M A (p, d, q)$ 模型，其中 $p, d, q$ 分别表示 $A R$ 的项数、差分阶数、 $M A$ 的项数。

Python 代码

二次指数平滑法

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

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

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:	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

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

DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

2026深度解读：Work Agent长程任务的技术机制与落地边界

Work Agent的出现，把交互形态推向新的阶段：接收一个目标，自主拆解步骤、调用工具、校验中间结果，持续完成跨时段、跨工具的长链任务。长程任务执行，正是Work Agent与传统聊天机器人最核心的区分标尺，本文从技术原理、能力场景、管控机制、边界约束等维度，拆解这类智能体的真实落地水平。在行业实践案例中，一份行业分析报告的生成，就完整覆盖上述环节：智能体先检索行业公开数据，抓取竞品公开资料，整