文章目录
ARIMA(自回归整合移动平均模型)模型是数学建模中常用到的时间序列模型,时间序列模型的数据特点是数据的观测值是与时间相关的,可能存在周期性或者季节性,通常还伴随着一定的随机性。ARIMA模型结构清晰,由自回归项(AR)、差分项(I)和移动平均项(MA)三部分组成,分别对应历史值的影响、随机误差项的影响以及过去随机误差项的影响,便于理解和解释。
一、什么是自回归项(AR)?
AR部分用于处理时间序列的自回归部分,它考虑的是过去若干个时期的观测值对当前观测值的影响。换种理解就是我们常说的滑动窗口。
以框选中的数据作为自变量从而来使用线性回归来预测下一项。
AR模型的数学表达式:
Y
t
=
c
+
φ
1
Y
t
−
1
+
φ
2
Y
t
−
2
+
.
.
.
+
φ
p
Y
t
−
p
+
ξ
t
Y_t=c+\varphi _1Y_{t-1}+\varphi _2Y_{t-2}+...+\varphi _pY_{t-p}+\xi _t
Yt=c+φ1Yt−1+φ2Yt−2+...+φpYt−p+ξt
- Y t Y_t Yt表示当前的时间序列模型
- c c c表是常数项
- ξ t \xi _t ξt表示误差项
使用AR模型表示了当前值 Y t Y_t Yt与它过去时刻的观测值有关,AR模型能很好的根据历史数据捕获趋势,但是也容易受到异常值的影响。
二、什么是移动平均(MA)?
MA部分是处理时间序列的移动平均部分,它考虑了过去的预测误差对当前值的影响。
MA模型的数学表达式:
Y
t
=
u
+
ϵ
t
+
θ
1
ϵ
t
−
1
+
θ
2
ϵ
t
−
2
+
.
.
.
+
θ
q
ϵ
t
−
q
Y_t=u+\epsilon _t+\theta _1\epsilon _{t-1}+\theta _2\epsilon _{t-2}+...+\theta _q\epsilon _{t-q}
Yt=u+ϵt+θ1ϵt−1+θ2ϵt−2+...+θqϵt−q
- Y t Y_t Yt表示当前的观测数据
- u u u表示常数项
- ϵ t \epsilon _t ϵt表示在t时间点的误差项
MA模型无法项AR模型那样捕捉数据的长期趋势,但是可以处理那些变化较大的异常值。
三、什么是差分(I)?
差分是一种数学操作,计算相邻数据的差值,能够将不平稳的数据转换成平稳的数据。如果有一个时间序列
Y
t
Y_t
Yt,那么该数据的一阶差分就为:
Δ
Y
=
Y
t
−
Y
t
−
k
\varDelta Y=Y_t-Y_{t-k}
ΔY=Yt−Yt−k
其中的k表示滞后阶数,当k=1时表示前后两个数据的差值,这样就会得到一个新的时间序列数据,这样通常得到的数据就平稳了,如果不平稳就在进行一次,进行几次差分一般我们称作几阶差分,ARIMA(p,d,q)三个参数中的d表示差分的阶数。
1、如何判断数据是否平稳?
判断数据平稳的方式有很多种,比如使用ADF检验是否存在单位根或者使用画图的方式看趋势,当然为了严谨一般使用前者。
代码实列如下:
import statsmodels.api as sm
print("ADF test: p=%f" % sm.tsa.stattools.adfuller(yourdatas)[1]) # ADF检验
2、为什么要使用差分?
在ARIMA模型中使用差分主要是为了满足模型的基本假设——时间序列数据的平稳性。平稳时间序列的特征是其均值、方差和协方差结构在时间上保持不变。对于非平稳时间序列,例如那些具有趋势、季节性或随机趋势成分的数据,直接应用自回归(AR)、移动平均(MA)或其他基于平稳序列的模型往往会得到不可靠的估计和预测结果。
四、超参数(p,d,q)的选取
1、自相关系数(ACF)和偏自相关系数(PACF)
自相关函数(ACF): 对于一个时间序列
Y
t
Y_t
Yt,t=1,2,…,n,自相关函数表达式:
A
C
F
(
k
)
=
C
o
v
(
Y
t
,
Y
t
−
k
)
V
a
r
(
Y
t
)
ACF\left( k \right) =\frac{Cov\left( Y_t,Y_{t-k} \right)}{Var\left( Y_t \right)}
ACF(k)=Var(Yt)Cov(Yt,Yt−k)
- C o v ( Y t , Y t − k ) Cov\left( Y_t,Y_{t-k} \right) Cov(Yt,Yt−k)表示 Y k Y_k Yk与 Y t − k Y_{t-k} Yt−k的协方差,其中 k k k表示滞后阶数
- V a r ( Y t ) Var\left( Y_t \right) Var(Yt)表示 Y t Y_t Yt的方差
偏自相关系数(PACF):偏自相关函数 表示在已知前 k-1 个滞后值的条件下,当前值与滞后 k 期值之间的线性关联度,可以避免中间滞后阶数的影响。PACF 的计算通常涉及递归公式或通过逆 Toeplitz 矩阵求解,其表达式相对复杂,不易直接给出封闭形式的公式,但在实际计算中,常采用Yule-Walker方程组或Levinson-Durbin算法等方法。
使用ACF和PACF确定ARIMA(p,d,q)参数中的p和q:
模型 | ACF | PACF |
---|---|---|
AR ( p p p) | 衰减趋于零(几何型或者振荡型) | p阶后截尾 |
MA( q q q) | q阶后截尾 | 衰减趋于零(几何型或振荡型) |
ARMA(p,q) | q阶后衰减趋于零(几何型或振荡型) | q阶后衰减趋于零(几何型或振荡型) |
- 截尾:落在置信区间内(95%的点都符合规则)
2、网格搜索
使用ACF和PACF后可以找到p和q的大致范围,一般只看图很难确定。一般找到使模型AIC或者BIC最小的p,q组合。
- AIC:赤池信息化准则
A I C = 2 k − 2 ln ( L ) AIC=2k-2\ln \left( L \right) AIC=2k−2ln(L) - BIC:贝叶斯信息准则
B I C = k ln ( n ) − 2 ln ( L ) BIC=k\ln \left( n \right) -2\ln \left( L \right) BIC=kln(n)−2ln(L)
其中 k k k表示参数的数量, n n n表示样本容量, L L L表示似然函数
示例代码
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy import stats
import statsmodels.api as sm
import warnings
from itertools import product
from datetime import datetime
warnings.filterwarnings('ignore')
plt.style.use('seaborn-poster')
# 中文乱码处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
data = pd.read_excel(r"C:\Users\28135\Desktop\美赛模拟\人数数据.xlsx")
# 时间戳转化为日常时间格式
data['Date_index'] = pd.to_datetime(data['Date'])
data = data.set_index('Date_index')
data.index = pd.to_datetime(data['Date'])
data = data.loc['2022-1-7':,:]
#test_data = data.loc['2022-12-1':,:]
#data = data.loc['2016-9-13':,:]
complete_index = pd.date_range(start=data.index.min(), end=data.index.max(), freq='D')
data = data.reindex(complete_index)
# 三次样条插值
# data = data.interpolate(method='spline', order=3)
# 比特币的价格变化趋势
fig = plt.figure(figsize=[15, 8])
plt.suptitle("比特币价格变化趋势, 单位 美元", fontsize=22)
plt.plot(data.Value, '-', label='By Days')
plt.legend()
plt.show()
# 单位根检验(未差分)
print("Dickey–Fuller test: p=%f" % sm.tsa.stattools.adfuller(data.Value)[1])
# 一阶差分
data['Value_diff'] = data.Value - data.Value.shift(1)
fig = plt.figure(figsize=[15, 8])
plt.suptitle("一阶差分比特币价格变化趋势, 单位 美元", fontsize=22)
plt.plot(data.Value_diff[2:], '-', label='By Days')
plt.legend()
plt.show()
# 单位根检验
print("Dickey–Fuller test: p=%f" % sm.tsa.stattools.adfuller(data.Value_diff[2:])[1])
# 自相关和偏自相关图
plt.figure(figsize=(15, 7))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(data.Value_diff[2:].values.squeeze(), lags=15, ax=ax)
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(data.Value_diff[2:].values.squeeze(), lags=15, ax=ax)
plt.tight_layout() # 自动调整当前图形的间距,以使其适合在屏幕上显示,优化布局显示
plt.show()
# 设置p和q的范围
p_range = range(10)
q_range = range(10)
# 初始化最小 AIC 和对应的 p、q
min_aic = float('inf')
best_p, best_q = 0, 0
# 循环尝试不同的 p 和 q 值
for p, q in product(p_range, q_range):
try:
# 拟合 ARIMA 模型
model = sm.tsa.ARIMA(data.Value, order=(p, 1, q))
results = model.fit()
except ValueError:
print('wrong parameters:', "p=%d, q=%d" % (p, q))
continue
# 计算 AIC
aic = results.aic
# 更新最小 AIC 和对应的 p、q
if aic < min_aic:
min_aic = aic
best_p, best_q = p, q
print(f"Best AIC: {min_aic}, Best p: {best_p}, Best q: {best_q}")
from statsmodels.tsa.arima.model import ARIMA
# Assuming 'data' is your DataFrame
model = ARIMA(data['Value'], order=(best_p, 1, best_q))
results = model.fit()
# Display the model summary
print(results.summary())
last_observation = data['Value'].iloc[-1]
# 假设您想预测未来60个时间点的值
forecast_steps = 60
# forecast_index = pd.date_range(start=data.index.min(), periods=len(data) + forecast_steps, freq='D')
# plt.figure(figsize=(15, 7))
forecast = results.get_forecast(steps=forecast_steps,alpha=0.05)
conf_int = forecast.conf_int()
# plt.plot(forecast_index, np.concatenate([data['Value'].values, forecast]), label='ARIMA预测结果', color='red', linewidth=1) # 预测结果
# plt.title('预测', fontsize=20)
# plt.legend()
# plt.show()
print(forecast)
print(conf_int)
# forecast.to_excel(r'C:\Users\28135\Desktop\美赛模拟\人数数据.xlsx')
plt.figure(figsize=(15, 7))
plt.plot(data['Value'], label='原始数据', color='blue')
plt.plot(forecast,label = '预测结果',color = 'red',linewidth = 1)
# 填充置信区间
# plt.fill_between(conf_int.index, conf_int.iloc[:, 0], conf_int.iloc[:, 1], color='r', alpha=0.05)
plt.legend()
plt.show()