🔷 4.1 Introduction

4.1 Introduction

This chapter presents the Monte Carlo simulation results of the hybrid pension participation model. The objective is to evaluate long-term retirement accumulation and pension benefit outcomes under stochastic contribution behaviour, actuarial mortality assumptions, and compound investment growth.

The simulation models individuals entering the pension system at age 18, contributing irregularly throughout their working life, and retiring at age 64. Mortality assumptions are derived from the Standard Ultimate Life Table, and an annual effective investment return of 8% is assumed.

4.2 Model Parameters

# Demographic settings
age_start  <- 18
age_retire <- 64
age_max    <- 100

# Financial assumption
i_annual <- 0.08
i_m <- (1 + i_annual)^(1/12) - 1
v_m <- 1 / (1 + i_m)

# Simulation size
N <- 10000

i_m
## [1] 0.00643403

4.3 Mortality Assumptions

Annual mortality probabilities are obtained from the Standard Ultimate Life Table. Since the table begins at age 20, ages 18 and 19 are approximated using the mortality rate at age 20.

Annual mortality rates are converted to monthly probabilities using:

\[ q_{x,m} = 1 - (1 - q_x)^{1/12} \]

mort <- data.frame(
  age = 20:100,
  qx = c(
    0.000250,0.000253,0.000257,0.000262,0.000267,0.000273,
    0.000280,0.000287,0.000296,0.000305,0.000315,0.000327,
    0.000341,0.000356,0.000372,0.000391,0.000412,0.000436,
    0.000463,0.000493,0.000527,0.000565,0.000608,0.000656,
    0.000710,0.000771,0.000839,0.000916,0.001003,0.001100,
    0.001209,0.001331,0.001469,0.001623,0.001797,0.001993,
    0.002212,0.002459,0.002736,0.003048,0.003398,0.003792,
    0.004234,0.004730,0.005288,0.005915,0.006619,0.007409,
    0.008297,0.009294,0.010413,0.011670,0.013081,0.014664,
    0.016440,0.018433,0.020668,0.023175,0.025984,0.029132,
    0.032658,0.036607,0.041025,0.045968,0.051493,0.057665,
    0.064554,0.072237,0.080798,0.090326,0.100917,0.112675,
    0.125708,0.140128,0.156052,0.173599,0.192887,0.214030,
    0.237134,0.262294,0.289584
  )
)

qx_lookup <- function(age){
  a <- floor(age)
  if (a %in% c(18,19)) a <- 20
  if (a > 100) a <- 100
  mort$qx[mort$age == a]
}

q_monthly <- function(age){
  qA <- qx_lookup(age)
  1 - (1 - qA)^(1/12)
}

# Check pattern
q_monthly(20)
## [1] 2.083572e-05
q_monthly(60)
## [1] 0.0002836086
q_monthly(100)
## [1] 0.02808998

4.4 Contribution Behaviour

Contribution activity is stochastic. For each age group, the probability of contributing in a given month is defined as follows:

pi_by_age <- function(age){
  if (age >= 18 && age <= 25) return(0.65)
  if (age >= 26 && age <= 35) return(0.75)
  if (age >= 36 && age <= 45) return(0.80)
  if (age >= 46 && age <= 55) return(0.78)
  if (age >= 56 && age <= 63) return(0.70)
  return(0)
}

# Visualise contribution probability by age
ages <- 18:63
probs <- sapply(ages, pi_by_age)

plot(ages, probs, type="l",
     main="Monthly Contribution Probability by Age",
     xlab="Age", ylab="Probability")

4.5 Contribution Amount Distribution

Contribution amounts are assumed to follow a Lognormal distribution.

mean_contrib <- 3000
sd_contrib   <- 2000

var_y <- sd_contrib^2
sigma2 <- log(1 + var_y/(mean_contrib^2))
mu <- log(mean_contrib) - 0.5*sigma2
sigma <- sqrt(sigma2)

# Plot distribution
sim_amounts <- rlnorm(10000, mu, sigma)

hist(sim_amounts, breaks=50,
     main="Simulated Monthly Contribution Distribution",
     xlab="Contribution (KSh)")

4.6 Simulation Procedure

simulate_one <- function(id){
  
  balance <- 0
  age <- age_start
  month_index <- 0
  
  total_contrib <- 0
  months_contributed <- 0
  
  retired_balance <- NA
  reached_retirement <- FALSE
  
  while(age < age_max){
    
    # Mortality
    if(runif(1) < q_monthly(age)){
      break
    }
    
    # Retirement
    if(age >= age_retire){
      retired_balance <- balance
      reached_retirement <- TRUE
      break
    }
    
    # Contribution decision
    contrib <- 0
    
    if(runif(1) < pi_by_age(age)){
      contrib <- rlnorm(1, mu, sigma)
      total_contrib <- total_contrib + contrib
      months_contributed <- months_contributed + 1
    }
    
    balance <- (balance + contrib) * (1 + i_m)
    
    month_index <- month_index + 1
    age <- age_start + month_index/12
  }
  
  data.frame(
    id = id,
    total_contributions = total_contrib,
    months_contributed = months_contributed,
    retirement_balance = retired_balance,
    reached_retirement = reached_retirement
  )
}
results_df <- do.call(rbind, lapply(1:N, simulate_one))

retirees <- subset(results_df, reached_retirement == TRUE)

length(retirees$retirement_balance)
## [1] 9478
annuity_due <- function(){
  surv <- 1
  a <- 0
  age <- age_retire
  
  for(k in 0:((age_max - age_retire)*12)){
    a <- a + (v_m^k) * surv
    surv <- surv * (1 - q_monthly(age))
    age <- age + 1/12
    if(surv < 1e-12) break
  }
  a
}
a64 <- annuity_due()

retirees$annual_pension <- 12 * (retirees$retirement_balance / a64)
set.seed(100)

sample_table <- retirees[
  sample(1:nrow(retirees), 10),
  c("id",
    "months_contributed",
    "total_contributions",
    "retirement_balance",
    "annual_pension")
]

knitr::kable(sample_table, digits = 2,
             caption = "Table 4.X: Sample of Simulated Individual Outcomes")
Table 4.X: Sample of Simulated Individual Outcomes
id months_contributed total_contributions retirement_balance annual_pension
4003 4003 373 1135019 10113344 980828.1
530 530 383 1208650 11622211 1127163.3
3630 3630 363 1041417 9373629 909087.9
3909 3909 357 1077269 8648986 838809.4
4322 4322 356 984510 9634130 934352.1
8312 8312 349 1027075 9642916 935204.3
3233 3233 365 1070283 10216210 990804.3
8632 8632 352 1024324 9778932 948395.6
9131 9131 363 1058448 10290297 997989.6
2095 2095 369 1114897 10147170 984108.6

`

# --- Summary Tables for Chapter 4 ---

# Vector of retirement balances (only those who reached retirement)
ret_bal <- retirees$retirement_balance

# Table 4.1: Survival to retirement
table_survival <- data.frame(
  Total_Simulated = N,
  Reached_Retirement = length(ret_bal),
  Died_Before_Retirement = N - length(ret_bal),
  Survival_Rate = length(ret_bal) / N
)

# helper for stats
stats_table <- function(x){
  data.frame(
    Statistic = c("N", "Mean", "Median", "Std Dev", "Min", "P10", "P25", "P75", "P90", "Max"),
    Value = c(
      length(x),
      mean(x),
      median(x),
      sd(x),
      min(x),
      as.numeric(quantile(x, 0.10)),
      as.numeric(quantile(x, 0.25)),
      as.numeric(quantile(x, 0.75)),
      as.numeric(quantile(x, 0.90)),
      max(x)
    )
  )
}

# Table 4.2: Retirement balances
table_retbal <- stats_table(ret_bal)

# Table 4.3: Annual pension
table_pension <- stats_table(retirees$annual_pension)

knitr::kable(table_survival, digits = 4, caption = "Table 4.1: Survival to Retirement")
Table 4.1: Survival to Retirement
Total_Simulated Reached_Retirement Died_Before_Retirement Survival_Rate
10000 9478 522 0.9478 
knitr::kable(table_retbal, digits = 2, caption = "Table 4.2: Summary Statistics of Retirement Balances (KSh)")
Table 4.2: Summary Statistics of Retirement Balances (KSh)
Statistic Value
N 9478.0
Mean 10335467.2
Median 10317399.3
Std Dev 639148.7
Min 8106685.4
P10 9524447.1
P25 9890345.2
P75 10764616.1
P90 11160169.1
Max 13137636.9 
knitr::kable(table_pension, digits = 2, caption = "Table 4.3: Summary Statistics of Annual Pension (KSh)")
Table 4.3: Summary Statistics of Annual Pension (KSh)
Statistic Value
N 9478.00
Mean 1002370.34
Median 1000618.05
Std Dev 61986.91
Min 786215.16
P10 923714.73
P25 959200.83
P75 1043990.72
P90 1082352.86
Max 1274134.71 
hist(ret_bal, breaks = 50,
     main = "Distribution of Retirement Balances",
     xlab = "Retirement Balance (KSh)")

hist(retirees$annual_pension, breaks = 50,
     main = "Distribution of Annual Pension",
     xlab = "Annual Pension (KSh/year)")

4.10 Interpretation of Results

The simulation results demonstrate substantial dispersion in retirement outcomes due to stochastic contribution behaviour and investment compounding.

Key observations include:

These findings confirm that irregular participation significantly affects long-term pension adequacy.