R/maximumSpendingForMinimumRuinTimeV2.R

In eaoestergaard/UNPIE: UNPIE Understanding Pensions In Europe

#' Calculates scenarios of future value of annuity payments (fv) with stochastic returns
#'
#' @param wealth The wealth at retirement. Must be entered as a positive number
#' @param minumumRuinTime Minimum time to ruin.  Must be entered as a positive integer
#' @param mu The expected interest real return per period. Default is zero. Must be entered as decimal
#' @param sigma Volatility of expected interest real return per period. Default is zero. Must be entered as decimal
#' @param nScenarios The total number of scenarios to be made. Default is one scenario
#' @param prob Probability to exceed minimum time to ruin. Must be entered as decimal.
#' @param seed Integer vector, containing the random number generator (RNG) state for random number generation in R
#' @import timeSeries
#' @export
#' @examples
#' maximumSpendingForMinimumRuinTimeV2(wealth=14000,minumumRuinTime = 16,mu=0.03,sigma=0.08,nScenarios=200, prob = 0.9, seed =NULL)
#'
maximumSpendingForMinimumRuinTimeV2 <- function(wealth=14000,
                                              minumumRuinTime=16,
                                              mu=0.03,
                                              sigma=0.08,
                                              nScenarios=200,
                                              prob=0.9,
                                              seed=NULL) {
  ##Type check
  if(!is.scalar(wealth)) return(stop("wealth must be of type scalar",call. = FALSE))
  if(!is.scalar(minumumRuinTime)) return(stop("minumumRuinTime must be of type scalar",call. = FALSE))
  if(!is.scalar(mu)) return(stop("mu must either be of type scalar", call. = FALSE))
  if(!is.scalar(sigma)) return(stop("sigma must either be of type scalar", call. = FALSE))
  if(!is.scalar(nScenarios)) return(stop("nScenarios must be of type scalar",call. = FALSE))
  if(!is.scalar(prob)) return(stop("prob must be of type scalar with value in [0,1]",call. = FALSE))
  if(prob>1 || prob<0) return(stop("prob must be of type scalar with value in [0,1]",call. = FALSE))
  if(is.numeric(seed)){
    set.seed(seed)
  }

  t <- minumumRuinTime + 50
  nper=t
  ruin <- minumumRuinTime
  n <- nScenarios
  p <- prob

  #if(sigma==0){
  #  nScenarios <- 1
  #  n <- 1
  #}

  r <- matrix(rnorm(nScenarios*t,mu,sigma),nScenarios,t)
  rf = exp(timeSeries::rowCumsums(r))

  rc = exp(timeSeries::rowCumsums(r[,1:(nper-1)]))
  rc = cbind(rep(0,nScenarios),rc)
  rc =1+ t(apply(rc,1,function(x) cumsum(x)))


  # Calculate saving.a
  Up <- rf*wealth

  scenariosCalculate <- function(x){
    Up - rc*x
  }

  nrs <- function(scenarios){
    signchanges <- scenarios <= 0
    nr <- apply(signchanges,1,function(x) which(diff(sign(x))!=0)[1])
    nr[is.na(nr)] <- t
    return(nr)
  }

  # Normal case where sigma > 0 and prob < 1.
  objectivFunction.Normal <- function(x){
    scenarios <- scenariosCalculate(x)
    nr <- nrs(scenarios)
    pm <- sum(sign(nr >= ruin))/n
    return(pm-p)
  }

  # Sigma = 0, simpler goalseeking.
  objectivFunction.SigmaZero <- function(x){
    scenarios <- scenariosCalculate(x)
    return(scenarios[1,ruin])
  }

  # Probability = 1, since 0 payment will achive 1 we need
  # approc
  objectivFunction.ProbOne <- function(x){
    scenarios <- scenariosCalculate(x)
    nr <- nrs(scenarios)
    return(min(nr/ruin)-1)
  }

  bounds <- c(0,wealth/ruin*2)
  objfun <- objectivFunction.Normal
  if(sigma==0 || minumumRuinTime <= 1 ){
    objfun <- objectivFunction.SigmaZero
  }else if(p == 1){
    objfun <- objectivFunction.ProbOne
  }

  res <- uniroot(objfun,bounds,extendInt = "yes")
  sce <- cbind(wealth,scenariosCalculate(res$root))


  for (i in 1:length(sce[,1])){
    if(any(sce[i,] <= 0))
      sce[i,which(sce[i,] <= 0)[1]:length(sce[1,])] = 0
  }

  return(list(res = res,scenarios = sce, nr = nrs(sce)))
}