BootMackChainLadder.R

library(ChainLadder)
library(mvtnorm)
##### Usefull functions                                                                #####
.diag <- function(Triangle) {
  # sympli returns the diagonal of the Triangle
  unname(rev(diag(Triangle[nrow(Triangle):1, ])))
}
.DF <- function(Triangle) {
  # Simply returns chain-ladder simples developpement factors.
  n <- dim(Triangle)[1]
  Triangle2 <- Triangle
  Triangle2[ col(Triangle2) == n - row(Triangle2) + sum(!is.na(Triangle2[, n]))  ] <- NA
  return(unname(c(colSums(Triangle[, 2:n], na.rm = TRUE) / colSums(Triangle2[, 1:(n - 1)], na.rm = TRUE), 1)))
}
.ultimates <- function(Triangle, df = .DF(Triangle)) {
  # simply returns chain-ladder ultimates.
  .diag(Triangle) * cumprod(rev(df))
}
.ibnr <- function(Triangle, df = .DF(Triangle)) {
  # simply returns the IBNR from classical chain-ladder model
  .diag(Triangle) * (cumprod(rev(df)) - 1)
}
.DFIndiv <- function(Triangle) {
  # Simply returns the Triangle of individual developpements factors.
  unname((cbind(Triangle, NA) / cbind(NA, Triangle))[, 2:(dim(Triangle)[[1]] + 1)])
}
.sigma <- function(Triangle, df=.DF(Triangle), DFIndiv = .DFIndiv(Triangle)) {
  # Returns Mack's estimate for the variances parameters sigma_j
  # Returns a line vector, indexed by Triangle columns.
  n <- dim(Triangle)[[1]]
  ecart <- DFIndiv - matrix(df, nrow = n, ncol = n, byrow = TRUE)
  rez <- Triangle * ecart^2
  sigma <- colSums(rez, na.rm = TRUE) / ((n - 2):(-1))
  
  # the last value is approximated by mack's sugestion :
  sigma <- c(
    sigma[1:(n - 2)],
    min(
      sigma[n - 2]^2 / sigma[n - 3],
      sigma[n - 2],
      sigma[n - 3]
    )
  )
  if (is.nan(sigma[n - 1])) {
    sigma[n - 1] <- 0
  }
  # we return not sigma^2 but just sigma
  return(unname(sqrt(sigma)))
}
.residuals <- function(Triangle, centered = FALSE, DF=.DF(Triangle), DFIndiv=.DFIndiv(Triangle), sigma=.sigma(Triangle)) {
  # estimation of mack's residuals, based on England & Verral 2007
  # returns a Triangle of same size as the input, but with NA on the diagonal (no residuals for the diagonal)
  residus <- Triangle
  n <- dim(Triangle)[1]
  
  for (i in 1:n) {
    for (j in 1:n) {
      residus[i, j] <- sqrt((n - j) / (n - j - 1)) * sqrt(Triangle[i, j]) * (DFIndiv[i, j] - DF[j]) / sigma[j]
    }
  }
  residus[is.nan(residus)] <- 0
  if (centered) {
    residus <- apply(residus, 2, function(x) {
      return(x - mean(x, na.rm = TRUE))
    })
  }
  return(residus)
}
.sampling <- function(Triangle, N=100, seuil = NA, zonnage=FALSE,morceaux=NULL) {
  # Fonction qui bootstrap un Triangle de positions
  # On fabrique une liste de Triangle de r?sidus bootstrapp?s :
  
  n <- length(Triangle)
  I <- dim(Triangle)[1]
  
  # Gestion du seuil d'exclusion des residus :
  Triangle2 <- Triangle
  if (!is.na(seuil)) {
    Triangle2[Triangle[(!is.na(Triangle))] > seuil] <- NA
  }
  if (zonnage) {
    # Il faut alors zonner les residus
    
    samples <- lapply(1:length(morceaux), function(i) {
      prob <- !is.na(Triangle2[, morceaux[[i]]])
      positions <- matrix(1:n, nrow = I)
      return(lapply(1:N, function(x) {
        matrix(sample(positions[, morceaux[[i]]],
                      replace = TRUE,
                      prob = prob,
                      size = length(prob)
        )
        ,
        nrow = I
        )
      }))
    })
    
    rez <- lapply(1:N, function(i) {
      return(do.call(cbind, lapply(samples, function(x) {
        x[[i]]
      })))
    })
  } else {
    prob <- !is.na(Triangle2)
    positions <- matrix(1:n, nrow = I)
    rez <- sample(positions, replace = TRUE, prob = prob, size = N * n)
    rez <- lapply(1:N, function(.x) {
      matrix(rez[(n * (.x - 1) + 1):(n * .x)], nrow = I)
    })
  }
  
  rez <- lapply(rez, function(x) {
    x[is.na(Triangle)] <- NA
    return(x)
  })
  return(rez)
}
.applysample <- function(sample, Triangle) {
  # Fonction d'application de cet sampling :
  return(ChainLadder::as.triangle(
    matrix(as.vector(Triangle)[as.vector(sample)],
           nrow = dim(Triangle)
    )
  ))
}
.DFPond <- function(DFIndiv, ponderation) {
  n <- dim(ponderation)[1]
  ponderation[col(ponderation) == n - row(ponderation) + 1] <- NA
  rez <- colSums(ponderation * DFIndiv, na.rm = TRUE) / colSums(ponderation, na.rm = TRUE)
  rez[n] <- 1
  return(rez)
}
.newDFPond <- function(NyDFIndiv, DFIndiv, Triangle) {
  n <- dim(Triangle)[1]
  DFIndiv[col(DFIndiv) == n - row(DFIndiv) + 1] <- c(NyDFIndiv, 1)
  rez <- colSums(Triangle * DFIndiv, na.rm = TRUE) / colSums(Triangle, na.rm = TRUE)
  rez[n] <- 1
  return(rez)
}
.formatOutput <- function(n, Triangle, diag, DF, sigma, residuals, DFBoot, Ultimates, IBNR, NyCum, NyInc, NyDF, NyUltimates, NyIBNR) {
  # output formatting :
  NyCum <- do.call(rbind, lapply(NyCum, rev))
  NyInc <- do.call(rbind, lapply(NyInc, rev))
  NyUltimates <- do.call(rbind, lapply(NyUltimates, rev))
  NyIBNR <- do.call(rbind, lapply(NyIBNR, rev))
  NyDF <- do.call(rbind, NyDF)
  DFBoot <- do.call(rbind, DFBoot)
  
  names(Ultimates) <- rownames(Triangle)
  names(IBNR) <- rownames(Triangle)
  names(diag) <- rownames(Triangle)
  colnames(NyCum) <- rownames(Triangle)[2:(n)]
  colnames(NyInc) <- rownames(Triangle)[2:(n)]
  colnames(NyUltimates) <- rownames(Triangle)
  colnames(NyIBNR) <- rownames(Triangle)
  colnames(NyDF) <- colnames(Triangle)
  colnames(DFBoot) <- colnames(Triangle)
  
  result <- list(
    Latest = rev(diag),
    DF = DF,
    sigma = sigma,
    residuals = residuals,
    DFBoot = DFBoot,
    Ultimates = Ultimates,
    IBNR = IBNR,
    NyCum = NyCum,
    NyInc = NyInc,
    NyDF = NyDF,
    NyUltimates = NyUltimates,
    NyIBNR = NyIBNR
  )
  class(result) <- c("BootMackChainLadder", class(result))
  return(result)
}
##### Boot Mack Chain Ladder                                                         #####
#' Boot Mack Chain Ladder model
#'
#' This function implement a simple bootstrap of the residuals from the mack model with a one-year reserving risk point of view.
#' Cf "One-year reserve risk including a tail factor : closed formula and bootstrapp aproach, 2011"
#'
#' @param Triangle A simple triangle from che Chain-ladder package.
#' @param B numeric. The number of bootstrap samples you want
#' @param distNy character. Distribution of next-year incremental payments. Either "normal" (default) or "residuals"
#' @param seuil numeric. Seuil for residuals exclusion. A value of NA (default) will prevent
#' @param zonnage logical. Do you want to zonne the residuals ?
#'
#' @return A BootMackChainLadder object with a lot of information about the bootstrapping. You can plot it, print it, str it to extract information.
#' @export
#'
#' @import magrittr
#'
#' @examples
#' data(ABC)
#' BootMackChainLader(Triangle = ABC, B = 100, distNy = "residuals", seuil = 2)
BootMackChainLadder <- function(Triangle, B=100, distNy="normal", seuil=NA, zonnage=FALSE,morceaux = NULL ) {
  if (!(distNy %in% c("normal", "residuals"))) {
    stop("DistNy Parameter must be 'normal' (classical MW) or 'residuals'")
  }
  # Cf "One-year reserve risk including a tail factor : closed formula and bootstrapp aproach, 2011"
  # First step : Mack model.
  n <- dim(Triangle)[1]
  diag <- rev(.diag(Triangle))
  DF <- .DF(Triangle)
  DFIndiv <- .DFIndiv(Triangle)
  sigma <- .sigma(Triangle, DF, DFIndiv)
  residuals <- .residuals(Triangle, centered = TRUE, DF, DFIndiv, sigma)
  Ultimates <- .ultimates(Triangle, DF)
  IBNR <- .ibnr(Triangle, DF)
  
  # Step 2 : Resampling residuals.
  samples <- .sampling(residuals, B, seuil = seuil, zonnage = zonnage,morceaux=morceaux)
  sampledResiduals <- lapply(samples, .applysample, residuals)
  # Step3 : Calculation of booted estimators
  DFIndivBoot <- lapply(sampledResiduals, function(.x) {
    t(t(.x * sqrt(t(c(sigma^2, 0) / t(Triangle)))) + DF)
  })
  DFBoot <- lapply(DFIndivBoot, .DFPond, Triangle) # ponderated by the original triangle !
  
  # Step 5 : Simulation of NY payments by a standard normal OR by the residuals... taking into acount process error
  if (distNy == "normal") {
    NyCum <- lapply(DFBoot, function(.x) {
      rnorm(n = n - 1, mean = (diag * .x)[1:(n - 1)], sd = sqrt((diag * c(sigma^2, 0))[1:(n - 1)]))
    })
  }
  if (distNy == "residuals") {
    # To simulate the same thing with residuals assuming they are all i.i.d between developpements years, we could do :
    if (!is.na(seuil)) { # FIrst, let's discard non-selected residuals.
      residuals[residuals[(!is.na(residuals))] > seuil] <- NA
    }
    if (!zonnage) {
      samples <- sample(residuals[!is.na(residuals)], size = (n - 1) * B, replace = TRUE)
      samples <- lapply(1:B, function(.x) {
        samples[((n - 1) * (.x - 1) + 1):((n - 1) * .x)]
      })
      NyCum <- mapply(function(.y, .x) {
        (.y * sqrt((diag * c(sigma^2, 0))[1:(n - 1)])) + (diag * .x)[1:(n - 1)]
      }, samples, DFBoot, SIMPLIFY = FALSE)
    } else {
      # If zonnig of residuals is needed, the following zonning will be implemented :
      
      samples <- list()
      
      samples <- lapply(1:length(morceaux), function(i) {
        x <- sample(residuals[, morceaux[[i]]][!is.na(residuals[, morceaux[[i]]])], size = length(morceaux[[i]]) * B, replace = TRUE)
        x <- lapply(1:B, function(.x) {
          x[(length(morceaux[[i]]) * (.x - 1) + 1):(length(morceaux[[i]]) * .x)]
        })
        return(x)
      })
      
      samples2 <- lapply(1:B, function(i) {
        plop <- samples[[1]][[i]]
        for (j in 2:length(morceaux)) {
          plop <- c(plop, samples[[j]][[i]])
        }
        return(plop)
      })
      
      NyCum <- mapply(function(.y, .x) {
        (.y * sqrt((diag * c(sigma^2, 0))[1:(n - 1)])) + (diag * .x)[1:(n - 1)]
      }, samples2, DFBoot, SIMPLIFY = FALSE)
    }
  }
  
  # NyCum is list of B vectors of size n-1, containing the next year cumulative payments, from origin n to 1 ( decreasing order)
  
  # Step 6 : Calculation of Ny estimators
  NyInc <- lapply(NyCum, function(.x) {
    .x - diag[1:(n - 1)]
  }) # Coresponding incremnts
  NyDFIndiv <- lapply(NyCum, function(.x) {
    .x / diag[1:(n - 1)]
  }) # Coresponding individual DF's
  NyDF <- lapply(NyDFIndiv, .newDFPond, DFIndiv, Triangle) # coredponing DF
  NyUltimates <- mapply(function(.x, .y) {
    c(.x * rev(cumprod(rev(.y[2:n]))), Triangle[1, n])
  }, NyCum, NyDF, SIMPLIFY = FALSE)
  NyIBNR <- lapply(NyUltimates, function(.x) {
    .x - diag
  })
  
  rez <- .formatOutput(n, Triangle, diag, DF, sigma, residuals, DFBoot, Ultimates, IBNR, NyCum, NyInc, NyDF, NyUltimates, NyIBNR)
  return(rez)
}

#' mean.BootMackChainLadder
#'
#' Calculate mean statiscics from the bootstraped mack model
#'
#' @param x A BootMackChainLadder Object
#'
#' @return Three data.frames containing data per origin year ("ByOrigin"), by developpement year ("ByDev) and globaly ("Totals)
#' @export
#'
#' @examples
#' data(ABC)
#' BMCL <- BootMackChainLader(Triangle = ABC, B = 100, distNy = "residuals", seuil = 2)
#' mean(BMCL)
mean.BootMackChainLadder <- function(x, ...) {
  
  # the purpose of this function is to return means of everything BMCL returned.
  ByOrigin <- data.frame(
    Latest = x$Latest,
    Ultimates = x$Ultimates,
    IBNR = x$IBNR,
    NyCum = c(x$Ultimates[1], colMeans(x$NyCum)),
    NyInc = c(0, colMeans(x$NyInc)),
    NyUltimates = colMeans(x$NyUltimates),
    NyIBNR = colMeans(x$NyIBNR)
  )
  row.names(ByOrigin) <- rev(row.names(ByOrigin))
  
  ByDev <- data.frame(
    DFBoot = colMeans(x$DFBoot),
    NyDF = colMeans(x$NyDF)
  )
  
  Totals <- colSums(ByOrigin)
  
  
  return(list(ByOrigin = ByOrigin, ByDev = ByDev, Totals = Totals))
}
#' CDR.BootMackChainLadder
#'
#' Calculate the one-year Claim developpement results from the bootstrapped mack model.
#'
#' @param BMCL A BootMackChainLadder Object
#'
#' @return A data.Frame with one-year results from the bootstrap.
#' @export
#'
#' @examples
#' data(ABC)
#' BMCL <- BootMackChainLader(Triangle = ABC, B = 100, distNy = "residuals", seuil = 2)
#' CDR(BMCL)
CDR.BootMackChainLadder <- function(BMCL) {
  Totals <- data.frame(IBNR = sum(BMCL$IBNR), `CDR(1)S.E.` = sd(rowSums(BMCL$NyIBNR)))
  row.names(Totals) <- "Totals"
  names(Totals) <- c("IBNR", "CDR(1)S.E.")
  rbind(setNames(data.frame(IBNR = BMCL$IBNR, `CDR(1)S.E.` = apply(BMCL$NyIBNR, 2, sd)), names(Totals)), Totals)
}
#' summary.BootMackChainLadder
#'
#'  Give summary statistics about the bootstrap.
#'
#' @param object A BootMackChainLadder Object
#'
#' @return NULL
#' @export
#'
#' @details
#' The function only print information about the model.
#'
#' @examples
#' data(ABC)
#' BMCL <- BootMackChainLader(Triangle = ABC, B = 100, distNy = "residuals", seuil = 2)
#' summary(BMCL)
summary.BootMackChainLadder <- function(object, ...) {
  mean <- mean(object)
  cat("This is a BootMackChainLadder model \n\n")
  cat("Detailed results by Origin years : \n")
  print(format(mean$ByOrigin, format = "i", nsmall = 0, big.mark = ","))
  # print(mean$ByDev)
  cat("\n Totals across origin years : \n")
  print(format(t(mean$Totals), format = "i", nsmall = 0, big.mark = ","))
  return(NULL)
}
#' print.BootMackChainLadder
#'
#' @param BMCL
#'
#' @return the BMCL object
#' @export
#'
#' @examples
#' data(ABC)
#' BMCL <- BootMackChainLader(Triangle = ABC, B = 100, distNy = "residuals", seuil = 2)
#' print(BMCL)
print.BootMackChainLadder <- function(x, ...) {
  print(summary(x))
  return(NULL)
}

##### Tests                                                                          #####


# test :
library(magrittr)

data(ABC)
data(auto)
MCL <- MackChainLadder(ABC,est.sigma="Mack")
BMCL1 <- BootMackChainLadder(ABC,B=10000,distNy = "normal")
BMCL2 <- BootMackChainLadder(ABC,B=10000,distNy = "residuals")
BMCL3 <- BootMackChainLadder(ABC,B=10000,distNy = "residuals",zonnage = TRUE,morceaux = list(1:2,3:5,6:10))

# r?sum? :
BMCL1
BMCL2
BMCL3
MCL

# cdr :
plot(as.numeric(rownames(CDR(MCL))[1:11]),CDR(MCL)$`CDR(1)S.E.`[1:11],type="l")
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL1)$`CDR(1)S.E.`[1:11],type="l",col=2)
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL2)$`CDR(1)S.E.`[1:11],type="l",col=3)
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL3)$`CDR(1)S.E.`[1:11],type="l",col=4)

# the convergence is neat :)


# lets try on other triangles :
data(MW2008)

MCL <- MackChainLadder(MW2008,est.sigma="Mack")
BMCL1 <- BootMackChainLadder(MW2008,B=10000,distNy = "normal")
BMCL2 <- BootMackChainLadder(MW2008,B=10000,distNy = "residuals")
BMCL3 <- BootMackChainLadder(MW2008,B=10000,distNy = "residuals",zonnage = TRUE,morceaux = list(1:3,4:8))
# r?sum? :
BMCL1
BMCL2
BMCL3
MCL

# cdr :
plot(as.numeric(rownames(CDR(MCL))[1:11]),CDR(MCL)$`CDR(1)S.E.`[1:11],type="l")
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL1)$`CDR(1)S.E.`[1:11],type="l",col=2)
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL2)$`CDR(1)S.E.`[1:11],type="l",col=3)
points(as.numeric(rownames(CDR(MCL))[1:11]),CDR(BMCL3)$`CDR(1)S.E.`[1:11],type="l",col=4)