多维奇异谱分析（Multivariate Singular Spectrum Analysis，MSSA）
1. HMSSA1. 1 分解1.1.1 第一步 嵌入1.1.2 第二步 SVD1.2 重构1.2.1 第三步 分组1.2.2 第四步 对角平均1.3 预测1.3.1 recurrent forecast1.3.2 vector forecast2. VMSSA2.1 分解2.1.1 第一步 嵌入2.1.2 第二步 SVD2.2 重构2.2.1 第三步 分组2.2.2 第四步 对角平均2.3 预测2.3.1 recurrent forecast2.3.2 vector forecast在实际运用中，通常需要对多个时间序列进行分析，而每个时间序列都有内部结构，且序列之间也会存在一定的依赖关系。因此，基本的单变量奇异谱分析(SSA)可以扩展到多元情况，也就是多维奇异谱分析(MSSA)。
轨迹矩阵是求解SSA的主要工具，对于多变量情况，轨迹矩阵可以有不同的定义方式。主要可以分为水平和垂直两种方式。而每一种定义方式下，又可以有recurrent和vector两种预测方法，如下图所示：
下面就主要梳理一下MSSA的内容。和SSA一样，MSSA也由两个互补阶段组成:分解和重建。第一阶段对序列进行分解，第二阶段对无噪声序列进行重构。
1. HMSSA
1. 1 分解
1.1.1 第一步 嵌入
考虑M个时间序列，先根据SSA的方法，生成第i(i=1,...,M)i(i=1,...,M)i(i=1,...,M)个序列的轨迹矩阵X(i)X^{(i)}X(i)，而后，水平地将这些轨迹矩阵拼接在一起：
X=[X(1):X(2):...:X(M)]X=[X^{(1)}:X^{(2)}:...:X^{(M)}]X=[X(1):X(2):...:X(M)]
要顺利的实现拼接，假设第iii个序列的轨迹矩阵X(i)X^{(i)}X(i)是Li∗KiL_i*K_iLi∗Ki的矩阵，其中，Li(2≤Li≤Ni−1)L_i(2\le L_i\le N_i-1)Li(2≤Li≤Ni−1)是每个长度为NiNiNi的时间序列的窗口长度，Ki=Ni−Li+1,i=1,...,MK_i=N_i-L_i+1, i=1,...,MKi=Ni−Li+1,i=1,...,M，则有L1=...=LM=L,L_1=...=L_M=L,L1=...=LM=L,而Ki,NiK_i,N_iKi,Ni可以不同。
1.1.2 第二步 SVD
这一步，对前面得到的XHX_HXH进行SVD操作，用λH1,...,λHL\lambda_{H_1},...,\lambda_{H_L}λH1,...,λHL表示XHXHTX_HX_H^TXHXHT的特征值，同样，有λH1≥,...,≥λHL≥0\lambda_{H_1} \ge ,...,\ge\lambda_{H_L}\ge0λH1≥,...,≥λHL≥0，UH1,...,UHLU_{H_1},...,U_{H_L}UH1,...,UHL则表示对应的特征向量。
对XHX_HXH进行SVD的过程可以表示为:XH=XH1+...+XHLX_H=X_{H_1}+...+X_{H_L}XH=XH1+...+XHL，其中
因此，XHXHTX_HX_H^TXHXHT的结构等价于：
1.2 重构
1.2.1 第三步 分组
与SSA类似，这个步骤相当于分解矩阵XH1,...,XHLX_{H_1},...,X_{H_L}XH1,...,XHL分成几个不相交的组，并对每组中的矩阵求和。下标集合1,…,L{1,…,L}1,…,L被划分成不相交的子集I1,……,ImI_1,……,I_mI1,……,Im,与XH=XI1+...+XImX_H=X_{I_1}+...+X_{I_m}XH=XI1+...+XIm相对应。
1.2.2 第四步 对角平均
对角平均的目的是将重构的hankel矩阵XIjX_{I_j}XIj变换为一维的序列。对于每个重构的轨迹矩阵而言，这一步与SSA相同。
1.3 预测
1.3.1 recurrent forecast
对于固定的LLL，为每个序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)分别构造轨迹矩阵X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,KiX^{(i)}= [X^{(i)}_{1},...,X^{(i)}_{K}]=(x_{mn})^{L,K_i}_{m,n=1}X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,Ki
构造块轨迹矩阵XVX_VXV如下：
定义UVj=(Uj(1)，…，Uj(M))T,U_{V_j}=(U_{j}^{(1)}，…，U_{j}^{(M)})^T,UVj=(Uj(1)，…，Uj(M))T, 作为一个XVXVTX_VX_V^TXVXVT的第jjj个特征向量。其中Uj(i)U_{j}^{(i)}Uj(i)的长度为L，对应于序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)
定义X^V=[X^1:...:X^K]=Σi=1rUViUViTXV\hat{X}_V=[\hat{X}_1:...:\hat{X}_K]=\Sigma_{i=1}^{r}U_{V_i}U^T_{V_i}X_VX^V=[X^1:...:X^K]=Σi=1rUViUViTXV,即
考虑矩阵
作为从上一步中获得的矩阵X^(i)\hat{X}^{(i)}X^(i)汉克尔化过程的结果
令Uj(i)∇U^{(i)\nabla}_{j}Uj(i)∇是特征向量Uj(i)U_{j}^{(i)}Uj(i)的前L−1L-1L−1个元素构成的向量。πj(i)\pi_{j}^{(i)}πj(i)表示特征向量Uj(i)U_{j}^{(i)}Uj(i)的最后一个元素。(i=1,...,M)(i=1,...,M)(i=1,...,M)
选择r个重建阶段的特征三元组用于预测。
定义矩阵U∇M=(U1∇M,...,Ur∇M)U^{\nabla M}=(U^{\nabla M}_1,...,U^{\nabla M}_r)U∇M=(U1∇M,...,Ur∇M),其中，Uj∇MU^{\nabla M}_jUj∇M为：
定义矩阵WWW:
如果矩阵(IM×M−WWT)−1(I_{M×M}−WW^T)^{−1}(IM×M−WWT)−1存在，且r≤Lsum−Mr≤L_{sum}−Mr≤Lsum−M，则存在h步VMSA预测，公式如下：
其中，y~ji(m)\tilde{y}^{(m)}_{j_i}y~ji(m)是重构序列mmm的第jij_iji个观察值，m=1,...,M.m=1,...,M.m=1,...,M.
1.3.2 vector forecast
HMSSA-V的过程非常类似于它的单变量版本SSA-V和HMSSA-R．算法的第一部分与HMSSA-R一样，继续考虑以下矩阵：
定义向量Zj(i),(i=1,...M)Z_j^{(i)},(i=1,...M)Zj(i),(i=1,...M)如下：
其中，X~j(i)\tilde{X}^{(i)}_{j}X~j(i)是分组并保留噪声分量后，第i个序列轨迹矩阵的重构列。通过构造矩阵Z(i)=[Z1(i),...,Zki+h+L−1(i)]Z^{(i)}= [Z^{(i)}_1,...,Z^{(i)}_{ k_i+h+L−1}]Z(i)=[Z1(i),...,Zki+h+L−1(i)]，并求对角平均，我们可以得到一个新的序列y^1(i),...,y^Ni+h(i)\hat{y}_1^{(i)},...,\hat{y}_{N_i+h}^{(i)}y^1(i),...,y^Ni+h(i)，其中，y^N+1(i),...,y^N+h(i)\hat{y}_{N+1}^{(i)},...,\hat{y}_{N+h}^{(i)}y^N+1(i),...,y^N+h(i)则是通过HMSSA-V方法获得的hhh个预测值。
2. VMSSA
2.1 分解
2.1.1 第一步 嵌入
与HMSSA一样，考虑M个时间序列，先根据SSA的方法，生成第i(i=1,...,M)i(i=1,...,M)i(i=1,...,M)个序列的轨迹矩阵X(i)X^{(i)}X(i)
而后，垂直地将这些轨迹矩阵拼接在一起，得到XVX_VXV：
要顺利的实现拼接，假设第iii个序列的轨迹矩阵X(i)X^{(i)}X(i)是Li∗KiL_i*K_iLi∗Ki的矩阵，其中，Li(2≤Li≤Ni−1)L_i(2\le L_i\le N_i-1)Li(2≤Li≤Ni−1)是每个长度为NiNiNi的时间序列的窗口长度，Ki=Ni−Li+1,i=1,...,MK_i=N_i-L_i+1, i=1,...,MKi=Ni−Li+1,i=1,...,M，与HMSSA不同的是，VMSSA要求K1=...=KM=K,K_1=...=K_M=K,K1=...=KM=K,而Li,NiL_i,N_iLi,Ni可以不同。
2.1.2 第二步 SVD
这一步，对前面得到的XVX_VXV进行SVD操作，用λV1,...,λVL\lambda_{V_1},...,\lambda_{V_L}λV1,...,λVL表示XVXVTX_VX_V^TXVXVT的特征值，同样，有λV1≥,...,≥λVL≥0\lambda_{V_1} \ge ,...,\ge\lambda_{V_L}\ge0λV1≥,...,≥λVL≥0，UV1,...,UVLU_{V_1},...,U_{V_L}UV1,...,UVL则表示对应的特征向量。
XVXVTX_VX_V^TXVXVT的结构等价于：
2.2 重构
2.2.1 第三步 分组
与HMSSA类似，这个步骤相当于分解矩阵XV1,...,XVLX_{V_1},...,X_{V_L}XV1,...,XVL分成几个不相交的组，并对每组中的矩阵求和。下标集合1,…,L{1,…,L}1,…,L被划分成不相交的子集I1,……,ImI_1,……,I_mI1,……,Im,与XV=XI1+...+XImX_V=X_{I_1}+...+X_{I_m}XV=XI1+...+XIm相对应。
2.2.2 第四步 对角平均
对角平均的目的是将重构的hankel矩阵XIjX_{I_j}XIj变换为一维的序列。对于每个重构的轨迹矩阵而言，这一步与SSA相同。
2.3 预测
2.3.1 recurrent forecast
考虑MMM个序列YNi(i)=(y1(i),...,yNi(i))Y^{(i)}_{N_i}=(y^{(i)}_{1},...,y^{(i)}_{N_i})YNi(i)=(y1(i),...,yNi(i))和对应的窗口长度Li,2≤Li≤Ni−1,i=1...,ML_i,2≤Li≤Ni−1,i=1...,MLi,2≤Li≤Ni−1,i=1...,M，VMSSA-R预测算法如下所示。
对于固定的KKK，为每个序列YNi(i)(i=1,...,M)Y^{(i)}_{N_i}(i=1,...,M)YNi(i)(i=1,...,M)分别构造轨迹矩阵X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,KiX^{(i)}= [X^{(i)}_{1},...,X^{(i)}_{K}]=(x_{mn})^{L,K_i}_{m,n=1}X(i)=[X1(i),...,XK(i)]=(xmn)m,n=1L,Ki
构造块轨迹矩阵XHX_HXH如下：
定义UHj=(u1j，…，uLj)T,U_{H_j}=(u_{1j}，…，u_{Lj})^T,UHj=(u1j，…，uLj)T, 长度为L，作为一个XHXHTX_HX_H^TXHXHT的第jjj个特征向量。
定义Hankel算子Xi^\hat{X_i}Xi^,考虑由r个分量获得的重构矩阵XH^=Σi=1rUHiUHiTXH\hat{X_H}=\Sigma_{i=1}^{r}U_{H_i}U^T_{H_i}X_HXH^=Σi=1rUHiUHiTXH
考虑矩阵
作为从上一步中获得的矩阵X^(i)\hat{X}^{(i)}X^(i)汉克尔化过程的结果
令UHj∇U^\nabla_{H_j}UHj∇是特征向量UHjU_{H_j}UHj的前L−1L-1L−1个元素构成的向量。πHj\pi_{H_j}πHj表示特征向量UHjU_{H_j}UHj的最后一个元素。(j=1,...,r)(j=1,...,r)(j=1,...,r)
定义ν2=Σj=1rπHj2\nu^2=\Sigma_{j=1}^{r}\pi_{H_j}^2ν2=Σj=1rπHj2
定义线性系数向量：
如果ν2<1\nu^2<1ν2<1，则存在HMSSA的H步预测，计算公式如下：
其中，i=1,...,M.i=1,...,M.i=1,...,M.
2.3.2 vector forecast
此算法的前10部分与VMSSA-R相同，继续考虑矩阵：
其中，
可以证明，矩阵Π\PiΠ是U∇U^{\nabla}U∇的列空间上的正交投影算子。令
其中，Π(i)\Pi^{(i)}Π(i)的维度为(Li−1)×(Lsum−M)(L_i−1)×(L_{sum}−M)(Li−1)×(Lsum−M),而R(i)(i=1,...,M)\mathscr{R}^{(i)}(i=1,...,M)R(i)(i=1,...,M)的长度为Lsum−ML_{sum}−MLsum−M,对应于序列YNi(i)Y_{N_i}^{(i)}YNi(i).然后，定义线性运算：
其中，YΔT=(YΔ(1),...,YΔ(M))Y_\Delta^T=(Y_\Delta^{(1)},...,Y_\Delta^{(M)})YΔT=(YΔ(1),...,YΔ(M)),而YΔ(i),(i=1,...,M)Y_\Delta^{(i)},(i=1,...,M)YΔ(i),(i=1,...,M)表示YiY_iYi(长度为LiL_iLi)的最后Li−1L_i-1Li−1个元素。使用上面的符号，VMSSA-V预测算法如下：
定义向量ZiZ_iZi如下：
其中，X~i\tilde{X}_{i}X~i是分组并保留噪声分量后，第i个序列轨迹矩阵的重构列。Lmax=max{L1,...,LM}L_{max}=max\{L_1,...,L_M\}Lmax=max{L1,...,LM}
2. 构造矩阵Z=[Z1,...,ZK+h+Lmax−1]Z= [Z_1,...,Z_{ K+h+L_{max}−1}]Z=[Z1,...,ZK+h+Lmax−1]，并求对角平均，我们可以得到一个新的序列y^1(i),...,y^N+h+Lmax(i)(i=1,...,M)\hat{y}_1^{(i)},...,\hat{y}_{N+h+L_{max}}^{(i)}(i=1,...,M)y^1(i),...,y^N+h+Lmax(i)(i=1,...,M)
3. 元素y^Ni+1(i),...,y^Ni+h(i)(i=1,...,M)\hat{y}_{N_i+1}^{(i)},...,\hat{y}_{N_i+h}^{(i)}(i=1,...,M)y^Ni+1(i),...,y^Ni+h(i)(i=1,...,M)则是通过VMSSA-V方法获得的hhh个预测值。
