Erlang混合分布广泛应用于保险损失数据的建模, 在保险破产理论和保险损失数据的拟合中都有良好的表现.保险破产理论中, 当利用混合Erlang分布对保险损失的严重程度建模时, 通常关注的一些指标将有明确的解析式, 比如无限破产概率, 随机破产时刻的拉普拉斯变换等, 这方面的研究可参考文献[3, 17, 21, 29]; 近几年, 学者更多关注于将Erlang混合分布用于拟合保险实际损失数据, 得到了很多令人满意的分布性质, 比如分布函数和矩都有解析式, 使得相关的风险测度value-at-risk (VaR)和tail VaR (TVaR)比较容易计算. Verbelen^[30]等将双边截断引入Erlang混合分布, 计算了再保险合同的纯保费.类似研究见文献[9, 10, 19, 26, 30]等. Lee和Lin^[20]提出Erlang混合分布的多元形式, 多元混合Erlang分布保留了一元Erlang混合分布的大部分有用的分布性质, 同时建模相依性, 与copula方法相呼应.关于多元Erlang混合分布的研究见文献[2, 16, 31, 32]等.

混合模型的首要问题是混合序的确定, Lee和Lin^{[19, 30]}等都利用BIC来确定Erlang混合分布的序, Yin和Lin^[33]提出了一种新的iSCAD惩罚函数, 建立惩罚似然函数, 运用EM算法给出参数的估计, 同时给出了混合序的估计.然而值得注意的是: Erlang分布是轻尾的, 用它来拟合重尾数据时可能很难达到预期效果.其次, 尾部数据相应的权重一般都很小, 公式(3.11) 可以看出, 权重小于阈值$\lambda$的相应Erlang分布都被删除, 这不利用保留拟合尾部数据的Erlang分布.为解决这些问题, 本文引入极值理论拟合尾部数据, 建立Erlang极值混合模型.

极值混合模型广泛应用于各领域的数据分析中, 尤其在保险、金融、水文和环境科学等领域.在保险领域, 大额索赔在保险公司的风险管理和产品定价, 尤其是再保险产品的定价方面, 有不可忽略的意义.文献[4, 12, 13, 23, 27]等将极值理论引入到保险的风险管理中.为使数据的主体和尾部都拟合的很好, Behrens^[5]提出单一参数分布与一个极值分布的混合模型, Carreau和Bengio^[7]}讨论混合参数分布与极值分布的混合模型, 类似的文献[5, 6, 15, 22, 24]等给出多种极值混合模型. Lee^[18]等最早将极值混合模型引入到保险数据中, 但是所有这些混合模型都没有考虑混合序的确定.

本文建立Erlang混合分布与广义帕累托(GPD)分布的混合模型, 广义帕累托(GPD)分布用于拟合数据的尾部, 而Erlang混合分布用于拟合数据的主体, 这样即有Erlang混合分布的优点, 同时保留了极值理论的长处.引入iSCAD惩罚来估计混合Erlang分布的参数, Yin和Lin^[33]已证明参数和混合序的估计都有一致性.

首先给出Erlang分布的密度函数为

$ f(x;\gamma, \theta)=\frac{x^{\gamma-1}e^{-x/\theta}}{\theta^\gamma(\gamma-1)!}, $

其中$\gamma$是取值为正整数的形状参数(shape parameter), $\theta>0$是尺度参数(scale parameter).

将$m$个不同的Erlang分布以权重$\boldsymbol{\alpha}=(\alpha_1, \cdots, \alpha_m)$混合, 则Erlang混合分布的密度函数为

$ f(x;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)=\sum\limits_{j=1}^m\alpha_j f(x;\gamma_j, \theta). $

相应的分布函数为$F(x;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)$, 其中权重参数$\boldsymbol{\alpha}=(\alpha_1, \cdots, \alpha_m)$满足$\alpha_j\ge 0$和$\sum\limits_{j=1}^{m}\alpha_j=1$, $\boldsymbol{\gamma}=(\gamma_1, \cdots, \gamma_m)$是形状参数, 为了可识别性的说明, 一般有$\gamma_1\le\cdots\le\gamma_m$, 而$\theta>0$是共有的尺度参数.

由于投保人的性别、车型、驾车经验和熟悉程度等的不同, 使得索赔数据一般有明显的异质性, 单一的Erlang分布可能很难给出好的拟合效果, 因此数据的主体部分本文仍然选用Erlang混合分布来拟合, 而尾部采用极值分布.故本文采用左右双边截断的Erlang混合分布, 大部分保险损失数据都是已知截断值, 比如保险中的免赔额和赔偿限额.以$\textit{l}$和$\mu$分别表示左右截断值(免赔额$\textit{l}$已知), 双边截断的Erlang混合分布的密度函数是

$ f(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta) =\sum\limits_{j=1}^{m}\alpha_j\frac{f(x;\gamma_j, \theta)}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(\textit{l};\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)} $

(2.1)

$ \;\;\;\;\;\;\;\;=\sum\limits_{j=1}^{m}\alpha_j\frac{F(\mu;\gamma_j, \theta)-F(\textit{l};\gamma_j, \theta)}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(\textit{l};\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}\frac{f(x;\gamma_j, \theta)}{F(\mu;\gamma_j, \theta)-F(\textit{l};\gamma_j, \theta)}\nonumber\\ \;\;\;\;\;\;\;\;=\sum\limits_{j=1}^{m}\pi_j f(x;\textit{l}, \mu, \gamma_j, \theta)\\ \;\;\;\;\;\;\;\;\triangleq f(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}, \theta)\nonumber, $

(2.2)

其中

$ \pi_j=\alpha_j\frac{F(\mu;\gamma_j, \theta)-F(\textit{l};\gamma_j, \theta)}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(\textit{l};\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}. $

(2.3)

显然, (2.2) 式是左右截断点为$\textit{l}$和$\mu$的Erlang分布$f(x;\textit{l}, \mu, \gamma_j, \theta)$的混合模型, 混合权重为$\boldsymbol{\pi}=(\pi_1, \cdots, \pi_m)$, 满足$\pi_j\ge 0$和$\sum\limits_{j=1}^{m}\pi_j=1$.密度函数(2.1) 相应的分布函数为$F(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)$.

在统计中, 广义帕累托分布(Generalized Pareto Distribution, i.e., GPD)经常被用于拟合其他分布或实际数据的尾部, 本文选用GPD拟合数据的尾部, 其密度函数是

$ g_\mu(x;\xi, \sigma)= \begin{cases} \frac{1}{\sigma}\{1+\frac{\xi(x-\mu)}{\sigma}\}^{-\frac{1}{\xi}-1}, \mbox{$x\in(\mu, \infty)\;\;$ if$\xi>0$}, \\ \frac{1}{\sigma}exp\{-\frac{\xi(x-\mu)}{\sigma}\}, \mbox{$x\in(\mu, \infty)\;\;$ if$\xi=0$}, \\ \frac{1}{\sigma}\{1+\frac{\xi(x-\mu)}{\sigma}\}^{-\frac{1}{\xi}-1}, \mbox{ $x\in(\mu, \mu-\frac{\sigma}{\xi})\;\;$if $\xi<0$}. \end{cases} $

(2.4)

广义帕累托分布(GPD)的生存函数为

$ \overline{G}_\mu(x;\xi, \sigma)= \begin{cases} \{1+\frac{\xi(x-\mu)}{\sigma}\}^{-\frac{1}{\xi}}, \mbox{$x\in(\mu, \infty)\;\;$if $\xi>0$}, \\ exp\{-\frac{\xi(x-\mu)}{\sigma}\}, \mbox{$x\in(\mu, \infty)\;\;$if $\xi=0$}, \\ \{1+\frac{\xi(x-\mu)}{\sigma}\}^{-\frac{1}{\xi}}, \mbox{ $x\in(\mu, \mu-\frac{\sigma}{\xi})\;\;$if $\xi<0$}. \end{cases} $

(2.5)

结合(2.1) 和(2.4) 式, 为弥补引言中提过的Erlang混合模型的不足, 本文建立的Erlang极值混合模型的密度函数为

$ h(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta, \xi, \sigma) = \begin{cases} (1-\psi_\mu)f(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta), \;\;\mbox{if $x\leq\mu$}, \\ \psi_\mu g_\mu(x;\xi, \sigma), \;\;\mbox{if $x>\mu$}, \end{cases} \label{mixgpd} $

(2.6)

其中$\mu$为阈值, $X$是服从$h(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta, \xi, \sigma)$分布的随机变量, 令$\psi_\mu=P(X>\mu)$, 其一般由大于$\mu$的样本比例来估计.

相应的生存函数为

$ \overline{H}(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta, \xi, \sigma) = \begin{cases} \psi_\mu+(1-\psi_\mu)(F(u;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)\\ -F(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)), \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mbox{if $x\leq\mu$}, \\ \psi_\mu \overline{G}_\mu(x;\xi, \sigma), \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mbox{if $x>\mu$}. \end{cases} $

(2.7)

风险测度就是各种风险度量指标的总称.现行的国际标准风险管理工具VaR最初由Morgan针对银行业务风险的需要提出的, 并很快被推广成为了一种产业标准.风险价值VaR是指在正常的市场条件、给定的置信水平以及给定的持有期间内, 投资组合所面临的潜在最大损失. VaR是一种分位数风险测度, 一般给定置信水平$p$, 典型的$p=95\%$或者$99\%$.但是, VaR作为风险测度只考虑了概率为$p$的事件的最大损失$\operatorname{VaR}_{p}$, 高于$\operatorname{VaR}_{p}$的损失并没有纳入风险测度, 为克服这个缺陷, Tail Value at Risk (or TVaR)被提出来.在给定置信水平$p$下, TVaR就是损失落入最糟的$1-p$部分的平均损失.下面给出Erlang极值混合模型关于风险指标VaR和TVaR的计算.

为便于风险指标${\rm VaR}_p$和${\rm TVaR}_p$的计算, 当$X\leq\mu$时, 将生存函数(2.7) 用Erlang密度函数分别表达为

$ \overline{H}(x)=\int_x^\mu(1-\psi_\mu)f(x;l, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)dx+\int_\mu^\infty \psi_\mu g_\mu(x;\xi, \sigma)dx\\ \;\;\;\;\;\;\;\;\;=\psi_\mu+(1-\psi_\mu)\int_x^\mu\sum\limits_{j=1}^{m}\alpha_j\cdot\frac{1}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(l;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}\cdot\frac{x^{j-1}e^{-x/\theta}}{\theta^j(j-1)!}dx\\ \;\;\;\;\;\;\;=\psi_\mu+C\sum\limits_{j=1}^{m}\frac{\alpha_j}{(j-1)!}(\sum\limits_{k=0}^{j-1}(\frac{x}{\theta})^{k}\frac{(j-1)!}{k!}e^{-x/\theta}-\sum\limits_{k=0}^{j-1}(\frac{\mu}{\theta})^{k}\frac{(j-1)!}{k!}e^{-\mu/\theta})\\ \;\;\;\;\;\;\;=\psi_\mu+C\sum\limits_{j=1}^{m}\alpha_j(\sum\limits_{k=1}^{j}(\frac{x}{\theta})^{k-1}\frac{1}{(k-1)!}e^{-x/\theta}-\sum\limits_{k=1}^{j}(\frac{\mu}{\theta})^{k-1}\frac{1}{(k-1)!}e^{-\mu/\theta})\\ \;\;\;\;\;\;\;=\psi_\mu+C\sum\limits_{k=1}^{m}\sum\limits_{j=k}^{m}\alpha_j((\frac{x}{\theta})^{k-1}\frac{1}{(k-1)!}e^{-x/\theta}-\sum\limits_{j=k}^{m}(\frac{\mu}{\theta})^{k-1}\frac{1}{(k-1)!}e^{-\mu/\theta})\\ \;\;\;\;\;\;\;=\psi_\mu+C\sum\limits_{j=1}^{m}\sum\limits_{k=j}^{m}\alpha_k((\frac{x}{\theta})^{j-1}\frac{1}{(j-1)!}e^{-x/\theta}-\sum\limits_{k=j}^{m}(\frac{\mu}{\theta})^{j-1}\frac{1}{(j-1)!}e^{-\mu/\theta})\\ \;\;\;\;\;\;\;=\psi_\mu+C\theta\sum\limits_{j=1}^{m}Q_j(f(x;j, \theta)-f(\mu;j, \theta)), \nonumber $

其中$C=\frac{(1-\psi_\mu)}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(l;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}$, $Q_j=\sum\limits_{k=j}^{m}\alpha_k^*, j=1, \cdots, m$, 其中当$k=\gamma_j$, $\alpha_k^*=\alpha_k$, 否则$\alpha_k^*=0$.生存函数(2.7) 也可以表示为

$ \overline{H}(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta, \xi, \sigma) = \begin{cases} \psi_\mu+(1-\psi_\mu)\frac{\theta \sum\limits_{j=1}^{m} Q_j (f(x;j, \theta)-f(\mu;j, \theta))}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(\textit{l};\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}, \mbox{if$x\leq\mu$}, \\ \psi_\mu \overline{G}_\mu(x;\xi, \sigma), \mbox{if$x>\mu$}. \end{cases} $

(2.8)

假设损失随机变量$X$服从Erlang极值混合分布(2.6), 给定置信水平$p$, 有

$ \overline{H}(x;\textit{l}, \mu, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta, \xi, \sigma) =1-p. \label{VaR2} $

(2.9)

方程(2.9) 的解即置信水平为$p$的${\rm VaR}_p$.

计算${\rm TVaR}_p$之前, 首先研究自付责任额为$R(>\textit{l})$的再保险的纯保费, 当$R\leq\mu$时,

$ \begin{align} & \prod (R)=E({{(x-R)}_{+}})=\int_{R}^{\mu }{(1-{{\psi }_{\mu }})(x-R)}\cdot f(x;l, \mu, \mathbf{\alpha }, \mathbf{\gamma }, \theta )dx\ \\ & \ \ \ \ \ +\int_{\mu }^{\infty }{{{\psi }_{\mu }}}(x-R)\cdot {{g}_{\mu }}(x;\xi, \sigma )dx \\ & \ \ \ \ \ =(1-{{\psi }_{\mu }})[\int_{R}^{\infty }{(x-R)}\cdot f(x;l, \mu, \mathbf{\alpha }, \mathbf{\gamma }, \theta )dx-\int_{\mu }^{\infty }{(x-\mu )}\cdot f(x;l, \mu, \mathbf{\alpha }, \mathbf{\gamma }, \theta )dx \\ & \ \ \ \ \int_{\mu }^{\infty }{(\mu-R)}\cdot f(x;l, \mu, \mathbf{\alpha }, \mathbf{\gamma }, \theta )dx]\ +{{\psi }_{\mu }}\int_{\mu }^{\infty }{(x-R)}\cdot {{g}_{\mu }}(x;\xi, \sigma )dx \\ & \ \ \ \ \ =\frac{(1-{{\psi }_{\mu }})}{F(\mu ;\mathbf{\alpha }, \mathbf{\gamma }, \theta )-F(l;\mathbf{\alpha }, \mathbf{\gamma }, \theta )}[{{\theta }^{2}}\sum\limits_{j=1}^{m}{Q_{j}^{*}}(f(R;j, \theta )f(\mu ;j, \theta )) \\ & \ \ \ \ (\mu-R)\sum\limits_{j=1}^{m}{{{\alpha }_{j}}}\bar{F}(\mu ;j, \theta )]\ \\ & \ \ \ \ \ \ +{{\psi }_{\mu }}[(\mu +\frac{\sigma }{1-\xi })I(0<\xi <1)-R], \\ \end{align} $

其中$I(\cdot)$是示性函数, $Q_j^*=\sum\limits_{k=j}^{m} Q_k, j=1, \cdots, m$.

当$R>\mu$时,

$ \begin{align} & \prod{(R)}=E({{(x-R)}_{+}})=\int_{R}^{\infty }{{{\psi }_{\mu }}}(x-R)\cdot {{g}_{\mu }}(x;\xi, \sigma )dx \\ & \ \ \ \ \ \ \ ={{\psi }_{\mu }}[\int_{R}^{\infty }{{{(1+\frac{\xi (x-\mu )}{\sigma })}^{-\frac{1}{\xi }}}}\frac{1}{\xi }dx+(\mu-R-\frac{\sigma }{\xi }){{(1+\frac{\xi (R-\mu )}{\sigma })}^{-\frac{1}{\xi }}}]\\ & \ \ \ \ \ \ \ ={{\psi }_{\mu }}\frac{(R-\mu )\xi +\sigma }{1-\xi }{{{\bar{G}}}_{\mu }}(R;\xi, \sigma ). \\ \end{align} $

综上, 自付责任额为$R(>\textit{l})$的再保险的纯保费为

$ E((X-R)_{+})= \begin{cases} \frac{1-\psi_\mu}{F(\mu;\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)-F(\textit{l};\boldsymbol{\alpha}, \boldsymbol{\gamma}, \theta)}[\theta^2\sum\limits_{j=1}^{m}Q_j^* (f(R; j, \theta)\\-f(\mu; j, \theta))- (\mu-R)\sum\limits_{j=1}^m \alpha_j\overline{F}(\mu, j, \theta)]\\ +\psi_\mu[(\mu+\frac{\sigma}{1-\xi})I(0<\xi<1)-R], &\mbox{if$R\leq\mu$}, \\ \frac{\sigma+(R-\mu)\xi}{1-\xi}I(0<\xi<1)\psi_\mu \overline{G}(R;\mu, \sigma, \xi), &\mbox{if$R>\mu$}. \end{cases} $

(2.10)

当自付责任额$R={\rm VaR}_p$时, 置信水平为$p$的${\rm TVaR}_p$为

$ \begin{split} {\rm TVaR}_p &=E(X|X>{\rm VaR}_p)=E((X-{\rm VaR}_p)|X>{\rm VaR}_p) +{\rm VaR}_p\\ &=\frac{E((X-R)_{+})}{1-p}+{\rm VaR}_p. \end{split} \label{TVaR2} $

(2.11)

文献[33]针对每一个分量权重参数$\pi_j$, $j=1, \cdots, m$, 提出的iSCAD惩罚函数为

$ \begin{aligned} P_{\varepsilon, \lambda}(\pi_j)=&\lambda\{\log\frac{a\lambda +\varepsilon}{\varepsilon}+\frac{a^2\lambda^2}{2}-\frac{a\lambda}{a\lambda+\varepsilon}\}I(\pi_j> a\lambda)\\ &+\lambda \{\log\frac{\pi_j+\varepsilon}{\varepsilon}-\frac{\pi_j^2}{2}+(a\lambda-\frac{1}{a\lambda+\varepsilon})\pi_j\}I(\pi_j\leq a\lambda), \end{aligned} \label{iSCAD} $

(3.1)

其中$I(\cdot)$是示性函数.本文建立的Erlang极值混合分布中Eralng混合分布的参数估计与新引入的极值分布的参数估计互不影响, 因此关于Eralng混合分布的极大惩罚似然估计仍然是一致的.

Expectation--Maximization (EM)算法最早由Dempster^[11]给出比较详细的说明, 当似然函数的最大值点不能直接得到时, EM算法通过迭代的方法找到最大值点. EM算法需引入隐变量, 隐变量可以是未知参数, 丢失的数据或者任何可以使模型简化的未观测数据量. EM算法分为$E$-step和$M$-step两步, 其中$E$-step计算目标函数关于隐变量$\mathbf{Z}$的条件期望, $M$-step是最大化目标函数, 求得参数的极大似然估计.王继霞等^[1]将EM算法用于有限混合Laplace分布的估计.

Erlang极值混合模型的所有待估参数是:拟合数据主体部分的Erlang混合分布的序$m$, 形状参数$\boldsymbol{\gamma}=(\gamma_1, \cdots, \gamma_m)$, 相应的权重参数$\boldsymbol{\alpha}=(\alpha_1, \cdots, \alpha_m)$, 所有Erlang分布共用的尺度参数$\theta$, 拟合数据尾部的广义帕累托分布的阈值$\mu$, 尺度参数$\sigma$, 形状参数$\xi$, 下面逐一介绍它们的估计.

由公式(2.2) 知, 密度函数(2.6) 也可以由新权重参数$\boldsymbol{\pi}$表示为

$ h(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}, \theta, \xi, \sigma) = \begin{cases} (1-\psi_\mu)f(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}, \theta), &\mbox{if$x\leq\mu$}, \\ \psi_\mu g_\mu(x;\xi, \sigma), &\mbox{if$x>\mu$}. \end{cases} \label{mixgpd1} $

(3.2)

假设$\boldsymbol{x}=(X_1, \cdots, X_n)$是独立同分布的随机变量, 服从密度函数$h(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}, \theta, \xi, \sigma)$, 即(3.2), 样本观测值为$\boldsymbol{x}=(x_1, \cdots, x_n)$, 相应有序样本观测值为$x_{(1)}\leq\cdots\leq x_{(n)}$, 记

$ \begin{equation} k=\sum\limits_{i=1}^{n}I(x_{(i)}>\mu). \end{equation} $

(3.3)

Pickands^[25]给出与阈值$\mu$相应的$k$的选择方法, 从1开始依次增加, 最大值为$[n/4]$, 而$\mu=x_{(n-k)}$, 本文最终由似然函数的大小选出$\mu$.为方便后面的说明, 重新表示$n'=n-k$和$x'=(x_{(1)}, \cdots, x_{(n')})$.

形状参数的估计采用Yin和Lin^[33]类似的方法, 即预先给定一个大的混合序$M$, 形状参数的所有可能取值是$\boldsymbol{\gamma}^0=(\gamma_1^0, \cdots, \gamma_M^0)$, 通过估计相应的权重参数, 来实现混合序的估计和形状参数的选择.

Erlang极值混合分布的密度函数$h(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}^0, \theta, \xi, \sigma)$中的部分未知参数记为$\boldsymbol{\phi}=(\pi_1, \cdots, \pi_M, \theta)$, 本文采用EM算法来估计$\boldsymbol{\phi}$.

样本$\boldsymbol{x}=(x_1, \cdots, x_n)$的对数似然函数为

$ \begin{align*} \ell_n(\boldsymbol{\phi}, \xi, \sigma;\boldsymbol{x})&=\sum\limits_{i=1}^{n}\ln h(x;\textit{l}, \mu, \boldsymbol{\pi}, \boldsymbol{\gamma}^0, \theta, \xi, \sigma)\\ &=\sum\limits_{i=1}^{n}\ln[(1-\psi_\mu)\sum\limits_{j=1}^M\pi_jf(x_i;\textit{l}, \mu, \gamma_j^0, \theta)I(x_i\leq\mu)\\ &+\psi_\mu g_\mu(x_i;\xi, \sigma)I(x_i>\mu)]\\ &=\sum\limits_{i=1}^{n'}\ln[(1-\psi_\mu)\sum\limits_{j=1}^M\pi_jf(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)] +\sum\limits_{i=n'+1}^{n}\ln[\psi_\mu g_\mu(x_{(i)};\xi, \sigma)]\\ &=\ell_{n'}(\boldsymbol{\phi};\boldsymbol{x}')+\sum\limits_{i=n'+1}^{n}\ln[\psi_\mu g_\mu(x_{(i)};\xi, \sigma)]. \end{align*} $

样本$\boldsymbol{x}=(x_1, \cdots, x_n)$的iSCAD惩罚对数似然函数, 其中与参数$\boldsymbol{\phi}=(\pi_1, \cdots, \pi_M, \theta)$有关的部分是

$ \ell_{n', P}(\boldsymbol{\phi}; \boldsymbol{x})=\ell_{n'}(\boldsymbol{\phi}; \boldsymbol{x}')-n'\sum\limits_{j=1}^{M}P_{\varepsilon, \lambda}(\pi_j). $

(3.4)

直接关于$\ell_{n', P}(\boldsymbol{\phi}; \boldsymbol{x})$求极大似然估计是困难的, 本文使用EM算法, 引入隐变量, 即$\boldsymbol{Z}=(\boldsymbol{Z}_1, \cdots, \boldsymbol{Z}_n)$, 其中$\boldsymbol{Z}_i=(Z_{ij}|i=1, \cdots, n, j=1, \cdots, M)$,

$ Z_{ij}= \begin{cases} 1, &\mbox{如果观测值$x_i$来自第$j$个分量分布$f(x_i;\textit{l}, \mu, \gamma_j^0, \theta)$}, \\ 0, &\mbox{否则}. \end{cases} $

(3.5)

那么完整样本($\boldsymbol{x}, \boldsymbol{Z}$)的似然函数为

$ \begin{split} L_n(\boldsymbol{\phi};\boldsymbol{x}, \boldsymbol{Z})=&\prod\limits_{i=1}^{n}\prod\limits_{j=1}^{M}\{[\pi_j(1-\psi_\mu)f(x_i;\textit{l}, \mu, \gamma_j^0, \theta)I(x_i\leq\mu)]^{z_{ij}}\\ &+\psi_\mu g_\mu(x_i;\xi, \sigma)I(x_i>\mu)\}\\ =&\prod\limits_{j=1}^{M}\prod\limits_{i=1}^{n'}[\pi_j(1-\psi_\mu)f(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)]^{z_{ij}}\prod\limits_{i=n'+1}^{n}[\psi_\mu g_\mu(x_{(i)};\xi, \sigma)]. \end{split} $

(3.6)

相应完整样本$(\boldsymbol{x}, \boldsymbol{Z})$的对数似然函数为

$ \begin{split} \ell_n(\boldsymbol{\phi};\boldsymbol{x}, \boldsymbol{Z}) =&\sum\limits_{j=1}^{M}\sum\limits_{i=1}^{n'}z_{ij}\{\ln(\pi_j)+\ln[(1-\psi_\mu) f(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)]\}\\ &+\sum\limits_{j=1}^{M}\sum\limits_{i=n'+1}^{n}\ln[\psi_\mu g_\mu(x_{(i)}, \xi, \sigma)]. \end{split} $

(3.7)

相应的完整样本$(\boldsymbol{x}, \boldsymbol{Z})$的iSCAD惩罚对数似然函数为

$ \begin{equation} \ell_{n, P}(\boldsymbol{\phi};\boldsymbol{x}, \boldsymbol{Z})=\ell_n(\boldsymbol{\phi};\boldsymbol{x}, \boldsymbol{Z})-n'\sum\limits_{j=1}^{M}P_{\varepsilon, \lambda}(\pi_j). \end{equation} $

(3.8)

EM算法是利用迭代过程来估计参数的方法, 假设已经完成第$k$次迭代, 获得的当前估计是$\boldsymbol{\phi}^{(k)}=(\pi_1^{(k)}, \cdots, \pi_M^{(k)}, \theta^{(k)})$, EM算法的$E$-step和$M$-step分别为

$E$-step $\ell_{n, P}(\boldsymbol{\phi};\boldsymbol{x}, \boldsymbol{Z})$关于隐变量$\boldsymbol{Z}$求条件期望, 得到关于可观测样本$\boldsymbol{x}$的边际似然函数, 即

$ \begin{equation} \begin{split} Q(\boldsymbol{\phi} \mid \boldsymbol{\phi}^{(k)})=&\sum\limits_{j=1}^{M}\sum\limits_{i=1}^{n'}q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})\{\ln(\pi_j)+\ln[(1-\psi_\mu) f(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)]\}\\ &+\sum\limits_{j=1}^{M}\sum\limits_{i=n'+1}^{n}\ln[\psi_\mu g_\mu(x_{(i)}, \xi, \sigma)]- n'\sum\limits_{j=1}^{M} P_{\varepsilon, \lambda}(\pi_j), \end{split} \label{Q2} \end{equation} $

(3.9)

其中$q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})$是观测值$x_{(i)}$(i=1, \cdots, n')$来自第$j$个分量分布的概率,

$ \begin{equation} \begin{split} q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})=\frac{\pi_j^{(k)}[(1-\psi_\mu)f(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)]}{\sum\limits_{j=1}^{M}\pi_j^{(k)}[(1-\psi_\mu)f(x_{(i)};\textit{l}, \mu, \gamma_j^0, \theta)]}. \end{split} \end{equation} $

(3.10)

$M$-step (3.9) 式是权重参数$\pi_j$(j=1, \cdots, M)$和尺度参数$\theta$的函数, 求函数(3.9) 的极大估计, 即

$ \hat{\boldsymbol{\phi}}^{(k+1)}=\mathop {{\text{arg}}\;{\text{max}}}\limits_{\phi \in \Phi } Q(\boldsymbol{\phi} \mid \boldsymbol{\phi}^{(k)})\\ \;\;\;\;=\mathop {{\text{arg}}\;{\text{max}}}\limits_{\phi \in \Phi}\{\sum\limits_{j=1}^{M}\sum\limits_{i=1}^{n'}q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})\ln(\pi_j) +\sum\limits_{j=1}^{M}\sum\limits_{i=1}^{n'}q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})\\ \;\;\;\;\cdot[\ln(1-\psi_\mu)-x_{(i)}\diagup \theta-\gamma_j^0\ln(\theta)-\ln(F(\mu;\gamma_j^0, \theta)-F(\textit{l};\gamma_j^0, \theta))]\\ \;\;\;\;+\sum\limits_{i=n'+1}^{n}\sum\limits_{j=1}^{M}[\ln(\psi_\mu)+\ln g_\mu(x_{(i)};\xi, \sigma)]- n'\sum\limits_{j=1}^{M} P_{\varepsilon, \lambda}(\pi_j)\}. $

权重参数$\pi_j$的第$(k+1)$次迭代的估计为

$ \begin{equation} \begin{split} \hat{\pi}_j^{(k+1)} = \bar{q}_{j}^{(k)} I(\bar{q}_{j}^{(k)} > a\lambda)+ \frac{M}{\lambda}(\bar{q}_{j}^{(k)}-\lambda)_{+}I(\bar{q}_{j}^{(k)}\leq a\lambda), \end{split} \label{pi} \end{equation} $

(3.11)

其中$\bar{q}_{j}^{(k)}\triangleq\frac{\sum\limits_{i=1}^{n'}q(\gamma_j^0 \mid x_{(i)}, \boldsymbol{\phi}^{(k)})}{n'}$.

尺度参数$\theta$的第$(k+1)$次迭代的估计为

$ \begin{equation} \hat{\theta}^{(k+1)} =\frac{\frac{1}{n'}\sum\limits_{i=1}^{n'} x_{(i)} -t^{(k)}}{\sum\limits_{j=1}^{M}\gamma_j^0\bar{q}_{j}^{(k)}}, \end{equation} $

(3.12)

其中

$ \begin{equation} \begin{split} t^{(k)}=\left.\sum\limits_{j=1}^{M}\bar{q}_{j}^{(k)}\frac{\textit{l}^{\gamma_j^0}e^{-\textit{l}/\theta}-\mu^{\gamma_j^0}e^{-\mu/\theta}}{\theta^{\gamma_j^0-1}(\gamma_j^0-1)![F(\mu;\gamma_j^0, \theta)-F(\textit{l};\gamma_j^0, \theta)]}\right|_{\theta=\theta^{(k)}}. \end{split} \end{equation} $

(3.13)

迭代过程一直持续到$|Q(\phi^{(k)})-Q(\phi^{(k-1)})|$小于某个既定的误差界.分别以$\hat{\boldsymbol{\pi}}=\{\hat{\pi}_j| \hat{\pi}_j\neq 0, j=1, \cdots, M\}$和$\hat{\theta}$表示EM迭代的最终结果.混合模型序的估计是

$ \hat{m}=\#\{\hat{\pi}_j| \hat{\pi}_j\neq 0, j=1, \cdots, M\}. $

为便于说明, 重新将$\hat{\boldsymbol{\gamma}}=\{\gamma_j^0|\hat{\pi}_j\neq 0, j=1, \cdots, M\}$表示为$\hat{\boldsymbol{\gamma}}=(\hat{\gamma}_1, \cdots, \hat{\gamma}_{\hat{m}})$, 对应的权重参数记为$\hat{\boldsymbol{\pi}}=(\hat{\pi}_1, \cdots, \hat{\pi}_{\hat{m}})$.原权重参数的估计记为$\hat{\boldsymbol{\alpha}}=(\hat{\alpha}_1, \cdots, \hat{\alpha}_{\hat{m}})$, 由方程(2.3) 可得

$ \begin{equation} \hat{\alpha}_j =c \frac{\hat{\pi}_j}{F(\mu;\hat{\gamma}_j, \hat{\theta})-F(\textit{l};\hat{\gamma}_j, \hat{\theta})}, \end{equation} $

(3.14)

其中$c$是常数, 选择合适的$c$进行标准化, 满足$\sum\limits_{j=1}^{\hat{m}} \hat{\alpha}_j =1$.

最后, 关于广义帕累托分布(GPD)的尺度参数$\sigma$和形状参数$\xi$的极大似然估计, Coles^[8]已经详细讨论过, 本文就不再作重复说明.

本文利用R软件进行计算, 基于Yin和Lin^[33]关于Erlang混合分布的R程序和软件包“ismev”, 编写本文Erlang混合分布和GPD分布混合模型的R程序, 完成模拟实验和实际数据中模型参数的估计.

为验证模型和估计的有效性, 本文给出一个模拟实验, 从密度函数(2.6) 中随机抽取了2500个随机数, 其中(2.6) 式中的所有参数见表 1中的真实参数.

参数的初始化主要参考文献[19, 28].事先给定$M=10$, 形状参数的备择范围即$\boldsymbol{\gamma}=(1, \cdots, 10)$, 以Tijms^[28]的方法初始化, 公式(3.11) 给出极大惩罚似然的权重参数估计, 其稀疏性实现了在形状参数备择范围$\boldsymbol{\gamma}=(1, \cdots, 10)$中进行合理选择.从表 1可以看出, 形状参数最终仅选中$\hat{\boldsymbol{\gamma}}=(2, 7)$, 只有这两个形状参数对应的权重参数估计为非零的, 即$(\hat{\alpha}_2, \hat{\alpha}_7)=(0.501, 0.499)$, 其它形状参数相应的权重参数估计均为零, 即$\hat{\alpha}_j=0, j=1, 3, 4, 5, 6, 8, 9, 10$.显然, 混合模型序的估计$\hat{m}=2$.

由本实验可以看出, 引入iSCAD惩罚的优势所在:通过对权重参数的估计, 同时实现了对形状参数的选择和混合模型序的估计. 表 1列出的所有参数估计值与真实值都很接近, 说明模型和算法都很有效, 能够反映出数据的特征. 图 1很好的反应了这一点, 图 1中的真实曲线和拟合曲线几乎是重合的.

丹麦火灾赔偿数据有2167个观测值, Embrechts^[14]和Mendes^[24]等都用极值理论研究过这组数据的尾部, 本文采用Erlang极值混合模型从总体上研究这组数据, 不再仅仅限于研究其尾部特征.

文献[33]讨论了带左截断点$l$的Erlang混合分布, 本文在其基础上提出了Erlang极值混合分布, 在本例中将利用这两种不同的分布分别拟合丹麦火灾赔偿数据, 比较两种分布的优劣.

表 2给出Erlang混合分布和Erlang极值混合分布(2.6) 拟合火灾损失数据得到的所有参数的估计值, 其中利用Erlang极值混合分布得到的结果说明拟合数据的主体部分采用了三个Erlang分布, 数据的尾部由广义帕累托分布来拟合, 两部分的阈值点为4.174, 尾部数据比例为0.152;而利用Erlang混合分布拟合同一组火灾数据则需要十个不同的Erlang分布的混合.

图 2是丹麦火灾数据的直方图、Erlang混合分布和Erlang极值混合分布的拟合曲线, 可以看出拟合效果较好.

图 3和4分别给出Erlang混合分布和Erlang极值混合分布的Q-Q图, 显然Erlang极值混合分布在尾部数据的拟合上更优.

本文给出VaR的非参数(nonparametric)法估计作为标杆, 在置信水平为$p$的条件下, ${\rm VaR}_p$的非参数估计是方程$F_n({\rm VaR}_p)=p$的解, 其中$F_n(x)= \frac{\sum\limits_{i=1}^n I(x_i\leq x)}{n}$.

表 3给出三种方法的$\operatorname{VaR}_{p}$估计值, 表 3可以看出, Erlang极限混合分布估计得到的$\operatorname{VaR}_{p}$与非参数法得到的$\operatorname{VaR}_{p}$非常接近, 估计效果很好.

表 4给出非参数法、Erlang混合分布和Erlang极值混合分布的$\operatorname{TVaR}_{p}$估计值, 其中TVaR的非参数估计为$ {\rm TVaR}_p = \frac{\sum\limits_{i=1}^n (x_i\cdot I(x_i> {\rm VaR}_p))}{\sum\limits_{i=1}^n I(x_i> {\rm VaR}_p)}$. Erlang混合分布的$\operatorname{\rm TVaR}_{p}$比非参数法的结果偏小, 这主要是因为Erlang混合分布对火灾损失数据的尾部拟合不足, 见图 3; 而Erlang极值混合分布的结果稍大, 而且越到尾部, 这种趋势越明显, 这主要是因为估计得到的$\hat{\xi}=0.661>0$, 即估计的极值分布为厚尾的, 而实际数据的尾部过于稀疏, 不足以表现这种厚尾性.