# Aging: Lee-Carter: Estimation

• The directory /tmp/drutex-598e650a042c91e7e969e49ca8ba2ed9-1 has been created.
• The directory /tmp/drutex-598e650a042c91e7e969e49ca8ba2ed9-2 has been created.

This Lee-Carter exercise is meant to give you a rough idea for the overall process used in fitting log-mortality rates with the SVD (dimensionality reduction), obtaining the representation in terms of the univariate time factor TeX Embedding failed! and the two age profiles TeX Embedding failed!.

```# use HMD central death rates, 5 year age groups, 1959-2002
file=url("http://coale.stanford.edu/data/HMD.USA19592002nMx.txt"),

# 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)```