缺失數(shù)據(jù)下滑動(dòng)平均模型的估計(jì)方法

陳博, 何幼樺

(上海大學(xué)理學(xué)院, 上海200444)

在對(duì)時(shí)間序列數(shù)據(jù)預(yù)處理的過程中,數(shù)據(jù)存在缺失值的情況屢見不鮮.這種情況的發(fā)生不僅會(huì)造成預(yù)處理時(shí)間的增加,還將嚴(yán)重影響數(shù)據(jù)質(zhì)量.此外,不同缺失數(shù)據(jù)的處理方法也會(huì)導(dǎo)致估計(jì)結(jié)果產(chǎn)生很大差異[1-2].因此,研究如何有效地處理缺失數(shù)據(jù)顯得十分必要.

對(duì)于含有缺失數(shù)據(jù)的時(shí)間序列的處理方法,大致可分為兩大類.第一類為借補(bǔ)法,即使用某種規(guī)則或方法來填補(bǔ)缺失數(shù)據(jù),如取平均值、線性插值法等.該方法易于實(shí)施,但沒有充分利用數(shù)據(jù)中一些關(guān)聯(lián)信息且不依賴于模型結(jié)構(gòu).基于此,將缺失數(shù)據(jù)與模型相結(jié)合的第二類方法[3]被提出.

龐新生[4]介紹了多重插補(bǔ)法的基礎(chǔ)理論與處理缺失數(shù)據(jù)時(shí)的基本思想.Junger等[5]將多重插補(bǔ)法應(yīng)用于空氣污染的時(shí)間序列中,并在缺失數(shù)據(jù)小于5%的情況下,取得了令人滿意的結(jié)果.Machado等[6]簡(jiǎn)述了存在數(shù)據(jù)缺失時(shí)自回歸模型的參數(shù)估計(jì)問題.Penzer等[7]運(yùn)用Kalman遞歸估計(jì)方法完成了對(duì)缺失數(shù)據(jù)下ARMA模型的似然方程擬合.仝倩等[8]通過經(jīng)驗(yàn)似然方法研究了含有缺失數(shù)據(jù)的半?yún)?shù)非線性模型的統(tǒng)計(jì)診斷問題.

除此之外,基于極大似然估計(jì)基礎(chǔ)上的EM算法[9]也是一種重要的處理缺失數(shù)據(jù)的統(tǒng)計(jì)方法.該算法通過期望(E)步來得到一個(gè)平均意義下的似然函數(shù),再通過最大化(M)步求得該次迭代似然函數(shù)的參數(shù)解,依此反復(fù)迭代得到最終的最優(yōu)估計(jì).Shumway等[10]為EM算法在時(shí)間序列上的應(yīng)用打下了基礎(chǔ).田萍等[11-12]應(yīng)用EM算法,針對(duì)含一個(gè)或連續(xù)兩個(gè)缺失數(shù)據(jù)的AR(p)模型給出了具體計(jì)算步驟,并針對(duì)含1～2個(gè)缺失數(shù)據(jù)的ARMA(1,1)模型做出了參數(shù)估計(jì),給出了其解析表達(dá)式.黃翔等[13]針對(duì)AR(p)模型的非左端缺失問題,且在缺失數(shù)據(jù)量較大時(shí),通過EM算法提出了有效的缺失值插補(bǔ)方案.

本工作針對(duì)存在連續(xù)缺失值的MA(q)模型,綜合了EM算法、極大似然原理及滑動(dòng)平均(moving average,MA)模型的相關(guān)理論,得到了模型參數(shù)的合理估計(jì),并在此基礎(chǔ)上進(jìn)一步完成了對(duì)缺失數(shù)據(jù)的填補(bǔ).通過數(shù)值模擬,研究了估計(jì)的有效性和精度,并將該算法應(yīng)用于實(shí)際數(shù)據(jù),得到了較好的估計(jì)結(jié)果.

1 模型參數(shù)估計(jì)

對(duì)于樣本序列X=(x1,x2,···,xn)′,滿足q階滑動(dòng)平均模型MA(q):

n?q,{εt}iid～N(0,σ2),可得到方程組

記

則原方程組可寫為

未知參數(shù)φ=(θ1,θ2,···,θq,σ2)的似然函數(shù)為

若樣本序列X中的(xs+1,xs+2,···,xs+r)為缺失數(shù)據(jù),記為

式中,T11,T22,T33分別為s,r和n-s-r階方陣,對(duì)數(shù)似然函數(shù)l為

估計(jì)參數(shù)φ的EM算法步驟如下.

(1)E步.在獲得第i步參數(shù)估計(jì)φ=φ(i)的條件下,對(duì)式(1)求關(guān)于Z2的條件期望:

式中,Rφ(i)=E(Z2Z′2)為缺失數(shù)據(jù)Z2的協(xié)方差矩陣.

(2)M步.極大化式(2),即第i+1步參數(shù)估計(jì)為

在給定參數(shù)φ的初始值φ(0)后,便可通過EM算法重復(fù)迭代,直至||φ(i+1)-φ(i)||小于給定精度τ,以此保證迭代的收斂性[9].

2 缺失數(shù)據(jù)估計(jì)

在給定參數(shù)φ下,缺失數(shù)據(jù)Z2的估計(jì)Z2有如下結(jié)論.

定理1 若含有缺失值的序列X為滑動(dòng)平均模型MA(q)的一組樣本,則在二次損失下缺失數(shù)據(jù)Z2的估計(jì)為

證明 考慮如下約束優(yōu)化問題:

該問題的拉格朗日乘子法目標(biāo)函數(shù)可寫為

式中,λ=(λ1,λ2,···,λn)′.計(jì)算

代入約束條件,可得

得到ε的最優(yōu)點(diǎn)為

最后極小化ε?′ε?,可得

式中,S11,S31,S13,S33分別為s×s,(n-s-r)×s,s×(n-s-r)和(n-s-r)×(n-s-r)階矩陣,則Z2的協(xié)方差矩陣為

推論1 當(dāng)缺失值在序列尾端,即s+r=n,則

證明 當(dāng)s+r=n時(shí),

因T與S/σ2互逆,由文獻(xiàn)[14]得

代入式(3),可得

推論2 當(dāng)缺失值在序列首端,即s=0,則

證明 同推論1.

而當(dāng)缺失值在非序列首尾部時(shí),類似地亦可完全用Tij直接表示出E(Z2Z′2),只是形式過于復(fù)雜,在此不再贅述.

3 數(shù)值模擬

設(shè)定序列長(zhǎng)度為50,σ=1,缺失位置在序列中部,參數(shù)初始值設(shè)置為θ(0)=0,EM算法精度τ=0.000 1,以均方誤差(mean square error,MSE)為評(píng)判標(biāo)準(zhǔn).本工作將通過數(shù)值模擬比較在相同情況下本算法與文獻(xiàn)[12]算法的參數(shù)估計(jì)效果,結(jié)果如表1所示.

表1 算法誤差對(duì)比Table 1 Algorithm error comparisons

由表1可知,本算法的精度明顯高于文獻(xiàn)[12],且文獻(xiàn)[12]算法只能解決單個(gè)或連續(xù)兩個(gè)數(shù)據(jù)缺失問題,也無法判斷MA系數(shù)的正負(fù)性,而本算法則沒有這些限制.

進(jìn)一步地,本工作將通過數(shù)值模擬研究樣本缺失比例、模型階數(shù)、序列長(zhǎng)度及模型特征根模長(zhǎng)等因素對(duì)參數(shù)估計(jì)整體均方誤差(MSEP)和缺失值估計(jì)整體均方差(MSEE)所造成的影響.記

3.1 模型特征根模長(zhǎng)、模型階數(shù)對(duì)整體均方誤差的影響

取序列長(zhǎng)度為50,樣本缺失比例為10%,觀察在相同序列長(zhǎng)度、樣本缺失比例情況下,模型特征根模長(zhǎng)與模型階數(shù)對(duì)估計(jì)整體均方誤差的影響,結(jié)果如表2,3所示.

表2 模型特征根模長(zhǎng)、模型階數(shù)與MSE P的關(guān)系Table 2 Relationships between length of model characteristic roots,order and MSE P

表3 模型特征根模長(zhǎng)、模型階數(shù)與MSE E的關(guān)系Table 3 Relationships between length of model characteristic roots,order and MSE E

由表2,3可得:當(dāng)序列長(zhǎng)度、模型特征根模長(zhǎng)、樣本缺失比例相同時(shí),MSEP與MSEE均與模型階數(shù)正相關(guān),即二者隨著模型階數(shù)的增加而增加;當(dāng)序列長(zhǎng)度、模型階數(shù)、樣本缺失比例相同時(shí),MSEP與MSEE同樣均與模型特征根模長(zhǎng)正相關(guān),即二者隨著模型特征根模長(zhǎng)的增加而增加.模型特征根模長(zhǎng)的大小反映了模型的平穩(wěn)程度.數(shù)值模擬結(jié)果顯示,模型越平穩(wěn),估計(jì)效果就越好.

3.2 序列長(zhǎng)度、樣本缺失比例對(duì)整體均方誤差的影響

設(shè)定模型為模型特征根模長(zhǎng)為0.1的MA模型,觀察在相同模型階數(shù)、模型特征根模長(zhǎng)的情況下,序列長(zhǎng)度與樣本缺失比例對(duì)估計(jì)整體均方誤差的影響,結(jié)果如表4,5所示.

表4 3階模型中序列長(zhǎng)度、樣本缺失比例與MSE P的關(guān)系Table 4 Relationships between length of sequence,proportion of missing data and MSE P in MA(3)

表5 3階模型中序列長(zhǎng)度、樣本缺失比例與MSE E的關(guān)系Table 5 Relationships between length of sequence,proportions of missing data and MSE E in MA(3)

本工作對(duì)不同模型階數(shù)、不同序列長(zhǎng)度和不同樣本缺失比例的情況進(jìn)行了數(shù)值模擬(表4,5僅為MA(3)的模擬結(jié)果).結(jié)果顯示:當(dāng)序列長(zhǎng)度、模型特征根模長(zhǎng)和模型階數(shù)保持固定時(shí),MSEP與樣本缺失比例正相關(guān),即MSEP隨著樣本缺失比例的增加而增加;當(dāng)模型特征根模長(zhǎng)、模型階數(shù)和樣本缺失比例相同時(shí),MSEP與序列長(zhǎng)度負(fù)相關(guān),即MSEP隨著的序列長(zhǎng)度的增加而減少.序列長(zhǎng)度及樣本缺失比例對(duì)MSEP均有明顯影響,MSEE但對(duì)此并不敏感.

4 實(shí)例計(jì)算

對(duì)2015年4月—2015年8月(共90日)的浦發(fā)銀行股票日對(duì)數(shù)收益率數(shù)據(jù)進(jìn)行建模,通過觀察序列的自相關(guān)函數(shù)圖可得序列的自相關(guān)函數(shù)是2階的.因此,在MA(2)模型的假設(shè)下,其參數(shù)的最大似然估計(jì)為θ1=-0.026 1,θ2=0.299 8,σ2=0.000 8,并在5%的顯著性水平下,接受該模型是二階滑動(dòng)平均模型的假設(shè).

假設(shè)因某種原因,序列中的10%數(shù)據(jù)缺失,缺失數(shù)據(jù)起始位置位于2015年6月3日,連續(xù)缺失9日數(shù)據(jù).通過本算法,在精度τ=0.000 1下,經(jīng)過10步迭代得到模型的參數(shù)估計(jì)值為θ1=-0.018 3,θ2=0.294 1,σ2=0.000 9,與完整數(shù)據(jù)下的結(jié)果基本相同.記缺失數(shù)據(jù)估計(jì)Z=(z1,z2,···,z9),缺失數(shù)據(jù)估計(jì)的方差Var(Z)=(0.000 078,0.000 077,0,0,0,0,0,0.000 077,0.000 078),相關(guān)系數(shù)矩陣為

計(jì)算結(jié)果表明,在缺失數(shù)據(jù)下的模型參數(shù)估計(jì)與完整數(shù)據(jù)的結(jié)果之間僅存在很小的誤差,證明本算法對(duì)參數(shù)估計(jì)是十分有效的.但是缺失數(shù)據(jù)估計(jì)的標(biāo)準(zhǔn)差與序列狀態(tài)近似處于同一數(shù)量級(jí)上,說明缺失數(shù)據(jù)估計(jì)效果不如參數(shù)估計(jì).由于當(dāng)MA(q)模型的預(yù)測(cè)點(diǎn)與樣本之間的距離大于q時(shí),其最佳線性預(yù)報(bào)值為模型的均值,因此實(shí)例中缺失數(shù)據(jù)估計(jì)的方差在與樣本距離大于2時(shí)均為0.同理,對(duì)于缺失數(shù)據(jù)估計(jì)的相關(guān)系數(shù)矩陣而言,當(dāng)缺失數(shù)據(jù)與左側(cè)樣本之間的距離小于等于q時(shí),則該類缺失數(shù)據(jù)估計(jì)之間存在著相關(guān)關(guān)系,而與其他缺失數(shù)據(jù)估計(jì)不相關(guān).類似地,當(dāng)缺失數(shù)據(jù)與右側(cè)樣本之間的距離小于等于q時(shí)亦是如此.