奇异谱分析SSA

这两个月跟着老师做项目，接触了一点时间序列分析和预测的内容，主要是基于矩阵分解的奇异谱分析（Singular Spectrum Analysis，SSA），趁国庆放假，整理一下做个笔记。

1. SSA （Singular Spectrum Analysis）

奇异谱分析(SSA)是一种处理非线性时间序列数据的方法，可以对时间序列进行分析和预测。它基于构造在时间序列上的特定矩阵的奇异值分解(SVD)，可以从一个时间序列中分解出趋势、振荡分量和噪声。SSA具有非常广泛的适用性，对于时间序列，既不需要假设参数模型，也不需要假设平稳性条件。

1.1 SSA算法的基本流程

考虑一个长度为N的时间序列$X^{\cdot }=X_N^{\cdot }=(x_1，\dots ，x_N)$。$N>2$，且$X$是一个非零序列，即，至少存在一个$i$使得$x_i\ne 0$。令整数$L(1<L<N)$为窗口长度，且$K=N_L+1$
SSA算法的过程由分解和重构两个互补的阶段组成。

1.2 分解

1. 第一步　嵌入
   我们将原始时间序列映射成一个长度为$L$的向量序列，形成$K=N-L+1$个长度为$L$的向量：

这些向量组成轨迹矩阵：

2. 第二步:奇异值分解
   在这一步，对轨迹矩阵$X$进行奇异值分解(SVD)。令$S=X{X}^T$,${\lambda }_1,...,{\lambda }_L$为$S$的特征值，且${\lambda }_1\ge ...\ge {\lambda }_L\ge 0$；而$U_1,...,U_L$是矩阵$S$对应于这些特征值的标准正交向量。
   令 d = r a n k ( X ) = m a x { i , λ i > 0 } d=rank(X)=max\{i,λ_i>0\} d=rank(X)=max{i,λi​>0}(注意，在实际序列中，我们通常有 d = L ∗ , L ∗ = m i n { L , K } d=L^∗,L^∗=min\{L,K\} d=L∗,L∗=min{L,K})。$V_i=X^T{U}_i/\sqrt{{\lambda }_i}(i=1，\dots ，d)$
   这种情况下，轨迹矩阵$X$的SVD可以写成:
   

其中，$X_i=\sqrt{{\lambda }_i}{U}_i{V}_i^T$

1.3 重构

3. 第三步 分组(grouping)
   先将下标集合 { 1 , . . . , d } \{1,...,d\} {1,...,d}划分成$m$个互不相交的子集$I_1,...I_m$，令 I = { i 1 , . . . , i p } I=\{i_1,...,i_p\} I={i1​,...,ip​}，则对应于$I$的合成矩阵$X_I=X_{i_1}+...+X_{i_p}$。则有：
   

4. 第四步:对角线平均
   在这一步中，我们将(2.3)中的每个矩阵$X_{I_j}$变换为一个长度为N的新序列，即得到分解后的序列。令$Y$为一个$L\ast K$的矩阵，元素为$y_{ij},1\le i\le L,1\le j\le K.$令$L^{\ast }=min(L,K),K^{\ast }=max(L,K),N=L+K-1$.如果$L<K$,$y_{ij}^{\ast }=y_{ij}$,否则,$y_{ij}^{\ast }=y_{ji}$
   我们利用下面的公式2.4进行对角线求平均,将矩阵$Y$转换为序列 $y_1，\dots ,y_N$
   
   即根据下图所示对角线求平均，将二维矩阵转换为一维序列：
   

1.4 SSA预测

SSA的预测，可以简单理解为一种线性递归过程，即：

$y_{N+1}=a_1{y}_N+a_2{y}_{N-1}+...+a_{L-1}{y}_{N-L+1}$

其中，系数$a_1,...a_{L-1}$根据SVD获得的特征值计算得到。

SSA的预测方法有递归(recurrent)和矩阵(vector)两种。

1.4.1 Recurrent Forecast

时间序列$Y_{N+M}=(y_1,...,y_{N+M})$的递归预测公式如下：

其中,$\widetilde{x_i}为根据公式2.4重构出来的时间序列值。$$y_{N+1},...,y_{N+M}$为M个预测值。
则，我们只需要求出系数$a_1,...a_{L-1}$即可。计算方法如下：
记向量$R=(a_{L-1},...,a_1{)}^T$,且有

式3.1中，${\nu }^2={\pi }_1^2+...+{\pi }_r^2$; ${\pi }_i$是向量$P_i(i=1,...,r)$的最后一个分量,$P_i$则是在SVD分解过程中得到的标准正交向量。

1.4.2 Vector Forecast

SSA预测的另一种方法是向量预测法。总的来说，向量预测比递归预测更稳定，特别是当序列中异常值较多的时候。考虑下面的矩阵：

定义线性算子如下：

其中，$Y_{\Delta }$是$Y$的后$L-1$个元素构成的向量。定义向量$Z_j$如下：

其中，$\widetilde{X_j}$为轨迹矩阵经过分组并剔除噪声分量后重构的第$j$列。现在，通过构造矩阵$Z=[Z_1,...,Z_{K+h+L-1}]$，并进行对角线平均，可以得到一个新的序列${\hat{y}}_1,...,{\hat{y}}_{K+h+L-1}$，其中，${\hat{y}}_{N+1},...,{\hat{y}}_{N+h}$则是通过Vector方法获得预测值。
Vector方法预测过程如下图所示：

2. python实现

SSA的代码实现，还是比较多的，基本上都是实现了分解和重构，预测一般也都是recurrent方法，如
pyts中的SSA
GitHub上也有比较多：
pssa
ssa-py

我基于GitHub上面的pssa，修改了一下，粗略实现了vector预测方法，有空的时候再整理传到GitHub吧。