Universal Life Policy Account Mechanics as a Differential Equation

JuliaActuary is an ecosystem of packages that makes Julia the easiest language to get started for actuarial workflows.
modeling
diffeq
mortalitytables
actuaryutilities
tutorial

Introduction

This demonstrates an example of defining an universal life policy value roll like a discrete differential equation.

It then uses the SciML DifferentialEquations package to “solve” the policy projection given a single point, but also to see how the policy projection behaves under different premium and interest rate conditions.

The first time this is run, it will download and precompile some large packages.

using Dates
using MortalityTables
using DifferentialEquations
using Plots
using ColorSchemes # for Turbo colors, which emphasize readability
using ActuaryUtilities
using LaTeXStrings

Let’s use the 2001 CSO table as the basis for cost of insurance charges:

cso = MortalityTables.table("2001 CSO Super Preferred Select and Ultimate - Male Nonsmoker, ANB")
MortalityTable (CSO/CET):
   Name:
       2001 CSO Super Preferred Select and Ultimate - Male Nonsmoker, ANB
   Fields: 
       (:select, :ultimate, :metadata)
   Provider:
       American Academy of Actuaries
   mort.SOA.org ID:
       1076
   mort.SOA.org link:
       https://mort.soa.org/ViewTable.aspx?&TableIdentity=1076
   Description:
       2001 Commissioners Standard Ordinary (CSO) Super Preferred Select and Ultimate Table – Male Nonsmoker. Basis: Age Nearest Birthday. Minimum Select Age: 0. Maximum Select Age: 99. Minimum Ultimate Age: 16. Maximum Ultimate Age: 120

Next, policy mechanics are coded. It’s essentially a discrete differential equation, so it leverages DifferentialEquations.jl

The projection is coded in the Discrete DifferentialEquation format:

\[u_{n+1} = f(u,p,t_{n+1})\]

In the code below, this translates to:

  • u is the state of the system. We will track three variables to represent the state:
    • state[1] is the account value
    • state[2] is the premium paid
    • state[3] is the policy duration
  • p are the parameters of the system.
  • t is the time, which will represent days since policy issuance
function policy_projection(state, p, t)
    # grab the state from the inputs
    av = state[1]

    # calculated variables
    cur_date = p.issue_date + Day(t)
    dur = duration(p.issue_date, cur_date)
    att_age = p.issue_age + dur - 1

    # lapse if AV <= 0 
    lapsed = (av <= 0.0) & (t > 1)

    if !lapsed

        monthly_coi_rate = (1 - (1 - p.mort_assump[att_age])^(1 / 12))

        ## Periodic Policy elements

        # annual events
        if Dates.monthday(cur_date) == Dates.monthday(p.issue_date) ||
           cur_date == p.issue_date + Day(1) # OR first issue date

            premium = p.annual_prem
        else
            premium = 0.0
        end

        # monthly_events
        if Dates.day(cur_date) == Dates.day(p.issue_date)
            coi = max((p.face - av) * monthly_coi_rate, 0.0)
        else
            coi = 0.0
        end

        # daily events
        int(av) = av * ((1 + p.int_rate)^(1 / 360) - 1.0)



        # av
        new_av = max(0.0, av - coi + premium + int(av - coi))

        # new state
        return [new_av, premium, dur] # AV, Prem, Dur

    else
        # new state
        return [0.0, 0.0, dur] # AV, Prem, Dur

    end


end
policy_projection (generic function with 1 method)

The following function will create a named tuple of parameters given a varying prem (premium) and int (credit rate).

params(prem, int) = (
    int_rate=int,
    issue_date=Date(2010, 1, 1),
    face=1e6,
    issue_age=25,
    mort_assump=MortalityTables.table("2001 CSO Super Preferred Select and Ultimate - Male Nonsmoker, ANB").ultimate,
    projection_years=75,
    annual_prem=prem,
)
params (generic function with 1 method)

Runing the system

This results in the following plot. The tracked output variables u1 and u2 represent the two vars that we tracked above: account value and cumulative premium.

begin
    p = params(
        8000.0, # 8,000 annual premium
        0.08    # 8% interest
    )

    # calculate the number of days to project
    projection_end_date = p.issue_date + Year(p.projection_years)
    days_to_project = Dates.value(projection_end_date - p.issue_date)

    # the [0.0,..] are the initial conditions for the tracked variables
    prob = DiscreteProblem(policy_projection, [0.0, 0.0, 0], (0, days_to_project), p)
    proj = solve(prob, FunctionMap())

    plot(proj)
end