一个时间序列往往是以下几类变化形式的叠加或耦合:
其中循环波动和季节变动一般可整合分解为有规律的周期序列,而长期趋势则可整合分解为趋势序列。分解出这两种序列是时序分解中的重要课题。
本次文章为大家整理分解周期序列和趋势序列的方法。
这里主要使用的技术是奇异谱分析(SSA),其是根据观测到的时间序列构造轨迹矩阵,并对轨迹矩阵进行分解和重构,从而提取出代表原时间序列不同成分的信号,如长期趋势信号、周期信号、噪声信号等,从而进一步对分解得到的信号进行分析。
算法流程如下:
输入:原始时间序列y,窗口长度L
显然矩阵X是一个汉克尔矩阵(每一条副对角线的元素都相等),矩阵的行数为窗口长度L,列数为N-L+1
其中:
奇异值分解将轨迹矩阵 X 分解为酉矩阵 U,对角阵 ∑ 和酉矩阵V的线性组合。这意味着:
r表示矩阵X的非零特征根数也即矩阵的秩。
不妨设根据某一分组原则将子矩阵分为了trend、periodic、noise 3组,对应的矩阵X分解得到的子序列将被组合为3部分。
其中:
'''数据处理'''
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()
数据来源:
https://www.kaggle.com/code/naveenkonam1985/crude-oil-price-prediction/data