在眾多時(shí)間序列模型中，SARIMA（seasonal autoregressive integrated moving average，季節(jié)性自回歸積分滑動(dòng)平均模型）能夠有效處理時(shí)間序列中的季節(jié)性成分。但是在實(shí)際應(yīng)用中，如何準(zhǔn)確識(shí)別和提取這些季節(jié)性模式一直是一個(gè)挑戰(zhàn)。

傳統(tǒng)上，識(shí)別季節(jié)性模式往往依賴(lài)于數(shù)據(jù)的可視化分析。但是我們可以使用傅里葉變換以及周期圖（Periodogram）這一強(qiáng)大工具，用一種更系統(tǒng)的方法來(lái)解決這個(gè)問(wèn)題。

本文將詳細(xì)介紹一些基礎(chǔ)但重要的概念，這些方法對(duì)于每位研究時(shí)間序列的數(shù)據(jù)科學(xué)家都具有實(shí)用價(jià)值。

1. 傅里葉變換的基本原理
2. Python實(shí)現(xiàn)傅里葉變換
3. 周期圖分析方法

問(wèn)題概述

我們以AEP（American Electric Power）能源消耗數(shù)據(jù)集為例（數(shù)據(jù)集使用CC0許可）：

import pandas as pd
 import matplotlib.pyplot as plt
 
 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()

從數(shù)據(jù)的初步可視化中可以明顯觀察到季節(jié)性模式的存在，但要精確捕獲所有這些模式并非易事。

傳統(tǒng)的手動(dòng)分析方法通常需要多個(gè)時(shí)間尺度的可視化，如下所示：

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])
 for date in df_3y[[True if x % (24 * 365.25 / 2) == 0 else False for x in range(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])
 for date in df_3M[[True if x % (24 * 7) == 0 else False for x in range(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])
 for date in df_7d[[True if x % 24 == 0 else False for x in range(len(df_7d))]].index.tolist():
     ax[2].axvline(date, color = 'r', alpha = 0.5)
 
 plt.tight_layout()
 plt.show()

不同時(shí)間尺度下的AEP能源消耗模式

通過(guò)多尺度分析，我們可以觀察到以下主要周期：

• 半年度周期（約180天）
• 周度周期（7天）
• 日度周期（24小時(shí)）

雖然對(duì)于能源消耗數(shù)據(jù)而言，這些季節(jié)性模式可以通過(guò)領(lǐng)域知識(shí)推斷，但僅依靠人工檢查存在以下局限性：

• 主觀性：容易忽略不明顯的模式
• 低效率：需要逐一檢查不同時(shí)間窗口
• 擴(kuò)展性差：難以應(yīng)用于大規(guī)模數(shù)據(jù)分析

作為數(shù)據(jù)科學(xué)家，我們需要一個(gè)能夠快速、準(zhǔn)確識(shí)別時(shí)間序列中關(guān)鍵頻率的工具。這正是傅里葉變換發(fā)揮作用的地方。

1、傅里葉變換的基本原理

傅里葉變換是一種將信號(hào)從時(shí)域轉(zhuǎn)換到頻域的數(shù)學(xué)工具。在時(shí)域中，我們觀察數(shù)據(jù)隨時(shí)間的變化；而在頻域中，我們可以看到構(gòu)成信號(hào)的各個(gè)頻率及其相對(duì)重要性。

根據(jù)傅里葉理論，任何滿(mǎn)足一定條件的函數(shù)f(x)都可以表示為不同頻率、振幅和相位的正弦函數(shù)的疊加。換言之，每個(gè)時(shí)間序列都可以分解為一系列基本波形的組合。

其中：

• F(f)表示頻域中的函數(shù)
• f(x)表示時(shí)域中的原始函數(shù)
• exp(-i2πf(x))是復(fù)指數(shù)函數(shù)，用作頻率分析器

函數(shù)F(f)的值表明頻率f在原始信號(hào)中的貢獻(xiàn)程度。

考慮一個(gè)由三個(gè)不同頻率（2Hz、3Hz和5Hz）的正弦波組成的復(fù)合信號(hào)：

應(yīng)用傅里葉變換后，我們可以提取這些基本頻率：

頻域表示清晰地顯示出信號(hào)中存在的三個(gè)基本頻率（2Hz、3Hz和5Hz）。將信號(hào)分解為基本波形：

原始信號(hào)（藍(lán)色）是三個(gè)基本波形（紅色）的疊加。這種分解方法可以應(yīng)用于任何時(shí)間序列，用于識(shí)別其主要頻率成分。

2、Python中的傅里葉變換實(shí)現(xiàn)

現(xiàn)在我們將這個(gè)理論應(yīng)用到之前的AEP能源消耗數(shù)據(jù)中。

Python的numpy.fft模塊提供了計(jì)算離散信號(hào)傅里葉變換的工具。FFT（快速傅里葉變換）是一種高效的算法，用于將離散信號(hào)分解為頻率成分：

from numpy import fft
 
 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
 
 # 繪制傅里葉變換結(jié)果
 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()

頻域可視化顯示了某些頻率具有較高的幅值，表明這些頻率在原始信號(hào)中具有重要作用。

3、周期圖分析

周期圖是信號(hào)功率譜密度（Power Spectral Density, PSD）的頻域表示。相比簡(jiǎn)單的傅里葉變換，周期圖能更好地量化各頻率分量的強(qiáng)度，并且可以有效降低次要頻率的噪聲影響。

周期圖的數(shù)學(xué)定義為：

其中：

• P(f)是頻率f處的功率譜密度
• X(f)是信號(hào)的傅里葉變換
• N是樣本數(shù)量

Python實(shí)現(xiàn)如下：

power_spectrum = np.abs(X)**2 / N  # 計(jì)算每個(gè)頻率的功率
 
 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()

從周期圖分析中，我們可以清晰地識(shí)別出以下主要頻率：

• 24Hz：對(duì)應(yīng)24小時(shí)周期
• 4,380Hz：對(duì)應(yīng)6個(gè)月周期
• 168Hz：對(duì)應(yīng)每星期

除此之外，還發(fā)現(xiàn)了三個(gè)次要周期：

• 12Hz周期
• 84Hz周期（對(duì)應(yīng)半周）
• 8,760Hz周期（對(duì)應(yīng)年度周期）

我們也可以使用scipy.signal模塊中的periodogram函數(shù)獲得相同的結(jié)果：

from scipy.signal import periodogram
 
 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()

總結(jié)

在時(shí)間序列分析中，準(zhǔn)確識(shí)別季節(jié)性模式是一個(gè)關(guān)鍵任務(wù)。本文介紹的周期圖分析方法為我們提供了一個(gè)系統(tǒng)、高效的工具，可以準(zhǔn)確識(shí)別時(shí)間序列中的各種周期性成分。這種方法不僅克服了傳統(tǒng)視覺(jué)分析的局限性，還能發(fā)現(xiàn)可能被忽略的次要周期模式，為時(shí)間序列分析提供了更全面的視角。