張文清， 錢夕元

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

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

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

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

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

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

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

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

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

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

2 參數(shù)估計

2.1 矩估計法

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

設隨機變量Y～GED(y;v) ，則其k階原點矩為

那么 |Y| 的k階原點矩為

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

因此 |X|k的k階原點矩為

類似的，可以推導出X的k階原點矩

因此X的期望為

X的方差為

假設X為來自AGED的一個樣本，X1為X≤0的部分，X2為X>0 的部分，X為X的期望，那么從原點矩的推導過程可知

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

矩估計法簡單、計算速度快、計算難度低，但是往往不能考慮到樣本中的所有相關信息，甚至會得出參數(shù)空間以外的估計值。所以一般把矩估計的結果作為極大似然估計法或者貝葉斯方法M-H鏈的初值。

2.2 極大似然估計法

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

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

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

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

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

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

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

2.3 貝葉斯估計法

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

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

比如

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

步驟1：對 α ，v1和v2的建議分布均采用截斷正態(tài)分布。初始化馬爾科夫鏈 θ0=(α0,v10,v20)T，(α0,v10,v20) 可采用矩估計的結果。設 θt?1是第t?1 次的迭代值。

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

3 模擬實驗

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

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

4 實例分析

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

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

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

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

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

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

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

比較AGED分布、SGED分布、GED分布和正態(tài)分布對該數(shù)據(jù)的擬合效果，結果如圖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ù)

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

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

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

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

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

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

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

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

5 結 語

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