精华内容
下载资源
问答
  • 背景任何事物在两个不同时刻都不可能保持完全相同的状态,但很多变化往往存在着一定的规律,例如 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蜗牛车交流群
    
    
    
    
    展开全文
  • 这是当初刚进公司时,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。

    展开全文
  • 时间序列的平稳性、随机性检验 在拿到时间序列数据后,首先要对数据的随机性和平稳性进行检测, 这两个检测时间序列预测的重要部分。根据不同检测结果需要采取不同的分析方法。 为什么时间序列要求平稳性呢?平稳...

    时间序列系列文章:

    时间序列(一):时间序列数据与时间序列预测模型
    时间序列(二):时间序列平稳性检测
    时间序列(三):ARIMA模型实战

    在上一篇文章时间序列(一):时间序列数据与时间序列预测模型中我们介绍了时间序列及一些时间序列预测模型。我们可以看到在进行预测时有一些模型表现较好,而另一些模型的预测结果却不尽人意。这是因为不同的时间序列模型对原始数据的要求是不同的,例如之前提到的ARIMA模型,要求时间序列数据平稳,否则得出的预测结果就会相差较大。本篇文章我们介绍时间序列的平稳性、随机性检验及相关时间序列数据处理方法。

    时间序列的平稳性、随机性检验

    在拿到时间序列数据后,首先要对数据的随机性和平稳性进行检测, 这两个检测是时间序列预测的重要部分。根据不同检测结果需要采取不同的分析方法。

    为什么时间序列要求平稳性呢?平稳性就是要求由样本拟合出的曲线在未来一段时间内仍然能够以现有的形态和趋势发展下去,这样预测结果才会有意义。

    对于平稳声序列, 它的均值和方差是常数, 现已有一套非常成熟的平稳序列的建模方法。 通常是建立一个线性模型来拟合该序列的发展 借此提取该序列的有用信息。

    对于非平稳序列, 由于它的均值和方差不稳定, 处理方法一般是将其转变为平稳序列,这样就可以应用有关平稳时间序列的分析方法, 如建立 ARIMA模型来进行相应的研究,或者分解趋势与季节性等并根据情况应用指数平滑模型等。

    对于纯随机序列, 又称为白噪声序列, 序列的各项之间没有任何相关关系, 序列在进行完全无序的随机波动, 可以终止对该序列的分析。 白噪声序列是没有信息可提取的平稳序列。

    在讲解平稳性和随机性的定义之前,我们先介绍一下时间序列中常用的几个特征统计量。

    时间序列的特征统计量

    对于一个时间序列任意时刻的序列值 { X t , t ∈ T } \left\{ X _ { t } , t \in T \right\} {Xt,tT},任意时刻的序列值 X t X _ { t } Xt都是一个随机变量,记其分布函数为 F t ( x ) F _ { t } ( x ) Ft(x),则其特征统计量均值、方差、自协方差函数、自相关系数的定义分别如下:

    均值: 表示时间序列在各个时刻取值的平均值,其定义如下:
    μ t = E X t = ∫ − ∞ ∞ x d F t ( x ) \mu _ { t } = E X _ { t } = \int _ { - \infty } ^ { \infty } x \mathrm { d } F _ { t } ( x ) μt=EXt=xdFt(x)

    方差: 表示时间序列在各个时刻围绕其均值波动的平均程度,其定义如下:
    σ t 2 = D X t = E ( X t − μ t ) 2 = ∫ − ∞ ∞ ( x − μ t ) 2 d F t ( x ) \sigma _ { t } ^ { 2 } = D X _ { t } = E \left( X _ { t } - \mu _ { t } \right) ^ { 2 } = \int _ { - \infty } ^ { \infty } \left( x - \mu _ { t } \right) ^ { 2 } \mathrm { d } F _ { t } ( x ) σt2=DXt=E(Xtμt)2=(xμt)2dFt(x)

    自协方差 : 表示时间序列任意两个时刻直接的相关性,任取 t , s ∈ T t , s \in T t,sT,则其定义如下:
    γ ( t , s ) = E [ ( X t − μ t ) ( X s − μ s ) ] \gamma ( t , s ) = E \left[ \left( X _ { t } - \mu _ { t } \right) \left( X _ { s } - \mu _ { s } \right) \right] γ(t,s)=E[(Xtμt)(Xsμs)]

    自相关系数: 同自协方差函数,其定义如下:

    ρ ( t , s ) = γ ( t , s ) D X t ⋅ D X s \rho ( t , s ) = \frac { \gamma ( t , s ) } { \sqrt { D X _ { t } \cdot D X _ { s } } } ρ(t,s)=DXtDXs γ(t,s)

    平稳时间序列的定义与检验

    平稳时间序列的定义

    平稳时间序列按照限定条件的严格程度可以分为以下两种类型:

    严平稳时间序列: 指时间序列的所有统计性质不会随着时间的推移而发生变化,即其联合概率分布在任何时间间隔都是相同的。设 { X t } \left\{ X _ { t } \right\} {Xt}为一时间序列,对任意的正整数 m m m,任取 t 1 , t 2 , ⋯   , t m ∈ T t _ { 1 } , t _ { 2 } , \cdots , t _ { m } \in T t1,t2,,tmT,对任意整数 τ \tau τ,有:

    F t 1 , t 2 , ⋯   , t m ( x 1 , x 2 , ⋯   , x m ) = F t 1 + τ , t 2 + τ , ⋯   , t m + τ ( x 1 , x 2 , ⋯   , x m ) F _ { t _ { 1 } , t _ { 2 } , \cdots , t _ { m } } \left( x _ { 1 } , x _ { 2 } , \cdots , x _ { m } \right) = F _ { t _ { 1 + \tau } , t _ { 2 + \tau } , \cdots , t _ { m + \tau } } \left( x _ { 1 } , x _ { 2 } , \cdots , x _ { m } \right) Ft1,t2,,tm(x1,x2,,xm)=Ft1+τ,t2+τ,,tm+τ(x1,x2,,xm)

    则称时间序列 { X t } \left\{ X _ { t } \right\} {Xt}为严平稳时间序列。

    宽平稳时间序列: 宽平稳时间序列则认为只要时间序列的低阶距(二阶)平稳,则该时间序列近似平稳。如果时间序列 { X t } \left\{ X _ { t } \right\} {Xt}满足以下三个条件:

    • 任取 t ∈ T t \in T tT,有 E X t 2 < ∞ EX _ { t }^ { 2 } <∞ EXt2<
    • 任取 t ∈ T t \in T tT,有 E X t = μ E X _ { t } = \mu EXt=μ,其中 μ \mu μ为常数;
    • 任取 t , s , k ∈ T t , s , k \in T t,s,kT, k + s − t ∈ T k + s - t \in T k+stT,有 γ ( t , s ) = γ ( k , k + s − t ) \gamma ( t , s ) = \gamma ( k , k + s - t ) γ(t,s)=γ(k,k+st)

    在现实生活中,时间序列是很难满足严平稳时间序列的要求的,因此,一般所讲的平稳时间序列在默认情况下都是指宽平稳时间序列。根据宽平稳时间序列的条件,我们可以容易得到宽平稳时间序列所具有的性质:

    • 均值为常数,即: E X t = μ , ∀ t ∈ T E X _ { t } = \mu , \quad \forall t \in T EXt=μ,tT
    • 方差也为均值,即: D X t = γ ( t , t ) = γ ( 0 ) , ∀ t ∈ T D X _ { t } = \gamma ( t , t ) = \gamma ( 0 ) , \quad \forall t \in T DXt=γ(t,t)=γ(0),tT
    • 自协方差函数和自相关系数只依赖于时间的平移长度,而与时间的起点无关。即: γ ( t , s ) = γ ( k , k + s − t ) , ∀ t , s , k ∈ T \gamma ( t , s ) = \gamma ( k , k + s - t ) , \quad \forall t , s , k \in T γ(t,s)=γ(k,k+st),t,s,kT
      因此,可以记 γ ( k ) \gamma ( k ) γ(k)为时间序列${ X _ { t } $的延迟k自协方差函数。

    由于平稳时间序列具有这些优良性质,因此,对于一个平稳时间序列来说,其待估计的参数量就变得少了很多,因为他们的均值、方差都是一样的,因此,可以利用全部的样本来估计总体的均值和方差,即:
    μ ^ = x ‾ = ∑ i = 1 n x i n γ ^ ( 0 ) = ∑ t = 1 n ( x t − x ‾ ) 2 n − 1 \widehat { \mu } = \overline { x } = \frac { \sum _ { i = 1 } ^ { n } x _ { i } } { n } \\ \widehat { \gamma } ( 0 ) = \frac { \sum _ { t = 1 } ^ { n } \left( x _ { t } - \overline { x } \right) ^ { 2 } } { n - 1 } μ =x=ni=1nxiγ (0)=n1t=1n(xtx)2
    这也是为什么说当拿到一个时间序列后,需要对其进行平稳性检验。

    平稳时间序列的检验

    那么,当拿到一个时间序列后,应该如何对其进行平稳性的检验呢?目前,对时间序列的平稳性检验主要有两种方法,一种是图检法,即根据时序图和自相关图进行直观判断,另一种是构造检验统计量的方法,有单位根检验法等方法。

    图检法

    对于图检法,我们一般绘制时间序列的时序图,考虑以下三个图形:

    在这里插入图片描述
    在第一幅图中,我们可以清楚地看到,均值随时间而变化(增加),呈现上升的趋势。因此,这是一个非平稳序列。平稳序列不应该呈现出随时间变化的趋势。

    第二幅图显然看不到序列的趋势,但序列的变化是一个时间的函数。正如前面提到的,平稳序列的方差必须是一个常数。

    再来看第三幅图,随着时间的增加,序列传播后变得更近,这意味着协方差是时间的函数。

    所以上述三个例子均是非平稳时间序列.

    再看下面的时序图:
    在这里插入图片描述
    在这张图中,均值、方差和协方差都是常数,这就是平稳时间序列

    另一方面,我们也可以通过自相关图来进行检验,对于平稳时间序列,其自相关图一般随着阶数的递增,自相关系统会迅速衰减至0附近,而非平稳时间序列则可能存在先减后增或者周期性波动等变动。如下图所示,该时间序列随着阶数的递增,自相关系数先减后增,因此,可以判断该时间序列不是平稳时间序列。
    在这里插入图片描述

    统计检验

    可以利用统计检验来代替目视检验:比如单位根平稳检验。单位根表名给定序列的统计特性(均值,方差和协方差)不是时间的常数,这是平稳时间序列的先决条件。下面是它的数学解释:

    假设我们有一个时间序列:
    y t = a ∗ y t − 1 + ε t y_t=a*y_{t-1 }+ε_t yt=ayt1+εt
    其中 y t y_t yt是t时刻的数据值, ε t ε_t εt 是误差项。
    仙子我们需要利用 y t − 1 y_{t-1 } yt1的值来计算 y t y_t yt,即:
    y t − 1 = a ∗ y t − 2 + ε t − 1 y_{t-1}=a*y_{t-2}+ε_{t-1} yt1=ayt2+εt1
    如果利用所有的观察值, y t y_{t} yt的值将是:
    y t = a n ∗ y t − n + ∑ ε t − i ∗ a i y_{t}=a^n*y_{t-n}+ \sumε_{t-i}*a^i yt=anytn+εtiai

    假设在上述方程中a的值为1(单位),则预测值将等于 y t − n y_{t-n} ytn 和从 t − n t-n tn t t t的所有误差之和,这意味着方差将随着时间的推移而增大,这就是时间序列中的单位根。众所周知,平稳时间序列的方差不能是时间的函数。单元根检验通过检查a=1的值来检查序列中是否存在单位根。以下是两个最常用的单位根平稳检测方法:

    ADF(增补迪基-福勒)检验

    迪基-福勒(Dickey Fuller)检验是最流行的统计检验方法之一,可以用它来确定序列中单位根的存在,从而帮助判断序列是否是平稳。ADF检验是对DF检验的扩展。这一检验的原假设与备选假设如下:

    • 原假设: 序列有一个单位根(序列非平稳)
    • 备选假设: 该序列没有单位根。(序列平稳)

    单位根是什么呢?当一个自回归过程中: y t = b y t − 1 + a + ϵ t y_{t} = by_{t-1} + a + \epsilon _{t} yt=byt1+a+ϵt,如果滞后项系数b为1,就称为单位根。当单位根存在时,自变量和因变量之间的关系具有欺骗性,因为残差序列的任何误差都不会随着样本量(即时期数)增大而衰减,也就是说模型中的残差的影响是永久的。这种回归又称作伪回归。如果单位根存在,这个过程就是一个随机漫步(random walk)。ADF检验就是判断序列是否存在单位根:如果序列平稳,就不存在单位根;否则,就会存在单位根.

    在Python中使用ADF检验可以在statsmodels中使用adfuller函数。在下面的代码中我们加上标题与输出值对应的名称:

    #定义ADF输出格式化函数
    from statsmodels.tsa.stattools import adfuller
    def adf_test(timeseries):
        print ('ADF检验结果:')
        dftest = adfuller(timeseries, autolag='AIC')
        dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','Number of Lags Used','Number of Observations Used'])
        for key,value in dftest[4].items():
           dfoutput['Critical Value (%s)'%key] = value
        print (dfoutput)
    
    #对数据集使用ADF检验
    adf_test(train['#Passengers'])
    

    检验结果如下:

    ADF检验结果:
    Test Statistic                   0.815369
    p-value                          0.991880
    Number of Lags Used             13.000000
    Number of Observations Used    130.000000
    Critical Value (1%)             -3.481682
    Critical Value (5%)             -2.884042
    Critical Value (10%)            -2.578770
    dtype: float64
    

    ADF的结果主要看以下两个方面:

    • Test Statistic的值如果比Critical Value (5%)小则满足稳定性需求.
    • p-value越低(理论上需要低于0.05)证明序列越稳定。

    在上面的例子中,test statistic > Critical Value (5%) ,这意味着序列不是平稳的。同时p值为0.99,这证实了我们最初在目视检测中观察的结果。

    随机性(白噪声)的定义与检验

    随机性时间序列的定义

    通过对时间序列进行平稳性检验后,我们可以将时间序列分为平稳时间序列和非平稳时间序列,对于非平稳时间序列,一般需要将其转化为平稳时间序列再进行分析,具体的转化方法随后再讲。而对于平稳时间序列,我们知道其有一个性质,即自协方差函数和自相关系数只依赖于时间间隔,而与起点无关,对于相同的时间间隔,其自协方差函数和自相关系数为一个常数,那么,就存在一种情况,当该常数为0时,照样满足平稳时间序列的条件,而此时序列之间的相关性则为0,即序列之间不相关,那么,这时我们的分析即可结束,因为对于一个毫无相关的序列,我们没法从中挖掘出可用的规律,此时的序列即为随机性时间序列,也称为白噪声序列。
        
    对于时间序列 X t X _ { t } Xt,如果满足:

    • 任取 t ∈ T t \in T tT,有 E X t = μ E X _ { t } = \mu EXt=μ
    • 任取 t , s ∈ T t , s \in T t,sT,有
      γ ( t , s ) = σ 2 ; t = s 0 ; t ≠ s \gamma ( t , s ) = \begin{array} { l l } { \sigma ^ { 2 }}; { t = s } \\ { 0 }; { t \neq s } \end{array} γ(t,s)=σ2;t=s0;t=s

    则称该时间序列为纯随机序列或白噪声序列,简记为 X t ∼ W N ( μ , σ 2 ) X _ { t } \sim W N \left( \mu , \sigma ^ { 2 } \right) XtWN(μ,σ2)。我们可以发现,其实白噪声序列的性质与平稳时间序列的性质一样,其均值和方差均为常数,只是自协方差函数或自相关系数为0,因此,该序列的任何两项之间不存在相关性,无法从中得到任何有用的信息,此时分析可以停止。

    纯随机性(白噪声)检验

    对于纯随机性序列,一般通过构建统计量的方法来检验。我们知道,白噪声序列除了0阶自相关系数外,即方差,其他阶的自相关系数应该均为0,因此,我们可以提出下面这样一个假设:

    H 0 : ρ 1 = ρ 2 = ⋯ = ρ m = 0 , ∀ m ⩾ 1 H 1 : 至 少 存 在 某 个 ρ k ≠ 0 , ∀ m ⩾ 1 , k ⩽ m H _ { 0 } : \rho _ { 1 } = \rho _ { 2 } = \cdots = \rho _ { m } = 0 , \quad \forall m \geqslant 1\\ H _ { 1 } :至少存在某个\rho _ { k } \neq 0 , \quad \forall m \geqslant 1 , k \leqslant m H0:ρ1=ρ2==ρm=0,m1H1:ρk=0,m1,km

    因此,围绕该假设,我们可以构建统计量进行检验,常用的统计量有Q统计量和LB统计量,其计算公式分别如下:

    Q = n ∑ k = 1 m ρ ^ k 2 L B = n ( n + 2 ) ∑ k = 1 m ( ρ ^ k 2 n − k ) Q = n \sum _ { k = 1 } ^ { m } \widehat { \rho } _ { k } ^ { 2 }\\L B = n ( n + 2 ) \sum _ { k = 1 } ^ { m } \left( \frac { \widehat { \rho } _ { k } ^ { 2 } } { n - k } \right) Q=nk=1mρ k2LB=n(n+2)k=1m(nkρ k2)

    其中, n n n为序列的观察期数, m m m为指定延迟期数, k k k为延迟阶数,Box和Pierce证明这两个统计量均服从自由度为m mm的卡方分布,当统计量大于 χ 1 − α 2 ( m ) \chi _ { 1 - \alpha } ^ { 2 } ( m ) χ1α2(m)或者P值小于 α α α 时,则认为可以拒绝原假设,即认为该序列是非随机序列

    我们可以使用acorr_ljungbox函数进行数据的纯随机性检验.语法为:

    acorr_ljungbox(x, lags=None, boxpierce=False) # 数据的纯随机性检验函数
    
    lags:为延迟期数,如果为整数,则是包含在内的延迟期数,如果是一个列表或数组,那么所有时滞都包含在列表中最大的时滞中
    
    boxpierce:为True时表示除开返回LB统计量还会返回Box和Pierce的Q统计量
    
    返回值:
    
    lbvalue:测试的统计量
    
    pvalue:基于卡方分布的p统计量
    
    bpvalue:((optionsal), float or array) – 基于 Box-Pierce 的检验的p统计量
    
    bppvalue:((optional), float or array) – 基于卡方分布下的Box-Pierce检验的p统计量
    

    代码实现:

    from statsmodels.stats.diagnostic import acorr_ljungbox
    acorr_ljungbox(train['#Passengers'])
    

    输出检验结果中会返回两个值:
    lbvalue: 测试的统计量 和 pvalue: 基于卡方分布的p统计量。如果p-value>0.05则可判断为白噪声序列

    时间序列的平稳化

    在熟悉了平稳性的概念及其不同的类型之后,接下来可以对序列进行平稳化操作。平稳化的方法有以下几种:

    差分法

    在该方法中,计算序列中连续项的差值。执行差分操作通常是为了消除均值的变化。从数学角度,差分可以写成:

    y t = y t – y t − 1 y_t = y_t – y_{t-1} yt=ytyt1

    其中 y t y_t yt 是t时刻的数值。相减数值之间的间隔数即为阶数。例如上述公式即为一阶差分。

    我们可以直接相减对序列差分,并绘制出对应线图:

    train['#Passengers_diff'] = train['#Passengers'] - train['#Passengers'].shift(1)
    train['#Passengers_diff'].dropna().plot()
    

    或者使用diff方法进行差分:

    train['#Passengers_diff'] = train['#Passengers'].diff(1)#一阶差分
    train['#Passengers_diff2'] = train['#Passengers_diff'].diff(1)#二阶差分
    

    在这里插入图片描述

    当数据存在季节性趋势时,我们可以利用季节性差分来消除季节性的不平稳因素。例如,星期一的观察值将与上星期一的观察值相减。从数学角度,它可以写成:

    y t = y t – y t − n y_t = y_t – y_{t-n} yt=ytytn

    n=7
    train['#Passengers_diff'] = train['#Passengers'] - train['#Passengers'].shift(n)
    

    变换

    变换用于对方差为非常数的序列进行平稳化。常用的变换方法包括幂变换、平方根变换和对数变换。

    train['#Passengers_log'] = np.log(train['#Passengers'])
    train['#Passengers_log_diff'] = train['#Passengers_log'] - train['#Passengers_log'].shift(1)
    train['#Passengers_log_diff'].dropna().plot()
    

    在这里插入图片描述

    在上面的变化过程中,我们首先对原始数据取对数,主要有两个用处:(1)将指数增长转为线性增长(2)可以平稳序列的方差。随后进行差分从而消除趋势的影响使序列平稳。

    分解

    对于有明显趋势或者周期性的时间序列二,我们也可以对其进行分解。分解需要用到statsmodels.tsa.seasonal.seasonal_decompose函数,可以将时间序列的数据分解为趋势(trend),季节性(seasonality)和残差(residual)三部分。

    from statsmodels.tsa.seasonal import seasonal_decompose
    decomposition = seasonal_decompose(train['#Passengers']).plot()#画出分解后时序图
    plt.show()
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    

    在这里插入图片描述

    这样趋势和季节性,还有残差值都被分解出来,之后我们就可以计算残差值的稳定性,从而得到一个平稳的时间序列

    参考文章:
    https://www.biaodianfu.com/arima.html

    在这里插入图片描述

    展开全文
  • 此代码是 [1] 中建议的 PMUCOS 方法的简单实现。 PMUCOS 解决了发现视频的所有周期片段并估计它们的周期的问题一种完全不受监督的方式。... [2] 中提出了一种更新的周期性检测方法。 [2] C.Panagiotakis和
  • 周期性时间序列的预测

    万次阅读 2018-08-17 18:00:00
    女主宣言AIOps 从立项到现在已经半年有余,从最开始的 LVS 异常检测,到如今的实时关联分析,智能运维已经渗透到我们日常运维中的许多场景,之后我们会将积累的经验分享出...

    女主宣言

    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
    展开全文
  • 把坑填上,用前后值的均值把缺失值补上,再看一眼: 填充好缺失值的序列 发现这份数据有这样几个特点,在模型设计和数据预处理的时候要考虑到: 1、这是一个周期性时间序列,数值有规律的以天为周期上下波动,图...
  • 时间序列的平稳性检验与随机性检验

    万次阅读 多人点赞 2019-04-05 14:42:56
        对于一个时间序列,在进行建模之前,首先需要进行平稳性检验和纯随机性检验,然后根据检验的结果再选择适合的模型。在讲解平稳性和随机性的定义之前,我们先介绍一下时间序列中常用的几个特征统计量。 2.1...
  • 异常点检测__ARIMA模型__时间序列中的4种常见异常:考虑 4 种特定的异常,分别是 innovational outlier (IO),additive outlier (AO),level shift (LS) 以及 temporary change (TC)。
  • 文章目录读取数据并绘制时序图绘制自相关图和偏自相关图进行季节差分模型拟合 读取数据并绘制时序图 install.packages("astsa") ...通过对相同数据的分析,可尝试S=12为季度周期。 进行季节差分 #计算月度均值 e1x2m
  • 上文我们已经知道了什么是时间序列的平稳性,也见到了一些平稳时间序列和非平稳的时间序列,那么当我们有一个新的时间序列数据时,怎么判断它是否是平稳的呢?时间序列平稳性检验方法,可分为三类:图形...
  • 用 Python 检验时间序列的平稳

    千次阅读 多人点赞 2020-10-27 08:40:00
    在做时间序列分析时,我们经常要对时间序列进行平稳性检验,而我们常用的软件是SPSS或SAS,但实际上python也可以用来做平稳性检验,而且效果也非常好,今天笔者就讲解一下如何用pyth...
  • 时间序列异常点检测

    万次阅读 2019-11-05 17:31:23
    时间序列异常检测--简单的语言介绍当前时间序列异常检测方法 在Statsbot,我们不断地回顾异常检测方法,并在此基础上改进我们的模型。 本文概述了目前最流行的时间序列异常检测算法及其优缺点。 这篇文章是写给...
  • 目录时间序列预测时间序列的平稳严平稳宽平稳三级目录 时间序列预测 按照时间的顺序把随机事件变化发展的过程记录下来就构成了一个时间序列。 x1,x2,x3,x4,...,xt 对时间序列进行观察、研究,找寻他变化发展的规律...
  • 时间序列预测/相似搜索/异常检测

    千次阅读 2020-10-19 19:43:55
    时间序列异常检测(一)—— 算法综述 时间序列异常检测(二)—— 基于KDD99数据集的实战 【论文分享】–多维时间序列异常检测 时序异常检测算法概览 AutoEncoder 是一种典型的无监督方法,可以将其扩展为 ...
  • P4J是一个python软件包,用于基于信息论目标函数对不规则采样和异方差时间序列进行周期检测。 P4J是为天文光曲线,恒星大小或通量的不规则采样时间序列开发的。 此程序包的核心是一个称为“周期图”的类,该类扫描...
  • matlab开发-应用于检测最新周期序列周期功率谱。周期功率谱及其在DNA序列潜在周期检测中的应用
  • 一文了解时间序列异常检测

    千次阅读 2020-12-08 20:19:36
    时间序列」是指某一个指标按照时间的统计或者观测而成的数列。比如,在运维的领域中,某主机每秒的CPU使用率、某业务每分钟的请求数量等,都可以形成一条时间序列;「异常检测」是指对反常的、和历史不同的行为...
  • 序列分解1、非季节性时间序列分解 移动平均MA(Moving Average)①SAM(Simple Moving Average) 简单移动平均,将时间序列上前n个数值做简单的算术平均。 SMAn=(x1+x2+…xn)/n②WMA(Weighted Moving Average) ...
  • Matlab时间序列分析

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

    万次阅读 多人点赞 2019-06-17 22:32:38
    SPSS(十九)SPSS之时间序列模型(图文+数据集) 时间序列是指将同一统计指标的数值按其发生的时间先后顺序排列而成的数列。正如人们常说,人生的出场顺序很重要,时间序列中隐藏着一些过去与未来的关系。时间序列...
  • 本文作者:张颖莹,文青松,孙亮,何凯,柯旻,蒋君伟近日,由阿里云计算平台和阿里云达摩院合作的时序多周期检测相关论文RobustPeriod: Robust Time-Frequency Mining for Multiple Periodicity Detection被...
  • 最近一直在接触时间序列,所以打算写一些有关时间序列的文章,预测部分会从规则开始、到传统模型、到机器学习、再到深度学习,此外也会介绍一些时间序列的基本概念,包括自相关、平稳、滞后、季节...
  • 浅谈LSTM对于周期时间序列数据的预测

    万次阅读 多人点赞 2019-03-06 12:33:40
    注:本文章主要针对的是长周期时间序列数据(10000-40000条为一个周期的数据)预测 产生训练和测试数据 我们需要做的是产生周期为20000条/周期的sin函数时间序列(用于训练) 以及周期为40000条/周期的sin...
  • 时间序列可预测度量

    千次阅读 2020-10-10 14:04:31
    时间序列可预测度量,讲解从序列长度到平稳,排序熵等指标
  • 主要介绍了python实现时间序列自相关图(acf)、偏自相关图(pacf)教程,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
  • 2.1 平稳性检验 一 、概率分布与特征统计量 Xt,t=1,2,⋅⋅⋅,tX_t,t=1,2,···,tXt​,t=1,2,⋅⋅⋅,t 在描述一个随机变量时是用 分布函数F(x)F(x)F(x) 特征统计量: 期望E(Xt)E(X_t)E(Xt​),方差D(Xt)D(X_t...
  • 时间序列--平稳介绍

    千次阅读 2021-04-08 20:51:26
    平稳时间序列中最重要的概念之一。 一个平稳的序列意味着它的均值、方差和协方差不随时间变化。 图一:均值是变化的(增长),整体是向上增长的趋势。在- 个平稳的序列里,它不应该有任何的变化趋势。 图二:...
  • 时间序列及异常检测综述(资料)

    千次阅读 多人点赞 2019-03-06 17:11:28
    文章目录1. 背景2. 时间序列预测方法3. ARIMA3.1 ARIMA模型预测的流程3.2 学习资料4. Prophet4.1 Prophet流程4.2 Prophet注意4.3 学习资料5. 其他时序方法6. 异常诊断相关方法7. 异常检测参考资料 1. 背景 时间...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 55,361
精华内容 22,144
关键字:

时间序列周期性检验