R/present_values.R

In o1i/UWS: Umwandlungssatz in der zweiten Säule

#' Discount rate matrix
#'
#' Creates the matrix of discount rates based on an interest rate curve and
#' a mortality table.
#'
#' @param qx Matrix with probabilities of death to be used.
#' @param i Named numeric vector with the interest rates to be used in given
#' years.
#' @return Matrix with the Discount rates at the right places.
make_v <- function(qx, i){
  v <- 1 / (1 + i)
  v_mat <- matrix(v[as.character(outer(as.numeric(rownames(qx)),
                                       as.numeric(colnames(qx)), `+`))],
                  nrow=nrow(qx), ncol = ncol(qx))
  return(v_mat)
}

#' Annuity
#'
#' Present value of annuity until death, given a mortality table and an interest
#' rate curve.
#'
#' @param qx Matrix with probabilities of death to be used. (Generation table)
#' @param i Named numeric vector with the interest rates to be used in given
#' years.
#' @param m Number of payments per year
#' @param k costs to be applied
#' @return Matrix with the present values of the annuities.
#'
#' @examples
#' äx(matrix(rep(0.5, 4), 2, 2, dimnames=list(108:109, 2015:2016)),
#' set_names(1:3/100, 2123:2125))
äx <- function(qx, i, m, k){
  out <- qx * 0
  v_mat <- make_v(qx, i)
  out[nrow(out), ] <- 1 + v_mat[nrow(out), ] * (1 - qx[nrow(out), ])
  for (j in (nrow(out) - 1):1) {
    out[j, ] <- 1 + v_mat[j, ] * out[j + 1, ] * (1 - qx[j, ])
  }
  out <- (out - (m - 1) / (2*m)) * (1 + k)
  rownames(out) <- rownames(qx)
  colnames(out) <- colnames(out)
  return(out)
}


#' Orphan's pension
#'
#' Present value of an orphan's pension.
#' Assumes no mortality of the orphan.
#' kx znd zx are very simple extrapolations from the numbers given in Koller.
#'
#' @param qx Generation Life Table to be used
#' @param i Fixed interest rate.
#' @param m Number of payments per year
#' @param k costs to be applied
#' @return Value of the orphan's pensions
#'
#' @examples
#' z(17, 0.02)
äk <- function(qx, i, m, k){
  kx <- function(x) pmin(0.24, pmax(0, exp(-(x - 60)/4) * 0.24))
  zx <- function(x) pmax(17.4, pmin(25, (x - 60)/2 + 17.4))
  äz <- function(z, i){
    normal <- -((1 + i) ^ (z - 24) - 1) / log(1 + i)
    normal[i == 0] <- 24 - z
    return(normal)
    }

  out <- qx * 0
  v_mat <- make_v(qx, i)
  ages <- as.numeric(rownames(out))
  out[nrow(out), ] <- qx[nrow(out), ] *
    v_mat[nrow(out), ] ^ 0.5 *
    kx(ages[length(ages)]) *
    äz(zx(ages[length(ages)]), 1 / v_mat[nrow(out), ] - 1)
  for (j in (nrow(out) - 1):1) {
    age <- ages[j]
    out[j, ] <- qx[j, ] *
      v_mat[j, ] ^ 0.5 *
      kx(age) *
      äz(zx(age), 1 / v_mat[j, ] - 1) +
      (1 - qx[j, ]) * v_mat[j, ] *  out[j + 1, ]
  }
  out <- out * (1 + k)
  rownames(out) <- rownames(qx)
  colnames(out) <- colnames(out)
  return(out)
}


#' Widow's pension
#'
#' Calculates the present value of widows pensions. Works just as well
#' for widower's pensions.
#'
#' @param qx Matrix with probabilities of death of the person on whose live
#' the start of the widow's pension depends.
#' @param äyw Matrix of present value of widow's pensions, once they have
#' started. In a first approximation, normal annuities can be used, although
#' this may slightly overestimate the present values for widow(er)s.
#' @param i Named numeric vector with the interest rates to be used in given
#' years.
#' @param yx a function yielding the age of the widow(er) if the insured dies at
#' age y.
#' @param hx probability of leaving an entitled widow(er) if dying at age x.
#' @param m Number of payments per year
#' @param k costs to be applied
#' @return Matrix with the present values of the widow(er)'s annuities.
#'
#' @examples
#' äx(matrix(rep(0.5, 4), 2, 2, dimnames=list(108:109, 2015:2016)),
#' set_names(1:3/100, 2123:2125))
äxw <- function(qx, äyw, i, yx, hx, m, k){
  out <- qx * 0
  v_mat <- make_v(qx, i)
  omega <- as.numeric(rownames(qx)[nrow(qx)])
  get_indices <- function(i, mat, d){
    row <- min(nrow(mat), max(1, i + d))
    cols <- pmin(ncol(mat), pmax(1, 1:ncol(mat) - d))
    return(cbind(row, cols))
  }
  get_d <- function(x) yx(x) - x
  out[nrow(out), ] <- qx[nrow(out), ] * v_mat[nrow(out), ] *
    äyw[get_indices(nrow(out), qx, get_d(omega))]
  for (x in (nrow(qx) - 1):1) {
    out[x, ] <- (1 - qx[x, ]) * v_mat[x, ] * out[x + 1, ] +
      qx[x, ] * v_mat[x, ] * äyw[get_indices(x, qx, get_d(as.numeric(rownames(qx)[x])))]
  }
  out <- (out - (m - 1) / (2 * m)) * (1 + k)
  return(out)
}