洪水頻率分布的尾指數(shù)估計

歐陽資生，甘 柳

0 引言

在各種水利規(guī)劃設計中，都需要進行水文頻率分析工作，以保證達到工程要求的設計值。而在洪水頻率分析中，無非是討論像洪水或超出某一警戒水平的水位或流量的分布情況。在極值統(tǒng)計中，像洪水等這類事件發(fā)生頻率很低，但是一經(jīng)發(fā)生就容易造成較大損失的事件叫極值事件。事實上，在洪災風險管理中，如何發(fā)現(xiàn)這些極值事件的發(fā)生概率和某個極值分位數(shù)點對風險管理者是相當重要的。我們知道，從極值理論的角度來看，這些極值事件的概率和極值分位數(shù)被一種稱為極值分布的尾指數(shù)所控制。對于正態(tài)分布，它的尾部呈指數(shù)函數(shù)衰減，其尾指數(shù)為零。當尾指數(shù)大于零時，分布尾部呈冪函數(shù)衰減，我們稱其為厚尾分布。厚尾模型在諸如金融、保險、水利等很多場合都是一個應用非常廣泛的分布模型。厚尾分布中，尾指數(shù)越大，其尾部越厚。因此尾指數(shù)可作為衡量洪水分布規(guī)律的重要指標。

假設X1,X2,???,Xn是一列正的，獨立同分布的隨機變量序列，具有共同的分布函數(shù)F(x)：

則分布函數(shù)F(x)被稱為厚尾的。模型（1）的F(x)也稱為Pareto型分布。?F是無窮遠處的緩變函數(shù)，滿足對所有的η＞0：

這里，γ稱為極值指數(shù)，而如果令α=1/γ，則α稱為尾指數(shù)。尾指數(shù)或極值指數(shù)的估計目前仍然是極值統(tǒng)計的一個熱點問題。如Danielsson(2001)利用自助法研究了尾指數(shù)的估計問題。Beirlant et al.(2008)提出了一種修正的偏差減少方法對尾指數(shù)和極值分位數(shù)進行估計。Beran and Schell(2010)在小樣本情況下構造了一個穩(wěn)健的M估計方法對尾指數(shù)進行了估計。Brito and Freitas(2010)研究了相依數(shù)據(jù)尾指數(shù)估計的相合性問題。但總的說來，對尾指數(shù)或極值指數(shù)的的估計，不外乎是Pickands估計、Hill估計和矩估計三種估計及其拓展估計方法。在這三種估計中，Hill估計和矩估計實際應用中相對較多，而Pickands估計實用性并不強，基本不用。

如果設：

為 n個觀測值 X1,X2,???,Xn的順序統(tǒng)計量，Hill估計定義為：

Hill估計雖然在理論上具有很好的大樣本性質，但是在實際應用中，并不好操作，我們可以從圖1就可很容易理解。

圖1 極值指數(shù)的Hill估計圖

圖1 是我們利用學生-t4分布對Hill估計結果作的一個隨機模擬，我們作了200次隨機模擬，每次模擬的樣本量是500，我們給出了模擬的200次估計的1/4分位數(shù)，中位數(shù)和3/4分位數(shù)，在學生-t4分布中，極值指數(shù)的真值γ=0.25，但是極值指數(shù)的真值到底取多少，我們從圖1中很難做出判斷。換句話說，我們不知道門限值取多大時才能對樣本進行有效分割。事實上，如何選取合適的門限值是估計極值指數(shù)或尾指數(shù)的基礎，也是我們進行洪水頻率分析的必要程序。

本文中，我們將基于指數(shù)回歸模型，給出矩估計的門限值和樣本點分割的選取原理和方法，然后利用MC方法進行模擬說明門限值選取的合理性，最后利用所構建的模型對湖南省四個水文觀測站的水文數(shù)據(jù)進行實證分析。

1 尾指數(shù)估計的統(tǒng)計建模

通過前面的模擬，我們看到Hill估計并不好操作，門限值選擇不當將導致極值指數(shù)估計的較大偏差。事實上，極值指數(shù)三種估計方法中，矩估計相對來說較為穩(wěn)健，實際應用中也相對較多。因此，本文我們采用矩估計來估計極值指數(shù)和尾指數(shù)。

矩估計由Dekkers、Einmahl&de Hall(1989)提出，其定義為：對 k∈{3,???,(n-1)}

盡管矩估計較Hill估計穩(wěn)健，在理論上也具有更好的大樣本性質，但是在實際應用中，也同樣存在門限值選取的問題。我們可以從圖2就很容易發(fā)現(xiàn)，雖然相對Hill估計而言，矩估計要穩(wěn)健得多，但是門限值的選取問題仍然是一個需要解決的問題。

圖2 極值指數(shù)的矩估計圖

圖2 是我們利用學生-t1分布對矩估計結果作的一個隨機模擬，我們同樣作了200次隨機模擬，每次模擬的樣本量大小為500，我們給出了模擬的500次估計的中位數(shù)。在學生-t1分布中，極值指數(shù)的真值為γ=1。雖然比起Hill估計而言要穩(wěn)健得多，但是，極值指數(shù)的真值到底取多少，我們從圖2中也不好判斷。換句話說，我們不知道門限值取多少時，才能對樣本實現(xiàn)最優(yōu)分割。

和其他估計一樣，在用矩估計對尾指數(shù)進行估計時，首先是確定門限值，找出超出門限值以上的觀察數(shù)據(jù)；也就是對所觀察到的樣本值的順序統(tǒng)計量進行有效分割，得到用于估計的觀察數(shù)據(jù)，然后才能進行估計。但需要指出的是，門限值的選取問題卻一直是困擾極值工作者的一個難題。門限值越大，可以分析的數(shù)據(jù)越少，這時，被分析的數(shù)據(jù)比較接近分布的極端，分析的偏差減少，但由于數(shù)據(jù)過少，估計的方差增加；反之，門限值過小，被分析的數(shù)據(jù)增加，分析的方差減少，但偏差卻增加了。對這個問題的研究，統(tǒng)計工作者提出了許多方案。如Dupuis(1998)建議從參數(shù)的穩(wěn)健性出發(fā)來確定門限值；Guillou(2001)、Matthys&Beirlant(2003),Beirlant et al.(1996,2004)，歐陽資生(2008)等建議使用最小化均方誤差或漸近二階矩來獲得門限值；Gomes et al.(2008)建議使二階參數(shù)估計的偏差達到最小從而通過一個啟發(fā)式適應過程得到門限選擇方法；Vandewalle et al.(2008)通過使PDC估計(partial density component estimation)的積分均方誤差達到最小來獲得門限值，對樣本進行分割。

采用SPSS 19.0軟件對數(shù)據(jù)進行分析處理，計量資料以（均數(shù)±標準差）表示，采用t檢驗；計數(shù)資料以（n，%）表示，采用χ2檢驗，以P＜0.05表示差異具有統(tǒng)計學意義。

下面：我們將基于指數(shù)回歸模型，在漸近最小均方誤差的準則下，給出矩估計的門限值和樣本點分割的選取原理和方法，并提出極值指數(shù)和尾指數(shù)估計的算法。

對于隨機變量序列X1,X2,???,Xn的分布函數(shù)F(x)，如前所述，我們假設F(x)是Pareto型的，其原因主要是基于Pareto型在極值分布中的地位和作用。我們知道，這種分布在金融、保險、水利中都被廣泛應用。例如，在巨災統(tǒng)計數(shù)據(jù)中，有一個廣為人知的事實，即巨災統(tǒng)計數(shù)據(jù)是厚尾的，因此，可以直接假設巨災統(tǒng)計數(shù)據(jù)分布服從Pareto分布。

對于模型(1)中的緩變函數(shù)，有一個被廣泛接受的假設：

假設?λ:存在一個實常數(shù)ρ＜0和一個正的比率函數(shù)b(x)，滿足當x→∞時，b(x)→0，且使得對所有的η≥1,

其中 kρ(η)= ∫1ηvρ-1d v=(ηρ-1)/ρ ，若 ρ=0 ，則 kρ(η)=log(η)。需要說明的是，假設?λ條件并不苛刻，一般的緩變函數(shù)均能滿足這個條件。

在假設?λ下，我們按照Beirlant et al.(2004)，歐陽資生(2008)的指數(shù)回歸模型方法來選取k，從而進一步確定參數(shù) γ?k,b?n,k,ρ?k。為此，建立如下指數(shù)回歸模型：

這里，f1,f2,???,fk是一列獨立的，服從標準指數(shù)分布的隨機變量。在式(6)中，利用最大似然估計，得到參數(shù)γ,bn,k,ρ的估計值：

類似于Beirlant et al.(2004)，歐陽資生(2008)，我們可得在指數(shù)回歸模型中，極值指數(shù)用矩估計作為估計量時的AMSE為：

因此，樣本的最優(yōu)分割k?optn為：

因此，根據(jù)以上原理，我們可得基于指數(shù)回歸模型的樣本分割方法，進而得到尾指數(shù)的矩估計的算法，算法如下：

(1)對指數(shù)回歸模型式(6)，利用極大似然估計，對k∈{3,???,(n-1)}計算參數(shù) γ ,bn,k,ρ 的估計值{(γ?k,b?n,k,ρ?k),k∈{3,???,(n-1)}

(2)對 k∈{3,???,(n-1)}計算 AMSE(γ?Mk)

(3)利用

獲得 k?optn

(4)根據(jù)矩估計式(4)和步驟3的最優(yōu)k?optn，可得極值指數(shù)的最優(yōu)估計和尾指數(shù)估計α?。

2 隨機模擬

現(xiàn)在，為驗證我們的模型，我們對下列極值分布進行蒙特卡洛模擬：

(1)Burr(1,1,1)分布。 Burr(θ,τ,λ)分布的分布函數(shù)滿足：

(2)Burr(1,0.5,2)分布

(3)Frechet(1)分布.Fre ch et(γ)分布的分布函數(shù)滿足：

(4)Frechet(2)分布

(5)學生-t4分布

表1 極值指數(shù)估計及其誤差估計模擬結果表

在蒙特卡洛模擬時，我們對每一種分布作了500次模擬，每次模擬的樣本量均為1000。表1分別給出了500次模擬中相應的最優(yōu)k值、γ的估計值的平均及其標準差、AMSE的平均。從表1可以看出，在矩估計中，借助于指數(shù)回歸模型獲得門限值、樣本點分割方法和極值指數(shù)估計值，其結果是令人非常滿意的。

3 湖南省洪水頻率分析的實證分析

3.1 數(shù)據(jù)描述

作為模型的一個應用，我們對洞庭湖周邊的桃源、津市、沙頭、石龜山等四個站點的水文數(shù)據(jù)中的水流量進行實證分析。數(shù)據(jù)跨度為1998年元月1日至2010年4月1日共4316個日數(shù)據(jù)。為對數(shù)據(jù)的基本情況有一直觀了解，我們在表2中列出了相應的統(tǒng)計量。同時，也繪畫了其相應的時間序列圖(圖3)。從表2可以發(fā)現(xiàn)，這四個站點的數(shù)據(jù)均呈現(xiàn)明顯的厚尾現(xiàn)象，同時，從圖3也可看出其波動明顯。

表2 洞庭湖周邊的四個水文觀測站水流量基本統(tǒng)計特征

圖3 （1998.1.1-2010.4.1）洞庭湖周邊四個水文觀測站水流量時間序列圖

3.2 尾指數(shù)的估計

根據(jù)前文的極值指數(shù)估計模型和計算方法，我們首先可以得到樣本的最優(yōu)分割方法從而得到超出門限值的樣本個數(shù)，然后得到極值指數(shù)的估計值，最后利用α?=1/γ?即可得湖南省四個水文觀測站點的水流量分布的尾指數(shù)估計。

在圖4中，我們給出了k=1,2,...,4310時極值指數(shù)的估計圖,其中虛線代表了最優(yōu)的k值。而表3分別給出了這四個站點在最優(yōu)的k值下的極值指數(shù)和尾指數(shù)最優(yōu)估計值，從表3可以看出，桃源站的尾指數(shù)是最大的，津市的尾指數(shù)是最小的。

兩點說明：(1)從前文對學生-t1分布的極值指數(shù)的隨機模擬圖2可看出，即使是估計時相對穩(wěn)健的矩估計也不好直接判斷尾指數(shù)的真值，這點在圖4中也得到了印證。從圖4中，如果沒有一定的準則，我們是無法獲知在何時對樣本實施分割，也就無法得到各個站點水流量數(shù)據(jù)的尾指數(shù)的真值了。但是，如果我們借助于指數(shù)回歸模型，采用使矩的AMSE達到最小作為評價標準，我們就可以很好地解決一個問題。(2)正如Beirlant et al.(1996,2004)在借助指數(shù)回歸模型對尾指數(shù)進行Hill估計時作的評述，指數(shù)回歸模型相對來說，較為穩(wěn)健。同時，由于有效的利用了極大似然估計，因此計算的速度也較快，這也是我們在進行矩估計時，借助指數(shù)回歸模型進行建模的主要原因。

圖4 洞庭湖周邊的桃源等四個水文站點水流量數(shù)據(jù)的極值指數(shù)估計圖

4 結論

以上四個水文觀測站點都是位于洞庭湖地區(qū)周邊站點，由前面的分析結果可以發(fā)現(xiàn)，四水系流域的水位變化情況均為厚尾分布，都可以通過極值分布加以較好地擬合。當然，流經(jīng)不同站點的水流量是不一樣的，且不同規(guī)模洪水流量的變化幅度亦有所區(qū)別，因此在實施防洪措施時應實事求是，依據(jù)不同的情況有區(qū)別的對待，這樣才能既做到全面有效防洪減災又能盡可能的降低不必要投入，減少浪費。

[1]Beirlant,J.,Figueiredo,F.,Gomes,M.I.,Vandewalle,B.Improved Re?duced-bias Tail Index and Quantile Estimators[J].J.Statist.Plann.and Inference，2008,138.

[2]Beirlant,J.,Goegebeur,Y.,Segers,J.,Teugels,J.Statistics of Ex?tremes.Theory and Applications[M].NewYork:Wiley,2004.

[3]Beirlant,J.,Vynckier,P.,Teugels,J.L.Tail Index Estimation,Pareto Quantile Plots,and Regression Diagnostics[J].J.Amer.Statist.Assoc，1996,91.

[4]Beran,J.,Schell,D.On Robust Tail Index Estimation.Computational Statisticsand Data Analysis,doi:10.1016/j.csda[J].2010.

[5]Brito,M.,Freitas,A.C.Consistent Erstimation of the Tail Index for De?pendent Data[J].Statistics and Probability Letters，2010,(80).

[6]Danielsson,J.Using a Bootstrap Method to Choose the Sample Fraction in Tail Index Estimation[J].Journal of Multivariate Analysis，2001,76.

[7]Dekkers,A.,de Haans,L.A Moment Estimator for the Index of an Ex?treme-value Distribution[J].Ann Statist,1989,17(4).

[8]Dupuis,D.J.Exceedances over High Thresholds:a Guide to Thresh?old Selection[J].Extremes,1998,3(1).

[9]Gomes,M.I.,Henriques Rodrigues,L,Vandewalle,B.，Viseu,C.A Heu?ristic Adaptive Choice of the Threshold for Bias-corrected Hill Esti?mators[J].J.Statist.Comput.And Simulation,2008,78(2).

[10]Guillou,A.,Hall,P.A Diagnostic for Selecting the Threshold in Ex?treme Analysis[J].J.R.Statist.Soc.Ser B,2001,63.

[11]Matthys,G.,Beirlant,J.Estimating the Extreme Value Index and High Quantiles with Exponential Regression Models[J].Statistica Si?nica,2003,13.

[12]Vandewalle,B.,Beirlant,J.,Christmann,A.,Hubert,M.A Robust Estimator for the Tail Index of Pareto-type Distributions[J].Compu?tational Statistics&Data Analysis,2007,51.

[13]歐陽資生.厚尾分布的極值分位數(shù)估計與極值風險測度研究[J].數(shù)理統(tǒng)計與管理,2008,27.