基于多尺度小波分解的時間序列預測方法研究

摘要：運用Mallat算法和Daubechies小波分解技術，把時間序列分解為比原始時間序列更單一的細節部分和概貌部分，然后把分解后的細節部分和概貌部分重構回原尺度，對重構后的各個時間序列用傳統時間序列模型進行預測，由此建立高階AR模型，最后累加各個時間序列預測結果得到原始時間序列的預測結果。通過對某地區工業總產值數據的分析和驗證，表明AR-wavelets模型與傳統單一模型相比可大大提高精度。

關鍵詞：多分辨率分析；Mallat算法；Daubechies小波；ARMA(P，Q)模型；時間序列分析

0 引言

傳統的預測方法如自回歸模型(AR)，滑動平均模型(MA)，自回歸滑動平均模型(ARMA)是在時間序列平穩的假設下，建立線性模型，然后采用模型外推的方法預測未來值。因此，這些方法只適用于平穩時間序列的預測。但現實生活中的時間序列數據往往是高度非平穩的，傳統的時間序列預測方法無法取得高精度的預測結果。

小波分解與重構實質上是通過不同的帶通濾波器將含有綜合信息的一組原始信號V(k)分解成了N＋I組特征不同的時間序列信號。其中一組近似信號反映了該時間序列內在的變換趨勢，而N組細節信號反映的是隨機擾動帶來的影響，二者的規律是不同的，對特征不同的信號選擇不同的參數進行預測，這樣分別預測的結果再合成，效果會比整體做預測的精度高。

1 預測原理

1.1 小波分解與重構預測原理簡介

小波分解與重構原理及預測過程可用圖1描述。

若V₀代表原始信號的集合，把第i組截止到采樣周期為止得到的原始時間序列信號記為V^o_i(K)，上標表示分解尺度。這里把原始信號視為。尺度上的信號，顯然有V^o_i(k)∈V_o。對V^o_i(k)進行分解，利用小波分解公式N尺度分解，得到一組基本時間序列信號V^N_i(k)和N組干擾信號ω^j_i(k)，(j＝1，2，…，N)。對分解得到的N＋1組時間序列信號分別單獨用Mallat算法重構到原尺度上，得到N＋I組在原始尺度上的經過分解重構處理的時間序列序號V^N_i(k)和ω^j_i(k)，j＝1，2，…，N)。然后再分別對每一組時間序列信號用ARMA模型進行預測，得到N＋1個預測值V^N_i(k)和ω^j_i(k)，(j＝1，2，…，N)。最后依據式(1)就得到預測結果。

1.2 對分解后的時間序列信號進行ARMA模型預測

小波分析工具可以使時間序列數據信號變得平滑，但這并不意味著分解后的概貌信號和細節信號變得平穩，因此還需要對分解后的信號進行分析。為此，首先要對數據進行初步的整理和必要的檢驗，這些工作稱之為預處理。

(1)測量數據的均值、方差和概率直方圖分析

測量數據的均值、方差和概率直方圖分析可以直接調用MATLAB函數：mean.m、vaLm、histfit.m。

(2)提取趨勢項

在Matlab中提供了提取趨勢項的函數detrend.m和dtrend.m。使用它們可以方便地提取趨勢項。

(3)分解后時間序列平穩化方法

一般來說，時間序列的不平穩表現在以下幾個方面：均值不平穩、方差不平穩、或均值、方差都不平穩。可以采取以下措施使時間序列平穩化。

A．方差平穩化

一般的，為了使方差平穩化可以用指數變換。見(2)式。式中，s(λ)稱為殘差平方和，μ是相對應的樣本均值。在指數變換中的λ可以作為模型參數由觀測序列去估計。使殘差平方和最小的λ相對應的變換即所需要的變換。見表1。

B.均值平穩化

非平穩時間序列可以通過對原始序列取適當階數差分而化為平穩序列。一般的時間序列差分兩次就可以轉化成平穩時間序列。考慮到原始時間序列差分后可能產生負數，方差平穩化操作必須在均值平穩化操作之前。

(4)模型類型和階次辨識

Matlab中提供了自相關函數autocoorr()和偏相關函數parcorr()。自相關函數和偏相關函數用于時間序列的模型識別，其準則是：如果時間序列的自相關函數是拖尾的，而偏相關函數是截尾的，那么這個時間序列符合AR(P)模型；如果時間序列的自相關函數是截尾的，而偏相關函數是拖尾的，那么這個時間序列符合MA(Q)模型；如果自相關函數和偏相關函數都是拖尾的，那么這個時間序列符合ARMA(P，Q)模型。關于階次，本文中是利用樣本自相關函數和樣本偏相關函數來推斷的。

2 實驗

2.1 基本資料

分析序列采用某地區1998-2007年生產總產值(單位：萬元)的月度資料，共有120個數據，圖2是這120個數據的折線圖。我們對1998-2006年數據建模，2007年的數據留做檢驗模型。

從圖2可以看出時間序列具有明顯的增長趨勢，并且含有周期為12個月的季節波動，即序列是非平穩的。

2.2 原始時間序列的分解與重構

本文對1998-2006年數據利用Daubechies小波系N＝3，即db3作為尺度函數進行多尺度分析。為了便于比較，將分解得到的各個尺度上的概貌信號、細節信號分別重構回原尺度。圖3是分解后各尺度上的近似信號的重構結果。

2.2.1 概貌信號的特征分析

從前期工作的成果中已經知道原始信號具有年度周期特征。即具有12個月的周期特征。概貌信號，去除了干擾信號，繼承了原始信號的主要特征和趨勢，那么理論上概貌信號應該和原始信號有同樣的周期特征。

從圖3中看出重構后的概貌信號是不平穩的。根據前面介紹，對序列進行方差變換和零均值變換，并進行T_a3＝II或T_a3＝12的季節差分。然后再看平穩化的自相關函數和偏相關函數。如圖4、圖5所示。

在進行季節差分時，根據T_a3＝II時和T_a3＝12時的自相關函數圖像和偏相關函數圖像判斷，T_a3＝12時能很好地反映時間序列的信息。圖4和圖5就是T_a3＝12時的自相關函數和偏相關函數圖像。由此可以確定平穩化后的概貌信號符合ARMA模型。有限階的ARMA(P，Q)序列可以轉化為無限階的AR(P)序列。因此可以用高階AR(P)模型代替ARMA(P，Q)。

2.2.2 細節信號的特征分析

細節信號可以被看作干擾信號，但是低頻干擾信號也有可能會受到與原始信號同樣的周期性因素影響。對于第三層細節分量d3，是所有細節分量中頻率最低的，也是各細節分量中噪聲成分最少的，它最可能和概貌信號有同樣的時序特征。為了判斷細節信號的時序特征，本文首先計算將細節分量在原尺度上重構后的原序列的自相關系數，零均值化后的細節信號d3的自相關函數如圖6所示。

由圖6可以看出，未考慮周期性而使用原分量重構序列計算的自相關系數呈現明顯的周期性，周期T_a3和T_a3相等，即T_a3＝12。我們把該序列按照時段化分為12個子時間序列，考察不同年同月的子時間序列的自相關系數，可以發現自相關系數呈現緩慢下降趨勢。這個特征表明該子序列屬于不穩定序列，且屬高階AR，階數p≥2模型，可以使用高階AR模型進行預測。

同理，可對第二層細節信號d2和第一層細節信號d1進行如上分析，建立高階AR模型。

2.2.3 子時間序列AR模型曲建立

各分量均使用高階AR模型逼近，因此首先要確定模型階數，不同的分量子序列分別用不同階數的AR模型。高階AR模型的關鍵問題是參數的確定(未知參數個數的控制就顯得很重要)。本文利用遞推最小二乘方法確定模型的參數，且使用常用的AIC(Akaike Information Criterion)方法確定模型的階數。AIC的準則的一般形式見(4)式：(4)式中，，N為樣本大小，L為預先給定的模型最高階數。當準則函數取最小值時的模型為適用模型。

2.3 序列預測

適用AIC確定概貌信號和細節信號的最佳回歸階數的時間序列，可以利用已有的ARMA方法或All、MA方法進行建模和預測，最后再將各序列的預測結果合成為對原始信號的預測結果。這一建模預測過程是先重構后預測的方法。

與確定分解尺度相似，將V^N_k和W^j_k(j＝1，2，…，N)N＋1個序列單獨重構到原尺度上，得到N＋1個在原始尺度上的時間序列V^N_k和W^j_k(j＝1，2，…，N)，這N＋1個時間序列具有相同的樣本數。它們的代數和等于原信號，見(5)式。即可以在原尺度上對N＋1個重構后的時間序列分別進行建模，得到模型(6)、(7)：

3 結束語

利用上述(6)，(7)模型進行預測，并將預測結果累加，就得到了原信號的預測結果。表2即預測的結果。我們可以看到誤差均在5％內，可以說預測的效果比較好。