1、ARIMA模型理论基础
ARIMA是差分自回归移动平均模型的引文缩写,其中AR表示的是自回归模型,MA表示的是移动平均模型,I表示的是差分。一般写成ARIMA(p,d,q),p是自回归阶数,q是移动平均阶数,d表示差分的次数。
它针对的是 平稳的时间序列模型.然而在现实生活中绝大多 数时间序列都是非平稳的。因此可以对数据进行差分,使其转化为平稳的时间序列,再用 ARIMA模型对其数据进行建模和预测。
ARIMA模型是根据过去不同时期数据的相关性,可以进行有效和精准的短期预测,它弥补了AR和MA进行预测出现的参数过多问题,在短期预测领域具有广泛的应用。
模型具体的数学表达式为:
$$ X_t=\varphi _iX_{t-1}+\varphi _2X_{t-2}+...+\varphi_pX_{t-p}+\varepsilon _t-\theta _1\varepsilon _{t-1}-...-\theta _q\varepsilon _{t-q} $$
式中,
- $X_t$(t=1,2,…,n)为信号的原始时间序列; p、q为模型的阶次;
- φ(i=1,2,…,q)为模型的参数系数;
- $\varepsilon_t$为误差值;
- 阶数i表示使非平稳序列平稳的差分阶数;
- 模型的阶数p、q则可以通过 序列间的自相关函数(ACF)以及偏相关函数(PACF)确定。
2、ARIMA建模步骤
- 首先需要对观测值序列进行平稳性检测,如果不平稳,则对其进行差分运算直到差分后的数据平稳;
- 在数据平稳后则对其进行白噪声检验,白噪声是指零均值常方差的随机平稳序列;
- 如果是平稳非白噪声序列就计算ACF(自相关系数)、PACF(偏自相关系数),进行ARMA等模型识别;
- 对已识别好的模型,确定模型参数,最后应用预测并进行误差分析。
3、ARIMA建模实战
数据集
3.1 导入模块
import sys
import os
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
from arch.unitroot import ADF
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import style
style.use('ggplot')
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.graphics.api import qqplot
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
pd.set_option('display.float_format', lambda x: '%.5f' % x)
np.set_printoptions(precision=5, suppress=True)
"""中文显示问题"""
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']
3.2 加载数据
#加载数据
data = pd.read_excel("data/data_arima.xls",index_col="年份",parse_dates=True)
data.head()
3.3 平稳性检验
#时序图:全靠肉眼的判断和人的经验。
#一阶差分
data["diff1"] = data["xt"].diff(1).dropna()
#二阶差分
data["diff2"] = data["diff1"].diff(1).dropna()
data1 = data.loc[:,["xt","diff1","diff2"]]
data1.plot(subplots=True, figsize=(18, 12),title="差分图")
时序图需要靠肉眼分辨,所以不是很有说服力
所以这里我们使用单位根检验来判断平稳性。
3.4 单位根检验
单位根检验有许多方法,最常用的是ADF检验。
ADF(Augmented Dickey-Fuller test)是Dickey-Fuller检验的增广形式,DF检验只能应用于一阶情况,当序列存在高阶的滞后相关时,可以使用ADF检验,所以说ADF是DF检验的扩展。
ADF检验的原理:ADF检验就是判断序列是否存在单位根,如果序列平稳,就不存在单位根;否则,就会存在单位根。
ADF检验的假设:H0假设就是存在单位根,如果得到的显著性检验统计量P值小于三个置信度(10%,5%,1%),则对应有(90%,95%,99%)的把握来拒绝原假设。
# 单位根检验
print("单位根检验:\n")
print(ADF(data.diff1.dropna()))
对一阶差分进行单位根检验,得到得到:1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较。
本数据中,P-value 为 0.023,接近0,ADF Test result同时小于5%、10%即说明很好地拒绝该假设,本数据中,ADF结果为-3.156,拒绝原假设,即一阶差分后数据是平稳的。
3.4 白噪声检验
白噪声检验也成为纯随机性检验,当数据是纯随机数据时,在对数据进行分析就没有任何意义了,所以拿到数据后最好对数据进行一个纯随机性检验。
可以通过statsmodels.stats.diagnostic
里面的acorr_ljungbox
实现
函数签名如下:
rom statsmodels.stats.diagnostic import acorr_ljungbox
lbvalue,pvalue= acorr_ljungbox(data, lags=1)
- lbvalue:测试的统计量
- pvalue:基于卡方分布的p统计量,pvalue<0.05,说明不是白噪声。
#白噪声检验:判断序列是否为非白噪声序列
from statsmodels.stats.diagnostic import acorr_ljungbox
# lbvalue,pvalue=acorr_ljungbox(data.diff1.dropna(),lags=[i for i in range(1,12)],boxpierce=True)
lbvalue,pvalue=acorr_ljungbox(data.diff1.dropna(),lags=1).values[0]
print('差分序列的白噪声检验结果如下:')
print('测试的统计量lbvalue:',lbvalue)
print('基于卡方分布的p统计量,若小于0.05,说明不是白噪声:',pvalue)
通过P<α,拒绝原假设,故差分后的序列是平稳的非白噪声序列,可以进行下一步建模
3.5 模型定阶
我们此时已经有了一个平稳的时间序列,然后就要构建合是的ARIMA模型,即确定ARIMA模型中的参数p和q,参数d=1是由于我们做了一阶差分序列就平稳了。
绘制平稳序列的自相关图(ACF)和偏自相关图(PACF)。
#绘制ACF和PACF图像
def show_acf_pacf(data_diff):
acf=plot_acf(data_diff,lags=20)
plt.title("ACF")
acf.show();
pacf=plot_pacf(data_diff,lags=17)
plt.title('PACF')
pacf.show()
show_acf_pacf(data.diff1.dropna())
从ACF图和PACF图中可以发现:
可以发现自相关图拖尾或一阶截尾;偏自相关图一阶截尾。
关于ARMA通用判断标准说明如下:
拖尾和截尾说明如下:
拖尾: 始终有非零取值,不会在大于某阶后就快速趋近于0(而是在0附近波动),可简单理解为无论如何都不会为0,而是在某阶之后在0附近随机变化。
截尾: 在大于某阶(k)后快速趋于0为k阶截尾,可简单理解为从某阶之后直接就变为0。
通常情况下:
- 如果说自相关图拖尾,并且偏自相关图在p阶截尾时,此模型应该为AR(p)。
- 如果说自相关图在q阶截尾并且偏自相关图拖尾时,此模型应该为MA(q)。
- 如果说自相关图和偏自相关图均显示为拖尾,那么可结合ACF图中最显著的阶数作为q值,选择PACF中最显著的阶数作为p值,最终建立ARMA(p,q)模型。
- 如果说自相关图和偏自相关图均显示为截尾,那么说明不适合建立ARMA模型。
虽然通过拖尾和截尾可以确定参数,但是具有很强的主观性。回想我们之前在参数预估的时候往往是怎么做的,不就是损失和正则项的加权么?我们这里能不能结合最终的预测误差来确定p,q的阶数呢?在相同的预测误差情况下,根据奥斯卡姆剃刀准则,模型越小是越好的。那么,平衡预测误差和参数个数,我们可以根据信息准则函数法,来确定模型的阶数。预测误差通常用平方误差即残差平方和来表示。
常用的信息准则函数法有下面几种:
AIC准则:AIC准则是由日本统计学家Akaike与1973年提出的,全称是最小化信息量准则(Akaike Information Criterion)。它是拟合精度和参数个数的加权函数
$$ AIC=2(模型参数的个数)-2ln(模型的极大似然函数) $$
BIC准则:AIC为模型选择提供了有效的规则,但也有不足之处。当样本容量很大时,在AIC准则中拟合误差提供的信息就要受到样本容量的放大,而参数个数的惩罚因子却和样本容量没关系(一直是2),因此当样本容量很大时,使用AIC准则选择的模型不收敛与真实模型,它通常比真实模型所含的未知参数个数要多。BIC(Bayesian InformationCriterion)贝叶斯信息准则是Schwartz在1978年根据Bayes理论提出的判别准则,称为SBC准则(也称BIC),弥补了AIC的不足。BIC的定义为:
$$ BIC = ln(n)(模型中参数的个数) - 2ln(模型的极大似然函数值) $$
进行ARMA参数的选择是,AIC准则和BIC准则的提出可以有效弥补根据自相关图和偏自相关图定阶的主观性,在有限的阶数范围内帮助我们寻找相对最优拟合模型。
#信息准则定阶:AIC、BIC、HQIC
def get_pq(data):
#AIC
AIC = sm.tsa.arma_order_select_ic(data, max_ar=6, max_ma=4, ic='aic')['aic_min_order']
#BIC
BIC = sm.tsa.arma_order_select_ic(data, max_ar=6, max_ma=4, ic='bic')['bic_min_order']
print('the AIC is{},\nthe BIC is{}\n'.format(AIC,BIC))
get_pq(data["xt"])
p=2,q=0,d=1。参数已经确定了。
3.6 参数估计
#参数估计
model=sm.tsa.ARIMA(data["xt"],order=(2,1,0))
result=model.fit()
print(result.summary())
现在已经建立了一个时序模型,现在我们需要找到一种方法来计算它。
3.7 模型的显著性检验
检验对象:残差序列
一个好的拟合模型应该能够提取观察值中几乎所有的样本相关信息,即残杀序列应该为白噪声序列。反之,如果残差序列为非白噪声序列,那就意味着残差序列中还残留着相关信息未被提取,这就说明拟合模型不够有效。
#模型的显著性检验
resid = result.resid#残差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
qq图显示,我们看到红色的KDE线与N(0,1)平行,这是残留物正太分布的良好指标,说明残差序列是白噪声序列,模型的信息的提取充分。
3.8 模型预测
#模型预测
pred = result.predict(start=str('1988'),end=('2000'),dynamic=True, typ='levels')
print (pred)
3.8 模型拟合效果展示
#模型拟合效果展示
plt.figure(figsize=(12, 8))
plt.xticks(rotation=45)
plt.plot(pred)
plt.plot(data.xt)
plt.show()
参考文献
论文:
[1]朱剑飞,陈文刚,宰洪涛,张轲,王新瑞.基于ARIMA与LSTM在电力负荷预测中的对比讨论[J].电气应用,2022,41(02):27-31.
[2]查华,石舢.基于ARIMA模型对江苏省GDP的预测[J].兰州文理学院学报(自然科学版),2022,36(03):33-36+54.DOI:10.13804/j.cnki.2095-6991.2022.03.010.
[3]魏书禾.基于ARIMA-SVR模型的中国生猪价格预测分析[J].科技创新与应用,2022,12(13):20-24+29.DOI:10.19981/j.CN23-1581/G3.2022.13.004.