robjhyndman / demography

demography package for R
https://pkg.robjhyndman.com/demography
73 stars 25 forks source link

LCA Function #19

Closed PalDaw closed 8 years ago

PalDaw commented 8 years ago

I am currently researching on the Lee-Carter(1992) model and hoping to make contribution to knowledge in the area of mortality forecasting. I went through the LCA function in R and I was able to produce a basic R summary codes for estimating the components of the Lee-Carter model (ax, bx and kt) . However , I am wondering whether there is anything I am missing. My results for "ax" seem to match with values from LCA function, however the values for ,"bx" and "kt" have significant differences . Could you kindly lend me some time to run the attached R codes. Is there anything I am missing? . I would appreciate very well your input and suggestions.

REPLICATION_LCA FUNC.docx

PalDaw commented 8 years ago

usa<- hmd.mx("USA", "USERNAME", "PASSWORD","USA") ######################################### usa1 <- extract.years(usa, years =1950:1959) mx <- usa1$rate$female mx <- t(mx) logrates <- log(mx) ax <- apply(logrates,2,mean) # ax is mean of logrates by column clogrates <- sweep(logrates,2,ax) # central log rates (with ax subtracted) (dimensions mn) svd.mx <- svd(clogrates) sumv <- sum(svd.mx$v[,1]) bx <- svd.mx$v[,1]/sumv kt <- svd.mx$d[1] * svd.mx$u[,1] \ sumv ##################################################

Fitting the model

clogratesfit <- outer(kt, bx) logratesfit <- sweep(clogratesfit,2,ax,"+") par(mfrow=c(1,2)) fit <- usa1 fit$y <- matrix(sapply(logratesfit, exp), ncol=2, byrow=TRUE) fit$x <- 0:110 plot(fit) #################################################

USING THE LCA FUNCTION

lc.female <- lca(usa1, series="female",ages=0:110) plot(lc.female$fitted) ###############################################

Comparing the output

lc.female$kt kt lc.female$bx bx lc.female$ax ax

Running the code above produced the following results for "kt"

lc.female$kt #from the "lca function"# Time Series: Start = 1950 End = 1959 Frequency = 1 [1] 11.386593 9.677363 6.597937 3.805083 -4.028068 -4.078909 -5.359104 [8] -3.214466 -5.607462 -8.974393 kt #from my code# [1] 5.415389 6.685072 6.689227 -1.464312 -1.885153 -2.889145 -1.989212 [8] -2.997783 -3.446872 -4.117211

Could you advise on the reason for the difference.

robjhyndman commented 8 years ago

You haven't applied the deaths adjustment described by Lee and Carter.