如何將一個時間序列分解為周期序列和趨勢序列的和？-魔扣目錄

時間序列分解是時序分析中的重要方法，廣泛應用于時間序列預測，時間序列異常檢測，時間序列聚類等場景，在工業界有很多的落地應用。

一個時間序列往往是以下幾類變化形式的疊加或耦合：

長期趨勢(Secular trend, T)：長期趨勢指現象在較長時期內持續發展變化的一種趨向或狀態。
季節變動(Seasonal Variation, S)：季節波動是由于季節的變化引起的現象發展水平的規則變動
循環波動(Cyclical Variation, C)：循環波動指以若干年為期限，不具嚴格規則的周期性連續變動
不規則波動(Irregular Variation, I): 不規則波動指由于眾多偶然因素對時間序列造成的影響

其中循環波動和季節變動一般可整合分解為有規律的周期序列，而長期趨勢則可整合分解為趨勢序列。分解出這兩種序列是時序分解中的重要課題。

本次文章為大家整理分解周期序列和趨勢序列的方法。

概述

這里主要使用的技術是奇異譜分析（SSA)，其是根據觀測到的時間序列構造軌跡矩陣，并對軌跡矩陣進行分解和重構，從而提取出代表原時間序列不同成分的信號，如長期趨勢信號、周期信號、噪聲信號等，從而進一步對分解得到的信號進行分析。

算法流程如下：

根據原始時間序列構建軌跡矩陣 X
對矩陣X進行奇異值分解
按奇異值生成r個子矩陣
根據某一分組原則將子矩陣 Xi 分為 m個組
對子矩陣 Xi 進行對角均值化處理得到子序列!
對m個組中的子序列相加得到分組子序列。

矩陣分解

01 時間序列嵌入形成軌跡矩陣

輸入：原始時間序列y，窗口長度L

顯然矩陣X是一個漢克爾矩陣（每一條副對角線的元素都相等），矩陣的行數為窗口長度L，列數為N-L+1

02 奇異值分解

其中：

奇異值分解將軌跡矩陣 X 分解為酉矩陣 U，對角陣 ∑ 和酉矩陣V的線性組合。這意味著：

r表示矩陣X的非零特征根數也即矩陣的秩。

03 特征分組

不妨設根據某一分組原則將子矩陣分為了trend、periodic、noise 3組，對應的矩陣X分解得到的子序列將被組合為3部分。

04 對角平均化

其中：

代碼實例

'''數據處理'''
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()