时间序列模型ARIMA
ARIMA模型是时间序列建模的一种方法,ARIMA模型的基本思想是将非平稳时间序列经过多个差值转换成平稳时间序列,再根据差值确定预测值。
时间序列:就是将某一个指标在不同时间上的不同数值,按照时间的先后次序排列成数列。这种数列由于受到各种偶然因素的影响,往往表现出某种随机性,彼此之间存在着统计上的依赖关系。
白噪声序列,是指白噪声过程的样本实称,简称白噪声。白噪声序列的特点表现在任何两个时点的随机变量都不相关,序列中没有任何可以利用的动态规律,因此不能用历史数据对未来进行预测和推断。
一、适用范围:
- 数据非平稳的
- 不具有季节性的变化趋势
- 时间序列
二、建模步骤
1)平稳性检验。通过对时间序列的时序图、自相关系数(ACF)图或ADF单位根检验,评判其平稳性。对非平稳序列我们可通过差分或取对数方法进行处理。
要求序列的均值与方差不发生明显变化
平稳性就是要求经由样本时间序列所得到的拟合曲线,在未来的一段时间内,仍能顺着现有的形态惯性地延续下去
严平稳:数据分布不随时间变化,均值与方差一直不变
弱平稳:期望值(均值)与相关系数的(依赖性)不变,也就是说数据前后存在某种联系,不然数据没有任何规律,研究也就没有意义。
方法一:(图检验相对主观)
自相关图检验:在某一范围内波动,且波动范围有限
时序图检验 :平稳序列具有短期相关性,伴随着延迟期数增减,自相关系数会趋于零
方法二:(统计学方法,推荐)
单位根检验:不存在单位根即使平稳序列
2)平稳序列还需进行纯随机性检验,即白噪声检验(White Noise)。如果一个序列是非白噪声序列,则说明序列之间存在信息的传递,这样才有做时间序列分析的意义。白噪声检验通常使用LB统计量对序列进行卡方检验。
原假设:H0: P1= P2=… Pm=0, ∀m>=1(Pk为序列的自相关系数,1≤k≤m);
备择假设:H1:至少存在某个Pk≠0,∀m>=1,k<=m
残差正态性检验
方法一:图检验(检验纯随机性)
自相关图检验:自相关系数等于零或者接近于零
QQ图检验:大部分点都在直线上,则数据符合正态分布
方法二:D-W(杜宾沃森)检验 或 LB统计量检验
3)建立模型。
4)估计模型中未知参数的值。
5)残差检验:当P-value都显著性大于0.05,则认为残差序列为白噪声序列,即模型拟合效果较好。
6)利用拟合模型,预测序列未来的走势[8]。
确定p,d,q参数,ARIMA(p,d,q)
比如这里q=2,p=10就差不多
不过主观看图确定未免过于主观,所以引入AIC和BIC来判断各种(p,d,q)组合,从中选择最有的组合。
AIC与BIC越小越好
参数确定与模型检验代码
顺便说一句,不要过于主观的选择参数
import numpy as np, pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
#建立ARIMA模型主要的工具,如下:
from statsmodels.graphics.tsaplots import plot_pacf,plot_acf#偏自相关与自相关图
from statsmodels.tsa.stattools import adfuller as ADF#平稳性看眼
from statsmodels.stats.diagnostic import acorr_ljungbox#白噪声检验
import statsmodels.api as sm#D-W检验,一阶自相关检验
from statsmodels.graphics.api import qqplot#画QQ图,检验数据是否正太分布
from statsmodels.tsa.arima_model import ARIMA#导入模型
from scipy import stats#用于检验正太分布 K-S检验
np.set_printoptions(precision=2,suppress=True)
plt.rcParams['font.family'] = 'YaHei Consolas Hybrid' # 设置字体样式
# plt.rcParams['font.size'] = '16' # 设置字体大小
# 通过pd.resample()函数可以进行按月统计、按年统计
# Series.resample('Y').sum()#按年求和
# Series.resample('M').sum()#按月求和
# Series.resample('D').sum()#按日求和
# Series.rolling(window=2).mean()滑动窗口:每widows个数计算一下他们的均值mean
# sd.rolling(window=2).mean().plot(linestyle='-.')
# sd.plot(linestyle='--')
# plt.show()
#第一列为时间,第二列为GDP数据
fname = r'D:\Python Project_1\Data\ARIMA.txt'
data = np.loadtxt(fname=fname,delimiter='\t',skiprows=1,dtype=np.float)[:200]
sd = pd.Series(data=data,index=range(1,201))
#diff
sd = sd.diff().dropna()##在这里进行差分
#平稳性检验 sig<0.05说明序列平稳 ;反之不平稳
res = ADF(sd);print(f"序列平稳性检验Pvalue:{res[1]}");
#白噪声检验 sig<0.05说明lags阶序列是非白噪声,则满足要求
reb = acorr_ljungbox(sd,lags=1);print(f"白噪声检验Pvalue:{reb['lb_pvalue'].values}");
# 通过以上两个检验从而确定d
#绘制初始数据图和pacf、acf,分布直方图
fig,ax = plt.subplots(2,2)
##查看原始序列
sd.plot(ax=ax[0][0],style='--')
#原始数据的频数分布直方图
sd.hist(ax=ax[0][1])
##使用差分使数据平稳
sd = sd.diff().diff().dropna()#差分后去除空值,这里二阶差分,所以d=2
## 自相关系数ACF 画图:X代表阶数, 阶数=2 代表y与y-2的自相关性
plot_acf(sd,lags=30,ax=ax[1][0])
##偏自相关系数剔除其他变量的影响下的ACF 比如y可能与y-1由关系,在计算y与y-2的自相关性系数时候,要剔除
plot_pacf(sd,lags=30,ax=ax[1][1])
# sns.despine(right=True,trim=True)
plt.tight_layout()
plt.savefig(r'D:\Python Project_1\Graph\acfpacf.svg')
"""
拖尾:相关系数指数衰减,始终由非零取值,几何型或震荡型
"""
"""confirm q and p by BIC and AIC"""
#继续,根据AIC和BIC准则确定pq参数, AIC和BIC越小越好
pmax = 4
qmax = 4
bic_matrix = np.zeros((pmax,qmax))
aic_matrix = np.zeros((pmax,qmax))
for p in range(pmax):
tmp = []
for q in range(qmax):
try:
model = sm.tsa.SARIMAX(sd.values,order=(p,1,q)).fit()
bic_matrix[p,q] = model.bic
aic_matrix[p, q] = model.aic
except:
bic_matrix[p,q] = None
aic_matrix[p,q] = None
bic = np.array(bic_matrix)
aic = np.array(aic_matrix)
fig,(ax1,ax2) = plt.subplots(1,2,figsize=(10,4))
np.savetxt(r'D:\Python Project_1\Result\bic.txt',X=bic,fmt='%.5e',delimiter=',')#保存的是BIC
sns.heatmap(bic,cmap='Reds',annot=True,mask=None,fmt='.2f',ax=ax1)
ax1.set_xlabel('p');ax1.set_ylabel('q');ax1.set_title('BIC');
sns.heatmap(aic,cmap='Reds',annot=True,mask=bic==None,fmt='.2f',ax=ax2)
plt.xlabel('p');plt.ylabel('q');ax2.set_title('AIC');plt.tight_layout();plt.show();#根据热力图选定参数qp
"""create a model to predict"""
print('='*40+'模型参数'+'='*40)
p = int(input("p = "));d = int(input("d = "));q = int(input("q = "));
Model = sm.tsa.SARIMAX(data,order=(p,d,q)).fit()
print('*'*50,'\n',Model.summary())#模型报告
"""check out ARIMA Model"""
print('='*40+'模型残差检验'+'='*40)
# 检验残差正太分布,说明残差不存在相关性,并且满足正态分布,说明模型预测良好
#方法一:QQ图
qqplot(data=Model.resid,line='q',fit=True);plt.show()
#方法二:KS检验
KS_res = stats.kstest(Model.resid, 'norm', (Model.resid.mean(), Model.resid.std()))
print(f'-->>[方法二]:K-S正太分布检验P-Value{KS_res[1]}')#p值大于0.05,为正态分布。
#方法三D-W检验
# D-W检验,是目前检验自相关最常用的方法,但是只适用于检验一阶自相关,因为一阶自相关系数pearson介于-1与1之间,所以0<=DE<=4。
# 并且DW=0 <--> 相关系数= 1
# DW=4 <--> 相关系数= -1
# DW=2 <--> 不存在(一阶自相关) 只有满足这个条件才可能
DW = sm.stats.durbin_watson(Model.resid)#接近2即满足要求
print((f"-->>[方法三]:D-W检验值:{DW}"))
# 方法四LB检验,当P-value都显著性大于0.05,则认为残差序列为白噪声序列,即模型拟合效果较好。
print(f"-->>[方法四]:白噪声检验:\n{acorr_ljungbox(Model.resid,lags=d)}")#看第一行第二个p值即可
"""Prediction """
print('='*40+'模型预测数据'+'='*40)
pre = Model.forecast(50)
print("-->>[预测数据]:\n",pre)
df = np.hstack((data,pre))
plt.plot(df)
plt.show()