アクチュアリーのためのPython入門（予定死亡率の補整１）

アクチュアリーのためのPython入門（数学編第1回）

予定死亡率の補整１

📚 アクチュアリーのためのPython入門
この記事はアクチュアリー数学編の一部です。
▶ 目次はこちら
▶ 逆引きガイドはこちら

はじめに

今回は、これまでの内容とは異なり、
生命保険数学というよりは確率・統計に近い部分を再現します。

生命保険会社では、自社の過去データを使って、
予定死亡率を作成することがあります。

自社の過去データから作成した粗死亡率は、
そのままでは、凸凹があったり、
高齢のデータがなかったりということで、
平滑化や補外を行います。

また、過去データの平均をそのまま使うのではなく、
ある程度の安全を見込んで、死亡率を割増もします。

予定死亡率の一般的なものとして、
標準生命表がありますので、
今回はその補整をPythonで再現します。

標準生命表の補整

まず、直近の標準生命表2018の作成過程は、
アクチュアリー会に資料がありますので、
それを参照ください。

その作成過程を確認すると、死亡率の補整は3回行っています。

• 1次補整
  正規近似により安全割増を付与
• 2次補整
  Greville補整による平滑化
• 3次補整
  Gompertz-Makehamによる高齢外挿

アクチュアリー会の資料と同じ数値を使って、
途中の計算結果をprintで確認しながら進めていきます。

1次補整の正規近似

まずは、前回まで使用している四捨五入のroundhuを使用。
指数や対数を使うのでmathライブラリをインポート。
mathは標準ライブラリなのでインストールは必要ありません。
途中の計算内容を分かりやすくするためにDataFrameを使用したいので、
pandasもインポートしておきます。

また、q0は粗死亡率、Lxは標本分散を計算するための標本数です。

mortality.py

from core import Decimal, roundhu
import math
import pandas as pd

# 補整前死亡率
q0_float = [
    0.00076,0.00032,0.00022,0.00015,0.00011,0.00009,0.00009,0.00008,0.00007,0.00007,
    0.00007,0.00009,0.00009,0.00011,0.0001,0.00018,0.0002,0.00027,0.00046,0.00036,
    0.0005,0.00038,0.00055,0.00058,0.00051,0.00053,0.00046,0.00052,0.00046,0.00047,
    0.00057,0.00054,0.0006,0.00049,0.00053,0.0006,0.00066,0.00073,0.00071,0.0008,
    0.00096,0.00103,0.00103,0.00117,0.00127,0.00134,0.00148,0.00166,0.00175,0.00203,
    0.00234,0.00241,0.00262,0.00298,0.00319,0.00344,0.00355,0.00413,0.0044,0.00506,
    0.00535,0.00592,0.00641,0.00705,0.00796,0.00857,0.00923,0.00986,0.01052,0.01205,
    0.01284,0.01417,0.01574,0.01765,0.01943,0.02234,0.0249,0.02815,0.03224,0.03678,
    0.04228,0.0477,0.05488,0.05996,0.06857,0.07548,0.08601,0.09882,0.11026,0.11164,
    0.11916,0.14499,0.15218,0.17747,0.15468,0.1832,0.16253,0.10308,0.11305,0.00802
]
# 補整用標本数
Lx_float = [
    515,609,718,844,988,1152,1338,1548,1785,2050,
    2346,2674,3037,3436,3873,4349,4865,5422,6020,6659,
    7338,8056,8811,9600,10421,11270,12142,13032,13935,14845,
    15754,16657,17545,18411,19247,20045,20799,21499,22140,22713,
    23215,23638,23979,24233,24397,24471,24452,24342,24142,23852,
    23478,23023,22492,21891,21226,20503,19731,18917,18068,17192,
    16297,15391,14481,13573,12674,11791,10927,10089,9281,8505,
    7764,7062,6399,5776,5194,4654,4154,3693,3272,2888,
    2539,2224,1941,1687,1461,1261,1084,928,792,673,
    570,481,404,338,282,234,194,160,131,108
]
# Decimal型に変換
q0 = [Decimal(str(x)) for x in q0_float]
Lx = [Decimal(str(x)) for x in Lx_float]

1次補整は安全割増です。
詳細はアクチュアリー会の資料にありますので省略します。
30%の増加を条件として、
2項分布の上側約$\varepsilon$点を正規近似して、
$$\overline{q}=q+u(\varepsilon)\sqrt{\frac{q(1-q)}{n}}$$
で計算します。

mortality.py

# 1次補整（安全割増）
# 1次補整死亡率
q1 = []
# 上側ε点
epsilon = Decimal("2")
# 上限(30%増が上限の意味)
lim = Decimal("0.3")
# 1次補整の実行
for q, L in zip(q0, Lx):
    sigma = Decimal(str(math.sqrt(q * (Decimal("1") - q) / L)))
    s = roundhu(q + min(epsilon * sigma, lim * q), 5)
    q1.append(s)

ここまでの途中結果q1を表示すると、
Decimal型なので、

[Decimal('0.00099'), Decimal('0.00042'), Decimal('0.00029'), …

となっています。

正規分布の関数

なお、今回は2$\sigma$と指定されていますが、
2項分布の上側約$\varepsilon$点の$\varepsilon$での指定だった場合、
累積正規分布の逆関数norm.ppf(a)を使うことができます。

次の例の場合は上側2.275%点を計算しています。
Decimal(str(u))が約2として計算しています。

なお、scipy.statsは確率分布を使うためのライブラリで、
scipyのインストールが必要です。

from scipy.stats import norm
# 上側ε点
epsilon = 0.02275
# 1次補整の実行
for q, L in zip(q0_float, Lx_float):
for q, L in zip(q0, Lx):
    sigma = Decimal(str(math.sqrt(q * (1 - q) / L)))
    u = norm.ppf(1 - epsilon)
    s = q + Decimal(str(u)) * sigma
    q1.append(roundhu(s, 5))

scipy.statsは、正規分布norm以外にも
二項分布binomやポアソン分布poissonなど、
アクチュアリーの確率統計でおなじみの確率分布があります。

2次補整の平滑化

続いて2次補整の平滑化です。
Greville補整は係数が決まっていて、
それを前後に乗じて滑らかにするわけですが、
0歳を補整するには、その下の仮数値を作ってあげる必要があります。
0～5歳を使って-1歳を作成、
-1～4歳を使って-2歳を作成と$\dots$、-6歳まで作成します。
insertはリストに挿入の意味です。
ここではリストの先頭に挿入します。

mortality.py

# 2次補整（平滑化）
# Greville係数
c_float = [0.240058, 0.214337, 0.147356, 0.065492, 0, -0.027864, -0.01935]
a_float = [1.016301, 0.360880, -0.021625, -0.160909, -0.138330, -0.056317]
ci = [Decimal(str(x)) for x in c_float]
aj = [Decimal(str(x)) for x in a_float]
# 0歳未満の外挿
for x in range(0, 6):
    q = 0
    for j in range(0, 6):
        q += q1[j] * aj[j]
    q = roundhu(q, 5)
    q1.insert(0, q) # q1の先頭に挿入

-6歳まで作成しますが、このままリストを使っていると、
何歳を指定したのか分からなくなります。
そこで一旦DataFrameを使い、インデックスを年齢にします。

mortality.py

# 外挿用死亡率に年齢を結合
age = list(range(-6, 100))
df = pd.DataFrame({
    "age": age,
    "qx": q1
})
# インデックス番号を年齢に
df = df.set_index("age")

DataFrameの中身を確認すると、
このようになっています。

          qx
age         
-6   0.00243
-5   0.00217
-4   0.00191
-3   0.00165
-2   0.00137
..       ...
 95  0.23378
 96  0.21129
 97  0.13400
 98  0.14697
 99  0.01043

2次補整は定数を乗じて、計算していきます。
DataFrame化したのは、次の式を分かりやすくするためです。
ただ、DataFrameのままでは、
単一のデータを取得するには式が長くなってしまうので、
計算結果はリストに戻しています。

mortality.py

# 2次補整死亡率
q2 = []
# Grevilleの平滑化
for t in range(0, 94):
    q = (df.loc[t, "qx"] * ci[0]
    + (df.loc[t - 1, "qx"] + df.loc[t + 1, "qx"]) * ci[1]
    + (df.loc[t - 2, "qx"] + df.loc[t + 2, "qx"]) * ci[2]
    + (df.loc[t - 3, "qx"] + df.loc[t + 3, "qx"]) * ci[3]
    + (df.loc[t - 4, "qx"] + df.loc[t + 4, "qx"]) * ci[4]
    + (df.loc[t - 5, "qx"] + df.loc[t + 5, "qx"]) * ci[5]
    + (df.loc[t - 6, "qx"] + df.loc[t + 6, "qx"]) * ci[6])
    q = roundhu(q, 5)
    q2.append(q)

q2が途中結果です。
これを確認すると、

[Decimal('0.00081'), Decimal('0.00056'), Decimal('0.00036'), …

となっています。

これで2次補整が終わりました。
ここまでは問題なく再現できたのですが、
3次補整では、少し問題が発生していますので、
説明のため次回とさせてください。

まとめ

今回は、標準生命表の2次補整までをPythonで再現しました。
同じことをすれば、自社データの粗死亡率も平滑化までできます。

ポイントは、0歳未満の外挿を分かりやすくするため、
DataFrame化したところです。
自分が分かるだけでいいなら、ここまでしなくていいのですが、
分かりやすさを重視しました。

また、計算では使っていませんが、
正規分布のnormも保険数学ではいろいろ応用できますので、
いろいろと試算してみてください。

---

📚 ナビゲーション

▶ 次の記事
予定死亡率の補整２

📚 目次
アクチュアリーのためのPython入門