1. 引言
近十年来,我国地震频发,汶川地震、雅安地震等地震灾害严重影响着广大民众的日常生活以及国民经济的健康发展。截至目前,我国的灾后经济重建主要依靠政府救济和社会捐助,与很多发达国家相比,我国的保险业赔付率非常低。长期以来,中国巨灾保险的发展举步维艰,为数不多的一些的积极尝试,多数以保险公司难以承担高额赔付而以失败告终。事实上,制约巨灾保险发展的主要原因在于我国缺乏有效的巨灾风险分散机制。而如何科学地评估巨灾风险,如何合理地分担巨灾风险,以及如何发挥金融市场在巨灾风险分散中的作用,这些是巨灾风险分散机制构建过程中需解决的重要问题。
现有研究表明,相对其他的厚尾分布,运用POT模型及帕累托分布对巨灾损失分布进行拟合的效果是最优的 [1] [2] ,而由于数据统计原因,我国地震灾害年损失数额值的数量目前较少,比如,郝军章、崔玉杰(2016) [2] 只选取了1961~2011年51个地震损失额损失。综合很多学者的研究 [2] [3] [4] [5] [6] ,可以看到,现有研究的中国地震巨灾风险评估数额不尽相同,尚未得到合理有效的结论。究其原因,小样本会对POT模型阈值的选取以及帕累托分布的拟合的准确性产生一定的影响。本文为了更好地刻画巨灾风险,利用了蒙特卡洛方法对数据进行扩充数据,试图避免由于样本的数据量较小对结果准确性的影响。在POT模型的运用过程中,阈值的选取是一个很重要的问题,耿贵珍、朱钰(2016) [3] 选取了较为直观的超额均值函数图与Hill图法,通过图像趋势进行判断,而大多数学者 [4] [5] [6] 都采用这两种方式或者其中的一种,进行阈值的选取。由于常用的超额均值函数图法与Hill图法较为主观,本文也同时运用了峰度法对阈值进行选取,试图得到更加客观的结果。
在巨灾债券方面,近年来,学者也做了一些研究。张节松和肖庆宪(2016) [7] 研究了巨灾保险公司的分保形式和自留风险的优化;程铖等(2014) [8] 使用Esscher变换为巨灾指数期权定价。可以看到国内相关研究已经取得了一些进展,但如何进一步合理描述巨灾风险的特性以及参数的选取是很重要的,这也正是巨灾债券在中国发展将要面临的问题,因而值得进一步探索。极值理论中的POT模型能够在总体分布未知的情况下,依靠样本数据推得总体极值分布的有关性质,故其在拟合巨灾风险厚尾部分的方面存在着相对优势。本文基于POT模型对地震巨灾风险进行了评估,更加合理地描述巨灾风险的特征,并依据所得的评估数据设计了一个一年期的零息地震巨灾风险债券,此债券考虑到一旦达到地震债券的触发值,债券发行方会面临较大的资金短缺的问题,同时也考虑到投资者不同的风险偏好,合理地对债券价格进行了分层,助力我国多层次地震巨灾风险保险体系的构建。
2. 极值损失统计模型
2.1. 极值理论
如果把总体分布函数记为,设为独立同分布的随机变量序列,个随机变量的最大值为,最小值为,则可得出最大值和最小值的分布函数分别为
在实际应用中,分布函数常常是未知的,如果先通过样本数据对其进行估计,再用进行极值分布的估计,那么估计中产生的细微误差将在中明显放大。目前可行的方式为承认未知,从出发,考虑的渐进模型,利用极值理论就解决问题。
2.2. POT模型的理论基础
广义帕累托分布(Generalized Pareto Distribution, GDP)被定义为:
式中:,是形状参数,是尺度参数。当时;当时;当时,广义帕累托分布为厚尾分布。
基于上面的广义帕累托分布,可以建立POT模型。
记为分布函数,设为独立同分布的随机变量序列,为阈值,记为超量损失,超量损失分布近似服从广义帕累托分布,其分布表示为:
对于足够大的阈值,超量损失分布函数可以用广义帕累托分布近似。即:
因此,超量损失分布函数可与广义帕累托分布近似。
阈值确定后,通常会使用历史模拟法对函数进行估计,即:记样本数为,记为超过阈值的样本数的个数,用近似表示。于是有:
3. 中国地震灾害损失概率的测算
3.1. 灾害损失数据的获得与处理
3.1.1. 数据来源与初步处理
查阅《中国地震统计年鉴》得到1990~2015年中国地震损失数额数据,共有26个数据,由于数据的时间跨度比较大,为了消除物价水平带来的影响,本文以2015年的CPI为定基指数,对1990~2015年的中国地震损失数额进行调整。经过物价调整后,得到以2015年为基期的1990~2015年的地震损失数据见图1。如图1所示,在2008年出现了明显的极端值状况,这是由于2008年发生了汶川大地震,为了消除这种大地震带来的影响,本文对数据进行对数处理,对数处理可以有效消除地震损失数据之间年度的巨大差异,处理结果如图2所示。
3.1.2. 扩充样本
本文选取了1990~2015年共26个数据,为了解决数据数量不足而导致模型的参数估计出现严重的偏差的问题,本文将通过模特卡罗模拟的方法扩大中国地震损失数据的样本空间。
首先,选取国内外相关研究中使用较多的3种厚尾分布模型作为候选模型 [9] ,分别为对数正态分布Lognormal、韦伯分布Weibull、伽马分布Gamma。之后,本文运用了参数方法对其进行拟合,利用极大似然估计法估计了相应分布的参数,见表1。接下来,需要进一步筛选哪一种分布拟合效果最优,本文选取了常用的一种拟合效果检验统计量,K-S检验法,如果得到的K-S检验统计量对应的P值较大,无法拒绝原假设,可认为拟合越接近原分布,结果如表1所示。
Figure 2. Logarithmic seismic loss value
图2. 经过对数处理的历年地震损失值
由表1的结果可以看到,Log-normal、Weibull、Gamma这三个分别对应的P值分别为0.9173、0.6282和0.8517,就K-S检验结果来说,在5%的显著性水平下,三个分布均通过了K-S检验,即认为数据拟合效果接近原始分布,三者比较,Log-normal分布的K-S检验的p值要高于其他分布的,可以认为在参数拟合法中,Log-normal分布是拟合地震损失数据的最优分布。同时,综合图3概率密度分布图与图4积累概率密度图来看,Log-normal分布的拟合效果是最好的,因此选取Log-normal分布来进行蒙特卡洛模拟。
本文根据已知的对数正态分布进行随机抽样1000次,即模拟1000次。由此产生已知的Log-normal分布的随机变量与之前的数据一起产生了新的样本,共1026个,如图5、图6所示。
3.2. 地震灾害数据的厚尾性检验
在构建模型之前,首先需要对地震损失数据进行厚尾性的检验。通常会将分布函数的尾部分布比指
Table 1. Maximum likelihood estimation of the corresponding parameter list and KS-test p table
表1. 极大似然估计法估计相应分布的参数表及KS-检验的p值表
Figure 3. The fitting effect and the probability density map of each distribution
图3. 各分部拟合效果对比图
Figure 4. The probability density map of each distribution
图4. 各分布积累概率密度图
数分布的尾部更厚的分布定义为厚尾分布,本文通过损失数据的概率密度函数图和关于标准指数分布的Q-Q图加以说明,如图7、图8。
从图7的概率密度函数图,可以看到损失数据顶部比正态分布更尖,尾部比正态分布更厚,呈现出明显的尖峰厚尾的性质,并且图像是向右偏斜的。通过图8关于标准指数分布的Q-Q图,可以看出上端向下倾斜,下端向上翘起,这表明数据的分布是尖峰态的,其尾部比正态分布的尾部厚。综上可见,用正态分布不能很好地拟合损失额数据,主要是极值数据的影响,因此,本文使用广义帕累托分布对这些极值数据进行拟合。
3.3. 阈值选取
阈值的选取是POT模型中非常关键的问题。阈值选择的准确性,对于参数和的估计起到重要的。如果阈值太大,会使得超量损失数据很少,将导致尺度参数偏高;如果阈值太小,会导致估计参数出现有偏性,无法达到应该的拟合效果。本文运用了三个阈值选取的方法:
3.3.1. 超额均值函数图
超额均值函数(Mean Excess Function, MEF)的定义为:
式中:表示超过阈值的样本观测量。MEF图是点的集合,横轴为,
Figure 7. Probability density function graph
图7. 概率密度函数图
Figure 8. Probability density function graph
图8. 概率密度函数图
纵轴为,得到相应的函数图。如果在某个观测值后函数曲线趋于线性,那么可以确定阈值为该观测值。图9为我国地震损失数据的样本平均超出量函数图及相应的95%置信区间,由图像可以看出阈值约在区间在(8,9)的一个范围内图形是近似线性的。
3.3.2. Hill图法
设表示独立同分布的次序统计量,可以得出:
横轴为,纵轴为,即Hill图为坐标的点构成的曲线,选取图中尾部指数稳定区域的起始点的横坐标对应的数值作为阈值。由图10的Hill图可以看到,尾部指数在超过大约31个数据后变得比较平稳,因此本文将阈值初步确定为8.5。
3.3.3峰度法
峰度的计算公式为
其中,
得到峰度后进行进一步计算,当时,将值最大的移出样本,循环该算法,直到计算出停止循环,从余下的观测值中选择最大的作为阈值。运用此方法确定阈值,结果为,接近于8.5。于是本文选取阈值,低于该阈值的地震损失属于正常范围,超过该阈值的地震损失额则认定为极值。
3.4. 参数估计及拟合效果检验
选定的阈值,运用极大似然估计方法进行参数估计,结果见表2。
利用估计结果,可以进一步得到拟合分布的诊断图。结果如图11所示,其中图a为P-P图,反映了各变量累计概率与其对应分布的累计概率之间的趋势关系,图b为Q-Q图,反映了所有数据变量分布的分位数和其对应分布的分位数之间的趋势关系。若P-P图和Q-Q图上的所有的点都近似在某一条直线上,则可以认为所选数据符合广义Pareto分布。图c为重现水平图,反映了重现期的对数数据与其对应的重现水平之间的关系,如图c中所示,所选的样本数据落在了指定分布的分位数置信区间内部,则可以认为所选取的数据符合广义Pareto分布。图d为尾部密度曲线的估计和直方图。从以上四个分布拟合诊断图可以看出,各数据点基本上分布在其参考线两侧,证明了对扩充过的地震年损失数据运用POT模型及广义Pareto分布拟合是合理的,同时拟合效果是很好的。
估计得出的广义帕累托分布函数为:
Table 2. Generalized Pareto distribution parameter estimation table
表2. 广义Pareto分布参数估计表
因此中国地震灾害损失概率分布函数为:
3.5. 巨灾风险水平的测算
高分位的估计式为:
由上述公式可以得到下表3。
从分位数的统计角度看,中国地震灾害损失额小于等于160,252.0万元的概率为90%,小于等于
Figure 11. Generalized Pareto distribution fitting diagnostic map: (a) P-P graph; (b) Q-Q graph; (c) reproducing horizontal graph; (d) histogram and density function estimation graph
图11. 广义Pareto分布拟合诊断图:(a) P-P图;(b) Q-Q图;(c) 重现水平图;(d) 直方图与密度函数估计图
Table 3. Generalized Pareto distribution parameter estimation table
表3. 广义Pareto分布参数估计表
512,894.9万元的概率为95%,小于等于9,281,176.0万元的概率为99%,小于等于35,305,717.1万元的概率为99.5%,小于等于982,081,246.9万元的概率为99.9%。
4. 地震巨灾债券的设计
4.1. 地震巨灾风险分担层次的划分
根据我国的国情,适宜建立政府主导的巨灾保险制度。因此,本文设计如下的地震巨灾风险损失的分担层次:
1) 地震灾害造成的0~16.025亿元之间的直接经济损失,可以通过我国国内的各大保险机构直接提供保险进行承担,根据相关数据,16.025亿元相对于2015年全国非寿险公司保费总收入的7994.97亿元而言,仅仅占了0.2004%,所以说国内保险市场完全有能力化解这个程度的风险;
2) 地震灾害造成的16.025~51.289亿元之间直接经济损失,可以通过国内主要的再保险市场使用分保等方式将部分风险和损失进行转移;
3) 地震灾害造成的51.289~928.118亿元之间的直接经济损失,可以通过成立地震巨灾保险基金,将地震巨灾风险证券化,在证券市场筹集资金,以期权等金融衍生品等形式将巨额风险打包转移到国内及国外的资本市场,以承担巨灾风险带来的巨额损失;
4) 地震灾害造成的超过928.118亿元的直接经济损失,可以通过政府以国家财政救援的方式为巨灾损失作最后的买单。
如图12所示。
4.2. 地震巨灾债券的定价模型
巨灾债券是通过发行收益与指定的巨灾损失相连结的债券,将保险公司部分巨灾风险转移给债券投资者。在资本市场上,需要通过专门中间机构(SPRVS)来确保巨灾发生时保险公司可以得到及时的补偿,以及保障债券投资者获得与巨灾损失相连结的投资收益。巨灾债券是指债券公开发行后,未来债券本金及债息的偿还与否,完全根据巨灾损失发生情况而定。即买卖双方通过资本市场债券发行的方式,一方支付债券本金作为债券发行的承购,另一方则约定按期支付高额的债息给另一方,并根据未来巨灾损失
Figure 12. Design scheme of earthquake catastrophic risk diversification mechanism
图12. 地震巨灾风险分担层次设计图
发生与否,作为后续付息与否及期末债券清偿与否的根据。
巨灾风险债券的有效运行有赖于设计合适的期限结构和选取合适的触发机制,这涉及利率是否随机、触发条件是只由单一事件触发还是包含跨期多事件等问题,而常用的触发机制包括实际损失触发机制、指数触发机制和混合触发机制等多种 [5] 。在本文的地震巨灾风险研究中,地震风险是以年为周期,因此如利用随机利率模型设计跨期多事件触发等复杂的巨灾风险债券,显然是不太合适的。而触发机制也难以简单划分到上述三种机制当中,根据本文采用的地震年损失数据,触发条件应为地震一年内造成的累计损失。
鉴于以上原因,本文决定选择结构较为简单的一年期零息债券来设计地震巨灾风险债券,其运作机制为:某地震巨灾风险债券发行方通过发行相关债券筹集一定规模的资金,若在债券期限内,地震的累计损失率达到或超过约定触发值,则在债券到期时,发行方有权没收债券面值一定比例的资金;若没有达到触发条件,则在债券到期时,发行方按照债券面值向投资者支付金额,因此投资者在债券到期日可获得的资金可表示为:
其中,T为债券到期时间,为债券到期偿还值,F为债券面值,D为触发值,A为在触发条件发行方的面值偿还比例,为损失率指数,I为指示函数,表示当时I取值为1,否则为0,表示当时I取值为1,否则为0。根据零息债券的一般定价方法,该地震巨灾风险债券的定价公式为:
由该公式可以看出,在无套利机会的条件下,地震巨灾风险债券的价格受到债券的相关参数:债券期限T、票面价值F、触发条件下的面值偿还比例A、触发值D、损失率指数的分布以及预期收益率r的影响。其中,由于本文设计的是一年期的零息债券,因此T = 1;根据本文设定的地震巨灾风险分担层次与之前求得的损失率分布函数,本文选择的触发值为5%,即51.289亿元;而损失率指数的分布已由前述中的中国地震灾害损失概率分布函数确定。对于其它影响因素,本文结合实际情况,选取了合适的参数范围以计算不同参数值组合下地震巨灾风险债券的价格。
参数设置如下:
1) 为了计算的方便,选取的债券面值F为10000;
2) 由于地震巨灾风险债券与其它金融产品的相关性较低,从而流动性较差,因此投资者一般要求较高的预期回报率,同时,综合考虑到投资者不同的风险偏好,应当对债券价格进行分层,本文设定的债券收益率r为5%、10%和15%;
3) 由于一旦发生地震巨灾,债券发行方会面临较大的资金短缺问题,因此触发条件下的面值偿还比例一般不会太高,考虑到债券发行方,本文的偿还比例A选择0、0.25和0.5;
不同参数值组合下的地震巨灾风险债券价格如表4所示。
从表4的数据可以看出,地震巨灾风险债券价格与各个参数间的关系为:在面值偿还比例一定时,债券价格随着收益率的上升而下降,这与传统的债券定价理论是一致的;第二,在收益率一定时,债券价格随着面值偿还比例的上升而上升,因为面值偿还比例上升会导致发行方在固定触发条件下的支付上升,投资者的回报上升,因此会使债券价格上升。
5. 总结
本文选取了1990~2015年的地震年损失数据,为了避免小样本对估计结果产生影响,运用了多种厚
Table 4. Generalized Pareto distribution parameter estimation table
表4. 广义Pareto分布参数估计表
尾分布对数据进行拟合,选取拟合效果最优的分布并运用了蒙特卡洛方法对样本进行扩充,从而对中国地震巨灾风险进行了更加合理的测算。接下来,引入极值理论中的POT模型,并选取了峰度法,综合利用常用的超额均值函数图法以及Hill图法,选取了合适的阈值8.5。之后,运用帕累托分布对新的样本空间进行了拟合,并对模型的形状参数以及尺度参数进行估计,拟合效果良好。下一步,本文根据测算的相关结果,对地震巨灾风险损失进行了简单的分担分层,其中考虑到了保险公司的担负能力。本文进一步根据测算的相关结果,设计了一个一年期的零息地震巨灾债券,一是基于POT模型测算出的相关数据更加切合地震巨灾风险的特点;二是考虑到一旦达到地震债券的触发值,债券发行方会面临较大的资金短缺的问题;三是考虑到利益投资者不同的风险偏好,对债券价格进行了分层,为债券发行方提供参考。本文在以往基础上做如下创新:1) 利用蒙特卡洛模拟方法实现了数据的扩充,避免了由于样本量过小导致的结果的偏差,从而增强了结果的准确性;2) 在阈值的选取过程中,除了常用的超额均值函数图与Hill图法,还引入了峰度法,避免由于过于主观的判断对阈值选取的影响,从而使结果更加准确;3) 给出了地震巨灾债券的定价模型,定价过程增加了巨灾风险分担和投资者风险偏好的考量,使定价结论更加契合市场需求。从而,对比以往的研究,本文对中国地震灾害损失进行了更加合理的测算,为地震巨灾风险估计的建模及巨灾债券的定价提供了更加准确的理论依据。
本文的研究结论将为保险公司、再保险公司和资本市场和国家财政各方合理分担巨灾风险提供理论框架,助力我国多层次地震巨灾风险保险体系的构建。
基金项目
辽宁省社会科学规划基金项目(L16BGL012);辽宁省教育厅科学研究项目(L2014024)。