面板數(shù)據(jù)下序貫?zāi)Ｐ妥凕c(diǎn)的漸近檢驗(yàn)法

陳卓恒,胡亦鈞

(1.華僑大學(xué)數(shù)學(xué)科學(xué)學(xué)院,福建泉州 362021)

(2.武漢大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,湖北武漢 430072)

1 引言

變點(diǎn)問題起源于工業(yè)質(zhì)量控制的領(lǐng)域,在質(zhì)量控制中非常重要的一點(diǎn)是如何快速檢測出生產(chǎn)線上不合格產(chǎn)品比例的增加.在建模和數(shù)據(jù)挖掘中,一個很常見的問題是,如何根據(jù)現(xiàn)有的樣本數(shù)據(jù),來判斷我們感興趣的一些量是否發(fā)生變化.這種問題在統(tǒng)計中我們稱其為變點(diǎn)問題.關(guān)于變點(diǎn)問題的建模和分析始于Page 1954年發(fā)表在《Biometrika》上的文獻(xiàn)[1],文獻(xiàn)主要考慮了利用分位數(shù)構(gòu)造的簡單檢測方法去檢驗(yàn)生產(chǎn)線上觀察到數(shù)據(jù)的穩(wěn)定性.在此之后,關(guān)于變點(diǎn)問題的研究變得日益活躍,各種研究文獻(xiàn)也日益增多.而且隨著學(xué)科的發(fā)展和深入,變點(diǎn)問題在經(jīng)濟(jì)、金融學(xué)、生物醫(yī)學(xué)、氣候?qū)W、導(dǎo)航系統(tǒng)、圖像處理、信號探測、計算機(jī)等很多領(lǐng)域都有廣泛的應(yīng)用背景.

很久以來大多數(shù)的變點(diǎn)檢測的文獻(xiàn)都是基于這樣的背景:對于一組固定容量的歷史數(shù)據(jù)集,我們設(shè)計一些檢驗(yàn)方法來判斷這組數(shù)據(jù)內(nèi)我們感興趣的參數(shù)是否發(fā)生變化.這種檢驗(yàn)即為我們常說的“事后檢驗(yàn)”.在這種檢驗(yàn)中,樣本是靜態(tài)的,這種數(shù)據(jù)我們稱為是離線數(shù)據(jù).離線數(shù)據(jù)分析主要在譬如歷史文本分析,圖像分析等領(lǐng)域有應(yīng)用.但是在更多的領(lǐng)域譬如質(zhì)量控制,醫(yī)療監(jiān)測,金融風(fēng)險控制,我們使用的都是在線數(shù)據(jù),也就是樣本采取連續(xù)抽樣的方式得到的.采取連續(xù)抽樣的方式來得到樣本進(jìn)行檢驗(yàn),這種檢驗(yàn)稱為“序貫檢驗(yàn)”.序貫檢驗(yàn)也是最近幾年變點(diǎn)檢測中的研究熱點(diǎn).具體來說,在序貫檢驗(yàn)中,一般假定靜止期的長度為m,也就是在m個觀察數(shù)據(jù)X1,X2,···,Xm內(nèi)不存在變點(diǎn),這m個觀察數(shù)據(jù)也常稱為歷史數(shù)據(jù).在考慮漸近性質(zhì)的時候一般令m→∞.在原假設(shè)成立即認(rèn)為沒有變點(diǎn)存在的情況下Xm+1,Xm+2,···的參數(shù)是相同的;若備則假設(shè)成立即認(rèn)為變點(diǎn)存在的情況下,則存在一個整數(shù) k?≥ 0,使得 Xm+1,Xm+2,···,Xm+k? 的參數(shù)與歷史數(shù)據(jù)相同,而 Xm+k?+1,Xm+k?+2,···的參數(shù)不同,這里k?為未知變點(diǎn).接下來要做的就是如何構(gòu)造合適的檢驗(yàn)統(tǒng)計量并制定合理的停止規(guī)則,從而判斷變點(diǎn)是否存在以及在認(rèn)為變點(diǎn)存在的情況下推斷出變點(diǎn)的位置,使得變點(diǎn)從產(chǎn)生到被檢測出需要的時間盡量的短,并盡量消除誤報,即快速準(zhǔn)確檢測出變點(diǎn).

文獻(xiàn)[2]是較早提出并研究變點(diǎn)檢測中序貫檢驗(yàn)的代表性文獻(xiàn)之一,文中提出了兩種檢測方案,一種利用波動差,一種利用殘差的累積和,并拓寬了經(jīng)典的不變原理得到所需要的收斂結(jié)果,最后確定了臨界值和停止規(guī)則,找到最優(yōu)停時,并在模擬中得到了檢驗(yàn).文獻(xiàn)[3]進(jìn)一步推廣了文獻(xiàn)[2]的結(jié)果,同樣對于線性回歸模型,文獻(xiàn)[3]提出了兩類檢驗(yàn)統(tǒng)計量,第一類統(tǒng)計量是建立在回歸系數(shù)β取最小二乘估計時殘差的加權(quán)累積和,第二類是利用回歸系數(shù)構(gòu)造的迭代殘差的部分和得到的統(tǒng)計量.并制定了合適的檢驗(yàn)規(guī)則,使得在歷史數(shù)據(jù)個數(shù)m→∞時錯誤預(yù)報率在約定的水平內(nèi),且檢驗(yàn)功效趨于1.學(xué)者Aue在變點(diǎn)的序貫檢驗(yàn)中也做了一系列工作,譬如率先開展研究停止時刻τn的極限分布,得到了文獻(xiàn)[4],后續(xù)工作有文獻(xiàn)[5–7],另外文獻(xiàn)[8]得到了在對ARMA時間序列中結(jié)構(gòu)變點(diǎn)進(jìn)行序貫檢驗(yàn)時,停止時刻的極限分布.在近幾年中,也有把自助抽樣技術(shù)引入到序貫?zāi)Ｐ椭械淖凕c(diǎn)檢測問題中,可參考文獻(xiàn)[9,10]和[11].

目前來說,現(xiàn)有的文獻(xiàn)主要研究單序列情形下的序貫?zāi)Ｐ?那么,如何在面板數(shù)據(jù)下研究序貫?zāi)Ｐ偷淖凕c(diǎn),這是一個嶄新的問題.在本文中,我們主要對這個問題進(jìn)行研究.本文剩余部分的結(jié)構(gòu)是這樣的:第2節(jié)將描述模型并得到相關(guān)的理論結(jié)果,即檢驗(yàn)統(tǒng)計量的大樣本性質(zhì);第3節(jié)基于前面的理論結(jié)果提出面板數(shù)據(jù)下關(guān)于變點(diǎn)的漸近檢驗(yàn)法,并進(jìn)行了Monte-Carlo模擬計算;第4節(jié)給出了前面引理和定理的證明過程;第5節(jié)對前面的結(jié)論進(jìn)行了總結(jié).

2 模型及相關(guān)理論結(jié)果

考慮面板均值變點(diǎn)模型

這里 μij和 εij分別表示均值和面板擾動項(xiàng).{μij,j=1,2,···,m+Tm}i和 {εij,j=1,2,···,m+Tm}i是相互獨(dú)立的.Tm與m 有關(guān)而且假定Tm≤ ∞,序列{εij}關(guān)于i相互獨(dú)立,εij～ (0,σ2)且在概率空間(?,F,P)上關(guān)于j是獨(dú)立同分布的(i.i.d.).假定0<σ2<∞而且對于某個υ>2有E|εij|υ<∞.

假定對于每個1≤i≤N,已經(jīng)觀察到長度為m的歷史數(shù)據(jù),它沒有變點(diǎn)發(fā)生,即對于1≤i≤N,

現(xiàn)在我們感興趣的是接下來到來的數(shù)據(jù)里是否存在均值發(fā)生變化的共同變點(diǎn),即想檢驗(yàn)原假設(shè)

和備則假設(shè)

這里σ2,μi0,δi0,k0m的值均為未知.檢測變點(diǎn)的方案是基于CUSUM型檢測統(tǒng)計量Γ(m,k,N)和一個邊界值,這里

此處

此外,基于歷史數(shù)據(jù)集 {Yi1,Yi2,···,Yim}i=1,2,···,N使用以下的方差估計

關(guān)于檢驗(yàn)法則,在

處停止并拒絕H0.τN(m)也被稱為停時,而且這里的c應(yīng)該被選擇成使我們能夠控制錯誤預(yù)報率.即在原假設(shè)H0下,對于某個給定水平0<α<1,

在備則假設(shè)H1下,要求

定理2.1對于面板數(shù)據(jù)模型(2.1),在原假設(shè)(i.e.(2.3)式)的假定下,當(dāng)m,N→∞時,有

這里{W(t),0≤t<∞}是一個Wiener過程,γ見(2.6)式.

3 漸近檢驗(yàn)法及其模擬計算

根據(jù)定理2.1知,在原假設(shè)成立即沒有變點(diǎn)時,若m,N→∞時,

的極限分布為

這里{W(t),0≤t<∞}是一個Wiener過程.因此可以通過此處已知的極限分布求出在給定的顯著性水平下檢驗(yàn)的臨界值c,然后利用這里的臨界值c以及檢驗(yàn)法則來確定是拒絕還是接受原假設(shè),并在原假設(shè)被拒絕,即認(rèn)為變點(diǎn)存在的情形下估計出變點(diǎn)的位置.這種方法我們稱它為“漸近檢驗(yàn)法”.

關(guān)于漸近檢驗(yàn)法,首先要計算檢驗(yàn)的臨界值.由于臨界值只與γ有關(guān),所以可記作cγ.為了計算經(jīng)驗(yàn)檢驗(yàn)水平,考慮一個不帶變點(diǎn)的模型,并且計算

的次數(shù),這樣就得到了這個次數(shù)和Monte-Carlo模擬次數(shù)之比.進(jìn)一步,為了計算經(jīng)驗(yàn)檢驗(yàn)功效,考慮帶一個變點(diǎn)的模型,并進(jìn)行M 次Monte-Carlo模擬.在這M 次模擬中,計算出檢驗(yàn)統(tǒng)計量

超過臨界值cγ的頻率,然后來計算經(jīng)驗(yàn)檢驗(yàn)功效.同時,|Γ(m,k,N)|超過臨界值cγ的最早出現(xiàn)的那個k被認(rèn)為是變點(diǎn)位置的估計值.

一般來說,作為一個好的檢驗(yàn)方法,經(jīng)驗(yàn)檢驗(yàn)水平應(yīng)該小于或等于給定的水平,而且經(jīng)驗(yàn)檢驗(yàn)功效應(yīng)該足夠大.另外,在變點(diǎn)后的停時應(yīng)該足夠的短才好.因此首先按照定理2.1來計算臨界值cγ.

第二步 在這M 次模擬的基礎(chǔ)上,計算臨界值cγ,使得P[Uγ>cγ]=α.

表1:由50000次Monte-Carlo模擬計算出的臨界值cγ

接下來,在表2–8中,利用表1中計算出的臨界值cγ,分別計算了在γ,Tm,N,m 取不同值,變點(diǎn)k0m在不同位置(左端,中間或右端)時,漸近方法的經(jīng)驗(yàn)檢驗(yàn)水平,檢驗(yàn)功效以及停止時刻的各數(shù)字特征.具體來說:

表2:利用cγ經(jīng)1000次Monte-Carlo模擬計算出的經(jīng)驗(yàn)檢驗(yàn)水平,N=200,Tm=500

表3:利用cγ經(jīng)1000次Monte-Carlo模擬計算出的經(jīng)驗(yàn)檢驗(yàn)水平,Tm=500,m=100

表4:利用cγ經(jīng)1000次Monte-Carlo模擬計算出的檢驗(yàn)功效,N=200,Tm=500,k0m=25

表5:利用cγ經(jīng)1000次Monte-Carlo模擬計算出的檢驗(yàn)功效,N=200,Tm=500,m=100

從表2可以看出,當(dāng)m固定時,隨著γ的增大,經(jīng)驗(yàn)檢驗(yàn)水平越小.當(dāng)γ固定時,隨著m的增大,經(jīng)驗(yàn)檢驗(yàn)水平越小.從第一行和第二行的數(shù)據(jù)可以看到,當(dāng)γ=0和0.25時,若歷史數(shù)據(jù)集m較少,得到的經(jīng)驗(yàn)檢驗(yàn)水平是會超過給定水平的.從而可以考慮在[0,1/2)的范圍內(nèi)選取較大的γ,來降低檢驗(yàn)犯第一類錯誤的概率.而m越大,即沒有變點(diǎn)的歷史數(shù)據(jù)越多,得到的經(jīng)驗(yàn)檢驗(yàn)水平越小,這與我們的直觀感覺也是相符合的.

從表3可以看出,γ的大小對經(jīng)驗(yàn)檢驗(yàn)水平影響很大,類似表2的表現(xiàn).當(dāng)N 固定時,隨著γ的增大,經(jīng)驗(yàn)檢驗(yàn)水平是降低的.在γ=0時,經(jīng)驗(yàn)檢驗(yàn)水平甚至超過了給定的水平.同時注意到,N的大小對經(jīng)驗(yàn)檢驗(yàn)水平的影響不大,不過比較而言,中等大的N似乎能產(chǎn)生相對較低的經(jīng)驗(yàn)檢驗(yàn)水平.

通過模擬得到的列出的,譬如表4,以及一些未列出的表格發(fā)現(xiàn),當(dāng)變點(diǎn)發(fā)生在m個經(jīng)驗(yàn)數(shù)據(jù)之后不久時,則無論m和N 值如何變化,對于不同的α和γ,檢驗(yàn)的功效都達(dá)到了1.

從表5可以看出,在變點(diǎn)k0m在m之后不久或者是更遠(yuǎn)一些,檢驗(yàn)的功效都是1.但是,在變點(diǎn)靠近右端點(diǎn)時,檢驗(yàn)的功效會降低,最高不超過0.35.在k0m固定時,隨著γ的增大,功效是越來越低的.因此,基于檢驗(yàn)水平和檢驗(yàn)功效的考慮,應(yīng)該選取不大不小的γ,譬如取γ=0.25,從而使兩者能達(dá)到一種平衡.

表6:N=200,Tm=500,=25,模擬次數(shù)1000次計算出的停止時刻的各數(shù)字特征

表6:N=200,Tm=500,=25,模擬次數(shù)1000次計算出的停止時刻的各數(shù)字特征

γ m=30 m=100 τm ↓;α → 0.025 0.05 0.10 0.025 0.05 0.10 0 Min 9.00 8.00 6.00 26.0 16.00 13.00 Q1 28.00 27.00 27.00 28.0 28.00 27.00 Median 28.00 28.00 28.00 28.0 28.00 28.00 Mean 28.44 28.06 27.58 28.5 28.22 27.92 Q3 29.00 29.00 29.00 29.0 29.00 29.00 Max 32.00 32.00 31.00 32.0 31.00 31.00 0.25 Min 5.00 4.00 4.00 26.00 26.00 11.00 Q1 28.00 28.00 27.00 28.00 28.00 27.00 Median 28.50 28.00 28.00 29.00 28.00 28.00 Mean 28.51 28.17 27.73 28.53 28.27 27.94 Q3 29.00 29.00 29.00 29.00 29.00 29.00 Max 33.00 32.00 31.00 32.00 31.00 31.00 0.45 Min 26.00 20.00 10.00 27.00 26.00 26.00 Q1 28.00 28.00 28.00 28.00 28.00 28.00 Median 29.00 29.00 29.00 29.00 29.00 29.00 Mean 29.21 28.84 28.5 29.12 28.83 28.55 Q3 30.00 30.00 29.00 30.00 29.00 29.00 Max 33.00 33.00 32.00 32.00 32.00 31.00 0.49 Min 27.00 26.00 19.00 27.00 26.00 26.00 Q1 29.00 28.00 28.00 29.00 28.00 28.00 Median 30.00 29.00 29.00 29.00 29.00 29.00 Mean 29.57 29.17 28.88 29.47 29.12 28.83 Q3 30.00 30.00 30.00 30.00 30.00 29.00 Max 34.00 33.00 33.00 32.00 32.00 32.00

從表6看出,如果以停止時刻τm的中位數(shù)和均值作為判斷標(biāo)準(zhǔn)的話,在γ固定時,m值的變化對停止時刻的影響不大,也就是說,對變點(diǎn)的估計值影響不大.但是,在m值固定時,隨著γ的增大,變點(diǎn)估計值與真實(shí)值的偏差是越來越大的.

從表7可以看出,在N 固定時,隨著γ的增大,變點(diǎn)值的估計偏差也是越來越大的,不過這些偏差的相差量不會很大,不超過1或者2.若γ固定,當(dāng)N 從50變到300時,變點(diǎn)估計值的精確度明顯提高了,在真實(shí)變點(diǎn)=25時,精確度提高了4～5個值.從這里可以看出,增大面板數(shù)據(jù)的橫截面?zhèn)€數(shù)N,是可以明顯提高變點(diǎn)估計的精確度的.

表7:Tm=500,m=100,=25,模擬次數(shù)1000次計算出的停止時刻的各數(shù)字特征

表7:Tm=500,m=100,=25,模擬次數(shù)1000次計算出的停止時刻的各數(shù)字特征

γ N=50 N=300 τm ↓;α → 0.025 0.05 0.10 0.025 0.05 0.10 0 Min 23.00 8.00 8.00 26.00 12.00 10.00 Q1 30.00 30.00 29.00 27.00 27.00 27.00 Median 32.00 31.00 31.00 28.00 28.00 28.00 Mean 31.79 31.13 30.56 27.94 27.66 27.44 Q3 33.00 32.00 32.00 28.00 28.00 28.00 Max 41.00 40.00 39.00 30.00 30.00 30.00 0.25 Min 26.00 23.00 18.00 26.00 26.00 12.00 Q1 31.00 30.00 30.00 27.00 27.00 27.00 Median 32.00 31.00 31.00 28.00 28.00 28.00 Mean 32.08 31.52 30.98 27.99 27.78 27.58 Q3 33.00 33.00 32.00 28.00 28.00 28.00 Max 42.00 41.00 41.00 30.00 30.00 30.00 0.45 Min 27.00 27.00 26.00 26.00 26.00 26.00 Q1 32.00 31.00 30.00 28.00 28.00 27.00 Median 33.00 32.00 32.00 28.00 28.00 28.00 Mean 33.29 32.59 32.01 28.44 28.18 27.98 Q3 35.00 34.00 33.00 29.00 29.00 28.00 Max 45.00 44.00 42.00 31.00 31.00 30.00 0.49 Min 27.00 26.00 26.00 27.00 27.00 26.00 Q1 32.00 32.00 31.00 28.00 28.00 28.00 Median 34.00 33.00 33.00 29.00 28.00 28.00 Mean 34.14 33.30 32.65 28.75 28.45 28.21 Q3 36.00 35.00 34.00 29.00 29.00 29.00 Max 46.00 45.00 43.00 31.00 31.00 31.00

最后看下在N,Tm,m都固定,而變點(diǎn)在不同位置時,變點(diǎn)位置估計值的情況.在表8中,若=25,這種情形在前面已經(jīng)探討過.我們來看下=390,即變點(diǎn)處于右端點(diǎn)的情形.從表格中的數(shù)據(jù)可以發(fā)現(xiàn),對于這種變點(diǎn)最好以變點(diǎn)的中位數(shù)而不是均值作為變點(diǎn)的估計值,因?yàn)橹形粩?shù)會離真實(shí)變點(diǎn)接近的多,估計的偏差在5到6之間.同時γ的大小與估計量的精確度的關(guān)系類似表7中的表現(xiàn),在N固定時,γ越大,估計的偏差也會相對大一些.

通過以上表格的具體分析,發(fā)現(xiàn)參數(shù)γ該如何選取是個很有趣的事情.從表2看出,較小或者中等大的γ(γ=0,0.25)且m取較小或中等大(m=30,100)時,數(shù)據(jù)表明經(jīng)驗(yàn)檢驗(yàn)水平與給定的水平接近.類似地,在表3中,γ=0,N=100同樣可使得經(jīng)驗(yàn)檢驗(yàn)水平接近給定的真實(shí)值α.我們再來看看參數(shù)γ與檢驗(yàn)功效的關(guān)系.表4中的檢驗(yàn)功效對應(yīng)著變點(diǎn)較早發(fā)生的情形,這時,無論γ在區(qū)間取何值,對應(yīng)的檢驗(yàn)功效均為1.從表5看,在變點(diǎn)靠近右端點(diǎn)時,γ越小(其它參數(shù)都固定)時,檢驗(yàn)功效越大.同樣,對表6–8的模擬結(jié)果進(jìn)行分析可知,越小的γ可使變點(diǎn)的估計值與真實(shí)值的偏離程度越小.根據(jù)以上這些分析,若從經(jīng)驗(yàn)檢驗(yàn)水平,檢驗(yàn)功效以及變點(diǎn)估計的精確度三方面來綜合考慮,我們傾向于選取區(qū)間中較小的γ,從而達(dá)到最佳的檢驗(yàn)和估計效果.

表8:N=200,Tm=500,m=100,模擬次數(shù)1000次計算出的停止時刻的各數(shù)字特征

4 定理的證明

在這一節(jié)中,將給出所有引理和定理的證明.

引理2.1的證明 注意到

那么

而且

經(jīng)過計算可得

于是有

和

以及

由(4.2)–(4.6)式以及Chebyshev不等式,得到

定理2.1的證明 借用文獻(xiàn)[3]中的思想來證明定理2.1.對于面板數(shù)據(jù)模型(2.1),在原假設(shè)(i.e.(2.3)式)下,注意到

和

此處{W1,m(t)}和{W2,m(t)}是相互獨(dú)立的Wiener過程,且v>2.從而推得

接下來,利用

由文獻(xiàn)[3]中的定理2.1的證明,有

以及

于是定理2.1得到了證明.

由定理2.1以及Chebyshev不等式,得到

5 結(jié)論

本文針對面板數(shù)據(jù)的序貫?zāi)Ｐ拖碌目赡茏凕c(diǎn),提出了CUSUM型檢驗(yàn)統(tǒng)計量并制定了檢驗(yàn)規(guī)則.隨后得到了檢驗(yàn)統(tǒng)計量的大樣本性質(zhì),并基于相關(guān)理論結(jié)果構(gòu)造了一種檢驗(yàn)方法:漸近檢驗(yàn)法.接下來,對漸近檢驗(yàn)法進(jìn)行了Monte-Carlo數(shù)值模擬,在模擬中對檢驗(yàn)方法的經(jīng)驗(yàn)檢驗(yàn)水平,檢驗(yàn)功效進(jìn)行了考察,并在變點(diǎn)存在的情況下估計出變點(diǎn)的位置.模擬顯示漸近檢驗(yàn)法是一種優(yōu)良的檢驗(yàn)估計方法.最后給出了理論結(jié)果的證明.