标签归档:学习经历

COS访谈第十一期:郁彬教授

【COS编辑部按】:受访者:郁彬             采访者:施涛

原文刊登于ICSA。本文由COS翻译组策划翻译。译者是密西根大学的冷静、新加坡国立大学的尤晓斌和中国人民大学的霍志骥,全文最终由采访者施涛和被访者郁彬审核、修改、定稿,个别地方对英文原文作了补充。本翻译征得了ICSA、郁彬和施涛的同意和支持,在此表示诚挚的谢意。此外,陈丽云、高涛、肖楠、牟官迅、邓一硕、姜晓东、邱怡轩、魏太云对译文也提出了一些修正建议,在此一并表示感谢。

郁彬是加州大学伯克利分校统计系和电子工程与计算机科学系的Chancellor’s Professor。她曾在威斯康星麦迪逊和耶鲁大学都任过教,并且曾经是贝尔实验室的技术研究成员。她在2009年到2012年间担任加州大学伯克利分校统计系系主任,还是北大微软统计和信息技术实验室的创办者和主任之一。

她在顶尖的科学期刊上发表了70余篇论文,涉及统计、机器学习、信息论、信号处理、遥感、神经科学和网络研究等领域。她还在许多期刊中担任编委,比如统计年刊(Annals of Statistics)、美国统计学会会刊(Journal of American Statistical Association)、机器学习研究期刊(Journal of Machine Learning Research)和技术计量学(Technometrics)。

她是美国艺术与科学学院(American Academy of Arts and Science)的院士。2006年当选Guggenheim Fellow,2012年作了伯努利协会的图基纪念演讲(Tukey Memorial Lecturer)。她还是泛华统计协会2012年首届许宝騄奖的三位获得者之一。她也是AAAS(American Association for the Advancement of Science,美国科学促进会)、IEEE(Institute of Electrical and Electronics Engineers,电气和电子工程师协会)、IMS(Institute of Mathematical Statistics,数理统计协会)和ASA(American Statistical Association,美国统计协会)的会士。

她是IMS的主席。她担任过Statistical and Applied Mathematical Sciences Institute的国家科学委员会的联合主席,现在是Institute for Pure and Applied Mathematics的科学顾问组,以及布朗大学Institute for Computational and Experimental Research in Mathematics的执政委员会。

2013年2月13日,郁彬在她位于伯克利埃文斯教学楼的办公室里接受她以前的学生、现于俄亥俄州立大学任教的施涛的采访。以下是采访的全部内容。

继续阅读COS访谈第十一期:郁彬教授

COS访谈第六期:张健(微软)

zhangjian_photo【COS编辑部按】:受访人:张健,微软公司担任data scientist。

写在前面的话:前面小编采访了微软的数据科学家谢梁,当时同小编一同吃酒的还有微软的另一位数据科学家张健,巧在张健兄乃小编的师兄,毕业于Ames村办大学(又名爱荷华州立大学),当年我前脚到村,他后脚离村,所以之前也不认识。敝村可能名气不大,张健兄在统计界可能也不会有太多人知晓,但俗话说(好吧,我承认我瞎编的)“村长亦干部,凡夫即圣人”,小编很好奇一个物理博士在统计行当里捣鼓什么,于是发去这次采访,希望对外专业的同行们有所启示。
继续阅读COS访谈第六期:张健(微软)

从线性模型到广义线性模型(2)——参数估计、假设检验

1.GLM参数估计——极大似然法

为了理论上简化,这里把GLM的分布限定在指数分布族。事实上,实际应用中使用最多的分布就是指数分布族,所以这样的简化可以节省很多理论上的冗长论述,也不会限制实际应用。
如前文如述,指数分布族的概率密度函数可以统一地写为:

$f_Y(y;\theta,\Psi)=exp[(y\theta – b(\theta))/{\Psi} + c(y;\Psi)]$

这里为了在模型中体现散布参数(dispersion parameter)$\phi$,把上述密度函数中的$\Psi$记做
$\Psi=a_i(\phi)={\phi}/w_i$
从而响应变量的单个观测值的(加权)对数似然函数可以表示为:

$logL({\theta}_i,\phi;y_i)=w_i[(y_i{\theta}_i-b({\theta}_i))/{\phi}]+c(y_i,\phi)$

再结合观测值之间的独立性,全体观测值的对数似然函数可记做:$\sum_i logL({\theta}_i,\phi;y_i)$
一般情况下最大化上述的对数似然函数很难找到解析解(正态分布是特例之一),因而必须使用数值方法求解。McCullagh和Nelder(1989)证明了使用Newton-Raphson方法,结合Fisher scoring算法,上述对数似然函数的最大化等价于连续迭代的加权最小二乘法(iteratively weighted least squares, or IRWLS)。

广义线性模型的IRWLS算法如下:
1.设置线性估计量和响应变量的均值的初始估计值: $\hat {\eta}_0$和$\hat {\mu}_0$
这里$\hat {\mu}_0$是根据经验或是专家意见等信息对$\mu=E(Y)$的一个估计值,而$\hat {\eta}_0$可以利用模型建立时选用的联接函数来获得,即$\hat {\eta}_0=g(\hat {\mu}_0)$。这一函数关系也用于计算步骤2和3中$\eta$对$\mu$一阶导数。
2.构造调整的因变量(adjusted dependent variable):$z_0=\hat {\eta}_0+(y-{\hat \mu}_0){d\eta \over d\mu}|_{\hat {\eta}_0}$
3.构造权重:$w^{-1}_0={({d\eta \over d\mu})}^2|{\hat {\eta}_0V(\hat {\mu}_0)}$
这里$V(\hat {\mu}_0)$是利用方差函数(variance function)和$\hat {\mu}_0$构造的$Var(Y)$的估计值。
4.利用步骤2和3构造的调整的因变量和权重,拟合普通线性模型(ordinary linear model),预测/拟合(predict)新的线性估计量和均值: $\hat {\eta}_1$和$\hat {\mu}_1$
5.重复步骤2-4直到收敛(满足一定的迭代步数或是精度要求)。
此时得到的模型就是极大似然估计方法下的广义线性模型。IRWLS的算法思路也从另一个方面说明了广义线性模型是普通线性模型的推广。在广义线性模型的实际应用中,IRWLS算法是最常用的极大似然估计求解方法。对于特殊的案例,也有其他的特殊的参数估计方法。比如对于在精算学科中最常用的列联表(contigency table)数据或案例就有Bailey-Simon法、边际总和法(marginal totals)、最小二乘法(least squares)、直接法(direct method)等。

2.假设检验

2.1 空模型和全模型

一个极端的情况,所有自变量$x_i$对于响应变量$Y$都没有影响,也即是为所有的响应变量$Y$拟合一个共同的均值,即只有一个参数。这样的模型称为空模型(null model)。对于普通线性模型(正态分布下的GLM)而言,空模型的具体形式就是$y=\mu + \epsilon$。对于特殊的数据或案例类型,可能存在着其他的限制条件(constraints)从而空模型的参数个数大于1。比如非寿险精算中经常用到的列联表(contigency table)数据,其空模型就可能包含了行号、列号、对角线序号等限制。

相反的一个极端情况就是,所有自变量$x_i$的每一个观测值或称为数据的样本点(data points)对于响应变量$Y$都有影响,这样的模型称为全模型(full or saturated model)。一般可以通过构造阶数足够高的多项式或者把所有的量化观测值(quantitative)视为质化观测值(qualitive),并且引入适当数量的交叉项(interactions)来构造全模型。

统计建模的目的之一就是把样本数据划分为随机成分和系统成分两大部分。在这一点上,空模型认为响应变量的变动(variation)完全由随机性(random variation)造成,而全模型则认为响应变量的变动完全来自于系统成分(systematic)。一个直观地理解就是全模型是在现有的数据或样本的条件下,针对某一种分布所能拟合的最优模型,因而可以做为检验目标模型拟合优度的一个标准(measure)。

2.2 偏差(Deviance)

如果把全模型的对数似然函数记为$l(y,\phi|y)$,把目标模型的对数似然函数记为$l({\hat {\mu}},\phi|y)$,那么目标模型与全模型在拟合优度上的偏离的定义可写成$2(l(y,\phi|y)-l({\hat {\mu}},\phi|y))$。再结合观测值的独立性假设和指数散布族的假设,那么上述偏离的定义可以简化为:

$\sum_i 2w_i(y_i({\hat {\theta}_i} – {\tilde {\theta}_i}) – b({\tilde {\theta}_i}) + b({\hat {\theta}_i})) /{\phi}$

其中$a_i(\phi)={\phi}/w_i$,$\tilde {\theta}$是全模型下的参数估计值,$\hat {\theta}$是目标模型下的参数估计值。如果把上式写成$D(y,\hat {\mu})/{\phi}$,那么$D(y,\hat {\mu})$称为偏差(Deviance),$D(y,\hat {\mu})/{\phi}$则称为标准化偏差(scaled deviace)。
此外,皮尔逊卡方统计量(Pearson’s chi-square statistics):

$X^2={\sum_i (y_i – {{\hat \mu}_i})^2 \over Var({\hat {\mu}}_i)}$

也是衡量模型偏离程度(discrepancy)的统计量之一,在一些场合可以做为偏差的替代选择。

2.3 拟合优度检验

广义线性模型的假设检验可以分为两种:一是检验目标模型相对于数据或预测值的拟合有效性的检验(goodness of fit test);另外一种则是对“大”模型以及对“大”模型的参数施加一定的线性约束(linear restrictions)之后得到的“小”模型之间的拟合优度比较检验。直观上的理解就是,“大”模型具有更多的参数,即从参数的线性约束总可把一个或多个参数用其他参数的线性组合来表示,然后代入“大”模型,从而参数的个数减少,派生出所谓的“小”模型,也就是说“大”和“小”并非任意的,而是具有一种派生关系(nested models)。如果把全模型认为是“大”模型,而目标模型是“小”模型,那么上述两种检验的本质是相同的。因而假设检验的零假设(null hypothsis)可以统一且直观地设定为:“小”模型(目标模型)是正确的模型。

如果把大模型记做$\Omega$,把小模型记做$\omega$,其标准化偏差之差记做$D_{\omega} – D_{\Omega}$,其自由度之差记做$df_{\omega}-df_{\Omega}$,则构造如下的统计量:${(D_{\omega} – D_{\Omega})/(df_{\omega}-df_{\Omega})} \over {\phi}$。

当$\phi$是已知常数时,比如泊松和二项分布的情况下$\phi=1$,上述统计量在零假设下渐近地(asymptotically)服从卡方分布(正态分布时正好是卡方分布)。当$\phi$未知时,通常需要用估计值代替。最常用的估计值是$\hat {\phi}=X^2/(n-p)$这里n是数据中观测值的数量,p是目标模型的参数个数。此时上述的统计量在零假设下近似地(approximately)服从F分布(正态分布时严格服从F分布)。注意上述两种情况下,渐近和近似的区别。

对于某一个参数,可以使用其估计值的标准误(standard error)来构造一个z统计量来检验其显著性,即$z=\hat {\beta}/se(\hat {\beta})$。在零假设下,z统计量在普通线性模型,也就是正态分布下的广义线性模型中就是我们熟知的t统计量,严格服从t分布。在其他分布下的广义线性模型中,渐近地服从正态分布。z检验也称为Wald检验,在广义线性模型中效果不如上述的偏差检验,因而较少使用。

从线性模型到广义线性模型(1)——模型假设篇

在统计学里,对特定变量之间的关系进行建模、分析最常用的手段之一就是回归分析。回归分析的输出变量通常记做$ Y$,也称为因变量(dependent)、响应变量(response)、被解释变量(explained)、被预测变量(predicted)、从属变量(regressand);输入变量通常记做$ x_1$,…,$x_p$,也称为自变量(independent)、控制变量(control&controlled)、解释变量(explanatory)、预测变量(predictor)、回归量(regressor)。本文根据作者自己的一些学习心得和理解,简单且不严格地介绍在模型假设方面普通线性模型和广义线性模型的区别和联系/推广(generalization)。广义线性模型的拟合检验、推断、诊断等方面的方法和手段依赖于模型所采用的分布类型,难以一概而论,将在作者后续的学习心得文章里具体介绍。

1.普通线性模型的简单回顾

普通线性模型(ordinary linear model)可以用下式表示:

$ Y={\beta}_0+{\beta}_1x_1+{\beta}_2x_2+…+{\beta}_{p-1}x_{p-1}+\epsilon$                          (1.1)

这里$ {\beta}_i$,$ i=1$,…,$p-1$称为未知参数,$ {\beta}_0$称为截矩项。

普通线性模型的假设主要有以下几点:

1.响应变量$ Y$和误差项$\epsilon$正态性:响应变量$ Y$和误差项$\epsilon$服从正态分布,且$\epsilon$是一个白噪声过程,因而具有零均值,同方差的特性。

2.预测量$x_i$和未知参数${\beta}_i$的非随机性:预测量$x_i$具有非随机性、可测且不存在测量误差;未知参数${\beta}_i$认为是未知但不具随机性的常数,值得注意的是运用最小二乘法或极大似然法解出的未知参数的估计值$\hat{\beta}_i$则具有正态性。

3.研究对象:如前所述普通线性模型的输出项是随机变量$ Y$。在随机变量众多的特点或属性里,比如分布、各种矩、分位数等等,普通线性模型主要研究响应变量的均值$E[ Y]$。

4.联接方式:在上面三点假设下,对(1.1)式两边取数学期望,可得

$ E[Y]={\beta}_0+{\beta}_1x_1+{\beta}_2x_2+…+{\beta}_{p-1}x_{p-1}$                                        (1.2)

从 (1.2)式可见,在普通线性模型里,响应变量的均值$E[ Y]$与预测量的线性组合${\beta}_0+{\beta}_1x_1+{\beta}_2x_2+…+{\beta}_{p-1}x_{p-1}$通过恒等式(identity)联接,当然也可认为通过形为$f(x)=x$的函数(link function)联接二者,即

$ E[Y]=f({\beta}_0+{\beta}_1x_1+{\beta}_2x_2+…+{\beta}_{p-1}x_{p-1})={\beta}_0+{\beta}_1x_1+{\beta}_2x_2+…+{\beta}_{p-1}x_{p-1}$

2.广义线性模型的简单介绍

广义线性模型(generalized linear model)正是在普通线性模型的基础上,将上述四点模型假设进行推广而得出的应用范围更广,更具实用性的回归模型。

1.响应变量的分布推广至指数分散族(exponential dispersion family):比如正态分布、泊松分布、二项分布、负二项分布、伽玛分布、逆高斯分布。exponential dispersion family的详细定义限于篇幅这里不做详细介绍。

2.预测量$x_i$和未知参数${\beta}_i$的非随机性:仍然假设预测量$x_i$具有非随机性、可测且不存在测量误差;未知参数${\beta}_i$认为是未知且不具有随机性的常数。

3.研究对象:广义线性模型的主要研究对象仍然是响应变量的均值$E[ Y]$。

4.联接方式:广义线性模型里采用的联连函数(link function)理论上可以是任意的,而不再局限于$f(x)=x$。当然了联接函数的选取必然地必须适应于具体的研究案例。同时存在着与假设2.1里提及的分布一一对应的联接函数称为标准联接函数(canonical link or standard link),如正态分布对应于恒等式,泊松分布对应于自然对数函数等。标准联接函数的推导及其应用上的优点涉及到指数分散族的标准化定义,这里不做详述。

3.简单的例子

考虑这样一个简单的退保案例:一个保险产品一共卖出12份保单(当然了这在现实中不可能,这里仅为示例),在保单期限内一共有6人退保。那么采用这12个投保人的特征数据如收入、职业、年龄等做为预测变量对退保/退保率进行回归分析时,普通线性模型不再适用,因为这里退保这一事件不再服从正态分布,而是二项分布(当然了如果观测值的个数足够大,比如大于30,正态分布是一个很好的近似)。此时就可采用广义线性模型(目标分布采用二项分布)进行建模。

4.补充:指数分布族的简介

指数分布族(exponential dispersion family)实质上是对一类具有以下形式的概率密度函数或具有此类密度函数的分布的总括:

$f_Y(y;\theta,\Psi)=exp[(y\theta – b(\theta))/{\Psi} + c(y;\Psi)]$

这里$\Psi$和$\theta$是实参数,$b(.)$和$c(.;.)$是实函数,该密度函数的支集(support)$D_{\Psi}$是$R$的子集,且不依赖于$\theta$。满足$\theta=\eta=g(\mu)$的联接函数$g(\mu)$称为标准联接函数(standard or canonical link)。

一般情况下参数$\Psi$的值是未知常数(fixed and unknown),因此在许多GLM文献里指数分布族又被称为单参数指数族(one-parameter exponential family)。对于比较常用的分布,$\Psi$和$\theta$的取值具有特殊的形式:

正态分布$N(\mu,{\sigma}^2)$:$\Psi={\sigma}^2$和$\theta=\mu$

泊松分布$Poisson(\lambda)$:$\Psi=1$和$\theta=log\lambda$

二项分布$Binomial(m,p)$:$\Psi=1$和$\theta=log(p/(p-1))$

负二项分布$Negative Binomial(r,p)$:$\Psi=1$和$\theta=log(1-p)$

伽玛分布$Gamma(\alpha,\beta)$:$\Psi=1/{\alpha}$和$\theta=-{\beta}/{\alpha}$

R软件在精算教学中的应用案例

本文作者为张缔香,文章由COS编辑部审核发表,略有修改。点击此处下载/阅读本文PDF版本


R软件做为一种统计软件,因其开源、免费、灵活的诸多优点得到越来越多的关注,无论网络上还是实体书店,关于R的教程铺天盖地,不甚枚举。因此,本文的目标不是做R的教程,而是将R和保险、精算教学结合起来,通过几个案例来说明R在保险、精算专业日常的教学和研究中可用之处。

作者在保险、精算的理论、专业知识方面水平有限,也没有编程基础,使用R也只2个月左右的时间,所以以下几个案例,无论在理论还是R编程技巧方面都有所欠缺。但做为加入COS的首篇文章,当求真求实,不做文抄公。

最后,本文的几个案例,均来自于作者在南开大学保险系读硕期间的一些课程作业,主要目的在于完成某一项任务,而没有过多地考虑通用性、用户界面友好性等,比如所有代码中几乎都没有写入报错、警告语句等。当然,作者也假定本文读者不单熟悉R的常用语法、函数等,而且具有相应的统计学、保险、精算专业知识,案例中的术语、结论等都直接引用,不做解释。

案例一 简单寿险产品的保费、准备金等

《南开大学经济学系列实验教材》丛书里有一本李秀芳老师主编《寿险精算实务实验教程》,是李秀芳老师2010年秋学期主讲的《寿险精算实务实验》课程的指定教材,作者也选修了这一门课。因此本案例主要来自这门课的实验内容,但在实现方法上,不再使用换算函数法,而转向概率的思路,利用精算等价原理来求得寿险产品的净保费、毛保费、准备金等。

本案例中的R脚本文件命名为li.R,存放生命表的CSV文件命名为clt03.csv。

1. 寿险计价函数 li()

描述:以列表的形式返回目标寿险产品的毛保费、理论责任准备金、费用。

参数:

  • 投保人年龄age:单个数值型,取值范围0-105;
  • 生存给付bf.s:年初被保险人生存条件下的给付,非负数值型向量;
  • 死亡给付bf.d:被保险人死亡的年度末的给付,非负数值型向量;
  • 定价利率ir:按年计量的复利,数值型向量;
  • 费用率feerate:预定费用占毛保费的比例,非负数值型向量;
  • 保费费率prerate:均衡、递增、增减、或是其他形式的,数值型向量;
  • 定价生命表file:字符串型对象,存放生命表的.CSV文件的文件路径;
  • 产品类型type:取值范围{1、2、3、4},分别对应男性、女养老金业务,男性、女性非养老金业务。

注释:

1. 生存给付、死亡给付的设置方式

以保额1000元,20年定期死亡保险为例,那么死亡给付 bf.d=rep(1000,20),生存给付 bf.s=0;
若是保额1000元的20年定期年金,则 bf.d=0,bf.s=rep(1000,20);
若保险给付不是有规律的序列,则必须以 c(x,y,…) 的形式具体地给出。

2. 费用率、保费费率的设置方式

费用率、保费费率必须设置为数值型向量,且长度必须与目标寿险产品的缴费期限相同。
以10年缴费为例,假定第一年的费用率0.6,第二年0.4,第三年0.2,第四年及以后为0.1,那么feerate=c(0.6,0.4,0.2,rep(0.1,7))。若费用率为0,则必须指明为0 ,不可省略。
同样以10年缴费为例,若是均衡保费,即每一年的保费都相同,则保费费率 prerate=rep(1,10)。递增型 prerate=c(10:1) 表示第一年缴10单位保费,第二年9单位,以此类推。

3. 生命表 CSV 文件的格式

生命表部分数据如下所示,存放在 CSV 文件里:

age	 male	    female      malep	     femalep
0	0.000722	0.000661	0.000627	0.000575
1	0.000603	0.000536	0.000525	0.000466
2	0.000499	0.000424	0.000434	0.000369
3	0.000416	0.000333	0.000362	0.00029
4	0.000358	0.000267	0.000311	0.000232
5	0.000323	0.000224	0.000281	0.000195

4. 本案例中采用的生命表的范围是0至105岁,所以在设置年龄、生存、死亡给付、费用率、保费费率时,年龄和其它四个参数长度最大者之和不能超过106。
示例:

目标产品:20岁投保、20年死亡保险、定价利率2.5%、10年均衡缴费、男性非养老金业务的净保费

source(“li.R”)
li(age=20,bf.s=0,bf.d=rep(1000,20),ir=0.025,feerate=rep(0,10),prerate=rep(1,10),file="clt03.csv",type=1)

输出结果如下:

$premium
 [1] 1.623953 1.623953 1.623953 1.623953 1.623953 1.623953 1.623953
 [8] 1.623953 1.623953 1.623953
$reserve
 [1]  1.036481  2.059095  3.076657  4.096164  5.119730  6.148513
 [7]  7.183697  8.229508  9.282238 10.335088  9.710075  9.017654
[13]  8.244942  7.390757  6.447828  5.406646  4.256400  2.982906
[19]  1.569502  0.000000  0.000000  0.000000  0.000000  0.000000
$fee
 [1] 0 0 0 0 0 0 0 0 0 0

附录:附录一

2. 应用案例:

做为《寿险实务实验课程》的结课作业,作者设计了一款少儿保险产品,并用R和li() 函数计算净保费、毛保费、现金价值、法定责任准备金、利润测试等,由于篇幅过大,待整理完毕之后再发上来。

案例二 BMS奖罚系统

本案例的R脚本文件命名为bms.R,索赔次数的CSV格式数据文件命名为BMSclaim.csv。

1. 目标BM系统

本案例的目标BMS系统采用中国保险行业协会2007年4月制定的奖罚系统,具体如下所示:

表1 中国保险行业协会制奖罚系统表(2007.04)
规则 保费等级系数 等级
连续3年及3年以上无赔款 0.7 1
连续2年无赔款 0.8 2
上年无赔款 0.9 3
新保或上年发生3次以下赔款 1.0 4
上年发生3次赔款 1.1 5
上年发生4次赔款 1.2 6
上年发生5次及5次以上赔款 1.3 7

2. 转移概率矩阵

$
M=\left( \begin{array}{ccc} {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; p_0 \; 0 \; 0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; p_0 \; 0 \; 0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; p_0 \; 0 \; 0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; p_0 \; 0 \; 0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; 0 \; p_0 \; 0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; 0 \; 0 \; p_0 \\ {1-\sum_{i=0}^4p_i} \; p_4 \; p_3 \; p_1+p_2 \; 0 \; 0 \; p_0\\ \end{array} \right)
$

其中$p_k=\frac{\lambda^k}{k!}e^{-\lambda}$,表示投保人在一个保险年度中恰好发生 k 次索赔的概率。

3. 计算稳态分布

在稳定状态的处于各个级别的概率$\pi=(\pi_6,\pi_7,\dots,\pi_1)$满足:

$\sum_i\pi_i=1$, \ $\pi=\pi M$

解得:

$\pi_7=1-\sum_0^4p_i; \ \pi_6=p_4; \ \pi_5=p_3; \ \pi_4=p_2+p_1; \ \pi_3=p_0-p_0^2; \ \pi_2=p_0^2-p_0^3; \ \pi_1=p_0^3$

4. 利用索赔次数的实际数据求出具体的转移概率矩阵

本案例的索赔次数的观察值数据以如下所示,存放在CSV文件里。

claim	exposure
0	   27141
1	   5789
2	   1443
3	   457
4	   155
5	   56
6	   27
7	   2
8	   2
9	   1
10	   0

这里采用泊松逆高斯分布来拟合索赔次数分布,采用极大似然估计的方法来估计参数(在做这个案例时,作者对R了解甚少,很多的功能、函数其实是可以在各种专业package里找到的,这里不详述。)。结果可通过以下的命令查看:

source("bms.R")
gamma_pig
beta_pig

按此估计的参数值,恰好发生k次索赔的概率可通过以下的命令查看:

print(pr)

按以上的索赔概率值,和事先设定好的转移规则,可得到转移概率矩阵,通过以下命令查看:

tmat

为求稳定状态下处于各保费等级的概率值,将转移矩阵取30次方(实际上只需要5次方左右即可达到稳定的极限分布),这个极限矩阵的各行已经趋于相同,也即所求的稳定状态下的保费等级分布。可能通过以下命令查看极限分布:

tmat_lim

5. 计算评价BMS系统优劣的指标

5.1 稳定状态下的平均保费率(新保费率为1.0)

稳定状态下的平均保费率即是以稳定状态下处于各保费等级的概率分布为权重,各状态下的保费费率的加权平均值,计算结果为0.8196。可以通过以下命令查看:

avp
5.2 相对稳定状态平均水平(RSAL)

RSAL=(稳定平均水平 – 最低水平)/(最高水平 – 最低水平)=0.1993,可以通过以下命令查看:

rsal
5.3 对新司机的隐性惩罚(ECL)

ECL=(新保费率 – 稳定平均费率)/稳定平均费率=0.2201,可以通过以下命令查看:

ecl
5.4 变异系数COD=标准差/均值

各年的平均保费以及变异系数可以通过以下命令查看:

cod_mat

平均保费和变异系数逐年变化曲线可以通过以下的命令查看:

graph1()
5.5 风险区分度

由于对索赔频率的数据采用的是泊松逆高斯拟合的,所以每个保费等级内不同保单的索赔频率服从参数为$\lambda$的泊松分布(条件分布);在所有保费等级中$\lambda$又服从逆高斯分布,其参数的极大似然估计已在前文交待。

对于各保费等级的保单,其索赔频率和均值和方差即为以上的复合分布的无条件均值和方差。均值和方差的计算均要计算无穷积分。在R里,采用离散化的方式来近似计算,即采取足够大的区间、足够小的间隔,通过求和的方式来计算无穷积分的近似值。各保费等级的内的索赔频率和方差的估计值以及相应结构 的参数可以通过以下命令来查看:

fin_mat

各个保费等级下的索赔频率曲线图可以通过以下的命令查看:

graph2()
5.6 BMS的弹性

BMS的弹性公式为:$elac=(dp/p)/(d^{\lambda}/\lambda)=d \ln p /d \ln \lambda$,该公式为连续形式下的弹性公式,在实验中稳定状态下平均保费的函数形式的精确表达式难以计算出,只能采取离散化的方式来近期计算。具体思路就是把微分用差分代替,取间隔够小、数量够多的$\lambda$及对应的稳定状态下的平均保费,得到弹性的曲线图。可以通过以下的命令查看:

graph3()

代码:附录二

案例三 利用LeeCarter模型进行死亡率预测

本案例中的脚本文件命名为forecast.R,两个CSV文件分别命名为93pm.csv 和 03pm.csv。

1. 数据来源及初步整理

选用中国生命表90-93养老金男表和00-03养老金男表做为本案例的基础数据。把101至105岁划分为一组,100岁及其以上每一整数年龄分为一组,共102组。100岁及其以上组死亡率直接采用死亡率qx,100岁以上组采用中心死亡率。采用线性插值法得出1994-2002年间的死亡率。90-93养老金男表和00-03养老金男表各存放于一个CSV文件中,按如下图所示的格式存放数据:

在R的命令窗口输入以下命令来观察1993-2003 死亡率矩阵:

source(“forecast.R”)
tab9303

把死亡率取自然对数,然后减去每一年龄的死亡率的对数在1993-2003年间的均值,得到残差矩阵A。在R的命令窗口输入以下命令来观察均值和A:

ax
amat


各年龄的均值如上图所示。

2. SVD分解

2.1 k值

利用 R 里的 svd() 命令对死亡率的残差矩阵A做SVD分解。通过以下命令可观察分解得出的U、S、V矩阵:

svda$u
svda$d
svda$v

把S矩阵右乘V的转置矩阵,取其结果的第一行,即为Lee-Carter 模型中的kt。按1993-2003的顺序,kt的各个值如下:

-2.2891315, -1.9135574, -1.5181148, -1.1002924, -0.6570094,-0.1844099,
0.3224553, 0.8701444, 1.4677250, 2.1284910, 2.8736996

可以在R里输入以下命令观察kt的值及其总和:

kt
sum(kt)
2.2 b值

在SVD分解得出的U矩阵里,取其第一列,即可得到未经中心化的bx的值。本人采用保险研究增刊里的方法,采用线性拟合的方法来求bx的值。在R里输入以下命令可观察bx的拟合值以及拟合结果:

bmat

在bmat矩阵中其中,如各列名称所示,第一列为年龄,第二列为bx,第三至六列分别为线性拟合bx时的t值、t临界值、F值、F临界值和拟合优度值,可以看出所有的bx都是显著的。

2.3 k值的调整

采用最小化RMSPE的方法来得到最优的kt值,RMSPE的公式如下:

采用非线性最小值的算法来求解使得RMSPE最小的kt值,具体如下:

-2.3445627, -1.9458174, -1.5271026, -1.0869498, -0.6237423, -0.1356942

0.3791744, 0.9230672, 1.4984474, 2.1080823, 2.7550978

通过如下命令来观察求出的k值:

kt_adj
2.4 预测

对kt值和时间变量1-11(1993年做为1,2003做为11,或者直接用1993,1994,……2003区别只在系数的值,对预测没有影响),得出kt的向后15期的预测值如下:

3.044373, 3.551768, 4.059164, 4.566559, 5.073955, 5.581350, 6.088746, 6.596141, 7.103537, 7.610932, 8.118328, 8.625723, 9.133118, 9.640514, 10.147909

可以通过如下命令查看kt的15期预测值及其95%置信区间:

k.lm

使用以上的kt预测值来预测后15年的出生平均余命,

80.95256 ,81.35479 ,81.74651 ,82.12812, 82.49997 ,82.86240, 83.21572,83.56023,

83.89621, 84.22392 ,84.54362, 84.85554, 85.15991, 85.45694, 85.74685

下图为向后预测39期得出的出生平均余命的预测值,在2050年左右,中国男性的平均出生余命达到90岁。这似乎有些过高,所以本模型只适合做15年以内的预测。

代码:附录三

附录

[box]

R软件在精算教学中的应用案例(附录)(包含附录一、二、三)

[/box]