# Author: tim
###############################################################################

stand <- function(x){
	x/sum(x)
}
mx2lx <- function(mx, trunc = 0){
	# pretend mx is continuous hazard, cheap trick
	# ought to work fine if we left truncate too
	x <- 1:length(mx) - 1
	mx <- mx[x >= trunc]

	c(1,exp(-cumsum(mx)))[1:length(mx)]
}

mx2dx <- function(mx, trunc = 0){
	x <- 1:length(mx) - 1
	mx <- mx[x >= trunc]
	# use previous function, first diff gives dx,
	# just need to peg on the zero to close out
	lx <- mx2lx(mx)
	-diff(c(lx,0))
}

mx2ex <- function(mx, trunc = 0){
	x <- 1:length(mx) - 1
	mx <- mx[x >= trunc]
	lx <- mx2lx(mx)
	lx <- c(lx,0)
	# spline interpolate to get Lx. force monotonic decline.
	# works best in log, then unlog :-)
	Lx <- exp(splinefun(log(lx)~c(x,x[length(x)+1]),method="monoH.FC")((x+.5)))
	# Tx are all the years lived above age x
	Tx <- rev(cumsum(rev(Lx)))
	# conditioned on survival, this give ex
	Tx / lx[-length(lx)]
}

# essentially same as excercise edagger, but we need it to start from mx
# this is why we needed all the above functions
mx2edagger <- function(mx, trunc = 0){

	# this is the new part
	dx      <- mx2dx(mx)
	ex      <- mx2ex(mx)

	age     <- 1:length(mx)-1
	ex 		<- ex[age >= trunc]
	dx 		<- dx[age >= trunc]
	dx 		<- stand(dx)
	age     <- age[age >= trunc]
	# assume constant hazard after 110 = same ex
	ex   	<- c(ex,ex[length(ex)])
	age2 	<- c(age,age[length(age)]+1)
	# linear interp would give the same thing.
	ex.5 	<- splinefun(ex~age2,method="monoH.FC")((age+.5))
	sum(ex.5 * dx)
}

library(DecompHoriuchi)
DecompContinuousOrig()

library(HMDHFDplus)
SWE <- readHMDweb("SWE","mltper_1x1",username=us,password=pw)
RUS <- readHMDweb("RUS","mltper_1x1",username=us,password=pw)

# OK, sort of experimental here to select years, but let's just compare
# for years of similar e(0). Could also compare within years...

e0swe <- SWE$ex[SWE$Age == 0]
e0rus <- RUS$ex[RUS$Age == 0]

yrswe <- sort(unique(SWE$Year))
yrrus <- sort(unique(RUS$Year))

# problem, there is basically no overlap!!!
plot(yrrus, e0swe[yrswe %in% yrrus], ylim = c(55,80))
lines(yrrus, e0rus)

# literally, we'd need to go back to close to 1900 to find
# a SWE e(0) close to Russia's

# so here's a decomp within year:
Year <- 2000
mxrus <- RUS$mx[RUS$Year == Year]
mxswe <- SWE$mx[SWE$Year == Year]


# compare edaggers
(edswe <- mx2edagger(mxswe))
(edrus <- mx2edagger(mxrus))
edrus - edswe # by this measure, deaths of Swedes result in 5 fewer years of life lost
# than those of Russians (year 2000)
# N  = time steps over which to integrate. Bigger isn't necessarily better...
Diff <- DecompContinuousOrig(mx2edagger, mxswe, mxrus, N = 20)

sum(Diff); edrus - edswe # pretty close. close enough says I.

barplot(Diff, main = "Age pattern to difference in edagger", sub =
"(where differences are assigned to ages in
which mortality rate differences occur)")