第⼀次尝试使⽤Python创建季节性ARIMA模型
Python构建时间序列ARIMA模型
前⾔
时间：19/04/25，天⽓晴，⼼情愉悦。
⼤学⽣活不知不觉已经快要结束了，我才第⼀次认真的开始检查、审视⾃⼰过去两年中究竟学到了些什么？仔细回想，学的东西好像还蛮多的，从我开始懵懂学习C语⾔后⼊坑编程系列学习以来，⾃⼰慢慢的喜欢上这⼀领域，开始⽤C实现⼀些⼩⼩的功能，例如12306抢票、⼀些经典的⼩游戏等。后续除了掌握⼀些基本数据库知识之外我⼜接触了Java这种oo语⾔，并且我在⾃⼰的课程设计⾥⾯独⽴编写了基于MINA框架的商城Shopping⼩程序，虽然写的有点烂，但是也基本实现了商城的⼤部分功能（对⾃⼰要求有点低）。再后来⼀次侥幸的机会跟着导师做了⼀个项⽬，是基于计算机视觉的智能⼩车导航避障的，为此，我专门去学习了在MFC框架下⽤opencv实现图像处理，嗯，opencv图像处理是我学了最长时间也是最系统的学习，从简单的图像形态学处理、图像去燥、图像分割到复杂颇具挑战的模式识别我都有了解过研究过，在做图像处理的过程中，我慢慢的发现⾃⼰对图像数据（数据）的各式各样处理计算颇感兴趣。那时候正值⼤数据的热潮，我从⽹上了⼀些⼤数据的视频学习了Hadoop、Spark等，掌握了⼤批量数据下并⾏处理的⽅式，但是⾃⼰并不知道⽤什么⽅法去处理数据，于是，18年的下半年，
我开始正式接触机器学习领域，⾄今学习机器学习时常也有半年的时间了，在这半年⾥也学习了很多公式推导也有很多算法流程，如今让我回忆⼀下⾃⼰曾经学的这些机器学习的东西，⾃⼰的印象其实已经没有那么深刻了，即使在学习当下对算法推导以及应⽤⾮常的熟悉，但是没有把⾃⼰学的内容整理下来，即使学过随着时间推移也差不多把⼤部分细节给忘记了。我挺后悔⾃⼰当初学习的时候没有做笔记，不过悔在过去、⾏在当下，从此刻开始记录我的学习，我相信我可以坚持记录下去。
ARIMA模型简介
具有如下结构的模型称为求和⾃回归移动平均(Autoregressive Integrated Moving Average)，简记为ARIMA(p,d,q)模型：
式中：
ARIMA模型的数学原理和公式我就不⼀⼀解释了，如果有不懂的地⽅可以参考，在本篇中我们主要探讨如何⽤Python对时序进⾏分析建模。
Python实现ARIMA模型预测
数据的获取与准备
Python3：
import csv
data_files = ader(open('data.csv','r'))
dataSet =[] # 训练数据集
dateSet =[] # 训练年份集
for data in data_files:
if data[0]=='200301': # 将csv⽂件中验证数据去除
break
if data[0]!='date':
dateSet.append(data[0])
dataSet.append(float(data[1]))
绘制1995-2002年时间序列趋势图
显⽰⼀⾏⼀⾏的数据是难以观察数据的趋势⾛向，这样不利于我们观测数据之间的相互联系，因此画出趋势图有利于数据分析。
Python3：
# 利⽤Pandas库来绘制时间序列趋势图
import matplotlib.pyplot as plt
import pandas as pd
df = pd.DataFrame(dataSet,index=dateSet,columns=['real value'])
df.plot()
plt.show()
result：
去均值化后ADF平稳性检验以及差分
①为什么要去均值化？做数据的去均值化其实就是对数据进⾏标准化，在很多情况下，我们对数据本⾝⼤⼩并不感兴趣，我们感兴趣的是数据之间的联系，去均值化并不会影响这种联系，因为去均值化后得到的时间序列趋势和没有去均值化得到的时间序列趋势是⼀样的。Python3：
# 对时间序列进⾏去均值化
mean =sum(dataSet)/len(dataSet) # 计算均值
dataSet_kill_mean =[data - mean for data in dataSet] # 得到去均值后的序列
②序列平稳是做ARIMA回归的前提条件，所以我们得到的序列必须是平稳的。观测上⾯的序列趋势图，我们可以很明显的发现数据有逐渐上升的趋势，因此我们判断该段时间序列为⾮平稳序列。为了得到平稳的时间序列，我们要做差分，因为差分⽬的是使变量序列平稳。Python3：
# 进⾏⼀阶差分并绘制趋势图
df_kill_mean = pd.DataFrame(dataSet_kill_mean,index=dateSet,columns=['kill mean value'])
df_kill_mean_1 = dataSet_kill_mean.diff(1).dropna()
df_kill_mean_1.plot()
plt.show()
import statsmodels.tsa.stattools as ts
adf_summary = ts.adfuller(np.array(df_kill_mean_1).reshape(-1)) # 进⾏ADF检验并打印结果
print(adf_summary)
>>>(-3.3635866783352,0.012264631212932628,12,82,{'1%':-3.512738056978279,'5%':-2.8974898650628984,'10%':-2.585948732897085},237.8 874357065477)
如何观察上⾯的统计结果来判断序列是否为平稳呢？
1. 1%、5%、10%不同程度拒绝原假设的统计值和ADF Test result的⽐较，ADF Test result同时⼩于1
%、5%、10%即说明⾮常好
地拒绝该假设，本数据中，adf结果为-3.3635866783352， ⼤于1% level的-3.512738056978279，但⼩于5% level的-
2.8974898650628984，如果建模对序列平稳性有相当⾼的要求，则认为该序列不平稳，若要求并不严格也可认为该序列是平稳
序列。如果⼀阶差分序列不平稳，则继续做⼆阶差分并检验。
2. P-value是否⾮常接近0。本结果中，P-value 为 0.012264631212932628 < 0.05并接近0。
由此，通过检验观测值，我们认为上述⼀阶差分后得到的数据满⾜平稳性检验要求。（⼀阶差分后的序列平稳性不太好，有可能通不过⽩噪声检验，在这⾥忽略⽩噪声检验环节，若⽩噪声检验得到的P值⼤于0.05，那么我们就得对时间序列进⾏⼆阶差分）
绘制⾃相关函数以及偏相关函数图确定参数p、q
Python3：
aphics.tsaplots import plot_acf, plot_pacf
# 绘制⾃相关图
plot_acf(df_kill_mean_1,lags=24).show() # 其中lags参数是指横坐标最⼤取值
# 绘制偏相关图
plot_pacf(df_kill_mean_1,lags=24).show()
plt.show()
rusult：
我们观测ACF函数图和PACF函数图，我们发现并没有呈现很好的拖尾或截尾情况。出现这个情况的原因有很多种情况，我们观测函数图发现，每隔⼀个周期，ACF和PACF都出现“尖峰“，例如上述图每隔12个⽉出现⼀个"尖峰"，由此我们可以很容易的判断，该序列可能存在季节性影响的因素（其实由开始的序列趋势图我们也可以看出序列存在季节性影响）。
我们可以通过分解的⽅式将时序数据分离成不同的成分，它主要将时序数据分离为Trend（成长趋势）、seasonal（季节性趋势）、Residuals（随机成分）。然后我们分别对这三个分离的序列进⾏ARIMA建模得到较好的模型，最后再将模型相加便可以得到最后的ARIMA模型。
Python3：
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition =seasonal_decompose(df_kill_mean_1, freq=12)
trend = d  # 趋势部分
seasonal = decomposition.seasonal # 季节性部分
residual = sid # 残留部分
decomposition.plot()
result：
java调用python模型这种分解建模的⽅式有些复杂，接下来我们采⽤季节性时间序列来对上述具有明显季节性的序列进⾏建模。
建⽴季节性时间序列模型ARIMA(k,D,m)S×(p,d,q)
我们称ARIMA(k,D,m)S×(p,d,q)为乘积季节模型，也可以写成ARIMA(p,d,q)(k,D,m)S模型，其中S为季节性周期。如果将模型中的AR因⼦和MA因⼦分别展开，可以得到类似的ARMA(kS+p,mS+q)的模型。当我们考虑⽤ARIMA(p,d,q)(k,D,m)S模型的时候，我们需要优化感兴趣度量的是ARIMA(p,d,q)(k,D,m)s各个参数的值。接下来，我们将使⽤“⽹格搜索”来迭代地探索参数的不同组合。 对于参数的每个组合，我们使⽤statsmodels模块的SARIMAX()函数拟合⼀个新的季节性ARIMA模型，并评估其整体质量。
Python3：
# 这段代码借鉴了其他博⽂的做法
import itertools
p = q =range(0,2) # p、q⼀般取值不超过2
d =range(1,2)
pdq =list(itertools.product(p, d, q))
seasonal_pdq =[(x[0], x[1], x[2],12)for x in list(itertools.product(p, d, q))]
warnings.filterwarnings("ignore") # 忽略警告信息
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(df_mean_1,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print('ARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic))
except:
continue
>>>
ARIMA(0,1,0)x(0,1,0,12)-AIC:322.6198841487564
ARIMA(0,1,0)x(0,1,1,12)-AIC:270.20449054798723
ARIMA(0,1,0)x(1,1,0,12)-AIC:279.7073077170952
ARIMA(0,1,0)x(1,1,1,12)-AIC:272.20429383065317
ARIMA(0,1,1)x(0,1,0,12)-AIC:248.2004574618413
ARIMA(0,1,1)x(0,1,1,12)-AIC:209.46795258585362
ARIMA(0,1,1)x(1,1,0,12)-AIC:219.76010995507397
ARIMA(0,1,1)x(1,1,1,12)-AIC:213.63266407486356
ARIMA(1,1,0)x(0,1,0,12)-AIC:292.5090537220262
ARIMA(1,1,0)x(0,1,1,12)-AIC:250.76816771015618
ARIMA(1,1,0)x(1,1,0,12)-AIC:251.43830162843227
ARIMA(1,1,0)x(1,1,1,12)-AIC:251.1402407735711
ARIMA(1,1,1)x(0,1,0,12)-AIC:243.34941237032209
ARIMA(1,1,1)x(0,1,1,12)-AIC:206.44248559031112
ARIMA(1,1,1)x(1,1,0,12)-AIC:211.64269806160195
ARIMA(1,1,1)x(1,1,1,12)-AIC:210.80404367428187
在这⾥，我们是通过最佳准则函数法来确定具体的参数值的，模型选择的AIC准则越⼩，我们就⼤概可
以认为该模型越优。通过上述结果我们可以容易得到ARIMA(0, 1, 1)(0, 1, 1,)12模型相对最优。因此我们选择ARIMA(0, 1, 1)(0, 1, 1,)12模型作为我们预测时间序列的最佳模型。
模型预测
我们使⽤我们选取的最佳模型进⾏预测。
Python3：
import statsmodels.api as sm
df = pd.DataFrame(dataSet,index=dateSet,columns=['real value']) # 原始时间序列
mod = sm.tsa.statespace.SARIMAX(df,order=(0,1,1),seasonal_order=(1,1,1,12),enforce_stationarity=False,enforce_invertibility=False)
result = mod.fit()
predict_sunspots = result.forecast(12)
forcast = np.array(predict_sunspots[:]).reshape(-1)
print(forcast)
>>>
[19.4291208517.4061460621.5329018622.2789123.0691896423.69105707
22.5403614223.4921939724.5610886324.6996556226.023*******.41517033]
模型指标MAPE
计算预测值和真实值得MAPE指标，得到MAPE指标值为3.14%，通过这个指标我们可以认为该模型是⼀个好模型。
预测值和真实值趋势对⽐图
结尾
这是我第⼀次建⽴季节性ARIMA模型，如果中间出现错误或者有不合理的地⽅，还望各位海涵并指正出来，⼤家⼀起学习就是要相互纠错，只有这样⼤家才可以在学习上相互收益、相互进步。同时我也希望这篇⽂章可以给苦恼于如何⽤Python建⽴时间序列模型的朋友指明⼀条明路。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。