精华内容
下载资源
问答
  • 周期性时间序列的预测

    万次阅读 2018-08-17 18:00:00
    最近在研究时间序列的时候,发现很多序列具有很强的周期性,那如何对此类序列进行预测呢? 1 数据处理 挑选一个如下图的具有周期性时间序列。该序列是取得是过去7天的数据,每小时一个点,一共7*24个点。   2 ...

    女主宣言

    AIOps 从立项到现在已经半年有余,从最开始的 LVS 异常检测,到如今的实时关联分析,智能运维已经渗透到我们日常运维中的许多场景,之后我们会将积累的经验分享出来,供大家学习参考,本文最先发布于 OpsDev,转载已获取作者授权。

    PS:丰富的一线技术、多元化的表现形式,尽在“HULK一线技术杂谈”,点关注哦!

    640?wx_fmt=jpeg

    At Tranquility Base, 1969

    by NASA IOTD

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


    1

    数据处理

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

    640?wx_fmt=png  

    2

    划分数据集

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


    3

    平滑处理

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

    640?wx_fmt=jpeg

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

    640?wx_fmt=png  

    4

    周期性分解

    具有周期性特征的序列需要将周期性特征提取出来。python里面的statsmodels工具包里面有针对周期性分解的函数seasonal_decompose,我们可以将序列进行分解。seasonal_decompose这个函数里面有个two_sided的参数,默认是True。Trend处理的时候用到移动平均的方法,熟悉此方法的读者就会发现,经过该方法处理以后,序列收尾两段有一部分数据缺失了,但是如果该参数为FALSE,则只有开始的时候有一段缺失值。

    640?wx_fmt=jpeg
    640?wx_fmt=png

    图3中的第一张图是observed,体现的原始数据;第二张是trend,体现的是分解出来的趋势部分;第三张是seasonal,体现的是周期部分;最后是residual,体现的是残差部分。


    本文采用的是seasonal_decompose的加法模型进行的分解,即 observed = trend + seasonal + residual,另还有乘法模型。在建模的时候,只针对trend部分学习和预测,如何将trend的预测结果加工成合理的最终结果?后面会有介绍。


    5

    预测

    我们对trend部分进行预测,最后再加上seasonal部分。对trend的预测,我们采用ARIMA模型。熟悉该模型的都知道,需要确定三个参数p,q和d,可以使用aic和bic的方法进行定阶,可以查阅相关的文献。

    640?wx_fmt=jpeg


    得到模型以后,就可以进行预测。

    640?wx_fmt=jpeg


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

    640?wx_fmt=png
     

    6

    评估

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


    7

    总结

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

    HULK一线技术杂谈

    由360云平台团队打造的技术分享公众号,内容涉及云计算数据库大数据监控泛前端自动化测试等众多技术领域,通过夯实的技术积累和丰富的一线实战经验,为你带来最有料的技术分享

    640?wx_fmt=gif
    展开全文
  • P4J是一个python软件包,用于基于信息论目标函数对不规则采样和异方差时间序列进行周期检测。 P4J是为天文光曲线,恒星大小或通量的不规则采样时间序列开发的。 此程序包的核心是一个称为“周期图”的类,该类扫描...
  • 有的开始比较平缓,后面开始增长… … 但是都属于典型的周期性时间序列,它的核心思想很简单:做好分解,做好预测结果的还原,和置信区间的设置,具体操作可根据具体业务逻辑做调整,祝大家建模愉快。 原文链接: ...

    公司平台上有不同的api,供内部或外部调用,这些api承担着不同的功能,如查询账号、发版、抢红包等等。

    日志会记录下每分钟某api被访问了多少次,即一个api每天会有1440条记录(1440分钟),将每天的数据连起来观察,有点类似于股票走势的意思。

    我想通过前N天的历史数据预测出第N+1天的流量访问情况,预测值即作为合理参考,供新一天与真实值做实时对比。当真实流量跟预测值有较大出入,则认为有异常访问,触发报警。

    数据探索

    我放了一份样例数据在data文件夹下,看一下数据大小和结构。

    data = pd.read_csv(filename)
    print('size: ',data.shape)
    print(data.head())
    

    file

    数据大小:

    共10080条记录,即10080分钟,七天的数据。

    字段含义:

    date:时间,单位分钟

    count:该分钟该api被访问的次数

    画图看一下序列的走势:(一些画图等探索类的方法放在了test_stationarity.py 文件中,包含时间序列图,移动平均图,有兴趣的可以自己尝试下)。

    def draw_ts(timeseries):
        timeseries.plot()
        plt.show()
    ​
    data = pd.read_csv(path)
    data = data.set_index('date')
    data.index = pd.to_datetime(data.index)
    ts = data['count']
    draw_ts(ts)
    

    file

    序列

    看这糟心的图,那些骤降为0的点这就是我遇到的第一个坑,我当初一拿到这份数据就开始做了。后来折腾了好久才发现,那些骤降为0的点是由于数据缺失,ETL的同学自动补零造成的,沟通晚了(TДT)。

    把坑填上,用前后值的均值把缺失值补上,再看一眼:

    file

    填充好缺失值的序列

    发现这份数据有这样几个特点,在模型设计和数据预处理的时候要考虑到:

    1、这是一个周期性的时间序列,数值有规律的以天为周期上下波动,图中这个api,在每天下午和晚上访问较为活跃,在早上和凌晨较为稀少。在建模之前需要做分解。

    2、我的第二个坑:数据本身并不平滑,骤突骤降较多,而这样是不利于预测的,毕竟模型需要学习好正常的序列才能对未知数据给出客观判断,否则会出现频繁的误报,令气氛变得十分尴尬( ´Д`),所以必须进行平滑处理。

    3、这只是一个api的序列图,而不同的api的形态差距是很大的,毕竟承担的功能不同,如何使模型适应不同形态的api也是需要考虑的问题。

    预处理

    1、划分训练测试集

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

    class ModelDecomp(object):
        def __init__(self, file, test_size=1440):
            self.ts = self.read_data(file)
            self.test_size = test_size
            self.train_size = len(self.ts) - self.test_size
            self.train = self.ts[:len(self.ts)-test_size]
            self.test = self.ts[-test_size:]
    

    2、对训练数据进行平滑处理

    消除数据的毛刺,可以用移动平均法,我这里没有采用,因为我试过发现对于我的数据来说,移动平均处理完后并不能使数据平滑,我这里采用的方法很简单,但效果还不错:

    把每个点与上一点的变化值作为一个新的序列,对这里边的异常值,也就是变化比较离谱的值剃掉,用前后数据的均值填充,注意可能会连续出现变化较大的点。

    def _diff_smooth(self, ts):
        dif = ts.diff().dropna() # 差分序列
        td = dif.describe() # 描述性统计得到:min,25%,50%,75%,max值
        high = td['75%'] + 1.5 * (td['75%'] - td['25%']) # 定义高点阈值,1.5倍四分位距之外
        low = td['25%'] - 1.5 * (td['75%'] - td['25%']) # 定义低点阈值,同上
    ​
        # 变化幅度超过阈值的点的索引
        forbid_index = dif[(dif > high) | (dif < low)].index 
        i = 0
        while i < len(forbid_index) - 1:
            n = 1 # 发现连续多少个点变化幅度过大,大部分只有单个点
            start = forbid_index[i] # 异常点的起始索引
            while forbid_index[i+n] == start + timedelta(minutes=n):
                n += 1
            i += n - 1
    ​
            end = forbid_index[i] # 异常点的结束索引
            # 用前后值的中间值均匀填充
            value = np.linspace(ts[start - timedelta(minutes=1)], ts[end + timedelta(minutes=1)], n)
            ts[start: end] = value
            i += 1
    ​
    self.train = self._diff_smooth(self.train)
    draw_ts(self.train)
    

    平滑后的训练数据:

    file

    3、将训练数据进行周期性分解

    采用statsmodels工具包:

    from statsmodels.tsa.seasonal import seasonal_decompose
    ​
    decomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)
    # self.ts:时间序列,series类型; 
    # freq:周期,这里为1440分钟,即一天; 
    # two_sided:观察下图2、4行图,左边空了一段,如果设为True,则会出现左右两边都空出来的情况,False保证序列在最后的时间也有数据,方便预测。
    ​
    self.trend = decomposition.trend
    self.seasonal = decomposition.seasonal
    self.residual = decomposition.resid
    decomposition.plot()
    plt.show()
    

    file

    分解图

    第一行observed:原始数据;第二行trend:分解出来的趋势部分;第三行seasonal:周期部分;最后residual:残差部分。

    我采用的是seasonal_decompose的加法模型进行的分解,即 observed = trend + seasonal + residual,另还有乘法模型。

    在建模的时候,只针对trend部分学习和预测,如何将trend的预测结果加工成合理的最终结果?当然是再做加法,后面会详细写。

    模型

    1、训练

    对分解出来的趋势部分单独用arima模型做训练:

    def trend_model(self, order):
        self.trend.dropna(inplace=True)
        train = self.trend[:len(self.trend)-self.test_size]
        #arima的训练参数order =(p,d,q),具体意义查看官方文档,调参过程略。
        self.trend_model = ARIMA(train, order).fit(disp=-1, method='css')
    

    2、预测

    预测出趋势数据后,加上周期数据即作为最终的预测结果,但更重要的是,我们要得到的不是具体的值,而是一个合理区间,当真实数据超过了这个区间,则触发报警,误差高低区间的设定来自刚刚分解出来的残差residual数据:

    d = self.residual.describe()
    delta = d['75%'] - d['25%']
    self.low_error, self.high_error = (d['25%'] - 1 * delta, d['75%'] + 1 * delta)
    

    预测并完成最后的加法处理,得到第七天的预测值即高低置信区间:

        def predict_new(self):
            '''
            预测新数据
            '''
            #续接train,生成长度为n的时间索引,赋给预测序列
            n = self.test_size
            self.pred_time_index= pd.date_range(start=self.train.index[-1], periods=n+1, freq='1min')[1:]
            self.trend_pred= self.trend_model.forecast(n)[0]
            self.add_season()
    ​
    ​
        def add_season(self):
            '''
            为预测出的趋势数据添加周期数据和残差数据
            '''
            self.train_season = self.seasonal[:self.train_size]
            values = []
            low_conf_values = []
            high_conf_values = []
    ​
            for i, t in enumerate(self.pred_time_index):
                trend_part = self.trend_pred[i]
    ​
                # 相同时间点的周期数据均值
                season_part = self.train_season[
                    self.train_season.index.time == t.time()
                    ].mean()
    ​
                # 趋势 + 周期 + 误差界限
                predict = trend_part + season_part
                low_bound = trend_part + season_part + self.low_error
                high_bound = trend_part + season_part + self.high_error
    ​
                values.append(predict)
                low_conf_values.append(low_bound)
                high_conf_values.append(high_bound)
    ​
            # 得到预测值,误差上界和下界
            self.final_pred = pd.Series(values, index=self.pred_time_index, name='predict')
            self.low_conf = pd.Series(low_conf_values, index=self.pred_time_index, name='low_conf')
            self.high_conf = pd.Series(high_conf_values, index=self.pred_time_index, name='high_conf')
    

    3、评估:

    对第七天作出预测,评估的指标为均方根误差rmse,画图对比和真实值的差距:

        md = ModelDecomp(file=filename, test_size=1440)
        md.decomp(freq=1440)
        md.trend_model(order=(1, 1, 3)) # arima模型的参数order
        md.predict_new() 
        pred = md.final_pred
        test = md.test
    ​
        plt.subplot(211)
        plt.plot(md.ts) # 平滑过的训练数据加未做处理的测试数据
        plt.title(filename.split('.')[0])
    ​
        plt.subplot(212)
        pred.plot(color='blue', label='Predict') # 预测值
        test.plot(color='red', label='Original') # 真实值
        md.low_conf.plot(color='grey', label='low') # 低置信区间
        md.high_conf.plot(color='grey', label='high') # 高置信区间
    ​
        plt.legend(loc='best')
        plt.title('RMSE: %.4f' % np.sqrt(sum((pred.values - test.values) ** 2) / test.size))
        plt.tight_layout()
        plt.show()
    

    file

    预测结果

    可以看到,均方根误差462.8,相对于原始数据几千的量级,还是可以的。测试数据中的两个突变的点,也超过了置信区间,能准确报出来。

    结语

    前文提到不同的api形态差异巨大,本文只展示了一个。

    我在该项目中还接触了其他形态的序列,有的有明显的上升或下降趋势;有的开始比较平缓,后面开始增长… …

    但是都属于典型的周期性时间序列,它的核心思想很简单:做好分解,做好预测结果的还原,和置信区间的设置,具体操作可根据具体业务逻辑做调整,祝大家建模愉快。

    原文链接:
    https://www.jianshu.com/p/09e5218f58b4

    文源网络,仅供学习之用,侵删。

    在学习Python的道路上肯定会遇见困难,别慌,我这里有一套学习资料,包含40+本电子书,800+个教学视频,涉及Python基础、爬虫、框架、数据分析、机器学习等,不怕你学不会!
    https://shimo.im/docs/JWCghr8prjCVCxxK/ 《Python学习资料》

    关注公众号【Python圈子】,优质文章每日送达。

    file

    展开全文
  • 这是当初刚进公司时,leader给的一个独立练手小项目,关于时间序列预测,情景比较简单,整个过程实现下来代码也仅100多行,但完成过程中踩了很多坑,觉的有必要分(tu)享(cao)一下。完整代码和样例数据放到了我的...

    这是当初刚进公司时,leader给的一个独立练手小项目,关于时间序列预测,情景比较简单,整个过程实现下来代码也仅100多行,但完成过程中踩了很多坑,觉的有必要分(tu)享(cao)一下。完整代码和样例数据放到了我的github上(文章仅粘贴部分): https://github.com/scarlettgin/cyclical_series_predict

    1、背景

    公司平台上有不同的api,供内部或外部调用,这些api承担着不同的功能,如查询账号、发版、抢红包等等。日志会记录下每分钟某api被访问了多少次,即一个api每天会有1440条记录(1440分钟),将每天的数据连起来观察,有点类似于股票走势的意思。我想通过前N天的历史数据预测出第N+1天的流量访问情况,预测值即作为合理参考,供新一天与真实值做实时对比。当真实流量跟预测值有较大出入,则认为有异常访问,触发报警。

    2、数据探索

    我放了一份样例数据在data文件夹下, 看一下数据大小和结构

    data = pd.read_csv(filename)
    print('size: ',data.shape)
    print(data.head())
    复制代码

    数据大小: 共10080条记录,即10080分钟,七天的数据。 字段含义: date:时间,单位分钟 count:该分钟该api被访问的次数

    画图看一下序列的走势:(一些画图等探索类的方法放在了test_stationarity.py 文件中,包含时间序列图,移动平均图,有兴趣的可以自己尝试下)。

    def draw_ts(timeseries):
        timeseries.plot()
        plt.show()
    
    data = pd.read_csv(path)
    data = data.set_index('date')
    data.index = pd.to_datetime(data.index)
    ts = data['count']
    draw_ts(ts)
    复制代码

    看这糟心的图,那些骤降为0的点这就是我遇到的第一个坑,我当初一拿到这份数据就开始做了。后来折腾了好久才发现,那些骤降为0的点是由于数据缺失,ETL的同学自动补零造成的,沟通晚了(TДT)。

    把坑填上,用前后值的均值把缺失值补上,再看一眼:

    发现这份数据有这样几个特点,在模型设计和数据预处理的时候要考虑到:

    1、这是一个周期性的时间序列,数值有规律的以天为周期上下波动,图中这个api,在每天下午和晚上访问较为活跃,在早上和凌晨较为稀少。在建模之前需要做分解。

    2、我的第二个坑:数据本身并不平滑,骤突骤降较多,而这样是不利于预测的,毕竟模型需要学习好正常的序列才能对未知数据给出客观判断,否则会出现频繁的误报,令气氛变得十分尴尬( ´Д`),所以必须进行平滑处理。

    3、这只是一个api的序列图,而不同的api的形态差距是很大的,毕竟承担的功能不同,如何使模型适应不同形态的api也是需要考虑的问题。

    3、预处理

    3.1 划分训练测试集

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

    class ModelDecomp(object):
        def __init__(self, file, test_size=1440):
            self.ts = self.read_data(file)
            self.test_size = test_size
            self.train_size = len(self.ts) - self.test_size
            self.train = self.ts[:len(self.ts)-test_size]
            self.test = self.ts[-test_size:]
    复制代码

    3.2 对训练数据进行平滑处理

    消除数据的毛刺,可以用移动平均法,我这里没有采用,因为我试过发现对于我的数据来说,移动平均处理完后并不能使数据平滑,我这里采用的方法很简单,但效果还不错:把每个点与上一点的变化值作为一个新的序列,对这里边的异常值,也就是变化比较离谱的值剃掉,用前后数据的均值填充,注意可能会连续出现变化较大的点:

    def _diff_smooth(self, ts):
        dif = ts.diff().dropna() # 差分序列
        td = dif.describe() # 描述性统计得到:min,25%,50%,75%,max值
        high = td['75%'] + 1.5 * (td['75%'] - td['25%']) # 定义高点阈值,1.5倍四分位距之外
        low = td['25%'] - 1.5 * (td['75%'] - td['25%']) # 定义低点阈值,同上
    
        # 变化幅度超过阈值的点的索引
        forbid_index = dif[(dif > high) | (dif < low)].index 
        i = 0
        while i < len(forbid_index) - 1:
            n = 1 # 发现连续多少个点变化幅度过大,大部分只有单个点
            start = forbid_index[i] # 异常点的起始索引
            while forbid_index[i+n] == start + timedelta(minutes=n):
                n += 1
            i += n - 1
    
            end = forbid_index[i] # 异常点的结束索引
            # 用前后值的中间值均匀填充
            value = np.linspace(ts[start - timedelta(minutes=1)], ts[end + timedelta(minutes=1)], n)
            ts[start: end] = value
            i += 1
    
    self.train = self._diff_smooth(self.train)
    draw_ts(self.train)
    复制代码

    平滑后的训练数据:

    3.3 将训练数据进行周期性分解

    采用statsmodels工具包:

    from statsmodels.tsa.seasonal import seasonal_decompose
    
    decomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)
    # self.ts:时间序列,series类型; 
    # freq:周期,这里为1440分钟,即一天; 
    # two_sided:观察下图2、4行图,左边空了一段,如果设为True,则会出现左右两边都空出来的情况,False保证序列在最后的时间也有数据,方便预测。
    
    self.trend = decomposition.trend
    self.seasonal = decomposition.seasonal
    self.residual = decomposition.resid
    decomposition.plot()
    plt.show()
    复制代码

    第一行observed:原始数据;第二行trend:分解出来的趋势部分;第三行seasonal:周期部分;最后residual:残差部分。 我采用的是seasonal_decompose的加法模型进行的分解,即 observed = trend + seasonal + residual,另还有乘法模型。在建模的时候,只针对trend部分学习和预测,如何将trend的预测结果加工成合理的最终结果?当然是再做加法,后面会详细写。

    4、模型

    4.1 训练

    对分解出来的趋势部分单独用arima模型做训练:

    def trend_model(self, order):
        self.trend.dropna(inplace=True)
        train = self.trend[:len(self.trend)-self.test_size]
        #arima的训练参数order =(p,d,q),具体意义查看官方文档,调参过程略。
        self.trend_model = ARIMA(train, order).fit(disp=-1, method='css')
    复制代码

    4.2 预测

    预测出趋势数据后,加上周期数据即作为最终的预测结果,但更重要的是,我们要得到的不是具体的值,而是一个合理区间,当真实数据超过了这个区间,则触发报警,误差高低区间的设定来自刚刚分解出来的残差residual数据:

    d = self.residual.describe()
    delta = d['75%'] - d['25%']
    self.low_error, self.high_error = (d['25%'] - 1 * delta, d['75%'] + 1 * delta)
    复制代码

    预测并完成最后的加法处理,得到第七天的预测值即高低置信区间:

            def predict_new(self):
                '''
                预测新数据
                '''
                #续接train,生成长度为n的时间索引,赋给预测序列
                n = self.test_size
                self.pred_time_index= pd.date_range(start=self.train.index[-1], periods=n+1, freq='1min')[1:]
                self.trend_pred= self.trend_model.forecast(n)[0]
                self.add_season()
    
    
            def add_season(self):
                '''
                为预测出的趋势数据添加周期数据和残差数据
                '''
                self.train_season = self.seasonal[:self.train_size]
                values = []
                low_conf_values = []
                high_conf_values = []
                
                for i, t in enumerate(self.pred_time_index):
                    trend_part = self.trend_pred[i]
                    
                    # 相同时间点的周期数据均值
                    season_part = self.train_season[
                    	self.train_season.index.time == t.time()
                    	].mean()
        
                    # 趋势 + 周期 + 误差界限
                    predict = trend_part + season_part
                    low_bound = trend_part + season_part + self.low_error
                    high_bound = trend_part + season_part + self.high_error
                    
                    values.append(predict)
                    low_conf_values.append(low_bound)
                    high_conf_values.append(high_bound)
    
                # 得到预测值,误差上界和下界
                self.final_pred = pd.Series(values, index=self.pred_time_index, name='predict')
                self.low_conf = pd.Series(low_conf_values, index=self.pred_time_index, name='low_conf')
                self.high_conf = pd.Series(high_conf_values, index=self.pred_time_index, name='high_conf')
    复制代码

    4.3 评估:

    对第七天作出预测,评估的指标为均方根误差rmse,画图对比和真实值的差距:

    	md = ModelDecomp(file=filename, test_size=1440)
    	md.decomp(freq=1440)
    	md.trend_model(order=(1, 1, 3)) # arima模型的参数order
    	md.predict_new() 
    	pred = md.final_pred
    	test = md.test
    
    	plt.subplot(211)
    	plt.plot(md.ts) # 平滑过的训练数据加未做处理的测试数据
    	plt.title(filename.split('.')[0])
    
    	plt.subplot(212)
    	pred.plot(color='blue', label='Predict') # 预测值
    	test.plot(color='red', label='Original') # 真实值
    	md.low_conf.plot(color='grey', label='low') # 低置信区间
    	md.high_conf.plot(color='grey', label='high') # 高置信区间
    
    	plt.legend(loc='best')
    	plt.title('RMSE: %.4f' % np.sqrt(sum((pred.values - test.values) ** 2) / test.size))
    	plt.tight_layout()
    	plt.show()
    复制代码

    可以看到,均方根误差462.8,相对于原始数据几千的量级,还是可以的。测试数据中的两个突变的点,也超过了置信区间,能准确报出来。

    5、结语

    前文提到不同的api形态差异巨大,本文只展示了一个,我在该项目中还接触了其他形态的序列,有的有明显的上升或下降趋势;有的开始比较平缓,后面开始增长... ... ,但是都属于典型的周期性时间序列,它的核心思想很简单:做好分解,做好预测结果的还原,和置信区间的设置,具体操作可根据具体业务逻辑做调整,祝大家建模愉快:-D。

    展开全文
  • 在我们的研究中,我们考虑了对印度旁遮普省降雨数据进行统计分析的季节性和周期性时间序列模型。 在本研究论文中,我们应用季节性自回归综合移动平均和周期自回归模型来分析旁遮普省的降雨数据。 为了评估模型识别...
  • GPS坐标时间序列周期性分析.pdf
  • 背景任何事物在两个不同时刻都不可能保持完全相同的状态,但很多变化往往存在着一定的规律,例如 24 小时日出日落,潮起潮落,这些现象通常称为「周期」。周期性,指时间序列中呈现出来的围绕长期趋...

    背景

    任何事物在两个不同时刻都不可能保持完全相同的状态,但很多变化往往存在着一定的规律,例如 24 小时日出日落,潮起潮落,这些现象通常称为「周期」

    周期性,指时间序列中呈现出来的围绕长期趋势的一种波浪形或振荡式变动。准确提取周期信息,不仅能反映当前数据的规律,应用于相关场景,还可以预测未来数据变化趋势。

    时间序列示例

    一般而言,时间序列周期性分为三种:

    • 「符号性周期」,例如序列 fbcnfkgbfopsf 周期为 4;

    • 「部分周期性」,例如序列 ansdcdmncdcacdascdmccd 周期为 4;

    • 「分段周期性」,例如上面给定的时间序列即为分段周期性;

    针对时间序列的周期性检测问题,这篇文章主要介绍「傅里叶变换」「自相关系数」两种方法及其在实际数据中的效果;

    傅里叶变换

    傅里叶变换是一种将时域、空域数据转化为频域数据的方法,任何波形(时域)都可以看做是不同振幅、不同相位正弦波的叠加(频域),详细介绍可以参考:An Interactive Guide To The Fourier Transform,???????? 此处放上经典图镇场!

    傅里叶变换

    对于一条具备周期性的时间序列,它本身就很接近正弦波,所以它包含一个显著的正弦波,周期就是该正弦波的周期,而这个正弦波可以通过傅里叶变换找到,它将时序数据展开成三角函数的线性组合,得到每个展开项的系数,就是傅里叶系数。傅里叶系数越大,表明它所对应的正弦波的周期就越有可能是这份数据的周期。

    自相关系数

    自相关系数(Autocorrelation Function)度量的是同一事件不同时间的相关程度,不同相位差(lag)序列间的自相关系数可以用 Pearson 相关系数计算。其数学表达如下:

    其中 表示相位 的数据延迟 「lag operator」

    自相关系数

    当序列存在周期性时,遍历足够多的相位差,一定可以找到至少一个足够大的自相关系数,而它对应的相位差就是周期。所以对于检测时序周期来说,只需找到两个自相关系数达到一定阈值的子序列,它们起始时间的差值就是我们需要的周期。

    实例说明

    为了保证结果的可靠性,可以将傅里叶分析和自相关系数结合起来判断周期性。主要思路是:先通过傅里叶变换找到可能的周期,再用自相关系数做排除,从而得到最可能的周期。

    给定一份周期性数据,时间间隔为 5 min。从这份数据中可以看出数据大体上具有周期为 1 day。

    示例数据

    下面使用傅里叶变换估计周期,代码如下所示

    from scipy.fftpack import fft, fftfreq
    
    fft_series = fft(data["value"].values)
    power = np.abs(fft_series)
    sample_freq = fftfreq(fft_series.size)
    
    pos_mask = np.where(sample_freq > 0)
    freqs = sample_freq[pos_mask]
    powers = power[pos_mask]
    
    top_k_seasons = 3
    # top K=3 index
    top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:]
    top_k_power = powers[top_k_idxs]
    fft_periods = (1 / freqs[top_k_idxs]).astype(int)
    
    print(f"top_k_power: {top_k_power}")
    print(f"fft_periods: {fft_periods}")
    

    取 top-3 振幅值为top_k_power: [ 614.8105282 890.33273899 1831.167168 ] 及其对应的周期 fft_periods: [ 72 278 292] 。???? 数据间隔为 5 min 所以真实周期应为 288,从傅里叶变换即可看出估计值 292 已经非常接近真实值。

    现在来计算自相关系数,代码如下所示:

    from statsmodels.tsa.stattools import acf
    
    # Expected time period
    for lag in fft_periods:
        # lag = fft_periods[np.abs(fft_periods - time_lag).argmin()]
        acf_score = acf(data["value"].values, nlags=lag)[-1]
        print(f"lag: {lag} fft acf: {acf_score}")
    
    
    expected_lags = np.array([timedelta(hours=12)/timedelta(minutes=5), timedelta(days=1)/timedelta(minutes=5), timedelta(days=7)/timedelta(minutes=5)]).astype(int)
    for lag in expected_lags:
        acf_score = acf(data["value"].values, nlags=lag, fft=False)[-1]
        print(f"lag: {lag} expected acf: {acf_score}")
    
    

    对应的输出如下:

    lag: 72 fft acf: 0.07405431832776994
    lag: 278 fft acf: 0.7834457453491087
    lag: 292 fft acf: 0.8259822269757922
    lag: 144 expected acf: -0.5942986094704665
    lag: 288 expected acf: 0.8410792774898174
    lag: 2016 expected acf: 0.5936030431473589
    

    通过自相关系数来得到显著分数最大值对应的周期,得出的结果为 292;

    此处实验补充了预设的三个周期值:12 hour、1 day、7 day,发现算出来还是周期 288 对应的相关分数最大,但是傅里叶变换没有估计出周期值 ????

    综上,这个小故事告诉我们:你算出来的还不如我预设的值呢!直接根据先验知识「预设周期」然后计算自相关系数就行了!


    建议阅读:

    高考失利之后,属于我的大学本科四年

    【资源分享】对于时间序列,你所能做的一切.

    【时空序列预测第一篇】什么是时空序列问题?这类问题主要应用了哪些模型?主要应用在哪些领域?

    【AI蜗牛车出品】手把手AI项目、时空序列、时间序列、白话机器学习、pytorch修炼

    公众号:AI蜗牛车
    
    保持谦逊、保持自律、保持进步
    
    
    
    个人微信
    备注:昵称+学校/公司+方向
    如果没有备注不拉群!
    拉你进AI蜗牛车交流群
    
    
    
    
    展开全文
  • 我们表明,可以通过几乎周期性相关(APC)时间序列的离散频谱以非参数方式对商业周期的某些特征进行建模。 根据该频谱的估计特征,通过过滤来提取业务周期。 作为说明,我们在波兰经济的工业生产指数中描述了商业...
  • 浅谈LSTM对于周期时间序列数据的预测

    万次阅读 多人点赞 2019-03-06 12:33:40
    注:本文章主要针对的是长周期时间序列数据(10000-40000条为一个周期的数据)预测 产生训练和测试数据 我们需要做的是产生周期为20000条/周期的sin函数时间序列(用于训练) 以及周期为40000条/周期的sin...
  • 元胞自动机是结构简单但行为复杂多样的离散动力系统.本文以模拟y0号初等元胞自动机为基础,证明了126号初等元胞自动机下任意有限初始条件迭代产生的时间序列是非周期的.
  • 序列分解1、非季节性时间序列分解 移动平均MA(Moving Average)①SAM(Simple Moving Average) 简单移动平均,将时间序列上前n个数值做简单的算术平均。 SMAn=(x1+x2+…xn)/n②WMA(Weighted Moving Average) ...
  • 时间序列算法分解一般是指把一个时间序列分解成:趋势序列,周期序列,残差序列。 时间序列分解算法最广为人知的可能是STL算法。它只能分解出一个周期序列。有很多博客和文章叙述了STL分解的原理,例如博客 时间...
  • 对商业周期的分析需要提取时间序列周期性成分,该时间序列通常也受到诸如潜在趋势或噪声等其他因素的影响。本文介绍了一些在最近的文献中用于从给定系列中提取商业周期的方法。它基于Stock and Watson(1999)在...
  • Matlab时间序列分析

    万次阅读 多人点赞 2018-11-13 18:53:46
    在引入时间序列前,先介绍几个matlab函数 matlab中的gallery函数简析 Matlab 中的 gallery 函数是一个测试矩阵生成函数。当我们需要对某些算法进行测试的时候,可以利用gallery函数来生成各种性质的测试矩阵。其用法...
  • R-时间序列-分解季节性时间序列

    千次阅读 2018-10-16 17:06:08
    1.季节性时间序列 包含:长期趋势Trend,季节趋势Seasonal,周期循环Circle,随机项Random 这里分解为相加模型X=T+S+C+R   在对时间序列进行分解之前,应该对序列进行检验:(下次写) 2.decompose()函数 将...
  • 时间序列用户生命周期的聚类方法

    千次阅读 2015-08-05 14:11:07
    这是一篇旧文:时间序列用户生命周期的聚类方法
  • 使用SPSS做时间序列预测时,如何自定义日期型数据周期/时间数据周期?问题的产生解决方案的来源:SPSS自带的语法参考DATE命令的基本使用方法两个示例讲解写好命令以后如何运行需要注意的点如果帮到了你,能留言鼓励...
  • 我们考虑一些时间序列,例如道路上的交通流量, > plot(T,X,type="l") > reg=lm(X~T) > abline(reg,col="red") 如果存在趋势,我们应该将其删除,然后处理残差 > Y=residuals(reg) > acf(Y...
  • 对商业周期的分析需要提取时间序列周期性成分,该时间序列通常也受到诸如潜在趋势或噪声等其他因素的影响。本文介绍了一些在最近的文献中用于从给定系列中提取商业周期的方法。它基于Stock and Watson(1999)在...
  • 时间序列模型的三个重要概念

    千次阅读 2017-08-02 14:57:09
    时间序列在量化投资中具有广泛的...探索金融时间序列的趋势和周期性 时间序列与其他变量的内在关系,为策略提供辅助和增强 不同时间序列之间的关系,发现新的策略 波动率建模,期权相关的策略 在时间序
  • 伪周期时间序列是一种广泛存在的数据形式,它具有伪周期性、非平稳性和特征值等特征。对这类时间序列进行预测,具有很强的研究和应用意义。然而,目前的相关研究对伪周期时间序列的关注度不足,一些已有的时间序列...
  • 经济时间序列的分析通常需要提取其周期性成分。这篇文章介绍了一些方法,可用于将时间序列分解为它们的不同部分。它基于《宏观经济学手册》中Stock和Watson(1999)关于商业周期的章节,但也介绍了一些较新的方法,...
  • 我的同僚经济学家的“错误信念”,我们称之为轰动效应的“神话”,被显示出来:有关商业周期时间序列和预测的“神话”。 尽管从长远来看,我们都死了,但经济历史……记住。 “贸易周期”理论最终成为“商业周期...
  • 使用时间延迟飞秒激光双脉冲序列在6H-SiC单晶上产生激光诱导的周期性表面结构
  • 季节变动(S):是时间序列在一年内重复出现的周期性波动。 循环波动©:是时间序列呈现出得非固定长度的周期性变动。 随机因素(I):是时间序列中除去长期趋势、季节变动和循环波动之后的随机波动。不规则波动通常总是...
  • 时间序列预测算法总结

    万次阅读 多人点赞 2018-10-18 10:30:48
    time series data mining 主要包括decompose(分析数据的各个成分,例如趋势,周期性),prediction(预测未来的值),classification(对有序数据序列的feature提取与分类),clustering(相似数列聚类)等。...
  • 时间序列预测,非季节ARIMA及季节SARIMA

    万次阅读 多人点赞 2019-03-24 21:55:00
    我们将首先介绍和讨论自相关,平稳和季节的概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。 介绍 时间序列提供了预测未来价值的机会。基于以前的价值观,可以使用时间序列来预测经济,天气和...
  • 时间序列相似度量-DTW

    万次阅读 多人点赞 2019-01-25 16:44:59
    最近项目中遇到求解时间序列相似问题,这里序列也可以看成向量。在传统算法中,可以用余弦相似度和pearson相关系数来描述两个序列的相似度。但是时间序列比较特殊,可能存在两个问题: 两段时间序列长度不同。...
  • 时间序列

    千次阅读 2018-09-02 00:33:43
    时间序列的定义 所谓时间序列就是按照时间的顺序记录的一列有序数据。对时间序列进行观察、研究、找寻他发展变化的规律,预测他将来的走势就是时间序列分析,时间序列分析方法只适用于近期与短期的预测。 相关特征...
  • 时间序列分析 - ARMA, ARIMA, SARIMA

    千次阅读 2019-03-06 17:19:26
    ARIMA: 针对非平稳非周期性时间序列分析 SARIMA: 针对非平稳周期性时间序列分析。 【自协方差与自相关系数】 时间序列在t时刻记作Xt,在s时刻记作Xs,那么这两个时刻对应的时间序列的自协方差的计算公式为: ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 123,232
精华内容 49,292
关键字:

周期性时间序列