Markov链随机利率下寿险精算函数的分布模拟

张连增，段白鸽，卜 林

0 引言

当前，随着计算技术的快速发展和统计软件的普及，为数值模拟现值随机变量的概率分布提供了技术支持。一方面，在多数情况下，现值随机变量的解析计算比较复杂，不如随机模拟方法直接方便。而且在实务中，对大多数寿险和生命年金产品的评估（如对现值随机变量的概率分布的研究），可以预期考虑解析解是相当困难的，甚至是不可行的。另一方面，随机模拟结果也可以作为验证解析计算的一种途径。鉴于此，Goovaerts(2000)[1]研究了随机模拟在现值随机变量精算函数中的应用。目前，随机模拟已经成为金融工程与统计研究的常用方法之一，在精算学中必然也将得到广泛的应用。为此，本文将在Markov链随机利率下，借助于随机模拟的方法来得到三类有代表性的精算函数的数值解。而Markov链作为最基本的随机过程，已经广泛应用于工程、经济科学、社会科学等领域，它具有易于理解、符合直观、编程方便、自动处理变量的不独立性问题等优点。另一个值得一提的优点在于，由Markov链模拟的随机利率不会取负值或较大值，这弥补了其它一些随机利率模型的不足之处。

本文将给出精算学中三个简单且有代表性的数值实例：20年期年金现值随机变量与积累值随机变量、20年期两全保险的现值随机变量、20年期期初付生命年金的现值随机变量，并借助于当前国际上日益流行的R软件编程实现。解决这些例子的思路与方法可以推广到精算学与金融学的进一步研究中。

1 理论模型

1.1 Markov链

关于Markov链的内容几乎在所有的随机过程教材中都有介绍，这里只给出最基本的介绍，参见Grinstead和Snell(1997)[2]。考虑离散时间的有限状态Markov链，记为{Xn:n=0,1,2,…}，状态空间为{s1,s2,…,sn}，转移概率矩阵可以表示为：

假设矩阵P对应于特征值为1的特征向量是一维的，则该Markov链存在唯一的平稳分布。如果记平稳分布为π=(π1,…,πn)，那么它满足如下方程：

平稳分布的含义可以表述为：如果设X0的概率分布为π，那么对任何n≥1，Xn的概率分布也为π。

另外，在某些条件下，Markov链存在极限分布，即当n→∞时，Xn的分布收敛于一个极限，该极限与初始状态X0无关。而且在很多情况下，极限分布与平稳分布是相同的。

1.2 确定性年金的现值与积累值

考虑期末付年金。对确定性n年期的年金，如果假设在时刻{1,2,…,n}的支付额分别为{b1,b2,…,bn}，在每年内的利率分别为{i1,i2,…,in}[9]，那么该年金在时刻0的现值可以表示为：

为了体现出n个利率{i1,i2,…,in}的影响，在计算年金积累值时，考虑期初付年金。如果假设在时刻{0,1,…,n-1}的支付额分别为{b0,b1,…,bn-1}，那么该年金在时刻n的积累值可以表示为：

1.3 两全保险

对于(x)岁的n年期两全保险，在每年内的死亡给付可以不同，而且死亡给付也可以与满期给付不同。为方便起见，假设死亡给付、满期给付都是1个单位。对给定的利率轨道{i1,i2,…,in}，该两全保险的现值随机变量Z的概率分布可以表示为：

其中，npx、n-1|qx等是寿险精算学中标准的生命表函数（生存概率、死亡概率）。

1.4 定期生命年金

对于(x)岁的n年期期初付生命年金，在每年初的生存给付可以不同。为方便起见，假设给付都是1个单位。对给定的利率轨道{i1,i2,…,in}，该定期生命年金的现值随机变量Z的概率分布[3]可以表示为：

其中，npx、n-1|qx等是寿险精算学中标准的生命表函数（生存概率、死亡概率）。

2 在Markov链随机利率下分布模拟的思路与方法

2.1 产生利率轨道

这里分两种情况展开讨论。第一种情况假设初始利率为常数，第二种情况假设初始利率的概率分布为Markov链的平稳分布。在第二种情况下，每个周期内的利率变量是同分布的（但并不独立）。对每一条利率轨道的产生，都应用Markov链的定义。利率轨道条数视两种不同的情况而定，通常对第一种情况1000条轨道即可，对第二种情况10000条轨道足够多了。

2.2 确定性年金的现值与积累值

针对每条利率轨道，计算确定性年金的现值[4]与积累值。对所有的利率轨道，重复计算后得到一系列年金现值与积累值的样本数据，进而通过对这些样本数据特征的描述便可得到现值随机变量和积累值随机变量的概率分布。

2.3 两全保险

由中国人寿保险业经验生命表（2000～2003），得到各个年度内的死亡概率，以及保险期限届满时的生存概率。针对每条利率轨道，得到一个现值随机变量的概率分布。对所有的利率轨道，重复计算后得到一系列现值随机变量的概率分布。这些现值随机变量的概率分布的特点是：概率值相同，而随机变量的取值可以不同。

下一步对得到的一系列概率分布进行混合（mixture），最终得到一个完整的概率分布。其原理可以表述如下：

如果设X|Λ=λ的条件概率分布函数为Fλ(x)，而且Λ的概率密度函数为u(λ)，那么X的无条件概率分布函数F(x)可以表示为：

假设模拟的利率轨道共有m条，在对上式进行数值计算时，首先要进行离散化处理。如果设每条利率轨道的权重均为，而且每个Fλ都由离散分布来代替，那么上述积分就可以变为离散求和。具体来说，在本文的两全保险数值实例中，每一条利率轨道对应一个λ*值，则有，而分布函数Fλ*的形式如式（5）所示。值得注意的是，因为本文数值实例部分采用离散时间的有限状态Markov链，所以不同的利率轨道的现值随机变量的取值是有可能重复的。

在得到分布函数F(x)之后，即可获得概率密度函数f(x)，进而可以计算随机变量X的均值、方差、进一步也可得到各个分位数以及相关的分布度量等。

2.4 定期生命年金

由中国人寿保险业经验生命表（2000～2003），得到各个年度内的死亡概率，以及生命年金期限届满时的生存概率。针对每条利率轨道，便可得到一个现值随机变量的概率分布[5]。对所有的利率轨道，重复计算后便可得到一系列现值随机变量的概率分布。最后通过对一系列概率分布进行混合后得到一个完整的概率分布。其计算原理与两全保险的情形相同，这里不再赘述。

3 数值实例及结果分析

假设Markov链利率取值为{0.02,0.03,0.04}，转移概率矩阵为：

下面通过数值实例进行分析，按照本文第2节的思路，基于Markov链随机利率对三个有代表性的精算函数进行了随机模拟，得到相应精算函数变量的概率分布及相关的分布特征。这里采用R软件对其进行算法实现。

3.1 在Markov链随机利率下，模拟确定性年金现值与积累值随机变量的概率分布

为了简便，假设年金的支付额都为1个单位。如果利率i为常数，那么这里考虑的期末付年金现值就是通常意义下的其计算公式为：

对于20年定期年金来说，根据初始利率的不同假设，下面分两种情况进行讨论。

（1）假设初始利率i1=0.03，每条利率轨道包含20个利率。通过模拟1000条利率轨道，分别得到包括1000个①为了统一，这里采用与本文后续两个实例相同的处理方式，在年金现值和积累值各自的1000次模拟结果中，不同结果的个数恰好都为1000个。定期年金现值与积累值的样本数据，进而可以描述该年金现值与积累值随机变量的概率分布。图1给出了初始利率为常数0.03情况下，基于Markov链随机利率模拟得到的年金现值与积累值随机变量的概率分布，其相应的分布特征如表1所示。

图1 初始利率为常数0.03时Markov链模拟的年金现值与积累值随机变量的概率分布

表1 初始利率为常数0.03时年金现值、积累值随机变量的概率分布的分布特征

为了便于比较，表2给出了利率分别为常数0.02、0.03、0.04情形下，计算的确定性年金现值和积累值。注意到表1中模拟得到的年金现值随机变量的均值为14.93，积累值随机变量的均值为27.62，它们明显处于这三个常数利率计算的现值、积累值的范围内。这表明在有限的利率波动范围内，采用确定性利率计算年金现值和积累值是合理的。另外，这两个均值都与利率为0.03时计算的年金现值和积累值非常接近，这与Markov链模拟时初始利率的选择有关。

若模拟次数为10000次，则得到的年金现值随机变量的均值为14.93，标准差为0.33；积累值随机变量的均值为27.57，标准差为0.77。与表1中的均值和标准差相比，差异不大，从一定意义上说，采用1000条利率轨道模拟的结果已相当可靠。

表2 三个常数利率下计算的20年定期年金现值和积累值

（2）假设初始利率的概率分布为Markov链的平稳分布，每条利率轨道包含20个利率。经计算，该Markov链的平稳分布为 (π1, π2, π3)=(0.2308,0.5769,0.1923)。通过模拟10000条利率轨道，分别得到包括10000个②与注释①类似，在年金现值和积累值各自的10000次模拟结果中，不同结果的个数都为9991个。定期年金现值与积累值的样本数据，进而可以描述该年金现值与积累值随机变量的概率分布。图2给出了初始利率为平稳分布情况下，基于Markov链随机利率模拟得到的年金现值与积累值随机变量的概率分布，其相应的分布特征如表3所示。

表3 初始利率为平稳分布时年金现值、积累值随机变量的概率分布的分布特征

注意到表3中模拟得到的年金现值随机变量的均值为14.94，积累值随机变量的均值为27.56，它们也明显处于这三个常数利率计算的现值、积累值的范围内，且与利率为0.03时计算的年金现值和积累值最为接近。另外，若模拟次数为1000次，则得到的年金现值随机变量的均值为14.93，标准差为0.37；积累值随机变量的均值为27.59，标准差为0.77。与表3中的均值和标准差相比，差异也不大，从一定意义上说，这里采用1000条利率轨道进行模拟也是相当可靠的。

3.2 在Markov链随机利率下，模拟两全保险现值随机变量的概率分布

考虑关于(30)岁男性的20年期两全保险，保额为1000元。引用中国人寿保险业经验生命表非养老金业务表（2000～2003,CL1）。对应于不同的利率假设，其趸缴净保费及相应的标准差的计算结果如表4所示。

表4 不同的利率假设下两全保险现值随机变量的均值和标准差

现模拟计算在Markov链随机利率下，该两全保险现值随机变量的概率分布。类似于前面的讨论，根据初始利率的不同假设，下面分两种情况进行讨论。

（1）假设初始利率i1=0.03，每条利率轨道包含20个利率。模拟1000条利率轨道，得到1000个概率分布，对这1000个概率分布混合后，得到该两全保险现值随机变量的概率分布，如图3中的左图所示，其相应的分布特征如表5第二列所示。值得注意的是，这里选定了随机模拟的“种子”数为set.seed(2525)③应用R软件进行随机模拟时，可以设定不同的“种子”数。选择同一个“种子”数，一方面可以唯一确定模拟结果，另一方面有助于对各种模拟方法的结果进行比较。。在这种情况下，该现值随机变量的不同取值共有902个。按照从小到大顺序排列这些现值，其中第500个现值为655.01，而取值大于655.01的概率仅为1.93%。

（2）假设初始利率的概率分布为Markov链的平稳分布，每条利率轨道包含20个利率。模拟10000条利率轨道，得到10000个概率分布，对这10000个概率分布混合后，得到该两全保险现值随机变量的概率分布，如图3中的右图所示，其相应的分布特征如表5第三列所示。在选定了随机模拟的“种子”数（set.seed(2525)）后，该现值随机变量的不同取值共有1240个。按照从小到大顺序排列这些现值，其中第700个现值为660.73，而取值大于660.73的概率仅为2.18%。的

图3 两种情况下两全保险现值随机变量的概率分布

表5 两种情况下两全保险现值随机变量的概率分布的分布特征

与表4相比，两种情况下模拟得到的两全保险现值随机变量的概率分布的均值都明显处于这三个常数利率计算的趸缴净保费的范围内，且与利率为0.03时的趸缴净保费最为接近。而两种情况下模拟得到的概率分布的标准差都高于这三个常数利率计算的标准差，这可以直观地理解为Markov链随机利率带来了更大的波动性。另外，当初始利率为常数0.03时，若模拟次数为10000次，则得到的两全保险现值随机变量的均值为563.29，标准差为41.67；当初始利率为平稳分布时，若模拟次数为1000次，则得到的两全保险现值随机变量的均值为562.21，标准差为42.43。从一定意义上说，采用1000条利率轨道对两种情况进行模拟都是可行的。

从图3可以看出，两种情况得到的两全保险现值随机变量的概率取值忽大忽小，其概率分布存在多个峰值，并非传统意义上的单峰分布。当然，从理论上讲，这种分布也是存在的，我们可以将其看作是一类特殊的奇异分布。

3.3 在Markov链随机利率下，模拟定期生命年金现值随机变量的概率分布

考虑关于(30)岁男性的20年期期初付生命年金，为了简便，假设每次支付额都为1个单位。引用中国人寿保险业经验生命表非养老金业务表（2000～2003,CL1）。如果利率i为常数，那么这里考虑的精算现值就是通常意义下的其计算公式为：

其中，式（11）中的k表示整数生存年数。

对应于不同的利率假设，该定期生命年金的精算现值及相应的标准差的计算结果如表6所示。

表6 不同利率假设下定期生命年金现值随机变量的均值和标准差

现模拟计算该定期生命年金现值随机变量的概率分布。类似于前面的讨论，根据初始利率的不同假设，下面分两种情况进行讨论。

（1）假设初始利率i1=0.03，每条利率轨道包含20个利率。模拟1000条利率轨道，得到1000个概率分布，对这1000个概率分布混合后，得到该定期生命年金现值随机变量的概率分布，如图4中的左图所示，其相应的分布特征如表7第二列所示。在选定了随机模拟的“种子”数（set.seed(2525)）后，该现值随机变量的不同取值共有10846个。按照从小到大顺序排列这些现值，其中第5000个现值为12.10，而取值小于12.10的概率仅为1.96%。实际上，经过细致的研究发现，在第5000个现值以后有1000个不同的现值对应的概率都相等，加总后的概率为0.9679，对应的现值随机变量的取值范围为[14.56,16.38]。为了更清晰地展示这一特征，图4中的右图给出了这1000个不同现值对应的概率值。

图4 初始利率为常数0.03时Markov链模拟的定期生命年金现值随机变量的概率分布

（2）假设初始利率的概率分布为Markov链的平稳分布，每条利率轨道包含20个利率。模拟1000条利率轨道，得到1000个概率分布，对这1000个概率分布混合后，得到该定期生命年金现值随机变量的概率分布，如图5中的左图所示，其相应的分布特征如表7第三列所示。在选定了随机模拟的“种子”数（set.seed(2525)）后，该现值的不同取值共有11849个。按照从小到大顺序排列这些现值，其中第6000个现值为12.07，而取值小于12.07的概率仅为1.97%。实际上，经过细致的研究发现，在第6000个现值以后也有1000个不同的现值对应的概率都相等，加总后的概率为0.9679，对应的现值随机变量的取值范围为[14.42,16.49]。为了更清晰地展示这一特征，图5中的右图给出了这1000个不同现值对应的概率值。

图5 初始利率为平稳分布时Markov链模拟的定期生命年金现值随机变量的概率分布

表7 两种情况下定期生命年金现值随机变量的概率分布的分布特征

与表6相比，两种情况下模拟得到的定期生命年金现值随机变量的概率分布的均值、标准差都明显处于这三个常数利率计算的精算现值、相应的标准差的范围内，且与利率为0.03时的精算现值、相应的标准差最为接近。另外，当初始利率为常数0.03时，若模拟次数为100次，则得到的定期生命年金现值随机变量的均值为15.18，标准差为1.23；当初始利率为平稳分布时，若模拟次数为100次，则得到的定期生命年金现值随机变量的均值为15.17，标准差为1.23。从一定意义上，采用100条④这里模拟100次之所以可行，是因为100次模拟可以得到100个概率分布，每个概率分布包括20个现值随机变量的取值，共可得到2000个现值随机变量的取值，且此时不同结果的个数已足够多了。利率轨道对两种情况进行模拟也是可行的。

从以上结论，初步猜测，不论初始利率如何选择，若利率轨道数为m，则该定期生命年金现值随机变量的取值都将以很大的概率均匀分布在m个较大的取值上。进一步注意到，在该定期生命年金的数值实例中，这个很大的概率正是npx=20p30=0.9679。正如图4、5所示，这两种情况模拟得到的定期生命年金现值随机变量的概率分布也不是传统意义上的单峰分布。从理论上讲，这种分布也是存在的，可视为一类特殊的奇异分布。

4 研究结论与方法推广

4.1 研究结论

（1）Markov链随机利率模型具有易于理解、符合直观、编程方便、能自动处理变量的不独立性问题、模拟的利率不会取负值或较大值等优点。在R软件下的编程实现运算时间少、效率较高。本文在R软件算法实现中，采用了同一个“种子数”，这样有助于进一步对两种利率轨道的结果进行比较。

（2）两种利率轨道模拟结果的一致性。从图1和图2、图3、图4和图5可以看出，两种利率假设下，基于Markov链随机利率模拟的三个精算函数随机变量相应的概率分布的图形都很相似；同时，从表1和表3、表5、表7也可以看出，其相应的分布特征也都很接近。

（3）实际上，对Markov链来说，如果存在极限（平稳）分布，那么该分布不依赖于初始状态。另外，一般来说，经过较少的转移次数，Markov链即可收敛于极限（平稳）分布。很容易验证，在本文选择的Markov链的转移概率下，不论初始状态如何，Markov链在不同时刻的概率分布一般都较快地收敛到极限（平稳）分布。因此，本文三个数值实例中，初始利率为常数和初始利率的概率分布为平稳分布两种情况得到的结论差异不大。

（4）从表1、表3可以看出，在实例1中，两种情况下模拟得到的确定性年金现值、积累值随机变量的概率分布的均值和中位数几乎相同；从表5、表7可以看出，在实例2和实例3中，两种情况下模拟得到的现值随机变量的概率分布的均值和中位数有一定差异。这可以从确定性年金、两全保险、定期生命年金本身的差异得到解释。其中，确定性年金与死亡概率无关，因而每次模拟都可以得到一个确定的年金支付额，多次模拟的结果一般呈现出接近正态分布的单峰分布；而两全保险、定期生命年金的支付额是不确定性的，受死亡概率的影响，每次模拟只能得到一个具体的概率分布，每个概率分布又包括若干个现值随机变量，其波动性很大，多次模拟的结果一般都不能呈现出正常意义下的分布，我们将其称为是一类特殊的奇异分布。

4.2 方法推广

（1）本文在Markov链随机利率模型假设下，应用随机模拟方法，对精算函数随机变量的概率分布进行了数值求解。具体地讲，本文采用了两种产生利率轨道的思路；一种是基于初始利率为常数的假设；另一种是基于初始利率的概率分布为Markov链的平稳分布的假设。并通过寿险精算学中的三个数值实例，分别给出了这两种假设下模拟得到的精算函数随机变量的概率分布以及相关的分布特征。解决这些例子的思路与方法可以推广到精算学和金融学的进一步研究中。这是一个有待深入研究的新方向。

（2）寿险利率市场化对我国精算技术提出了很大的挑战，将更符合市场变化规律的随机利率模型应用到寿险公司的定价、准备金评估以及随机资产负债管理中，具有重要的现实意义。本文对随机利率下寿险精算中一些基本的随机变量的概率分布进行了研究，作为后续研究，将进一步对Vasicek随机利率模型、Cox-Ingersoll-Ross（CIR）随机利率模型、Multiplicative Shock随机利率模型等进行系统研究，并应用R软件进行算法实现，以期对各种寿险产品的保费、准备金、现金价值、部分利润指标等的变动进行分析，并进一步比较随机利率下保险合同的价值与固定利率下保险合同的价值等。

（3）不断更新的计算技术和日益复杂的金融市场风险，以及对金融数学认识的深入和学科交叉，促进了保险公司利率市场化的进程。进一步的研究内容是：在对随机利率下寿险精算理论进行深入研究的基础上，针对我国寿险业的利率市场化问题，提出可行性建议，为我国寿险利率市场化的进程提供理论支持和技术保证。

[1] Goovaerts,M.,Dhaene,J.,De Schepper,A.Stochastic Upper Bounds for Present Value Functions[J].Journal of Risk and Insurance,2000，(67).

[2] Grinstead,C.M.,Snell,J.L.Introduction to Probability(2ndEdition)[M].American Mathematical Society,1997.

[3] Bühlmann,H.Stochastic Discounting[J].Insurance:Mathematics and Economics,1992，(11).

[4] Beekman,J.A.,Fuelling,C.E.Interest and Mortality Randomness in Some Annuities[J].Insurance:Mathematics and Economics,1990，(9).

[5] De Schepper,A.,De Vylder,E.,Goovaerts,M.,Kaas,R.Interest Ran⁃domness in Annuities Certain[J].Insurance:Mathematics and Econom⁃ics,1992,11(3).

[6] Vanneste,M.,M.Goovaerts,E.Labie.The Distribution of Annuities[J].Insurance:Mathematics and Economics,1994,15(1).