張亞寧，馬軍海

（天津大學(xué) 管理與經(jīng)濟(jì)學(xué)部，天津 300072）

0 引言

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

1 數(shù)學(xué)模型

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

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

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

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

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

2 算法原理與實現(xiàn)

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

2.1 算法原理

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

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

（2）計算殘差r與矩陣R各列的相關(guān)系數(shù)，并找出相關(guān)系數(shù)最大的列對應(yīng)的列標(biāo)ik。

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

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

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

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

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

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

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

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

2.2 預(yù)測

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

圖1 對上證指數(shù)的預(yù)測

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

圖2 對我國第三產(chǎn)業(yè)產(chǎn)值的預(yù)測

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

表1 本文方法與ARMA模型預(yù)測結(jié)果對比 （單位：億元）

表2 對2012年到2016年我國第三產(chǎn)業(yè)產(chǎn)值的預(yù)測 （單位：億元）

3 結(jié)論

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

[1]李子奈,葉阿忠.高級計量經(jīng)濟(jì)學(xué)[M].北京：清華大學(xué)出版社,2003.

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