首先要讲的是，在短时交通流预测领域，ARIMA模型是一种相当粗糙的基于时间序列的预测方法，相应地，这种预测方法的准确性也不高。同时，ARIMA方法也不在我的研究范围内，所以本次实验就直接使用网络上现成的python第三方库进行实现，不再重新造轮子。

ARIMA模型表达为ARIMA(p,d,q)，是在ARMA模型的基础上新引入了'I'也就是差分的思想。ARIMA模型由三个部分组成，分别是AR(p)模型、I(d)差分方法、MA(q)模型。

AR(p)模型(AutoRegressive)：p即为时间序列的长度，AR(p)模型的预测值是将过去p长度的时间序列的各真实值进行加权平均得到的。

I(d)差分：d为差分的阶数。差分的目的是将原本不平稳的时间序列经过差分，转化为更加平稳的序列。d的取值取决于原序列经过几次差分后可以转化为平稳或者弱平稳的序列。

MA(q)模型(MovingAverage):MA 模型适用于平稳的序列上，q也表示了过去的供预测的时间序列的长度，而通过MA模型进行的预测是将过去q个时间段的误差值取不同权重进行累加，再同平稳序列的均值加和得到的。

基于上述的思想，ARIMA(p,d,q)模型先经过d阶差分将原序列转化为较为平稳的时间序列，再通过计算确定效果最好的p、q的值，最后通过p、d、q的取值建立一个ARIMA模型。

实现代码如下：

import pandas as pd

from statsmodels.stats.diagnostic import acorr_ljungbox

from statsmodels.tsa.arima_model import ARIMA

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

import matplotlib.pyplot as plt

from statsmodels.tsa.stattools import adfuller

##通过数据的前288行计算p、d、q取值

filename = r'train13720.csv'

data = pd.read_csv(filename, index_col = 'timestamp').iloc[:288,:]

#画出折线图

#data.plot()

#plt.show()

#画出自相关性图

#plot_acf(data)

#plt.show()

#原始序列的平稳性检测

print('原始序列的检验结果为：',adfuller(data['hourly_traffic_count']))

#返回值依次为：adf, pvalue p值， usedlag, nobs, critical values临界值 , icbest, regresults, resstore

#单位检测统计量对应的p 值显著大于 0.05 ， 说明序列可以判定为 非平稳序列

#对数据进行差分后得到 自相关图和 偏相关图

D_data = data.diff().dropna()

D_data.columns = [u'流量差分']

#画出差分后的时序图

#D_data.plot()

#plt.show()

#画出自相关图

#plot_acf(D_data)

#plt.show()

#画出偏相关图

#plot_pacf(D_data)

#plt.show()

print(u'差分序列的ADF 检验结果为： ', adfuller(D_data[u'流量差分'])) #平稳性检验

#一阶差分后的序列的时序图在均值附近比较平稳的波动， 自相关性有很强的短期相关性， 单位根检验 p值小于 0.05 ，所以说一阶差分后的序列是平稳序列#

print(u'差分序列的白噪声检验结果：',acorr_ljungbox(D_data, lags= 1)) #返回统计量和 p 值

#对模型进行定阶

print(D_data)

#一般阶数不超过 length /10，但由于计算性能限制，故暂取leneth/30

pmax = int(len(D_data) / 30)

qmax = int(len(D_data) / 30)

bic_matrix = []

#计算不同p、q值对应的准确性收益

#print(bic_matrix)

for i in range(pmax+1):

temp=[]

for j in range(qmax+1):

try:

a=ARIMA(data,(i,1,j)).fit().bic

temp.append(a)

except:

temp.append(None)

finally:

bic_matrix.append(temp)

print(bic_matrix)

bic_matrix = pd.DataFrame(bic_matrix) #将其转换成Dataframe 数据结构

p,q = bic_matrix.stack().idxmin() #使用 idxmin 找出最小值的位置

print(u'BIC 最小的p值 和 q 值：%s,%s' %(p,q)) # BIC 最小的p值 和 q 值：0,1

#建立ARIMA 模型，ARIMA(0,1,1)

model = ARIMA(data, (p,1,q)).fit()

# model.summary2() #生成一份模型报告

# model.forecast(1) #进行1个时间单位的预测， 返回预测结果， 标准误差， 和置信区间

print('predict data:')

print((model.forecast(1)[0][0]))

通过上述的代码，我们就建立了一个关于该样本的ARIMA模型，其模型的p、d、q参数分别为0、1、1。通过下述的代码，使用该模型进行100步的预测：

p=0

d=1

q=1

result1=[]

for i in range(100):

temp=[]

data = pd.read_csv(filename, index_col = 'timestamp').iloc[i:288+i,:]

model = ARIMA(data, (p,d,q)).fit()

a = pd.read_csv(filename, index_col = 'timestamp').iloc[288+i:289+i,:].values[-1][0]

b = model.forecast(1)[0][0]

temp.append(a)

temp.append(b)

result1.append(temp)

result1=pd.DataFrame(result1)

print(result1)

得到实际观察值和预测值如下，左侧为观察值，右侧为预测值：

0 1

0 9 8.923906

1 12 8.957891

2 9 9.509028

3 6 9.455487

4 7 8.869332

.. .. ...

95 10 8.668516

96 5 8.893771

97 10 8.183186

98 15 8.506519

99 9 9.691928

[100 rows x 2 columns]

在论文《基于 ARIMA 模型的短时交通流实时自适应预测》(韩超，宋苏 ，王成红)中，作者根据他们的数据集求得了ARIMA(2,1,2)的模型，将其简单挪用过来如下：

p=2

d=1

q=2

result2=[]

for i in range(100):

temp=[]

data = pd.read_csv(filename, index_col = 'timestamp').iloc[i:288+i,:]

model = ARIMA(data, (p,d,q)).fit()

a = pd.read_csv(filename, index_col = 'timestamp').iloc[288+i:289+i,:].values[-1][0]

b = model.forecast(1)[0][0]

temp.append(a)

temp.append(b)

result2.append(temp)

result2=pd.DataFrame(result2)

print(result2)

得到实际观察值和预测值如下，左侧为观察值，右侧为预测值：

0 1

0 9 8.787160

1 12 8.915429

2 9 9.263472

3 6 9.627119

4 7 9.185042

.. .. ...

95 10 9.379399

96 5 9.061870

97 10 8.501785

98 15 8.161913

99 9 9.204021

[100 rows x 2 columns]

还可以使用刚才提到的文章定义的评价指标进行结果分析，代码如下：

def RE(r,p):

if r>p:

return (r-p)/r

else:

return (p-r)/r

rt1=[]

for i in result1.values:

temp=[]

for j in i:

temp.append(j)

temp.append(RE(i[0],i[1]))

rt1.append(temp)

c=0

sum=0

for i in rt1:

if i[-1]!=float('inf'):

s=float(i[-1])

sum+=s

c+=1

print(c)

print(sum/c)

rt1=pd.DataFrame(rt1)

print(rt1)

rt2=[]

for i in result2.values:

temp=[]

for j in i:

temp.append(j)

temp.append(RE(i[0],i[1]))

rt2.append(temp)

c=0

sum=0

for i in rt2:

if i[-1]!=float('inf'):

s=float(i[-1])

sum+=s

c+=1

print(c)

print(sum/c)

rt2=pd.DataFrame(rt1)

print(rt2)

得到两个模型的预测评价如下：

95

0.48252386991776014

0 1 2

0 9.0 8.923906 0.008455

1 12.0 8.957891 0.253509

2 9.0 9.509028 0.056559

3 6.0 9.455487 0.575914

4 7.0 8.869332 0.267047

.. ... ... ...

95 10.0 8.668516 0.133148

96 5.0 8.893771 0.778754

97 10.0 8.183186 0.181681

98 15.0 8.506519 0.432899

99 9.0 9.691928 0.076881

[100 rows x 3 columns]

95

0.46626540909769215

0 1 2

0 9.0 8.923906 0.008455

1 12.0 8.957891 0.253509

2 9.0 9.509028 0.056559

3 6.0 9.455487 0.575914

4 7.0 8.869332 0.267047

.. ... ... ...

95 10.0 8.668516 0.133148

96 5.0 8.893771 0.778754

97 10.0 8.183186 0.181681

98 15.0 8.506519 0.432899

99 9.0 9.691928 0.076881

[100 rows x 3 columns]

可以看到，即便是针对另外一个数据集，上述文献构建的预测模型取得了约0.47的MAPE值，而自行计算得到的模型取得了约0.48的MAPE值，这说明文献中构建的模型的预测结果总体上还是右于本次自行构建的模型的。