多维奇异谱分析（Multivariate Singular Spectrum Analysis，MSSA)

在实际运用中，通常需要对多个时间序列进行分析，而每个时间序列都有内部结构，且序列之间也会存在一定的依赖关系。因此，基本的单变量奇异谱分析(SSA)可以扩展到多元情况，也就是多维奇异谱分析(MSSA)。
轨迹矩阵是求解SSA的主要工具，对于多变量情况，轨迹矩阵可以有不同的定义方式。主要可以分为水平和垂直两种方式。而每一种定义方式下，又可以有recurrent和vector两种预测方法，如下图所示：

下面就主要梳理一下MSSA的内容。和SSA一样，MSSA也由两个互补阶段组成:分解和重建。第一阶段对序列进行分解，第二阶段对无噪声序列进行重构。

1. HMSSA

1. 1 分解

1.1.1 第一步 嵌入

考虑M个时间序列，先根据SSA的方法，生成第$i(i=1,...,M)$个序列的轨迹矩阵$X^{(i)}$，而后，水平地将这些轨迹矩阵拼接在一起：

$X=[X^{(1)}:X^{(2)}:...:X^{(M)}]$

要顺利的实现拼接，假设第$i$个序列的轨迹矩阵$X^{(i)}$是$L_i\ast K_i$的矩阵，其中，$L_i(2\le L_i\le N_i-1)$是每个长度为$Ni$的时间序列的窗口长度，$K_i=N_i-L_i+1,i=1,...,M$，则有$L_1=...=L_M=L,$而$K_i,N_i$可以不同。

1.1.2 第二步 SVD

这一步，对前面得到的$X_H$进行SVD操作，用${\lambda }_{H_1},...,{\lambda }_{H_L}$表示$X_H{X}_H^T$的特征值，同样，有${\lambda }_{H_1}\ge ,...,\ge {\lambda }_{H_L}\ge 0$，$U_{H_1},...,U_{H_L}$则表示对应的特征向量。
对$X_H$进行SVD的过程可以表示为:$X_H=X_{H_1}+...+X_{H_L}$，其中

因此，$X_H{X}_H^T$的结构等价于：

1.2 重构

1.2.1 第三步 分组

与SSA类似，这个步骤相当于分解矩阵$X_{H_1},...,X_{H_L}$分成几个不相交的组，并对每组中的矩阵求和。下标集合$1,\dots ,L$被划分成不相交的子集$I_1,\dots \dots ,I_m$,与$X_H=X_{I_1}+...+X_{I_m}$相对应。

1.2.2 第四步 对角平均

对角平均的目的是将重构的hankel矩阵$X_{I_j}$变换为一维的序列。对于每个重构的轨迹矩阵而言，这一步与SSA相同。

1.3 预测

1.3.1 recurrent forecast

1. 对于固定的$L$，为每个序列$Y_{N_i}^{(i)}(i=1,...,M)$分别构造轨迹矩阵$X^{(i)}=[X_1^{(i)},...,X_K^{(i)}]=(x_{mn}{)}_{m,n=1}^{L,K_i}$

2. 构造块轨迹矩阵$X_V$如下：
   

3. 定义$U_{V_j}=(U_j^{(1)}，\dots ，U_j^{(M)}{)}^T,$ 作为一个$X_V{X}_V^T$的第$j$个特征向量。其中$U_j^{(i)}$的长度为L，对应于序列$Y_{N_i}^{(i)}(i=1,...,M)$

4. 定义${\hat{X}}_V=[{\hat{X}}_1:...:{\hat{X}}_K]={\Sigma }_{i=1}^r{U}_{V_i}{U}_{V_i}^T{X}_V$,即
   

5. 考虑矩阵
   作为从上一步中获得的矩阵${\hat{X}}^{(i)}$汉克尔化过程的结果

6. 令$U_j^{(i)\nabla }$是特征向量$U_j^{(i)}$的前$L-1$个元素构成的向量。${\pi }_j^{(i)}$表示特征向量$U_j^{(i)}$的最后一个元素。$(i=1,...,M)$

7. 选择r个重建阶段的特征三元组用于预测。

8. 定义矩阵$U^{\nabla M}=(U_1^{\nabla M},...,U_r^{\nabla M})$,其中，$U_j^{\nabla M}$为：
   

9. 定义矩阵$W$:
   

10. 如果矩阵$(I_{M\times M}-W{W}^T{)}^{-1}$存在，且$r\le L_{sum}-M$，则存在h步VMSA预测，公式如下：
    
    其中，${\tilde{y}}_{j_i}^{(m)}$是重构序列$m$的第$j_i$个观察值，$m=1,...,M.$

1.3.2 vector forecast

HMSSA-V的过程非常类似于它的单变量版本SSA-V和HMSSA-R．算法的第一部分与HMSSA-R一样，继续考虑以下矩阵：

1. 定义向量$Z_j^{(i)},(i=1,...M)$如下：
   
   其中，${\tilde{X}}_j^{(i)}$是分组并保留噪声分量后，第i个序列轨迹矩阵的重构列。
2. 通过构造矩阵$Z^{(i)}=[Z_1^{(i)},...,Z_{k_i+h+L-1}^{(i)}]$，并求对角平均，我们可以得到一个新的序列${\hat{y}}_1^{(i)},...,{\hat{y}}_{N_i+h}^{(i)}$，其中，${\hat{y}}_{N+1}^{(i)},...,{\hat{y}}_{N+h}^{(i)}$则是通过HMSSA-V方法获得的$h$个预测值。

2. VMSSA

2.1 分解

2.1.1 第一步 嵌入

与HMSSA一样，考虑M个时间序列，先根据SSA的方法，生成第$i(i=1,...,M)$个序列的轨迹矩阵$X^{(i)}$

而后，垂直地将这些轨迹矩阵拼接在一起，得到$X_V$：

要顺利的实现拼接，假设第$i$个序列的轨迹矩阵$X^{(i)}$是$L_i\ast K_i$的矩阵，其中，$L_i(2\le L_i\le N_i-1)$是每个长度为$Ni$的时间序列的窗口长度，$K_i=N_i-L_i+1,i=1,...,M$，与HMSSA不同的是，VMSSA要求$K_1=...=K_M=K,$而$L_i,N_i$可以不同。

2.1.2 第二步 SVD

这一步，对前面得到的$X_V$进行SVD操作，用${\lambda }_{V_1},...,{\lambda }_{V_L}$表示$X_V{X}_V^T$的特征值，同样，有${\lambda }_{V_1}\ge ,...,\ge {\lambda }_{V_L}\ge 0$，$U_{V_1},...,U_{V_L}$则表示对应的特征向量。

$X_V{X}_V^T$的结构等价于：

2.2 重构

2.2.1 第三步 分组

与HMSSA类似，这个步骤相当于分解矩阵$X_{V_1},...,X_{V_L}$分成几个不相交的组，并对每组中的矩阵求和。下标集合$1,\dots ,L$被划分成不相交的子集$I_1,\dots \dots ,I_m$,与$X_V=X_{I_1}+...+X_{I_m}$相对应。

2.2.2 第四步 对角平均

对角平均的目的是将重构的hankel矩阵$X_{I_j}$变换为一维的序列。对于每个重构的轨迹矩阵而言，这一步与SSA相同。

2.3 预测

2.3.1 recurrent forecast

考虑$M$个序列$Y_{N_i}^{(i)}=(y_1^{(i)},...,y_{N_i}^{(i)})$和对应的窗口长度$L_i,2\le Li\le Ni-1,i=1...,M$，VMSSA-R预测算法如下所示。

1. 对于固定的$K$，为每个序列$Y_{N_i}^{(i)}(i=1,...,M)$分别构造轨迹矩阵$X^{(i)}=[X_1^{(i)},...,X_K^{(i)}]=(x_{mn}{)}_{m,n=1}^{L,K_i}$

2. 构造块轨迹矩阵$X_H$如下：
   

3. 定义$U_{H_j}=(u_{1j}，\dots ，u_{Lj}{)}^T,$ 长度为L，作为一个$X_H{X}_H^T$的第$j$个特征向量。

4. 定义Hankel算子$\widehat{X_i}$,考虑由r个分量获得的重构矩阵$\widehat{X_H}={\Sigma }_{i=1}^r{U}_{H_i}{U}_{H_i}^T{X}_H$
   

5. 考虑矩阵
   作为从上一步中获得的矩阵${\hat{X}}^{(i)}$汉克尔化过程的结果

6. 令$U_{H_j}^{\nabla }$是特征向量$U_{H_j}$的前$L-1$个元素构成的向量。${\pi }_{H_j}$表示特征向量$U_{H_j}$的最后一个元素。$(j=1,...,r)$

7. 定义${\nu }^2={\Sigma }_{j=1}^r{\pi }_{H_j}^2$

8. 定义线性系数向量：
   

9. 如果${\nu }^2<1$，则存在HMSSA的H步预测，计算公式如下：
   
   其中，$i=1,...,M.$

2.3.2 vector forecast

此算法的前10部分与VMSSA-R相同，继续考虑矩阵：

其中，

可以证明，矩阵$\Pi$是$U^{\nabla }$的列空间上的正交投影算子。令

其中，${\Pi }^{(i)}$的维度为$(L_i-1)\times (L_{sum}-M)$,而${\mathcal{R}}^{(i)}(i=1,...,M)$的长度为$L_{sum}-M$,对应于序列$Y_{N_i}^{(i)}$.然后，定义线性运算：

其中，$Y_{\Delta }^T=(Y_{\Delta }^{(1)},...,Y_{\Delta }^{(M)})$,而$Y_{\Delta }^{(i)},(i=1,...,M)$表示$Y_i$(长度为$L_i$)的最后$L_i-1$个元素。使用上面的符号，VMSSA-V预测算法如下：

1. 定义向量$Z_i$如下：
   

其中，${\tilde{X}}_i$是分组并保留噪声分量后，第i个序列轨迹矩阵的重构列。 L m a x = m a x { L 1 , . . . , L M } L_{max}=max\{L_1,...,L_M\} Lmax​=max{L1​,...,LM​}
2. 构造矩阵$Z=[Z_1,...,Z_{K+h+L_{max}-1}]$，并求对角平均，我们可以得到一个新的序列${\hat{y}}_1^{(i)},...,{\hat{y}}_{N+h+L_{max}}^{(i)}(i=1,...,M)$
3. 元素${\hat{y}}_{N_i+1}^{(i)},...,{\hat{y}}_{N_i+h}^{(i)}(i=1,...,M)$则是通过VMSSA-V方法获得的$h$个预测值。