許一涵,尤添革,溫芫姚,寧 靜,肖揚嵐,尤學(xué)敏

(1. 福建農(nóng)林大學(xué) 計算機與信息學(xué)院,福州 350002;2. 福建省統(tǒng)計信息研究中心,福州 350002;3. 國網(wǎng)信通億力科技有限責(zé)任公司,福州 350002)

我國是一個極端災(zāi)害多發(fā)的國家,巨大自然災(zāi)害不但會影響民眾的生命安全,還會對我國政府造成極大的經(jīng)濟負擔(dān),影響國家經(jīng)濟發(fā)展.2003年非典,2008年四川8.0級大地震,2010年西南大旱災(zāi),2021年鄭州暴雨等等巨災(zāi)事件,都對國家經(jīng)濟造成巨大的損失.

巨災(zāi)的特點是極少發(fā)生,但只要發(fā)生都將會對各方面或者某一方面造成巨大損失,因此,有必要對巨災(zāi)進行相關(guān)分析.巨災(zāi)事件通常都有“尖峰厚尾”的特性,此時極端數(shù)據(jù)比中間數(shù)據(jù)更具有分析價值,研究這部分數(shù)據(jù)對于研究巨災(zāi)風(fēng)險更具實際意義,極值模型正是為此提供了一個合適的模型.

極值模型在水文、降水、地震、工程、環(huán)境、金融及保險等方面都有重要應(yīng)用.近年來,極值理論在巨災(zāi)風(fēng)險建模方面得到了廣泛的應(yīng)用.在地震災(zāi)害方面,主要運用極值理論分析其損失數(shù)據(jù)、震級數(shù)據(jù),計算地震最大震級的分布[1],復(fù)發(fā)規(guī)律并預(yù)測,研究其震級上限及重現(xiàn)水平[2-5],分析地震災(zāi)害損失及其尾部特征[6-7],并計算了VaR、ES等值度量地震災(zāi)害風(fēng)險[8-9],以此得到地震保險純保費和人均保費[10].基于VaR值,建立了風(fēng)險分散層級[11],還利用了CVaR、PML這兩種風(fēng)險測度指標改進VaR的測度[12],近年來,許多學(xué)者還采用混合模型分析我國地震損失[13].以上對地震的研究大多基于其極端值,也有學(xué)者認為中間數(shù)據(jù)對于GPD模型的建立也十分重要,建立了四種混合模型,并提出用貝葉斯法估計參數(shù)和閾值[14].本文將基于極值理論,綜合運用MCMC、極大似然估計等研究方法,分析我國地震損失分布和巨災(zāi)保險政策,為研究我國地震損失提供理論基礎(chǔ).

1 模型與方法

1.1 模型

極值理論是統(tǒng)計建模中的一個重要理論,常常用于研究和分析樣本數(shù)據(jù)中罕見情況.極值理論是由Fisher和Tippet[15](1928)提出的,隨后Gnedenko[16](1943)對此深入分析,之后Gumbel[17](1958)將其標準化.

極值模型主要由兩種方法組成:一種是分塊極值方法(Block Maxima Method,BMM);另一種是超越閾值方法(Peak Over threshold,POT).POT模型由Smith[18](1986)提出,是一種更為有效、應(yīng)用更加廣泛的模型,主要分析數(shù)據(jù)尾部的漸近分布,也就是廣義Pareto分布(Generalized Pareto Distribution,GPD).相比BMM模型,POT模型更充分利用數(shù)據(jù),計算也更準確,是BMM模型的改進.

對某固定的大值μ小于x*,成為閾值,若Xi大于μ,則稱它為超閾值,此時Xi-μ表示超出量,得到

Fμ(x)=Pr(X-μ≤x|X>μ)=

(1)

稱Fμ(x)為超出量分布[19].對應(yīng)的密度函數(shù)為

(2)

(3)

稱為超閾值分布函數(shù),對應(yīng)的密度函數(shù)為

(4)

e(μ)=E(X-μ|x>μ)

(5)

稱為X的平均超出量函數(shù).

考慮它的極限分布,如果X的分布函數(shù)為

(6)

則稱X服從廣義Pareto分布.其中:μ∈R代表位置參數(shù),σ>0代表尺度參數(shù),ξ∈R代表形狀參數(shù).

1.2 研究方法

1.2.1 厚尾性檢驗

通過繪制QQ圖檢驗數(shù)據(jù)的尾部特征,先將數(shù)據(jù)升序排序,然后利用數(shù)據(jù)和標準正態(tài)分布的各個分位點繪制散點圖.分位點計算通常使用:

(7)

如果符合正態(tài)分布,圖像中的大部分點應(yīng)該圍繞著一條直線波動;否則,圖像將在一端或者兩端有擺動.當(dāng)經(jīng)驗分位數(shù)增速較理論分位數(shù)快時,圖像將向上擺動,為厚尾分布;相反,圖像將向下擺動,為短尾分布.

1.2.2 閾值的選取

POT模型主要考察超閾值的樣本數(shù)據(jù),閾值的取值大小會影響研究結(jié)果.如果閾值太大,可用數(shù)據(jù)量太少,會導(dǎo)致方差增大;如果閾值太小,樣本數(shù)據(jù)太多,將產(chǎn)生有偏的參數(shù)估計,不服從GPD分布.因此,合理選擇閾值,既要保證閾值充分大,又要保證數(shù)據(jù)量足夠.本文使用兩種圖示法確定閾值:1)Hill圖法,繪制次序統(tǒng)計量與Hill估計值的圖像,Hill估計式為:

(8)

Hill圖中尾部穩(wěn)定區(qū)域開始時所對應(yīng)的橫坐標Xk就為所要選取的閾值;2)均值超額函數(shù)圖法,取超出量為Yi=Xi-μ,均值超額函數(shù)是閾值μ的線性函數(shù),函數(shù)為:

(9)

在圖中X≥μ時,如果大部分點圍繞在一條直線附近波動,就可以選取這個值為閾值.

1.2.3 參數(shù)估計

1)極大似然估計

(10)

2)Gibbs抽樣

設(shè)θ=(θ1,…,θp)′是p維參數(shù)向量,π(θ|D)是觀察到數(shù)據(jù)集D后θ的后驗分布,則基本Gibbs抽樣方法如下:

第0步,任意選取一個初始點θ(0)=(θ1,0,θ2.0,…,θp,0)′,并置i=0;

第1步,按下列方法生成θ(i+1)=(θ1,i+1,θ2,i+1,…,θp,i+1)′:

生成θ1,i+1～π(θ1|θ2,i…,θp,i,D),

生成θ2,i+1～π(θ2|θ1,i+1,Q3,i,…,θp,i,D),

……

生成θp,i+1～π(θp|θ1,i+1,θ2,i+1,…,θp-1,i+1,D);

第2步,置i=i+1,并返回到第1步.

在這個算法過程中,θ的每一個分量按照自然順序生成,每一個循環(huán)需要生成p個隨機變量.

2 實證分析

2.1 數(shù)據(jù)來源與預(yù)處理

由于地處板塊的相互作用、自然環(huán)境的破壞和人類活動的增加,我國經(jīng)常受到地震的侵擾.地震發(fā)生后,民眾傷亡慘重、經(jīng)濟損失巨大,政府部門和保險機構(gòu)都需要投入巨大的資金,例如2008年汶川8.0級大地震,造成約七萬人死亡,損失金額八千多億.隨著我國現(xiàn)代化進程的不斷加快,交通、建筑等城市設(shè)施不斷地優(yōu)化,人口越來越密集,地震的發(fā)生只會使經(jīng)濟損失和人員傷亡越來越大.

以我國省份空間布局視角,截止2019年,我國32個省市自治區(qū)中有27個發(fā)生過地震,其中四川、云南、新疆、青海省、西藏這五個省市自治區(qū)發(fā)生頻率較高,是我國發(fā)生地震的主要省份.

考慮到1990年以前數(shù)據(jù)誤差較大,決定以中國大陸1990～2019年地震數(shù)據(jù)作為樣本數(shù)據(jù).數(shù)據(jù)中震級均為4級以上;記相同地區(qū)短時間發(fā)生的全部地震和余震的總損失為該次地震經(jīng)濟損失,記最大震級為該次地震震級.本文中所有數(shù)據(jù)均是國家標準統(tǒng)計數(shù)據(jù),來源為各年《中國地震年鑒》等,數(shù)據(jù)處理軟件為R、Openbugs.

由于損失數(shù)據(jù)都是根據(jù)當(dāng)年經(jīng)濟水平記錄的,為了盡量減少經(jīng)濟發(fā)展對研究結(jié)果的影響,使數(shù)據(jù)更具說服力,本文按照2019年的GDP對數(shù)據(jù)進行調(diào)整,即每次的地震損失數(shù)據(jù)乘以2019年的GDP,再除以每年的GDP.

2.2 描述性分析

通過表1描述性統(tǒng)計可以看到,數(shù)據(jù)的最小值與最大值相差巨大;75%分位數(shù)遠小于均值,標準差較大;偏度系數(shù)背離正態(tài)分布0特征值,右偏嚴重;峰度系數(shù)明顯背離正態(tài)分布3特征值.由QQ圖1可見,圖像尾部呈現(xiàn)下凸的形狀,下尾偏離明顯,初步認為數(shù)據(jù)有“尖峰、厚尾”特征.利用散點圖2能夠更直觀看到圖中有一極強影響點發(fā)生于2013年四川省蘆山縣的7級特大地震,本次特大地震直接經(jīng)濟損失1 102.35億元,死亡人數(shù)196人,受傷人數(shù)13 019人.此次地震的特點是震級較大;地震后余震發(fā)生十余次且震級較高,對震區(qū)造成持續(xù)破壞,交通、通信等設(shè)施都有較大程度的毀壞;相比同震級地震,雖然此次地震造成經(jīng)濟損失較大,但是人員傷亡情況較輕,主要由于2008年汶川特大地震災(zāi)后重建、加固工作較好,大部分建筑的主體結(jié)構(gòu)并沒有倒塌,人群密集的醫(yī)院、學(xué)校等大型建筑毀壞、倒塌程度較輕.

圖1 地震災(zāi)害損失數(shù)據(jù)QQ圖Figure 1 QQ chart of earthquake disaster damage data

圖2 地震災(zāi)害損失數(shù)據(jù)散點圖Figure 2 Scatterplot of earthquake disaster damage data

表1 中國大陸地震災(zāi)害損失數(shù)據(jù)描述性統(tǒng)計(億元)

2.3 模型估計與分析

以上結(jié)果可以得出正態(tài)分布并不能很好的擬合地震損失數(shù)據(jù),考慮使用GPD分布來擬合.首先對數(shù)據(jù)進行Box-Cox變換得到變換參數(shù)λ近似為0,即對數(shù)據(jù)進行對數(shù)變換.通過圖3中的Hill圖和均值超額函數(shù)圖確定閾值,比較不同閾值下GPD分布參數(shù)(σ,ξ)的極大似然估計及相應(yīng)的95%置信區(qū)間,發(fā)現(xiàn)閾值μ=3.4較合理,此時該閾值附近呈現(xiàn)平穩(wěn)狀態(tài),建立POT模型得到參數(shù)結(jié)果如表2所示,此時超閾值樣本數(shù)為48個,數(shù)據(jù)量滿足模型的需要,也就是說這48個數(shù)據(jù)屬于高損數(shù)據(jù).并且通過GPD分布擬合圖和殘差QQ圖檢驗GPD分布擬合情況,由圖4可以看出擬合效果較好.使用極大似然估計法對廣義帕累托分布的參數(shù)進行估計,得到結(jié)果表2、診斷圖5,在P-P圖和分位數(shù)圖中只有幾個強影響點偏離直線;重現(xiàn)水平圖中幾乎所有的點都落在了置信區(qū)間內(nèi);密度曲線的估計和直方圖基本一致,說明極大似然法可以估計該模型參數(shù).

圖3 Hill圖和均值超額函數(shù)圖Figure 3 Hill chart and mean-excess function chart

圖4 GPD分布擬合圖和殘差QQ圖Figure 4 GPD distribution fit plots and residual QQ plots

表2 POT模型參數(shù)估計表

圖6 核密度圖、迭代軌跡圖和自相關(guān)系數(shù)圖Figure 6 Nuclear density diagram, iterative trajectory diagram and auto-correlation graph

比較極大似然法和MCMC法估計不同閾值時的參數(shù),結(jié)果如表3所示,當(dāng)超閾值樣本較多時,兩者估計結(jié)果相近;當(dāng)超閾值樣本較少時,極大似然法變得很不穩(wěn)定,此時兩者的參數(shù)ξ估計值相差0.807 7,因此,小樣本情況下選擇MCMC法估計參數(shù)更加理想.

表3 極大似然法與 MCMC 法參數(shù)估計結(jié)果比較Table 3 Comparison of parameter estimation results between the maximum likelihood method and MCMC method

表4 我國地震災(zāi)害損失金額的高分位數(shù)點估計值(億元)Table 4 High-quartile point estimates of the amount of earthquake disaster damage in China (billion yuan)

3 結(jié) 語

本文基于GPD分布與Gibbs抽樣法對我國地震損失數(shù)據(jù)展開實證研究,研究結(jié)果表明:1)我國地震損失數(shù)據(jù)具有明顯的“尖峰、厚尾”特性,經(jīng)分布擬合圖、殘差QQ圖等驗證了GPD分布能夠很好擬合此特性,并通過Hill圖、均值超額函數(shù)圖等確定最優(yōu)閾值μ=3.4;2)基于該閾值建立POT模型,使用極大似然法估計參數(shù),并通過P-P圖、分位數(shù)圖、重現(xiàn)水平圖等檢驗;3)在參數(shù)估計中,考慮到極大似然法在小樣本情形下可能失效的問題,再次采用蒙特卡洛法進行參數(shù)估計,并通過調(diào)整閾值、控制樣本量來比較不同數(shù)據(jù)量下兩者的參數(shù)估計結(jié)果,發(fā)現(xiàn)在小樣本情形下蒙特卡洛法比極大似然法估計效果更好;4)通過高分位數(shù)點估計值,根據(jù)我國實際國情,簡單說明適合我國的巨災(zāi)保險制度.本文運用POT模型充分利用了樣本數(shù)據(jù),保留了巨災(zāi)風(fēng)險的厚尾特征,能更好地把握巨災(zāi)數(shù)據(jù),對于研究巨災(zāi)損失具有重要意義.