谱密度估计与周期分析：用傅里叶把隐藏的周期从噪声里拎出来 • halo 的技术博客

先说结论：谱密度估计能把时域里完全看不见的周期从噪声里拎出来——合成实验里信噪比只有 0.19（噪声方差是信号的 5 倍），肉眼看序列就是随机抖动，Welch 谱上 20 日和 63 日两个周期的峰清清楚楚。但这门技术有三个大坑，每一个都足以让你把噪声当成周期去交易：一是原始周期图的方差不随样本量收敛，毛刺永远满屏；二是随机游走的谱天然带 1/f² 斜坡，最大伪峰能达到中位数的 375 倍，不差分就做谱分析必然”发现”假周期；三是频率分辨率——本次实验中 Welch 估出 19.69 日、真实周期 20 日，误差只有 1.5%，外推 512 日后相位漂移 144°，择时方向准确率掉到 0.393，比抛硬币还差。把频率精调到 20.06 日后，准确率回升到 0.584，日均收益从 -0.34 变成 +0.32。

谱密度估计：时域看不见的周期，频域一眼看穿

一、为什么频域能看见时域看不见的东西#

时域里的一条序列，是所有频率成分叠加后的总和。当一个振幅 0.6 的 20 日正弦被方差 1.4 的 AR(1) 噪声淹没时，任何一天的观测值里周期成分只占很小一部分，肉眼和大多数时域统计量（均值、自相关的单点值）都很难分辨。

频域的杀手锏是能量集中：周期信号的能量全部堆在一个频点上，而噪声的能量摊在整个频率轴上。1024 个样本里，20 日正弦的全部能量集中到 1/20 = 0.05 这一个频点，噪声则被摊薄到 512 个频点——局部信噪比被放大了两个数量级。这就是”时域看不见、频域一眼看穿”的数学原理。

功率谱密度（PSD） $S(f)$ 的定义是自协方差函数的傅里叶变换（Wiener–Khinchin 定理）：

$S(f) = \sum_{k=-\infty}^{\infty} \gamma(k) e^{-2\pi i f k}$

它回答的问题是：序列的总方差在各个频率上是怎么分配的。谱上的峰 = 该频率贡献了超额方差 = 可能存在周期。

二、实验设置：把周期埋进噪声里#

合成一条 1024 日的序列，模拟”某个日度收益溢价”：

import numpy as np
from scipy import signal as sps

rng = np.random.default_rng(42)
n = 1024
t = np.arange(n)
P1, P2 = 20, 63    # 隐藏周期：约一个月、一个季度的交易日数
sig_clean = 0.6*np.sin(2*np.pi*t/P1) + 0.4*np.sin(2*np.pi*t/P2 + 0.7)

# AR(1) 噪声，phi=0.5，方差约为信号的 5 倍
eps = np.zeros(n)
e = rng.normal(0, 1.0, n)
for i in range(1, n):
    eps[i] = 0.5*eps[i-1] + e[i]
x = sig_clean + eps
# 信号方差 0.261，噪声方差 1.377，SNR = 0.19

python

选 AR(1) 而不是白噪声是刻意的：金融序列的噪声几乎都有正自相关（波动聚集的残留、微观结构效应），AR(1) 噪声的谱在低频端天然抬高——这正是后面显著性检验必须用”红噪声基线”而不是”白噪声基线”的原因。

三、原始周期图 vs Welch：方差不收敛是周期图的原罪#

最朴素的谱估计是周期图（periodogram）：对整条序列做一次 FFT，取模平方。

freqs_pg, pxx_pg = sps.periodogram(x, fs=1.0, detrend='constant')

python

周期图是渐近无偏的，但有个致命缺陷：它的方差不随样本量增加而收敛。每个频点的周期图值渐近服从 $S(f) \cdot \chi^2_2 / 2$ 分布——自由度永远是 2，标准差永远和均值同阶。样本从 1024 加到 10 万，毛刺一根不少。

Welch 方法用”分段 + 加窗 + 平均”换取方差收敛：

freqs_w, pxx_w = sps.welch(x, fs=1.0, nperseg=256,
                            noverlap=128, detrend='constant')

python

1024 个样本切成 256 长的段、50% 重叠，共 7 段，平均后方差降为原来的约 1/7。代价是频率分辨率从 1/1024 降到 1/256——这个代价后面会要了择时策略的命。

原始周期图毛刺满屏 vs Welch 平均后双峰清晰

左图原始周期图上真实峰完全埋在毛刺里；右图 Welch 谱上 0.05（20 日）和 0.016（63 日）两个峰清晰可辨。峰值检测报出的周期是 19.69 日和 64 日——注意都不精确等于真值，这不是 bug，是 Welch 的频率格点只有 1/256 = 0.0039 的间距，真实频率落在两个格点之间时只能就近取整。记住这个 19.69，第五节它会闯祸。

四、显著性：峰高不算数，越过红噪声包络才算数#

谱上有峰 ≠ 有周期。噪声自己也会随机堆出峰来，而且有色噪声堆出的峰还有系统性偏好——AR(1) 噪声的理论谱是：

$S_{AR(1)}(f) = \frac{\sigma^2}{1 + \phi^2 - 2\phi\cos(2\pi f)}$

$\phi > 0$ 时低频端单调抬高。如果拿白噪声的平坦基线做检验，低频区的一切都”显著”——这是气候学和金融里最常见的假周期来源。

正确做法：用数据拟合 AR(1) 参数，蒙特卡洛模拟红噪声的谱分布，取分位数做包络。

phi_hat = np.corrcoef(x[:-1], x[1:])[0, 1]
s2 = np.var(x) * (1 - phi_hat**2)

nsim = 500
sim = np.zeros((nsim, len(freqs_w)))
for k in range(nsim):
    ee = rng.normal(0, np.sqrt(s2), n)
    xx = np.zeros(n)
    for i in range(1, n):
        xx[i] = phi_hat*xx[i-1] + ee[i]
    _, sim[k] = sps.welch(xx, fs=1.0, nperseg=256, noverlap=128,
                          detrend='constant')

thr95 = np.quantile(sim, 0.95, axis=0)
thr99 = np.quantile(sim, 0.99, axis=0)

python

显著性检验：只有越过红噪声包络的峰才算数

结果很有信息量：20 日周期越过了 99% 包络，63 日周期只越过 95% 包络。63 日成分振幅更小（0.4 vs 0.6），又落在 AR(1) 噪声天然抬高的低频区，检验功效被双重削弱。实务含义：低频周期（月度、季度、年度效应）天然更难与红噪声区分，声称发现了低频周期的研究，要求的证据强度应该更高，而市面上大多数”季节性效应”研究恰恰在这里最草率。

顺带一提：这里用逐频点分位数做包络，存在多重比较问题——512 个频点每个都给 5% 的假阳性率，全局至少一个假阳性的概率远超 5%。更严格的做法是用”模拟谱的全局最大值分布”做阈值（类似 White’s Reality Check 的思路），代价是功效下降。本次实验里 20 日峰高出包络一个数量级，两种标准下都显著。

五、反面教材：随机游走的谱里全是假周期#

把谱分析直接用在价格（而不是收益率）上，是这个领域第一大事故现场。随机游走 $p_t = p_{t-1} + \varepsilon_t$ 的理论谱是 $S(f) \propto 1/f^2$ ——低频端能量爆炸。

实测：对一条纯随机游走做 Welch 谱，最大伪峰是谱中位数的 375 倍。任何”找峰”算法都会兴奋地报告一个低频”周期”，而这条序列里连一丝周期都没有。

随机游走的 1/f² 假峰 vs 差分后的平坦白噪声谱

差分（取收益率）之后，谱立刻变平——白噪声，无周期，干干净净。规则很简单但值得写进 checklist：

永远对收益率（或其他平稳序列）做谱分析，不要对价格做；
即使是收益率，也先 detrend 去掉均值和线性趋势；
看到低频大峰，第一反应是怀疑非平稳残留，不是庆祝发现周期。

顺便拆一个流传很广的说法：“股价有 XX 年康波周期/太阳黑子周期”。这类研究大多直接对价格指数（甚至对数价格都不取）做谱分析，低频伪峰要多少有多少——1/f² 斜坡上随便切一段都能讲出周期故事。

六、择时实验：频率差 0.4%，收益差出天壤#

检出了显著的 20 日周期，能不能用？前 512 日做训练（估周期 + 拟合相位），后 512 日做样本外择时：预测的周期成分为正就做多，为负就做空。

第一版：直接用 Welch 峰值频率 19.69 日。 结果惨烈——方向准确率 0.393（比抛硬币还差），日均收益 -0.34。

问题出在相位漂移。频率误差看起来只有 $|1/19.69 - 1/20| = 0.08\%$ 每日，但相位误差是累积的：

$\Delta\phi = T \cdot \left|\frac{1}{\hat{P}} - \frac{1}{P}\right| \cdot 360° = 512 \times 0.00079 \times 360° \approx 144°$

外推 512 日后相位漂了 144°——接近反相。前半段预测还大致对齐，后半段完全反着来，做多的日子恰好是周期低谷。方向准确率低于 0.5 不是”没信号”，是信号被反着用了。

第二版：把周期作为自由参数精调。 用 Welch 峰值做初值，对训练段做非线性最小二乘：

from scipy.optimize import curve_fit

def sin_free(tt, A, P, ph, c):
    return A*np.sin(2*np.pi*tt/P + ph) + c

popt, _ = curve_fit(sin_free, np.arange(512), train,
                    p0=[1, 19.69, 0, 0])   # Welch 峰值做初值
# 拟合结果：P = 20.056（真值 20），振幅 0.557

python

正弦拟合的频率精度不受 FFT 格点限制（本质是把 512 个样本的信息全部用于估一个连续参数），拟合出 20.056 日，误差 0.28%。样本外：方向准确率 0.584，日均 +0.32，t 统计量 5.7。

粗频率相位漂移 vs 精调频率择时的天壤之别

同一个”发现了 20 日周期”的结论，频率估计差 0.4%（19.69 vs 20.06），策略从日均 -0.34 到 +0.32。周期交易的生死线不在”有没有周期”，而在”频率估得够不够准、能撑多久不漂”。 实务上的对策有三个：用正弦拟合或 MUSIC 类高分辨率方法精调频率；缩短外推窗口、滚动重估相位；或者干脆放弃外推，改用带通滤波器实时提取周期成分（代价是滤波器的群延迟）。

七、诚实的部分：真实市场里这套东西的胜率#

上面是合成数据，周期是真实存在且频率恒定的。真实市场残酷得多：

第一，真实周期几乎都是时变的。 就算某段时间存在 20 日左右的周期（期权到期、月度调仓资金流），它的频率和振幅都会漂移。固定频率的正弦外推在合成数据上撑 512 日，在真实数据上能撑 50 日就不错了。对策是滚动窗口重估，但窗口短了频率精度又不够——第六节说了，频率精度就是生死线。这是个真实的两难，不要假装它不存在。

第二，可交易的周期会被交易掉。 稳定、显著、频率恒定的价格周期是最容易被套利的模式之一。谱分析在金融里更现实的用途不是直接择时，而是：检验策略残差里有没有周期性结构（有则说明模型漏了日历效应）、辨认数据管道的伪周期（比如每周五的结算价异常会在谱上打出 1/5 频率的峰）、以及给波动率建模提供周期性成分（日内 U 型、周内效应）。

第三，样本外这次实验赢了，部分靠运气。 AR(1) 噪声的样本外段恰好没有出现长的负漂移。把随机种子换 20 个重跑，精调版的方向准确率在 0.55–0.60 之间波动，粗频率版在 0.35–0.50 之间——精调版稳定更好，但”日均 +0.32”这个具体数字不要外推。

八、Checklist：谱分析用于量化研究的最小纪律#

只对平稳序列做谱分析——收益率可以，价格不行；先 detrend；
用 Welch（或 multitaper）而不是原始周期图——方差不收敛的估计量没有实用价值；
显著性对着红噪声基线检验——拟合 AR(1)，蒙特卡洛做包络，白噪声基线会在低频区大量误报；
低频峰要求更高的证据标准——那里是红噪声和非平稳残留的重灾区；
要交易就必须精调频率——FFT 格点精度不够外推用，正弦拟合或滚动重估相位；算一下 $T \cdot \Delta f \cdot 360°$ ，相位漂移超过 90° 的外推等于自杀；
对”发现周期”保持默认怀疑——先排除数据管道伪周期、日历效应、非平稳残留，剩下的才配讨论。

谱密度估计是把时间序列换一个坐标系观察的技术。它不生产 alpha，它生产的是”这条序列的方差花在了哪些频率上”这个事实。事实本身没有方向，用对了是显微镜，用错了是万花筒——每转一下都能看见新的”周期”。