Most popular Environment to run R

http://www.rstudio.com/products/rstudio/

Free & Open-Source Integrated Development Environment (IDE) for R

Features:

1. Built-in R console
2. R, LaTeX, HTML, C++ syntax highlighting, R code completion
3. Easy management of multiple projects
4. Integrated R documentation
5. Interactive debugger
6. Package development tools

Note: R must be installed first!

Basic quantities in the analysis of mortality

• Survival function

\(s_{x} = Pr (X > x)\)

• Probability that (x) will survive for another t

\(p_{xt} = s_{x+t}/s_{x}\)

• Probability that (x) will die within t years

\(q_{xt} = [s_{x} - s_{x+t}]/s_{x}\)

• Mortality intensity (hazard function or force of mortality)

\(\mu_{x} = lim_{h \to 0} 1/h * q_{xh}\)

Probability that (x) will die within h

Some special parametric laws of mortality

• De Moivre

\(\mu_{x} = 1/ (ω – x)\) subject to \(0 \leq x < ω\)

• Gompertz

\(\mu_{x} = Bc^{x}\) subject to \(x \geq 0, B>0, c>1\)

• Makeham

\(\mu_{x} = A + Bc^{x}\) subject to \(x \geq 0, B>0, c>1, A>=-B\)

• Thiele

\(\mu_{x} = B_{1}C_{1}^{-x}+B_{2}C_{2}^{[-1/2(x-k)^2]}+B_{3}C_{3}^{x}\) subject to \(x\geq 0, B_{1}, B_{2}, B_{3}>0, C_{1}, C_{2}, C_{3}>1\)

Generalised Age-Period-Cohort Stochastic Mortality Models

• Lee-Carter (LC) \(\hat{\mu_{xt}} = \alpha_{x} + \beta^{(1)}_{x} \kappa^{(1)}_{t}\)
• Age-Period-Cohort (APC) \(\hat{\mu_{xt}} = \alpha_{x} + \kappa^{(1)}_{t} +\gamma_{t-x}\)
• Cairns-Blake_Dowd (CBD) \(\hat{\mu_{xt}} = \kappa^{(1)}_{t} + (x-\bar{x})\kappa^{(2)}_{t}\)
• Quadratic CBD with cohort effects M6 \(\hat{\mu_{xt}} = \kappa^{(1)}_{t} + (x-\bar{x})\kappa^{(2)}_{t} +\gamma_{t-x}\)
• Quadratic CBD with cohort effects M7 \(\hat{\mu_{xt}} = \kappa^{(1)}_{t} + (x-\bar{x})\kappa^{(2)}_{t} + ((x-\bar{x})^2-\hat{\sigma^{2}_{x}}) \kappa^{3}_{t}\)
• Quadratic CBD with cohort effects M8 \(\hat{\mu_{xt}} = \kappa^{(1)}_{t} + (x-\bar{x})\kappa^{(2)}_{t} + (x_{c}-x)\gamma_{t-x}\)
• Plat \(\hat{\mu_{t}} = \alpha_{x} + \kappa^{(1)}_{t} + (x-\bar{x})\kappa^{(2)}_{t} + \gamma_{t-x}\)

## LC model under a Binomial setting - M1
LC <- lc(link = "logit")
LCfit <- fit(LC, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## Initialising
## Running start-up iterations..
## Running main iterations...............
## Done
## StMoMo: Finish fitting with gnm

LCres <- residuals(LCfit)

## Collect the fit statistics
modelFit[1, 1] <- AIC(LCfit)
modelFit[1, 2] <- BIC(LCfit)

plot(LCfit, nCol = 3)

## Forecast
LCfor <- forecast(LCfit, h = forecastTime)
LCqxt <- cbind(LCfor$fitted, LCfor$rates)

## Fitting APC Currie (2006) - M3
APC <- apc("logit")
APCfit <- fit(APC, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
              ages.fit = ages.fit, wxt = weights, start.ax = LCfit$ax,
              start.bx = LCfit$bx, start.kt = LCfit$kx)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

## Collecting fit results
modelFit[2, 1] <- AIC(APCfit)
modelFit[2, 2] <- BIC(APCfit)

## Forecasting
APCfor <- forecast(APCfit, h = forecastTime, gc.order = c(1, 1, 0))
APCqxt <- cbind(APCfor$fitted, APCfor$rates)

CBD model Cairns (2009) under a Binomial distribution of deaths Haberman and Renshaw (2011) - M5

CBD <- cbd()
CBDfit <- fit(CBD, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
              ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

modelFit[3, 1] <- AIC(CBDfit)
modelFit[3, 2] <- BIC(CBDfit)

CBDfor <- forecast(CBDfit, h = forecastTime)
CBDqxt <- cbind(CBDfor$fitted, CBDfor$rates)

## M6
M6 <- m6()
M6fit <- fit(M6, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

modelFit[4, 1] <- AIC(M6fit)
modelFit[4, 2] <- BIC(M6fit)

M6for <- forecast(M6fit, h = forecastTime, gc.order = c(2, 0, 0))
M6qxt <- cbind(M6for$fitted, M6for$rates)

M7 <- m7(link = "logit")
M7fit <- fit(M7, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

modelFit[5, 1] <- AIC(M7fit)
modelFit[5, 2] <- BIC(M7fit)

M7for <- forecast(M7fit, h = forecastTime, gc.order = c(2, 0, 0))
M7qxt <- cbind(M7for$fitted, M7for$rates)

M8 <- m8(link = "logit", xc = 65)
M8fit <- fit(M8, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

modelFit[6, 1] <- AIC(M8fit)
modelFit[6, 2] <- BIC(M8fit)

M8for <- forecast(M8fit, h = forecastTime, gc.order = c(2, 0, 0))
M8qxt <- cbind(M8for$fitted, M8for$rates)

f2 <- function(x, ages) mean(ages) - x
constPlat <- function(ax, bx, kt, b0x, gc, wxt, ages){
  nYears <- dim(wxt)[2]
  x <- ages
  t <- 1:nYears
  c <- (1 - tail(ages, 1)):(nYears - ages[1])
  xbar <- mean(x)
  phiReg <- lm(gc ~ 1 + c + I(c^2), na.action = na.omit)
  phi <- coef(phiReg)
  gc <- gc - phi[1] - phi[2] * c - phi[3] * c^2
  kt[2, ] <- kt[2, ] + 2 * phi[3] * t
  kt[1, ] <- kt[1, ] + phi[2] * t + phi[3] * (t^2 - 2 * xbar * t)
  ax <- ax + phi[1] - phi[2] * x + phi[3] * x^2
  ci <- rowMeans(kt, na.rm = TRUE)
  ax <- ax + ci[1] + ci[2] * (xbar - x)
  kt[1, ] <- kt[1, ] - ci[1]
  kt[2, ] <- kt[2, ] - ci[2]
  list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)
}

PLAT <- StMoMo(link = "logit", staticAgeFun = TRUE,
               periodAgeFun = c("1", f2), cohortAgeFun = "1",
               constFun = constPlat)
PLATfit <- fit(PLAT, Dxt = Dxt, Ext= E0xt, ages = ages, years = years, 
               ages.fit = ages.fit, wxt = weights)

## StMoMo: The following cohorts have been zero weigthed: 1860 1861 1862 1982 1983 1984 
## StMoMo: Start fitting with gnm
## StMoMo: Finish fitting with gnm

modelFit[7, 1] <- AIC(PLATfit)
modelFit[7, 2] <- BIC(PLATfit)

PLATfor <- forecast(PLATfit, h = forecastTime, gc.order = c(2, 0, 0))
PLATqxt <- cbind(PLATfor$fitted, PLATfor$rates)

What model has the best fit?

• Limited data, often poor quality at high ages
• Kannisto extrapolation

extrapolate <- kannistoExtrapolation(qx=PLATqxt, ages=ages.fit, 
                                     years=years_chart, max_age=120, nObs=15)
PLATqxtExtr <- extrapolate$qxt

Curious about the function?

Just type kannistoExtrapolation in the console.

• Assumption: we have imaginary portofolio
• Premium read from input data
• Youngest entry age 25
• Retirement age 65
• Annuity is paid until age 120
• Annuity annual amount 1000 EUR
• Valuation year 2015

## Read in portfolio data
portfolio <- read.csv('portfolio.csv')

## Assumptions
ages.fit <- 25:120
valyear <- 2015
pension <- 1000

The expected or mean value of the present value of future cash flows for current obligations, projected over the contract’s run-off period, taking into account all up-to-date financial market and actuarial information.

• Calculate BEL for all models: LC, APC, CBD, M6, M7, M8 and PLAT.
• Use our function DFcashflowfunction

DFcashflow = function(qxt, ageStart, omegaAge, pensionAge, valyear, ir, type = 1)

• Function arguments:
  • qxt - mortality rates
  • ageStart - age of the insured
  • omegaAge - omega age
  • pensionAge - retirement age
  • valyear - valuation year
  • ir - interest rate value (for simplifaction constant)
  • type - indicator for discounting premium or annuity

BEL <- array(NA, c(7,1))
rownames(BEL) <- c("LC", "APC", "CBD", "M6", "M7", "M8", "PLAT")
colnames(BEL) <- "BEL"

for (m in 1:length(models)){
output <- list()
output2 <- list()
for (i in 1:nrow(portfolio)){
  output[[i]] <- DFcashflow(models[[m]]*expF, ageStart = portfolio$age[i], 
                            omegaAge = 120, pensionAge = 65, valyear = valyear, 
                            ir = 0.02, type = 1)*pension
  output2[[i]] <- DFcashflow(models[[m]]*expF, ageStart = portfolio$age[i], 
                             omegaAge = 120, pensionAge = 65, valyear = valyear, 
                             ir = 0.02, type = 2)*portfolio$Premium[i]
}
BEL[m, 1] <- round(do.call(sum, output)-do.call(sum, output2),2)
}

modelSim <- list()
selectBEL <- list()
for (m in 1:models2run){
  modelSim[[m]] <- simulate(modelsFitted[[m]], nsim = nsim, h = forecastTime)
  collectBEL <- array(NA, c(200,1))
  for (s in 1:nsim){
    prem_s <- 0
    ben_s <- 0
    qx <- cbind(modelSim[[m]]$fitted[, , s], modelSim[[m]]$rates[, , s])
    extrapolate <- kannistoExtrapolation(qx, ages.fit, years_chart)
    for (i in 1:nrow(portfolio)){
      ben_s <- ben_s+DFcashflow(extrapolate$qxt*expF, ageStart = portfolio$age[i], 
              omegaAge = 120, pensionAge = 65, valyear = valyear, 
              ir = 0.02, type = 1)*pension
      prem_s <- prem_s+DFcashflow(extrapolate$qxt*expF, ageStart = portfolio$age[i], 
                omegaAge = 120, pensionAge = 65, valyear = valyear, 
                ir = 0.02, type = 2)*portfolio$Premium[i]
    }
    collectBEL[s, 1] <- round(ben_s - prem_s, 2)
  }
  selectBEL[[m]] <- quantile(collectBEL, probs = 0.995, type = 1)
}
SCR <- round(as.numeric(selectBEL) - BEL, 2)
colnames(SCR) <- "SCR"

"StMoMO: An R Package for Stochastic Mortality Modelling" A.M. Villegas et al.

"Stochastic Modelling of Mortality Risks" Frankie Gregorkiewicz

"Solvency II Glossary" CEA and the Groupe Consultatif

"R: The most powerful and most widely used statistical software" Revolution Analytics

An independent and innovative consultancy firm specialized in the field of actuarial science and risk management

• Founded in 2006
• 80 FTE in 2015
• Headquarters in Amsterdam
• Local, dedicated office in Warsaw
• Highly qualified and trained actuaries and risk professionals

Comprehensive, hands-on services for insurance companies

• Actuarial modeling
• Solvency II and regulation
• Risk management function
• Actuarial function
• Governance
• Strategy and profitability
• Investment consulting and ALM