基于POT模型的我国地震损失分布研究

许一涵,尤添革,温芫姚,宁 静,肖扬岚,尤学敏

(1. 福建农林大学 计算机与信息学院,福州 350002;2. 福建省统计信息研究中心,福州 350002;3. 国网信通亿力科技有限责任公司,福州 350002)

我国是一个极端灾害多发的国家,巨大自然灾害不但会影响民众的生命安全,还会对我国政府造成极大的经济负担,影响国家经济发展.2003年非典,2008年四川8.0级大地震,2010年西南大旱灾,2021年郑州暴雨等等巨灾事件,都对国家经济造成巨大的损失.

巨灾的特点是极少发生,但只要发生都将会对各方面或者某一方面造成巨大损失,因此,有必要对巨灾进行相关分析.巨灾事件通常都有“尖峰厚尾”的特性,此时极端数据比中间数据更具有分析价值,研究这部分数据对于研究巨灾风险更具实际意义,极值模型正是为此提供了一个合适的模型.

极值模型在水文、降水、地震、工程、环境、金融及保险等方面都有重要应用.近年来,极值理论在巨灾风险建模方面得到了广泛的应用.在地震灾害方面,主要运用极值理论分析其损失数据、震级数据,计算地震最大震级的分布[1],复发规律并预测,研究其震级上限及重现水平[2-5],分析地震灾害损失及其尾部特征[6-7],并计算了VaR、ES等值度量地震灾害风险[8-9],以此得到地震保险纯保费和人均保费[10].基于VaR值,建立了风险分散层级[11],还利用了CVaR、PML这两种风险测度指标改进VaR的测度[12],近年来,许多学者还采用混合模型分析我国地震损失[13].以上对地震的研究大多基于其极端值,也有学者认为中间数据对于GPD模型的建立也十分重要,建立了四种混合模型,并提出用贝叶斯法估计参数和阈值[14].本文将基于极值理论,综合运用MCMC、极大似然估计等研究方法,分析我国地震损失分布和巨灾保险政策,为研究我国地震损失提供理论基础.

1 模型与方法

1.1 模型

极值理论是统计建模中的一个重要理论,常常用于研究和分析样本数据中罕见情况.极值理论是由Fisher和Tippet[15](1928)提出的,随后Gnedenko[16](1943)对此深入分析,之后Gumbel[17](1958)将其标准化.

极值模型主要由两种方法组成:一种是分块极值方法(Block Maxima Method,BMM);另一种是超越阈值方法(Peak Over threshold,POT).POT模型由Smith[18](1986)提出,是一种更为有效、应用更加广泛的模型,主要分析数据尾部的渐近分布,也就是广义Pareto分布(Generalized Pareto Distribution,GPD).相比BMM模型,POT模型更充分利用数据,计算也更准确,是BMM模型的改进.

对某固定的大值μ小于x*,成为阈值,若Xi大于μ,则称它为超阈值,此时Xi-μ表示超出量,得到

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

(1)

称Fμ(x)为超出量分布[19].对应的密度函数为

(2)

(3)

称为超阈值分布函数,对应的密度函数为

(4)

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

(5)

称为X的平均超出量函数.

考虑它的极限分布,如果X的分布函数为

(6)

则称X服从广义Pareto分布.其中:μ∈R代表位置参数,σ>0代表尺度参数,ξ∈R代表形状参数.

1.2 研究方法

1.2.1 厚尾性检验

通过绘制QQ图检验数据的尾部特征,先将数据升序排序,然后利用数据和标准正态分布的各个分位点绘制散点图.分位点计算通常使用:

(7)

如果符合正态分布,图像中的大部分点应该围绕着一条直线波动;否则,图像将在一端或者两端有摆动.当经验分位数增速较理论分位数快时,图像将向上摆动,为厚尾分布;相反,图像将向下摆动,为短尾分布.

1.2.2 阈值的选取

POT模型主要考察超阈值的样本数据,阈值的取值大小会影响研究结果.如果阈值太大,可用数据量太少,会导致方差增大;如果阈值太小,样本数据太多,将产生有偏的参数估计,不服从GPD分布.因此,合理选择阈值,既要保证阈值充分大,又要保证数据量足够.本文使用两种图示法确定阈值:1)Hill图法,绘制次序统计量与Hill估计值的图像,Hill估计式为:

(8)

Hill图中尾部稳定区域开始时所对应的横坐标Xk就为所要选取的阈值;2)均值超额函数图法,取超出量为Yi=Xi-μ,均值超额函数是阈值μ的线性函数,函数为:

(9)

在图中X≥μ时,如果大部分点围绕在一条直线附近波动,就可以选取这个值为阈值.

1.2.3 参数估计

1)极大似然估计

(10)

2)Gibbs抽样

设θ=(θ1,…,θp)′是p维参数向量,π(θ|D)是观察到数据集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步.

在这个算法过程中,θ的每一个分量按照自然顺序生成,每一个循环需要生成p个随机变量.

2 实证分析

2.1 数据来源与预处理

由于地处板块的相互作用、自然环境的破坏和人类活动的增加,我国经常受到地震的侵扰.地震发生后,民众伤亡惨重、经济损失巨大,政府部门和保险机构都需要投入巨大的资金,例如2008年汶川8.0级大地震,造成约七万人死亡,损失金额八千多亿.随着我国现代化进程的不断加快,交通、建筑等城市设施不断地优化,人口越来越密集,地震的发生只会使经济损失和人员伤亡越来越大.

以我国省份空间布局视角,截止2019年,我国32个省市自治区中有27个发生过地震,其中四川、云南、新疆、青海省、西藏这五个省市自治区发生频率较高,是我国发生地震的主要省份.

考虑到1990年以前数据误差较大,决定以中国大陆1990～2019年地震数据作为样本数据.数据中震级均为4级以上;记相同地区短时间发生的全部地震和余震的总损失为该次地震经济损失,记最大震级为该次地震震级.本文中所有数据均是国家标准统计数据,来源为各年《中国地震年鉴》等,数据处理软件为R、Openbugs.

由于损失数据都是根据当年经济水平记录的,为了尽量减少经济发展对研究结果的影响,使数据更具说服力,本文按照2019年的GDP对数据进行调整,即每次的地震损失数据乘以2019年的GDP,再除以每年的GDP.

2.2 描述性分析

通过表1描述性统计可以看到,数据的最小值与最大值相差巨大;75%分位数远小于均值,标准差较大;偏度系数背离正态分布0特征值,右偏严重;峰度系数明显背离正态分布3特征值.由QQ图1可见,图像尾部呈现下凸的形状,下尾偏离明显,初步认为数据有“尖峰、厚尾”特征.利用散点图2能够更直观看到图中有一极强影响点发生于2013年四川省芦山县的7级特大地震,本次特大地震直接经济损失1 102.35亿元,死亡人数196人,受伤人数13 019人.此次地震的特点是震级较大;地震后余震发生十余次且震级较高,对震区造成持续破坏,交通、通信等设施都有较大程度的毁坏;相比同震级地震,虽然此次地震造成经济损失较大,但是人员伤亡情况较轻,主要由于2008年汶川特大地震灾后重建、加固工作较好,大部分建筑的主体结构并没有倒塌,人群密集的医院、学校等大型建筑毁坏、倒塌程度较轻.

图1 地震灾害损失数据QQ图Figure 1 QQ chart of earthquake disaster damage data

图2 地震灾害损失数据散点图Figure 2 Scatterplot of earthquake disaster damage data

表1 中国大陆地震灾害损失数据描述性统计(亿元)

2.3 模型估计与分析

以上结果可以得出正态分布并不能很好的拟合地震损失数据,考虑使用GPD分布来拟合.首先对数据进行Box-Cox变换得到变换参数λ近似为0,即对数据进行对数变换.通过图3中的Hill图和均值超额函数图确定阈值,比较不同阈值下GPD分布参数(σ,ξ)的极大似然估计及相应的95%置信区间,发现阈值μ=3.4较合理,此时该阈值附近呈现平稳状态,建立POT模型得到参数结果如表2所示,此时超阈值样本数为48个,数据量满足模型的需要,也就是说这48个数据属于高损数据.并且通过GPD分布拟合图和残差QQ图检验GPD分布拟合情况,由图4可以看出拟合效果较好.使用极大似然估计法对广义帕累托分布的参数进行估计,得到结果表2、诊断图5,在P-P图和分位数图中只有几个强影响点偏离直线;重现水平图中几乎所有的点都落在了置信区间内;密度曲线的估计和直方图基本一致,说明极大似然法可以估计该模型参数.

图3 Hill图和均值超额函数图Figure 3 Hill chart and mean-excess function chart

图4 GPD分布拟合图和残差QQ图Figure 4 GPD distribution fit plots and residual QQ plots

表2 POT模型参数估计表

图6 核密度图、迭代轨迹图和自相关系数图Figure 6 Nuclear density diagram, iterative trajectory diagram and auto-correlation graph

比较极大似然法和MCMC法估计不同阈值时的参数,结果如表3所示,当超阈值样本较多时,两者估计结果相近;当超阈值样本较少时,极大似然法变得很不稳定,此时两者的参数ξ估计值相差0.807 7,因此,小样本情况下选择MCMC法估计参数更加理想.

表3 极大似然法与 MCMC 法参数估计结果比较Table 3 Comparison of parameter estimation results between the maximum likelihood method and MCMC method

表4 我国地震灾害损失金额的高分位数点估计值(亿元)Table 4 High-quartile point estimates of the amount of earthquake disaster damage in China (billion yuan)

3 结 语

本文基于GPD分布与Gibbs抽样法对我国地震损失数据展开实证研究,研究结果表明:1)我国地震损失数据具有明显的“尖峰、厚尾”特性,经分布拟合图、残差QQ图等验证了GPD分布能够很好拟合此特性,并通过Hill图、均值超额函数图等确定最优阈值μ=3.4;2)基于该阈值建立POT模型,使用极大似然法估计参数,并通过P-P图、分位数图、重现水平图等检验;3)在参数估计中,考虑到极大似然法在小样本情形下可能失效的问题,再次采用蒙特卡洛法进行参数估计,并通过调整阈值、控制样本量来比较不同数据量下两者的参数估计结果,发现在小样本情形下蒙特卡洛法比极大似然法估计效果更好;4)通过高分位数点估计值,根据我国实际国情,简单说明适合我国的巨灾保险制度.本文运用POT模型充分利用了样本数据,保留了巨灾风险的厚尾特征,能更好地把握巨灾数据,对于研究巨灾损失具有重要意义.