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

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

予定死亡率の補整３

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

はじめに

前回までで、標準生命表をある程度は再現できたのですが、
今回は、今までの内容を少し修正して、
決算数理編で使用する予定死亡率を作成します。

具体的には、

• 男性と女性の予定死亡率を計算
• 30%上限の2$\sigma$ではなく、25%上限の1.5$\sigma$で計算

とします。Pythonの難易度は特に変化ありませんが、
Gompertz-Makehamのパラメータ推計の方法や初期値によって、
微妙に死亡率が変わってしまうので、
修正後のコードを載せておきます。

1次補整 安全割増

まずは、女性の粗死亡率と補整用標本数を加えます。
これらの数値は前回に引き続き標準生命表の作成過程からです。

mortality15.py

# 1.5σの死亡率作成
from core import Decimal, roundhu
import math
import pandas as pd

# 補整前死亡率（男性）
q0M_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
]
# 補整前死亡率（女性）
q0F_float = [
    0.00074,0.00029,0.0002,0.00013,0.0001,0.00008,0.00007,0.00007,0.00006,0.00005,
    0.00005,0.00005,0.00006,0.00007,0.00009,0.00011,0.00011,0.00018,0.00014,0.00018,
    0.00019,0.00022,0.0002,0.00021,0.00022,0.00023,0.00024,0.00023,0.00024,0.00027,
    0.00025,0.00031,0.00039,0.00033,0.00042,0.00047,0.00048,0.00054,0.00062,0.00064,
    0.00065,0.00075,0.00074,0.00082,0.00085,0.00092,0.00099,0.00118,0.00128,0.00143,
    0.00158,0.00161,0.00164,0.00189,0.00189,0.00224,0.0022,0.00223,0.00242,0.00249,
    0.0029,0.00309,0.00326,0.00337,0.00359,0.00385,0.00385,0.00434,0.00471,0.00528,
    0.00574,0.00628,0.00724,0.00822,0.00904,0.01072,0.01162,0.013,0.01482,0.0165,
    0.01981,0.02295,0.02626,0.02975,0.03487,0.03948,0.04511,0.052,0.05913,0.07086,
    0.07992,0.08549,0.09833,0.10293,0.1167,0.13895,0.12354,0.14921,0.09122,0.07899
]
# 補整用標本数（男性）
LxM_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
]
# 補整用標本数（女性）
LxF_float = [
    715,828,956,1100,1262,1443,1645,1869,2116,2389,
    2689,3016,3372,3759,4177,4626,5107,5620,6165,6742,
    7348,7984,8648,9336,10047,10778,11525,12285,13053,13825,
    14596,15361,16114,16851,17565,18250,18903,19516,20084,20604,
    21069,21477,21822,22103,22315,22458,22530,22530,22458,22315,
    22103,21822,21477,21069,20604,20084,19516,18903,18250,17565,
    16851,16114,15361,14596,13825,13053,12285,11525,10778,10047,
    9336,8648,7984,7348,6742,6165,5620,5107,4626,4177,
    3759,3372,3016,2689,2389,2116,1869,1645,1443,1262,
    1100,950,828,715,615,528,451,385,327,277
]
# Decimal型に変換
q0M = [Decimal(str(x)) for x in q0M_float]
LxM = [Decimal(str(x)) for x in LxM_float]
q0F = [Decimal(str(x)) for x in q0F_float]
LxF = [Decimal(str(x)) for x in LxF_float]

何度も計算しなおすわけではないので、
関数化までしていません。
sexを男性なら"M"、女性なら"F"を入力します。
死亡率、標本数、最終年齢は辞書化で性別をキーにしておきます。
なお、最終年齢は本来ならば3次補正後に残存表を作成して、
確認する必要がありますが省略します。

1次補整では安全割増を変更しています。

mortality15.py

# 死亡率の作成
sex = "F"
q0 = {"M":q0M, "F":q0F}
Lx = {"M":LxM, "F":LxF}
omega = {"M":110, "F":113}

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

2次補整 平滑化

前回までは男性の補整だったために抜いていた部分を
今回は女性も対応するため追加します。

男性はGompertz-Makehamのパラメータ推計の際に、
81歳～92歳の死亡率を使うので、
2次補整は98歳まで十分でした。
女性は81～94歳の死亡率を使うので、
2次補整を100歳までする必要があります。

粗死亡率は99歳までしかないので、
0歳未満で使った外挿を100歳以上にも使用します。

mortality15.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]
# 99歳超の外挿
for x in range(0, 6):
    q = 0
    for j in range(0, 6):
        q += q1[99 + x - j] * aj[j]
    q = roundhu(q, 5)
    q1.append(q)
# 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の先頭に挿入
# 外挿用死亡率も合わせて年齢を結合
age = list(range(-6, 106))
df = pd.DataFrame({
    "age": age,
    "qx": q1
})
# インデックス番号を年齢に
df = df.set_index("age")

外挿して-6歳から106歳までの死亡率ができたので、
99歳まで2次補整の平滑化を行います。

mortality15.py

# 2次補整死亡率
q2 = []
# Grevilleの平滑化
for t in range(0, 100):
    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)

3次補整 高齢外挿

まず、誤差を取る上限の年齢を変数にして、男女で別にします。
パラメータの初期値も男女で別にします。
あとは大きな変更がありません。
最初の方のsexを変更すれば、男性または女性の死亡率が計算できます。

mortality15.py

# 計算年齢の上限
max_age ={"M":92, "F":94}
# 3次補整（高齢外挿）
# 生存者数
lx = [1]
for t in range(0, 100):
    q = float(q2[t])    # しばらく端数処理をしないので、float型
    lx_next = lx[t] * (1 - q) 
    lx.append(lx_next)
# 死力をLagrangeの補間公式で近似
mux = []
for t in range(0, max_age[sex]+1):
    mux.append((8 * (lx[t - 1] - lx[t + 1]) - (lx[t - 2] - lx[t + 2])) / (12 * lx[t]))
# Gompertz-Makehamの最小2乗法の誤差
xmin = 81   # 年齢範囲下限
xmax = max_age[sex]   # 年齢範囲上限
def sse(params):
    parA, parB, parC = params
    err = 0
    for t in range(xmin, xmax + 1):
        err += (parA + parB * math.exp(parC * (t - xmin)) - mux[t]) ** 2
    # 精度を上げたいのでlogを取る
    return math.log(err)
   
# 誤差の最小化
from scipy.optimize import minimize

# パラメータの初期値と最小化
def_val = {"M":[-0.017, 0.07, 0.1], "F":[-0.01, 0.035, 0.1]}
result = minimize(sse, x0 = def_val[sex], method="Nelder-Mead", tol=1e-10)
parA, parB, parC, = result.x

# Gompertz-Makehamで84歳で接続
q3 = []
for t in range(0, omega[sex]):
    if t < 84:
        q = q2[t]
    else:
        q = 1 - math.exp(-(parA + parB / parC * (math.exp(parC) - 1) * math.exp(parC * (t - xmin))))
    q = roundhu(q, 5)
    q3.append(q)

df2 = pd.DataFrame({
    "age":range(0, omega[sex]),
    "q":q3
    })

このあと、実際はExcelファイルに出力して、
決算数理編で使用しています。

まとめ

今回は標準生命表の作成過程を再現して、
そこから、標準生命表より安全割増の少ない死亡率を作成しました。

実務では、引受基準を厳しくするかわり、
今回のように死亡率の安全割増を少なくしたり、
逆に引受基準を緩くする代わりに、
死亡率の安全割増を多めにするということが考えられます。

---

📚 ナビゲーション
◀ 前の記事
予定死亡率の補整２

▶ 次の記事

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