|
|
# 代码9-6
|
|
|
#ARIMA模型
|
|
|
import pandas as pd
|
|
|
import matplotlib.pyplot as plt
|
|
|
import warnings
|
|
|
import math
|
|
|
import itertools
|
|
|
import numpy as np
|
|
|
import statsmodels.api as sm
|
|
|
from statsmodels.graphics.tsaplots import plot_acf
|
|
|
from statsmodels.stats.diagnostic import acorr_ljungbox
|
|
|
from statsmodels.tsa.stattools import adfuller as ADF
|
|
|
from statsmodels.tsa.arima_model import ARIMA
|
|
|
import matplotlib.ticker as ticker
|
|
|
|
|
|
|
|
|
on_h = pd.read_csv('../tmp/on_h.csv', index_col=0
|
|
|
,encoding='utf-8')
|
|
|
train = pd.DataFrame(on_h.iloc[0:426, 0])
|
|
|
test = pd.DataFrame(on_h.iloc[426:, 0])
|
|
|
#train.plot(title = '训练集时序图') # 画训练集时序图
|
|
|
x = train.index
|
|
|
y = train.on_man
|
|
|
fig, ax = plt.subplots(1,1)
|
|
|
ax.plot(x, y)
|
|
|
ticker_spacing = 70
|
|
|
ax.xaxis.set_major_locator(ticker.MultipleLocator(ticker_spacing))
|
|
|
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
|
|
|
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('客流量')
|
|
|
plt.title('训练集时序图')
|
|
|
plt.tight_layout() # tight_layout()方法可以保证图像的完整度
|
|
|
plt.show()
|
|
|
|
|
|
plot_acf(train,lags=400)
|
|
|
plt.xlabel('日期索引')
|
|
|
plt.ylabel('自相关系数')
|
|
|
plt.title('训练集自相关图')
|
|
|
plt.tight_layout()
|
|
|
plt.show() # 画训练集自相关图
|
|
|
|
|
|
print('原始序列的ADF检验结果为：', ADF(train['on_man'])) # 检验训练集平稳性
|
|
|
print('白噪声检验结果为：', acorr_ljungbox(train['on_man'], lags=1))
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-7
|
|
|
import statsmodels.api as sm
|
|
|
train['on_man'] = train['on_man'].astype(float)
|
|
|
# 定阶
|
|
|
bic_matrix = [] # bic矩阵，选择bic矩阵中最小值对应的行(p),列(q)
|
|
|
# 存在部分报错，所以用try来跳过报错
|
|
|
for p in range(11):
|
|
|
tmp = []
|
|
|
for q in range(5):
|
|
|
try:
|
|
|
tmp.append(sm.tsa.ARIMA(train, order=(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('BIC最小的p值和q值为：%s、%s' % (p, q))
|
|
|
|
|
|
model = sm.tsa.ARIMA(train, order=(p, 1, q)).fit() # 建立ARIMA(p, 1, 1)模型
|
|
|
summary = model.summary() # 给出一份模型报告
|
|
|
forecast = model.forecast(10)
|
|
|
print('10天的预测结果、标准误差和置信区间分别为：\n', forecast)
|
|
|
|
|
|
# pre = pd.DataFrame(forecast, columns = ['predict'])
|
|
|
pre = pd.DataFrame({'predict':forecast})
|
|
|
pre.index = test.index #如果索引不同将pre加到test中时会出错
|
|
|
test['pre'] = pre
|
|
|
plt.plot(test.index, test.on_man)
|
|
|
plt.plot(test.index, test.pre, linestyle=':')
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.legend(['真实', '预测'])
|
|
|
plt.title('预测结果和实际结果的对比')
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('上车人数')
|
|
|
plt.tight_layout() # tight_layout()方法可以保证图像的完整度
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-8
|
|
|
# 剔除节假日之后的时序图
|
|
|
on_nh = on_h[on_h.iloc[:,2] != '小长假'] # 剔除假期
|
|
|
x = on_nh.index
|
|
|
y = on_nh.on_man
|
|
|
fig, ax = plt.subplots(1,1)
|
|
|
ax.plot(x, y)
|
|
|
plt.xticks(range(len(x)), on_nh['date.1'])
|
|
|
ticker_spacing = 70
|
|
|
ax.xaxis.set_major_locator(ticker.MultipleLocator(ticker_spacing))
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.title('剔除节假日后ST111-01站点客流量')
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('客流量（人）')
|
|
|
plt.tight_layout() # tight_layout()方法可以保证图像的完整度
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-9
|
|
|
# 剔除节假日及前一天后的时序图
|
|
|
new_h = on_h #存放节假日更新的数据
|
|
|
new_h.index = range(len(new_h))
|
|
|
for i in range(1, len(new_h)):
|
|
|
if new_h.iloc[i, 2] == '小长假' :
|
|
|
new_h.iloc[i-1, 2] = '小长假'
|
|
|
new_nh = new_h[new_h.iloc[:,2] != '小长假'] # 剔除节假日及其前一天
|
|
|
|
|
|
plt.plot(new_nh['date.1'], new_nh['on_man'],color='green')
|
|
|
plt.gca().xaxis.set_major_locator(ticker.MultipleLocator(70))
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.title('剔除节假日及其前一天后ST111-01站点客流量')
|
|
|
plt.legend(['on_man'], loc=2)
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('客流量（人）')
|
|
|
plt.tight_layout() # tight_layout()方法可以保证图像的完整度
|
|
|
plt.savefig('../tmp/9-10.png', dpi=1080)
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
# 代码9-10
|
|
|
new_nh.index = range(len(new_nh))
|
|
|
train1 = pd.DataFrame(new_nh.iloc[0:378, 0])
|
|
|
test1 = pd.DataFrame(new_nh.iloc[378:, 0])
|
|
|
p = q = range(4)
|
|
|
d = range(2)
|
|
|
pdq = list(itertools.product(p, d, q))
|
|
|
seasonal_pdq = [(x[0], x[1], x[2], 7)for x in list(itertools.product(p, d, q))]
|
|
|
print('Examples of parameter combinations for Seasonal ARIMA...')
|
|
|
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
|
|
|
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
|
|
|
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
|
|
|
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
|
|
|
warnings.filterwarnings("ignore") # specify to ignore warning messages
|
|
|
sa =[]
|
|
|
for param in pdq:
|
|
|
for param_seasonal in seasonal_pdq:
|
|
|
try:
|
|
|
mod = sm.tsa.statespace.SARIMAX(train1,
|
|
|
order=param,
|
|
|
seasonal_order=param_seasonal,
|
|
|
enforce_stationarity=False,
|
|
|
enforce_invertibility=False)
|
|
|
results = mod.fit()
|
|
|
print('ARIMA{}x{}7 - AIC:{}'.format(param, param_seasonal, results.aic))
|
|
|
sa.append(param)
|
|
|
sa.append(param_seasonal)
|
|
|
sa.append(results.aic)
|
|
|
except:
|
|
|
continue
|
|
|
AIC = [i for i in sa if type(i) == np.float64]
|
|
|
AIC_min = min(AIC)
|
|
|
for i in np.arange(2,len(sa),3):
|
|
|
if sa[i] == min(AIC):
|
|
|
param = sa[i-2]
|
|
|
param_seasonal = sa[i-1]
|
|
|
mod = sm.tsa.statespace.SARIMAX(train1,
|
|
|
order=(param),
|
|
|
seasonal_order=(param_seasonal),
|
|
|
enforce_stationarity=False,
|
|
|
enforce_invertibility=False)
|
|
|
print('模型最终定阶为：', (param, param_seasonal))
|
|
|
results = mod.fit()
|
|
|
print(results.summary().tables[1])
|
|
|
fig = plt.figure(figsize=(15, 12))
|
|
|
results.plot_diagnostics(figsize=(15, 12), fig=fig)
|
|
|
plt.show()
|
|
|
pre_10 = results.predict(start=378, end=387,dynamic=True)
|
|
|
out_pre = pd.DataFrame(np.zeros([10,3]),columns = ['real', 'pre', 'error'])
|
|
|
out_pre['real'] = list(test1['on_man'])
|
|
|
out_pre['pre'] = list(pre_10)
|
|
|
# 计算相对误差
|
|
|
error_seasonal = (out_pre.loc[:, 'pre']-out_pre.loc[:,'real'])/out_pre.loc[:,'real']
|
|
|
# 平均相对误差
|
|
|
error_mean = abs(error_seasonal).mean()
|
|
|
print('预测平均相对误差为：', error_mean)
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-11
|
|
|
# 节假日客流规律
|
|
|
# 由于之前步骤对on_h有修改，所以此处重新载入
|
|
|
holiday = pd.read_csv('../tmp/holiday.csv', index_col=0
|
|
|
,encoding = 'utf-8')
|
|
|
Train_Station = pd.read_csv('../tmp/Train_Station.csv', index_col=0
|
|
|
,encoding='utf-8')
|
|
|
Train_ST111_01 = Train_Station[Train_Station.iloc[:, 0] == 'ST111-01']
|
|
|
on_h = Train_ST111_01.groupby('date')['on_man'].sum()
|
|
|
on_h = pd.DataFrame(on_h)
|
|
|
on_h['date'] = 0
|
|
|
on_h['holiday'] = 0
|
|
|
# 添加日期和类型（工作日或者小长假）
|
|
|
for i in range(len(holiday)):
|
|
|
for j in range(len(on_h)):
|
|
|
if holiday.iloc[i,0] == on_h.index[j]:
|
|
|
on_h.loc[on_h.index[j], 'holiday'] = holiday.iloc[i,1]
|
|
|
on_h.loc[on_h.index[j], 'date'] = holiday.iloc[i,0]
|
|
|
# 2015春节
|
|
|
fig = plt.figure(figsize=(12, 6)) # 设置画布
|
|
|
ax = fig.add_subplot(1, 1, 1)
|
|
|
ax.plot(on_h.loc['2015-01-19':'2015-02-18', 'on_man'], color = 'blue')
|
|
|
ax.plot(on_h.loc['2015-02-18':'2015-02-25', 'on_man'], color = 'red', linestyle=':')
|
|
|
ax.plot(on_h.loc['2015-02-25':'2015-03-01', 'on_man'], color = 'blue')
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('上车人数')
|
|
|
plt.title('2015春节客流量')
|
|
|
plt.legend(['工作日','节假日'])
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.show()
|
|
|
|
|
|
# 2015劳动节
|
|
|
fig1 = plt.figure(figsize=(12, 6)) # 设置画布
|
|
|
ax1 = fig1.add_subplot(1, 1, 1)
|
|
|
ax1.plot(on_h.loc['2015-04-27':'2015-05-01', 'on_man'], color = 'blue')
|
|
|
ax1.plot(on_h.loc['2015-05-01':'2015-05-04', 'on_man'], color = 'red', linestyle=':')
|
|
|
ax1.plot(on_h.loc['2015-05-04':'2015-05-11', 'on_man'], color = 'blue')
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('上车人数')
|
|
|
plt.title('五一客流量')
|
|
|
plt.legend(['工作日','节假日'])
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.show()
|
|
|
|
|
|
# 国庆和中秋
|
|
|
fig2 = plt.figure(figsize=(12, 6)) # 设置画布
|
|
|
ax2 = fig2.add_subplot(1, 1, 1)
|
|
|
ax2.plot(on_h.loc['2015-09-21':'2015-09-26', 'on_man'],color = 'blue')
|
|
|
ax2.plot(on_h.loc['2015-09-26':'2015-09-28', 'on_man'],color = 'red', linestyle=':')
|
|
|
ax2.plot(on_h.loc['2015-09-28':'2015-09-30', 'on_man'],color = 'blue')
|
|
|
ax2.plot(on_h.loc['2015-09-30':'2015-10-08', 'on_man'],color = 'red', linestyle=':')
|
|
|
ax2.plot(on_h.loc['2015-10-08':'2015-10-12', 'on_man'],color = 'blue')
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('上车人数')
|
|
|
plt.title('中秋和国庆客流量')
|
|
|
plt.legend(['工作日', '节假日'])
|
|
|
plt.xticks(rotation = 45)
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-12
|
|
|
# 2015和2016春节客流量比较
|
|
|
compare = pd.DataFrame(on_h.loc['2015-02-05':'2015-02-26', 'on_man'])
|
|
|
compare['2016'] = list(on_h.loc['2016-01-25':'2016-02-15', 'on_man'])
|
|
|
compare.columns = ['2015', '2016']
|
|
|
compare.index = range(len(compare))
|
|
|
plt.plot(compare.index, compare['2015'], linestyle=':')
|
|
|
plt.plot(compare.index, compare['2016'])
|
|
|
plt.legend(['2015', '2016'])
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('客流量')
|
|
|
plt.title('2015和2016春节客流量比较')
|
|
|
|
|
|
|
|
|
|
|
|
# 代码9-13
|
|
|
# 2015春节节假日波动系数
|
|
|
M = on_h.loc['2015-01-01':'2015-02-17']
|
|
|
M = M[M.loc[:,'holiday'] != '小长假']
|
|
|
M1 = on_h.loc['2015-01-01':'2015-02-07']
|
|
|
M1 = M1[M1.loc[:,'holiday'] != '小长假']
|
|
|
B_coef = []
|
|
|
# 春节前10天
|
|
|
for i in on_h.loc['2015-02-08':'2015-02-17',:].index:
|
|
|
B_coef.append('%.2f' % (on_h.loc[i, 'on_man']/math.ceil(M1.iloc[-30:].mean())))
|
|
|
# 春节及春节后两天
|
|
|
for i in on_h.loc['2015-02-18':'2015-02-26',:].index:
|
|
|
B_coef.append('%.2f' % (on_h.loc[i, 'on_man']/math.ceil(M.iloc[-30:].mean())))
|
|
|
B_coef = pd.DataFrame(B_coef)
|
|
|
B_coef.columns = ['on_man']
|
|
|
B_coef[u'on_man'] = B_coef[u'on_man'].astype(float)
|
|
|
fig3 = plt.figure(figsize=(8,6)) # 设置画布
|
|
|
ax3 = fig3.add_subplot(1, 1, 1)
|
|
|
ax3.plot(B_coef.iloc[:,0], color='blue')
|
|
|
plt.xlabel('Index')
|
|
|
plt.ylabel('系数')
|
|
|
plt.title('2015春节客流量波动系数')
|
|
|
# 设置数字标签
|
|
|
for a, b in zip(B_coef.index, B_coef.iloc[:, 0]):
|
|
|
plt.text(a, b, b, ha='center', va='bottom', fontsize=20)
|
|
|
plt.legend()
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
# 预测2016年春节客流量
|
|
|
MM = on_h.loc['2015-01-01':'2016-02-06',:]
|
|
|
MM = MM[MM.loc[:, 'holiday'] != '小长假']
|
|
|
MM_mean = math.ceil(MM.iloc[-30:].mean())
|
|
|
MM1 = on_h.loc['2015-01-01':'2016-01-27',:]
|
|
|
MM1 = MM1[MM1.loc[:, 'holiday'] != '小长假']
|
|
|
MM1_mean = math.ceil(MM1.iloc[-30:].mean())
|
|
|
pre_2016_b = B_coef.iloc[0:10, 0] * MM1_mean
|
|
|
pre_2016_a = B_coef.iloc[10:, 0] * MM_mean
|
|
|
pre_2016 = pd.DataFrame(on_h.loc['2016-01-28':'2016-02-15', 'on_man'])
|
|
|
pre_2016['pre'] = 0
|
|
|
pre_2016.loc[0:10,'pre'] = list(pre_2016_b)
|
|
|
pre_2016.loc[10:,'pre'] = list(pre_2016_a)
|
|
|
pre_2016.columns=['real', 'pre']
|
|
|
plt.plot(pre_2016.index, pre_2016.real)
|
|
|
plt.plot(pre_2016.index, pre_2016.pre, linestyle=':')
|
|
|
plt.xticks(rotation=45)
|
|
|
plt.legend(['real', 'pre'])
|
|
|
plt.xlabel('日期')
|
|
|
plt.ylabel('客流量')
|
|
|
plt.title('预测2016年春节客流量')
|
|
|
|
|
|
# 计算相对误差
|
|
|
error_pre = (pre_2016.loc[:, 'pre'] - pre_2016.loc[:, 'real'])/pre_2016.loc[:, 'real']
|
|
|
# 平均相对误差
|
|
|
error_pre_mean = abs(error_pre).mean()
|
|
|
print('预测的平均相对误差为：', error_pre_mean)
|
