时间序列丨周期性检测方法及其 Python 实践

背景

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

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

一般而言，时间序列周期性分为三种：

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

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

傅里叶变换

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

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

自相关系数

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

$\text { autocorr }(X, t)= corr[X, lag(X,t)]=\frac{\operatorname{cov}[X, \operatorname{lag}(X, t)]}{\operatorname{std}[X] \operatorname{std}[\operatorname{lag}(X, t)]}=\frac{\operatorname{cov}[X, \operatorname{lag}(X, t)]}{\operatorname{var}(X)}= \\ =\frac{\sum_{i=1}^{N}\left[X_{i}-\operatorname{mean}(X)\right]\left[X_{i-t}-\operatorname{mean}(X)\right]}{\sum_{i=1}^{N}\left[X_{i}-\operatorname{mean}(X)\right]^{2}}$

其中 $lag(X_i,t) = X_{i-t}$ 表示相位 $t$ 的数据延迟 lag operator。

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

Python 实践

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

给定一份周期性数据，时间间隔为 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