奇异谱分析 SSA 周期提取：把价格嵌入成轨迹矩阵做 SVD，看不见的周期自己浮现 • halo 的技术博客

你有一段价格，怀疑里面藏着「约 37 天的短周期」和「约 120 天的中周期」，想把它们单独拎出来做轮动或者确认趋势里有没有周期性扰动。傅里叶变换能做，但它要整段周期稳定（非平稳的价格里基不合适）；小波能做，但要选母小波；VMD 能做，但要拍 K 和 α。

结论先放这：奇异谱分析 SSA（Vautard & Ghil 1989；Golyandina & Zhigljavsky 2013）是这几类里「最不预设结构」的一个。 它的核心想法出奇简单——把一维时间序列延迟嵌入成一个矩阵，对这个矩阵做 SVD，SVD 出来的本征分量天然就是「趋势 / 周期 / 噪声」的候选，你只要决定怎么分组。 不挑基、不拍 K、不需要周期先验。

在我们合成的价格（慢饱和趋势 + 周期 120 + 周期 37 + 噪声）上，SSA 趋势相对真实趋势的 RMSE≈2.05，优于 HP 滤波的 2.34；前 6 个本征分量吃掉了 100% 的方差；隐藏的周期 37 用一对本征分量还原后，与真实正弦的相关高达 0.934。诚实地说：SSA 也有真问题——嵌入窗口 L 要拍、分量怎么配对成周期靠经验、弱周期会被噪声淹没、头尾 L−1 个点补不回来，后文五类陷阱要一一拆穿。但它确实是「让看不见的周期自己从 SVD 里浮现」最干净的工具。

SSA 把价格嵌入成轨迹矩阵、做 SVD，最平滑的本征分量就是趋势

1. 为什么傅里叶/小波不够「自适应」#

傅里叶：把信号投影到整段正弦基。非平稳价格（趋势 + 突变 + 周期漂移）里，一个频率会「漏」到很多基上，频谱糊成一片。
小波：要选母小波和尺度，选错基就全错，而且金融序列的多尺度结构很难先验匹配。
VMD：要拍模态数 K 和带宽 α（上篇聊过）。

SSA 反其道而行：不引入任何外部基。它只用序列自身的延迟结构构造矩阵，然后让 SVD 自己长出「基」来——这些基是数据驱动的，且天然按「解释方差」排序。

2. 第一步：延迟嵌入（Embedding）#

给定长度为 N 的序列 $x_1,\dots,x_N$ ，选一个窗口长度 $L$ （通常取 $1 < L < N/2$ ），构造轨迹矩阵 $X$ ：

X = \begin{bmatrix} x_1 & x_2 & \cdots & x_{N-L+1} \\ x_2 & x_3 & \cdots & x_{N-L+2} \\ \vdots & \vdots & \ddots & \vdots \\ x_L & x_{L+1} & \cdots & x_N \end{bmatrix}

形状是 $L \times K$ ，其中 $K = N - L + 1$ 。第 j 列就是「从第 j 个样本开始的、长度为 L 的延迟切片」。这个矩阵有个重要性质：它是「轨迹矩阵」（Hankel 型）——每条反对角线上的元素相等。后面还原时要利用这一点。

import numpy as np
from scipy.linalg import svd

def embed(x, L):
    N = len(x)
    K = N - L + 1
    X = np.column_stack([x[i:i + K] for i in range(L)])   # 形状 (L, K)
    return X

# 例: N=600, L=120 -> X 形状 (120, 481)

python

3. 第二步：SVD 把方差拆成本征分量#

对轨迹矩阵做 SVD：

X = \sum_{i=1}^{d} \sqrt{\lambda_i}\, U_i V_i^T, \qquad d = \mathrm{rank}(X)

$\sqrt{\lambda_i}$ 是第 i 个奇异值（ $\lambda_i$ 是方差的代理，按从大到小排）。
$U_i$ 是左奇异向量（长度 L，描述「分量在时间窗内的形状」）。
$V_i$ 是右奇异向量（长度 K，描述「分量沿原序列的演化」）。

每一对 $(U_i, V_i)$ 加上奇异值，就是一个本征三元组（eigentriple）。前几个本征三元组解释绝大部分方差。

SSA 谱（scree）：少数本征分量吃掉大部分方差

在我们的合成价格上（L=120），前 4 个奇异值约为 28997 / 778 / 724 / 366，前 6 个本征分量就吃掉了 100% 的方差——说明信号结构高度低秩（一个趋势 + 两个周期 + 少量噪声，正好对应少数几个强分量）。

4. 第三步：对角平均还原（Reconstruction）#

SVD 把矩阵拆了，但我们要的是一维序列。对每个（或每组）本征三元组，先构造秩-1 矩阵 $X_i = \sqrt{\lambda_i} U_i V_i^T$ ，再用**对角平均（Hankelization）**把它映射回一维：

\tilde{x}_t = \frac{1}{\#\text{反对角}} \sum_{i+j=t+1} (X_i)_{ij}

也就是对矩阵的每条反对角线取平均。这一步强制还原出的序列是「干净的」（Hankel 结构），消除嵌入时引入的冗余。

def diagonal_average(D):
    """对角平均：把 (L,K) 矩阵还原成 1-D 序列（Hankelization）。"""
    L, K = D.shape
    N = L + K - 1
    rec = np.zeros(N)
    counts = np.zeros(N)
    for i in range(L):
        for j in range(K):
            rec[i + j] += D[i, j]
            counts[i + j] += 1
    return rec / counts

def ssa_components(x, L, groups):
    """groups: 每组要合并的本征分量下标列表；返回每组还原后的 1-D 序列。"""
    X = embed(x, L)
    U, s, VT = svd(X, full_matrices=False)
    out = []
    for g in groups:
        D = np.zeros_like(X)
        for k in g:
            D += s[k] * np.outer(U[:, k], VT[k])
        out.append(diagonal_average(D))
    return out

python

5. 怎么分组：趋势 / 周期 / 噪声#

SVD 出来的本征分量按方差排序，但**怎么把它们映射到「趋势 / 周期 / 噪声」**是 SSA 的关键手艺：

趋势：最平滑的分量。判断方法——看左奇异向量 $U_i$ 的「曲率」（二阶差分平方和），曲率最小的就是最平滑、对应趋势。
周期：成对出现的分量。一个正弦振荡会拆成两个相邻的本征三元组，它们的左奇异向量是「正弦 / 余弦」的正交对（ quadrature）。把这对合并，就还原出一个干净正弦。
噪声：剩余的高频、低方差分量，丢弃即可。

def smoothest_indices(U, s, top=2):
    """返回最平滑的 top 个本征分量下标（趋势候选）。"""
    smooth = []
    for k in range(len(s)):
        u = U[:, k]
        curv = np.sum(np.diff(u, 2) ** 2)        # 二阶差分平方和 = 曲率
        smooth.append((curv, k))
    smooth.sort()
    return [idx for _, idx in smooth[:top]]

# 趋势 = 两个最平滑分量之和
trend_idx = smoothest_indices(U, s, top=2)
trend = ssa_components(price, L=120, groups=[trend_idx])[0]

python

6. 看结果：趋势提取与基准对比#

SSA 趋势 vs HP 滤波：无端点漂移、低频干净

相对真实趋势的 RMSE：

方法	趋势 RMSE
SSA（L=120）	2.05
HP 滤波	2.34

SSA 略优于 HP，且有个 HP 没有的优点：SSA 的趋势是「从数据自身最平滑的本征结构里长出来的」，不依赖惩罚项 λ，端点也不会被拽向最后一点（HP 的端点漂移是老问题）。代价是你要拍窗口 L。

7. SSA 的隐藏能力：把看不见的周期还原出来#

这是 SSA 最迷人的地方。我们注入的周期 37，在合成价格里被噪声盖着、肉眼难辨。但 SVD 把它的能量集中到了一对本征分量上（下标 4 和 5）。把这对合并、还原，就得到一个干净正弦：

SSA 把看不见的周期还原出来 + 趋势对噪声的鲁棒性

左图：还原出的振荡（分量 4+5）与真实周期 37 正弦几乎重合，相关系数高达 0.934——SSA 真的把「看不见的周期」从噪声里抠了出来。右图：噪声从 0.02 加到 0.20，趋势 RMSE 只从 2.08 温和升到 2.21，对噪声鲁棒。

# 自动找「最像某个周期」的相邻分量对
best_corr, best_pair = -1, None
for k in range(0, 14, 2):                 # 相邻对
    comp = ssa_components(price, L=120, groups=[[k, k + 1]])[0]
    comp = comp[:len(true_cycle37)] - comp[:len(true_cycle37)].mean()
    corr = np.corrcoef(comp, true_cycle37)[0, 1]
    if corr > best_corr:
        best_corr, best_pair = corr, (k, k + 1)
# 结果: 最佳对 = (4, 5), 与真实周期37相关 = 0.934

python

8. 五类真实陷阱（不拆穿就是自欺）#

陷阱 1：窗口 L 主观，且决定了一切。 L 太小，轨迹矩阵太「瘦」，周期长于 L 的分量根本进不去；L 太大，计算量爆炸、且把多个周期搅进同一分量。经验法则： $L \approx N/4$ 起步，且 L 要大于你想找的最长周期。没有自动最优 L。

陷阱 2：分量配对（哪些 i 合成一个周期）靠经验。 一个真实周期拆成两个相邻本征三元组，但噪声也会产生相邻小分量。配对错了（把噪声当周期、或拆错对）就还原出假振荡。实务上靠「看左奇异向量的正弦形状 + 验证还原序列的周期」双保险。

陷阱 3：弱周期会被噪声淹没。 当某个周期振幅远小于噪声，它的能量分散到很多小本征分量上，SVD 里不再突出，你还原不出来。这不是 bug，是信噪比的物理下限——别强行宣称「找到了」一个相关 0.3 的弱周期。

陷阱 4：端点缺失，头尾 L−1 个点补不齐。 嵌入用长度为 L 的窗，序列头尾各有 L−1 个位置无法被完整窗覆盖，还原后这些点要么缺失要么方差小。对「最新一根 K 线」做实时周期检测时要留缓冲，别直接用最后 L−1 步。

陷阱 5：计算成本随 L·K·min(L,K) 涨。 SVD 本身是 $O(\min(L,K) \cdot L \cdot K)$ 。N=600、L=120 是毫秒级，但日线十年（~2500 点）配 L=600，SVD 会很慢；反复网格搜索 L 更慢。生产里建议固定 L、缓存 SVD、或只对感兴趣的子段做。

9. 小结：SSA 在分解家族里的位置#

方法	要预设什么	频率/周期	端点	混叠风险
傅里叶	整段周期	给出（但非平稳会糊）	周期延拓	低
小波	母小波 + 尺度	多尺度	边界处理	中
EMD	无（筛）	不显式给	端点效应	高（混叠）
VMD	K + α	显式给出	较稳	低
SSA	窗口 L	从 SVD 浮现	头尾缺失	低

一句话：SSA = 用「延迟嵌入 + SVD」把时间序列的自相关结构变成可排序的本征分量，趋势 / 周期 / 噪声靠你怎么分组来决定。 它不挑基、不拍 K/α，代价是要拍窗口 L、要靠经验配对分量、头尾补不齐。

下一篇我们对比过 VMD（把频谱当可估计量），而 SSA 走的是「矩阵分解」路线——两者殊途同归：都是把「趋势 + 周期 + 噪声」从一维信号里自适应地剥开。你的数据非平稳、想要频率自现，SSA 值得进工具箱。

代码与数据完全可复现：本文所有图由 gen_ssa_decomposition.py 生成，仅依赖 numpy + scipy.linalg（SVD），趋势 RMSE、scree 奇异值、周期还原相关、抗噪曲线均直接打印在脚本输出里。