# cassandra.R - a program that computes the value of certain pension
#               guarantee structures.
#
# You must fill in the blanks for a bunch of constants, which are at the
# top of the file, which pertain to your country and your beliefs about
# the most sensible values for these constants.
# After that, when you run the program, you get a bunch of useful results.
#
# The program embeds a range of alternative investment strategies and
# guarantee structures. If you want to make changes to those, then you
# will need to delve deeper into the program. But if you merely want
# numerical values that pertain to your setting, it's enough to tweak
# the numerical values at the top and run.
#
# Ajay Shah
# http://www.mayin.org/ajayshah
# Version 0.1, Sat Jul  1 22:44:23 IST 2006

# ---------------------------------------------------------------------------
# Assumptions
# ---------------------------------------------------------------------------

# A. Country characteristics
     # 1. starting wage, in local currency units, expressed per working day
     wage.range <- seq(20,200,10)
     # 2. Wage growth anticipated over lifetime - a mix of GDP growth and the
     #    wage-experience profile. Expressed in percent in real terms per year.
     wage.growth.per.year <- 3
     # 3. Your view about the riskless rate in real terms (in percent per year)
     real.riskless <- 2.5
     # 4. Your view about the equity premium (in percentage points per year)
     equity.premium <- 7
     # 5. Your view about the DAILY volatility of the equity index (percent)
     sigma.M <- 1.3
     # 6. Your view about the price of an annuity deep in the future
     #    This is expressed as the lumpsum required at retirement date
     #    in order to buy an annuity of 1 LCU per day.
     annuity.price <- 3842
     # 7. Post-retirement income - LCU/day - below which it is poverty
     poverty.income <- 46

# B. Other assumptions which are not country-specific.
#    From here on, you are tweaking assumptions about the program which
#    are generic and don't particularly refer to any one country.
     # Number of working days per year
     days.per.year <- 250                    # working days only, of course
     # Working life (in years)
     workinglife.years <- 40
     # Pension contribution rate (fixed for life), in percent
     contribution.rate <- 8.33
     # How much compute resources you have to splurge :-)
     N.simulations <- 500           # 5000 is nice, 10000 is infinity.

# ---------------------------------------------------------------------------
# You shouldn't need to edit anything from here on.
# ---------------------------------------------------------------------------

# A function that converts a number like "15% per year" into a number
# like 1.000559 - in 1+r format - for the per-day returns.
annual.to.daily <- function(r) {
  (1 + r/100)^(1/days.per.year)
}

# The function simengine does the heavy lifting. It is given 3 facts:
# the mean and sigma of the pension portfolio, and a `cvector' object
# which is a column from the contributions matrix.
# It pumps through N.simulations outcomes for the terminal pension
# wealth. This generates a vector W of all the pension wealth outcomes.
simengine <- function(mean, sd, cvector) {
    W <- as.numeric(N.simulations)
    for (simulation in 1:N.simulations) {
      r <- mean + (Z[,simulation] * sd)
      piece1 <- c(1, cumprod(r))
      ## Not c(1, cumprod(rev(r))), as rev(r) and r are equally i.i.d. N() :)
      W[simulation] <- piece1 %*% cvector
    }
    W
}

# This function is given a long vector with N.simulations outcomes
# from the risk-neutral distribution of terminal pension wealth. It
# applies various kinds of guarantee structures, and reports the price
# of each of 'em.
facets.rn <- function(W, lastwage, sumcontributions) {
  ## Guarantee 1: we will top you up beyond poverty
  target <- poverty.wealth
  new <- W
  new[new < target] <- target
  g1.payoff <- sum(new-W)/N.simulations

  # Guarantee 2: we will top you up to 50% replacement rate
  target <- 0.5*lastwage*annuity.price # The 0.5 IS TWEAKABLE 
  new <- W
  new[new < target] <- target
  g2.payoff <- sum(new-W)/N.simulations

  # Guarantee 3: we will top you off to sumcontributions
  target <- sumcontributions
  new <- W
  new[new < target] <- target
  g3.payoff <- sum(new-W)/N.simulations

  discounting <- (1+(real.riskless/100))^workinglife.years
  c(g1.payoff, g2.payoff, g3.payoff)/discounting
}

# This function is given a long vector with N.simulations outcomes
# from the true distribution of terminal pension wealth. It reports useful
# things about this such as quantiles, Pr(below poverty), etc.
facets.true <- function(W, lastwage, sumcontributions) {
  pr.poverty <- sum(W<poverty.wealth)/N.simulations
  v <- quantile(W, probs=c(.25,.5,.75))/(lastwage*annuity.price)
  rr.median <- v[2]
  rr.IQR <- v[3]-v[1]

  # Guarantee 1: we will top you up beyond poverty
  target <- poverty.wealth
  new <- W
  new[new<target] <- target
  v <- quantile(new, probs=c(.25,.5,.75))/(lastwage*annuity.price)
  g1.pr.poverty <- sum(new<poverty.wealth)/N.simulations
  g1.median <- v[2]
  g1.IQR <- v[3]-v[1]

  # Guarantee 2: we will top you up to 50% replacement rate
  target <- 0.5*lastwage*annuity.price # TWEAKABLE 0.5 HERE
  new <- W
  new[new<target] <- target
  v <- quantile(new, probs=c(.25,.5,.75))/(lastwage*annuity.price)
  g2.pr.poverty <- sum(new<poverty.wealth)/N.simulations
  g2.median <- v[2]
  g2.IQR <- v[3]-v[1]

  # Guarantee 3: we will top you off to sumcontributions
  target <- sumcontributions
  new <- W
  new[new<target] <- target
  v <- quantile(new, probs=c(.25,.5,.75))/(lastwage*annuity.price)
  g3.pr.poverty <- sum(new<poverty.wealth)/N.simulations
  g3.median <- v[2]
  g3.IQR <- v[3]-v[1]

  toreturn <-        c(pr.poverty,    rr.median, rr.IQR,
                       g1.pr.poverty, g1.median, g1.IQR,
                       g2.pr.poverty, g2.median, g2.IQR,
                       g3.pr.poverty, g3.median, g3.IQR)
  names(toreturn) <- c("Pr(poverty)", "Median(rr)", "IQR(rr)",
                       "g1.Pr(Pov)",  "g1.Median", "g1.IQR",
                       "g2.Pr(Pov)",  "g2.Median", "g2.IQR",
                       "g3.Pr(Pov)",  "g3.Median", "g3.IQR")
  toreturn
}

# Compute a few variants of the above assumptions...
poverty.wealth <- poverty.income * annuity.price
workinglife.days <- workinglife.years * days.per.year # switch to life in days
wage.growth.per.day <- annual.to.daily(wage.growth.per.year)

# Prestored N(0,1) random numbers.
set.seed(1001)
Z <- matrix(rnorm(workinglife.days*N.simulations), nrow=workinglife.days)

# Make a matrix containing the vector of contributions for all values
# of initial wage. This matrix has one row per day of working life and
# one column per initial-wage scenario that we want to analyse.
time <- 1:workinglife.days
contributions <- matrix(NA, nrow=(workinglife.days+1), ncol=length(wage.range))
for (i in 1:length(wage.range)) {
    initial.wage <- wage.range[i]
    ts.wage <- initial.wage * wage.growth.per.day^time
    c.ts <- contribution.rate * ts.wage / 100
    contributions[,i] <- c(rev(c.ts), 1)
}
lastwage <- wage.range*((1+(wage.growth.per.year/100))^workinglife.years)
summation <- colSums(contributions)

# So now the matrix contributions is ready to roll. We will repeatedly
# use it in a large number of simulations. And, the vector `summation'
# stores the sum total of moneys that went in (zero real interest) under
# each income assumption.
#
# Aside. Two things are precomputed: Z and contributions. Of these,
# the real cost is of making Z. Compared with that, computing
# `contributions' is low.

# Okay, now we're ready to do some real work.
# For every level of equity exposure, we have two classes of
# questions.
# The first is to envision the outcomes in realistic terms - correct
# mean, correct sd.
# The second is to price derivatives, which requires shifting to a
# risk-neutral measure: correct sd, fake mean.

mean.rf <- annual.to.daily(real.riskless)
all <- NULL
for (equityfraction in seq(0,1,.05)) { # TWEAKABLE - YOU COULD USE SOME OTHER GRID
    for (relevant in 1:length(wage.range)) {
      cat("\n\nWorking for equity ", equityfraction, " and initial wage ", wage.range[relevant], "\n")

      # Assemble materials
      mean.true <- annual.to.daily(
                     equityfraction*(equity.premium+real.riskless) +
                     ((1-equityfraction)*(real.riskless))
                   )
      sd <- equityfraction*sigma.M/100
      tmp.c <- contributions[,relevant]
      tmp.s <- summation[relevant]
      tmp.w <- lastwage[relevant]

      # Compute
      W.true <- simengine(mean.true, sd, tmp.c)
      W.rn   <- simengine(mean.rf,   sd, tmp.c)

      # Summarise & print.
      m <- matrix(facets.true(W.true, tmp.w, tmp.s), ncol=4)
      m <- rbind(m, c(0, facets.rn(W.rn, tmp.w, tmp.s)))
      colnames(m) <- c("g: None", "g: Poverty", "g: 50%", "g: 0% real")
      rownames(m) <- c("Pr(Poverty)", "Median(rr)", "IQR(rr)", "Price")
      print(m)
      all <- rbind(all, c(equityfraction, wage.range[relevant], as.vector(m)))
   }
}

colnames(all) <- c("equityfraction", "initial.wage",
                   "g0.Pr", "g0.rr.median", "g0.rr.IQR", "g0.price",
                   "g1.Pr", "g1.rr.median", "g1.rr.IQR", "g1.price",
                   "g2.Pr", "g2.rr.median", "g2.rr.IQR", "g2.price",
                   "g3.Pr", "g3.rr.median", "g3.rr.IQR", "g3.price")
tail(all)
all <- as.data.frame(all)
save(all, file="results.rda")