ARIMA时间序列预测模型详细讲解+Python案例演示
🤵♂️ 个人主页：@艾派森的个人主页
✍🏻作者简介：Python学习者
🐋 希望大家多多支持，我们一起进步！😄
如果文章对你有帮助的话，
欢迎评论 💬点赞👍🏻 收藏 📂加关注+
目录
模型引出
ARIMA模型
自回归模型（AR）
移动平均模型（MA）
自回归移动平均模型（ARMA）
差分自回归移动平均模型（ARIMA）
ARIMA模型的建模步骤
股票价格预测
源代码
模型引出
时间序列
时间序列也称动态序列，是指将某种现象的指标数值按照时间顺序排列而成的数值序列。时间序列分析大致可分成三大部分，分别是描述过去、分析规律和预测未来，本文将主要介绍时间序列分析中常用ARIMA模型。
时间序列数据
对同一对象在不同时间连续观察所取得的数据，它具备两个要素，第一个要素是时间要素，第二个要素是数值要素。例如：
- 从出生到现在，你的体重的数据(每年生日称重一次)
- 中国历年来GDP的数据
- 在某地方每隔一小时测得的温度数据
时间序列根据时间和数值性质的不同，可以分为时期时间序列和时点时间序列。
- 时期序列中，数值要素反映现象在一定时期内发展的结果
- 时点序列中，数值要素反映现象在一定时点上的瞬间水平
区分时期和时点时间序列
- 1)从出生到现在，你的体重的数据(每年生日称重一次)
- 2)中国历年来GDP的数据
- 3)在某地方每隔一小时测得的温度数据
1)和3)是时点时间序列;2)是时期时间序列
时期序列可加，时点序列不可加
时期序列中的观测值反映现象在一段时期内发展过程的总量，不同时期的观测值可以相加相加结果表明现象在更长一段时间内的活动总量；而时点序列中的观测值反映现象在某一瞬间上所达到的水平，不同时期的观测值不能相加，相加结果没有实际意义。
ARIMA模型
ARIMA模型全称叫差分整合移动平均自回归模型，AR代表自回归模型，MA代表移动平均模型，中间的I代表差分次数。
自回归模型（AR）
描述当前值和历史值之间的关系，用变量自身的历史数据对自身进行预测，其必须要满足平稳性要求，只适用于预测与自身前期相关的现象(时间序列的自相关性)
移动平均模型（MA）
移动平均模型关注的是自回归模型中误差项的累计
即时间序列当前值与历史值没有关系，而只依赖于历史白噪声的线性组合
移动平均法能有效地消除预测中的随机波动
自回归移动平均模型（ARMA）
自回归与移动平均的结合
该式表明：
- 一个随机时间序列可以通过一个自回归移动平均模型来表示，即该序列可以由其自身的过去或滞后值以及随机扰动项来解释。
- 如果该序列是平稳的，即它的行为并不会随着时间的推移而变化，那么我们就可以通过该序列过去的行为来预测未来。
差分自回归移动平均模型（ARIMA）
将自回归模型(AR)、移动平均模型(MA)和差分法结合，我们就得到了差分自回归移动平均模型
ARIMA(p、d、q)
p是自回归项，q为移动平均项数，d为时间序列成为平稳时所做的差分次数
原理：将非平稳时间序列转化为平稳时间序列然后将因变量仅对它的滞后值以及随机误差项的
现值和滞后值进行回归所建立的模型。
ARIMA模型的建模步骤
- 对序列绘图，进行平稳性检验，观察序列是否平稳；对于非平稳时间序列要先进行d阶差分，转化为平稳时间序列；
- 经过第一步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数(ACF)和偏自相关系数(PACF)，通过对自相关图和偏自相关图的分析，得到最佳的阶数p、q;
- 由以上得到的d、q、p，得到 ARIMA 模型。然后开始对得到的模型进行模型检验。
股票价格预测
已知一个上市公司（阿里巴巴）一段时期的开盘价，最高价，最低价，收盘价等信息，要求建立模型，预测股价。
Python版本:3.9
代码编辑器：jupyter notebook
首先导入我们的数据集
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')df = pd.read_csv('BABA.csv',index_col='Date')
df
在股票价格预测问题中，一般是以收盘价为准，这里我们先绘制一下原始数据集中的收盘价线图
df.index = pd.to_datetime(df.index)
plt.figure(figsize=(15,6))
plt.plot(df['Close'])
plt.show()
从图中可以看出，在2022年之前，价格变动的幅度较大，为了提高模型的准确率，这里我们只选取2022年以后的数据进行研究。
new_df = df.loc['2022':]
new_df
提取出收盘价Close列
# 提取Close列
close_data = new_df['Close']
close_data
拆分数据为训练集和测试集，其中将2022年和2023年作为训练集，2024年的数据作为测试集
# 划分训练集和测试集
train = close_data.loc['2022':'2024'] # 2022年-2023年作为训练集
test = close_data.loc['2024':] # 2024年之后作为测试集绘制一下训练集和测试集的线图 ，其中蓝色的是训练集，橙色的是测试集
plt.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示
plt.rcParams['axes.unicode_minus'] = False #解决符号无法显示
plt.figure(figsize=(15,6))
plt.plot(train,label='训练集')
plt.plot(test,label='测试集')
plt.legend()
plt.show()
平稳性
平稳性就是要求经由样本时间序列所得到的拟合曲线在未来的一段时间内仍然能够按照现有的形态延续下去。
平稳性要求序列的均值和方差不发生明显变化。
- 严平稳:序列所有的统计性质(期望，方差)都不会随着时间的推移而发生变化
- 宽平稳:期望与相关系数(依赖性)不变，就是说t时刻的值X依赖于过去的信息
实际数据大致上都是宽平稳.
平稳性对于我们分析时间序列至关重要。如果一个时间序列不是平稳的，通常需要通过差分的方式将其转化为平稳时间序列。
差分法实现
时间序列在t和t-1时刻的差值，将非平稳序列变平稳。Δyx= y(x+ 1)-y(x),(x= 0,1,2,..)
比如一组数列[0,1,2,3,4,5,6,7]，进行差分后就会得到新数列[1,1,1,1,1,1,1]
new_df['diff_1'] = new_df['Close'].diff(1) # 1阶差分
new_df['diff_2'] = new_df['diff_1'].diff(1) # 2阶差分
fig = plt.figure(figsize=(12,10))
# 源数据
ax1 = fig.add_subplot(311)
ax1.plot(new_df['Close'])
# 1阶差分
ax2 = fig.add_subplot(312)
ax2.plot(new_df['diff_1'])
# 2阶差分
ax3 = fig.add_subplot(313)
ax3.plot(new_df['diff_2'])
plt.show()
从图中看出，一阶差分的结果就已经很不错了，接近于平稳序列。但如果你觉得2阶差分更好，那你也可以做2次差分。
自相关系数（ACF）
有序的随机变量序列与其自身相比较。自相关系数反映了统一序列在不同时序的取值之间的相关性，对于时间序列yt，yt与yt-k的相关系数称为y间隔k的自相关系数。
偏自相关系数（PACF）
对于一个平稳AR(p)模型，求出滞后k自相关系数p(k)时，实际上得到并不是x(t)与x(t-k)之间单纯的相关关系。
因为x(t)同时还会受到中间k-1个随机变量x(t-1),x(t-2),…,x(t-k+1)的影响，而这k-1个随机变量又都和x(t-k)具有相关关系，所以自相关系数里面实际掺杂了其他变量对x(t)与x(t-k)的影响。
为了能单纯测度x(t-k)对x(t)的影响，引进偏自相关系数(PACF)的概念。对于平稳时间序列{x(t)}，所谓滞后k偏自相关系数指在剔除了中间k-1个随机变量x(t-1),x(t-2),.,x(t-k+ 1)的干扰之后，x(t-k)对x(t)影响的相关程度。
ADF检验
对于一个时间序列，如何确定它是否满足平稳性要求?通常采用图检验法(通过时间序列趋势图或者自相关函数图判断)或ADF检验
ADF大致的思想就是基于随机游走(不平稳的一个特殊序列)的，对其进行回归，如果发现p=1，说明序列满足随机游走，就是非平稳的。
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller as ADF# 计算原始序列、一阶差分序列、二阶差分序列的单位根检验结果
new_df['diff_1'] = new_df['diff_1'].fillna(0)
new_df['diff_2'] = new_df['diff_2'].fillna(0)timeseries_adf = ADF(new_df['Close'].tolist())
timeseries_diff1_adf = ADF(new_df['diff_1'].tolist())
timeseries_diff2_adf = ADF(new_df['diff_2'].tolist())# 打印单位根检验结果
print('timeseries_adf：',timeseries_adf)
print('timeseries_diff1_adf：',timeseries_diff1_adf)
print('timeseries_diff2_adf：',timeseries_diff2_adf)运行结果如下：
timeseries_adf： (-2.896948320582682, 0.045704018855239004, 0, 583, {'1%': -3.441616425652826, '5%': -2.866510292004876, '10%': -2.5694170294713863}, 2916.325614092606)
timeseries_diff1_adf： (-24.914648751916296, 0.0, 0, 583, {'1%': -3.441616425652826, '5%': -2.866510292004876, '10%': -2.5694170294713863}, 2926.0720185762493)
timeseries_diff2_adf： (-9.67035487508312, 1.281822983964346e-16, 18, 565, {'1%': -3.4419770063102213, '5%': -2.866669060591297, '10%': -2.5695016312945413}, 2977.896861773869)ADF检验的结果共有五个参数：
第一个值:表示Test Statistic，即T检验，表示T统计量，假设检验值
第二个值:p-value，即p值，表示T统计量对应的概率值
第三/四个值:Lags Used，即表示延迟和测试的次数
第五个参数{10%: xxx,'1%':xxx,'5%:xxx}:不同程度拒绝原假设的统计值
如何确定该序列是否平稳呢?
- 1%、5%、10%不同程度拒绝原假设的统计值和 ADF 假设检验值比较，ADF 假设检验值同时小于1%、5%、10%即说明非常好地拒绝该假设
- P-value是否非常接近0
这里我们重点看第二个值，p值越接近0越好，如果为0则表明拒绝随机游走，说明序列是平稳的。观察结果，我们发现一阶差分的结果p值已经为0，所有我们只需要一次差分即可。
平稳性检验
import statsmodels.api as sm
# 绘制
fig = plt.figure(figsize=(12,7))ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(train,lags=40,ax=ax1)
ax1.xaxis.set_ticks_position('bottom')ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(train,lags=40,ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
plt.show()
上图是训练集的ACF和PACF图，由图形可以看出，大部分的值都没有落在置信区间内，也间接说明数据集需要进行差分处理。接着我们继续绘制一阶差分后的ACF和PACF图
import statsmodels.api as sm
# 绘制
fig = plt.figure(figsize=(12,7))ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(new_df['diff_1'],lags=20,ax=ax1)
ax1.xaxis.set_ticks_position('bottom')ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(new_df['diff_1'],lags=20,ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
plt.show()
上图是一阶差分后训练集的ACF和PACF图，由图形可以看出，大部分的值都落在了置信区间内，说明一阶差分序列为平稳序列，无需再差分。
确定p、q
拖尾和截尾:拖尾指序列以指数率单调递减或震荡衰减，而截尾指序列从某个时点变得非常小
截尾(出现以下情况，通常视为(偏)自相关系数d阶 截尾)
- 1)在最初的d阶明显大于2倍标准差范围
- 2)之后几乎95%的(偏)自相关系数都落在2倍标准差范围以内
- 3)且由非零自相关系数衰减为在零附近小值波动的过程非常突然
拖尾(出现以下情况，通常视为(偏)自相关系数拖尾)
- 1)如果有超过5%的样本(偏)自相关系数都落入2倍标准差范围之外
- 2)或者是由显著非0的(偏)自相关系数衰减为小值波动的过程比较缓慢或非常连续
从训练集的ACF和PACF图确认，ACF拖尾，PACF1阶截尾，所以选择ARIMA(1，1，0)模型。
通过拖尾和截尾对模型定阶，具有很强的主观性。回顾一下我们对于模型参数估计得方法，是通过对损失和正则项的加权评估。我们在参数选择的时候，需要平衡预测误差与模型复杂度。我们可以根据信息准则函数法，来确定模型的阶数。这里介绍AIC、BIC准则。
AIC 准则全称是最小化信息量准则(Akaike Information Criterion):
AIC=-2ln(L)+2K ，其中L 表示模型的极大似然函数，K表示模型参数个数·
AIC准则存在一定的不足。当样本容量很大时，在 AIC 准则中拟合误差提供的信息就要受到样本容量的放大，而参数个数的惩罚因子却和样本容量没关系(一直是2)，因此当样本容量很大时使用 AIC 准则的模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多.
BIC(Bayesian InformationCriterion)贝叶斯信息准则弥补了 AIC 的不足:
BIC=-2ln(L)+Kln(n)其中 n 表示样本容量。
显然，这两个评价指标越小越好。
我们通过网格搜索，确定 AIC、BIC 最优的模型(p、q)
我们以BIC准则为例，确定p,q的取值范围为[0,5]，通过循环网格搜索所有组合的BIC的值
# 遍历寻找适当的参数
import itertools
import numpy as np
import seaborn as sns
# 确定pq的取值范围
d = 1
p_min = 0
q_min = 0
p_max = 5
q_max = 5# 初始化一个DataFrame来存储结果,以BIC准则
results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
for p,q in itertools.product(range(p_min,p_max+1),range(q_min,q_max+1)):if p==0 and q==0:results_bic.loc['AR{}'.format(p),'MA{}'.format(q)] = np.nancontinuetry:model = sm.tsa.ARIMA(train,order=(p,d,q))results = model.fit()results_bic.loc['AR{}'.format(p),'MA{}'.format(q)] = results.bicexcept:continueresults_bic
绘制成热力图
results_bic = results_bic[results_bic.columns].astype(float)
# 绘制热力图
fig,ax = plt.subplots(figsize=(10,8))
ax = sns.heatmap(results_bic,mask=results_bic.isnull(),ax=ax,annot=True,fmt='.3f')
ax.set_title('BIC')
plt.show()
从图中可以看到，BIC最小值的组合为('AR1'，'MA0')
也可以通过如下代码进行查看
我们也可以使用 AIC和 BIC 准则对训练数据 train 进行 ARMA 模型阶数的选择，这里限制了自回归项(AR)和移动平均项(MA)的最大阶数，将其设置为8
train_results = sm.tsa.arma_order_select_ic(train,ic=['aic','bic'],trend='n',max_ar=8,max_ma=8)
print('AIC',train_results.aic_min_order)
print('BIC',train_results.bic_min_order)
最终的结果都是在说明p取1，q取0
模型检验
检验参数估计的显著性(t检验)
检验残差序列的随机性，即残差之间是独立的
残差序列的随机性可以通过自相关函数法来检验，即做残差的自相关函数图
# 根据以上求得
p = 1
d = 1
q = 0
# 初始化模型
model = sm.tsa.ARIMA(train,order=(p,d,q))
# 训练模型
results = model.fit()
# 获取残差
resid = results.resid
# 绘制残差的ACF图
fig,ax = plt.subplots(figsize=(12,5))
ax = sm.graphics.tsa.plot_acf(resid,lags=40,ax=ax)
plt.show()
从ACF图中可以看出，大部分值都非常接近0，说明残差之间独立性比较高，模型效果不错
打印模型信息
results.summary()
模型预测
我们以测试集的时间轴进行预测
predict_ = results.predict(start='2024',end='2024-04-30')
predict_
将训练集、测试集、预测值进行对比
plt.figure(figsize=(15,6))
plt.plot(train,label='训练集')
plt.plot(test,label='测试集')
plt.plot(predict_,label='预测值')
plt.legend()
plt.show()
蓝色是训练集，橙色是测试集，绿色是预测值，可以看到预测值与测试集基本吻合，我们放大看看。
plt.figure(figsize=(12,6))
plt.plot(test,label='真实值')
plt.plot(predict_,label='预测值')
plt.legend()
plt.show()
可以看到预测值与真实值是相接近的，趋势也基本保持一致、
最后我们预测未来5天股票收盘价
results.forecast(5) # 预测未来5天
源代码
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')df = pd.read_csv('BABA.csv',index_col='Date')
df
df.index = pd.to_datetime(df.index)
plt.figure(figsize=(15,6))
plt.plot(df['Close'])
plt.show()
new_df = df.loc['2022':]
new_df
# 提取Close列
close_data = new_df['Close']
close_data
# 划分训练集和测试集
train = close_data.loc['2022':'2024'] # 2022年-2023年作为训练集
test = close_data.loc['2024':] # 2024年之后作为测试集
plt.rcParams['font.sans-serif'] = ['SimHei'] #解决中文显示
plt.rcParams['axes.unicode_minus'] = False #解决符号无法显示
plt.figure(figsize=(15,6))
plt.plot(train,label='训练集')
plt.plot(test,label='测试集')
plt.legend()
plt.show()
差分法
new_df['diff_1'] = new_df['Close'].diff(1) # 1阶差分
new_df['diff_2'] = new_df['diff_1'].diff(1) # 2阶差分
fig = plt.figure(figsize=(12,10))
# 源数据
ax1 = fig.add_subplot(311)
ax1.plot(new_df['Close'])
# 1阶差分
ax2 = fig.add_subplot(312)
ax2.plot(new_df['diff_1'])
# 2阶差分
ax3 = fig.add_subplot(313)
ax3.plot(new_df['diff_2'])
plt.show()
ADF检验
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller as ADF# 计算原始序列、一阶差分序列、二阶差分序列的单位根检验结果
new_df['diff_1'] = new_df['diff_1'].fillna(0)
new_df['diff_2'] = new_df['diff_2'].fillna(0)timeseries_adf = ADF(new_df['Close'].tolist())
timeseries_diff1_adf = ADF(new_df['diff_1'].tolist())
timeseries_diff2_adf = ADF(new_df['diff_2'].tolist())# 打印单位根检验结果
print('timeseries_adf：',timeseries_adf)
print('timeseries_diff1_adf：',timeseries_diff1_adf)
print('timeseries_diff2_adf：',timeseries_diff2_adf)
参数确定
import statsmodels.api as sm
# 绘制
fig = plt.figure(figsize=(12,7))ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(train,lags=40,ax=ax1)
ax1.xaxis.set_ticks_position('bottom')ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(train,lags=40,ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
plt.show()
import statsmodels.api as sm
# 绘制
fig = plt.figure(figsize=(12,7))ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(new_df['diff_1'],lags=20,ax=ax1)
ax1.xaxis.set_ticks_position('bottom')ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(new_df['diff_1'],lags=20,ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
plt.show()
模型建立
# 遍历寻找适当的参数
import itertools
import numpy as np
import seaborn as sns
# 确定pq的取值范围
d = 1
p_min = 0
q_min = 0
p_max = 5
q_max = 5# 初始化一个DataFrame来存储结果,以BIC准则
results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
for p,q in itertools.product(range(p_min,p_max+1),range(q_min,q_max+1)):if p==0 and q==0:results_bic.loc['AR{}'.format(p),'MA{}'.format(q)] = np.nancontinuetry:model = sm.tsa.ARIMA(train,order=(p,d,q))results = model.fit()results_bic.loc['AR{}'.format(p),'MA{}'.format(q)] = results.bicexcept:continue
results_bic
results_bic = results_bic[results_bic.columns].astype(float)
# 绘制热力图
fig,ax = plt.subplots(figsize=(10,8))
ax = sns.heatmap(results_bic,mask=results_bic.isnull(),ax=ax,annot=True,fmt='.3f')
ax.set_title('BIC')
plt.show()
results_bic.stack().idxmin()
利用模型取p和q的最优值
train_results = sm.tsa.arma_order_select_ic(train,ic=['aic','bic'],trend='n',max_ar=8,max_ma=8)
print('AIC',train_results.aic_min_order)
print('BIC',train_results.bic_min_order)
模型检验
# 根据以上求得
p = 1
d = 1
q = 0
# 初始化模型
model = sm.tsa.ARIMA(train,order=(p,d,q))
# 训练模型
results = model.fit()
# 获取残差
resid = results.resid
# 绘制残差的ACF图
fig,ax = plt.subplots(figsize=(12,5))
ax = sm.graphics.tsa.plot_acf(resid,lags=40,ax=ax)
plt.show()
results.summary()
模型预测
predict_ = results.predict(start='2024',end='2024-04-30')
predict_
plt.figure(figsize=(15,6))
plt.plot(train,label='训练集')
plt.plot(test,label='测试集')
plt.plot(predict_,label='预测值')
plt.legend()
plt.show()
plt.figure(figsize=(12,6))
plt.plot(test,label='真实值')
plt.plot(predict_,label='预测值')
plt.legend()
plt.show()
results.forecast(5) # 预测未来5天
资料获取，更多粉丝福利，关注下方公众号获取
