goalseek.py

from functools import wraps
from collections import defaultdict
import pandas as pd
import torch
from dataclasses import dataclass
from torch import optim

@dataclass(eq=False)
class Ctx:
    net_premium_pp: torch.Tensor


torch.set_default_dtype(torch.float64)

# constants
max_proj_len = 12 * 20 + 1

mp = pd.read_csv("BasicTerm_M/model_point_table.csv")
disc_rate = torch.tensor(pd.read_csv("BasicTerm_M/disc_rate_ann.csv")['zero_spot'].values)
mort_np = torch.tensor(pd.read_csv("BasicTerm_M/mort_table.csv").drop(columns=["Age"]).values)
sum_assured = torch.tensor(mp["sum_assured"].values)
issue_age = torch.tensor(mp["age_at_entry"].values)
policy_term = torch.tensor(mp["policy_term"].values)

# classes
class Cash:
    def __init__(self):
        self.reset()

    def reset(self):
        self.caches = defaultdict(dict)

    def __call__(self, func):
        @wraps(func)
        def wrapper(*args, **kwargs):
            key = (args, frozenset(kwargs.items()))
            if key not in self.caches[func.__name__]:
                self.caches[func.__name__][key] = func(*args, **kwargs)
            return self.caches[func.__name__][key]

        return wrapper

cash = Cash()

@cash
def get_annual_rate(duration: int):
    return mort_np[issue_age + duration - 18, min(duration, 5)]
@cash
def get_monthly_rate(duration: int):
    return 1 - (1 - get_annual_rate(duration)) ** (1/12)
@cash
def duration(t: int):
    return t // 12
@cash
def pols_death(t: int):
    return pols_if(t) * get_monthly_rate(duration(t))
@cash
def pols_if(t: int):
    if t == 0:
        return 1
    return pols_if(t - 1) - pols_lapse(t - 1) - pols_death(t - 1) - pols_maturity(t)

@cash
def lapse_rate(t: int):
    return max(0.1 - 0.02 * duration(t), 0.02)
@cash
def pols_lapse(t: int):
    return (pols_if(t) - pols_death(t)) * (1 - (1 - lapse_rate(t)) ** (1/12))
@cash
def pols_maturity(t: int):
    if t == 0:
        return 0
    return (t == 12 * policy_term) * (pols_if(t - 1) - pols_lapse(t - 1) - pols_death(t - 1))

@cash
def discount(t: int):
    return (1 + disc_rate[duration(t)]) ** (-t/12)
@cash
def claims(t: int):
    return pols_death(t) * sum_assured
@cash
def inflation_rate():
    return 0.01
@cash
def inflation_factor(t):
    return (1 + inflation_rate()) ** (t/12)
@cash
def expense_acq():
    return 300
@cash
def expense_maint():
    return 60
@cash
def pv_pols_if():
    return sum(pols_if(t) * discount(t)  for t in range(max_proj_len))
@cash
def pv_claims():
    return sum(claims(t) * discount(t) for t in range(max_proj_len))
@cash
def net_premium_pp():
    return pv_claims() / pv_pols_if()
@cash
def loading_prem():
    return 0.5
@cash
def expenses(t):
    return (t == 0) * expense_acq() * pols_if(t) \
           + pols_if(t) * expense_maint()/12 * inflation_factor(t)
@cash
def premium_pp(ctx: Ctx):
    return (1 + loading_prem()) * ctx.net_premium_pp
@cash
def premiums(t, ctx):
    return premium_pp(ctx) * pols_if(t)
@cash
def pv_premiums(ctx):
    return sum(premiums(t, ctx) * discount(t) for t in range(max_proj_len))
@cash
def pv_expenses():
    return sum(expenses(t) * discount(t) for t in range(max_proj_len))
@cash
def commissions(t, ctx):
        return (duration(t) == 0) * premiums(t, ctx)
@cash
def pv_commissions(ctx):
    return sum(commissions(t, ctx) * discount(t) for t in range(max_proj_len))
@cash
def net_cf(t):
    return premiums(t) - claims(t) - expenses(t) - commissions(t)
@cash
def pv_net_cf(ctx):
    return pv_premiums(ctx) - pv_claims() - pv_expenses() - pv_commissions(ctx)

ctx = Ctx(net_premium_pp=torch.zeros(len(mp), requires_grad=True))

optimizer = optim.LBFGS([ctx.net_premium_pp], lr=1)
def closure():
    cash.caches.clear()
    optimizer.zero_grad()
    loss = torch.sum(pv_net_cf(ctx) ** 2)
    print(loss)
    loss.backward()
    return loss
optimizer.step(closure)