an_introduction_to_lifecontingencies_package.R
In lifecontingencies: Financial and Actuarial Mathematics for Life Contingencies

### R code from vignette source 'an_introduction_to_lifecontingencies_package.Rnw'

###################################################
### code chunk number 1: setup
###################################################
	options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
	set.seed(123)
	numSim=100


###################################################
### code chunk number 2: load
###################################################
library("lifecontingencies")


###################################################
### code chunk number 3: structure1
###################################################
showClass("lifetable")


###################################################
### code chunk number 4: structure2
###################################################
showClass("actuarialtable")


###################################################
### code chunk number 5: structure3
###################################################
showMethods(classes=c("actuarialtable","lifetable"))


###################################################
### code chunk number 6: ir1
###################################################
interest2Discount(0.03)
discount2Interest(interest2Discount(0.03))


###################################################
### code chunk number 7: ir1
###################################################
convertible2Effective(i=0.10,k=4)


###################################################
### code chunk number 8: npv1
###################################################
capitals <- c(-1000,200,500,700)
times <- c(0,1,2,5)
presentValue(cashFlows=capitals, timeIds=times, interestRates=0.03)


###################################################
### code chunk number 9: npv3
###################################################
presentValue(cashFlows=capitals, timeIds=times, 
interestRates=c( 0.04, 0.02, 0.03, 0.05), 
probabilities=c(1,1,1,0.5))


###################################################
### code chunk number 10: npv3
###################################################
getIrr <- function(p) (presentValue(cashFlows=capitals, timeIds=times, 
interestRates=p) - 0)^2
nlm(f=getIrr, p=0.1)$estimate


###################################################
### code chunk number 11: ann1
###################################################
100 * annuity(i=0.03, n=5)


###################################################
### code chunk number 12: ann2
###################################################
100 * accumulatedValue(i=0.03, n=5)


###################################################
### code chunk number 13: ann3
###################################################
ann1 <- annuity(i=0.03, n=5, k=1, type="immediate")
ann2 <- annuity(i=0.03, n=5, k=12, type="immediate")
c(ann1,ann2)


###################################################
### code chunk number 14: ann5
###################################################
incrAnn <- increasingAnnuity(i=0.03, n=10, type="due")
decrAnn <- decreasingAnnuity(i=0.03, n=10, type="immediate")
c(incrAnn, decrAnn)


###################################################
### code chunk number 15: ann6
###################################################
annuity(i=((1+0.04)/(1+0.03)-1), n=10)


###################################################
### code chunk number 16: capAmort1
###################################################
capital <- 100000
interest <- 0.05 
payments_per_year <- 2
rate_per_period <- (1+interest)^(1/payments_per_year)-1
years <- 30
R <- 1/payments_per_year * 
capital/annuity(i=interest, n=years, 
		k=payments_per_year)
R


###################################################
### code chunk number 17: capAmort2
###################################################
balanceDue <- numeric(years * payments_per_year)
balanceDue[1] <- capital * (1+rate_per_period) - R
for(i in 2:length(balanceDue)) balanceDue[i]<-
  balanceDue[i-1] * (1+rate_per_period) - R


###################################################
### code chunk number 18: figBalanceDue
###################################################
plot(x=c(1:length(balanceDue)), y=balanceDue, main="Loan amortization", 
		ylab="EoP balance due", xlab="year", type="l",col="steelblue")


###################################################
### code chunk number 19: BPFun1
###################################################
bond<-function(faceValue, couponRate, couponsPerYear, yield,maturity)
{
	out <- numeric(1)
	numberOfCF <- maturity * couponsPerYear
	CFs <- numeric(numberOfCF)
	payments <- couponRate * faceValue / couponsPerYear 
	cf <- payments * rep(1,numberOfCF)
	cf[numberOfCF] <- faceValue + payments 
	times <- seq.int(from=1/couponsPerYear, to=maturity, 
                by=maturity/numberOfCF)
	out <- presentValue(cashFlows=cf, interestRates=yield, 
                   timeIds=times)
	return(out)
}

perpetuity<-function(yield, immediate=TRUE)
{
	out <- numeric(1)
	out <- 1 / yield
	out <- ifelse(immediate==TRUE, out, out*(1+yield))
	return(out)
}



###################################################
### code chunk number 20: BPFun2
###################################################
bndEx1 <-bond(1000, 0.06, 2, 0.05, 3)
bndEx2 <-bond(1000, 0.06, 2, 0.06, 3)
ppTy1 <-perpetuity(0.1)
c(bndEx1, bndEx2, ppTy1)


###################################################
### code chunk number 21: durationAndConvexity
###################################################
cashFlows <- c(100,100,100,600,500,700)
timeVector <- seq(1:6)
interestRate <- 0.03
dur1 <-duration(cashFlows = cashFlows, timeIds = timeVector, 
		i = interestRate, k = 1, macaulay = TRUE)
dur2 <-duration(cashFlows = cashFlows, timeIds = timeVector, 
		i = interestRate, k = 1, macaulay = FALSE)
cvx1 <-convexity(cashFlows = cashFlows, timeIds = timeVector, 
		i = interestRate, k = 1)
c(dur1, dur2, cvx1)


###################################################
### code chunk number 22: alm1
###################################################
GTCFin<- 10000 * (1 + 0.05)^7
GTCFin


###################################################
### code chunk number 23: alm2
###################################################
yieldT0 <- 0.04
durLiab <- 7
pvLiab <- presentValue(cashFlows = GTCFin,timeIds = 7,
		interestRates = yieldT0)
convLiab <- convexity(cashFlows=GTCFin, timeIds = 7, 
		i=yieldT0)
pvBond <- bond(100,0.03,1,yieldT0,5)
durBond <- duration(cashFlows=c(3,3,3,3,103), 
		timeIds=seq(1,5), i = yieldT0)
convBond <- convexity(cashFlows=c(3,3,3,3,103), 
		timeIds=seq(1,5), i = yieldT0)
pvPpty <- perpetuity(yieldT0)
durPpty <- (1+yieldT0)/yieldT0
covnPpty <- 2/(yieldT0^2)


###################################################
### code chunk number 24: alm3
###################################################
a <- matrix(c(durBond, durPpty,1,1), nrow=2, 
		byrow=TRUE)
b <- as.vector(c(7,1))
weights <-solve(a,b)
weights


###################################################
### code chunk number 25: alm4
###################################################
bondNum <- weights[1] * pvLiab / pvBond
pptyNum <- weights[2] * pvLiab / pvPpty	
bondNum
pptyNum


###################################################
### code chunk number 26: alm5
###################################################
convAsset <- weights[1] * convBond + weights[2] * covnPpty
convAsset>convLiab


###################################################
### code chunk number 27: alm6
###################################################
yieldT1low <- 0.03
immunizationTestLow <- (bondNum * bond(100,0.03,1,yieldT1low,5) + 
			pptyNum * perpetuity(yieldT1low)> 
			GTCFin / (1+yieldT1low)^7)
yieldT1high <- 0.05
immunizationTestHigh <- (bondNum * bond(100,0.03,1,yieldT1high,5) + 
			pptyNum * perpetuity(yieldT1high)>
			GTCFin/(1+yieldT1high)^7)
immunizationTestLow
immunizationTestHigh


###################################################
### code chunk number 28: createALifecontingenciesObject
###################################################
x_example <- seq(from=0,to=9, by=1)
lx_example <- c(1000,950,850,700,680,600,550,400,200,50)
exampleLt <- new("lifetable", x=x_example, lx=lx_example, 
		name="example lifetable")


###################################################
### code chunk number 29: printShow
###################################################
print(exampleLt)


###################################################
### code chunk number 30: headAndTail
###################################################
head(exampleLt)


###################################################
### code chunk number 31: fromDataFrame1
###################################################
data("demoUsa")
data("demoIta") 
usaMale07 <- demoUsa[,c("age", "USSS2007M")]
usaMale00 <- demoUsa[,c("age", "USSS2000M")]
names(usaMale07) <- c("x","lx")
names(usaMale00) <- c("x","lx")
usaMale07Lt <-as(usaMale07,"lifetable")
usaMale07Lt@name <- "USA MALES 2007"
usaMale00Lt <-as(usaMale00,"lifetable")
usaMale00Lt@name <- "USA MALES 2000"


###################################################
### code chunk number 32: fromDataFrame2
###################################################
lxIPS55M <- with(demoIta, IPS55M)
pos2Remove <- which(lxIPS55M %in% c(0,NA))
lxIPS55M <-lxIPS55M[-pos2Remove]
xIPS55M <-seq(0,length(lxIPS55M)-1,1)
ips55M <- new("lifetable",x=xIPS55M, lx=lxIPS55M, 
		name="IPS 55 Males")
lxIPS55F <- with(demoIta, IPS55F)
pos2Remove <- which(lxIPS55F %in% c(0,NA))
lxIPS55F <- lxIPS55F[-pos2Remove]
xIPS55F <- seq(0,length(lxIPS55F)-1,1)
ips55F <- new("lifetable",x=xIPS55F, lx=lxIPS55F, 
		name="IPS 55 Females")


###################################################
### code chunk number 33: createFromSurvivalRates
###################################################
data("demoIta")
itaM2002 <- demoIta[,c("X","SIM92")]
names(itaM2002) <- c("x","lx")
itaM2002Lt <- as(itaM2002,"lifetable")
itaM2002Lt@name <- "IT 2002 Males"
itaM2002 <- as(itaM2002Lt,"data.frame")
itaM2002$qx <- 1-itaM2002$px
for(i in 20:60) itaM2002$qx[itaM2002$x==i] = 0.2 * itaM2002$qx[itaM2002$x==i]
itaM2002reduced <- probs2lifetable(probs=itaM2002[,"qx"], radix=100000,
		type="qx",name="IT 2002 Males reduced")


###################################################
### code chunk number 34: createAnActuarialtableObject
###################################################
exampleAct <- new("actuarialtable",x=exampleLt@x, lx=exampleLt@lx, 
interest=0.03, name="example actuarialtable")


###################################################
### code chunk number 35: methods1
###################################################
getOmega(exampleAct)


###################################################
### code chunk number 36: methods2
###################################################
print(exampleLt)
print(exampleAct)


###################################################
### code chunk number 37: methods3
###################################################
exampleActDf <- as(exampleAct, "data.frame")


###################################################
### code chunk number 38: toMarkovChain
###################################################
data(soa08)
require(markovchain)
soa08Mc<-as(soa08,"markovchainList")


###################################################
### code chunk number 39: figSurvivalFunction0
###################################################
data("soa08Act")
soa08ActDf <- as(soa08Act, "data.frame")


###################################################
### code chunk number 40: figSurvivalFunction
###################################################
plot(soa08Act, type="l",col="steelblue")


###################################################
### code chunk number 41: probabilityAndDemographics
###################################################
demoEx1<-pxt(ips55M,20,1)
demoEx2<-qxt(ips55M,30,2) 
demoEx3<-exn(ips55M, 50,20,"complete") 
c(demoEx1,demoEx2,demoEx3)


###################################################
### code chunk number 42: mxAndqx
###################################################
mx20t1 <- mxt(ips55M,20,1)
qx20t1 <- mx2qx(mx20t1)
c(mx20t1,qx20t1)


###################################################
### code chunk number 43: fractionalAges
###################################################
data("soa08Act")
pxtLin <- pxt(soa08Act,80,0.5,"linear") 
pxtCnst <- pxt(soa08Act,80,0.5,"constant force") 
pxtHyph <- pxt(soa08Act,80,0.5,"hyperbolic") 
c(pxtLin,pxtCnst,pxtHyph)


###################################################
### code chunk number 44: moreThanOneHead
###################################################
tablesList <- list(ips55M, ips55F)
jsp <- pxyzt(tablesList, x=c(65,63), t=2)
lsp <- pxyzt(tablesList, x=c(65,63), t=2, status="last") 
jelt <- exyzt(tablesList, x=c(65,63), status="joint") 
c(jsp,lsp,jelt)


###################################################
### code chunk number 45: lifeIns1
###################################################
data(soa08Act)
UComm <- Axn(actuarialtable=soa08Act, x=25, n=65-25, k=12)
UCpt <- ((soa08ActDf$Mx[26]-soa08ActDf$Mx[66])/soa08ActDf$Dx[26]) * 
		0.06/real2Nominal(i=0.06,k=12)
c(UComm, UCpt)


###################################################
### code chunk number 46: lifeIns2
###################################################
P <- UCpt/axn(actuarialtable=soa08Act,x=25,n=10)
P


###################################################
### code chunk number 47: lifeIns3
###################################################
(10 + 1 ) * Axn(actuarialtable=soa08Act, x=25, n=10) 
DAxn(actuarialtable = soa08Act, x=25, n=10) + 
IAxn(actuarialtable = soa08Act, x=25, n=10)


###################################################
### code chunk number 48: annuity1
###################################################
UCpt <- axn(actuarialtable=soa08Act, x=75, m=10)
UComm <- with(soa08ActDf,Nx[86]/Dx[76])
c(UCpt,UComm)


###################################################
### code chunk number 49: annuity2
###################################################
P=axn(actuarialtable=soa08Act, x=75, m=10) / 
		axn(actuarialtable=soa08Act, x=75, n=5)
P

PComm <- with(soa08ActDf,(Nx[86]/Dx[76]) / 
				((Nx[76]-Nx[81])/Dx[76]))
PComm


###################################################
### code chunk number 50: annuity3
###################################################
U <- axn(actuarialtable=soa08Act, x=75, m=10, k=12)
P <- axn(actuarialtable=soa08Act, x=75, m=10, k=12) / 
		axn(actuarialtable=soa08Act, x=75, n=5)
c(U,P)


###################################################
### code chunk number 51: lifeInsuranceBenefitReserve
###################################################
P=100000 * Axn(soa08Act,x=25,n=40)/axn(soa08Act,x=25,n=40)
reserveFun = function(t) return(100000*Axn(soa08Act,x=25+t,n=40-t)-P*
					axn(soa08Act,x=25+t,n=40-t))
for(t in 0:40) {if(t%%5==0) cat("At time ",t,
				" benefit reserve is ", 
				reserveFun(t),"\n")}


###################################################
### code chunk number 52: annuityReserve
###################################################
yearlyRate <- 12000
irate <- 0.02
APV <- yearlyRate*axn(soa08Act, x=25, i=irate,m=65-25,k=12)
levelPremium <- APV/axn(soa08Act, x=25,n=65-25,k=12)

annuityReserve<-function(t) {
	out<-NULL
	if(t<65-25) out <- yearlyRate*axn(soa08Act, x=25+t, 
    i=irate, m=65-(25+t),k=12)-levelPremium*axn(soa08Act, 
              x=25+t, n=65-(25+t),k=12) else {
		out <- yearlyRate*axn(soa08Act, x=25+t, i=irate,k=12)
	}
	return(out)
}

years <- seq(from=0, to=getOmega(soa08Act)-25-1,by=1)
annuityRes <- numeric(length(years))
for(i in years) annuityRes[i+1] <- annuityReserve(i)
dataAnnuityRes <- data.frame(years=years, reserve=annuityRes)


###################################################
### code chunk number 53: annuityReserveGraph
###################################################
plot(y=dataAnnuityRes$reserve, x=dataAnnuityRes$years,
col="steelblue", main="Deferred annuity benefit reserve",
ylab="amount",xlab="years",type="l")


###################################################
### code chunk number 54: expAugmented
###################################################
G <- (100000*Axn(soa08Act, x=35) + (2.5*100000/1000 + 25)*
			axn(soa08Act,x=35))/((1-.1)*axn(soa08Act,x=35))
G


###################################################
### code chunk number 55: twoHeadsAnnuitImmediate
###################################################
twoLifeTables <- list(maleTable=soa08Act, femaleTable=soa08Act)
axn(soa08Act, x=65,m=1)+axn(soa08Act, x=70,m=1)-
axyzn(tablesList=twoLifeTables,	x=c(65,y=70),status="joint",m=1) 
axyzn(tablesList=twoLifeTables, x=c(65,y=70), status="last",m=1)


###################################################
### code chunk number 56: revesionaryAnuity
###################################################
axn(actuarialtable = soa08Act, x=60,m=1)-
		axyzn(tablesList = twoLifeTables, 
				x=c(65,60),status="joint",m=1)


###################################################
### code chunk number 57: rLife1
###################################################
rLife(n = 5, object = soa08Act, x = 45, type = "Kx")


###################################################
### code chunk number 58: rLife2
###################################################
futureLifetimes <- as.data.frame(rLifexyz(n=numSim, 
				tablesList=list(husband=ips55M,wife=ips55F),
				x=c(68,65), type="Tx"))
names(futureLifetimes) <- c("husband","wife")
temp <- futureLifetimes$wife - futureLifetimes$husband
futureLifetimes$widowance  <-  sapply(temp, max,0)
mean(futureLifetimes$widowance)


###################################################
### code chunk number 59: widowanceFig
###################################################
hist(futureLifetimes$widowance, freq=FALSE, main="Distribution of widowance yars",
			xlab="Widowance years", col="steelblue", nclass=100);abline(v=mean(futureLifetimes$widowance), 
			col="red", lwd=2)


###################################################
### code chunk number 60: AxnAPVAndStochastic
###################################################
APVAxn <- Axn(soa08Act,x=25,n=40,type="EV")
APVAxn
sampleAxn <- rLifeContingencies(n=numSim, lifecontingency="Axn",
		object=soa08Act,x=25,t=40,parallel=FALSE)
tt1 <-t.test(x=sampleAxn,mu=APVAxn)$p.value

APVIAxn <- IAxn(soa08Act,x=25,n=40,type="EV")
APVIAxn
sampleIAxn <- rLifeContingencies(n=numSim, lifecontingency="IAxn",
		object=soa08Act,x=25,t=40,parallel=FALSE)
tt2 <-t.test(x=sampleIAxn,mu=APVIAxn)$p.value

APVaxn <- axn(soa08Act,x=25,n=40,type="EV")
APVaxn
sampleaxn <- rLifeContingencies(n=numSim, lifecontingency="axn",
		object=soa08Act,x=25,t=40,parallel=FALSE)
tt3 <- t.test(x=sampleaxn,mu=APVaxn)$p.value

APVAExn <- AExn(soa08Act,x=25,n=40,type="EV")
APVAExn
sampleAExn <- rLifeContingencies(n=numSim, lifecontingency="AExn",
		object=soa08Act,x=25,t=40,parallel=FALSE)
tt4 <-t.test(x=sampleAExn,mu=APVAExn)$p.value
c(tt1, tt2,tt3, tt4)


###################################################
### code chunk number 61: figsim
###################################################
	par(mfrow=c(2,2))
	hist(sampleAxn, main="Term Insurance", xlab="Actuarial present value",nclass=50, col="steelblue",freq=FALSE);abline(v=APVAxn, col="red", lwd=2)
	hist(sampleIAxn, main="Increasing Life Insurance", xlab="Actuarial present value",nclass=50, col="steelblue",freq=FALSE);abline(v=APVIAxn, col="red", lwd=2)
	hist(sampleaxn, main="Temporary Annuity Due", xlab="Actuarial present value",nclass=50, col="steelblue",freq=FALSE);abline(v=APVaxn, col="red", lwd=2)
	hist(sampleAExn, main="Endowment Insurance", xlab="Actuarial present value",nclass=50, col="steelblue",freq=FALSE);abline(v=APVAExn, col="red", lwd=2)


###################################################
### code chunk number 62: randomMultipleLifeCon
###################################################
tablesList=list(soa08Act,soa08Act);x=c(60,60);m=0;status="last";t=30;k=1
APVAxyz<-Axyzn(tablesList=tablesList,x=x,n=t,status=status,type="EV")
samplesAxyz<-rLifeContingenciesXyz(n=numSim,lifecontingency = "Axyz",
		tablesList = tablesList,x=x,t=t,m=m,k=k,status=status,
		parallel=FALSE)
tt5<-t.test(x=samplesAxyz, mu=APVAxyz)$p.value
APVaxyz<-axyzn(tablesList=tablesList,x=x,n=t,m=m,k=k,status=status,type="EV")
samplesaxyz<-rLifeContingenciesXyz(n=numSim,lifecontingency = "axyz",
		tablesList = tablesList,x=x,t=t,m=m,k=k,status=status,
		parallel=FALSE)
tt6<-t.test(x=samplesaxyz, mu=APVaxyz)$p.value
c(tt5,tt6)


###################################################
### code chunk number 63: variance
###################################################
var(sampleAxn)
Axn(soa08Act, x=25,n=40, power=2)-Axn(soa08Act, x=25,n=40, power=1)^2


###################################################
### code chunk number 64: benefitPremium1
###################################################
APV <- Axn(actuarialtable = soa08Act, x=25, n=40)
APV


###################################################
### code chunk number 65: benefitPremium2
###################################################
samples <- rLifeContingencies(n=numSim, lifecontingency = "Axn", 
		object= soa08Act, x=25,t=40,parallel=FALSE)
pct90Pr <- as.numeric(quantile(samples,.90))
pct90Pr


###################################################
### code chunk number 66: benefitPremium3
###################################################
pct90Pr2 <- qnorm(p=0.90,mean=APV, sd=sd(samples)/sqrt(1000))
pct90Pr2


###################################################
### code chunk number 67: stochasticExampleFull1
###################################################
nsim <- 50
employees <- 100
salaryDistribution <- rlnorm(n=employees,m=10.77668944,s=0.086177696)
ageDistribution <- round(runif(n=employees,min=25, max=65))
policyLength <- sapply(65-ageDistribution, min, 1)

getEmployeeBenefit<-function(index,type="EV") {
	out <- numeric(1)
	out <- salaryDistribution[index]*Axn(actuarialtable=soa08Act, 
			x=ageDistribution[index],n=policyLength[index], 
			i=0.02,m=0,k=1, type=type)
	return(out)
}

require(parallel)
cl <- makeCluster(1) 
worker.init <- function(packages) {
	for (p in packages) {
		library(p, character.only=TRUE)
	}
	invisible(NULL)
}
clusterCall(cl, 
		worker.init, c('lifecontingencies'))
clusterExport(cl, varlist=c("employees","getEmployeeBenefit",
				"salaryDistribution","policyLength",
				"ageDistribution","soa08Act"))


###################################################
### code chunk number 68: stochasticExampleFull2
###################################################
employeeBenefits <- numeric(employees)
employeeBenefits <- parSapply(cl, 1:employees,getEmployeeBenefit, type="EV")
employeeBenefit <- sum(employeeBenefits)



benefitDistribution<-numeric(nsim)
yearlyBenefitSimulate<-function(i)
{
	out <- numeric(1)
	expenseSimulation <- numeric(employees)
	expenseSimulation <- sapply(1:employees, getEmployeeBenefit, type="ST")
	out <- sum(expenseSimulation)
	return(out)
}

benefitDistribution <- parSapply(cl, 1:nsim,yearlyBenefitSimulate )
stopCluster(cl)

riskMargin <- as.numeric(quantile(benefitDistribution,.75)-employeeBenefit)
totalBookedCost <- employeeBenefit+riskMargin

employeeBenefit
riskMargin
totalBookedCost


###################################################
### code chunk number 69: mdt1
###################################################
valdezDf<-data.frame(
		x=c(50:54),
		lx=c(4832555,4821937,4810206,4797185,4782737),
		heart=c(5168, 5363, 5618, 5929, 6277),
		accidents=c(1157, 1206, 1443, 1679,2152),
		other=c(4293,5162,5960,6840,7631)
)
valdezMdt<-new("mdt",name="ValdezExample",table=valdezDf)


###################################################
### code chunk number 70: md3a (eval = FALSE)
###################################################
## print(valdezMdt)


###################################################
### code chunk number 71: md3b
###################################################
valdezDf<-as(valdezMdt,"data.frame")
require(markovchain)
valdezMarkovChainList<-as(valdezMdt,"markovchainList")


###################################################
### code chunk number 72: mdt4
###################################################
getOmega(valdezMdt)
getDecrements(valdezMdt)


###################################################
### code chunk number 73: summary.mdt
###################################################
summary(valdezMdt)


###################################################
### code chunk number 74: mdt.dx1
###################################################
dxt(valdezMdt,x=51,decrement="other")
dxt(valdezMdt,x=51,t=2, decrement="other")
dxt(valdezMdt,x=51)


###################################################
### code chunk number 75: mdt.dx2
###################################################
dxt(valdezMdt,x=51,t=2, decrement="other")
pxt(valdezMdt,x=50,t=3)
qxt(valdezMdt,x=53,t=2,decrement=1)


###################################################
### code chunk number 76: mdt.randomSamples
###################################################
rmdt(n = 2,object = valdezMdt,x = 50,t = 2)


###################################################
### code chunk number 77: mdt.udd1
###################################################
qxt.prime.fromMdt(object = valdezMdt,x=53, decrement="accidents")


###################################################
### code chunk number 78: mdt.udd2
###################################################
qxt.fromQxprime(qx.prime = 0.01,other.qx.prime = c(0.03,0.06))


###################################################
### code chunk number 79: mdt.act1
###################################################
myTable<-data.frame(x=c(16,17,18),
  lx=c(20000,17600,14520),
  da=c(1300,1870,2380),
  doc=c(1100,1210,1331)
)
myMdt<-new("mdt",table=myTable,name="Sample")


###################################################
### code chunk number 80: mdt.act2
###################################################
Axn.mdt(object=myMdt,x=16,i=.1,decrement="da")


###################################################
### code chunk number 81: deadifalco.1a
###################################################
axnmdt.firsttype<-function (object, x, n, i , payment="advance", delta=0) {
  #delta is the annuity indexing
    out <- numeric(1)
    if (!(class(object) %in% c("lifetable", "actuarialtable", "mdt"))) 
      stop("Error! Only lifetable, actuarialtable or mdt classes are accepted")
    if (missing(object)) 
      stop("Error! Need a Multiple decrement table")
    if (missing(x)) 
      stop("Error! Need age!")
    if (x > getOmega(object)) {
      stop("Age greater than Omega")
    }
    if(class(object)=="mdt"){
    if (x < min(object@table$x)) {
      stop("Age lower than minimum age")
    }}
    if(class(object)=="actuarialtable"){
      if (x < min(object@x)) {
        stop("Age lower than minimum age")
      }}
    if(!(missing(i))){
      interest <- i
    }else{
      if(class(object)=="actuarialtable"){
        interest=object@interest
      }else{
        stop("Needed Interest Rate ")
      }
    }
    if (missing(n)) 
      n <- (getOmega(object)  - x)
    if (n == 0) {
      stop("Contract duration equal to zero")
    }
    probs = numeric(n)
    times = seq(from = 0, to = n-1, by = 1)
    if (payment == "arrears") times = times + 1
    for (j in 1:length(times)) probs[j] = pxt(object, x, times[j])
    out <- sum(apply(cbind(probs,((1 + interest)/(1+delta))^-times),1,prod))
  return(out)
}


###################################################
### code chunk number 82: deadifalco.1b
###################################################
data("de_angelis_di_falco")
HealthyMaleTable2013 <- de_angelis_di_falco$HealthyMaleTable2013
DAT<-new("actuarialtable", x=de_angelis_di_falco$DisabledMaleLifeTable$age, lx=de_angelis_di_falco$DisabledMaleLifeTable$'2013',
         name="DisabledTable",i=0.03)
axnmdt.firsttype(DAT,x=65,n=10,i=0.03,payment="arrears",delta=0.02)
axnmdt.firsttype(DAT,65,10,payment="arrears",delta=0.02)
axnmdt.firsttype(DAT,65,10,payment="arrears",i=0.03,delta=0.02)
#Last case equal to axn
axnmdt.firsttype(DAT,65,10,payment="arrears",delta=0)
axn(DAT,65,10,payment="arrears")