CT4 Survival Model 精算模型学习笔记

Jackie

根据赫瑞瓦特的 CT4 精算模型课程对全书的知识点做的梳理，和协会的CT4考试内容不完全一致，适合预习或者复习CT4的时候理一遍思路。

把一学期老师上课内容梳理一遍+自己理解～很基础但很有用～

如果有问题请指出～希望能有所收获～

本文架构：

生存模型引入
估计生存时间分布
马尔可夫模型基础
马尔可夫模型数据及估计
死亡率的二项和柏松分布
修匀（Graduation）方法及统计检验
风险暴露

生存模型引入

讨论生存模型时很多都是在医学统计上常用的概念。

这个模型，运用很多不同的方法，想得出的就是生存函数和死亡率的估计来解决实际问题。

这篇文章需要部分人寿保险数学的基础，但是没有也没关系，简单的～

首先要知道这几个概念：

$$
F(t)=P(T\leq t)=\qx{t}{x}
$$

其中，$T$ 指存活时间，$F(t)$ 就是 $T$ 的分布函数，即在 $t$ 时刻前死亡的概率分布。

$$
S(t)=P(T\geq t)=1-F(t)=\px{t}{x}
$$

$S(t)$ 为生存函数，年龄为 $x$ 岁的人，活过 $t$ 时刻的概率。

$\mu_x$ 为 $x$ 时刻的死亡力（一瞬间的事儿）。

其中还涉及到一些变换技巧，刷刷题练练手感就好。

估计生存时间分布（K-M方法，Cox模型）

1.Cohort

罗马古代军团的意思，在这里代表拿出一群人，从出生到死亡一一记录，形成生命表。

但是，这个太费时太费力，所以，我们一般用 mini-cohort，也就是比如一群在2017年30岁的人，在2017年到2019年，死亡和存活的人数，然后计算，这样比较方便简洁，但是误差肯定比前面的大。

为了估计死亡或者生病的时间，引入 censor 的概念：right-censoring, left-censoring 和 interval-censoring，一般来说用第一个。

举个 right-censoring 的例子：a 时刻，观察一个人，好，他活着，说明他在 a 时间点之后死亡。

对于 left-censoring：同样 a 时刻，那个人生病了，说明 a 点之前某点是健康的状态。

对于interval-ensoring：假设 a,b 两点，a 点在 b 点前，a 时刻健康，b 时刻死亡，说明 a,b 中某点是生病的状态。

还有 random censoring, non-information censoring, type 1, type 2 censoring。

详细的 censoring 的知识点总结可以看 Jackie 的推文：CENSORED DATA：当注射格列宁的小白鼠逃出实验室

2. 用 Kaplan-Meier 估计方法估计生存函数。

首先引入概念：

$$
\px{t}{x}=S(t)=\px{}{x} \times \px{}{x+1} \times \px{}{x+2} \times \px{}{x+3} \times \ldots \times \px{}{x+t-1}
$$

然后，我们在每个时间点 $t$ 都有：

$n_t$：刚好在 $t$ 点之前（$t^-$）所有存活的人数
$dt$：在 $t$ 点的死亡人数

所以说 $t$ 点的存活率为：
$$
\frac{(n_t-dt)}{n_t}
$$

根据最开始的概念得到：

$$
S(t)=\prod_{j=1}^{k}(n j-d j) / nj
$$

然后根据我们的数据，就可以把生存函数算出来啦～

3.Cox Regression Model（Cox回归模型）

但是呢，KM方法有一个大缺点，就是在很多协变量（covariate）的情况下很麻烦。协变量不是x,y但是影响实验的结果，例如我们现在找时间和死亡率的关系，性别就是一个协变量。如果协变量少，我们就可以分组来用 KM，但是，现实是残酷的，有很多协变量（是否吸烟等），所以引入了Cox Regression Model（Cox回归模型）。

定义了 Cox proportional hazards model，因为其是一个半参数模型，已知 baseline hazard，运用 partial likelihood function（部分似然函数）算出 $\beta$ 期望和方差。

关于为什么用 partial likelihood function，老师上课没说，我搜到的解释是因为在 Cox proportional hazards model 里面 baseline hazard 是早就定义好的，所以是一个半参数模型，用 partial likelihood function。

之后对 $\beta$ 做假设检验，$H_0$为 $\beta=0$ 可以用 Wald Test, Likelihood-Ratio Test( $-2\times (l_0-l_1)$ ), Score test( $\frac{U\left(\beta_{0}\right)^{2}}{I\left(\beta_{0}\right)}$)（检验所有回归的经典算法）去算，差别很小，都用 likelihood，就是估计模型的多少的差别，都是套路。

马尔科夫模型基础

对于学过随机变量的筒子应该对这个不陌生，但没有学过也没有关系，其实很简单易懂。

小白我的理解：

马尔可夫模型其实是多个状态的转换。

可以想象你在一个地方呆了一会儿，然后以一个超快的速度移动到另一个点，之后再移动到另外一个状态再移动......（想想很早很早很火的表情包，哎呀找不到了哭）直到到了一个吸收态出不去了（比如死亡），那就一直待到那里了。

我们把移动的速率叫做 transition intensity。

在生存模型中，我们最初可以考虑只有两种状态的情况，状态1是活，2是死，我们移动的速度就是我们的死亡力。之后考虑1是健康，2是生病，3是死亡等等的推广。

为了理解和计算，对于马尔可夫模型，我们给出三个假设：

第一：是马尔可夫模型的定义就是不考虑之前的状态，对于未来，只管现在的年龄和状态。

比如A之前生病很长时间，然后蹦到健康的状态，之后蹦到什么状态只跟现在有关，虽然有点不切实际，但是如果在开始就考虑所有的话会混乱。

第二：在短时间内，我们只能转移一次，转移的几率很小（o(dt)）。

比如给定很短很短的时间，我们只能从a跳到b，而不能中途再来一个a到c到d到b这种情况。其实是符合实际情况的，很短的时间我们不能很快健康到生病再到健康。

第三：短时间内，两个不同状态a,b转换的概率为：

$$
d t \actsymb{}{}{p}{ab}{x+t}=\mu_{x+t}^{a b} d t+o(d t)
$$

也就是在二的基础上，说直接两状态的 transition intensity 乘以很短的时间，再加一个转移多次的 $o(dt)$.

对于多状态马尔可夫模型，有了transition intensity，就可以根据假设和公式把转移进入某状态和转移出某状态的概率求出来了。

有一个很直观和简单的方法。

比如你在 $x$ 岁时是状态1，如果要求从1到2的概率。

首先列出 $x+t$ 岁，状态1，2，3......，然后，在时间点 $x+t+dt$，这些状态到状态2的 transition intensity 乘以 $dt$，最后转换一下，求出从1到2的，关于时间 $t$ 的倒数，概率之后自会有方法解决～

马尔可夫模型数据及估计

啦啦啦，终于进入实战了，首先是最简单的两个状态（生或死）模型。

首先来几个概念，之前介绍了censor，也就是保险公司（或者其他什么的）的调查时间点（不可能一直调查，考虑实际情况）。

现在呢，对于某人 $i$，我们把开始调查的时间记作 $a$，结束调查的时间（没死）为 $b$，调查终止的时间（死或 censor）记为 $t$（可能为 $b$）。

然后，我们的等待（调查）时间就是 $t-a$，在等待时间内 $D=0$ 表示活着，$D=1$ 表示死。

接下来呢，推导两种情况：
$$
\begin{aligned}
&P\left(D_{i}=0, V_{i}=b_{i}-a_{i}\right)\\
&=P\left(a_{i} 到 b_{i} 存活\right)\\
&=\exp \left(u \times v_{i}\right)(因为恒定 u)\\
&=e^{-\mu v_{i}} \times \mu^{d i}
\end{aligned}
$$

$$
\begin{aligned}
&P\left(D_{i}=1, v_{i}<V_{i}<v_{i}+d_{i}\right)\\
&=v_{i} P_{a_{i}+x} \times\left(\mu d v_{i}+o\left(d v_{i}\right)\right)\\
&=e^{-\mu v_{i}} \times \mu+o\left(d v_{i}\right)
\end{aligned}
$$

无论在D=0或1的情况下，等待时间的概率都可以用一个函数（包含d与v的混合分布）代表（假定在 $t$ 属于$[0,1]$内，$\mu_{x+t}=\mu_x$ 恒定）:

$$
f_{i}\left(d_{i}, v_{i}\right)=e^{-\mu v_{i}} \times \mu^{d i}
$$

然后，把所有人乘起来，就是极大似然函数，最后求得 MLE，得到：

死亡力=所有死亡数量/等待时间的加和。

$$
\tilde{\mu}=D / V
$$

我的直观理解：

假定一天内，恒定死亡力，10个人，1个逝去，那么死亡力就是1/10。

然后，再常规的 score function, information function 用极大似然函数求到，最后算出方差标准差。

类似的，可以推广到多状态模型，先算出一个的概率，之后到极大似然，不同的transition intensity 都可以用MLE算出。

死亡率的二项和柏松分布

因为简化下来，死亡人数就应该是个二项分布（生or死），所以死亡人数$D\sim B(n,\qx{}{x})$，然后我们已知二项的 $f(x)$，近似在（0，1）死亡率不变，算个 MLE 就可以得到$\bar{\qx{}{x}}=d/n$ .

但是 $\qx{}{x}$ 是会随着时间变化而改变的，所以我们现在就需要引入假设。

我们引入三种假设（都是在 $t \in [0,1] $ 的近似假设， $\qx{t}{x}$ 指现在 $x$ 岁，在 $x$ 岁和$x+1$ 岁之间 $t$ 刻死亡的概率）：

Uniform distribution of deaths (UDD)

$$
\qx{t}{x}=t\cdot \qx{}{x}
$$

所谓 uniform 就是指 $\qx{}{x}$ 在这个区间不变。

Constant force of mortality (CFM)

$$
\mu_{x+t}=\mu_x
$$

不用说。

Balducci assumption (BA)

就是跟 UDD 假设相反，UDD 是从 $x$ 开始看，BA 是从 $x+1$ 开始看。所以：

$$
\qx{1-t}{x+t}=(1-t)\cdot \qx{}{x}
$$

之后就进入实战了，像上一章所说：

$$
\qx{}{i}=(b_i-a_i)\qx{}{x+a_i}
$$

然后运用上面随意一个假设求出近似值然后运用 MLE 就好了～

Poisson Model 就是因为上一章我们 markov model 写出来的 maximum likelihood 可以用 poisson 的形式表示，所以认为 $D\sim Poisson(\mu E_x^c)$,其中 $E_x^c$ 就是之前的 $V$（observed waiting time）。

然后我们就可以根据不同情况运用这三种模型了。

Method of Graduation and Statistical Tests

Gompertz law

根据 Gompertz law: $\mu_x=BC^x$ （$x$ 为年龄，B，C为常数），那么 $ln(\mu_x)$ 就应该跟 $x$ 线性相关。

接下来，就是linear regression的战场了。

对于 normal error 的话，用 simple 就好。

但涉及到poisson和bionomial（五讲的知识哦），就得用 generalised linear model 各自的 link。

$$
D_{x} \sim B\left(n, q_{x}\right)
$$

$$
\bar{\qx{}{x}}=\frac{D_{x}}{E_{x}}
$$

$$
D_{x} \sim Poisson \left(\mu_{x}E_{x}^{c}\right)
$$
$$
\log \left[E\left(D_{x}\right)\right]=\log \left(E_{x}^{c}\right)+\log \left(\mu_{x}\right)
$$

然后计算residual。

还可以用 standard table，运用的 weighted least squares和generalised linear model 来 graduate，把数据转换一下就好。

例如运用 possion，直接把后面 $\log \left(\mu_{x}\right)=b+a \times \log \left(\mu_{x}^{t}\right)$（之前 $\log \left(\mu_{x}\right)=a+b \times x(a g e)$）。

但是运用 standard table 有个问题，他没有算进一个人持有多份保单的情况，数据上报的时候保险公司也只会提供数量而不会把人头也报上去。

所以这个 duplicate 的问题，因为一人持有多份保单使得方差变大，大概估算每个年龄段 dupplicate 的比例，成比例缩小即可。

我们就可以开始 statistical tests 他是否合适了。

熟悉的 Chi Square test、standarlised deviation test，到对 residual 进行的 sign test（因为residual 太小，就对符号进行检验）还有 change of sign test、grouping of signs test，从+-符号数量，改变频次，有多少 group 来确定是否合适；最后再来一个 serial correlation test来检验 residual 是否 positive correlated，运用 R 就很简单了。

风险暴露

计算 $\tilde{\mu_{x}}=\frac{d_{x}}{E_{x}}$, $\tilde{\qx{}{x}}=\frac{d_{x}}{E_{x}^{c}}$，我们需要知道：第一是 central exposed to risk $E_{x}$，第二是initial exposed to risk $E_{x}^{c}$, 关系是：$E_{x}^{c} \approx E_{x}+0.5 \cdot d_{x}$.

然后我们探讨 principle of correspondence，也就是 exposure data death data 必须是同一时间界定的，不然 $dx$ 我界定为 nearest birthday，而 exposure 我按 calendar year censor，两个算起来就不对。

那么问题来了，年内死亡的时间该归为多少岁？

所以，为了确定时间界定，接下来引入age label。

分为 age last birthday, age nearest birthday。

前一就是在62.8岁死亡，归为62岁，后一就归为63岁。

然后看到 rate interval（没查到中文是啥），就是怎么区分 interval，分为 life year, calendar year 和 policy year，然后在每个 interval 里面，引入上面的时间界定。

也就是 rate interval 确定每个 interval 的年份，再由 age label 确定到底归为哪个年纪。

来一个例子，就如 rate interval: policy year; age label: age last birthday at the nearest policy anniversary。

一个人，birth 09/09/1910，died 10/05/2000， policy start 11/06/1998，死时算多少岁？

好我们看看，nearest policy anniversary 11/06/2000，此时的last birthday：09/09/1999，为89岁，然后我们就算为89岁。

之前都是知道进入和离开的具体时间，现在问题来了，如果没有记录时间呢？还有如果census和death data的age definition不同呢？

现在我们从census data里估算出 $E_{x}^{c}$。保险公司把数据提交的时候，一般是在一个特定的时间，我们就把 age label 定义为 age x last birthday。

然后，$t_1$,$t_2$ 之间的observation就是 $P_{x, t}, E_{x}^{c}=\int_{t 1}^{t 2} P_{x, t} d t$，我们可以用 trapezium rule 来近似，例如：

$$
\int_{0}^{2} P_{x, t} d t \approx 0.5 *\left(P_{x, 0}+P_{x, 1}\right)+0.5 *\left(P_{x, 1}+P_{x, 2}\right)=0.5 * P_{x, 0}+P_{x, 1}+0.5 * P_{x, 2}
$$

得到了central exposure，算出 $\mu_x$, $\qx{}{x}$.

之后涉及到 age label 的问题，我们就需要 adjust census data。

例如 $dx$ 是 age last birthday,而 $P_{x, t}^C$ 为 age nearest birthday测量的，我们就需要 adjust $P_{x, t}^{D} \approx 0.5 * P_{x, t}^{C}+0.5 * P_{x+1, t}^{C}$.

所以：

$$
E_{x}^{c}=\int_{t 1}^{t 2} P_{x, t}^{D} d t \approx 0.5 * \int_{t 1}^{t 2} P_{x, t}^{C}+P_{x+1, t}^{C} d t
$$

这个画图比较好理解，主要确定何时 $x$ 转变为 $x+1$ 。

就是这样啦～希望大家又学习到一点东西～有不对的请指出，谢谢啦～