非對稱三參數(shù)廣義誤差分布的參數(shù)估計(jì)及應(yīng)用

張文清， 錢夕元

（ 華東理工大學(xué)數(shù)學(xué)學(xué)院，上海 200237）

隨著科學(xué)技術(shù)的高速發(fā)展，各種各樣的數(shù)據(jù)分析方法和統(tǒng)計(jì)方式被廣泛地應(yīng)用于金融經(jīng)濟(jì)和社會活動(dòng)的各個(gè)方面。越來越多的研究表明這些數(shù)據(jù)并不滿足對稱分布，它們往往呈現(xiàn)出帶偏、尖峰厚尾的特性，其中厚尾通常指的是數(shù)據(jù)集不滿足正態(tài)假設(shè)，并且較正態(tài)分布的尾部更厚，比如標(biāo)準(zhǔn)t分布和廣義誤差分布，金融和生物學(xué)領(lǐng)域的數(shù)據(jù)更加呈現(xiàn)這一特點(diǎn)。在這種情況下，如果假定其服從某一對稱分布往往會帶來錯(cuò)誤的分析結(jié)果，從而誤導(dǎo)決策，產(chǎn)生重大損失。因而非對稱及尖峰厚尾分布的研究引起了廣泛關(guān)注，并構(gòu)建了一系列模型和分布族來擬合高度不對稱并且峰度大于正態(tài)分布的數(shù)據(jù)。

廣義誤差分布(Generalized Error Distribution,GED)是指數(shù)族函數(shù)中的一個(gè)對稱單峰函數(shù)，是一種靈活的概率分布函數(shù)，其峰度隨參數(shù)趨向于無窮大。比起標(biāo)準(zhǔn)正態(tài)分布，廣義誤差分布有更厚的尾部，因而被廣泛應(yīng)用于描述金融市場價(jià)格波動(dòng)情況。Subbotin[1]提出，廣義誤差分布包括作為特殊情況的正態(tài)分布、拉普拉斯分布和均勻分布等。Box等[2]在貝葉斯估計(jì)中采用廣義誤差分布來模擬先驗(yàn)密度。Nelson[3]用廣義誤差分布模擬股票市場收益率的分布。Hsieh[4]用廣義誤差分布模擬匯率的分布。同時(shí)廣義誤差分布也是Mcdonald等[5]提出的廣義t分布以及Theodossiou[6]提出的帶偏廣義t分布的特殊情況。與廣義t分布不同，廣義誤差分布的各階矩及矩母函數(shù)都存在。

研究人員一般通過對原有分布進(jìn)行改造來構(gòu)造新的帶偏分布。Cappuccio等[7]在廣義誤差分布中引入偏尾參數(shù)，構(gòu)造了帶偏廣義誤差分布模型，給出了其概率密度函數(shù)的表達(dá)式，并將其應(yīng)用于隨機(jī)波動(dòng)SV模型中。本文首先在廣義誤差分布的基礎(chǔ)上，采用兩個(gè)尾部參數(shù)分別控制左右尾部，并引入偏度參數(shù)，構(gòu)造了非對稱三參數(shù)廣義誤差分布，同時(shí)研究了該分布的基本性質(zhì)，包括累積分布函數(shù)、分位數(shù)函數(shù)及各階原點(diǎn)矩等，并給出了其隨機(jī)變量的抽樣方法；其次分別給出了用矩估計(jì)、極大似然方法和貝葉斯估計(jì)法來估計(jì)該模型參數(shù)的步驟，并通過馬爾科夫蒙特卡羅方法生成的模擬數(shù)據(jù)驗(yàn)證比較了這3種方法；最后將該分布應(yīng)用于兩組實(shí)際數(shù)據(jù)中，驗(yàn)證了非對稱三參數(shù)廣義誤差分布在擬合非對稱、尖峰厚尾數(shù)據(jù)方面優(yōu)于帶偏廣義誤差分布、廣義誤差分布和正態(tài)分布。

1 非對稱三參數(shù)廣義誤差分布

對于隨機(jī)變量X，如果它的概率密度函數(shù)(Probability Distribution Function, PDF)是

當(dāng) α=0.5 且v1=v2時(shí)，AGED的概率密度函數(shù)圖像左右對稱；當(dāng) α∈(0,0.5) 時(shí)，AGED的概率密度函數(shù)圖像向右偏；當(dāng) α∈(0.5,1) 時(shí)，AGED的概率密度函數(shù)圖像向左偏。當(dāng)v1>v2時(shí)，AGED在y軸右側(cè)的圖像尾部比左側(cè)厚；當(dāng)v1＜v2時(shí)，AGED在y軸左側(cè)的圖像尾部比右側(cè)厚。

AGED在不同參數(shù)取值下對應(yīng)的概率密度函數(shù)圖像如圖1所示。

性質(zhì)1如果隨機(jī)變量X～AGED(x;α,v1,v2) ，那么X的累積分布函數(shù)(Cumulative Distribution Function,CDF)為

性質(zhì)3若X為服從 AGED(x;α,v1,v2) 的隨機(jī)變量，Y1為服從 GED(y;v1) 的一個(gè)隨機(jī)變量，Y2為服從GED(y;v2) 的一個(gè)隨機(jī)變量，那么X可以通過式(4)用Y1，Y2進(jìn)行抽樣。

2 參數(shù)估計(jì)

2.1 矩估計(jì)法

矩估計(jì)法首先將總體矩(即所考慮的隨機(jī)變量的冪的期望值)表示為相關(guān)參數(shù)的函數(shù)，然后將這些表達(dá)式等于樣本矩，從而建立方程組，其中方程個(gè)數(shù)與待估計(jì)參數(shù)的個(gè)數(shù)相同。解出方程組即可得到這些參數(shù)的估計(jì)值。

設(shè)隨機(jī)變量Y～GED(y;v) ，則其k階原點(diǎn)矩為

那么 |Y| 的k階原點(diǎn)矩為

設(shè)隨機(jī)變量X～AGED(x;α,v1,v2) ，Y1～GED(y;v1)，Y2～GED(y;v2)，那么期望分別可以用下式表示

因此 |X|k的k階原點(diǎn)矩為

類似的，可以推導(dǎo)出X的k階原點(diǎn)矩

因此X的期望為

X的方差為

假設(shè)X為來自AGED的一個(gè)樣本，X1為X≤0的部分，X2為X>0 的部分，X為X的期望，那么從原點(diǎn)矩的推導(dǎo)過程可知

聯(lián)立方程(6)和(8)，(7)和(9)分別可得

矩估計(jì)法簡單、計(jì)算速度快、計(jì)算難度低，但是往往不能考慮到樣本中的所有相關(guān)信息，甚至?xí)贸鰠?shù)空間以外的估計(jì)值。所以一般把矩估計(jì)的結(jié)果作為極大似然估計(jì)法或者貝葉斯方法M-H鏈的初值。

2.2 極大似然估計(jì)法

極大似然估計(jì)(Maximum Likelihood Estimation,MLE)是一種通過最大化似然函數(shù)來估計(jì)概率分布函數(shù)參數(shù)的方法，其目標(biāo)是找到使似然函數(shù)在參數(shù)空間內(nèi)最大的模型參數(shù)值。

假設(shè)X1,X2,···,Xn獨(dú)立同分布且為來自AGED的一個(gè)樣本，x1,x2,···xn是X1,X2,···,Xn的觀測值。將x1,x2,···xn代入AGED中得到關(guān)于參數(shù) α 、v1和v2的似然函數(shù)

上式兩邊取對數(shù)得

將式(13)分別對 α ，v1和v2求一階偏導(dǎo)得到對數(shù)似然方程組

該對數(shù)似然方程組含有非線性方程，因而無法給出其顯式解，所以這里通過Newton-Raphson迭代法[9]計(jì)算。運(yùn)用Newton-Raphson迭代法需要計(jì)算其二階偏導(dǎo)，該對數(shù)似然函數(shù)的二階偏導(dǎo)分別為

記 θ=(θ1,θ2,θ3)T=(α,v1,v2)T，根據(jù)Newton-Raphson迭代法，迭代關(guān)系式為

極大似然方法邏輯簡單又方法靈活，因此已經(jīng)成為統(tǒng)計(jì)推斷的主要手段，其缺點(diǎn)是計(jì)算量較大。

2.3 貝葉斯估計(jì)法

貝葉斯方法在18世紀(jì)由Thomas Bayes提出。與頻率學(xué)派的觀點(diǎn)不同，貝葉斯方法中參數(shù)為在參數(shù)空間 Θ 內(nèi)取值的一個(gè)隨機(jī)變量 θ 。研究人員用先驗(yàn)分布 π(θ) 來概括觀測數(shù)據(jù)前 θ 的可能值[11]。當(dāng)觀測到數(shù)據(jù)X后，通過后驗(yàn)分布 π(θ|X) 將關(guān)于 θ 的樣本信息與先驗(yàn)信息 π(θ) 結(jié)合。

先驗(yàn)信息使得用貝葉斯估計(jì)進(jìn)行統(tǒng)計(jì)推斷的準(zhǔn)確性更高，因而需要合理選取先驗(yàn)分布。當(dāng)缺少參數(shù)值的分布規(guī)律相關(guān)信息時(shí)，先驗(yàn)分布通常采用Jeffrey先驗(yàn)，對于連續(xù)分布函數(shù)來說，Jeffrey先驗(yàn)是一種標(biāo)準(zhǔn)無信息先驗(yàn)，它在數(shù)值上正比于Fisher信息矩陣I的行列式的平方根。信息矩陣I可以通過對數(shù)似然函數(shù)的海森矩陣進(jìn)行計(jì)算，這里I(θ)=?H(θ) ，θ=(θ1,θ2,θ3)T=(α,v1,v2)T，信息矩陣I具體形式如下

比如

由于馬爾科夫鏈蒙特卡羅方法(Markov Chain Monte Carlo, MCMC)的發(fā)現(xiàn)，貝葉斯推論的研究和應(yīng)用在19世紀(jì)80年代有了巨大的增長。該方法解決了許多計(jì)算問題，使得研究人員對非標(biāo)準(zhǔn)的、復(fù)雜的應(yīng)用也越來越感興趣[12]。MCMC方法中M-H算法和Gibbs算法的使用最為廣泛，本文采用了M-H算法，流程如下：

步驟1：對 α ，v1和v2的建議分布均采用截?cái)嗾龖B(tài)分布。初始化馬爾科夫鏈 θ0=(α0,v10,v20)T，(α0,v10,v20) 可采用矩估計(jì)的結(jié)果。設(shè) θt?1是第t?1 次的迭代值。

步驟5：若算法收斂到參數(shù)的后驗(yàn)分布則停止迭代，否則繼續(xù)迭代該算法。

3 模擬實(shí)驗(yàn)

為了驗(yàn)證以上3種算法是否有效，通過接受拒絕算法，令 α=0.3 、v1=3 、v2=6 ，樣本容量分別取50、500、2000, 得到了3組服從AGED分布的模擬數(shù)據(jù)，然后分別采用這3種算法進(jìn)行參數(shù)估計(jì)，參數(shù)估計(jì)的結(jié)果如圖2所示，其中，ML、MLE、BAYESIAN分別為矩估計(jì)法、極大似然估計(jì)法、貝葉斯估計(jì)法，α、v1、v2均為參數(shù)。從圖中可以看出，隨著樣本容量的增加，參數(shù)估計(jì)的準(zhǔn)確性越來越高，估計(jì)值的離散程度也越來越低，集中分布在真值附近。

圖2 3種方法參數(shù)的估計(jì)結(jié)果Fig. 2 Parameter estimation results of three methods

4 實(shí)例分析

本文通過分析兩組實(shí)際數(shù)據(jù)驗(yàn)證了AGED模型對尖峰厚尾數(shù)據(jù)的擬合效果。使用統(tǒng)計(jì)軟件R（版本4.0.2）進(jìn)行實(shí)例分析，然后比較了模型AGED、SGED（Skew Generalized Error Distribution ）、GED和正態(tài)分布Normal的擬合結(jié)果，這里用貝葉斯方法估計(jì)模型AGED、SGED的參數(shù)，用矩方法估計(jì)GED和正態(tài)分布的參數(shù)。其中，SGED的概率密度函數(shù)是

其中， γ(b;w)=tw?1e?tdt。

4.1 火山高度數(shù)據(jù)

第一個(gè)實(shí)例分析采用全新世(大約過去10000年)期間爆發(fā)的1416座火山高度數(shù)據(jù)集，該數(shù)據(jù)集可以在網(wǎng)站dx.doi.org/10.5479/si.GVP.VOTW4-2013獲得。由于原數(shù)據(jù)最小值與最大值間跨度較大，最小值為?5700，最大值為6879，數(shù)據(jù)集中分布在1694.17附近，所以在用R軟件分析前對數(shù)據(jù)進(jìn)行了處理。表1所示為處理前后該數(shù)據(jù)集的描述統(tǒng)計(jì)量，偏度和峰度在處理前后一致，偏度為0.49，峰度為1.57，這表明非對稱和尖峰厚尾模型適用于分析該數(shù)據(jù)集。

表1 火山高度數(shù)據(jù)的描述統(tǒng)計(jì)量Table 1 Descriptive statistics for the volcano height data

將該數(shù)據(jù)集的經(jīng)驗(yàn)累積分布函數(shù)和擬合的AGED模型的累積分布函數(shù)進(jìn)行比較，結(jié)果如圖3所示。從圖中可以看出，兩條曲線高度重合，表明AGED模型對該數(shù)據(jù)集的擬合效果較好。

圖3 經(jīng)驗(yàn)累積分布函數(shù)和擬合的AGED模型累積分布函數(shù)比較Fig. 3 CDF comparision of empirical model and fitted AGED model

比較AGED分布、SGED分布、GED分布和正態(tài)分布對該數(shù)據(jù)的擬合效果，結(jié)果如圖4所示。從圖中可以看出，AGED模型的擬合效果最好，不僅擬合出了該數(shù)據(jù)集概率分布的形狀，而且很好地反映了數(shù)據(jù)集的概率分布趨勢；正態(tài)分布的擬合效果次之，但是沒能反映數(shù)據(jù)的偏態(tài)和尾部情況；SGED分布和GED分布擬合效果較差，無法擬合出其尖峰的特性，尾部的擬合效果也比較差。

圖4 火山高度數(shù)據(jù)擬合曲線Fig. 4 Fitting curve for the volcano height data

4.2 恒星豐度數(shù)據(jù)

第二個(gè)數(shù)據(jù)集是68個(gè)太陽這類恒星的測量數(shù)據(jù)，數(shù)據(jù)來源于R軟件astrodatR包，數(shù)據(jù)集名為censor_Be，一共68條數(shù)據(jù)，本文分析因變量lg N(Be)，lg N(Be)表示鈹豐度的對數(shù)值。此前Mattos等[13]運(yùn)用偏正態(tài)截尾回歸的尺度混合模型(Scale Mixture of Skew Normal Censored Regression，SMSNCR)分析了該數(shù)據(jù)集，Heleno等[14]在此基礎(chǔ)上運(yùn)用The Asymmetric Alpha-Power Skew-t Distribution進(jìn)行分析。

由于原數(shù)據(jù)過小且分布高度集中，在用R軟件分析前對數(shù)據(jù)進(jìn)行了處理。表2所示為處理前后相關(guān)描述性統(tǒng)計(jì)量，可以發(fā)現(xiàn)處理后數(shù)據(jù)均值減小，標(biāo)準(zhǔn)差變大，偏度與峰度不變，偏度為?1.51，峰度為2.3，表明AGED模型適用于分析該模型。

表2 恒星豐度數(shù)據(jù)的描述統(tǒng)計(jì)量Table 2 Descriptive statistics for the stellar abundances data

將該數(shù)據(jù)集的經(jīng)驗(yàn)累積分布函數(shù)和擬合的AGED模型的累積分布函數(shù)進(jìn)行比較，結(jié)果如圖5所示。從圖中可以看出，在 (?2,?1) 這一區(qū)間上，兩條曲線有些偏差，其他區(qū)間都非常接近甚至重合，表明AGED模型對該數(shù)據(jù)集的擬合效果較好。

圖5 經(jīng)驗(yàn)累積分布函數(shù)和擬合的AGED模型累積分布函數(shù)圖Fig. 5 CDF comparision of empirical model and fitted AGED model

比較AGED分布、SGED分布、GED分布和正態(tài)分布對該數(shù)據(jù)的擬合效果，結(jié)果如圖6所示。

圖6 恒星豐度數(shù)據(jù)擬合曲線Fig. 6 Fitting curve for the stellar abundances data

與4.1節(jié)火山高度數(shù)據(jù)集的擬合效果類似，從圖中可以看出，AGED模型的擬合效果最好，同時(shí)擬合出了該數(shù)據(jù)集概率分布的形狀及其概率分布趨勢；正態(tài)分布的擬合效果僅次于AGED模型，但是沒能反映數(shù)據(jù)的偏態(tài)，該數(shù)據(jù)集左側(cè)厚尾右側(cè)薄尾的情況也未能反映出來；SGED分布對數(shù)據(jù)集小于0的部分?jǐn)M合效果尚可，但是對尖峰和右側(cè)尾部的擬合效果不足；GED分布的擬合效果最差，數(shù)據(jù)集尖峰厚尾非對稱的特性均未能體現(xiàn)。

5 結(jié) 語

本文針對實(shí)際數(shù)據(jù)的尖峰厚尾和非對稱的特性提出了非對稱三參數(shù)廣義誤差分布，該分布在廣義誤差分布的基礎(chǔ)上，通過左尾參數(shù)和右尾參數(shù)分別控制左右兩側(cè)的尖峰厚尾情況，并引入偏度參數(shù)控制偏度。新分布在擬合對稱性和尾部方面有更大的靈活性，便于擬合帶偏厚尾數(shù)據(jù)。文中研究了新分布的理論性質(zhì)和參數(shù)估計(jì)方法，給出了其矩估計(jì)、極大似然估計(jì)和貝葉斯估計(jì)的具體步驟，并通過模擬數(shù)據(jù)檢驗(yàn)了這3種方法的有效性。在火山高度數(shù)據(jù)和恒星豐度數(shù)據(jù)上的應(yīng)用表明，該分布能更好地描述數(shù)據(jù)尖峰厚尾和非對稱的特性，并且貝葉斯估計(jì)對該分布的參數(shù)估計(jì)效果較好。