在众多时间序列模型中，SARIMA（seasonal autoregressive integrated moving average，季节性自回归积分滑动平均模型）能够有效处理时间序列中的季节性成分。但是在实际应用中，如何准确识别和提取这些季节性模式一直是一个挑战。
传统上，识别季节性模式往往依赖于数据的可视化分析。但是我们可以使用傅里叶变换以及周期图（Periodogram）这一强大工具，用一种更系统的方法来解决这个问题。
本文将详细介绍一些基础但重要的概念，这些方法对于每位研究时间序列的数据科学家都具有实用价值。
- 傅里叶变换的基本原理
- Python实现傅里叶变换
- 周期图分析方法
问题概述
我们以AEP（American Electric Power）能源消耗数据集为例（数据集使用CC0许可）：
importpandasaspd
importmatplotlib.pyplotasplt
df=pd.read_csv("data/AEP_hourly.csv", index_col=0)
df.index=pd.to_datetime(df.index)
df.sort_index(inplace=True)
fig, ax=plt.subplots(figsize=(20,4))
df.plot(ax=ax)
plt.tight_layout()
plt.show()
从数据的初步可视化中可以明显观察到季节性模式的存在，但要精确捕获所有这些模式并非易事。
传统的手动分析方法通常需要多个时间尺度的可视化，如下所示：
fig, ax=plt.subplots(3, 1, figsize=(20,9))
df_3y=df[(df.index>='2006–01–01') & (df.index<'2010–01–01')]
df_3M=df[(df.index>='2006–01–01') & (df.index<'2006–04–01')]
df_7d=df[(df.index>='2006–01–01') & (df.index<'2006–01–08')]
ax[0].set_title('AEP energy consumption 3Y')
df_3y[['AEP_MW']].groupby(pd.Grouper(freq='D')).sum().plot(ax=ax[0])
fordateindf_3y[[Trueifx% (24*365.25/2) ==0elseFalseforxinrange(len(df_3y))]].index.tolist():
ax[0].axvline(date, color='r', alpha=0.5)
ax[1].set_title('AEP energy consumption 3M')
df_3M[['AEP_MW']].plot(ax=ax[1])
fordateindf_3M[[Trueifx% (24*7) ==0elseFalseforxinrange(len(df_3M))]].index.tolist():
ax[1].axvline(date, color='r', alpha=0.5)
ax[2].set_title('AEP energy consumption 7D')
df_7d[['AEP_MW']].plot(ax=ax[2])
fordateindf_7d[[Trueifx%24==0elseFalseforxinrange(len(df_7d))]].index.tolist():
ax[2].axvline(date, color='r', alpha=0.5)
plt.tight_layout()
plt.show()
不同时间尺度下的AEP能源消耗模式
通过多尺度分析，我们可以观察到以下主要周期：
- 半年度周期（约180天）
- 周度周期（7天）
- 日度周期（24小时）
虽然对于能源消耗数据而言，这些季节性模式可以通过领域知识推断，但仅依靠人工检查存在以下局限性：
- 主观性：容易忽略不明显的模式
- 低效率：需要逐一检查不同时间窗口
- 扩展性差：难以应用于大规模数据分析
作为数据科学家，我们需要一个能够快速、准确识别时间序列中关键频率的工具。这正是傅里叶变换发挥作用的地方。
1、傅里叶变换的基本原理
傅里叶变换是一种将信号从时域转换到频域的数学工具。在时域中，我们观察数据随时间的变化；而在频域中，我们可以看到构成信号的各个频率及其相对重要性。
根据傅里叶理论，任何满足一定条件的函数f(x)都可以表示为不同频率、振幅和相位的正弦函数的叠加。换言之，每个时间序列都可以分解为一系列基本波形的组合。
其中：
- F(f)表示频域中的函数
- f(x)表示时域中的原始函数
- exp(-i2πf(x))是复指数函数，用作频率分析器
函数F(f)的值表明频率f在原始信号中的贡献程度。
考虑一个由三个不同频率（2Hz、3Hz和5Hz）的正弦波组成的复合信号：
应用傅里叶变换后，我们可以提取这些基本频率：
频域表示清晰地显示出信号中存在的三个基本频率（2Hz、3Hz和5Hz）。将信号分解为基本波形：
原始信号（蓝色）是三个基本波形（红色）的叠加。这种分解方法可以应用于任何时间序列，用于识别其主要频率成分。
2、Python中的傅里叶变换实现
现在我们将这个理论应用到之前的AEP能源消耗数据中。
Python的numpy.fft模块提供了计算离散信号傅里叶变换的工具。FFT（快速傅里叶变换）是一种高效的算法，用于将离散信号分解为频率成分：
fromnumpyimportfft
X=fft.fft(df['AEP_MW'])
N=len(X)
frequencies=fft.fftfreq(N, 1)
periods=1/frequencies
fft_magnitude=np.abs(X) /N
mask=frequencies>=0
# 绘制傅里叶变换结果
fig, ax=plt.subplots(figsize=(20, 3))
ax.step(periods[mask], fft_magnitude[mask]) # 仅绘制正频率
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
ax.set_title('AEP energy consumption - Frequency-Domain')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Magnitude')
plt.show()
频域可视化显示了某些频率具有较高的幅值，表明这些频率在原始信号中具有重要作用。
3、周期图分析
周期图是信号功率谱密度（Power Spectral Density, PSD）的频域表示。相比简单的傅里叶变换，周期图能更好地量化各频率分量的强度，并且可以有效降低次要频率的噪声影响。
周期图的数学定义为：
其中：
- P(f)是频率f处的功率谱密度
- X(f)是信号的傅里叶变换
- N是样本数量
Python实现如下：
power_spectrum=np.abs(X)**2/N # 计算每个频率的功率
fig, ax=plt.subplots(figsize=(20, 3))
ax.step(periods[mask], power_spectrum[mask])
ax.set_title('AEP energy consumption Periodogram')
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.show()
从周期图分析中，我们可以清晰地识别出以下主要频率：
- 24Hz：对应24小时周期
- 4,380Hz：对应6个月周期
- 168Hz：对应每星期
除此之外，还发现了三个次要周期：
- 12Hz周期
- 84Hz周期（对应半周）
- 8,760Hz周期（对应年度周期）
我们也可以使用scipy.signal模块中的periodogram函数获得相同的结果：
fromscipy.signalimportperiodogram
frequencies, power_spectrum=periodogram(df['AEP_MW'], return_onesided=False)
periods=1/frequencies
fig, ax=plt.subplots(figsize=(20, 3))
ax.step(periods, power_spectrum)
ax.set_title('Periodogram')
ax.set_xscale('log')
ax.xaxis.set_major_formatter('{x:,.0f}')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.show()
总结
在时间序列分析中，准确识别季节性模式是一个关键任务。本文介绍的周期图分析方法为我们提供了一个系统、高效的工具，可以准确识别时间序列中的各种周期性成分。这种方法不仅克服了传统视觉分析的局限性，还能发现可能被忽略的次要周期模式，为时间序列分析提供了更全面的视角。
https://avoid.overfit.cn/post/e459dbb3b4f8461c9224eefee2672d77
作者：Lorenzo Mezzini
