时间序列 Python

背景

任何事物在两个不同时刻都不可能保持完全相同的状态，但很多变化往往存在着一定的规律，例如 24 小时日出日落，潮起潮落，这些现象通常称为「周期」。
周期性，指时间序列中呈现出来的围绕长期趋势的一种波浪形或振荡式变动。准确提取周期信息，不仅能反映当前数据的规律，应用于相关场景，还可以预测未来数据变化趋势。
2021-07-31-08-27-14-359295.png
时间序列示例
一般而言，时间序列周期性分为三种：

  • 「符号性周期」，例如序列 fbcnfkgbfops中 f 周期为 4；
  • 「部分周期性」，例如序列 ansdcdmncdcacdascdmc 中 cd 周期为 4；
  • 「分段周期性」，例如上面给定的时间序列即为分段周期性；

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

傅里叶变换

傅里叶变换是一种将时域、空域数据转化为频域数据的方法，任何波形（时域）都可以看做是不同振幅、不同相位正弦波的叠加（频域）👇🏻 此处放上经典图！
2021-07-31-08-27-14-477296.jpeg
傅里叶变换
对于一条具备周期性的时间序列，它本身就很接近正弦波，所以它包含一个显著的正弦波，周期就是该正弦波的周期，而这个正弦波可以通过傅里叶变换找到，它将时序数据展开成三角函数的线性组合，得到每个展开项的系数，就是傅里叶系数。傅里叶系数越大，表明它所对应的正弦波的周期就越有可能是这份数据的周期。

自相关系数

自相关系数（Autocorrelation Function）度量的是同一事件不同时间的相关程度，不同相位差（lag）序列间的自相关系数可以用 Pearson 相关系数计算。其数学表达如下：
其中 表示相位 的数据延迟 「lag operator」。
2021-07-31-08-27-14-601309.jpeg
自相关系数
当序列存在周期性时，遍历足够多的相位差，一定可以找到至少一个足够大的自相关系数，而它对应的相位差就是周期。所以对于检测时序周期来说，只需找到两个自相关系数达到一定阈值的子序列，它们起始时间的差值就是需要的周期。

实例说明

为了保证结果的可靠性，可以将傅里叶分析和自相关系数结合起来判断周期性。主要思路是：先通过傅里叶变换找到可能的周期，再用自相关系数做排除，从而得到最可能的周期。
给定一份周期性数据，时间间隔为 5 min。从这份数据中可以看出数据大体上具有周期为 1 day。
2021-07-31-08-27-14-861759.png
示例数据
下面使用傅里叶变换估计周期，代码如下所示

  1. from scipy.fftpack import fft, fftfreq
  3. fft_series = fft(data["value"].values)
  4. power = np.abs(fft_series)
  5. sample_freq = fftfreq(fft_series.size)
  7. pos_mask = np.where(sample_freq > 0)
  8. freqs = sample_freq[pos_mask]
  9. powers = power[pos_mask]
  11. top_k_seasons = 3
  12. # top K=3 index
  13. top_k_idxs = np.argpartition(powers, -top_k_seasons)[-top_k_seasons:]
  14. top_k_power = powers[top_k_idxs]
  15. fft_periods = (1 / freqs[top_k_idxs]).astype(int)
  17. print(f"top_k_power: {top_k_power}")
  18. 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 已经非常接近真实值。
现在来计算自相关系数，代码如下所示：

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

对应的输出如下：

  1. lag: 72 fft acf: 0.07405431832776994
  2. lag: 278 fft acf: 0.7834457453491087
  3. lag: 292 fft acf: 0.8259822269757922
  4. lag: 144 expected acf: -0.5942986094704665
  5. lag: 288 expected acf: 0.8410792774898174
  6. lag: 2016 expected acf: 0.5936030431473589

通过自相关系数来得到显著分数最大值对应的周期，得出的结果为 292；
此处实验补充了预设的三个周期值：12 hour、1 day、7 day，发现算出来还是周期 288 对应的相关分数最大，但是傅里叶变换没有估计出周期值。
综上：算出来的还不如预设的值呢！直接根据先验知识「预设周期」然后计算自相关系数就行了！