###############################################################################
#
# Copyright (c) Longevtas Ltd. All rights reserved. This script is provided
# without warranty of any kind. Updated versions of this script can be found
# at https://www.longevitas.co.uk/site/informationmatrix/portfoliomortalitytrackingusav.uk.html
#
# This script calculates the real-time mortality tracker as per Richards, S. J.
# (2022) Real-time measurement of portfolio mortality levels in the presence of
# shocks and reporting delays, Annals of Actuarial Science,
# doi: 10.1017/S1748499522000021
#
# The CSV input file is presumed to follow the format generated by the 'Data
# Audit' feature of Longevitas, i.e. with three columns with a header row
# containing (i) calendar time, (ii) lives (or policies) and (iii) deaths
# (or cessations). Rows are assumed sorted by ascending time. The data format
# is shown in Figure 1 of the accompany blog.
#
# The generation of this CSV input file is computationally intensive, with the
# number of comparisons increasing with the square of the number of annuities.
# For this reason, the CSV input file is generated using custom C++ code using
# parallel processing. Background material on this topic can be found here:
#
# https://www.longevitas.co.uk/site/informationmatrix/followingthethread.html
#
# Note that this script does not implement the allowance for occurred-but-not-
# reported (OBNR) deaths discussed in the following blog:
#
# https://www.longevitas.co.uk/site/informationmatrix/hereisthenowcast.html
#
# However, such 'nowcasting' is available using the spreadsheet automatically
# generated by the 'Data Audit' feature in Longevitas.
# Smoothing constant to apply.
c <- 0.15
# Read in CSV data file.
d <- read.csv(file="C:/sandpit/US.csv", header=TRUE, skip=1)
# Calculate Nelson-Aalen cumulative hazard in time as per equation (1).
vdNelsonAalen <- cumsum(d$Deaths/d$Lives)
# Calculate smoothed hazard in time as per equation (2).
n <- length(d$Time)
vdMu <- rep(NA, n)
for (y in 2:n)
{
interval <- findInterval(c(d$Time[y]-c/2, d$Time[y]+c/2), d$Time)
if (interval[1]*interval[2] > 0)
vdMu[y] <- (vdNelsonAalen[interval[2]] - vdNelsonAalen[interval[1]]) / c
}
# Plot c-smoothed hazard. Adjust the xlim and ylim ranges according to your data.
par(mar=c(2, 3, 0, 0), col="darkblue", col.axis="darkblue", family="serif")
plot(d$Time, vdMu, xlim=range(2015, 2021.5), ylim=range(0.0, 0.10), type="n", axes=FALSE, xlab="", ylab="")
axis(1, las=1)
axis(2, las=1)
lines(d$Time, vdMu)
legend("topleft", legend=paste0("Smoothed with c=", c), bty="n", cex=1.2)