SSA时间序列分解

奇异谱分析（SSA)根据观测到的时间序列构造轨迹矩阵，并对轨迹矩阵进行分解和重构，从而提取出代表原时间序列不同成分的信号，如长期趋势信号、周期信号、噪声信号等，从而进一步对分解得到的信号进行分析。

算法流程

Step1:根据原始时间序列构建轨迹矩阵$X$

Step2:对矩阵X进行奇异值分解：$X={\sum }_{i=1}^r{\sigma }_i{U}_i{V}_i^T$

Step3:按奇异值生成$r$个子矩阵:$X_i={\sigma }_i{U}_i{V}_i^T$

Step4:根据某一分组原则将子矩阵$X_i$分为$m$个组

Step5:对子矩阵$X_i$进行对角均值化处理得到子序列$\widetilde{F_i}$

Step6:对$m$个组中的子序列相加得到分组子序列。

矩阵分解

时间序列嵌入形成轨迹矩阵

输入：原始时间序列$y$，窗口长度$L(2\le L\le N/2)$，$N=len(y)$
$$\mathbf{X}=\left[\begin{array}{cccccc}{y}_1& y_2& y_3& y_4& \dots & y_{N-L+1}\\ y_2& y_3& y_4& y_5& \dots & y_{N-L+2}\\ y_3& y_4& y_5& y_6& \dots & y_{N-L+3}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ y_L& y_{L+1}& y_{L+2}& y_{L+3}& \dots & y_N\end{array}\right]\text{\hspace{1em}(1)}$$
显然矩阵$X$是一个汉克尔矩阵（每一条副对角线的元素都相等），矩阵的行数为窗口长度$L$，矩阵的列数为$N-L+1$。

奇异值分解（SVD)

$$X=U\Sigma V^T\text{\hspace{1em}(2)}$$

其中：

• $U$：$L\times L$酉矩阵

• $\Sigma$：$L\times K$对角阵

• $V$：$K\times K$酉矩阵

奇异值分解将轨迹矩阵$X$分解为酉矩阵$U$，对角阵$\Sigma$和酉矩阵$V$的线性组合。这意味着：
$$\begin{array}{rl}\mathbf{X}& =\sum _{i=1}^r{\sigma }_i{U}_i{V}_i^{\text{T}}\\ & =\sum _{i=1}^r{\mathbf{X}}_i\end{array}\text{\hspace{1em}(3)}$$

$r$表示矩阵$X$的非零特征根数也即矩阵的秩。

序列重构

特征分组

不妨设根据某一分组原则将子矩阵分为了trend、periodic、noise 3组，对应地矩阵$X$分解得到的子序列将被组合为3部分。
$$\begin{array}{rlrl}{\tilde{\mathbf{X}}}^{\text{(trend)}}& =\sum _{t\in trend}{\tilde{\mathbf{X}}}_t& ⟹& {\tilde{F}}^{\text{(trend)}}=\sum _{t\in trend}{\tilde{F}}_t\\ {\tilde{\mathbf{X}}}^{\text{(periodic)}}& =\sum _{p\in periodic}{\tilde{\mathbf{X}}}_p& ⟹& {\tilde{F}}^{\text{(periodic)}}=\sum _{p\in periodic}{\tilde{F}}_p\\ {\tilde{\mathbf{X}}}^{\text{(noise)}}& =\sum _{n\in noise}{\tilde{\mathbf{X}}}_n& ⟹& {\tilde{F}}^{\text{(noise)}}=su{m}_{n\in noise}{\tilde{F}}_n\end{array}\text{\hspace{1em}(4)}$$

对角平均（化矩阵为序列）

$${\tilde{y}}_k=\{\begin{array}{lr}\frac{1}{k}{\sum }_{m=1}^k{x}_{m,k-m+1}& 1\le k<L^{\ast }\\ \frac{1}{L^{\ast }}{\sum }_{m=1}^{L^{\ast }}{x}_{m,k-m+1}& L^{\ast }\le k\le K^{\ast }\\ \frac{1}{N-k+1}{\sum }_{m=k-K^{\ast }+1}^{N-K^{\ast }+1}{x}_{m,k-m+1}& K^{\ast }<k\le N\end{array}\text{\hspace{1em}(5)}$$

其中，$L^{\ast }=min(L,N-L+1),K^{\ast }=max(L,N-L+1)$​

实例

数据来源于：Kaggle原油价格预测

'''数据处理'''
import pandas as pd
import numpy as np
'''数据可视化'''
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 800
#调整分辨率
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文
plt.rcParams['axes.unicode_minus']=False
#正常显示负号
###加载数据集###
file = r'D:/机器学习/crude-oil-price.csv'
oil_info = pd.read_csv(file,index_col='date')
price = oil_info['price']
# 算法封装
class SSA(object):
__supported_types = (pd.Series,np.ndarray,list) #限制时间序列的输入类型
def __init__(self,tseries,L):
'''
Args:
tseries:原始时间序列
L:窗口长度
'''
if not isinstance(tseries, self.__supported_types):
raise TypeError("请确保时间序列的数据类型为Pandas Series,NumpPy array 或者list")
else:
self.orig_TS = pd.Series(tseries)
self.N = len(tseries)
#原始时间序列长度
if not 2 <= L <= self.N / 2:
raise ValueError("窗口长度必须介于[2,N/2]")
self.L = L
#窗口长度，轨迹矩阵的行数
self.K = self.N - self.L + 1 #轨迹矩阵的列数
self.X = np.array([self.orig_TS.values[i:L+i] for i in range(0,self.K)]).T
#奇异值分解
self.U,self.Sigma,VH = np.linalg.svd(self.X)
self.r = np.linalg.matrix_rank(self.X)
#矩阵的秩等于非零特征值的数量
#每一个非零特征值都对应一个子矩阵，子矩阵对角平均化后得到原始时间序列的一个子序列
self.TS_comps = np.zeros((self.N,self.r))
#对角平均还原
for i in range(self.r):
X_elem = self.Sigma[i] * np.outer(self.U[:,i], VH[i,:])
X_rev = X_elem[::-1]
self.TS_comps[:,i] = [X_rev.diagonal(j).mean() for j in range(-X_rev.shape[0]+1, X_rev.shape[1])]
def comps_to_df(self):
'''将子序列数组转换成DataFrame类型
'''
cols = ["F{}".format(i) for i in range(self.r)]
return pd.DataFrame(data=self.TS_comps,columns=cols,index=self.orig_TS.index)
def reconsruct(self,indices):
'''重构，可以是部分重构（相当于子序列的分组合并），也可以是全部合并（重构为原序列）
Args:
indices
重构所选择的子序列
'''
if isinstance(indices,int):
indices = [indices]
ts_vals = self.TS_comps[:,indices].sum(axis=1)
return pd.Series(ts_vals,index=self.orig_TS.index)
def vis(self):
'''可视化子序列
'''
fig,axs = plt.subplots(self.r,sharex='all')
for i in range(self.r):
axs[i].plot(self.reconsruct(i),lw=1)
price_SSA = SSA(price,7)
comps_df = price_SSA.comps_to_df()

最后

以上就是顺心哑铃为你收集整理的SSA时间序列分解的全部内容，希望文章能够帮你解决SSA时间序列分解所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。