R/LCLS.R

In demofit: Parametric Mortality Curve Fitting and Mortality Forecasting Tools

#' Lee-Carter model for limited data
#'
#' Fits and forecasts mortality rates using Lee-Carter model with sparse data in irregular years.
#'
#' @param x vector of ages.
#' @param t vector of years.
#' @param M matrix of mortality rates (rows as years and columns as ages).
#' @param curve name of mortality curve for smoothing forecasted mortality rates (including gompertz, makeham, oppermann, thiele, wittsteinbumsted, perks, weibull, vandermaen, beard, heligmanpollard, rogersplanck, siler, martinelle, thatcher, gompertz2, makeham2, oppermann2, thiele2, wittsteinbumsted2, perks2, weibull2, vandermaen2, beard2, heligmanpollard2, rogersplanck2, siler2, martinelle2, thatcher2, where first 14 curves' parameters are unconstrained and last 14 curves' parameters are generally restricted to be positive).
#' @param h forecast horizon (default = 10).
#' @param jumpoff if 1, forecasts are based on estimated parameters only; if 2, forecasts are anchored to observed mortality rates in final year (default = 1). 
#'
#' @details
#' The Lee-Carter (LC) model is specified as 
#' 
#' \eqn{ln(m_{x,t}) = \alpha_x + \beta_x \kappa_t + \epsilon_{x,t}}.
#'
#' The model is estimated by singular value decomposition and is forecasted by random walk with drift applied to \eqn{\kappa_t}. Constraints include sum of \eqn{\beta_x} is one and sum of \eqn{\kappa_t} is zero. It can be applied to whole age range.
#'
#' @importFrom stats fitted prcomp sd rnorm
#' @importFrom graphics par lines legend points image 
#' @importFrom grDevices colorRampPalette
#'
#' @return
#' An object of class LCLS with associated S3 methods coef, forecast (which = 1 for smoothed (default); which = 2 for raw), plot, residuals, and simulate (nsim for setting number of simulations; seed for initialising random number generator).
#'
#' @references
#' Li, N., Lee, R., and Tuljapurkar, S. (2004). Using the Lee-Carter method to forecast mortality for populations with limited data. International Statistical Review, 72(1), 19-36.
#'
#' @examples
#' x <- 60:89
#' t <- c(1991,1996,2001,2006,2011:2020)
#' a <- c(-4.8499,-4.7676,-4.6719,-4.5722,-4.4847,-4.3841,-4.2813,-4.1863,-4.0861,-3.9962,
#' -3.8885,-3.7896,-3.6853,-3.5737,-3.4728,-3.3718,-3.2586,-3.1474,-3.0371,-2.9206,
#' -2.7998,-2.6845,-2.5653,-2.4581,-2.3367,-2.2159,-2.1017,-1.9941,-1.8821, -1.7697)
#' b <- c(0.0283,0.0321,0.0335,0.0336,0.0341,0.0358,0.0368,0.0403,0.0392,0.0395,
#' 0.0396,0.0399,0.0397,0.0386,0.039,0.0375,0.0367,0.0368,0.035,0.0354,
#' 0.0336,0.0323,0.0313,0.0295,0.0282,0.0265,0.024,0.0226,0.0219,0.0183)
#' k <- c(12.11,8.21,
#' 3.27,-1.03,
#' -5.18,-5.64,-6,-6.51,-6.91,-6.9,-8.32,-8.53,-9.69,-9.31)
#' set.seed(123)
#' M <- exp(outer(k,b)+matrix(a,nrow=14,ncol=30,byrow=TRUE)+rnorm(420,0,0.035))
#' fit <- LCLS(x=x,t=t,M=M,curve="makeham",h=30,jumpoff=2)
#' coef(fit)
#' forecast::forecast(fit)
#' plot(fit)
#' residuals(fit)
#'
#' @export
LCLS <- function(x,t,M,curve=c("gompertz","makeham","oppermann","thiele","wittsteinbumsted","perks","weibull","vandermaen","beard","heligmanpollard","rogersplanck","siler","martinelle","thatcher","gompertz2","makeham2","oppermann2","thiele2","wittsteinbumsted2","perks2","weibull2","vandermaen2","beard2","heligmanpollard2","rogersplanck2","siler2","martinelle2","thatcher2"),h=10,jumpoff=1) {
if (!is.numeric(x)||!is.numeric(t)||!is.numeric(M)) { stop("x and t and M must be numeric") }
if (!is.vector(x)||!is.vector(t)) { stop("x and t must be a vector") }
if (!is.matrix(M)) stop("M must be a matrix with its rows as years and columns as ages")
if (length(x)!=ncol(M)) stop("the number of ages must match the number of columns of M")
if (length(t)!=nrow(M)) stop("the number of years must match the number of rows of M")
if (is.unsorted(x,strictly=TRUE)||is.unsorted(t,strictly=TRUE)) { stop("x and t must be in ascending order") }
if (any(x<0)||any(t<0)) { stop("x and t must be non-negative") }
if (any(M<=0)) { stop("all M values must be positive") }
if (nrow(M)<3) stop("it requires at least 3 years of data for this forecast")
if (max(t)-min(t)+1<20) stop("it requires at least a period of 20 years of data for this forecast")
if (!is.numeric(h)||!is.numeric(jumpoff)) { stop("h and jumpoff must be numeric") }
if (h<1) { stop("h must be at least 1") }
if (jumpoff!=1&&jumpoff!=2) { stop("jump-off must be either 1 or 2") }
curve <- tryCatch(match.arg(curve),error = function(e) { stop("invalid curve choice") })
tryCatch({
nr <- nrow(M); nc <- ncol(M)
PCA <- prcomp(log(M),center=TRUE,scale=FALSE)
a <- numeric(); for (j in 1:nc) { a[j] <- mean(log(M[,j])) }
b <- PCA$rotation[,1]/sum(PCA$rotation[,1])
k <- PCA$x[,1]*sum(PCA$rotation[,1])
res <- array(NA,c(nr,nc))
for (i in 1:nr) { for (j in 1:nc) { res[i,j] <- log(M[i,j])-a[j]-b[j]*k[i] }}
res <- res/sd(res)
mu <- (k[nr]-k[1])/(t[nr]-t[1])
kf <- k[nr]+mu*(1:h)
Mf <- array(NA,c(h,nc)); Mfs <- array(NA,c(h,nc))
if (jumpoff==1) { for (i in 1:h) { for (j in 1:nc) { Mf[i,j] <- exp(a[j]+b[j]*kf[i]) }}}
if (jumpoff==2) { for (i in 1:h) { for (j in 1:nc) { Mf[i,j] <- M[nr,j]*exp(b[j]*(kf[i]-k[nr])) }}}
for (i in 1:h) { Mfs[i,] <- fitted(MC(x=x,m=Mf[i,],curve=curve)) }
invisible(structure(
list(curve=curve,x=x,t=t,M=M,h=h,jumpoff=jumpoff,alpha=a,beta=b,kappa=k,standardresiduals=res,forecast=Mf,smoothforecast=Mfs),
class="LCLS"
))
}, error = function(e) { stop(paste0("model fitting and forecasting are unsuccessful - please make sure the data and age range are suitable for the model and curve\n",e$message),call.=FALSE) })
}

#' @export
coef.LCLS <- function(object,...) {
list(alpha=object$alpha,beta=object$beta,kappa=object$kappa)
}

#' @export
forecast.LCLS <- function(object,which=1,...) {
if (length(which)!=1||!(which%in%c(1,2))) { stop("which must be 1 or 2") }
if (which==1) { object$smoothforecast }
else if (which==2) { object$forecast }
}

#' @export
plot.LCLS <- function(x,...) {
old <- par(no.readonly=TRUE)
on.exit(par(old))
par(mfrow=c(3,2),mar=c(3.5,2.5,1.5,0.5),mgp=c(1.5,0.5,0))
plot(x$x,x$alpha,xlab="age",ylab="alpha",bty="n",type="l")
plot(x$x,x$beta,xlab="age",ylab="beta",bty="n",type="l")
plot(x$t,x$kappa,xlab="year",ylab="kappa",bty="n",type="l")
colband <- colorRampPalette(c("black","grey","white"))
image(x=x$x,y=x$t,z=t(x$standardresiduals),col=colband(6),xlab="age",ylab="year",main="standardised residuals",cex.main=1,font.main=1)
plot(x$x,log(x$M[1,]),xlab="age",ylab="log death rate",pch=16,cex=0.5,bty="n",ylim=c(min(log(x$M),log(x$smoothforecast)),max(log(x$M),log(x$smoothforecast))))
points(x$x,log(x$M[nrow(x$M),]),pch=1,cex=0.5)
lines(x$x,log(x$smoothforecast[x$h,]))
temp <- paste("forecast",x$h,"years")
legend("bottomright",legend=c("observed first data","observed last data",temp),pch=c(16,1,NA),lty=c(NA,NA,1),pt.cex=0.5,cex=0.8,bty="n")
}

#' @export
residuals.LCLS <- function(object,...) {
object$standardresiduals
}

#' @export
simulate.LCLS <- function(object,nsim=10,seed=123,...) {
if (!is.numeric(nsim)||nsim!=floor(nsim)||nsim<1||!is.numeric(seed)||seed!=floor(seed)||seed<1) { stop("nsim and seed must be positive integers") }
t <- object$t; M <- object$M; h <- object$h; jumpoff <- object$jumpoff
a <- object$alpha; b <- object$beta; k <- object$kappa
nr <- nrow(M); nc <- ncol(M)
amat <- matrix(a,nr,nc,byrow=TRUE)
bmat <- matrix(b,nr,nc,byrow=TRUE)
kmat <- matrix(k,nr,nc,byrow=FALSE)
res <- log(M)-amat-bmat*kmat
s <- sd(res)
mu <- (k[nr]-k[1])/(t[nr]-t[1])
s2 <- sum((diff(k)-mu*diff(t))^2)/(t[nr]-t[1]-sum(diff(t)^2)/(t[nr]-t[1]))
Msim <- array(NA,c(h,nc,nsim))
am <- matrix(a,h,nc,byrow=TRUE)
bm <- matrix(b,h,nc,byrow=TRUE)
set.seed(seed)
if (jumpoff==1) {
for (z in 1:nsim) {
ksim <- k[nr]+cumsum(rnorm(h,mean=mu,sd=sqrt(s2)))
km <- matrix(ksim,h,nc,byrow=FALSE)
Msim[,,z] <- exp(am+bm*km+matrix(rnorm(h*nc,0,s),h,nc,byrow=TRUE))
}} 
if (jumpoff==2) {
for (z in 1:nsim) {
ksim <- k[nr]+cumsum(rnorm(h,mean=mu,sd=sqrt(s2)))
km <- matrix(ksim,h,nc,byrow=FALSE)
Msim[,,z] <- exp(matrix(log(M[nr,]),h,nc,byrow=TRUE)+bm*(km-matrix(k[nr],h,nc,byrow=FALSE))+matrix(rnorm(h*nc,0,s),h,nc,byrow=TRUE))
}}
list(Msim=Msim,sigma=s)
}