Aging: Lee-Carter: Estimation

# use HMD central death rates, 5 year age groups, 1959-2002
HMD.nMx<- read.table(
       file=url("http://coale.stanford.edu/data/HMD.USA19592002nMx.txt"),
       header=T, skip=2) ;
 
# We only work with combined data here. Avoid sexual dimorphism in mortality evolution
# Put the rates into a time x age data matrix
nMx <- NULL;
for( year in 1959:2000){
  nMx <- rbind( nMx, HMD.nMx[  HMD.nMx$Year == year,5])
}
 
# Lee-Carter with SVD
# Note: the svd from R returns Vi Uj such that sum( v[i,]^2)=1 and sum( u[j,]^2)=1
# which differs from many normalizations but is standard for the SVD routines.
# NB that bx can have negative elements - so normalizing by sum(bx)=1 will not work out well.
# NB that the sign of bx and kt can be switched from what you might expect, so we check for
# the first bx value and change signs if it is negative
f.LCsvd=function(lm.t.x ){
  ax = apply(lm.t.x,2,mean);
  Y = sweep(lm.t.x, 2, ax);             #mean centered data
  Y.svd = svd(Y);                       #returns U %*% diag(d) %*% t(V) = Y
  bx = Y.svd$v[,1];
  b1sign=sign(bx[1]);
  kt = Y.svd$d[1]*Y.svd$u[,1];
  bx=bx*b1sign; kt=kt*b1sign;
  predictions = outer(kt,bx);           # time by age k[t] * b[x]
  predictions = sweep(predictions,2,ax,"+"); # add back in column means
  dimnames(predictions) = dimnames(lm.t.x);
  return(list(ax=ax, bx=bx, kt=kt,svd=Y.svd,predictions=predictions ));
};
 
lc <- f.LCsvd(log(nMx))
 
# take a look at results of the procedure
## k_t  ; mortality index parameter which will be forecast
plot(1959:2000,lc$kt,type="b",xlab="Years", ylab="")
title("kt mortality parameter from Lee-Carter SVD") 
 
## a_x  
ages <- sort( unique(HMD.nMx$Age))
plot(ages,lc$ax,type="b",xlab="Age",ylab="")      
title("ax mean log mortality")
## b_x
 plot(ages,lc$bx,type="b",xlab="Age",ylab="")      
 title("bx  age-specific sensitivity to mortality change")
#Q: what do negative b(x) mean?
 
## where does this model not fit?
persp( log(nMx) - lc$predictions)