中国是世界上受风暴潮灾害影响最为严重的国家之一。近年来, 由于全球气候异常变化, 中国风暴潮灾害的成灾频率不断升高, 造成的人员与经济损失也逐年递增。中国海洋灾害公报显示, 2015年风暴潮灾害造成的直接经济损失占海洋灾害总损失的99.8%。风暴潮灾害不仅会使水产养殖设施、码头、海堤、护岸、渔船等遭到严重损毁, 还会造成房屋损坏和人员伤亡, 给人民生命财产安全和沿海地区经济社会发展造成巨大威胁。因此, 为了合理度量风暴潮灾害风险、估计风暴潮灾害造成的巨大损失缺口, 科学测算风暴潮灾害的风险价值成为当务之急。

风险价值VaR(value at risk)作为一种较理想的风险度量指标, 是指一定的持有期内给定置信水平下可能发生的最大损失。它最早应用于金融领域的风险测度, 现在逐渐向其他领域扩展。VaR最明显的优势在于可以将风险量化为简单明了的数值。由于不同置信水平下的VaR值代表了不同概率下的最大可能损失, 可据此判断出单次灾害造成损失的大致范围。由于中国风暴潮灾害损失数据存在明显的厚尾性, 传统VaR将正态分布作为基本假设的作法可能会导致极端条件下VaR的低估, 干扰风暴潮灾害风险度量精度。鉴于此, 本文引入极值理论来修正传统方法存在的VaR估计缺陷。

1   材料与方法

1.1   研究背景

目前, 灾害风险度量的研究方法大致可以归结为以下几类:GIS技术分析法、信息扩散理论、超概率曲线法、建立风险评估模型法、模糊数学法等。李喜仓等^[1]、Kourgialas^[2]等学者都运用GIS空间分析技术对暴雨洪水等灾害风险进行分析, 形成风险区划空间分布图。Liu^[3]、Hao^[4]等学者基于信息扩散理论实现对森林火灾、草原生物灾害的风险评估。许闲等^[5]、张琳等^[6]等学者利用超概率曲线方法对地震、森林火灾进行灾害损失评估。魏章进等^[7]、胡俊锋等^[8]等学者选择构建风险评估模型来实现对台风、洪涝等灾害风险的综合评价。He^[9]、Gao^[10]等学者基于模糊数学方法最终确立了旱灾、洪灾和冰雪灾害等的风险等级。李晓等^[11]、李远平等^[12]等学者将层次分析法、自然灾害风险指数法、因子叠加分析、指标归一法等同GIS技术分析相结合, 绘制出风险等级的地理分布图, 实现对灾害风险的综合评定。极值理论作为次序统计学的分支, 在描述极端事件特征方面具有较强优势。极值理论具有超越样本数据的分析能力, 在总体分布函数未知的情况下, 仅依靠少量数据就可以实现极端数据特征的刻画。基于此, 近年来学者们开始将极值理论运用到自然灾害的风险度量中。耿贵珍和王慧彦^[13]通过POT模型分析中国地震风险损失分布, 得出GPD分布和POT模型在巨灾风险厚尾分布的拟合中具有巨大的优势和其他分布无法取代的地位。邵腾伟和冉光和^[14]将极值理论中的POT模型应用于估算农业自然灾害带给粮食产出的风险价值。

中国风暴潮灾害的风险度量研究主要集中在灾害风险评估模型的构建和灾害风险等级的划分等方面。郜志超等^[15]利用GIS强大的技术分析绘制出更精细的台风风暴潮灾害风险等级区划图。郑慧和赵昕^[16]通过将数量级作为灾害损失定级标准, 依据模糊数学原理构建隶属度函数, 得出了风暴潮灾害损失等级模糊定级测度法。齐鹏等^[17]将信息扩散理论引入风暴潮灾害的风险计算模型中, 得到最高潮位的超越概率估计曲线。李国胜和李阔^[18]通过构建多年一遇风暴潮综合风险评估模型, 实现对风暴潮风险的评估。

纵观现有研究成果, 大多数研究方法无法较好地解决数据厚尾性对风险度量造成的约束, 因此只能通过构建风险评估指标体系、制定风险等级划分标准、得出风险等级区划等角度进行灾害风险度量, 缺少基于真实损失分布度量风险的相关研究。中国风暴潮灾害的风险度量研究也多集中于风险等级划分, 缺少对灾害风险进行具体地量化测度。而学者们通过实证发现:极值理论超越样本的分析能力和在厚尾估计中的优势可以弥补传统方法估计VaR的不足, 实现基于真实损失分布的风险价值测算。因此, 本文采用极值理论的POT模型对风暴潮灾害直接经济损失数据进行建模, 并得出基于广义帕累托分布的损失分布函数, 从而测算出不同置信水平下的VaR。根据VaR测算结果, 得出未来单次风暴潮灾害可能造成的最大损失以及该损失发生的概率, 从而有效估算风暴潮灾害的损失缺口, 实现对我国风暴潮灾害风险的合理度量, 为今后有效分散风险提供科学依据。

1.2   理论模型构建

1.2.1   极值理论定义

极值理论主要包括BMM(block maxima method)模型和POT(peaks over threshold)模型。BMM模型是将所有数据进行分组, 然后选取每组最大值构成一个新的极值数组, 并以此为样本进行建模。POT模型则是对所有超过某一较大阈值的样本数据进行建模, 并在建模过程将超额分布近似为广义帕累托分布。本文使用POT模型完成风暴潮灾害直接经济损失总分布估计, 主要原因是:风暴潮灾害风险价值VaR测算依赖于其损失分布函数的分位数, 因此在较少风暴潮灾害损失数据的基础上, 如何较准确地拟合风暴潮灾害的损失分布尤为重要。Zheng等^[19]、Singh and Allen^[20]、林龙腾和沈利生^[21]等通过实证得出, POT模型在观测值较少的情况下估计准确性较高。因此, 本文运用POT模型来拟合中国风暴潮灾害直接经济损失总分布。

1.2.2   基于广义帕累托分布的POT模型

假设风暴潮损失样本数据为独立同分布的随机变量{X_i}(i=1, 2, ……n, n为选取的样本总量), {X_i}的分布函数记作F(x)。对于一个给定的充分大的阈值μ, 若X_i>μ, 则称X_i为超阈值, 称Y_i=X_i-μ为超出量。我们通常称超出量Y_i=X_i-μ的分布为条件超额分布函数, 即:

通过变换式(1)可得出用条件超额分布函数F_μ(y)表示的总分布函数F(x):

根据式(2)可知, 我们可以先通过求出超过阈值μ的条件超额分布函数, 然后得出风暴潮灾害的损失分布函数F(x)。

Pickands和Balkema and de Haan提出对于充分大的阀值μ, 条件超额分布F_μ(y)近似服从广义帕累托分布(Generalized Pareto Distribution, GPD)。因此根据Pickands-Balkname-De Haan定理可知, 若F∈MDA(H), 即当分布函数F属于极值分布H(x)的最大值吸引场时, 对于充分大的阈值μ, 条件超额分布函数F_μ(y)可近似为广义帕累托分布(GPD分布), 即:

根据G_{ξ, σ(μ)}(y)对样本损失数据尾部分布的拟合, 将其嵌入到式(2)中可得出:

本文采用历史模拟方法对F(μ)进行合理估计, 即:

式中:N_μ为超阈值样本个数。综合式(3)~式(5)最后得出风暴潮灾害的损失分布函数F(x):

1.2.3   基于POT模型的VaR推导

从统计学来看, VaR可确切地表示为:

式中:ΔX为未来一段时间内风暴潮灾害的潜在损失; q为置信水平。由此我们可以看出, VaR实际上是损失ΔX分布函数F(x)在q分位数下的点估计值。即

式中:F^-1(q)为损失分布函数F(q)的逆函数, 因此根据式(6)可推出:

2   结果与讨论

2.1   指标选取与数据分析

本文以中国1989~2014年间风暴潮灾害直接经济损失数据作为样本, 数据来源于《中国海洋灾害公报》(1989~2014), 运用R语言和Eviews 8.0软件来进行实证分析。为使数据更具有可比性, 本文采用GDP调节法对直接经济损失数据进行了调整, 即以2014年名义GDP为基准。样本损失数据描述性统计, 见图 1和表 1。

图  1  风暴潮灾害直接经济损失数据的直方图

Fig.  1  The histogram of direct economic loss of storm surge disaster

下载: 全尺寸图片 幻灯片

表  1  风暴潮灾害直接经济损失数据的描述性统计(单位:亿元)

Tab.  1  The descriptive statistics of direct economic loss of storm surge disaster

下载: 导出CSV 

| 显示表格

从图 1我们可以看出, 样本损失数据的直方图呈现出右偏态。同时, 表 1所示的描述性统计中峰度(Kurtosis)为24.83390, 远大于3;偏度(Skewness)为3.99270大于零。因此我们可以得出风暴潮灾害的直接经济损失数据具有厚尾特征。接下来, 对损失数据的厚尾性进行检测。

2.2   厚尾检测

本文采用以下方法进行厚尾检测。第一, QQ图。当数据独立同分布且满足正态分布时, QQ图中的散点近似为一条直线; 当数据分布为厚尾分布时, 则QQ图向上凸。第二, 经验平均超出函数图。根据McNeil(1997)^[22]的结论, 由平均超出函数e(μ)=E(X-μ|X>μ)的性质可知, e(μ)呈水平直线分布时, X服从某一指数分布; e(μ)有向上倾斜的趋势时, 表示X为厚尾分布; e(μ)有向下倾斜趋势时, 表示X为短尾分布。

图 2和图 3分别为样本数据的QQ图和经验平均超出函数图。通过QQ图可以看出, 风暴潮损失数据样本点呈现明显的上凸趋势, 表明损失数据具有厚尾性; 同时, 经验平均超出函数图中的各点呈近似线性分布且具有明显的向上倾斜趋势, 再次证明了损失数据的确具有厚尾性。因此确定风暴潮灾害直接经济损失分布为厚尾分布, 可以运用GPD分布进行拟合。

图  2  中国风暴潮灾害损失数据的QQ图

Fig.  2  The QQ plot of storm surge disaster loss data in China

下载: 全尺寸图片 幻灯片

图  3  中国风暴潮灾害损失数据的经验平均函数图(EMEF)

Fig.  3  The EMEF plot of storm surge disaster loss data in China

下载: 全尺寸图片 幻灯片

2.3   阈值选取

选择合理的阈值是POT模型拟合的关键。阈值μ的选取不仅关乎超阈值样本数据的数量还关系到GPD分布中参数的估计。阈值μ选取过大, 则会导致超阈值样本数据量过少, 影响模型拟合效果; 阈值μ选取过小, 超出量分布与GPD分布相差较大, 则会产生有偏估计。运用HILL图选取阈值的原则为选取图形中尾部稳定区域的起始点横坐标k所对应的数据X_k作为阈值μ。图 4为损失数据的HILL图。从图中可以看出, 中尾部首次达到平稳状态的k值约为67至77之间, 此时确定阈值μ的范围为(5.12959, 7.64503)。

图  4  HILL图

Fig.  4  The HILL plot

下载: 全尺寸图片 幻灯片

由于HILL图只能提供阈值的大致范围, 因此为了更加准确地选取阈值μ, 本文还结合参数稳定性估计图来进一步判定合理阈值。在HILL图确定的阈值取值范围内, 若μ₀为一个适当的初始阈值, 且超阈值样本服从广义帕累托分布, 则对于大于μ₀的μ, 相应的形状参数和尺度参数应保持稳定。图 5为参数稳定性估计图。图中的尺度参数σ和形状参数ξ在(5.15625, 5.84375)内均保持相对稳定, 为保证POT模型拟合的准确性, 我们选取该相对稳定区间内较大值作为阈值, 最终阈值μ=5.84375。

图  5  参数稳定性估计图

Fig.  5  The stability figure for parameters

下载: 全尺寸图片 幻灯片

2.4   参数估计与诊断

在得到阈值μ的基础上, 我们利用极大似然估计法(MLE)对GPD分布中的尺度参数σ和形状参数ξ进行估计。根据G_{ξ, σ(μ)}(y)求其密度函数, 进而得到对数似然函数如下:

对该似然函数求极大值便可估计出尺度参数σ与形状参数ξ, 计算结果如表 2所示。

表  2  极大似然估计法的参数估计结果

Tab.  2  The parameters estimation results using maximum likelihood estimation method

下载: 导出CSV 

| 显示表格

通过表 2可知, 尺度参数σ为20.86929, 形状参数ξ为0.29040, 当阈值为5.84375时, 超出阈值的样本数据为75个。将上述参数分析结果:μ=5.84375、Nμ=75、σ=20.86929、ξ=0.29040代入式(6)中可得中国风暴潮灾害直接经济损失分布函数为:

基于上述参数估计的结果, 进行风暴潮灾害直接经济损失数据GPD拟合分布的效果诊断。

从图 6可以看出, 概率图和分位数图上所有点几乎在一条直线上, 并且重现水平图中的点也都在重现水平的置信区间内。最后, 密度函数图中密度曲线与直方图的走势也相吻合。因此, 4个诊断图均证明参数估计结果具有一定的有效性, 并且POT模型拟合效果较好。

图  6  风暴潮灾害直接经济损失数据GPD拟合诊断图

Fig.  6  The GPD diagnosis charts ofdirect economic loss of storm surge disaster

下载: 全尺寸图片 幻灯片

2.5   基于POT模型的VaR估计与分析

根据可估算出基于POT模型的VaR值, 结果如表 3所示。

表  3  基于POT-GPD模型不同置信水平下的VaR测度(单位:亿元)

Tab.  3  The VaR in different confidence levels based on POT-GPD model

下载: 导出CSV 

| 显示表格

不同置信水平下的风暴潮灾害风险价值VaR是指风暴潮灾害在不同损失发生概率下的可能最大损失。通过表 3可知, 中国发生一次风暴潮灾害造成的直接经济损失不超过50.02860亿元的可能性为90%, 不超过75.90515亿元的可能性为95%, 同时有99%的把握认为发生一次风暴潮灾害造成的最大直接经济损失不会超过160.46587亿元。

为了确保风暴潮灾害VaR估测的准确性和有效性, 本文运用Kupiec失败频率检验法对VaR测算结果进行检验。设考察的实际损失天数为T, 失败天数为N, 则实际失败概率, 而期望的失败概率应为p^*=1-q(q为置信水平), 因此假定零假设为VaR失败概率为p^*, 构建统计量LR=, 并根据实证结果计算出不同置信水平下的LR统计量。若LR统计量小于其在该置信水平下的临界值, 则接受原假设, 从而证明该模型测算结果较好。检验结果如下表所示:

通过表 4可知, 不同置信水平下的风险价值VaR测算结果全部通过了Kupiec失败频率检验, 由此可以说明运用POT模型测算的VaR结果较准确, 与实际损失情况贴合程度较好。

表  4  Kupiec失败频率检验法

Tab.  4  The Kupiec failure frequency test

下载: 导出CSV 

| 显示表格

3   结论

(1) 极值理论的POT模型可以较好地描述中国风暴潮灾害的尾部损失分布, 弥补传统方法对VaR估算的缺陷, 实现合理测算中国风暴潮灾害风险价值的研究目标。

(2) 基于GPD分布的POT模型对中国风暴潮灾害直接经济损失历史观测值拟合, 得到我国风暴潮灾害直接经济损失的总分布函数为:

(3) 不同置信水平下基于POT模型的VaR值, 即当置信水平为90%、95%和99%时, 估测的VaR_q^POT值分别对应为50.02860亿元、75.90515亿元和160.46587亿元。根据上述VaR测算结果, 有99%的把握认为未来单次风暴潮灾害造成的最大可能损失额不会超过160.46587亿元, 为中国今后有效分散风暴潮灾害风险提供了科学依据。