Nothing
# Model North --------------------------------------------------
# With correction by Jim Oeppen (September 2018)
cdmltn <- function(sex = "F"){
if (sex != "F" & sex !="M")
stop("sex must be either F or M!")
if(sex == "F") { # Female
# CD lifetables indexed by e_10 = eten
# Level 1 is ef=20, Level 25 is ef = 80
# JO optimisation results
etenf <- c( 22.40195, 26.40432, 29.99710, 33.25105, 36.21856,
38.94234, 41.45523, 43.78386, 45.95114, 47.97472,
49.87049, 51.64820, 52.94942, 54.08511, 55.28232,
56.53927, 57.85128, 59.21201, 60.61233, 62.04218,
63.49014, 64.96932, 67.39459, 70.66124, 74.81290)
eten <- etenf
# Unit vector of same length as eten
eee <- 1 + 0*eten
# ages at start of groups; later extended.
xx <- c(0, 1, seq(5,75, by=5) )
# FEMALE ax and bx are the coefficients of the simple regressions of
# the nqx's on e10
# JO correction age 5
ax <- c(0.47504, 0.45025, 0.19376, 0.10041, 0.10126, 0.11261,
0.13137, 0.15448, 0.17693, 0.1844, 0.1944, 0.22364,
0.30043, 0.41033, 0.56691, 0.77206, 0.96175)
bx <- c(-0.006923, -0.006805, -0.002928, -0.001497, -0.001480, -0.001618, -0.001893, -0.002239, -0.002566, -0.002612, -0.002712, -0.003011, -0.004053, -0.005394, -0.007187, -0.009334, -0.010681)
# a1x and b1x are the coefficients of the log-regressions of the nqx's on e10
a1x <- c(5.7332, 7.6298, 7.1271, 6.1089, 5.4984, 5.2649, 5.2547, 5.3691, 5.3186, 4.9099, 4.6164, 4.3673, 4.4363, 4.4163, 4.4030, 4.3826, 4.3108)
# JO correction age 20
b1x <- c(-0.05133, -0.08909, -0.08647, -0.07192, -0.05955,
-0.05372, -0.05236, -0.05339, -0.05136, -0.04261,
-0.03627, -0.02961, -0.02858, -0.02511, -0.02152,
-0.01784, -0.01355)
#####################################################################
##### # Calculate an array of probabilities of dying nqx up to age 75
#######################################################################
# The first subscript refers to e10; the second to the age x.
nqxa <- eee %o% ax + eten %o% bx # nqx from linear model
nqxl <- (1/10000)*10^( eee %o% a1x + eten %o% b1x ) # exponential model
#
# Choose among linear, mixed, and exponential models
# on basis of crossover points and slopes.
slopechk <- ( (0*eten + 1/log(10)) %o% (bx/b1x) ) * (1/nqxl ) # Slope check
nqx <- ( nqxa + nqxl) / 2
nqx <- ifelse( nqxa <= nqxl & slopechk < 1, nqxa, nqx )
nqx <- ifelse( nqxa <= nqxl & slopechk >= 1, nqxl, nqx )
####################################################################
#### Calculate the survivorship proportions lx up to age 75
##################################################################
u <- log( 1 - nqx)
uu <- t( apply( u, c(1), cumsum ) )
km1 <- length( uu[1,]) - 1 # number of age groups minus 1
lx <- cbind( 1+0*eee, exp( uu[ , 1:km1]) )
#####################################################################
# Special extension from ages 80 to 100 (i.e. for j> 18
###############################################################
nL75 <- 5*lx[,17]*(1-nqx[,17]) + (1/2)*5*lx[,17]*nqx[,17]
leighty <- lx[,17]*(1-nqx[,17]) # Survivorship to age 80
m77 <- lx[,17]*nqx[,17]/nL75 # m77 is hazard age 77.5.
m105 <- 0.613+ 1.75*nqx[,17]
kappa <- (1/27.5)*log( m105 / m77 ) # THIS CORRECTS COALE'S MISTAKE
m80 <- m77*exp(kappa*2.5) # Hazard at age 80
xlong <- seq(80,121, by = 0.2 ) - 80 # ages in 1/5 years -80
uu <- kappa %o% xlong
u <- (-m80/kappa) %o% ( 1 + 0*xlong)
llx <- leighty %o% ( 1 + 0*xlong)
llx <- llx * exp ( u*exp(uu) - u )
nL80 <- llx[, 1:26 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL85 <- llx[, 26:51 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL90 <- llx[, 51:76 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL95 <- llx[, 76:101 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
T100 <- as.vector( llx[, 101:201] %*% c( 1/2, rep(1,99), 1/2)*(1/5))
#################################################################
##### Compute other functions in the lifetable
#################################################################
xx <- c(0, 1, seq(5,95, by=5) ) # ages
nn <- c( diff(xx),8 ) # age group widths
lx <- cbind( lx[ , 1:17], llx[ , c(1,26,51,76)] )
qqx <- ( llx[ , c(1,26,51) ] - llx[ , c(26,51,76)] )/llx[, c(1,26,51)]
nqx <- cbind( nqx[ , 1:17], qqx, eee )
dimnames(nqx)[[1]] <- paste(seq(eten))
dimnames(nqx)[[2]] <- paste(xx)
dimnames(lx) <- dimnames(nqx)
na0 <- 0.050 + 3.00*nqx[,1]
na0 <- ifelse( nqx[,1]< 0.100, na0, 0.35 )
# JO corrected na1
na1 <- 1.733 - 1.627 * nqx[, 1]
na1 <- ifelse(nqx[, 1] < 0.1, na1, 1.570)
nax <- cbind( na0, na1, 2.25*eee)
nax <- cbind( nax, 2.6*eee %o% rep(1,13), 2.5*eee)
nax <- cbind( nax, 2.5*eee %o% rep(1,4) ) # Extend to length 21
dimnames(nax) <- dimnames(nqx)
ndx <- lx*nqx
nLx <- ( eee %o% nn )*( lx - ndx ) + nax*ndx
nLx <- cbind( nLx[ , 1:17], nL80, nL85, nL90, nL95 )
# JO add dimnames
dimnames(nLx)[[1]] <- paste(seq(eten))
dimnames(nLx)[[2]] <- paste(xx)
nmx <- ndx/nLx
K <- rev(seq(xx)) # Reversed indices for age groups
Tx <- t( apply( nLx[, K], c(1), cumsum) )
Tx <- ( T100 %o% rep(1,21) ) + Tx[,K]
ex <- Tx/lx
out <- list(age = xx, width = nn, e10 = eten,
lx=lx,nqx=nqx,nax=nax,ndx=ndx,nLx=nLx,nmx=nmx,Tx=Tx,ex=ex)
return(out)
} # end Female part
# Male
# etenm is a linear function of etenf from runs of Coale-Demeny female
# program etenm <- 2.298048 + 0.914067*etenf
# JO optimisation results
etenm <- c(20.82759, 24.69877, 28.17353, 31.32028, 34.19063,
36.82507, 39.25501, 41.50779, 43.60369, 45.56088,
47.39436, 49.11385, 50.37175, 51.47042, 52.62865,
53.84409, 55.11337, 56.42968, 57.78417, 59.16709,
60.56808, 61.99947, 64.34571, 67.50589, 71.52278)
eten <- etenm #
eee <- 1 + 0*eten # Unit vector of same length as eten
xx <- c(0, 1, seq(5,75, by=5) ) # ages at start of groups; later extended.
# MALE ax and bx are the coefficients of the simple regressions of the
# nqx's on e10
ax <- c(0.54327, 0.46169, 0.18983, 0.09551, 0.09666, 0.13472, 0.14325, 0.15280, 0.17535, 0.20924, 0.24673, 0.28578, 0.36171, 0.45849, 0.59986, 0.82662, 1.03681)
# JO correction age 60
bx <- c(-0.008251, -0.00729, -0.002974, -0.001476, -0.001422,
-0.001968, -0.002103, -0.002244, -0.002589, -0.003083,
-0.003605, -0.004016, -0.005037, -0.006124, -0.007677,
-0.010241, -0.011906)
# MALE a1x and b1x are the coefficients of the log-regressions of the
# nqx's on e10
a1x <- c(5.6151, 7.2025, 6.1947, 5.3488, 4.5662, 4.6970, 4.7661, 4.7248, 4.7568, 4.7280, 4.6020, 4.3499, 4.3718, 4.2977, 4.2858, 4.3482, 4.3197)
b1x <- c(-0.05022, -0.08475, -0.07195, -0.06047, -0.04322, -0.04277, -0.04372, -0.04236, -0.04197, -0.03986, -0.03578, -0.02857, -0.02682, -0.02244, -0.01913, -0.01710, -0.01357)
#####################################################################
##### # Calculate an array of probabilities of dying nqx up to age 75
#######################################################################
# The first subscript refers to e10; the second to the age x.
nqxa <- eee %o% ax + eten %o% bx # nqx from linear model
nqxl <- (1/10000)*10^( eee %o% a1x + eten %o% b1x ) # exponential model
#
# Choose among linear, mixed, and exponential models
# on basis of crossover points and slopes.
slopechk <- ( (0*eten + 1/log(10)) %o% (bx/b1x) ) * (1/nqxl ) # Slope check
nqx <- ( nqxa + nqxl) / 2
nqx <- ifelse( nqxa <= nqxl & slopechk < 1, nqxa, nqx )
nqx <- ifelse( nqxa <= nqxl & slopechk >= 1, nqxl, nqx )
####################################################################
#### Calculate the survivorship proportions lx up to age 75
##################################################################
u <- log( 1 - nqx)
uu <- t( apply( u, c(1), cumsum ) )
km1 <- length( uu[1,]) - 1 # number of age groups minus 1
lx <- cbind( 1+0*eee, exp( uu[ , 1:km1]) )
#####################################################################
# Special extension from ages 80 to 100 (i.e. for j> 18
###############################################################
nL75 <- 5*lx[,17]*(1-nqx[,17]) + (1/2)*5*lx[,17]*nqx[,17]
leighty <- lx[,17]*(1-nqx[,17]) # Survivorship to age 80
m77 <- lx[,17]*nqx[,17]/nL75 # m77 is hazard age 77.5.
# JO corrected
m105 <- 0.551 + 1.75 * nqx[, 17] # Male
kappa <- (1/27.5)*log( m105 / m77 ) # THIS CORRECTS COALE'S MISTAKE
m80 <- m77*exp(kappa*2.5) # Hazard at age 80
xlong <- seq(80,121, by = 0.2 ) - 80 # ages in 1/5 years -80
uu <- kappa %o% xlong
u <- (-m80/kappa) %o% ( 1 + 0*xlong)
llx <- leighty %o% ( 1 + 0*xlong)
llx <- llx * exp ( u*exp(uu) - u )
nL80 <- llx[, 1:26 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL85 <- llx[, 26:51 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL90 <- llx[, 51:76 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
nL95 <- llx[, 76:101 ] %*% c( 1/2, rep(1,24), 1/2)*(1/5)
T100 <- as.vector( llx[, 101:201] %*% c( 1/2, rep(1,99), 1/2)*(1/5))
#################################################################
##### Compute other functions in the lifetable
#################################################################
xx <- c(0, 1, seq(5,95, by=5) ) # ages
nn <- c( diff(xx),8 ) # age group widths
lx <- cbind( lx[ , 1:17], llx[ , c(1,26,51,76)] )
qqx <- ( llx[ , c(1,26,51) ] - llx[ , c(26,51,76)] )/llx[, c(1,26,51)]
nqx <- cbind( nqx[ , 1:17], qqx, eee )
dimnames(nqx)[[1]] <- paste(seq(eten))
dimnames(nqx)[[2]] <- paste(xx)
dimnames(lx) <- dimnames(nqx)
# JO corrected na0 and na1
na0 <- 0.0425 + 2.875 * nqx[, 1]
na0 <- ifelse(nqx[, 1] < 0.1, na0, 0.33)
na1 <- 1.859 - 3.013 * nqx[, 1]
na1 <- ifelse(nqx[, 1] < 0.1, na1, 1.558)
nax <- cbind( na0, na1, 2.25*eee)
nax <- cbind( nax, 2.6*eee %o% rep(1,13), 2.5*eee)
nax <- cbind( nax, 2.5*eee %o% rep(1,4) ) # Extend to length 21
dimnames(nax) <- dimnames(nqx)
ndx <- lx*nqx
nLx <- ( eee %o% nn )*( lx - ndx ) + nax*ndx
nLx <- cbind( nLx[ , 1:17], nL80, nL85, nL90, nL95 )
# JO add dimnames
dimnames(nLx)[[1]] <- paste(seq(eten))
dimnames(nLx)[[2]] <- paste(xx)
nmx <- ndx/nLx
K <- rev(seq(xx)) # Reversed indices for age groups
Tx <- t( apply( nLx[, K], c(1), cumsum) )
Tx <- ( T100 %o% rep(1,21) ) + Tx[,K]
ex <- Tx/lx
out <- list(age = xx,width = nn,lx = lx, e10 = eten,
nqx=nqx,nax=nax,ndx=ndx,nLx=nLx,nmx=nmx,Tx=Tx,ex=ex)
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.