life_modeling_problem.md

library(data.table)

q <- c(0.001,
       0.002,
       0.003,
       0.003,
       0.004,
       0.004,
       0.005,
       0.007,
       0.009,
       0.011)

w <- c(0.05,
       0.07,
       0.08,
       0.10,
       0.14,
       0.20,
       0.20,
       0.20,
       0.10,
       0.04)

P <- 100
S <- 25000
r <- 0.02

dt <- as.data.table(cbind(q,w))

npv <- function(cf, r, S, P) {
  cf[, inforce := shift(cumprod(1 - q - w), fill = 1)
  ][, lapses := inforce * w
  ][, deaths := inforce * q
  ][, claims := deaths * S
  ][, premiums := inforce * P
  ][, ncf := premiums - claims
  ][, d := (1/(1+r))^(.I)
  ][, sum(ncf*d)]
}

npv(dt,r,S,P)
#> [1] 50.32483

microbenchmark::microbenchmark(npv(dt,r,S,P))
#> Unit: milliseconds
#>              expr    min      lq     mean  median      uq    max neval
#>  npv(dt, r, S, P) 2.5791 2.71035 2.964293 2.85625 3.10385 6.0357   100

^{Created on 2021-01-15 by the reprex package (v0.3.0)}

Here is the Rust implementation for reference:

pub fn npv(mortality_rates: &Vec<f64>, lapse_rates: &Vec<f64>, interest_rate: f64, sum_assured: f64, premium: f64, init_pols: f64, term: Option<usize>) -> f64 {

    let term = term.unwrap_or_else(|| mortality_rates.len());
    let mut result = 0.0;
    let mut inforce = init_pols;
    let v: f64 = 1.0 / (1.0 + interest_rate);

    for (t, (q, w)) in mortality_rates.iter().zip(lapse_rates).enumerate() {
        let no_deaths = if t < term {inforce * q} else {0.0};
        let no_lapses = if t < term {inforce * w} else {0.0};
        let premiums = inforce * premium;
        let claims = no_deaths * sum_assured;
        let net_cashflow = premiums - claims;
        result += net_cashflow * v.powi(t as i32);
        inforce = inforce - no_deaths - no_lapses;
    }

    result
}

And the benchmark:

npv                     time:   [362.23 ns 364.30 ns 366.35 ns]
Found 9 outliers among 100 measurements (9.00%)
  7 (7.00%) low mild
  1 (1.00%) high mild
  1 (1.00%) high severe

The reason there are three measurements above is that it's showing the upper and lower ranges. The Rust benchmarking system is a statistical one.

I wrote this fast so I may have made a mistake. I also did not try to optimize it so there's some performance gain to be had for sure. Even though it looks serial I'm sure the compiler is auto-vectorizing it.

Nice, @paddyhoran. I have not used rust so played with it here: https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=255a33e84c70e3bcd43258ddbaa66f27 Looks like the answer matches the other above.

Also looks like your algorithm for the problem is more efficient overall. It's over 2x faster than the prior Julia fastest Julia version above. Code and benchmark:

function npv3(q,w,P,S,r,term=nothing)
	term = isnothing(term) ? length(q) + 1 : term + 1
	inforce = 1.0
	result = 0.0
	
	for (t,(q,w)) in enumerate(zip(q,w))
		deaths = t < term ? inforce * q : 0.0
		lapses = t < term ? inforce * w : 0.0
		premiums = inforce * P
		claims = deaths * S
		ncf = premiums - claims
		result += ncf * (1 / ( 1 + r)) ^ (t-1)
		inforce = inforce - deaths - lapses
	end
	
  	return result
end

BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     135.023 ns (0.00% GC)
  median time:      135.262 ns (0.00% GC)
  mean time:        137.513 ns (0.00% GC)
  maximum time:     1.327 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     873

Updated Rankings for the "Life Modeling Problem":

code	algorithm	relative	absolute timing (mean)
Julia	accumulating for loop	1x	135 nanoseconds
Rust	accumulating for loop	3x	360 nanoseconds
Julia	vector	4x	500 nanoseconds
Python numpy	vector	155x	21 microseconds
R vectorized	vector	520x	70 microseconds
R data.table	data.table vector	22000x	3 milliseconds