Last active
July 20, 2018 14:46
-
-
Save lrnv/01c37d66f7394e221a30d54873c37a3d to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment