开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> Python知识库 -> Python数学建模系列（十）：时间序列 -> 正文阅读

[Python知识库]Python数学建模系列（十）：时间序列

前言

Hello！小伙伴！
非常感谢您阅读海轰的文章，倘若文中有错误的地方，欢迎您指出～
?
自我介绍 ?(?ˊ?ˋ)?
昵称：海轰
标签：程序猿｜C++选手｜学生
简介：因C语言结识编程，随后转入计算机专业，有幸拿过一些国奖、省奖…已保研。目前正在学习C++/Linux/Python
学习经验：扎实基础 + 多做笔记 + 多敲代码 + 多思考 + 学好英语！
?
初学Python 小白阶段
文章仅作为自己的学习笔记用于知识体系建立以及复习
题不在多学一题懂一题
知其然知其所以然！

往期文章

Python数学建模系列（一）：规划问题之线性规划

Python数学建模系列（二）：规划问题之整数规划

Python数学建模系列（三）：规划问题之非线性规划

Python数学建模系列（四）：数值逼近

Python数学建模系列（五）：微分方程

Python数学建模系列（六）：蒙特卡洛算法

Python数学建模系列（七）：差分

Python数学建模系列（八）：图论

Python数学建模系列（九）：回归

1 JetRail高铁乘客量预测——7种时间序列方法

内容简介

时间序列预测在日常分析中常会用到，是重要的时序数据处理方法。

高铁客运量预测

假设要解决一个时序问题：根据过往两年的数据（2012 年 8 月至 2014 年 8月），需要用这些数据预测接下来 7 个月的乘客数量。

数据获取：获得2012-2014两年每小时乘客数量

import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt
 
df = pd.read_csv('../profile/train2.csv')
df.head()

运行结果
在这里插入图片描述
数据集处理：（以每天为单位构造和聚合数据集）

从2012年8月—2013年12月的数据中构造一个数据集
创建train and test文件用于建模。前14个月（2012年8月—2013年10月）用作训练数据，后两个月（2013年11月—2013年12月）用作测试数据。
以每天为单位聚合数据集

import pandas as pd
import matplotlib.pyplot as plt
 
df = pd.read_csv('../profile/train2.csv',nrows=11856)

train = df[0:10392]
test = df[10392:]
 
df['Timestamp'] = pd.to_datetime(df['Datetime'], format='%d-%m-%Y %H:%M')  # 4位年用Y，2位年用y
df.index = df['Timestamp']
df = df.resample('D').mean() #按天采样，计算均值
 
train['Timestamp'] = pd.to_datetime(train['Datetime'], format='%d-%m-%Y %H:%M')
train.index = train['Timestamp']
train = train.resample('D').mean() 
 
test['Timestamp'] = pd.to_datetime(test['Datetime'], format='%d-%m-%Y %H:%M')
test.index = test['Timestamp']
test = test.resample('D').mean()
 

train.Count.plot(figsize=(15,8), title= 'Daily Ridership', fontsize=14)
test.Count.plot(figsize=(15,8), title= 'Daily Ridership', fontsize=14)
plt.show()

运行结果

在这里插入图片描述

1.1 朴素法

如下图所示， y 轴表示物品的价格，x 轴表示时间（天）
在这里插入图片描述
如果数据集在一段时间内都很稳定，我们想预测第二天的价格，可以取前面一天的价格，预测第二天的值。这种假设第一个预测点和上一个观察点相等的预测方法就叫朴素法。即
$\quad \hat y_{t+1} = y_t$

Demo代码

dd = np.asarray(train['Count'])
y_hat = test.copy()
y_hat['naive'] = dd[len(dd) - 1]
plt.figure(figsize=(12, 8))
plt.plot(train.index, train['Count'], label='Train')
plt.plot(test.index, test['Count'], label='Test')
plt.plot(y_hat.index, y_hat['naive'], label='Naive Forecast')
plt.legend(loc='best')
plt.title("Naive Forecast")
plt.show()

运行结果

在这里插入图片描述
朴素法并不适合变化很大的数据集，最适合稳定性很高的数据集。我们计算下均方根误差，检查模型在测试数据集上的准确率：

from sklearn.metrics import mean_squared_error
from math import sqrt
 
rms = sqrt(mean_squared_error(test['Count'], y_hat['naive']))
print(rms)

# 43.91640614391676

运行结果

在这里插入图片描述
最终均方根误差RMS为：43.91640614391676

1.2 简单平均法

如下图所示， y 轴表示物品的价格，x 轴表示时间（天）
在这里插入图片描述
物品价格会随机上涨和下跌，平均价格会保持一致。我们经常会遇到一些数据集，虽然在一定时期内出现小幅变动，但每个时间段的平均值确实保持不变。这种情况下，我们可以预测出第二天的价格大致和过去天数的价格平均值一致。这种将预期值等同于之前所有观测点的平均值的预测方法就叫简单平均法。即 $\quad \hat y_{x+1} = \frac{1}{x} \sum^x_{i=1}y_i$

y_hat_avg = test.copy()
y_hat_avg['avg_forecast'] = train['Count'].mean()
plt.figure(figsize=(12,8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['avg_forecast'], label='Average Forecast')
plt.legend(loc='best')
plt.show()

运行结果

在这里插入图片描述
简单平均法用之前全部已知的值计算出它们的平均值，将它作为要预测的下一个值。当然这不会很准确，但这种预测方法在某些情况下效果是最好的。

from sklearn.metrics import mean_squared_error 
from math import sqrt 
rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['avg_forecast'])) 
print(rms)

# 109.88526527082863

运行结果

在这里插入图片描述
这种模型并没有改善准确率。因此我们可以从中推断出当每个时间段的平均值保持不变时，这种方法的效果才能达到最好。虽然朴素法的准确率高于简单平均法，但这并不意味着朴素法在所有的数据集上都比简单平均法好。

1.3 移动平均法

如下图所示， y 轴表示物品的价格，x 轴表示时间（天）
在这里插入图片描述
物品价格在一段时间内大幅上涨，但后来又趋于平稳。我们也经常会遇到这种数据集，比如价格或销售额某段时间大幅上升或下降。如果我们这时用之前的简单平均法，就得使用所有先前数据的平均值，但在这里使用之前的所有数据是说不通的，因为用开始阶段的价格值会大幅影响接下来日期的预测值。因此，我们只取最近几个时期的价格平均值。很明显这里的逻辑是只有最近的值最要紧。这种用某些窗口期计算平均值的预测方法就叫移动平均法。

计算移动平均值涉及到一个有时被称为“滑动窗口”的大小值p。使用简单的移动平均模型，我们可以根据之前数值的固定有限数p的平均值预测某个时序中的下一个值。这样，对于所有的 𝑖>𝑝：

$\quad \hat y_l = \frac{1}{p}(y_{i-1} + y_{i-2} + ... + y_{i-p})$

移动平均法实际很有效，特别是当你为时序选择了正确的p值时

y_hat_avg = test.copy() 
y_hat_avg['moving_avg_forecast'] = train['Count'].rolling(60).mean().iloc[-1] 
plt.figure(figsize=(16,8)) 
plt.plot(train['Count'], label='Train') 
plt.plot(test['Count'], label='Test') 
plt.plot(y_hat_avg['moving_avg_forecast'], label='Moving Average Forecast') 
plt.legend(loc='best') 
plt.show()

运行结果

在这里插入图片描述

from sklearn.metrics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['moving_avg_forecast']))
print(rms)

此方法计算出来的均方根差为：46.72840725106963

我们可以看到，对于这个数据集，朴素法比简单平均法和移动平均法的表现要好。此外，我们还可以试试简单指数平滑法，它比移动平均法的一个进步之处就是相当于对移动平均法进行了加权。在上文移动平均法可以看到，我们对“p”中的观察值赋予了同样的权重。但是我们可能遇到一些情况，比如“p”中每个观察值会以不同的方式影响预测结果。将过去观察值赋予不同权重的方法就叫做加权移动平均法。加权移动平均法其实还是一种移动平均法，只是“滑动窗口期”内的值被赋予不同的权重，通常来讲，最近时间点的值发挥的作用更大了。

$\quad \hat y_l = \frac{1}{m}(w_1y_{i-1} + w_2y_{i-2} + w_3y_{i-3} + ... + w_my_{i-m})$

这种方法并非选择一个窗口期的值，而是需要一列权重值（相加后为1）。例如，如果我们选择[0.40, 0.25, 0.20, 0.15]作为权值，我们会为最近的4个时间点分别赋给40%，25%，20%和15%的权重。

1.4 简单指数平滑法

我们注意到简单平均法和加权移动平均法在选取时间点的思路上存在较大的差异。我们就需要在这两种方法之间取一个折中的方法，在将所有数据考虑在内的同时也能给数据赋予不同非权重。例如，相比更早时期内的观测值，它会给近期的观测值赋予更大的权重。按照这种原则工作的方法就叫做简单指数平滑法。它通过加权平均值计算出预测值，其中权重随着观测值从早期到晚期的变化呈指数级下降，最小的权重和最早的观测值相关：

在这里插入图片描述
其中0≤α≤1是平滑参数。对时间点T+1的单步预测值是时序所有观测值的加权平均数。权重下降的速率由参数α控制，预测值 $\hat y_x$ 是时间点T的单步预测值与 $\alpha) \hat y_{x-1}$ 的和。因而写作

在这里插入图片描述

所以本质上，我们是用两个权重α和1?α得到一个加权移动平均值。我们可以看到 $\hat y_{x-1}$ 和1?α相乘，让表达式呈递进形式，这也是该方法被称为“指数”的原因。时间 t+1 处的预测值为最近观测值yt和最近预测值 $\hat y_{t|t-1}$ 之间的加权平均值。

from statsmodels.tsa.api import SimpleExpSmoothing

y_hat_avg = test.copy()
fit = SimpleExpSmoothing(np.asarray(train['Count'])).fit(smoothing_level=0.6, optimized=False)
y_hat_avg['SES'] = fit.forecast(len(test))
plt.figure(figsize=(16, 8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['SES'], label='SES')
plt.legend(loc='best')
plt.show()

运行结果

在这里插入图片描述

from sklearn.metrics import mean_squared_error
from math import sqrt

rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['SES']))
print(rms)

均方根差为：43.357625225228155

1.5 霍尔特线性趋势法

如下图所示， y 轴表示物品的价格，x 轴表示时间（天）

如果物品的价格是不断上涨的（见上图），我们上面的方法并没有考虑这种趋势，即我们在一段时间内观察到的价格的总体模式。在上图例子中，我们可以看到物品的价格呈上涨趋势。虽然上面这些方法都可以应用于这种趋势，但我们仍需要一种方法可以在无需假设的情况下，准确预测出价格趋势。这种考虑到数据集变化趋势的方法就叫做霍尔特线性趋势法。

每个时序数据集可以分解为相应的几个部分：趋势（Trend），季节性(Seasonal)和残差(Residual)。任何呈现某种趋势的数据集都可以用霍尔特线性趋势法用于预测。

Demo代码

import statsmodels.api as sm

sm.tsa.seasonal_decompose(train['Count']).plot()
result = sm.tsa.stattools.adfuller(train['Count'])
plt.show()

运行结果

在这里插入图片描述
我们从图中可以看出，该数据集呈上升趋势。因此我们可以用霍尔特线性趋势法预测未来价格。该算法包含三个方程：一个水平方程，一个趋势方程，一个方程将二者相加以得到预测值 $\hat y$

我们在上面算法中预测的值称为水平（level）。正如简单指数平滑一样，这里的水平方程显示它是观测值和样本内单步预测值的加权平均数，趋势方程显示它是根据 ?(t)??(t?1) 和之前的预测趋势 b(t?1) 在时间t处的预测趋势的加权平均值。

我们将这两个方程相加，得出一个预测函数。我们也可以将两者相乘而不是相加得到一个乘法预测方程。当趋势呈线性增加和下降时，我们用相加得到的方程；当趋势呈指数级增加或下降时，我们用相乘得到的方程。实践操作显示，用相乘得到的方程，预测结果会更稳定，但用相加得到的方程，更容易理解。
在这里插入图片描述
Demo代码

from statsmodels.tsa.api import Holt

y_hat_avg = test.copy()

fit = Holt(np.asarray(train['Count'])).fit(smoothing_level=0.3, smoothing_slope=0.1)
y_hat_avg['Holt_linear'] = fit.forecast(len(test))

plt.figure(figsize=(16, 8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['Holt_linear'], label='Holt_linear')
plt.legend(loc='best')
plt.show()

运行结果

在这里插入图片描述

from sklearn.metrics import mean_squared_error
from math import sqrt

rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['Holt_linear']))
print(rms)
# 43.056259611507286

均方根误差为：43.056259611507286

1.6 Holt-Winters季节性预测模型

在应用这种算法前，我们先介绍一个新术语。假如有家酒店坐落在半山腰上，夏季的时候生意很好，顾客很多，但每年其余时间顾客很少。因此，每年夏季的收入会远高于其它季节，而且每年都是这样，那么这种重复现象叫做“季节性”（Seasonality）。如果数据集在一定时间段内的固定区间内呈现相似的模式，那么该数据集就具有季节性。

在这里插入图片描述
我们之前讨论的5种模型在预测时并没有考虑到数据集的季节性，因此我们需要一种能考虑这种因素的方法。应用到这种情况下的算法就叫做Holt-Winters季节性预测模型，它是一种三次指数平滑预测，其背后的理念就是除了水平和趋势外，还将指数平滑应用到季节分量上。
在这里插入图片描述

其中 s 为季节循环的长度，0≤α≤ 1, 0 ≤β≤ 1 ， 0≤γ≤ 1。水平函数为季节性调整的观测值和时间点t处非季节预测之间的加权平均值。趋势函数和霍尔特线性方法中的含义相同。季节函数为当前季节指数和去年同一季节的季节性指数之间的加权平均值。

在本算法，我们同样可以用相加和相乘的方法。当季节性变化大致相同时，优先选择相加方法，而当季节变化的幅度与各时间段的水平成正比时，优先选择相乘的方法。

from statsmodels.tsa.api import ExponentialSmoothing

y_hat_avg = test.copy()
fit1 = ExponentialSmoothing(np.asarray(train['Count']), seasonal_periods=7, trend='add', seasonal='add', ).fit()
y_hat_avg['Holt_Winter'] = fit1.forecast(len(test))
plt.figure(figsize=(16, 8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['Holt_Winter'], label='Holt_Winter')
plt.legend(loc='best')
plt.show()

运行结果

在这里插入图片描述

from sklearn.metrics import mean_squared_error
from math import sqrt

rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['Holt_Winter']))
print(rms)

# 25.26453787766697

1.7 自回归移动平均模型（ARIMA）

另一个场景的时序模型是自回归移动平均模型（ARIMA）。指数平滑模型都是基于数据中的趋势和季节性的描述，而自回归移动平均模型的目标是描述数据中彼此之间的关系。ARIMA的一个优化版就是季节性ARIMA。它像Holt-Winters季节性预测模型一样，也把数据集的季节性考虑在内。

Demo代码

import statsmodels.api as sm

y_hat_avg = test.copy()
fit1 = sm.tsa.statespace.SARIMAX(train.Count, order=(2, 1, 4), seasonal_order=(0, 1, 1, 7)).fit()
y_hat_avg['SARIMA'] = fit1.predict(start="2013-11-1", end="2013-12-31", dynamic=True)
plt.figure(figsize=(16, 8))
plt.plot(train['Count'], label='Train')
plt.plot(test['Count'], label='Test')
plt.plot(y_hat_avg['SARIMA'], label='SARIMA')
plt.legend(loc='best')
plt.show()

运行结果

在这里插入图片描述

from sklearn.metrics import mean_squared_error
from math import sqrt

rms = sqrt(mean_squared_error(test['Count'], y_hat_avg['SARIMA']))
print(rms)

# 26.05142646431944