非參數模型的穩健跳點檢測估計

韓忠成林金官

(1.東南大學數學學院,江蘇 南京211189;2.南京審計大學統計與數學學院,江蘇 南京211815)

1.引言

觀測數據的曲線擬合具有廣闊的應用前景,非參數回歸模型為曲線擬合問題提供了一個主流的統計工具,其形式為

其中,m(·)是未知回歸函數有有界支撐U= [a,b],X是一維解釋變量,ε是獨立同分布的隨機誤差項.在某些情況下,回歸函數可能在某些未知位置存在跳點,表示相關過程的結構變化.比如,當生產線失控時,產品的質量指標可能在未知的時間點發生向下或向上的移動.在這種情況下,跳點的檢測對回歸函數結構的刻畫十分重要.

近年來,非參數模型跳點的估計已被廣泛研究.文[1]指出回歸函數可能存在不連續點,傳統光滑方法得到的擬合曲線在跳點處存在較大偏差.在跳點個數已知的假設下,文[2]提出了跳點和回歸函數的核估計方法.文[3]利用小波方法給出了跳點的檢測方法.文[4]利用回歸函數的單邊非參數回歸方法估計不連續點的位置.文[5]基于局部線性估計量構造跳點估計過程,證明了跳點估計過程的收斂性質.文[6]探討了不同方法下跳點估計問題的最優表現.在實際問題中,跳點的個數和位置通常是未知的.文[7]提出了一種不連續點的檢測方法.該方法通過比較任意給定點的三種估計量確定不連續點的位置.文[8]基于局部線性估計量提出了一種保跳曲線擬合方法.文[9]指出局部線性估計量不可避免地存在巨大的計算負擔,而B樣條在擬合不連續回歸函數時表現更好.

上述文獻的結果都是在最小二乘方法下得到的.然而,最小二乘方法對觀測數據存在異常點或重尾分布的情形十分敏感.眾所周知,M-估計常用來處理異常點的情形(見文[10]),但是當誤差項服從正態分布時M-估計會損失一些效率.因此,當帶跳非參數模型存在異常點時,需要發展一種合適的估計方法能同時獲得穩健性和有效性.但是,據知,目前還未有此類研究文獻出現.本文在跳點個數和位置未知的假設下,結合B樣條提出一個穩健有效的跳檢測方法,通過引入一個調節參數,改善回歸函數的估計效率.蒙特卡洛模擬和實例分析說明了提出的估計方法不僅在回歸函數的連續區間而且在跳點的鄰域內都有很好的表現.

本文結構如下: 第2節介紹估計方法;第3節通過數值模擬給出提出的方法在有限樣本下的表現;第4節用本文提出的方法處理上證指數數據.

2.估計方法

假設模型(1.1)中的回歸函數m(·)有如下表達式:

其中,g(x)是一元光滑函數,I(·)是示性函數當條件為真時取1,否則取0.q表示回歸函數中跳點的個數,dj和sj分別表示第j個跳點的幅度和位置.稱滿足式(2.1)的模型(1.1)為帶跳非參數模型.

Ⅰ眾數估計

假設{(Xi,Yi),i= 1,··· ,n}是來自模型(1.1)的一組樣本.為避免局部多項式估計的缺點,回歸函數m(x)可通過B樣條近似給出.令U= (u1,··· ,uK)表示支撐[a,b]上的內節點向量,對應的擴展節點向量記為則

其中,B(x) = (B1,p(x),··· ,BK+p+1,p(x))表示p階B樣條基函數,K表示內節點個數.根據眾數光滑思想,我們可通過最大化下式

估計α,其中,?h(t) =h?1?(t/h),h是需要選擇的帶寬,?(t)表示對稱核密度函數.?(t)的選擇不是非常嚴格,為了便于計算,本文?(t)取標準正態密度.

注意到最大化式(2.3)無法直接得到α的顯式解.為了估計α,給出如下的EM算法:

步0 計算α的初始值α(0).設置k=0.

步1 更新π(j|α(k)):

步2 更新α(k+1):

其中,MT= (B(X1),··· ,B(Xn)),Wk是以π(j|α(k))為元素的對角陣,Y= (Y1,··· ,Yn)T.設置k=k+1,并返回至步1.

步3 重復步1至步2,直到收斂.α的最終估計量,記作.回歸函數在點x處的估計量記為(x,U?)=B(x)T.

進一步,如果在U?內加入p+1個同樣的新節點x0∈(a,b),不失一般性,假設x0∈(ui,ui+1),則新的節點向量記為,即

和

類似(2.3)式,(2.5)式和(2.6)式的最優解可通過同樣的算法步驟獲得,分別記為和.則回歸函數在點x處的估計量記為令RSS0表示殘差平方和,即插入新節點之后的殘差平方和包含兩部分

和

注1步0中α(0)的計算可參見文[9]的方法.

Ⅱ跳點檢測估計

由文[11]可知,如果回歸函數m(x)在支撐[a,b]上是光滑的,則每個設計點(x;U?)是m(x)的相合估計;如果m(x)在支撐[a,b]上存在跳點,那么在跳點的鄰域內(x;U?)不是m(x)的相合估計.(x;)在區間[a,x0)和[x0,b]上也具有相同的性質.因此,為了提高回歸函數的估計精度,需要檢測觀測數據中的跳點.

為了檢測跳點,回歸函數估計量的距離函數定義如下:

直觀來說,若x0位于回歸曲線的連續區域在區間[a,x0) 和[x0,b]上與(x;U?)相差無幾,包括在跳點的鄰域內也是如此,所以接近很小;若x0位于跳點的鄰域內,僅在x0的左鄰域內相合,在x0的右鄰域內非相合,而在跳點兩側均是不相合的,因此,當x0接近跳點時,的差異十分顯著,D(x0)相應增加.特別地,如果x0與跳點重合,D(x0)可得到局部極大值點.

總體來說,當x0的鄰域內存在跳點,D(x0)變大且存在一個局部極大值點,否則D(x0)的值很小.根據D(x0)在跳點處的信息,我們提出下面的跳點檢測步驟:

第1 步: 對任一點x0,若滿足|D(x0)|≥?n,其中?n是非負閾值,則x0被標記為跳點.

第2 步: 假設{νi,i=1,··· ,q}是第一步檢測的跳點,且?n=Xi ?Xi?1均相等.若存在整數1≤i1?n且νi2+1?νi2>?n,則稱(νi1+νi2)/2是一個估計的跳點,代替跳點{νi1,··· ,νi2}.

利用上述程序可檢測出回歸函數中跳點的位置和個數,記作{ν?1,··· ,ν?q?}和q?.令ν?0=a,ν?q?+1=b,V={ν?0,··· ,ν?q?+1},不難發現,回歸函數在區間[ν?0,ν?1),··· ,[ν?q?,ν?q?+1]上是連續的.記新的節點向量為可通過最大化下式

進行估計,其中B?(x)是節點向量下的B樣條基函數向量.與(2.3)式類似,回歸函數在點x處的估計量為稱為穩健跳點檢測估計量.

Ⅲ參數選擇

在利用B樣條函數擬合回歸函數的過程中,有四個參數需要選擇: 內節點個數K,基函數階數p,帶寬h和閾值?n.首先討論參數K和p的選擇,通?？紤]以下二維交叉驗證準則

獲得.其次,由文[12]可知,基于B樣條函數的局部眾數估計量與最小二乘估計量的漸近方差之比如下所示:

其中σ2= E(ε2),F(h) = E(?′′h(ε)),G(h) = E(?′h(ε)2).比值R(h)僅依賴帶寬h,且在估計量的有效性和穩健性方面扮演著重要角色.因此,帶寬h的理想選擇為

由(2.7)式可知,hopt與樣本大小n無關,只與ε的條件誤差分布有關.

實際問題中,隨機誤差項的分布是未知的,因此F(h)和G(h)無法直接獲得.一個靈活的處理方法是通過

分別估計F(h)和G(h).則R(h)可利用來估計,其中表示基于初始估計得到的殘差項.利用格點搜索方法,很容易找到hopt最小化(h).

參數?n的選擇需要合適的跳點檢測準則,常用的評價準則為Hausdorff距離

其中J和分別表示真實的和估計的跳點集合.由于J未知,無法直接計算,故采用bootstrap方法.假設存在B個bootstrap樣本,根據第k個樣本檢測到的跳點記為則的估計為

?n的最優值可通過最小化獲得.

注2參數選擇的其他方法可參見文[9,11-12].

3.數值模擬

本節通過數值例子評價提出的跳點檢測方法和回歸函數估計量的有限樣本表現.考慮一組觀測值{(Xi,Yi),i=1,··· ,n}來自模型

其中Xi是來自[0,1]的均勻分布,回歸函數表達式如下函數m(x)有兩個跳點,分別位于0.3和0.7處,幅度分別是2.8和1.7.樣本量取n= 200和400,每次實驗重復N=200次.誤差分布考慮以下兩種不同情形:

情形1εi ～N(0,0.12),正態分布.

情形2εi ～0.95N(0,0.12)+0.05N(0,32),5%的數據可近似看作異常點.

首先,研究跳點檢測方法檢測跳點的能力.表3.1給出了不同情形下檢測到的跳點出現在真實跳點0.02范圍內的次數.與情形2相比,情形1中的跳點檢測方法的表現明顯更好.這一現象表明誤差分布的噪聲水平較小,跳點檢測方法的表現越好.進一步地,在情形2中,樣本量增加相應地提高了跳檢測方法檢測跳點的能力.同時,當跳點的幅度增加時有類似的結論.

表3.1 200次重復實驗下真實跳點0.02范圍內檢測出跳點的次數

其次,研究回歸函數估計量的有限樣本表現.在獲得跳點個數和位置的估計之后,使用提出的跳點檢測方法和眾數回歸方法(MPS)估計回歸函數.為了說明其有效性與穩健性,我們將該方法與基于分段樣條擬合和最小二乘提出的跳點檢測(LSPS)估計方法[9]進行比較,兩個曲線估計量分別記作在200次重復實驗下,對這兩個估計量計算相應的平均積分平方誤差(mean integral squared error,MISE)和跳點附近的局部MISE的值,結果如表3.2所示.

表3.2 回歸函數的MISE 和跳點附近的局部MISE 的模擬結果

4.實證分析

股票市場作為國民經濟的晴雨表，受到政府和投資者的高度重視。由于股票市場充滿了不確定性、機遇和風險,因此,挖掘有效信息可以幫助投資者抓住機遇并規避風險.

股票價格指標是度量金融市場信息的有效工具,從統計學角度分析股票價格指標對獲取信息十分重要.作為示例,我們收集了一組上海證券綜合指數從2014年1月2日至2016年12月30日的日收盤價數據(見http://q.stock.sohu.com).這三年中,股票市場經歷了幾次危機,稱為中國股市動蕩.從圖4.1可知,動蕩起始于2015年6月15日,于2016年2月早期終止.三個暴跌點出現在2015年6月,2015年8月,2016年1月.然而,由于噪聲的影響,跳點位置和幅度均是未知的.因此,跳點檢測以及收盤價曲線擬合需要格外關注.值得注意的是,在分析數據之前,有必要對數據進行歸一化處理.

圖4.1 2014年1月2日至2016年12月30日上海證券綜合指數的日收盤價數據

圖4.2 2014年1月2日至2016年12月30日上海證券綜合指數的擬合曲線

圖4.3 Y200 =5000作為異常點時,2014年1月2日至2016年12月30日上海證券綜合指數的擬合曲線

根據第2節的跳點檢測方法,從圖4.2中可觀測到三個跳點,分別位于0.483,0.548和0.667(對應日期2015年6月15日,2015年8月21日和2016年1月4日).檢測出的跳點位置與三個暴跌點的位置十分接近.同時,圖4.2中的擬合曲線與真實數據的變化趨勢保持一致,進一步說明提出的跳點檢測估計方法在跳點附近和連續區域內的表現良好.

為了檢驗本文提出的方法對異常點是否穩健,將第200個觀測值設為Y200=5000,見圖4.3.不難發現,本文提出的方法與最小二乘法的跳點檢測方法的跳點檢測結果與圖4.2中的結果保持一致.當存在異常點的時候,基于最小二乘的跳點檢測方法的回歸函數估計量(虛線)在異常點附近明顯偏離了真實曲線.然而,基于眾數的跳點檢測方法的回歸估計量(虛點線)與圖4.2中的結果保持一致.因此,基于眾數的跳點檢測方法是穩健的.