基于联合分位数回归的地震灾害二元损失评估模型及其应用

杨雅明 张 晴 李秀芳

一、引言

我国位于环太平洋地震带和欧亚地震带之间，受太平洋板块、印度洋板块和菲律宾海板块的挤压，地震灾害频发。20世纪以来，我国共发生6级以上的地震近800次，因地震造成的死亡人数约占所有自然灾害死亡总人数的54%。其中，2008年汶川地震造成8 352亿元直接经济损失，受灾总人口达4 625.6万人 ；2010年玉树地震造成228.5亿元直接经济损失，2 698人死亡。为降低地震对国家经济和生产生活的影响，政府需要建立起完善有效的巨灾风险管理制度。长期以来，我国巨灾风险损失补偿主要以国家财政救济和社会捐助为主，随着我国社会财富的集中程度和人口密度的上升，各种巨灾所带来的损失程度和影响范围不断扩大，传统基于政府救济式的巨灾管理模式已经无法满足需求。近些年来，政府针对巨灾风险管理问题提出了一系列的建议，2021年《中华人民共和国国民经济和社会发展第十四个五年规划和2035年远景目标纲要》中提出要“优化国家应急管理能力体系建设，提高防灾减灾抗灾救灾能力”并进一步“发展巨灾保险”。

进行巨灾风险管理的前提是对灾害造成的损失进行量化分析，建立精确且贴切实际的巨灾风险评估模型。目前评估地震经济损失的模型主要是基于分布拟合的方法(田玲和姚鹏，2013［1］；郝军章和崔玉杰，2016［2］)。由于地震巨灾经济损失数据具有厚尾特征，单一的损失分布模型难以对其尾部数据进行有效拟合。基于此，郭静和张连增(2021)［3］提出对阈值之前损失数据用Mixed Erlang分布进行拟合，阈值之后尾部数据用帕累托分布进行拟合，即构建组合分布模型。李云仙等(2017)［4］利用对数正态分布和广义帕累托分布的组合分布模型对1966—2014年间5.0～7.0级地震灾害经济损失进行拟合。与单分布相比，组合分布模型能更精确地拟合巨灾的经济损失数据。地震灾害生命损失数据一般为计数型变量，具有零膨胀、过度分散的特征，孟生旺和李政宵(2018)［5］分别用右截断负二项分布和右截断广义帕累托分布拟合死亡人数，用Panjer迭代法和快速傅里叶变换计算地震死亡人数在特定时期的分布以及风险度量值。

鉴于巨灾损失数据具有尖峰厚尾的统计特征，以VaR等测度为代表的尾部风险测度成为巨灾风险量化的重要工具，除以上基于分布拟合方法为基础的量化方法外，分位数回归模型则是另一类尾部风险量化的主要方法。自Koenker和Basset(1978)［6］的开创性工作以来，分位数回归模型不断发展，Koenker和Machado(1999)［7］研究了非对称拉普拉斯与分位数回归之间的联系，随后Yu(2001)［8］提出了利用贝叶斯方法对分位数回归时的参数进行估计，Benoit 和Vanden(2017)［9］系统总结了分位数回归的贝叶斯估计方法，并分别介绍了因变量为连续和离散情况下的模型及变量选择方法。如今，分位数回归方法在地震巨灾风险量化中已有众多应用。例如：田玲等(2019)［10］利用半参数分位数回归模型对1990—2015年地震经济损失数据进行建模。李云仙和孟生旺(2019)［11］分别利用函数系数的分位数回归模型、传统分位数回归模型对我国1990—2015年地震灾害直接经济损失数据进行建模。此外，分位数回归模型也被应用于地震生命损失量化中，如Jiang等(2018)［12］利用1969—2006年我国地震灾害死亡人数建立了针对计数数据的贝叶斯半参数分位数回归模型。

上述讨论的分布拟合与分位数回归的风险量化思路各具利弊。依托分布拟合方法可以方便计算更多风险测度指标并提供较多估计信息，但受限于巨灾损失的低频高损与厚尾特征，分布拟合过程中会存在分布选择主观性与阈值设定偏误等问题。依托分位数回归方法虽无需设定响应变量的分布且能将多个解释变量融入分位数回归，但该方法仅能得到分位点估计，而其他分布信息，例如另一常用尾部风险测度ES的评估结果则无法获取。除此之外，通过文献梳理发现当前上述两类方法在地震巨灾研究中的应用大多从财产或人身单一损失的角度出发，而忽略了二元损失之间的相依性。根据刘新红等(2019)［13］研究结论可知，地震直接经济损失和死亡人数之间存在正相依关系，若对两者进行独立性假设，风险度量值VaR会被低估，故考虑损失之间的相依关系是必要的。一些学者已经展开了基于分位数回归的多元变量之间相依性的研究，如Petrella和Raponi(2019)［14］将基于非对称拉普拉斯分布的分位数回归模型推广到多元的背景下，解释多个因变量之间的相依关系，并运用EM算法对相关系数矩阵进行估计。在此基础上，Merlo等(2021)［15］利用联合分位数回归模型对金融资产的VaR和ES进行估计，他们的研究为本文综合建立经济损失和生命损失的二元分位数回归模型提供了可能性。

综上所述，为刻画经济损失和生命损失之间的相依结构，量化不同震级下地震灾害的二元损失特征，本文将在联合分位数回归的框架下，运用EM算法联合估计震级对死亡人数和直接经济损失的回归系数，以验证震级与地震损失之间的高度相关性，最终通过计算VaR和ES对地震巨灾风险进行评估。在实证研究部分，本文收集了我国大陆地区1990—2020年间306次成灾地震的损失数据，建立联合分位数回归模型，并通过计算不同震级在不同概率水平下的VaR值和ES值进行风险评估，最后利用核估计方法对震级发生情况进行分布估计，从而实现对我国地震巨灾指数保险的总保费规模的测算。本文的主要贡献体现在两个方面：第一，在联合分位数回归模型的框架下评估地震的二元损失，综合评估直接经济损失和死亡人数的风险，并将两者的相依性纳入模型。第二，不同于大多数学者利用随机模拟法计算地震巨灾损失的ES(李云仙等，2019［11］；刘新红等，2019［12］)，本文基于多元非对称拉普拉斯尺度参数与ES之间的联系，可实现对地震巨灾损失的VaR和ES联合估计。

二、基于地震损失的联合分位数回归模型

本文用直接经济损失和死亡人数来评估地震灾害的二元损失，参考Machado和Silva(2005)［16］提出的方法对死亡人数进行连续化处理，并根据Petrella等(2019)［14］的研究，构建二元损失的联合分位数回归模型。

(一)单变量分位数回归模型与非对称拉普拉斯分布的初步研究

为了更好地解释多元非对称拉普拉斯分布与联合分位数回归之间的联系，本节简要阐述单变量分位数回归模型与非对称拉普拉斯分布(Asymmetric Laplace distribution，AL)的直接联系。随机变量Y服从位置参数μ、尺度参数δ>0和偏度参数τ∈(0，1)的非对称拉普拉斯分布，即AL(μ，δ，τ)，其概率密度函数为：

fAL(y;μ,δ,τ)=τ(1-τ)δ-1exp[-ρτ{(y-μ)/δ}]

上式中ρτ(t)=t{τ-1(t<0)}为损失(或检验)函数，其中1为指示函数，ρτ{(y-μ)/δ}服从参数为1/δ的指数分布。上式的AL分布还可以被表示为：

其中U服从参数为1/δ的指数分布，Z服从标准正态分布。另外，为了保证参数μ在给定的水平τ上与Y的分位数重合，ξτ和θτ必须满足以下条件：

ξτ=(1-2τ)/{τ(1-τ)},θτ=2/{τ(1-τ)}

对于i∈{1，…，n}，yi为被解释变量，xi为k×1维解释变量，Qyi(τ|xi)是在给定xi时yi的分位数回归函数，则Qyi(τ|xi)与xi的关系如下：

其中τ∈(0，1)为分位数水平，βτ为k×1维回归参数。分位数回归模型主要用于研究响应变量分位数与解释变量之间的关系，分位数回归模型的参数估计可通过最小化目标函数实现：

(1)

(2)

(二)联合分位数回归与多元非对称拉普拉斯分布

根据单变量分位数回归模型与非对称拉普拉斯分布的直接联系，本节利用多元非对称拉普拉斯(Multivariate Asymmetric Laplace，MAL)分布来联合建模多元响应变量的条件分位数。令Yi=(Yi1，…，Yip)Τ是对于i∈{1，2…n}的p维响应变量，Xi是k×1维自变量，对于j∈{1，…，p}，联合分位数回归模型如下：

其中βτ=(βτ1，…，βτp)Τ是p×k维未知系数矩阵，βτj=(β1，τj，…，βk，τj)，上式等价于如下多元线性回归模型：

Yi=βτXi+εi.

(3)

(三)基于二元损失的半参数分位数回归模型构建

在模型构建前，本节需要对响应变量进行处理，令二元变量Zi=(Zi1，Zi2)Τ分别为地震造成的直接经济损失和死亡人数，其中i∈1，2，…，n。由于直接经济损失数据常存在极端值，因此本节将对直接经济损失进行对数处理。另外，本节采用Machado和Silva(2005)［16］提出的方法，将死亡人数转化为连续变量，且该连续变量的分位数与原变量之间存在着对应关系。令Yi1=log(Zi1)为对数处理后的直接经济损失，Wi2=Zi2+Ui，其中Ui是服从区间[0，1]上的均匀分布随机变量。设τ=(τ1，τ2)为二元损失变量分位数水平，构造单调变换T(Wi2；τ2)如下：

其中，ζ为一个较小正数，沿用Machado和Silva(2005)［16］的设定，ζ取10-5。令Yi2=T(Wi2；τ2)，经过前述变换，二元损失变量的联合分位数回归模型可表示为：

其中，解释变量Xi为地震震级，βτ0=(βτ10,βτ20)′,βτ1=(βτ11, βτ21)′。

(四)基于EM算法的参数估计

本节利用EM算法对模型参数进行极大似然估计。EM算法本质上是在计算期望(E)步和最大化(M)步之间交替进行，期望(E)步计算了当前参数估计的对数似然函数的期望，最大化(M)步通过最大化E步中获得的对数似然函数的期望来进行参数估计。根据Petrella等(2019)［14］给出的两个命题可进行对数似然函数期望的计算和最优参数估计。

EM算法通过E步和M步的循环迭代计算来实现，该迭代过程一直进行直到满足收敛条件，本文将收敛判据设定为待估参数连续两次迭代差的二范数小于10-6。

三、地震风险测度及地震指数保险保费规模测算

(一)风险测度值

VaR和ES是广泛使用的风险度量工具。在本文中，VaR可以理解为在给定概率水平τ下，地震灾害可能造成的最大损失值。VaR的数学表达式为：

VaRτ=inf{l∈:p(Y>l)≤τ}

其中，Y表示地震灾害损失变量，p为概率水平，l为实数，为实数集。根据分位数回归的定义，直接经济损失的VaR可表示为：

死亡人数的VaR为：

由于VaR无法充分反映尾部风险信息且不满足次可加性要求，本文进一步采用了一致性风险测度ES刻画尾部风险。ES不仅继承了VaR的诸多优点，还提供了更多关于损失尾部的信息(Artzner等，1999［17］)。因此，ES作为风险度量的替代指标，越来越受风险管理者的关注，但在普通分位数回归框架下ES与VaR之间无法关联。Bassett等(2004)［18］提出了分位数回归在基于非对称拉普拉斯分布情况下尺度参数与ES之间的数量关系，因此将分位数回归方法的评估作用进行了提升，Taylor(2019)［19］则进一步将这一关联应用于测度金融资产组合收益Rt的左侧尾部风险ES，公式如下：

本文测度目标为损失变量的右尾较大损失分位数下的ES，且并非基于时变尺度参数假设，因此基于上述理论框架，测算过程中按照损失额的负数进行参数估计，且将δt的处理退化为常数δ。

(二)核密度估计

x1，x2，…，xn为收集的样本数据，设其密度函数为f(x)，则f(x)的非参数核密度估计为：

其中，k(·)为核函数，h为窗宽，n为样本容量。在本文中，核函数选取为高斯核函数。然而在不同窗宽下，核密度估计曲线形状差距比较大，一般最佳窗宽选择为h=cn-1/5(其中c为常数)，通过不断调整c，使得核密度估计达到满意的估计结果。常见的窗宽的计算公式为：

h=0.785(x0.75-x0.25)n-1/5

(4)

其中x0.75-x0.25为数据的0.75分位数与0.25分位数之差，本节选用此公式计算的窗宽h与MATLAB根据拇指经验法则计算的最优默认窗宽w进行比较，选择对密度函数估计效果更好的窗宽。

(三)地震巨灾保险保费规模测算思路

参考李云仙和孟生旺(2019)［11］关于地震指数保险保费规模测算方法，本文将震级设定为地震巨灾保险的“指数”，针对地震造成的直接经济损失，以事先约定的地震指数为基础，通过分位数回归模型量化地震指数与直接经济损失之间的关系，从而实现对总保费规模的测算。对于分位数回归模型而言，要先利用核密度法估计地震震级的概率密度函数，然后求出在给定概率水平下的VaR均值，最后将VaR均值与地震平均发生频次相乘计算总保费。VaR均值计算公式如下：

其中，a、b为地震震级上下限，f(x)为地震震级的概率密度函数。则每年的总保费规模为：

其中，P为每年保费规模，N为地震发生次数。

四、基于我国地震灾害数据的实证分析

(一)数据来源及描述性统计

本文的数据来自2005—2020年的《中国大陆地震灾害损失述评》［20-33］及中国地震信息网(http：//www.csi.ac.cn/)，从中获取了1990—2020年共计306次地震灾害的震级、直接经济损失(单位：亿元)和死亡人数数据。由于时间跨度较长，所以本节将直接经济损失数据调整到2020年的物价水平，以消除通货膨胀的影响。表1描述了直接经济损失和死亡人数的统计特征。通过观察发现，两个变量均具有尖峰厚尾的特点，且直接经济损失离散程度较大，本节将对其进行对数处理。针对具有零膨胀特征的死亡人数，需要根据第二节的方法进行连续化处理。通过图1的直方图以及QQ图可以更直观地观察变量的厚尾性。

图1 处理后的地震损失分布直方图和QQ图

表1 地震直接经济损失和死亡人数的描述性统计

从图2的散点图可以看出随着震级的增大，地震带来的经济损失和死亡人数也不断增大，但是并非单一的线性关系。在同一震级下，不同地震带来的损失之间的差异很大，仅研究震级与地震灾害损失均值之间的关系并不能很好地拟合和预测损失数据。地震巨灾损失数据的过离散和厚尾特征，使得传统的均值回归模型并不适用，而分位数回归模型在一定程度上可以克服均值回归模型的一些缺陷。图2给出了τ=0.5、0.75、0.95时处理后的二元损失变量与震级拟合的分位数回归线，从中可以发现，在各分位数水平下，处理后的二元损失变量与震级之间关系可以通过线性模型描述，因此分位数回归模型可适用于本文样本数据的研究。

图2 地震损失与震级的拟合分位数回归线

(二)参数估计结果

表2 不同分位数水平下的回归系数估计值

表2中直接经济损失对数的回归系数估计值为βτ11，转换后的死亡人数在不同分位数水平下的回归系数估计值为βτ21。结果表明，在0.1的显著水平下，各分位数处βτ11、βτ21是显著大于0的。这一结论验证了震级与直接经济损失之间的高度相关性，为以震级作为地震巨灾保险的触发指数的巨灾保险制度提供了实证支持。值得一提的是，随着分位数水平的不断攀升，βτ11虽有波动，但仍然存在一个上升的趋势，而βτ21本应该在0.99分位数处取得最大值，但表中并没有呈现出这样的规律，这可能是因为在不同分位数水平下，震级可解释的地震经济和生命损失是不同的。众所周知，对于地震损失而言，其受多种因素的影响，人口密度、居民防灾意识和地区经济水平均可影响地震灾害损失，在本文中，这些影响因素均被笼统地概括为常数βτ10、βτ20。可是，对于生命损失而言，通过对原始数据进行分析可以发现，我国地震灾害震级主要分布在4.5～6.5之间，云南、新疆等地区经常发生大地震，但由于人口稀少导致死亡人数较少，而对于四川和东部沿海各省份来说，由于人口密度上升，导致中等地震即可能造成重大伤亡。因此，除了震级之外，地震灾害死亡人数还受各地区人口密度的影响，并且人口密度是非常重要的一个影响因素。

对于相关系数ρ而言，表3给出了τ1=τ2=0.5，…，0.95时的估计值以及估计值对应的标准误。由估计结果可以发现，直接经济损失和死亡人数之间的相关系数显著大于零，即两者之间存在正相依关系。但随着分位数水平的不断上升，两者之间的相依性逐渐减弱，这是由于在高分位数水平处，死亡人数波动十分剧烈，波动幅度较大，减弱了其与经济损失之间的相关性。

表3 二元损失变量相关系数估计值

(三)地震灾害风险度量

在二元损失独立和相依的模型假设下，通过估计τ1=τ2=0.5，…，0.99时的回归系数，进一步测算直接经济损失在不同分位数水平下的VaR值(见表4)。结果表明，若假设直接经济损失和死亡人数相互独立，则VaR值会被低估，故为了测算得到更充足的保费规模，考虑二元损失的相依性是必要的。在二元损失相依的模型假设下，测算得到的VaR值随着震级和置信度的升高而增加，且在同一置信度下，VaR随震级的提高具有加速增加趋势。如在置信度取0.9时，4级地震的VaR为1.819亿，5级地震为11.225亿，增长9.406亿；6级地震为69.279亿，增长58.054亿；7级地震为427.557亿，增长358.278亿；另外，在同一震级下，VaR随概率水平增长而增加的速度也不断加快。

表4 直接经济损失的VaR值 (单位：亿元)

表5给出了在二元损失独立和相依两种情形下，死亡人数的VaR值。结果表明，在独立的情形下，死亡人数的VaR值同样会被低估。通过观察在二元损失相依的模型假设下测算的VaR值发现，当地震震级较低时，造成的生命损失比较小，如表中4、5级地震几乎不造成很大的人员伤亡。但随着震级的增大，生命损失会出现加速骤增，因此大地震造成生命损失不可小觑。

表5 死亡人数的VaR值 (单位：人)

表6给出了独立和相依情形下，处理后的二元损失在不同分位数水平下的ES值。与VaR相同的是，在模型假设为独立的情形下，测算的ES值也会被低估，再次证实考虑二元损失的正向相依性的必要性。在二元损失相依的模型假设下，取对数的直接经济损失的ES值随着τ1的增大而增大，即尾部风险极大。在对两变量均取对数的情况下，死亡人数ES值的变化率要高于直接经济损失的变化率。尤其是在τ2≥0.9后，ES值波动幅度非常大，尾部风险极高且难以把控，所以生命损失的风险管理非常值得探讨和研究。

表6 二元损失的ES值

(四)地震指数保险保费规模测算

前述结论已证明了当前地震巨灾保险以地震震级作为指数保险的触发条件的合理性。云南省大理州政策性地震指数保险将赔付触发点设置为5级，而通过对我国地震灾害损失数据进行分析发现，我国成灾地震中4.5～7级地震造成的灾害损失占所有地震灾害损失的比重约为90%，而4.5级以下地震形成的灾害损失仅占比不到6%，将触发点设置为4.5级可以涵盖90%以上的地震灾害，所以本文将触发点设置为4.5级。

图3 不同窗宽下的核密度估计

接下来根据τ1=τ2=0.5，…，0.99时回归系数的估计结果，计算4.5级以上的地震灾害经济损失的VaR均值。由于我国大陆1990—2020年地震灾害震级最大为8.1级，因此震级上限选取为8.1。计算公式如下：

(4)

本文运用随机模拟法求解式(4)，设定模拟次数为C=50 000，计算结果见表7。在90%的概率水平下，4.5级以上单次地震灾害的直接经济损失VaR值为73.776亿元。根据1990—2020年地震灾害数据可知，我国大陆4.5级以上地震灾害平均每年发生9.29次，因而预测若每年的总保费规模为685.4亿元，则该保费在90%的概率下能够满足地震灾害直接经济损失的赔付。

表7 不同概率水平下4.5级以上地震的直接经济损失VaR值 (单位：亿元)

在我国现行巨灾风险管理体系下，当地震灾害发生时，政府会对遇难者家属支付一定数额的死亡抚恤金，本节尝试测算我国每年地震死亡保险人数规模。利用式(4)可计算不同分位数水平下每次4.5级以上地震造成的最大可能生命损失。如表8显示，在50%的概率水平下，4.5级以上地震造成的死亡人数小于2人；在70%的概率水平下，4.5级以上地震造成的死亡人数小于8人；在90%的概率水平下，4.5级以上地震造成的死亡人数小于98人；在95%的概率水平下，4.5级以上地震造成的死亡人数小于472人。基于我国大陆平均每年发生地震灾害的次数预测可得，我国每年地震灾害造成的死亡人数在90%的分位数约为903人。

表8 不同概率水平下4.5级以上地震死亡人数的VaR值

五、结论与展望

评估地震灾害造成的生命和财产损失对地震保险的设计具有重要的现实意义。本文基于我国大陆地区1990—2020年306次成灾地震数据，建立了以震级为解释变量的财产损失与人身损失二元损失联合分位数回归模型。结果表明，震级对直接经济损失的影响较为显著，且其对高分位数处损失的影响更大，这一结论为以震级作为地震指数保险的触发参数提供了支持。由于死亡人数具有零膨胀的特性，较低震级造成的人员伤亡影响有限，但随震级升高，人身伤亡损失陡然上升。除了震级之外，地区人口密度也是影响地震灾害死亡人数的重要影响因素。本文利用多元非对称拉普拉斯分布中的相关系数矩阵来刻画二元损失之间的相关关系，对相关系数的估计结果表明两者之间存在正相关关系，特别是接近中位数区域的正相关关系较强。

结合非对称拉普拉斯分布的性质，本文计算了二元损失独立和相依情形下，我国地震灾害直接经济损失和死亡人数的VaR和ES。结果表明，在独立的模型假设下，VaR和ES值都会被低估，所以为了使得测算的保费规模更为充足，考虑二元损失的正向相依性是非常必要的。在相依情形下，本文计算了在不同震级和概率水平下，我国地震灾害损失的VaR和ES。结果表明，VaR的增长幅度会随着震级的增大而升高，所以大地震(7级以上)造成的损失是不可轻视的。在概率水平较高时，死亡人数的ES波动幅度非常大，即尾部风险极高且难以把控，因此生命损失的风险管理是非常值得深入探讨和研究的。考虑到数据的可获得性，本文仅选用了震级作为自变量。然而，还有很多其他因素也可能会影响地震的损失。后续研究可以考虑如何将这些数据纳入模型，从而得到更加全面的结论。