Question

本章容量预测模型的建模流程如图11-6所示。

图11-6　容量预测建模图

首先需要对观测值序列进行平稳性检验，如果不平稳，则对其进行差分处理直到差分后的数据平稳。在数据平稳后，对其进行白噪声检验。如果没有通过白噪声检验，就进行模型识别，识别其模型属于AR、MA和ARMA中的哪一种模型。并且通过BIC信息准则对模型进行定阶，确定ARIMA模型的p，q参数。在模型识别后需进行模型检验，检测模型残差序列是否为白噪声序列。如果模型没有通过检测，需要对其进行重新识别。对已通过检验的模型采用极大似然估计方法进行模型参数估计。最后，应用模型进行预测，将实际值与预测值进行误差分析。如果误差比较小（误差阈值需通过业务分析进行设定），表明模型拟合效果较好，则模型可以结束。反之需要重新估计参数。

模型构建的过程中需要用到以下方法。

1）平稳性检验：为了确定原始数据序列中没有随机趋势或确定趋势，需要对数据进行平稳性检验，否则将会产生“伪回归”的现象。本章采用单位根检验（ADF）的方法或者时序图的方法进行平稳性检验，其检验的结果见表11-6，时序图的方法见11.2.2小节数据探索分析。

表11-6　平稳性检验结果

平稳性检验的Python代码如代码清单11-2所示。

代码清单11-2　平稳性检验代码

#-*- coding: utf-8 -*-
#平稳性检验
import pandas as pd
#参数初始化
discfile = '../data/discdata_processed.xls'
data = pd.read_excel(discfile)
data = data.iloc[: len(data)-5] #不使用最后5个数据
#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
diff = 0
adf = ADF(data['CWXT_DB:184:D:\\'])
while adf[1] >= 0.05: # adf[1]为p值，p值小于0.05认为是平稳的。
  diff = diff + 1
  adf = ADF(data['CWXT_DB:184:D:\\'].diff(diff).dropna())
print(u'原始序列经过%s阶差分后归于平稳，p值为%s' %(diff, adf[1]))

代码详见：demo/code/stationarity_test.py

2）白噪声检验：为了验证序列中有用的信息是否已被提取完毕，需要对序列进行白噪声检验。如果序列检验为白噪声序列，就说明序列中有用的信息已经被提取完毕了，剩下的全是随机扰动，无法进行预测和使用。本章采用LB统计量的方法进行白噪声检验，其结果见表11-7。

表11-7　白噪声检验结果

白噪声检验的Python代码如代码清单11-3所示。

代码清单11-3　白噪声检验代码

#-*- coding: utf-8 -*-
#白噪声检验
import pandas as pd
#参数初始化
discfile = '../data/discdata_processed.xls'
data = pd.read_excel(discfile)
data = data.iloc[: len(data)-5] #不使用最后5个数据
#白噪声检测
from statsmodels.stats.diagnostic import acorr_ljungbox
[[lb], [p]] = acorr_ljungbox(data['CWXT_DB:184:D:\\'], lags = 1)
if p < 0.05:
  print(u'原始序列为非白噪声序列，对应的p值为：%s' %p)
else:
  print(u'原始该序列为白噪声序列，对应的p值为：%s' %p)
[[lb], [p]] = acorr_ljungbox(data['CWXT_DB:184:D:\\'].diff().dropna(), lags = 1)
if p < 0.05:
  print(u'一阶差分序列为非白噪声序列，对应的p值为：%s' %p)
else:
  print(u'一阶差分该序列为白噪声序列，对应的p值为：%s' %p)

代码详见：demo/code/whitenoise_test.py

3）模型识别：采用极大似然比方法进行模型的参数估计，估计各个参数的值。然后针对各个不同模型，采用BIC信息准则对模型进行定阶，确定p，q参数，从而选择最优模型。根据此方法选择的模型结果见表11-8。

表11-8　模型结果

模型识别代码见代码清单11-4。

代码清单11-4　模型识别代码

#-*- coding: utf-8 -*-
#确定最佳p、d、q值
import pandas as pd
#参数初始化
discfile = '../data/discdata_processed.xls'
data = pd.read_excel(discfile, index_col = 'COLLECTTIME')
data = data.iloc[: len(data)-5] #不使用最后5个数据
xdata = data['CWXT_DB:184:D:\\']
from statsmodels.tsa.arima_model import ARIMA
#定阶
pmax = int(len(xdata)/10) #一般阶数不超过length/10
qmax = int(len(xdata)/10) #一般阶数不超过length/10
bic_matrix = [] #bic矩阵
for p in range(pmax+1):
  tmp = []
  for q in range(qmax+1):
    try: #存在部分报错，所以用try来跳过报错
      tmp.append(ARIMA(xdata, (p,1,q)).fit().bic)
    except:
      tmp.append(None)
  bic_matrix.append(tmp)
bic_matrix = pd.DataFrame(bic_matrix) #从中可以找出最小值
p,q = bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置
print(u'BIC最小的p值和q值为：%s、%s' %(p,q))

代码详见：demo/code/find_optimal_pq.py

4）模型检验：模型确定后，检验其残差序列是否为白噪声。如果不是白噪声，说明残差中还存在有用的信息，需要修改模型或者进一步提取。本案例所确定的ARIMA（0，1，1）模型成功地通过了检验。模型检验代码如代码清单11-5所示。

代码清单11-5　模型检验代码

#-*- coding: utf-8 -*-
#模型检验
import pandas as pd
#参数初始化
discfile = '../data/discdata_processed.xls'
lagnum = 12 #残差延迟个数
data = pd.read_excel(discfile, index_col = 'COLLECTTIME')
data = data.iloc[: len(data)-5] #不使用最后5个数据
xdata = data['CWXT_DB:184:D:\\']
from statsmodels.tsa.arima_model import ARIMA #建立ARIMA(0,1,1)模型
arima = ARIMA(xdata, (0, 1, 1)).fit() #建立并训练模型
xdata_pred = arima.predict(typ = 'levels') #预测
pred_error = (xdata_pred - xdata).dropna() #计算残差
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
lb, p= acorr_ljungbox(pred_error, lags = lagnum)
h = (p < 0.05).sum() #p值小于0.05，认为是非白噪声。
if h > 0:
  print(u'模型ARIMA(0,1,1)不符合白噪声检验')
else:
  print(u'模型ARIMA(0,1,1)符合白噪声检验')

代码详见：demo/code/arima_model_check.py

5）模型预测：应用通过检验的模型进行预测，获取未来5天的预测值，并且与实际值作比较，也就是我们在建模型的时候所忽略的最后5个数据。为了方便比较，将单位换算成GB，其结果见表11-9。

表11-9　预测结果