# R/BootstrapReserve.R In edalmoro/ChainLadderQuantileV1: Statistical Methods and Models for Claims Reserving in General Insurance

```## Author: Markus Gesmann
## Copyright: Markus Gesmann, [email protected]
## Date:19/10/2008

BootChainLadder <- function(Triangle, R = 999, process.distr=c("gamma", "od.pois")){

if(!'matrix' %in% class(Triangle))
Triangle <- as.matrix(Triangle)

process.distr <- match.arg(process.distr)
weights <- 1/Triangle

triangle <- Triangle
if(nrow(triangle) != ncol(triangle))
stop("Number of origin periods has to be equal or greater then the number of development periods.\n")

triangle <- checkTriangle(Triangle)
output.triangle <- triangle

m <- dim(triangle)[1]
n <- dim(triangle)[2]
origins <- c((m-n+1):m)

## Obtain the standard chain-ladder development factors from cumulative data.

triangle <- array(triangle, dim=c(m,n,1))
weights <-  array(weights, dim=c(m,n,1))
inc.triangle <- getIncremental(triangle)
Latest <- getDiagonal(triangle,m)

## Obtain cumulative fitted values for the past triangle by backwards
## recursion, starting with the observed cumulative paid to date in the latest
## diagonal

dfs <- getIndivDFs(triangle)
avDFs <- getAvDFs(dfs, 1/weights)
ultDFs <- getUltDFs(avDFs)
ults <- getUltimates(Latest, ultDFs)
## Obtain incremental fitted values, m^ ij, for the past triangle by differencing.
exp.inc.triangle <- getIncremental(getExpected(ults, 1/ultDFs))
## exp.inc.triangle[is.na(inc.triangle[,,1])] <- NA
exp.inc.triangle[is.na(inc.triangle[origins,,1])] <- NA
## Calculate the unscaled Pearson residuals for the past triangle using:
inc.triangle <- inc.triangle[origins,,]
dim(inc.triangle) <- c(n,n,1)
unscaled.residuals  <- (inc.triangle - exp.inc.triangle)/sqrt(abs(exp.inc.triangle))

## Calculate the Pearson scale parameter
nobs  <- 0.5 * n * (n + 1)
scale.factor <- (nobs - 2*n+1)
scale.phi <- sum(unscaled.residuals^2,na.rm=TRUE)/scale.factor
## Adjust the Pearson residuals using

## Sample incremental claims
## Resample the adjusted residuals with replacement, creating a new
## past triangle of residuals.

## Fit the standard chain-ladder model to the pseudo-cumulative data.
## Project to form a future triangle of cumulative payments.

simCum <- makeCumulative(simClaims)
simLatest <- getDiagonal(simCum,nrow(simCum))
simDFs <- getIndivDFs(simCum)
simWghts <- simCum
simAvDFs <- getAvDFs(simDFs, simWghts)
simUltDFs <- getUltDFs(simAvDFs)
simUlts <- getUltimates(simLatest, simUltDFs)
## Get expected future claims
## Obtain the corresponding future triangle of incremental payments by
## differencing, to be used as the mean when simulating from the process
## distribution.

simExp <- getIncremental(getExpected(simUlts, 1/simUltDFs))
simExp[!is.na(simClaims)] <- NA
processTriangle <-array(NA,c(n,n,R))

if(process.distr=="gamma")
processTriangle[!is.na(simExp)] <- sign(simExp[!is.na(simExp)])*rgamma(length(simExp[!is.na(simExp)]), shape=abs(simExp[!is.na(simExp)]/scale.phi), scale=scale.phi)

if(process.distr=="od.pois")
processTriangle <-  apply(simExp,c(1,2,3), function(x)
ifelse(is.na(x), NA, sign(x)*rpois.od(1, abs(x), scale.phi)))
##if(process.distr=="od.nbionm")
##  processTriangle <-  apply(simExp,c(1,2,3), function(x)
##          ifelse(is.na(x), NA, sign(x)*rnbinom.od(1, mu=abs(x), size=sqrt(1+abs(x)),d=scale.phi)))

processTriangle[is.na(processTriangle)] <- 0
simExp[is.na(simExp)] <- 0

IBNR.Triangles <- processTriangle
IBNR <- makeCumulative(IBNR.Triangles)[,n,]
dim(IBNR)<-c(n,1,R)
ParamDist <- makeCumulative(simExp)[,n,]
dim(ParamDist)<-c(n,1,R)

# Giuseppe Crupi
# Generate and process Next Year triangle for one year re-reserving approach (Solvency 2 purpose)

NYCost<-getNYCost(triangle[origins,,,drop=F],IBNR.Triangles,R)
NYParamDist<-getNYCost(triangle[origins,,,drop=F],simExp,R)

if(m>n){
IBNR <- apply(IBNR, 3, function(x) c(rep(0, m-n),x))
dim(IBNR) <- c(m,1,R)

NYCost <- apply(NYCost, 3, function(x) c(rep(0, m-n),x))
dim(NYCost) <- c(m,1,R)

ParamDist <- apply(ParamDist, 3, function(x) c(rep(0, m-n),x))
dim(ParamDist) <- c(m,1,R)

NYParamDist <- apply(NYParamDist, 3, function(x) c(rep(0, m-n),x))
dim(NYParamDist) <- c(m,1,R)

IBNR.Triangles <- apply(IBNR.Triangles,3, function(x) rbind(matrix(0,m-n,n),x))
dim(IBNR.Triangles) <- c(m,n,R)

simClaims <- apply(simClaims,3, function(x) rbind(matrix(0,m-n,n),x))
dim(simClaims) <- c(m,n,R)
}

IBNR.Totals <- apply(IBNR.Triangles,3,sum)

output <- list( call=match.call(expand.dots = FALSE),
Triangle=output.triangle,
f=as.vector(avDFs),
simClaims=simClaims,
IBNR.ByOrigin=IBNR,
IBNR.Triangles=IBNR.Triangles,
IBNR.Totals = IBNR.Totals,
ParamDist.ByOrigin=ParamDist,
NYCost.ByOrigin = NYCost,
NYParamDist.ByOrigin=NYParamDist,
#NYIBNR.ByOrigin = NYIBNR,
process.distr=process.distr,
R=R)

return(output)

}

############################################################################
##

}
############################################################################
##
quantile.BootChainLadder <- function(x,probs=c(0.75, 0.95), na.rm = FALSE,
names = TRUE, type = 7,...){

m <- dim(x\$Triangle)[1]
n <- dim(x\$Triangle)[2]

ByOrigin <- t(apply(x\$IBNR.ByOrigin, 1, quantile, probs=probs, names=names, type=type,...))

if(length(probs)>1){
ByOrigin <- as.data.frame(ByOrigin)
}else{
ByOrigin <- as.data.frame(t(ByOrigin))
}
names(ByOrigin) <- paste("IBNR ", probs*100, "%", sep="")

origin <- dimnames(x\$Triangle)[[1]]

if(length(origin)==nrow(ByOrigin)){
rownames(ByOrigin) <- origin
}

Total.IBNR.q <- quantile(x\$IBNR.Totals, probs=probs, names=names, type=type, ...)

Totals <- as.data.frame(Total.IBNR.q)

colnames(Totals)=c("Totals")
rownames(Totals) <- paste("IBNR ", probs*100, "%:", sep="")

output <- list(ByOrigin=ByOrigin, Totals=Totals)

return(output)
}
############################################################################
##

m <- dim(x\$Triangle)[1]
n <- dim(x\$Triangle)[2]

ByOrigin <- apply(x\$IBNR.ByOrigin, 1, mean,...)

ByOrigin <- as.data.frame(ByOrigin)

names(ByOrigin) <- "Mean IBNR"

origin <- dimnames(x\$Triangle)[[1]]
if(length(origin)==nrow(ByOrigin)){
rownames(ByOrigin) <- origin
}

Total.IBNR <- mean(x\$IBNR.Totals,...)
Totals <- as.data.frame(Total.IBNR)

colnames(Totals)=c("Total")
rownames(Totals) <- "Mean IBNR:"

output <- list(ByOrigin=ByOrigin, Totals=Totals)

return(output)
}

############################################################################
##

.Triangle <- object\$Triangle
m <- dim(.Triangle)[1]
n <- dim(.Triangle)[2]

dim(.Triangle) <- c(dim(.Triangle),1)
Latest <- as.vector(getLatest(getIncremental(.Triangle)))

IBNR <- object\$IBNR.ByOrigin
dim(IBNR) <- dim(IBNR)[c(1,3)]
IBNR.q <- apply(IBNR, 1, quantile, probs=probs,...)
IBNR.mean <- apply(IBNR, 1, mean)
IBNR.sd <- apply(IBNR, 1, sd)
sumIBNR <- t(rbind(IBNR.mean, IBNR.sd, IBNR.q))
sumIBNR <- as.data.frame(sumIBNR)

## ByOrigin
ByOrigin <- data.frame(Latest,
Ult.mean=Latest+sumIBNR\$IBNR.mean,
sumIBNR)
names(ByOrigin) <- c("Latest", "Mean Ultimate", "Mean IBNR",
"SD IBNR", paste("IBNR ", probs*100, "%", sep=""))
ex.origin.period <- Latest!=0
ByOrigin <- ByOrigin[ex.origin.period,]

origin <- dimnames(object\$Triangle)[[1]]
if(length(origin)==nrow(ByOrigin)){
rownames(ByOrigin) <- origin
}

## Totals
Total.Latest <- sum(Latest,na.rm=TRUE)
Total.IBNR <- object\$IBNR.Totals
Total.IBNR.mean <-  mean(Total.IBNR)
Total.IBNR.sd <-  sd(Total.IBNR)
Total.IBNR.q <- quantile(Total.IBNR, probs=probs,...)

Totals <-  c(Total.Latest, Total.Latest+Total.IBNR.mean,
Total.IBNR.mean, Total.IBNR.sd, Total.IBNR.q)
Totals <- as.data.frame(Totals)

colnames(Totals)=c("Totals")
rownames(Totals) <- c("Latest:","Mean Ultimate:",
"Mean IBNR:","SD IBNR:",
paste("Total IBNR ", probs*100, "%:", sep="") )

output <- list(ByOrigin=ByOrigin, Totals=Totals)
return(output)
}
############################################################################
##
print(x\$call)
cat("\n")
summary.x <- summary(x,probs=probs)

names(summary.x\$ByOrigin)[4] <- "IBNR.S.E"
print(format(summary.x\$ByOrigin, big.mark = ",", digits = 3),...)

cat("\n")
Totals <- summary.x\$Totals
rownames(Totals)[4] <- "IBNR.S.E"
print(format(Totals, big.mark=",",digits=3), quote=FALSE)

}

makeCumulative <- function(tri){
## Author: Nigel de Silva - Giuseppe Crupi
for (i in 2:dim(tri)[2])
tri[,i,]<-tri[,i-1,]+tri[,i,]
return(tri)
}

getIncremental <- function(tri){
## Author: Nigel de Silva - Giuseppe Crupi
n<-dim(tri)[2]
tri[,2:n,]<-tri[,2:n,]-tri[,1:(n-1),]

return(tri)
}

getLatest <- function(incr.tri){
## Author: Nigel de Silva
out <- apply(incr.tri, c(1,3), sum, na.rm=T)
dim(out) <- c(1, dim(out))
out <- aperm(out, c(2,1,3))
return(out)
}

getIndivDFs <- function(tri){
## Author: Nigel de Silva
.nDevPrds <- dim(tri)[2]
.numerator <- tri[ , c(2:.nDevPrds, .nDevPrds), , drop=F]
tri <- .numerator / tri
tri[ , .nDevPrds, ] <- NA
return(tri)
}

getAvDFs <- function(dfs, wghts){
## Author: Nigel de Silva - Giuseppe Crupi
.include <- dfs
.include[!is.na(.include)] <- 1
.include[is.na(wghts)] <- NA
.include[is.na(.include)] <- 0
dfs[is.na(dfs)]<-0
wghts[is.na(wghts)]<-0

out<-colSums(dfs * wghts * .include)/colSums(wghts * .include)
out[is.nan(out)] <- 1

out <- array(out, c(1,dim(out)))
return(out)

}

getUltDFs <- function(avDFs){
## Author: Nigel de Silva - Giuseppe Crupi
for (i in seq.int(dim(avDFs)[2]-1,1,-1))
avDFs[,i,]<-avDFs[,i+1,]*avDFs[,i,]
return(avDFs)

}

getUltimates <- function(latest, ultDFs){
## Author: Nigel de Silva - Giuseppe Crupi
m<-dim(ultDFs)[1]
n<-dim(ultDFs)[2]
r<-dim(ultDFs)[3]

ultDFs <- ultDFs[,(dim(ultDFs)[2]):1,]
dim(ultDFs)<-c(n,1,r)
out <- latest * ultDFs

return(out)

}

expandArray <- function(x, d, new.size){
## Author: Nigel de Silva
## Expand array x, in d th dimension to new.size
## Permute dimension to set changing dimension last
.n <- length(dim(x))
.perm <- 1:.n
.perm[d] <- .n
.perm[.n] <- d
x <- aperm(x, .perm)

## Expand dimension
x <- array(x, dim=c(dim(x)[-.n], new.size))

x <- aperm(x, .perm)

return(x)

}

getExpected <- function(ults, ultDFs){
## Author: Nigel de Silva
## get backward chain ladder incremental triangles
ults <- expandArray(ults, 2, dim(ultDFs)[2])
ultDFs <- expandArray(ultDFs, 1, dim(ults)[1])

return(ults * ultDFs)
}

sampleResiduals <- function(resids, positions, n.sims){
## Author: Nigel de Silva
## Worry about excluding, zoning, etc. residuals later

resids <- as.vector(resids)
resids <- resids[!is.na(resids)]

.sample <- sample(resids, prod(dim(positions)[-3])*n.sims, replace=T)
.sample <- array(.sample, dim=c(dim(positions)[-3], n.sims))
.sample[is.na(positions)] <- NA

return(.sample)
}

randomClaims <- function(exp.clms, resids, n.sims){
## Author: Nigel de Silva
.residSample <- sampleResiduals(resids, exp.clms, n.sims)
exp.clms <- expandArray(exp.clms, 3, n.sims)
out <- .residSample * sqrt(abs(exp.clms)) + exp.clms

return(out)
}

# Mike Lonergan mel at mcs.st-and.ac.uk
# Fri Jun 14 19:31:24 CEST 2002 https://stat.ethz.ch/pipermail/r-help/2002-June/022425.html
# functions to generate random numbers following
# over-dispersed Poisson and negative binomial distributions
# based on simple cheat of using a standard negative binomial,
# but choosing the scale parameter to give the desired mean/variance
# ratio at the given value of the mean.

rpois.od<-function (n, lambda,d=1) {
if (d==1 | lambda<=0)
rpois(n, lambda)
else
rnbinom(n, size=(lambda/(d-1)), mu=lambda)
}

rnbinom.od<-function (n, size, prob, mu, d=1) {
if (!missing(prob)) {
if (!missing(mu))
stop("prob and mu both specified")
mu<-size*(1-prob)/prob
}
size2 <- mu/(d-1+(d*mu/size))
prob2 <- size2/(size2 + mu)
rnbinom(n, size2, prob2)
}

############################################################################
##
x, mfrow=NULL, title=NULL, log=FALSE,
which=1:4, ...){

if(is.null(title)){
myoma <- c(0,0,0,0)
}else{
myoma <- c(0,0,2,0)
}

Total.IBNR <- x\$IBNR.Total

if(is.null(mfrow)){
mfrow <- c(ifelse(length(which) < 2,1,
ifelse(length(which) < 3, 2,
ceiling(length(which)/2))),
ifelse(length(which)>2,2,1))
}

op=par(mfrow=mfrow, oma=myoma,...)

if(1 %in% which){
## Histogram
hist(Total.IBNR, xlab="Total IBNR")
lines(density(Total.IBNR))
rug(Total.IBNR)
}
if(2 %in% which){
## Empirical distribution
plot(ecdf(Total.IBNR), xlab="Total IBNR")
}
if(3 %in% which){
## Plot simultated Ultiamtes by Origin
plotBootstrapUltimates(x)
}
if(4 %in% which){
## Backtest last developemt year and check for cal. year trends
backTest <- backTestLatestIncremental(x)
plotLatestIncremental(backTest, log=log)
}
title( title , outer=TRUE)
par(op)
}

backTestLatestIncremental <- function(x){
stop("This function expects an object of class BootChainLadder")
# Simulated Latest Incremental
simLatest <- getLatest(getIncremental(x\$simClaims))

## Actual Latest Incremental
triangle <- x\$Triangle
dim(triangle) <- c(dim(triangle)[1:2],1)
actual.incr.Latest <- as.vector(getLatest(getIncremental(getIncremental(triangle))))
df.actual <- data.frame(actual.incr.Latest, origin=dimnames(x\$Triangle)[[1]])

Long.simLatest <- expand.grid(origin=dimnames(x\$Triangle)[[1]], sample=c(1:x\$R))
Long.simLatest\$sim.incr.Latest <- as.vector(simLatest)
LongSimActual <- merge(Long.simLatest, df.actual)
return(LongSimActual)
}

plotLatestIncremental <- function(LongSimActual,##sim.Latest,
##actual.Latest,
log=TRUE,
xlab="origin period",
ylab="latest incremental claims",
main="Latest actual incremental claims\nagainst simulated values",
col="bisque",
...){
## Plot the latest simulated incremental claims against actual observations to identify calendar year trends
## Example:
## x <- backTestLatestIncremental(B)
## plotLatestIncremental(x[[1]], x[[1]])

LongSimActual <- LongSimActual[LongSimActual\$sim.incr.Latest>0 & LongSimActual\$actual.incr.Latest>0,]
LongSimActual\$origin <- factor(LongSimActual\$origin)
if(log){
boxplot(LongSimActual\$sim.incr.Latest ~ LongSimActual\$origin, log="y",
xlab=xlab, ylab=paste("Log(",ylab,")", sep=""), main=main,col=col,...)
}else{
boxplot(LongSimActual\$sim.incr.Latest ~ LongSimActual\$origin,
xlab=xlab, ylab=ylab, main=main,col=col,...)
}
points(y=LongSimActual\$actual.incr.Latest,x=LongSimActual\$origin, col=2, pch=19)
legend("topleft","Latest actual", pch=19, col=2)
}

plotBootstrapUltimates <- function(x, xlab="origin period", ylab="ultimate claims costs",
main="Simulated ultimate claims cost",...){

triangle <- x\$Triangle
dim(triangle) <- c(dim(triangle),1)
Latest <- as.vector(getLatest(getIncremental(triangle)))

.origin <- dimnames(x\$Triangle)[[1]]
meanUlt <- data.frame(mean(x)\$ByOrigin+Latest, .origin)
names(meanUlt) <- c("meanUltimate", "origin")
Ultimate <- apply(x\$IBNR.ByOrigin, 3, function(.x) .x+Latest)

simUlt <- expand.grid(origin=.origin, sample=c(1:dim(x\$IBNR.ByOrigin)[3]))
simUlt\$Ultimate <- as.vector(Ultimate)

simUlt <- subset(merge(simUlt, meanUlt) , Ultimate > 0)
simUlt\$origin <- factor(simUlt\$origin)

boxplot(Ultimate ~ origin, data=simUlt,
xlab=xlab, ylab=ylab,main=main,col="bisque",...)

points(y=simUlt\$meanUltimate, x=simUlt\$origin,col=2, pch=19)
legend("topleft","Mean ultimate claim", pch=19, col=2)

}

getTriangleNextYear <- function(t.orig,t.ibnr){
## Author: Giuseppe Crupi

m<-dim(t.orig)[1]
n<-dim(t.orig)[2]
r<-dim(t.ibnr)[3]

out<-array(NA,c(m,n,r))
out[,,]<-getIncremental(t.orig)[,,]
IBNRIndexes<-getDiagonalIndexes(t.ibnr,m+1)
out[IBNRIndexes]<-t.ibnr[IBNRIndexes]
out<-makeCumulative(out)

return(out)

}

getDiagonalIndexes <-function(tri,d){
## Author: Giuseppe Crupi

n<-dim(tri)[1]
m<-dim(tri)[2]
r<-dim(tri)[3]
if (is.na(r)) r<-1

endRow<-min(d,n)
startRow<-max(1,d-m+1)

startCol<-min(d,m)
endCol<-max(1,d-n+1)

i<-startRow:endRow
j<-startCol:endCol
out<-array(0,c(length(i)*r,3))

out[,1]<-rep(i,r)
out[,2]<-rep(j,r)
out[,3]<-rep(1:r,each=length(i))

return(out)

}

getDiagonal <-function(tri,d){
## Author: Giuseppe Crupi

n<-dim(tri)[1]
m<-dim(tri)[2]
r<-dim(tri)[3]
if (is.na(r)) r<-1
dim(tri)<-c(n,m,r)

i<-getDiagonalIndexes(tri,d)
out<-tri[i]
dim(out)<-c(dim(i)[1]/r,1,r)

return(out)

}

getNYCost <-function(t.orig,t.ibnr,R){

m <- dim(t.orig)[1]
n <- dim(t.orig)[2]

triangleNY<-getTriangleNextYear(t.orig,t.ibnr)

NYLatest <- getDiagonal(triangleNY,m+1) ## Cumulative
NYPayments <- getDiagonal(t.ibnr,m+1) # Incremental

NYDFs <- getIndivDFs(triangleNY)

NYWghts <- triangleNY

NYAvDFs <- getAvDFs(NYDFs, NYWghts)

NYUltDFs <- getUltDFs(NYAvDFs)[,2:n,]
dim(NYUltDFs)<-c(1,n-1,R)

NYUlts <- getUltimates(NYLatest, NYUltDFs)

NYIBNR <- NYUlts - NYLatest

out <- array(0,c(m,1,R))

out[2:m,,] <- NYIBNR+NYPayments ## New rereserved ultimate

return(out)
}

CDR.BootChainLadder <- function(x, probs=c(0.75, 0.95), ...){
B <- x

IBNR <-  c(summary(B)\$ByOrigin[,3], summary(B)\$Totals[3,1])
IBNR.S.E <-  c(summary(B)\$ByOrigin[,4], summary(B)\$Totals[4,1])
CDR <- apply(B[["NYCost.ByOrigin"]],1, mean)
CDR.SD <- apply(B[["NYCost.ByOrigin"]],1, sd)
CDR.Totals <- apply(B[["NYCost.ByOrigin"]],3,sum)

CDR.Param.SD<-apply(B[["NYParamDist.ByOrigin"]],1, sd)
CDR.Param.Totals <-apply(B[["NYParamDist.ByOrigin"]],3,sum)

#CDR <- IBNR - c(CDR, mean(CDR.Totals))
CDR.Param.S.E <- c(CDR.Param.SD, sd(CDR.Param.Totals ))
CDR.S.E <- c(CDR.SD, sd(CDR.Totals ))
CDR.Process.S.E<-(CDR.S.E^2-CDR.Param.S.E^2)^0.5

if(length(probs)>1){
CDR.Q <-  t(cbind(apply(B[["NYCost.ByOrigin"]],1, quantile, probs),
quantile(CDR.Totals, probs))) #- IBNR
}else{
CDR.Q <-  c(apply(B[["NYCost.ByOrigin"]],1,quantile, probs),
quantile(CDR.Totals, probs)) #- IBNR
}

res <- data.frame(IBNR, IBNR.S.E,
#CDR.Param.S.E,
#CDR.Process.S.E,
CDR.S.E, CDR.Q)

names(res)[-c(1:2)] <- c("CDR(1)S.E",
paste0("CDR(1)", 100*probs,"%"))

rownames(res) <- c(rownames(B\$Triangle), "Total")
return(res)
}
```
edalmoro/ChainLadderQuantileV1 documentation built on Oct. 1, 2018, 12:23 a.m.