JuliaActuary
Practical, extensible, and open-source actuarial modeling and analysis.

Universal Life Policy Account Mechanics as a Differential Equation

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.

begin
    using Dates 
    using MortalityTables 
    using DifferentialEquations 
    using Plots 
    using ActuaryUtilities
end

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:

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

Moving up the ladder of abstraction

An excellent way to understand the behavior of a model is to move up the ladder of abstraction. Below, we will see what happens to the projection at varying levels of credit rates and annual premiums.

The full range of simulated outcomes can take a couple of minutes to run:

begin
    prem_range = 1000.0:100.0:9000.0
    int_range = 0.02:0.0025:0.08
    
    function ending_av(ann_prem,int,days_to_project)
        p = params(ann_prem,int)
        prob = DiscreteProblem(policy_projection,[0.0,0.0,0],(0,days_to_project),p)
        proj = solve(prob,FunctionMap())
        end_av = proj[end][1] 
        if end_av == 0.0
            lapse_time = findfirst(isequal(0.0),proj[1,2:end])
        else
            lapse_time = length(proj)
        end
        duration = proj[3,lapse_time]
        end_age = p.issue_age + duration - 1.0
        
        return end_av,end_age
        
    end
    
    end_age = zeros(length(prem_range),length(int_range))
    end_av = zeros(length(prem_range),length(int_range))
    
    # loop through each projection we did and fill our ranges with the ending AV and ending age
    for (i,vp) in enumerate(prem_range)
        for (j,vi) in enumerate(int_range)
            end_av[i,j],end_age[i,j] = ending_av(vp,vi,days_to_project)
        end
    end
end

Now let's plot the result. Not surprising, interest has a huge effect on the policy projection. Premium is also a major influence.

One thing that's remarkable is how going from 2000 premium to just ~2200 of premium results in about a 5m difference at 8% interest. The power of compound interest!

begin
    using ColorSchemes # for Turbo colors, which emphasize readability
    
    viz = plot(layout=2) # side by side plot
    
    # 
    contour!(viz[1],int_range,
        prem_range,
        end_av ./ 1e6, # scale to millions for readability
    
        contour_labels=true,
        c=cgrad(ColorSchemes.turbo.colors),
        fill=true,
        title="AV at age 100 (M)",
        ylabel=L"Annual Premium (\$)",
        xlabel="credit rate"
        
    
    )
    
    contour!(viz[2],int_range,
        prem_range,
        end_age, 
    
        contour_labels=true,
        c=cgrad(ColorSchemes.turbo.colors),
        fill=true,
        yaxis=false,
        title="Age at Lapse",
        xlabel="credit rate"
    )
    
    annotate!(viz[2],[0.055,7000,Plots.text("Doesn't lapse \nbefore age 100", 8, :white, :center)])
    
end
using LaTeXStrings

Conclusion

This shows how universal life mechanics are a dynamic system. the growth/decay is governed by two competing feedback loops:

Endnotes

This is not meant to represent any particular insurance product, nor fully replicate typical account mechanics.

Built with Julia 1.8.5 and

ActuaryUtilities 3.4.2
ColorSchemes 3.19.0
DifferentialEquations 7.2.0
LaTeXStrings 1.3.0
MortalityTables 2.3.0
Plots 1.31.2

Run this Pluto Notebook

To run this page locally, download this file and open it with Pluto.jl.

The packages in JuliaActuary are open-source and liberally licensed (MIT License) to allow wide private and commercial usage of the packages, like the base Julia language and many other packages in the ecosystem. See terms of this site.