一種基于壓縮感知的時間序列預測方法

張亞寧，馬軍海

（天津大學 管理與經濟學部，天津 300072）

0 引言

對于時間序列數據來說已知數據和未知數據間通常都存在某種線性或非線性的關系，時間序列預測方法就是通過將它們間的這種關系估計出來而進行預測的。在現有的預測方法中以參數模型估計法最為流行，其中又以AR模型、ARMA和ARIMA模型的研究最為廣泛。這三種模型都假定時間序列是由某個白噪聲E={ei}，i=1，2…，N激勵而來的。參數模型估計法的任務就是在噪聲未知的前提下，估計出模型參數。以AR模型為例，AR模型的參數估計往往需要假定觀測數據為平穩時間序列，然后根據自相關函數建立起Yule-Walker方程(或 Wiener-Hopf方程)，再通過 Levinson遞推算法(或最小二乘)估計出參數[1]。但對于很多經濟數據來說，采用AR模型預測存在一些問題，如經濟數據往往呈現出遞增、遞減或周期性的變化趨勢，通常都不能滿足平穩時間序列均值和方差恒定的前提，再比如所能獲得的經濟數據量經常有限，而自相關函數往往需要大量的數據才能較為準確的估計出來。除此之外，p階AR模型只能考慮到距離預測值最近的p個已知數據的貢獻，而忽略其他數據。ARMA和ARIMA模型都以AR模型為基礎，因此也存在類似的問題。為解決這些問題，本文提出了基于壓縮感知[2]的參數估計方法，該方法并不局限于前p個觀測值對預測值的貢獻，而是根據數據的特點在一定范圍內找到最優的選擇；并且該方法對時間序列的平穩性沒有要求，只要求數據間確實存在某種可以表達的關系。

1 數學模型

若預測值為它的前期值和隨機項的線性函數，則時間序列數據預測模型可由式(1)描述：

其中，xn+1為預測數據，xi，i=1，2，…，n為前期值，ai為模型的待估計參數，en+1為隨機項，服從相互獨立的均值為0、方差為σ2的正態分布。如果每一組觀測數據，均可由同一組模型參數表達，則式(1)表達為矩陣形式：

為便于表達，本文將式（2）表達為矩陣相乘的形式：p=Ru+e。矩陣向量p、R、u和e的定義見式(2)，其中R是一個M+1行、N+1列矩陣。根據極大似然估計法，似然函數L為：

由于隨機項向量e未知，模型參數向量u無法通過直接解方程的方法求出。由于似然函數L取極大值等價于，因此模型參數向量u的可通過式(4)求解：

下面利用正交匹配追蹤法實現式(4)所示的優化問題的求解。

2 算法原理與實現

AR模型和ARMA模型在經濟數據預測中的成功應用表明，未知的經濟數據可通過有限個已知的經濟數據比較好的預測出來。據此，本文做如下假設：在式(2)中，當M和N都比較大時，向量為稀疏向量。即，向量u中非零元素的個數遠小于向量u中元素的總數，但這些非零元素的個數、大小和位置仍未知。至此，本節將時間序列預測問題歸結為：在u為稀疏向量的前提下求解式(4)的問題。壓縮感知理論是求解稀疏信號有力工具，并且已經在信號處理領域當中取得了巨大的成功，其基本思想是在向量u為稀疏向量的前提下，式(2)存在唯一解，且該解可求。

2.1 算法原理

本節選用壓縮感知理論中經典的正交匹配追蹤(Orthogonal Matching Pursuit,OMP)為基礎來求解上述問題。OMP算法屬于貪婪算法的范疇，它的基本思想是，在對問題求解時，總是做出在當前看來是最好的選擇，其數學原理見文獻[2]。本文在OMP算法的基礎上，針對式(4),提出的參數模型求解流程如下：

（1）初始化殘差r=p、迭代次序k=0、支持集Λ為空集。

（2）計算殘差r與矩陣R各列的相關系數，并找出相關系數最大的列對應的列標ik。

（3）首先將列標ik存入支持集中Λ=Λ∪{ik}，再將矩陣R中列標不屬于支持集Λ的元素置為零得到RΛ，最后在新的系統矩陣RΛ上根據（5）式估計出最小二乘解ut：

（4）更新殘差 r=p-RΛut。

（5）判斷終止條件是否成立。成立，則轉步驟6。不成立，迭代次序自增1：k=k+1，并轉步驟2。

（7）選擇置信度η，判斷去均值后的殘差數據r在置信度為η的條件下是否服從0均值方差為的正態分布。是，則輸出模型參數向量ut；否，選擇其他參數重新計算。

顯然，OMP算法的基本思想是從矩陣R中選擇與預測向量p最接近的列來逼近預測向量p并求出殘差r，然后再從投影矩陣R中選擇其他列向量進一步消除殘差。根據式(5)可得：

可見，更新后觀測矩陣RΛ的列向量始終與殘差r=p-RΛut保持正交。殘差r的能量必定不包含觀測矩陣RΛ列向量的貢獻。隨著矩陣RΛ維數不斷的擴大，殘差r必然不斷衰減。可見，該方法最大的特點是在殘差最小的前提下，自適應的選擇參數向量u中的非零元素，而不是像AR模型那樣，事先指定非零元素的個數和位置。步驟5中的終止條件有多種選擇方式，本文僅給出兩種選擇方式作為參考：

方式2：多次實驗確定最佳終止條件，k≥P時停止，P通過多次實驗確定。

具體選擇何種方式可根據實際情況自行選擇。除此之外，影響預測結果的還有矩陣R的行數M+1和列數N+1。其中，N主要影響非零元素的搜索范圍，因此N應當在條件允許的情況下，選擇一個足夠大的數值。參數M主要影響方程組的個數，如果M選擇過小，則模型參數較易受誤差的影響；如果M選擇過大，則同樣可能導致模型參數的不精確。這主要是因為式(2)暗含一個假設，即所有方程均滿足同一組模型參數u，如果M選擇過大，這個假設可能不成立。在實際應用中，參數M應當根據經驗在一定范圍內逐次嘗試選出最優值。

2.2 預測

為了驗證本文方法的有效性，筆者在軟件MATLAB(2010b)上實現了上述方法，并應用于對上證指數和我國第三產業值的預測。第一組實驗選用1994年3月1日到2013年2月28日的日收盤價對上證指數進行預測。該實驗主要用于說明本文算法的執行過程。為了給予參數向量u較大的選取空間，本文式（2）中的參數M和N分別選為49和199，對上證指數進行預測。預測結果與真實值的對比見圖1（b）。關于步驟5的終止條件，本文選擇第二種方式來確定。具體來說就是，針對觀測到的數據進行預測，并逐次嘗試迭代次數k的所有可能取值，選定預測值平均絕對百分比（Mean Absolute Percentage Error,MAPE）最低的迭代次數來進行預測。MAPE在每次迭代后的值見圖1.（a）。可見，，隨著k的增加，MAPE則表現為先下降后上升，在k=10時最小，因此對上證指數數據的預測中，本文推薦將步驟5中的迭代終止條件設定為迭代次序大于等于10。據此，對上證指數中的序列進行預測，預測結果與真實值見圖2（b），此時的平均絕對百分比為2.00%。經檢驗，殘差通過了顯著性水平為0.05的高斯分布檢驗。

圖1 對上證指數的預測

第二組實驗選取我國第三產業產值1952～2011年的數據來進行預測。本文式（2）中的參數M取為6，N取為7。關于步驟5的終止條件，本文選擇第二種方式來確定。具體來說就是逐次嘗試迭代次數k的所有可能取值，選定預測值平均絕對百分比（Mean Absolute Percentage Error,MAPE）最低的迭代次數來進行預測。由圖2（a）可知，隨著迭代次數k的增加，MAPE與之呈現出正相關的關系，因此步驟5的終止條件為P=1，此時MAPE=0.0573，殘差通過了顯著性水平為0.05的高斯分布檢驗，預測結果與真實值見圖2（b）。

圖2 對我國第三產業產值的預測

在2007～2011年數據未知的情況下，表1給出了對我國第三產業產值樣本數據進行預測得到的預測值與真實值的對比，經檢驗，殘差同樣也通過了顯著性水平為0.05的高斯分布檢驗。另外還給出了本方法與ARMA（1,1）模型利用同樣樣本數據得到的2007～2011年的預測值和真實值及誤差百分比。通過以上對比可知本文方法在預測精度上是優于ARMA模型的。最后，用本文方法對2012～2016年的我國第三產業產值進行預測，預測結果見表2。

表1 本文方法與ARMA模型預測結果對比 （單位：億元）

表2 對2012年到2016年我國第三產業產值的預測 （單位：億元）

3 結論

本文提出了一種基于壓縮感知的時間序列預測方法，首先根據極大似然估計理論建立了數學模型，再在壓縮感知的理論框架下將參數估計歸結為一個壓縮感知中的經典問題，進而采用正交匹配追蹤法求解出參數，最后得到時間序列的預測值。利用該方法對我國第三產業產值的預測結果表明，本文方法與ARMA模型相比具有更高的預測精度，預測結果有很高的實用參考價值。將本文方法和AR模型做對比不難發現，如果將向量u的非零元素設定在前p個，則該方法將退化為p階AR模型。也就是說，AR模型是本文方法的一個特例。本文方法的預測精度高于ARMA模型，另外由于對數據的平穩性，季節性等沒有要求，因此也有著更廣的適用范圍。

[1]李子奈,葉阿忠.高級計量經濟學[M].北京：清華大學出版社,2003.

[2]Tropp J,Gilbert A.Signal Recovery from Random Measurements Via Orthogonal Matching Pursuit[J].IEEE Transactions on Information Theory,2007,53(12).