殷崔紅,林小東,2,袁海麗
(1.廈門大學數學科學學院,福建廈門361005)
(2.多倫多大學統(tǒng)計系,加拿大多倫多M5S 3G3)
(3.武漢大學數學與統(tǒng)計學院,湖北武漢430072)
厄蘭極值混合模型的有效估計及其在保險中的應用
殷崔紅1,林小東1,2,袁海麗3
(1.廈門大學數學科學學院,福建廈門361005)
(2.多倫多大學統(tǒng)計系,加拿大多倫多M5S 3G3)
(3.武漢大學數學與統(tǒng)計學院,湖北武漢430072)
本文研究了Erlang混合分布和廣義帕累托分布混合模型的估計問題.通過引入iSCAD懲罰函數,利用EM算法極大化iSCAD懲罰似然函數的方法,獲得了混合序和參數的估計值,計算出有效的度量風險指標value-at-risk(VaR)和tail-VaR(TVaR),通過模擬實驗和實際數據說明了模型和算法的有效性.推廣了有限Erlang極值混合模型在保險數據擬合中的應用.
極值理論;極值混合模型;iSCAD懲罰;EM算法;似然函數
Erlang混合分布廣泛應用于保險損失數據的建模,在保險破產理論和保險損失數據的擬合中都有良好的表現.保險破產理論中,當利用混合Erlang分布對保險損失的嚴重程度建模時,通常關注的一些指標將有明確的解析式,比如無限破產概率,隨機破產時刻的拉普拉斯變換等,這方面的研究可參考文獻[3,17,21,29];近幾年,學者更多關注于將Erlang混合分布用于擬合保險實際損失數據,得到了很多令人滿意的分布性質,比如分布函數和矩都有解析式,使得相關的風險測度value-at-risk(VaR)和tail VaR(TVaR)比較容易計算.Verbelen[30]等將雙邊截斷引入Erlang混合分布,計算了再保險合同的純保費.類似研究見文獻[9, 10,19,26,30]等.Lee和Lin[20]提出Erlang混合分布的多元形式,多元混合Erlang分布保留了一元Erlang混合分布的大部分有用的分布性質,同時建模相依性,與copula方法相呼應.關于多元Erlang混合分布的研究見文獻[2,16,31,32]等.
混合模型的首要問題是混合序的確定,Lee和Lin[19,30]等都利用BIC來確定Erlang混合分布的序,Yin和Lin[33]提出了一種新的iSCAD懲罰函數,建立懲罰似然函數,運用EM算法給出參數的估計,同時給出了混合序的估計.然而值得注意的是:Erlang分布是輕尾的,用它來擬合重尾數據時可能很難達到預期效果.其次,尾部數據相應的權重一般都很小,公式(3.11)可以看出,權重小于閾值λ的相應Erlang分布都被刪除,這不利用保留擬合尾部數據的Erlang分布.為解決這些問題,本文引入極值理論擬合尾部數據,建立Erlang極值混合模型.
極值混合模型廣泛應用于各領域的數據分析中,尤其在保險、金融、水文和環(huán)境科學等領域.在保險領域,大額索賠在保險公司的風險管理和產品定價,尤其是再保險產品的定價方面,有不可忽略的意義.文獻[4,12,13,23,27]等將極值理論引入到保險的風險管理中.為使數據的主體和尾部都擬合的很好,Behrens[5]提出單一參數分布與一個極值分布的混合模型, Carreau和Bengio[7]討論混合參數分布與極值分布的混合模型,類似的文獻[5,6,15,22, 24]等給出多種極值混合模型.Lee[18]等最早將極值混合模型引入到保險數據中,但是所有這些混合模型都沒有考慮混合序的確定.
本文建立Erlang混合分布與廣義帕累托(GPD)分布的混合模型,廣義帕累托(GPD)分布用于擬合數據的尾部,而Erlang混合分布用于擬合數據的主體,這樣即有Erlang混合分布的優(yōu)點,同時保留了極值理論的長處.引入iSCAD懲罰來估計混合Erlang分布的參數, Yin和Lin[33]已證明參數和混合序的估計都有一致性.
首先給出Erlang分布的密度函數為
其中γ是取值為正整數的形狀參數(shape parameter),θ>0是尺度參數(scale parameter).
將m個不同的Erlang分布以權重α=(α1,···,αm)混合,則Erlang混合分布的密度函數為
相應的分布函數為F(x;α,γ,θ),其中權重參數α=(α1,···,αm)滿足αj≥0和γ=(γ1,···,γm)是形狀參數,為了可識別性的說明,一般有γ1≤···≤γm,而θ>0是共有的尺度參數.
由于投保人的性別、車型、駕車經驗和熟悉程度等的不同,使得索賠數據一般有明顯的異質性,單一的Erlang分布可能很難給出好的擬合效果,因此數據的主體部分本文仍然選用Erlang混合分布來擬合,而尾部采用極值分布.故本文采用左右雙邊截斷的Erlang混合分布,大部分保險損失數據都是已知截斷值,比如保險中的免賠額和賠償限額.以l和μ分別表示左右截斷值(免賠額l已知),雙邊截斷的Erlang混合分布的密度函數是
其中
顯然,(2.2)式是左右截斷點為l和μ的Erlang分布f(x;l,μ,γj,θ)的混合模型,混合權重為π=(π1,···,πm),滿足πj≥0和=1.密度函數(2.1)相應的分布函數為F(x;l,μ,α,γ,θ).
在統(tǒng)計中,廣義帕累托分布(Generalized Pareto Distribution,i.e.,GPD)經常被用于擬合其他分布或實際數據的尾部,本文選用GPD擬合數據的尾部,其密度函數是
廣義帕累托分布(GPD)的生存函數為
結合(2.1)和(2.4)式,為彌補引言中提過的Erlang混合模型的不足,本文建立的Erlang極值混合模型的密度函數為
其中μ為閾值,X是服從h(x;l,μ,α,γ,θ,ξ,σ)分布的隨機變量,令ψμ=P(X>μ),其一般由大于μ的樣本比例來估計.
相應的生存函數為
風險測度就是各種風險度量指標的總稱.現行的國際標準風險管理工具VaR最初由Morgan針對銀行業(yè)務風險的需要提出的,并很快被推廣成為了一種產業(yè)標準.風險價值VaR是指在正常的市場條件、給定的置信水平以及給定的持有期間內,投資組合所面臨的潛在最大損失.VaR是一種分位數風險測度,一般給定置信水平p,典型的p=95%或者99%.但是,VaR作為風險測度只考慮了概率為p的事件的最大損失VaRp,高于VaRp的損失并沒有納入風險測度,為克服這個缺陷,Tail Value at Risk(or TVaR)被提出來.在給定置信水平p下,TVaR就是損失落入最糟的1-p部分的平均損失.下面給出Erlang極值混合模型關于風險指標VaR和TVaR的計算.
為便于風險指標VaRp和TVaRp的計算,當X≤μ時,將生存函數(2.7)用Erlang密度函數分別表達為
假設損失隨機變量X服從Erlang極值混合分布(2.6),給定置信水平p,有
方程(2.9)的解即置信水平為p的VaRp.
計算TVaRp之前,首先研究自付責任額為R(>l)的再保險的純保費,當R≤μ時,
當R>μ時,
綜上,自付責任額為R(>l)的再保險的純保費為
當自付責任額R=VaRp時,置信水平為p的TVaRp為
文獻[33]針對每一個分量權重參數πj,j=1,···,m,提出的iSCAD懲罰函數為
其中I(·)是示性函數.本文建立的Erlang極值混合分布中Eralng混合分布的參數估計與新引入的極值分布的參數估計互不影響,因此關于Eralng混合分布的極大懲罰似然估計仍然是一致的.
Expectation-Maximization(EM)算法最早由Dempster[11]給出比較詳細的說明,當似然函數的最大值點不能直接得到時,EM算法通過迭代的方法找到最大值點.EM算法需引入隱變量,隱變量可以是未知參數,丟失的數據或者任何可以使模型簡化的未觀測數據量. EM算法分為E-step和M-step兩步,其中E-step計算目標函數關于隱變量Z的條件期望, M-step是最大化目標函數,求得參數的極大似然估計.王繼霞等[1]將EM算法用于有限混合Laplace分布的估計.
Erlang極值混合模型的所有待估參數是:擬合數據主體部分的Erlang混合分布的序m,形狀參數γ=(γ1,···,γm),相應的權重參數α=(α1,···,αm),所有Erlang分布共用的尺度參數θ,擬合數據尾部的廣義帕累托分布的閾值μ,尺度參數σ,形狀參數ξ,下面逐一介紹它們的估計.
由公式(2.2)知,密度函數(2.6)也可以由新權重參數π表示為
假設X=(X1,···,Xn)是獨立同分布的隨機變量,服從密度函數h(x;l,μ,π,γ,θ,ξ,σ),即(3.2),樣本觀測值為x=(x1,···,xn),相應有序樣本觀測值為x(1)≤···≤x(n),記
Pickands[25]給出與閾值μ相應的k的選擇方法,從1開始依次增加,最大值為[n/4],而μ=x(n-k),本文最終由似然函數的大小選出μ.為方便后面的說明,重新表示n'=n-k和x'=(x(1),···,x(n')).
形狀參數的估計采用Yin和Lin[33]類似的方法,即預先給定一個大的混合序M,形狀參數的所有可能取值是γ0=(),通過估計相應的權重參數,來實現混合序的估計和形狀參數的選擇.
Erlang極值混合分布的密度函數h(x;l,μ,π,γ0,θ,ξ,σ)中的部分未知參數記為φ= (π1,···,πM,θ),本文采用EM算法來估計φ.
樣本x=(x1,···,xn)的對數似然函數為
樣本x=(x1,···,xn)的iSCAD懲罰對數似然函數,其中與參數φ=(π1,···,πM,θ)有關的部分是
直接關于?n',P(φ;x)求極大似然估計是困難的,本文使用EM算法,引入隱變量,即Z=(Z1,···,Zn),其中Zi=(Zij|i=1,···,n,j=1,···,M),
那么完整樣本(x,Z)的似然函數為
相應完整樣本(x,Z)的對數似然函數為
相應的完整樣本(x,Z)的iSCAD懲罰對數似然函數為
EM算法是利用迭代過程來估計參數的方法,假設已經完成第k次迭代,獲得的當前估計是φ(k)=EM算法的E-step和M-step分別為
E-step ?n,P(φ;x,Z)關于隱變量Z求條件期望,得到關于可觀測樣本x的邊際似然函數,即
M-step(3.9)式是權重參數πj(j=1,···,M)和尺度參數θ的函數,求函數(3.9)的極大估計,即
權重參數πj的第(k+1)次迭代的估計為
尺度參數θ的第(k+1)次迭代的估計為
其中
迭代過程一直持續(xù)到|Q(φ(k))-Q(φ(k-1))|小于某個既定的誤差界.分別以0,j=1,···,M}和表示EM迭代的最終結果.混合模型序的估計是
最后,關于廣義帕累托分布(GPD)的尺度參數σ和形狀參數ξ的極大似然估計,Coles[8]已經詳細討論過,本文就不再作重復說明.
本文利用R軟件進行計算,基于Yin和Lin[33]關于Erlang混合分布的R程序和軟件包“ismev”,編寫本文Erlang混合分布和GPD分布混合模型的R程序,完成模擬實驗和實際數據中模型參數的估計.
為驗證模型和估計的有效性,本文給出一個模擬實驗,從密度函數(2.6)中隨機抽取了2500個隨機數,其中(2.6)式中的所有參數見表1中的真實參數.
表1:真實參數與參數估計值的對比
參數的初始化主要參考文獻[19,28].事先給定M=10,形狀參數的備擇范圍即γ=(1,···,10),以Tijms[28]的方法初始化,公式(3.11)給出極大懲罰似然的權重參數估計,其稀疏性實現了在形狀參數備擇范圍γ=(1,···,10)中進行合理選擇.從表1可以看出,形狀參數最終僅選中=(2,7),只有這兩個形狀參數對應的權重參數估計為非零的,即=(0.501,0.499),其它形狀參數相應的權重參數估計均為零,即=0,j=1,3,4,5,6,8,9,10.顯然,混合模型序的估計
由本實驗可以看出,引入iSCAD懲罰的優(yōu)勢所在:通過對權重參數的估計,同時實現了對形狀參數的選擇和混合模型序的估計.表1列出的所有參數估計值與真實值都很接近,說明模型和算法都很有效,能夠反映出數據的特征.圖1很好的反應了這一點,圖1中的真實曲線和擬合曲線幾乎是重合的.
圖1:模擬數據的直方圖,真實曲線與擬合曲線
丹麥火災賠償數據有2167個觀測值,Embrechts[14]和Mendes[24]等都用極值理論研究過這組數據的尾部,本文采用Erlang極值混合模型從總體上研究這組數據,不再僅僅限于研究其尾部特征.
文獻[33]討論了帶左截斷點l的Erlang混合分布,本文在其基礎上提出了Erlang極值混合分布,在本例中將利用這兩種不同的分布分別擬合丹麥火災賠償數據,比較兩種分布的優(yōu)劣.
表2給出Erlang混合分布和Erlang極值混合分布(2.6)擬合火災損失數據得到的所有參數的估計值,其中利用Erlang極值混合分布得到的結果說明擬合數據的主體部分采用了三個Erlang分布,數據的尾部由廣義帕累托分布來擬合,兩部分的閾值點為4.174,尾部數據比例為0.152;而利用Erlang混合分布擬合同一組火災數據則需要十個不同的Erlang分布的混合.
表2 :參數估計值
圖2:丹麥火災數據的直方圖與擬合曲線
圖2是丹麥火災數據的直方圖、Erlang混合分布和Erlang極值混合分布的擬合曲線,可以看出擬合效果較好.
圖3和4分別給出Erlang混合分布和Erlang極值混合分布的Q-Q圖,顯然Erlang極值混合分布在尾部數據的擬合上更優(yōu).
本文給出VaR的非參數(nonparametric)法估計作為標桿,在置信水平為p的條件下, VaRp的非參數估計是方程Fn(VaRp)=p的解,其中Fn(x)=
表3:非參數法、Erlang混合分布和Erlang極值混合分布的VaRp值的比較
圖3:丹麥火災數據的Q-Q圖(Erlang混合分布)
圖4:丹麥火災數據的Q-Q圖(Erlang極值混合分布)
表3給出三種方法的VaRp估計值,表3可以看出,Erlang極限混合分布估計得到的VaRp與非參數法得到的VaRp非常接近,估計效果很好.
表4:非參數法、Erlang混合分布和Erlang極值混合分布的TVaRp值的比較
表4給出非參數法、Erlang混合分布和Erlang極值混合分布的TVaRp估計值,其中TVaR的非參數估計為TVaRp=.Erlang混合分布的TVaRp比非參數法的結果偏小,這主要是因為Erlang混合分布對火災損失數據的尾部擬合不足,見圖3;而Erlang極值混合分布的結果稍大,而且越到尾部,這種趨勢越明顯,這主要是因為估計得到的=0.661>0,即估計的極值分布為厚尾的,而實際數據的尾部過于稀疏,不足以表現這種厚尾性.
[1]王繼霞,汪春峰,苗雨.有限混合Laplace分布回歸模型局部估計的EM算法(英文)[J/OL].數學雜志, 2016,36(4):667-675.
[2]Badescu A L,Lan G,Lin X S,et al.Modeling correlated frequencies with application in operational risk management[J].J.Oper.Risk,2015,10(1):1-43.
[4]Beirlant J,Goegebeur Y,Segers J,et al.Statistics of extremes:theory and applications[M].Ltd England:John Wiley&Sons,2006.
[5]Behrens C N,Lopes H F,Gamerman D.Bayesian analysis of extreme events with threshold estimation[J].Stat.Model.,2004,4(3):227-244.
[6]Carreau J,Bengio Y.A hybrid Pareto model for asymmetric fat-tail data[R].Technical Report 1283, Canada:Dept.IRO,Universitde Montral,2006.
[7]Carreau J,Bengio Y.A hybrid pareto mixture for conditional asymmetric fat-tailed distributions[J]. Neural Net.,IEEE Trans.,2009,20(7):1087-1101.
[8]Coles S,Bawa J,Trenner L,et al.An introduction to statistical modeling of extreme values[M]. London:Springer,2001.
[9]Cossette H,Mailhot M,Marceau.TVaR-based capital allocation for multivariate compound distributions with positive continuous claim amounts[J].Insur.:Math.Econ.,2012,50(2):247-256.
[10]Cossette H,CotM P,Marceau E,et al.Multivariate distribution defined with Farlie Gumbel Morgenstern copula and mixed Erlang marginals:Aggregation and capital allocation[J].Insur.: Math.Econ.,2013,52(3):560-572.
[11]Dempster A P,Laird N M,Rubin D B.Maximum likelihood from incomplete data via the EM algorithm[J].J.Royal Stat.Soc.,Ser.B(Methodological),1977:1-38.
[12]Embrechts P,Klppelberg C,Mikosch T.Modelling extremal events for insurance and finance[J]. Springer,1997,71(2):183-199.
[13]Embrechts P,Resnick S I,Samorodnitsky G.Extreme value theory as a risk management tool[J]. North Amer.Actua.J.,1999,3(2):30-41.
[14]Embrechts P,Klppelberg C,Mikosch T.Modelling extremal events:for insurance and finance[M]. Germany:Springer Sci.Business Media,2013.
[15]Frigessi A,Haug O,Rue H.A dynamic mixture model for unsupervised tail estimation without threshold selection[J].Extremes,2002,5(3):219-235.
[16]Hashorva E,Ratovomirija G.On Sarmanov mixed Erlang risks in insurance applications[J].Astin Bull.,2015,45(01):175-205.
[17]Landriault D,Willmot G E.On the joint distributions of the time to ruin,the surplus prior to ruin, and the deficit at ruin in the classical risk model[J].North Amer.Actua.J.,2009,13(2):252-270.
[18]Lee D,Li W K,Wong T S T.Modeling insurance claims via a mixture exponential model combined with peaks-over-threshold approach[J].Insur.:Math.Econ.,2012,51(3):538-550.
[19]Lee S C K,Lin X S.Modeling and evaluating insurance losses via mixtures of Erlang distributions[J]. North Amer.Actua.J.,2010,14(1):107-130.
[20]Lee S C K,Lin X S.Modeling dependent risks with multivariate Erlang mixtures[J].Astin Bull., 2012,42(01):153-180.
[21]Lin X S,Willmot G E.The moments of the time of ruin,the surplus before ruin,and the deficit at ruin[J].Insur.:Math.Econ.,2000,27(1):19-44.
[22]MacDonald A,Scarrott C J,Lee D,et al.A flexible extreme value mixture model[J].Comput.Stat. Data Anal.,2011,55(6):2137-2157.
[23]McNeil A J.Estimating the tails of loss severity distributions using extreme value theory[J].Astin Bull.,1997,27(01):117-137.
[24]Melo Mendes B V,Lopes H F.Data driven estimates for mixtures[J].Comput.Stat.Data Anal., 2004,47(3):583-598.
[25]Pickands III J.Statistical inference using extreme order statistics[J].Ann.Stat.,1975:119-131.
[26]Porth L,Zhu W,Seng Tan K.A credibility-based Erlang mixture model for pricing crop reinsurance[J].Agricul.Finance Rev.,2014,74(2):162-187.
[27]Resnick S I.Discussion of the Danish data on large fire insurance losses[J].Astin Bull.,1997,27(01): 139-151.
[28]Tijms H C.A first course in stochastic models[M].UK:John Wiley and Sons,2003.
[29]Tsai C C L,Willmot G E.On the moments of the surplus process perturbed by diffusion[J].Insur.: Math.Econ.,2002,31(3):327-350.
[30]Verbelen R,Gong L,Antonio K,et al.Fitting mixtures of Erlangs to censored and truncated data using the EM algorithm[J].Astin Bull.,2015,45(03):729-758.
[31]Verbelen R,Antonio K,Claeskens G.Multivariate mixtures of Erlangs for density estimation under censoring[J].Life.Data Anal.,2015:1-27.
[32]Willmot G E,Woo J K.On some properties of a class of multivariate Erlang mixtures with insurance applications[J].Astin Bull.,2015,45(01):151-173.
[33]Yin C,Lin X S.Efficient estimation of Erlang mixtures using iSCAD penalty with insurance application[J].Astin Bull.,Available on CJO2016,doi:10.1017/asb.2016.14.
2010 MR Subject Classification:62E15;62F10
EFFICIENT ESTIMATION OF ERLANG AND GPD MIXTURES USING ISCAD PENALTY WITH INSURANCE APPLICATION
YIN Cui-hong1,LIN Xiao-dong1,2,YUAN Hai-li3
(1.School of Mathematical Sciences,Xiamen University,Xiamen 361005,China)
(2.Department of Statistical Sciences,University of Toronto,Ontario M5S 3G3,Canada)
(3.School of Mathematics and Statistics Sciences,Wuhan University,Wuhan 430072,China)
In this paper,we study efficient estimation of Erlang&GPD mixture model.By using a new thresholding penalty function and a corresponding EM algorithm,we estimate model parameters and determine the order of the mixture model.We obtain risk measure including VaR and TVaR and show efficiency of the new mixture model in simulation studies and a real data application,which improve Erlang&extreme value mixture model in modeling insurance losses.
extreme value theory;mixture model;iSCAD penalty;EM algorithm;likelihood function
MR(2010)主題分類號:62E15;62F10O212.1
A
0255-7797(2016)06-1315-13
?2016-04-09接收日期:2016-06-28
國家自然科學基金資助(11201352).
殷崔紅(1982-),女,山東濰坊,博士,主要研究方向:非壽險精算.