頻譜分析 (Spectral Analysis，以下想簡稱的時候會叫 SA) 是過去沒學過的東西，此外也意外發現 Singular Spectrum Analysis(簡稱 SSA )，紀錄一下。
文很長，尤其是 SSA 的部分，超長。
Spectral Analysis 的目標是藉由將具有週期性(seasonal)的時間序列資料從 time domain 轉換為 frequency domain 的形式，擷取出訊號的特徵，再就這些特徵做出更好的預測或解釋。
從一個簡單的例子開始：
y1 <- rep(c(1,0,0),80)
y2 <- rep(c(0,0,0,2),60)
y3 <- rep(c(0,0,0,0,0,5),40)
y <- y1+y2+y3
ts_len <- 1:240
plot(1:20, y1[1:20], type="h", xlab="Time", ylab="Value", main ="y1, y2, y3 time series(head 20)")
lines(y2[1:20], col="red")
lines(y3[1:20], col="green")
plot(ts_len, y, type="h", xlab="Time", ylab="Value", main ="y time series")
y1、y2、y3各表示一個振動幅度的訊號，結合出來的y是這三種不同週期波的疊加，3個週期及幅度如下表：
| 變數 | 週期 | 幅度 | 重複次數 | 元素數 |
|---|---|---|---|---|
| y1 | 3 | 1 | 80 | 240 |
| y2 | 4 | 2 | 60 | 240 |
| y3 | 6 | 5 | 40 | 240 |
實務上會看到的資料可能就是y的長相，那麼，我們該如何辨識出潛藏在y後面的訊號呢? 可以使用 Fast Fourier Transform ，並且畫出功率譜（power spectrum）圖：
I <- abs(fft(y))^2/240
P <- (4/240)*I[1:120]
f <- 0:119/240
plot(f[-1],P[-1], type="l", xlab="Frequency", ylab="Power")
尖峰越高表示該頻率的能量(power)越強，週期性越明顯。以這個例子來說是週期為 3 的部分較明顯。
以上的動作在 R 可以只用一列 code 完成，不過這是經過平滑的版本，採用更高階的 Power Spectrum Analysis 方法處理:
spectrum(y, spans = 2, log = "no")
與經典的統計模型 ARIMA 家族應用時可能會出現的問題類似，Spectral Analysis 在由固定的頻率訊號組成的序列可以做出良好的預測，不過一旦實際序列中出現新的頻率的訊號，預測效果就會變差，因為模型根本沒有在訓練階段看過它。因此 Spectral Analysis 只能在具有穩定的週期且不會被突發事件影響的時間序列上才會有良好的表現。
(雖然我個人覺得這是預測長時間的時間序列都會面臨的問題
Forecasting: Principles and Practice也提出類似的解法：將最早的資料丟掉，只用近期的資料，避免模型僵化。雖然這主要針對 ARIMA 家族，但也適用於 Spectral Analysis。
Singular spectrum analysis (SSA) 是一種無母數方法(non-parametric technique)，這種方法的特色是不需要滿足統計假設即可使用，而且 SSA 除了擁有 decompensation 的能力外，搭配 linear recurrent formula 可以拿來做時間序列預測。
從 Singular spectrum analysis 的名詞顧名思義，它會把時間序列做 SVD 分解，不過時間序列都是一維的資料，要如何達成呢？其實 SSA 會經過以下步驟：
Embedding
Singular Value Decomposition
Grouping
Diagonal Averaging
讓我們來一步步理解這些是甚麼。
講到 Embedding ，會先想到的是修 NLP 時看到的 word embedding，但兩者的概念完全無關！它的概念為：把一維時間序列，轉成一個高維的矩陣，使後續可以執行 SVD 分解。具體來說：
給定一個一維且長度為 \(n\)時間序列\(\mathbf{y} = (y_1, y_2, \dots, y_n)'\)，首先要先選擇 window length \(L\)，代表「一次要看多長的一段歷史」，所以\(L\)的數學條件是：
選擇 \(L\) 的大小效果為：
再來，要建立延遲向量（lagged vectors），意思是從原始序列中，切出 \(K\) 個長度為 \(L\) 的小序列，其中 \(K = n - L + 1\)
然後對每個 \(i = 1, \dots, K\)，定義：
\[\mathbf{Y}_i = \begin{pmatrix} y_i \\ y_{i+1} \\ \vdots \\ y_{i+L-1} \end{pmatrix}\]
這其實是一個 Moving Window 的概念：
最後排成一個 trajectory matrix：
\[ \mathbf{Y}_{L \times K} = \begin{pmatrix} | & | & & | \\ \mathbf{Y}_1 & \mathbf{Y}_2 & \cdots & \mathbf{Y}_K \\ | & | & & | \end{pmatrix} \]
寫成元素形式：
\[ \mathbf{Y}_{L \times K} = \begin{pmatrix} y_1 & y_2 & \cdots & y_K \\ y_2 & y_3 & \cdots & y_{K+1} \\ \vdots & \vdots & \ddots & \vdots \\ y_L & y_{L+1} & \cdots & y_n \end{pmatrix} \]
這個矩陣是一個Hankel matrix，有一個很重要的特性是它每條反對角線上的元素都一樣。
這個步驟比較簡單，就是對 trajectory matrix 做 SVD 分解：
\[Y=UDV'\]
其中 \(U\) 和 \(V\) 有 orthonormal columns ，也就是說它的數學性質是：
\[\mathbf{U}_i' \mathbf{U}_j = 0,\quad\mathbf{V}_i' \mathbf{V}_j = 0\]
與 PCA 類似，SVD 拆解出來的東西統稱為 Components，運用在時間序列上就是有包含(signal)與雜訊(noise)的集合，從訊號再細分則有趨勢( trend )與季節性 ( seasonal )的部分，雜訊則是胡亂抖動、不規則的部份。而 \(U\) 簡單來說就是判斷這些 Components 屬於哪一種類型( pattern )。
\(D\) 是 diagonal matrix，\(D\)，即只有特徵值( eigenvalue )在對角元素的矩陣，且對角元素按大至小排序：\(d_1 \ge d_2 \ge \dots d_p\)。
\(V\)則是表示\(U\)在時間軸上是怎麼被啟動的，代表某個 component 在每個時間點的強弱。
這個 SVD 分解也可以改寫成 sum of rank-one bi-dimensional matrices 的形式
\[ \mathbf{Y} = \sum_{i=1}^{d} \mathbf{Y}_i = \sum_{i=1}^{d} \sqrt{\lambda_i} \mathbf{U}_i \mathbf{V}_i^{\prime} \]
文字說明：
\(d\)：Y 的 rank
\(\lambda_i\)：第 \(i\) 個特徵值
\(\mathbf{U}_i\)：左特徵向量
\(\mathbf{V}_i^{\prime}\)：右特徵向量的轉置
\(\mathbf{Y}_i = \lambda_i \mathbf{U}_i \mathbf{V}_i^{\prime}\)：第 \(i\) 個分量矩陣
考慮到 \(Y\) 的對稱結構，取 rank 時可以根據以下性質
\[\mathrm{rank}(Y) = \mathrm{rank}(YY') = \mathrm{rank}(Y'Y)\]
改寫成以下：
\[\mathbf{S} = YY'\]
這樣 \(S\) 的 rank-one bi-dimensional matrices 會變成
\[ \mathbf{S} = \mathbf{YY'} = \sum_{i=1}^{d} \lambda_i \mathbf{U}_i \mathbf{V}_i^{\prime}\mathbf{V}_i\mathbf{U}_i^{\prime} = \sum_{i=1}^{d} \lambda_i\]
就是特徵值的總和，也是 components 的權重(weight)總和！
Grouping 目的是要從所有特徵值中，挑出具代表性的訊號(signal)，其他則是雜訊(noise)。換句話說，我們可以把\(Y\)寫成
\[ \mathbf{Y} = \sum_{i=1}^{d} \sqrt{\lambda_i} \mathbf{U}_i \mathbf{V}_i^{\prime} + \varepsilon \]
\(\sqrt{\lambda_i} \mathbf{U}_i \mathbf{V}_i^{\prime}\)就是訊號+權重。
這一步的意思是把重建好的訊號矩陣，還原成一條序列，具體來說是要前述提取的訊號+權重重建成新的 Hankel matrix。要達成這個目標，有時也會用到 Smoothing 的方法讓重建出來的訊號更加平滑(沒有多餘的雜訊)，常見的作法是將所有多個反對角線的元素進行平均值化。
SSA 是一種無母數方法(Nonparametric method)，研究者需要控制 2 參數：Window length \(L\)以及選擇\(r\)個特徵作為訊號。那該怎麼設定才算是好的呢？以下有一個簡單的方法－看相關係數，用數學說，對每一群重建的\(𝑌^{(1)}_T\)與\(𝑌^{(2)}_T\)，他們的相關係數可以寫成：
\[𝜌_{(12)}^w = \frac{(𝑌^{(1)}_T,𝑌^{(2)}_T)}{ \|𝑌^{(1)}_T\|_w \times\|𝑌^{(2)}_T\|_w}\]
\(\|𝑌^{(i)}_T\|=\sqrt{(𝑌^{(1)}_T,𝑌^{(2)}_T)_w}\)
\((𝑌^{(1)}_T,𝑌^{(2)}_T)_w=\sum \limits _{k=1}^{t}w_ky^i_ky^j_k, i=1,j=2\)
如果兩者接近 0 相關，表示兩者可以好好分離，如果很高，表示兩群訊號有不少個訊號應該要被分在同一群。我們可以將其應用在分離訊號部分與雜訊部分的結果。這些\(Y\)兩兩相關的程度可以組成一個 W-correlation matrix (這裡的 w= 相關度) ，想要一目了然的查看這些\(Y\)兩兩相關的程度，可以用 W-correlation matrix 繪製的 heatmap。
根據文獻提供的說法，最適\(L\)值接近時間序列長度的一半，並且與每個時期的觀察數成比例，如此可以擁有較好的預測效果。不過如果時間許可，多加嘗試不同組合會是較好的選擇。
前述介紹的 SA 與 SSA 本質上都是一種對時間序列做分解，進而提取特徵的技術，沒有辦法拿來直接進行預測。要進行預測需要搭配一些演算法。最 basic 的方法是搭配 linear recurrent formula (LRF)，它看起來就像是 \(AR(d)\) 模型：
\[y_t=a_1y_{t-1}+a_2{t-1}+\dots+a_dy_{t-d} \\ t=d+1,\dots ,T\]
其中\(a_i\)系列就是賦予不同時間步的\(y\)的權重，它們可以經由前面分解出的 components 算出：
\[a^T=(a_{L-1}, \dots , a_1 )= \frac{1}{1-v^2}\sum \limits _{j=1}^{r}\pi_j U_j^{\bigtriangledown}\]
其中：
\(U_j^{\bigtriangledown}\)：前 \(L-1\) 個 components of \(U_j\)，\(U_j\) 來自 \(U\)
\(\pi_j\)：\(U_j\)最後一個 component
\(v^2=\sum \limits _{j=1}^{r}\pi_j^2\)
這個方法實務上比較常用，R 也有對應的 package 可用。
在 R 裡，想要使用 SSA 可以使用 RSSA package 來處理。
# install.packages("Rssa")
library(Rssa)
# 沿用 y 當例子
s <- ssa(y, L = 110) # L = window length
plot(s) # 看 eigenvalues
可以觀察到在 window length 設置為 110 的情況下，約莫有 8 個左右的特徵有權重(eigenvalues)，進一步去看看這些特徵的波動：
# 設定subplot繪圖參數
par(mfrow = c(4, 3))
# 繪製特徵向量 (Eigenvectors)
plot(s, type = "vectors", idx = 1:10, main = "Top 10 Eigenvectors")
可以觀察到，前 8 個 components 的波型很規律，可以設為訊號，後面佔比 0% 的 components 都是不規則狀，屬於雜訊。
觀察分群結果可以看看熱度圖：
plot(s,type = "wcor") # 看熱度圖
顏色越深越一致就越代表是同一群，可以觀察到整體區塊可以分成 2 個：前1-8個 components 圖片看起來很乾淨，第 9 個 lag 以後出現很多在不同區域的相關，這些即為雜訊，呼應雜訊的長相都有點不規則。
我們來利用 1-8 個 components 重建時間序列。
# 重建時間序列
r <- reconstruct(s, groups = list(Trend1= 1,Trend2= 2,Trend3=3:4,Trend4=5:6,Trend5=7:8))
# 畫重建後的序列
par(mfrow = c(3,2))
plot(r$Trend1, type = "l",main = "Trend1")
plot(r$Trend2, type = "l",main = "Trend2")
plot(r$Trend3, type = "l",main = "Trend3")
plot(r$Trend4, type = "l",main = "Trend4")
plot(r$Trend5, type = "l",main = "Trend5")
combined_series_ex1 <- r$Trend1+r$Trend2+r$Trend3+r$Trend4+r$Trend5
plot(combined_series_ex1, type = "l",main = "Trends")
par(mfrow = c(1,1))
plot(combined_series_ex1, type = "l",main = "Trends and Original Time series")
# 疊加原始序列
lines(y, col = "red")
# 加上圖例
legend("topright",
legend = c("Reconstructed", "Original"),
col = c("black", "red"), lty = 1)
因為選擇效果太好，直接用 8 個components 100%還原了，太強了！不過這是可能因為這個例子沒有真正雜訊擾亂的緣故。
不過這個例子有點太簡單了，再換一個。
先換成一個組成更複雜的例子：
改寫自 github 上分享的 ipynb，也是用 python 處理的版本。
# The number of time 'moments'
N <- 240
# Time index
t <- 0:(N - 1)
# Components
trend <- 0.001 * (t - 120)^2
p1 <- 40
p2 <- 70
periodic1 <- 2 * sin(2 * pi * t / p1)
periodic2 <- 0.75 * sin(2 * pi * t / p2)
# Set seed for reproducibility
set.seed(123)
# Generate noise
noise <- 2 * (runif(N) - 0.5)
# Final time series
y_two <- trend + periodic1 + periodic2 + noise
# Plot
plot(t, y_two, type = "l", lwd = 2.5,
xlab = "t", ylab = "F(t)",
main = "The Time Series y_tow and its Components")
lines(t, trend, col = "blue", lty = 2)
lines(t, periodic1, col = "red", lty = 2)
lines(t, periodic2, col = "green", lty = 2)
lines(t, noise, col = "purple", lty = 3)
legend("topright",
legend = c("Time Series (y_two)", "Trend",
"Periodic #1", "Periodic #2", "Noise"),
col = c("black", "blue", "red", "green", "purple"),
lty = c(1, 2, 2, 2, 3),
lwd = c(2.5, 1, 1, 1, 1))
在這個例子裡，y_two 包含了趨勢、2種週期性波動(或稱季節性波動)，以及雜訊，結構如下：
\[ Y_2 = \underbrace{0.001 (t-100)^2}_{\text{trend}} + \underbrace{2\sin\left(\frac{2\pi t}{20}\right)}_{\text{periodic 1}} + \underbrace{0.75\sin\left(\frac{2\pi t}{30}\right)}_{\text{periodic 2}} + \underbrace{\varepsilon_t}_{\text{noise}} \]
其中\(\varepsilon_t \sim \text{Uniform}(-1,1)\)，為 white noise。但是 trend 是非線性(non-linear)，非一般的 SA 與 ARMA 等線性模型可正確估計。1
2種週期性波動則是：
| 數學表示 | 週期 (Period) | 振幅 (Amplitude) | 頻率 (Frequency) | |
|---|---|---|---|---|
| 第一個 | \(2\sin\left(\frac{2\pi t}{20}\right)\) | 20 | 2 | \(1/20\) |
| 第二個 | \(0.75\sin\left(\frac{2\pi t}{30}\right)\) | 30 | 0.75 | \(1/30\) |
因為 20 和 30 的最小公倍數是 60，所以疊加後整體 pattern 每 60 單位會重複一次。
# install.packages("Rssa")
# library(Rssa)
# 用 y_two 當例子
s <- ssa(y_two, L = 110) # L = window length
plot(s) # 看 eigenvalues
可以觀察到在 window length 設置為 110 的情況下，共有 50 個 eigenvalues，約莫有 8 個左右的特徵權重較高，後面逐漸逼近 0 ，進一步去看看這些特徵的占比與波動：
# 設定subplot繪圖參數
par(mfrow = c(10, 5))
# 繪製特徵向量 (Eigenvectors)
plot(s, type = "vectors", idx = 1:20, main = "Top 20 Eigenvectors")
這是 \(U\) 裡，前 20 個 Eigenvectors (即\(U_j,\ j=1,\dots 20\))的波形圖，可以依照這些特徵數的權重占比去決定訊號與雜訊，此外，也可以藉由 \(U_j\) 的波形來判斷，如果波型不規則，極有可能是雜訊。
這裡可以先猜：第1、2個占比較多，設為趨勢；第3、4次之，為週期性波動；其他占比較少，且不少波形不規則，設定為雜訊。
也可以用 w-correlation 繪製的 heatmap 決定，並觀察目前選擇的 \(L\) 能不能有效分解：
plot(s,type = "wcor") # 看熱度圖
方塊顏色越深表示越有相關，如果 window length \(L\) 設得較好，理想情況屬於訊號的部分應該要盡量只有自己跟自己(如 F1 跟 F1)顏色最深，其他都是淺的 。以這裡來說，可以看到 F1 與 F2那裡只有自己顏色最深，可以選為 trend；F3 與 F4 顏色一樣深且最深，可以視做同一個 component；F5 與 F6 類似 F1 與 F2 的情況，但他們之間的關係度稍高，考慮到前面還有 1% 左右的占比，視做同一個 component，列入訊號中；之後的東西除了因為占比少，不同區域的方塊顏色也是不規則分布。
因此這樣猜測：
F1 與 F2 為trend
F3 與 F4 為週期 1
F5 與 F6 為週期 2
其他視作雜訊
我們利用 reconstruct() 執行前面提過的 LRF 演算法，重建時間序列。
# 重建時間序列
r <- reconstruct(s, groups = list(Trend = 1:2, Period1 = 3:4,Period2 = 5:6))
# r 裡面會有Trend, Period1, Period2 , original time series 等多維度資料，不建議直接用plot(r)
# 分開畫
par(mfrow = c(2,2))
plot(r$Trend, type = "l",main = "Trend")
plot(r$Period1, type = "l",main = "Period1")
plot(r$Period2, type = "l",main = "Period2")
combined_series <- r$Trend + r$Period1 + r$Period2
plot(combined_series, type = "l", main = "Reconstructed Trend + Period")
# plot(r$Rest, type = "l",
# main = "Remain")
右下角即是重建後的時間序列，跟原本的資料相比：
plot(combined_series, type = "l", main = "Reconstructed Signal")
lines(y_two,col="red") # 紅線為原始資料
顯然效果不錯，擷取訊號部分重建的時間序列可視作原本時間序列的平滑版本。
在預測方面，網路上的討論幾乎沒有看到 SSA 的缺點，都是讚美。這主要是是因為 SSA 用上了 SVD 的技術，使得它不需要嚴格的統計假設，研究者只需調整 \(L\) 的個數、選擇多少個特徵作為訊號重建跟使用什麼樣的平滑技巧。不過我在去年(2025/11月)的演講訪問到的教授 Paulo Canas Rodrigues 倒是說明了 SSA 本身可能會出現的問題，2原文翻譯如下：
單獨使用 SSA 並不會導致更好的預測準確性。原因是某些時間序列的訊號可能無法被我們選擇的 components 捕捉（如果我們選擇太少的components），或者該訊號可能包括一些噪音（如果我們選擇太多的 components ）。
SSA by itself does not result in better forecasting accuracy. The reason is that some signal from the time series might not be captured by the components that we select (if we select too few components), or that the signal might include some noise (if we select too many components).
這點在第一個例子中不太明顯，但第二個例子確實感覺得到選擇 components 需要一點感覺跟猜測，而且 SSA 是按照資料本身去分解，在不知道原資料有哪些訊號的情況不好處理。
我們來就第二個例子來驗證老師的說法。
# 重建時間序列
# 只選一個 component
r <- reconstruct(s, groups = list(Trend = 1))
# 分開畫
par(mfrow = c(2,1))
plot(r$Trend, type = "l",main = "Trend(only 1 component )")
plot(r$Trend, type = "l",main = "Trend and Original time series")
lines(y_two, col="red")
這個例子驗證了 component 選太少作為訊號的問題。如果是選太多：
# 重建時間序列
# 選30個 component
r_30 <- reconstruct(s, groups = list(Trend = 1:2, Guess1=3:4,Guess2=5:6, Guess3=6:7,Guess3=8:9, Guess4=10:30 ))
# 分開畫
par(mfrow = c(3,2))
plot(r_30$Trend, type = "l",main = "Trend")
plot(r_30$Guess1, type = "l",main = "Guess1")
plot(r_30$Guess2, type = "l",main = "Guess2")
plot(r_30$Guess3, type = "l",main = "Guess3")
plot(r_30$Guess4, type = "l",main = "Guess4")
combined_series_30 <- r_30$Trend + r_30$Guess1 + r_30$Guess2 +r_30$Guess3+ r_30$Guess4
plot(combined_series_30 , type = "l",main = "Signal and Original time series")
lines(y_two, col="red")
明顯可以看出Guess3與Guess4是雜訊，因為它有不規則的波形，誤把雜訊加入的結果是估計出來的序列變得很不平滑，這種情況像是迴歸裡的 overfitting，模型對未來新資料的預測能力很可能會不佳，解釋性也不高。
# 用 y_two 當例子
s2 <- ssa(y_two, L = 2) # L = window length
s2
# 設定subplot繪圖參數
par(mfrow = c(2, 1))
# 繪製特徵向量 (Eigenvectors)
plot(s, type = "vectors", idx = 1:2, main = "2 Eigenvectors")
Call:
ssa(x = y_two, L = 2)
Series length: 240, Window length: 2, SVD method: eigen
Special triples: 0
Computed:
Eigenvalues: 2, Eigenvectors: 2, Factor vectors: 0
Precached: 0 elementary series (0 MiB)
Overall memory consumption (estimate): 0.004471 MiB
dim(s2$U)只剩下 2 個 components 可以用啦 XDD 因為選擇 \(L\) 便超短的代價是SSA 最多只能有 L 個非零特徵值，\(U\)的維度也變成\(2 \times 2\)。看起來第 1 個 component 占比很大，可以當趨勢，第 2 個 component 當雜訊。
看一下 W-correlation
plot(s2,type = "wcor") # 看熱度圖
看起來分得很好 XDD ，再來看看重建會發生甚麼事
# 重建時間序列
r2 <- reconstruct(s2, groups = list(Trend = 1))
# r 裡面會有Trend, Period1, Period2 , original time series 等多維度資料，不建議直接用plot(r)
# 分開畫
par(mfrow = c(2,2))
plot(r2$Trend, type = "l",main = "Trend")
plot(combined_series, type = "l", main = "Reconstructed Signal(L=110)")
plot(r2$Trend, type = "l",main = "Trends and Original Time series")
lines(y_two,col="red") # 紅線為原始資料
問題變得比較明顯了，SSA 的目標是要有效分離出訊號與雜訊，但是因為 components 只有兩個，即使挑其中一個占比大的，重建時資料也變得會捕捉到雜訊，不像 \(L=110\) 的版本要平滑。這種情況像是迴歸裡的 overfitting，模型對未來新資料的預測能力會不佳。
再來看看 \(L\) 設超大會怎樣：
# 用 y_two 當例子
s2 <- ssa(y_two, L = 239) # L = window length
s2
# 設定subplot繪圖參數
par(mfrow = c(2, 1))
# 繪製特徵向量 (Eigenvectors)
plot(s, type = "vectors", idx = 1:2, main = "2 Eigenvectors")
Call:
ssa(x = y_two, L = 239)
Series length: 240, Window length: 239, SVD method: eigen
Special triples: 0
Computed:
Eigenvalues: 2, Eigenvectors: 2, Factor vectors: 0
Precached: 0 elementary series (0 MiB)
Overall memory consumption (estimate): 0.008087 MiB
還是只剩兩個 components ，看來問題會與 \(L\) 設超小一樣。
SA 部分：
SSA 部分