最后，本文的几个案例，均来自于作者在南开大学保险系读硕期间的一些课程作业，主要目的在于完成某一项任务，而没有过多地考虑通用性、用户界面友好性等，比如所有代码中几乎都没有写入报错、警告语句等。当然，作者也假定本文读者不单熟悉R的常用语法、函数等，而且具有相应的统计学、保险、精算专业知识，案例中的术语、结论等都直接引用，不做解释。
《南开大学经济学系列实验教材》丛书里有一本李秀芳老师主编《寿险精算实务实验教程》，是李秀芳老师2010年秋学期主讲的《寿险精算实务实验》课程的指定教材，作者也选修了这一门课。因此本案例主要来自这门课的实验内容，但在实现方法上，不再使用换算函数法，而转向概率的思路，利用精算等价原理来求得寿险产品的净保费、毛保费、准备金等。
本案例中的R脚本文件命名为li.R，存放生命表的CSV文件命名为clt03.csv。
**描述：**以列表的形式返回目标寿险产品的毛保费、理论责任准备金、费用。
参数：
注释：
以保额1000元，20年定期死亡保险为例，那么死亡给付bf.d=rep(1000,20)，生存给付bf.s=0；
若是保额1000元的20年定期年金，则bf.d=0，bf.s=rep(1000,20)；
若保险给付不是有规律的序列，则必须以c(x,y,…)的形式具体地给出。
费用率、保费费率必须设置为数值型向量，且长度必须与目标寿险产品的缴费期限相同。
以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单位，以此类推。
生命表部分数据如下所示，存放在CSV文件里：
age male femalemalep femalep0 0.000722 0.000661 0.000627 0.0005751 0.000603 0.000536 0.000525 0.0004662 0.000499 0.000424 0.000434 0.0003693 0.000416 0.000333 0.000362 0.000294 0.000358 0.000267 0.000311 0.0002325 0.000323 0.000224 0.000281 0.0001954.本案例中采用的生命表的范围是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.6239531.6239531.6239531.6239531.6239531.6239531.623953[8]1.6239531.6239531.623953$reserve[1]1.0364812.0590953.0766574.0961645.1197306.148513[7]7.1836978.2295089.28223810.3350889.7100759.017654[13]8.2449427.3907576.4478285.4066464.2564002.982906[19]1.5695020.0000000.0000000.0000000.0000000.000000$fee[1]0000000000**附录：**附录一
做为《寿险实务实验课程》的结课作业，作者设计了一款少儿保险产品，并用R和li()函数计算净保费、毛保费、现金价值、法定责任准备金、利润测试等，由于篇幅过大，待整理完毕之后再发上来。
本案例的R脚本文件命名为bms.R，索赔次数的CSV格式数据文件命名为BMSclaim.csv。
本案例的目标BMS系统采用中国保险行业协会2007年4月制定的奖罚系统，具体如下所示：
表1中国保险行业协会制奖罚系统表(2007.04)
$$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次索赔的概率。
在稳定状态的处于各个级别的概率\(\pi=(\pi_6,\pi_7,\dots,\pi_1)\)满足：
$$\sum_i\pi_i=1,\\pi=\piM$$
解得：
$$\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$$
本案例的索赔次数的观察值数据以如下所示，存放在CSV文件里。
laim exposure0 271411 57892 14433 4574 1555 566 277 28 29 110 0这里采用泊松逆高斯分布来拟合索赔次数分布，采用极大似然估计的方法来估计参数(在做这个案例时，作者对R了解甚少，很多的功能、函数其实是可以在各种专业package里找到的，这里不详述。)。结果可通过以下的命令查看：
source("bms.R")gamma_pigbeta_pig按此估计的参数值，恰好发生k次索赔的概率可通过以下的命令查看：
print(pr)按以上的索赔概率值，和事先设定好的转移规则，可得到转移概率矩阵，通过以下命令查看：
tmat为求稳定状态下处于各保费等级的概率值，将转移矩阵取30次方(实际上只需要5次方左右即可达到稳定的极限分布)，这个极限矩阵的各行已经趋于相同，也即所求的稳定状态下的保费等级分布。可能通过以下命令查看极限分布：
tmat_lim5.计算评价BMS系统优劣的指标5.1稳定状态下的平均保费率（新保费率为1.0）稳定状态下的平均保费率即是以稳定状态下处于各保费等级的概率分布为权重，各状态下的保费费率的加权平均值，计算结果为0.8196。可以通过以下命令查看：
avp5.2相对稳定状态平均水平(RSAL)RSAL=(稳定平均水平–最低水平)/(最高水平–最低水平)=0.1993，可以通过以下命令查看：
rsal5.3对新司机的隐性惩罚(ECL)ECL=(新保费率–稳定平均费率)/稳定平均费率=0.2201，可以通过以下命令查看：
ecl5.4变异系数COD=标准差/均值各年的平均保费以及变异系数可以通过以下命令查看：
cod_mat平均保费和变异系数逐年变化曲线可以通过以下的命令查看：
graph1()5.5风险区分度由于对索赔频率的数据采用的是泊松逆高斯拟合的，所以每个保费等级内不同保单的索赔频率服从参数为\(\lambda\)的泊松分布(条件分布)；在所有保费等级中\(\lambda\)又服从逆高斯分布，其参数的极大似然估计已在前文交待。
对于各保费等级的保单，其索赔频率和均值和方差即为以上的复合分布的无条件均值和方差。均值和方差的计算均要计算无穷积分。在R里，采用离散化的方式来近似计算，即采取足够大的区间、足够小的间隔，通过求和的方式来计算无穷积分的近似值。各保费等级的内的索赔频率和方差的估计值以及相应结构的参数可以通过以下命令来查看：
fin_mat各个保费等级下的索赔频率曲线图可以通过以下的命令查看：
graph2()5.6BMS的弹性BMS的弹性公式为：\(elac=(dp/p)/(d^{\lambda}/\lambda)=d\lnp/d\ln\lambda\)，该公式为连续形式下的弹性公式，在实验中稳定状态下平均保费的函数形式的精确表达式难以计算出，只能采取离散化的方式来近期计算。具体思路就是把微分用差分代替，取间隔够小、数量够多的\(\lambda\)及对应的稳定状态下的平均保费，得到弹性的曲线图。可以通过以下的命令查看：
graph3()**代码：**附录二
本案例中的脚本文件命名为forecast.R，两个CSV文件分别命名为93pm.csv和03pm.csv。
选用中国生命表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：
各年龄的均值如上图所示。
利用R里的svd()命令对死亡率的残差矩阵A做SVD分解。通过以下命令可观察分解得出的U、S、V矩阵：
svda$usvda$dsvda$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的值及其总和：
ktsum(kt)2.2b值在SVD分解得出的U矩阵里，取其第一列，即可得到未经中心化的bx的值。本人采用保险研究增刊里的方法，采用线性拟合的方法来求bx的值。在R里输入以下命令可观察bx的拟合值以及拟合结果：
bmat在bmat矩阵中其中，如各列名称所示，第一列为年龄，第二列为bx，第三至六列分别为线性拟合bx时的t值、t临界值、F值、F临界值和拟合优度值，可以看出所有的bx都是显著的。
采用最小化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值：
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年以内的预测。