#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.