β

周期性时间序列的预测

奇虎360-addops 46 阅读
 

背景

最近在研究时间序列的时候,发现很多序列具有很强的周期性,那如何对此类序列进行预测呢?

数据处理

挑选一个如下图的具有周期性的时间序列。该序列是取得是过去7天的数据,每小时一个点,一共7*24个点。

alt

划分数据集

我们取前六天的数据做训练,第七天做测试集。

平滑处理

时间序列经常会出现毛刺的点,需要做平滑处理才能分析,类似上图中的数据。消除数据的毛刺,可以用移动平均法,但是移动平均有时候处理完后并不能使数据平滑,我这里采用的方法很简单,但效果还不错:把每个点与上一点的变化值作为一个新的序列,对这里边的异常值,也就是变化比较离谱的值剃掉,用前后数据的均值填充:

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

经过处理以后,上图的时间序列得到了平滑处理,效果如下图。

alt

周期性分解

具有周期性特征的序列需要将周期性特征提取出来。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

alt

图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')

下面是预测的结果,从图中可以看到预测的结果将周期性的特征完美地体现出来了。

alt

评估

对第七天作出预测,评估的指标为均方根误差rmse,本序列的rmse小于5,效果还是不错的。

总结

本文介绍了周期性序列的预测方法,你可能会问并不是所有的序列都具有周期性,事实确实如此,接下来几篇博客,我会重点介绍周期性检测的一些方法。希望此博客对您研究时间序列有所帮助。

 
作者:奇虎360-addops
应用运维|运维开发|opsdev|addops|虚拟化|openstack|docker|容器化|k8s|智能运维
原文地址:周期性时间序列的预测, 感谢原作者分享。

发表评论