周期性时间序列的预测
背景
最近在研究时间序列的时候,发现很多序列具有很强的周期性,那如何对此类序列进行预测呢?
数据处理
挑选一个如下图的具有周期性的时间序列。该序列是取得是过去7天的数据,每小时一个点,一共7*24个点。
划分数据集
我们取前六天的数据做训练,第七天做测试集。
平滑处理
时间序列经常会出现毛刺的点,需要做平滑处理才能分析,类似上图中的数据。消除数据的毛刺,可以用移动平均法,但是移动平均有时候处理完后并不能使数据平滑,我这里采用的方法很简单,但效果还不错:把每个点与上一点的变化值作为一个新的序列,对这里边的异常值,也就是变化比较离谱的值剃掉,用前后数据的均值填充:
def diff_smooth(ts, interval):
'''时间序列平滑处理'''
# 间隔为1小时
wide = interval/60
# 差分序列
dif = ts.diff().dropna()
# 描述性统计得到:min,25%,50%,75%,max值
td = dif.describe()
# 定义高点阈值,1.5倍四分位距之外
high = td['75%'] + 1.5 * (td['75%'] - td['25%'])
# 定义低点阈值
low = td['25%'] - 1.5 * (td['75%'] - td['25%'])
i = 0
forbid_index = dif[(dif > high) | (dif < low)].index
while i < len(forbid_index) - 1:
# 发现连续多少个点变化幅度过大
n = 1
# 异常点的起始索引
start = forbid_index[i]
while forbid_index[i+n] == start + datetime.timedelta(minutes=n):
n += 1
i += n - 1
# 异常点的结束索引
end = forbid_index[i]
# 用前后值的中间值均匀填充
value = np.linspace(ts[start - datetime.timedelta(minutes=wide)], ts[end + datetime.timedelta(minutes=wide)], n)
ts[start: end] = value
i += 1
return ts
经过处理以后,上图的时间序列得到了平滑处理,效果如下图。
周期性分解
具有周期性特征的序列需要将周期性特征提取出来。python里面的statsmodels工具包里面有针对周期性分解的函数seasonal_decompose,我们可以将序列进行分解。seasonal_decompose这个函数里面有个two_sided的参数,默认是True。Trend处理的时候用到移动平均的方法,熟悉此方法的读者就会发现,经过该方法处理以后,序列收尾两段有一部分数据缺失了,但是如果该参数为FALSE,则只有开始的时候有一段缺失值。
decomposition = seasonal_decompose(smooth_data, two_sided=False)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
图3中的第一张图是observed,体现的原始数据;第二张是trend,体现的是分解出来的趋势部分;第三张是seasonal,体现的是周期部分;最后是residual,体现的是残差部分。
本文采用的是seasonal_decompose的加法模型进行的分解,即 observed = trend + seasonal + residual,另还有乘法模型。在建模的时候,只针对trend部分学习和预测,如何将trend的预测结果加工成合理的最终结果?后面会有介绍。
预测
我们对trend部分进行预测,最后再加上seasonal部分。对trend的预测,我们采用ARIMA模型。熟悉该模型的都知道,需要确定三个参数p,q和d,可以使用aic和bic的方法进行定阶,可以查阅相关的文献。
trend.dropna(inplace=True)
# arima的训练参数order =(p,d,q),具体意义查看官方文档,调参过程略。
predict_model = ARIMA(trend, (1,1,2), freq='H').fit(disp=-1, method='css')
得到模型以后,就可以进行预测。
# 续接train,生成长度为n的时间索引,赋给预测序列
n = 24
predict_time_index = pd.date_range(start=trend.index[-1], periods=n+1, freq='H')[1:]
trend_predict = predict_model.forecast(n)[0]
'''
为预测出的趋势数据添加周期数据和残差数据
'''
values = []
low_conf_values = []
high_conf_values = []
# enumerate() 函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列,同时列出数据和数据下标,一般用在for循环当中。
for i, t in enumerate(predict_time_index):
trend_part = trend_predict[i]
# 相同时间点的周期数据均值
# t为2018-08-09 15:18:00类型的时间,t.time()为15:18:00类型的时间
season_part = seasonal[seasonal.index.time == t.time()].mean()
# 趋势 + 周期 + 误差界限
predict = trend_part + season_part
values.append(predict)
# 得到预测值
final_predict = pd.Series(values, index=predict_time_index, name='predict')
下面是预测的结果,从图中可以看到预测的结果将周期性的特征完美地体现出来了。
评估
对第七天作出预测,评估的指标为均方根误差rmse,本序列的rmse小于5,效果还是不错的。
总结
本文介绍了周期性序列的预测方法,你可能会问并不是所有的序列都具有周期性,事实确实如此,接下来几篇博客,我会重点介绍周期性检测的一些方法。希望此博客对您研究时间序列有所帮助。