吴世涵的博客
0%

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

发表于 分类于
本文字数： 14k 阅读时长 ≈ 13 分钟
本篇文章主要介绫了数学建模中的时间序列模型，包括移动平均法、指数平滑法、趋势移动平均法、差分指数平滑法、自适应滤波法、季节性时间序列预测、ARMA 时间序列、ARIMA 模型等内容。

24. 时间序列模型

定义

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

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

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

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

  • 加法模型 ​；
  • 乘法模型 ​；
  • 混合模型。

其中 是目标的观测记录，。

移动平均法

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

简单移动平均法

设观测序列为 ，取移动平均项数 ，一次简单移动平均值计算公式为

令 ，得到预测模型，并有误差

加权移动平均法

设时间序列为 ，加权移动平均公式为 其中 为 期加权移动平均数， 为 的权数，得预测公式

趋势移动平均法

一次移动得平均数为 在一次移动平均的基础上再进行一次移动平均，就是二次移动平均，计算公式为

设时间序列某时期开始具有直线趋势，则设 则有

故

即 同理有 故 同理有 于是有 一次移动平均实际上认为最近 N 期数据对未来影响相同，而 N 期以前的数据对未来没有影响。高阶移动平均则认为两端项权数小，中间项权数大，是对称的。

指数平滑法

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

一次指数平滑法

设时间序列 ， 为加权系数，一次指数平滑公式为

将式展开，有

表明 是全部历史数据的加权平均，系数分别为 ，符合指数规律。

有预测模型

二次指数平滑法

一次指数平滑法存在明显的滞后偏差，可以再作二次指数平滑，有

设时间序列从某时期开始具有直线趋势，可得

三次指数平滑法

时间序列变动表现为二次曲线趋势时，需要用三次指数平滑法，即

预测模型为

差分指数平滑法

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

同样的，当时间序列呈现二次曲线增长时，可用二阶差分指数平滑模型预测

自适应滤波法

基本预测公式为

调整权数的公式为 其中 为学习常数， 为第 期预测误差。

步骤：

  1. 设定初始权数、学习常数和开始时期 ；
  2. 根据预测公式计算 期的预测值 ；
  3. 计算预测误差 ；
  4. 根据调整公式更新权数；
  5. 进 1，回到步骤 2，直到 超出序列的最后一个时期，迭代结果。

季节性时间序列预测

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

  1. 计算每年所有季度的算数平均值 其中 分别表示年份、季度的序号；

  2. 计算同季度的算数平均值

  3. 计算季度系数 ​；

  4. 预测下一年的年加权平均

  5. 得到各季度的预测值

ARMA 时间序列

设随机序列 满足

  • 常数​，
  • 与 ​ 无关，

则称 为平稳序列。

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

  • AR 序列

其中 是均值为 0，方差为 的平稳白噪声，称 是阶数为 的自回归序列，记 。

  • MA 序列

其中 是均值为 0，方差为 的平稳白噪声，称 是阶数为 的移动平均序列，记 。

于是有

其中 是均值为 0，方差为 的平稳白噪声，称 是阶数为 的自回归移动平均序列，记 。

模型构建

AIC 准则：选 使得 其中 是样本容量， 是 的估计值。​

参数估计

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

模型预测

记 的估计量为 ，对于 AR 序列 对于 MA 序列，有 令 ，得 可取初值 ​，当 ​ 充分大后，初始误差的影响逐渐消失。

因此对于 ARMA 序列，令 有递推式

其中 满足 ，可令初值 。

ARIMA 模型

对差分运算后得到的平稳序列用 ARMA 模型拟合，称 模型，其中 分别表示 的项数、差分阶数、 的项数。

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

输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#! /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 年的太阳黑子数量，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#! /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()

输出如下：

1
2
3
4
5
6
7
8
9
10
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 模型作预测，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#! /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)

输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
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