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

```## Author: Markus Gesmann
## Copyright: Markus Gesmann, [email protected]
## Date:10/11/2007
## Date:08/09/2008
## Date:22/03/2009
## Date:25/05/2010 Dan

Triangle,
weights=1,
alpha=1,
est.sigma="log-linear",
tail=FALSE,
tail.se=NULL,
tail.sigma=NULL,
mse.method = "Mack")
{
## idea: have a list for tail factor
## tail=list(f=FALSE, f.se=NULL, sigma=NULL, F.se=NULL)
##
# 2013-02-25 Parameter risk recursive formula may have a third term per
#   Murphy and BBMW
if (! mse.method %in% c("Mack", "Independence")) stop("mse.method must be 'Mack' or 'Independence'")

Triangle <- checkTriangle(Triangle)
m <- dim(Triangle)[1]
n <- dim(Triangle)[2]

## Mack uses alpha between 0 and 2 to distinguish
## alpha = 0 straight averages
## alpha = 1 historical chain ladder age-to-age factors
## alpha = 2 ordinary regression with intercept 0

## However, in Zehnwirth & Barnett they use the notation of delta, whereby delta = 2 - alpha
## the delta is than used in a linear modelling context.
delta <- 2-alpha
alpha <- 2 - CL\$delta

# Estimate expected values and standard errors in four steps:
# 1) Squaring the Triangle: Expected values and f/F SE's from the data in the triangle
# 2) Expected values and f/F SE's from the tail factor specifications
# 3) Process Risk and Parameter Risk estimates of the squared triangle (incl tail column)
# 3) Expected values and SE's for the totals-across-origin-periods of the predicted values

## 1) Squaring the Triangle

## EXPECTED VALUES: Predict the chain ladder models
## f/F SE's
StdErr <- Mack.S.E(CL[["Models"]], FullTriangle, est.sigma = est.sigma,
weights = CL[["weights"]], alpha = alpha)

## 2) Tail
## Check for tail factor
if(is.logical(tail)){
if(tail){
tail <- tailfactor(StdErr\$f)
tail.factor <- tail\$tail.factor
StdErr\$f <- c(StdErr\$f, tail.factor = tail.factor)
}else{
#            tail.factor <- tail
# MunichChainLadder needs an 'f' vector as long as the triangle is wide
#   and adding on a harmless tail doesn't seem to cause any
#   Mack difficulties
# Default will not be named in the output
tail.factor <- 1.000
StdErr\$f <- c(StdErr\$f, tail.factor)
}
}else{
#        if(is.numeric(tail))
# Documentation says tail must be logic or numeric
tail.factor <- as.numeric(tail)
StdErr\$f <- c(StdErr\$f, tail.factor = tail.factor)
}
# Then finally, ...
if (tail.factor > 1) {
## EXPECTED VALUES
FullTriangle <- tail.E(FullTriangle, tail.factor)
## STANDARD ERRORS
## Estimate the standard error of f and F in the tail
##  If tail.se and/or tail.sigma provided, return those values
StdErr <- tail.SE(FullTriangle, StdErr, Total.SE, tail.factor,
tail.se = tail.se, tail.sigma = tail.sigma)
}

## 3) Calculate process and parameter risks of the predicted loss amounts
StdErr <- c(StdErr, MackRecursive.S.E(FullTriangle, StdErr\$f, StdErr\$f.se, StdErr\$F.se, mse.method = mse.method))

## 4) Total-across-origin-periods by development period
## EXPECTED VALUES
##   Not complicated. Not required at this time.
## STANDARD ERRORS
## Calculate process and parameter risk for the sum of the predicted loss amounts
Total.SE <- TotalMack.S.E(FullTriangle, StdErr\$f, StdErr\$f.se, StdErr\$F.se, StdErr\$FullTriangle.procrisk, mse.method = mse.method)

## Collect the output
output <- list()
output[["call"]] <-  match.call(expand.dots = FALSE)
output[["Triangle"]] <- Triangle
output[["FullTriangle"]] <- FullTriangle
output[["Models"]] <- CL[["Models"]]
output[["f"]] <- StdErr\$f
output[["f.se"]] <- StdErr\$f.se
output[["F.se"]] <- StdErr\$F.se
output[["sigma"]] <- StdErr\$sigma
output[["Mack.ProcessRisk"]]   <- StdErr\$FullTriangle.procrisk  # new dmm
output[["Mack.ParameterRisk"]] <- StdErr\$FullTriangle.paramrisk  # new dmm
output[["Mack.S.E"]] <- sqrt(StdErr\$FullTriangle.procrisk^2 + StdErr\$FullTriangle.paramrisk^2)
output[["weights"]] <- CL\$weights
output[["alpha"]] <- alpha
## total.procrisk <- apply(StdErr\$FullTriangle.procrisk, 2, function(x) sqrt(sum(x^2)))
output[["Total.Mack.S.E"]] <- Total.SE[1] # [1] removes attributes
output[["Total.ProcessRisk"]] <- attr(Total.SE, "processrisk")
output[["Total.ParameterRisk"]] <- attr(Total.SE, "paramrisk")
output[["tail"]] <- tail
return(output)
}

##############################################################################
## Calculation of the mean squared error and standard error
## mean squared error = stochastic error (process variance) + estimation error
## standard error = sqrt(mean squared error)
approx.equal <- function (x, y, tol=.Machine\$double.eps^0.5) abs(x-y)<tol
Mack.S.E <- function(MackModel, FullTriangle, est.sigma="log-linear", weights, alpha) {
n <- ncol(FullTriangle)
m <- nrow(FullTriangle)
f <- rep(1, n - 1)
f.se <- rep(0, n - 1)
sigma <- rep(0, n - 1)

## Extract estimated slopes, std. error and sigmas
## 2015-10-9 Replace lm's warning with more appropriate message
smmry <- suppressWarnings(lapply(MackModel, summary))
f <- sapply(smmry, function(x) x\$coef["x","Estimate"])
f.se <- sapply(smmry, function(x) x\$coef["x","Std. Error"])
sigma <- sapply(smmry, function(x) x\$sigma)
df <- sapply(smmry, function(x) x\$df[2L])
tolerance <- .Machine\$double.eps
perfect.fit <- (df > 0) & (f.se < tolerance)
w <- which(perfect.fit)
if (length(w)) {
warn <- "Information: essentially no variation in development data for period(s):\n"
nms <- colnames(FullTriangle)
periods <- paste0("'", paste(nms[w], nms[w+1], sep = "-"), "'")
warn <- c(warn, paste(periods, collapse = ", "))
# Print warning message
warning(warn)
}
#
isna <- is.na(sigma)

if(est.sigma[1] %in% "log-linear"){
if (sum(!isna) == 1) {
warning(paste("Too few (1) link ratios for fitting 'loglinear' model to estimate sigma_n.",
"est.sigma will be overwritten to 'Mack'.\n",
"Mack's estimation method will be used instead."))
est.sigma <- "Mack"
}
else {
## estimate sigma[n-1] via log-linear regression
sig.model <- suppressWarnings(estimate.sigma(sigma))
sigma <- sig.model\$sigma

p.value.of.model <- tryCatch(summary(sig.model\$model)\$coefficient[2,4],
error = function(e) e)
if (inherits(p.value.of.model, "error") |
is.infinite(p.value.of.model) |
is.nan(p.value.of.model)
) {
warning(paste("'loglinear' model to estimate sigma_n doesn't appear appropriate.\n",
"est.sigma will be overwritten to 'Mack'.\n",
"Mack's estimation method will be used instead."))
est.sigma <- "Mack"
}
else
if(p.value.of.model > 0.05){
warning(paste("'loglinear' model to estimate sigma_n doesn't appear appropriate.",
"\np-value > 5.\n",
"est.sigma will be overwritten to 'Mack'.\n",
"Mack's estimation method will be used instead."))

est.sigma <- "Mack"
}
else{
f.se[isna] <- sigma[isna]/sqrt(weights[1,isna]*FullTriangle[1,isna]^alpha[isna])
}
}
}
if(est.sigma[1] %in% "Mack"){
for(i in which(isna)){   # usually i = n - 1
ratio <- (sigma[i - 1]^4/sigma[i - 2]^2) # Bug fix: sigma[i - 2]^2 could be zero
if(is.nan(ratio) | is.infinite(ratio)){
sigma[i] <- sqrt(abs(min(sigma[i - 2]^2, sigma[i - 1]^2)))
}else{
sigma[i] <- sqrt(abs(min(ratio, min(sigma[i - 2]^2, sigma[i - 1]^2))))
}
f.se[i] <- sigma[i]/sqrt(weights[1,i]*FullTriangle[1,i]^alpha[i])
}
}
if(is.numeric(est.sigma)){
# Markus, I'm not sure what your intention is in this next loop when the lengths of est.sigma and sigma differ
for(i in seq(along=est.sigma)){
l <- length(est.sigma)
# BUG?! If length(est.sigma) > n, then n-i could be negative
sigma[n-i] <- est.sigma[l-i+1]
f.se[n-i] <- sigma[n-i]/sqrt(weights[1,n-i]*FullTriangle[1,n-i]^alpha[n-i])
}
}

W <- weights
W[is.na(W)] <- 1
F.se <- t(sigma/t(sqrt(W[,-n]*t(t(FullTriangle[,-n])^alpha[-n]))))

return(list(sigma = sigma,
f = f,
f.se = f.se,
F.se = F.se)
)
}
################################################################
MackRecursive.S.E <- function(FullTriangle, f, f.se, F.se, mse.method = "Mack"){
n <- ncol(FullTriangle)
m <- nrow(FullTriangle)

FullTriangle.procrisk <- FullTriangle[, 1:n] * 0
FullTriangle.paramrisk <- FullTriangle[, 1:n] * 0

## Recursive Formula
colindex <- 1:(n-1)
for (k in colindex) {
for (i in (m-k+1):m) {
FullTriangle.procrisk[i,k+1] <- sqrt(
FullTriangle[i,k]^2*(F.se[i,k]^2)
+ FullTriangle.procrisk[i,k]^2*f[k]^2
)
FullTriangle.paramrisk[i,k+1] <- sqrt(
FullTriangle[i,k]^2*(f.se[k]^2)
+ FullTriangle.paramrisk[i,k]^2*f[k]^2
# 2013-02-25 Parameter risk recursive formula may have a third term
+ ifelse(mse.method == "Mack", 0, FullTriangle.paramrisk[i, k]^2 * (f.se[k]^2))
)
}
}

return(list(FullTriangle.procrisk=FullTriangle.procrisk,
FullTriangle.paramrisk=FullTriangle.paramrisk))
}

################################################################################
## Total reserve SE

TotalMack.S.E <- function(FullTriangle, f, f.se, F.se, FullTriangle.procrisk, mse.method = "Mack") {
# 2013-02-25
# Parameter Risk and Process Risk components of Total Risk broken out

# The current program design expects a scalar, total.seR, from this
#   function equal to the total risk of the total reserve.
#   So as not to break existing code, the Parameter Risk and
#   Process Risk vectors (one element per column) are attached as attributes

C <- FullTriangle
n <- ncol(C)
m <- nrow(C)

total.paramrisk <- numeric(n)

# Assumption is that origin years are independent, therefore
total.procrisk <- sqrt(colSums(FullTriangle.procrisk^2, na.rm = TRUE))

# For parameter risk recursion, the sum of future loss expected values
#   plus the diagonal value, by development age, comes in handy.
# Currently, code relies on the fact that nrows(triangle) >= ncols
#   Actually, function 'checkTriangle' makes sure Triangle is square
M <- sapply(1:ncol(FullTriangle), function(k) sum(FullTriangle[(m + 1 - k):m, k], na.rm = TRUE))

for(k in c(1:(n-1))){
#        total.seR[k+1] <- sqrt(total.seR[k]^2 * f[k]^2 +
#                               sum(C[c((m+1-k):m),k]^2 *
#                                   (F.se[c((m+1-k):n),k]^2),na.rm=TRUE)
#                               + sum(C[c((m+1-k):m),k],na.rm=TRUE)^2 * f.se[k]^2 )
total.paramrisk[k + 1] <- sqrt(
sum(M[k], na.rm = TRUE)^2 * f.se[k]^2 +
total.paramrisk[k]^2 * f[k]^2 +
ifelse(mse.method == "Mack", 0, total.paramrisk[k]^2 * f.se[k]^2))
}

# The scalar returned is the total risk of total reserves in the rightmost column
total.seR <- sqrt(total.procrisk^2 + total.paramrisk^2)[n]
attr(total.seR, "processrisk") <- total.procrisk
attr(total.seR, "paramrisk") <- total.paramrisk

return(total.seR)
}

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

estimate.sigma <- function(sigma){
if(!all(is.na(sigma))){
n <- length(sigma)
dev <- 1:n
my.dev <- dev[!is.na(sigma) & sigma > 0]
my.model <- lm(log(sigma[my.dev]) ~ my.dev)
sigma[is.na(sigma)] <- exp(predict(my.model, newdata=data.frame(my.dev=dev[is.na(sigma)])))
}

return(list(sigma=sigma, model=my.model))
}

########################################################################
## Estimate expected value when tail

tail.E <- function(FullTriangle, tail.factor){
n <- ncol(FullTriangle)
m <- nrow(FullTriangle)

FullTriangle <- cbind(FullTriangle, FullTriangle[,n] * tail.factor)
dimnames(FullTriangle) <- list(origin=dimnames(FullTriangle)[[1]],
dev=c(dimnames(FullTriangle)[[2]][1:n], "Inf"))
return(FullTriangle)
}

########################################################################
## Estimate standard error for tail

tail.SE <- function(FullTriangle, StdErr, Total.SE, tail.factor, tail.se = NULL, tail.sigma = NULL) {
n <- ncol(FullTriangle)
m <- nrow(FullTriangle)

## Idea: linear model for f, estimate dev for tail factor
## linear model for f.se and sigma and put dev from above in
## 2013-02-28: The tail factor is given and stored in StdEff\$f[n-1];
##    want to relate the f.se's and sigma's from link ratios estimated
##    from the interior of the triangle to the given tail.
##    The link ratios from the interior of the triangle are stored in
##    StdEff\$f[1:(n-2)].
stopifnot(n > 2) # Must have at least two columns in the original triangle
# to impute estimates for the tail-driven,
# previously c-bind-ed Ultimate column
start <- 1
.f <- StdErr\$f[start:(n - 2)]
.dev <- c(start:(n - 2))
##    mf <- lm(log(.f[.f>1]-1) ~ .dev[.f>1])
mf <- lm(log(.f-1) ~ .dev)
tail.pos <- (log(StdErr\$f[n - 1] - 1) - coef(mf)[1]) / coef(mf)[2]

if(is.null(tail.se)){
.fse <- StdErr\$f.se[start:(n-2)]
mse <- lm(log(.fse) ~ .dev)
tail.se <- exp(predict(mse, newdata=data.frame(.dev=tail.pos)))
}
StdErr\$f.se <- c(StdErr\$f.se, tail.se = tail.se)

if(is.null(tail.sigma)){
.sigma <- StdErr\$sigma[start:(n-2)]
msig <- lm(log(.sigma) ~ .dev)
tail.sigma <- exp(predict(msig, newdata = data.frame(.dev = tail.pos)))
}
StdErr\$sigma <- c(StdErr\$sigma, tail.sigma = as.numeric(tail.sigma))

## estimate the standard error of the tail factor ratios
se.F.tail <- tail.sigma / sqrt(FullTriangle[, n - 1])
StdErr\$F.se <- cbind(StdErr\$F.se, se.F.tail)

return(StdErr)
}

##############################################################################
## Summary
##
## Summarise my results
Latest <- getLatestCumulative(object\$Triangle)

ex.origin.period <- Latest!=0

Ultimate <- object[["FullTriangle"]][,ncol(object[["FullTriangle"]])]
Dev.To.Date <- Latest/Ultimate
IBNR <- Ultimate-Latest
Mack.S.E <- object[["Mack.S.E"]][,ncol(object[["Mack.S.E"]])]
CV <- Mack.S.E/(Ultimate-Latest)

ByOrigin <- data.frame(Latest, Dev.To.Date, Ultimate, IBNR, Mack.S.E, CV)
names(ByOrigin)[6]="CV(IBNR)"
ByOrigin <- ByOrigin[ex.origin.period,]

Totals <-  c(sum(Latest,na.rm=TRUE),
sum(Latest,na.rm=TRUE)/sum(Ultimate,na.rm=TRUE),
sum(Ultimate,na.rm=TRUE),
sum(IBNR,na.rm=TRUE), object[["Total.Mack.S.E"]],
object[["Total.Mack.S.E"]]/sum(IBNR,na.rm=TRUE)
)
# Totals <- c(Totals, round(x[["Total.Mack.S.E"]]/sum(res\$IBNR,na.rm=TRUE),2))
Totals <- as.data.frame(Totals)

colnames(Totals)=c("Totals")
rownames(Totals) <- c("Latest:","Dev:","Ultimate:",
"IBNR:","Mack S.E.:",
"CV(IBNR):")

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

#getLatestCumulative <- function(cumulative.tri)
#  apply(cumulative.tri, 1L, function(x) ifelse(length(w <- which(!is.na(x))) > 0L, x[tail(w, 1L)], x[1L]))

##############################################################################
## print
##

summary.x <- summary(x)
print(x\$call)
cat("\n")
print(format(summary.x\$ByOrigin, big.mark = ",", digits = 3),...)

Totals <- summary.x\$Totals
rownames(Totals)[5] <- "Mack.S.E"
Totals[1:6,] <- formatC(Totals[1:6,], big.mark=",",digits=2,format="f")
cat("\n")
print(Totals, quote=FALSE)
#invisible(x)

}

################################################################################
## plot
##
x, mfrow=NULL, title=NULL,
lattice=FALSE, which=1:6, ...){

.myResult <-  summary(x)\$ByOrigin

.FullTriangle <- x[["FullTriangle"]]
.Triangle <- x[["Triangle"]]

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))
}
if(!lattice){
if(is.null(title)) myoma <- c(0,0,0,0) else myoma <- c(0,0,2,0)

op=par(mfrow=mfrow, oma=myoma, mar=c(4.5,4.5,2,2))

plotdata <- t(as.matrix(.myResult[,c("Latest","IBNR")]))
n <- ncol(plotdata)

if(1 %in% which){

if(getRversion() < "2.9.0") { ## work around missing feature

bp <- barplot(plotdata,
legend.text=c("Latest","Forecast"),
##    args.legend=list(x="topleft"), only avilable from R version >= 2.9.0
names.arg=rownames(.myResult),
xlab="Origin period",
ylab="Amount",#paste(Currency,myUnit),
ylim=c(0, max(apply(.myResult[c("Ultimate", "Mack.S.E")],1,sum),na.rm=TRUE)))

}else{
bp <- barplot(plotdata,
legend.text=c("Latest","Forecast"),
args.legend=list(x="topleft"),
names.arg=rownames(.myResult),
xlab="Origin period",
ylab="Amount",#paste(Currency,myUnit),
ylim=c(0, max(apply(.myResult[c("Ultimate", "Mack.S.E")],1,sum),na.rm=TRUE)))
}
## require("Hmisc")
.errbar(x=bp, y=.myResult\$Ultimate,
yplus=(.myResult\$Ultimate + .myResult\$Mack.S.E),
yminus=(.myResult\$Ultimate - .myResult\$Mack.S.E),
cap=0.05,
}
if(2 %in% which){

matplot(t(.FullTriangle), type="l",
main="Chain ladder developments by origin period",
xlab="Development period", ylab="Amount", #paste(Currency, myUnit)
)
}
Residuals=residuals(x)
if(3 %in% which){
plot(standard.residuals ~ fitted.value, data=Residuals,
ylab="Standardised residuals", xlab="Fitted")
lines(lowess(Residuals\$fitted.value, Residuals\$standard.residuals), col="red")
abline(h=0, col="grey")
}
if(4 %in% which){
plot(standard.residuals ~ origin.period, data=Residuals,
ylab="Standardised residuals", xlab="Origin period")
lines(lowess(Residuals\$origin.period, Residuals\$standard.residuals), col="red")
abline(h=0, col="grey")
}
if(5 %in% which){
plot(standard.residuals ~ cal.period, data=Residuals,
ylab="Standardised residuals", xlab="Calendar period")
lines(lowess(Residuals\$cal.period, Residuals\$standard.residuals), col="red")
abline(h=0, col="grey")
}
if(6 %in% which){
plot(standard.residuals ~ dev.period, data=Residuals,
ylab="Standardised residuals", xlab="Development period")
lines(lowess(Residuals\$dev.period, Residuals\$standard.residuals), col="red")
abline(h=0, col="grey")
}
title( title , outer=TRUE)
par(op)

}else{

## require(grid)
## Set legend
fl <-
grid.layout(nrow = 2, ncol = 4,
heights = unit(rep(1, 2), "lines"),
widths =
unit(c(2, 1, 2, 1),
c("cm", "strwidth", "cm",
"strwidth"),
data = list(NULL, "Chain ladder dev.", NULL,
"Mack's S.E.")))

foo <- frameGrob(layout = fl)

foo <- placeGrob(foo,
linesGrob(c(0.2, 0.8), c(.5, .5),
gp = gpar(col=1, lty=1)),
row = 1, col = 1)
foo <- placeGrob(foo,
linesGrob(c(0.2, 0.8), c(.5, .5),
gp = gpar(col=1, lty=3)),
row = 1, col = 3)
foo <- placeGrob(foo,
row = 1, col = 2)
foo <- placeGrob(foo,
textGrob(label = "Mack's S.E."),
row = 1, col = 4)

long <- expand.grid(origin=as.numeric(dimnames(.FullTriangle)\$origin),
dev=as.numeric(dimnames(.FullTriangle)\$dev))
long\$value <- as.vector(.FullTriangle)
long\$valuePlusMack.S.E <-  long\$value + as.vector(x\$Mack.S.E)
long\$valueMinusMack.S.E <- long\$value - as.vector(x\$Mack.S.E)
sublong <- long[!is.na(long\$value),]
xyplot(valuePlusMack.S.E + valueMinusMack.S.E + value ~ dev |
factor(origin), data=sublong, t="l", lty=c(3,3,1), as.table=TRUE,
main="Chain ladder developments by origin period",
xlab="Development period",
ylab="Amount",col=1,
legend = list(top = list(fun = foo)),...)
}
}
################################################################################
## residuals
##
m <- nrow(object[["Triangle"]])
n <- ncol(object[["Triangle"]])
myresiduals <- unlist(lapply(object[["Models"]], resid,...))

standard.residuals <- rep(NA, length(myresiduals))

## Identify if we have standardized residuals for all items,
## e.g. this is not the case if some weights have been set to zero
x=lapply(lapply(object\$Model, resid), function(x) as.numeric(names(x)))
y=lapply(lapply(object\$Model, rstandard), function(x) as.numeric(names(x)))
names(y)=1:length(y)
names(x)=1:length(x)
rst <- unlist(lapply(names(x), function(z) x[[z]] %in% y[[z]]))
## rst holds the indices with standardised residuals for weights > 0

standard.residuals[rst] <- unlist(lapply(object[["Models"]], rstandard,...))
fitted.value <- unlist(lapply(object[["Models"]], fitted))
origin.period <- as.numeric(unlist(lapply(lapply(object[["Models"]],residuals),names)))
dev.period <-rep(1:(n-1), sapply(lapply(object[["Models"]],residuals),length))
cal.period <- origin.period + dev.period - 1

myResiduals <- data.frame(origin.period,
dev.period,
cal.period,
residuals=myresiduals,
standard.residuals,
fitted.value)

return(na.omit(myResiduals))
}

Mack <- function(triangle, weights=1, alpha=1){
## triangle: cumulative claims triangle
## weights : weights
## alpha=0 : gives straight average of the chain ladder age-to-age factors
## alpha=1 : gives the historical chain ladder age-to-age factors
## alpha=2 : is the result of an ordinary regression of C_{i,k+1} against C_{i,k} with intercept 0.
weights <- checkWeights(weights, triangle)

m <- nrow(triangle)
n <- ncol(triangle)

## Individual chain ladder age-to-age factors:
F <- triangle[,-1]/triangle[,-n]
#    if(is.numeric(weights)){
#        weights2 <- triangle
#        weights2[!is.na(weights2)] <- weights
#        weights <- weights2
#    }

wCa <- weights[,-n] * triangle[,-n]^alpha
## Note NA^0=1, hence set all cells which are originally NA back to NA
wCa[is.na(triangle[,-n])] <- NA

wCaF <- wCa*F
## Set diagonal of wCa to NA
wCa[row(wCa)==n+1-col(wCa)] <- NA
f <- apply( wCaF, 2, sum, na.rm=TRUE)/apply( wCa, 2, sum, na.rm=TRUE)

## Get full triangle
avDFs <- c(f,1)
dim(avDFs) <- c(1,n,1)
ultDFs <- getUltDFs(avDFs)

Latest <- getLatestCumulative(triangle)
ults <- getUltimates(Latest, ultDFs)
FullTriangle <- getExpected(ults, 1/ultDFs)
dim(FullTriangle)=c(n,n)

## Estimate standard errors

k <- c(1:(n-2))
wCa_F_minus_f <- weights[,k] * triangle[,k]^alpha * t(t(F)-f)[,k]^2

wCa_F_minus_f <-  triangle[,k]^alpha * t(t(F)-f)[,k]^2

sigma <- sqrt( 1/c(n-k-1) * apply(wCa_F_minus_f, 2, sum, na.rm=TRUE) )
# Mack approximation for sigma_{n-1}
sigma_n1 <- sqrt( abs(min(sigma[n-2]^4/sigma[n-3]^2, min(sigma[n-3]^2, sigma[n-2]^2)) ))
sigma <- c(sigma, sigma_n1)

# Estimate f.se
wCa2 <-  weights[,c(1:(n-1))] * triangle[,c(1:(n-1))]^alpha
## Note NA^0=1, hence set all cells which are originally NA back to NA
wCa2[is.na(triangle[,-n])] <- NA

## Set diagonal of wCa to na
wCa2[row(wCa2)==n+1-col(wCa2)] <- NA
f.se <- sigma / sqrt( apply(wCa2, 2, sum, na.rm=TRUE))

# Estimate F.se
W <- weights
W[is.na(W)] <- 1
F.se <- t(sigma/sqrt(t(W[,c(1:(n-1))]*FullTriangle[,c(1:(n-1))]^alpha)))
#    F.se <- t(sigma/sqrt(t(FullTriangle[,c(1:(n-1))]^alpha)))

F.se[is.na(FullTriangle[,c(1:(n-1))])] <- NA

FullTriangle.se <- FullTriangle * 0
## Recursive Formula
rowindex <- 2:m
if(m>n)
rowindex <- c((m-n+1):m)
for(i in rowindex){
for(k in c((n+1-i):(n-1))){
if(k>0)
FullTriangle.se[i,k+1] = sqrt(
FullTriangle[i,k]^2*(F.se[i,k]^2+f.se[k]^2) #
+ FullTriangle.se[i,k]^2*f[k]^2
)
}
}

output <- list()
output[["call"]] <-  match.call(expand.dots = FALSE)
output[["Triangle"]] <- triangle
output[["Latest"]] <- Latest
output[["FullTriangle"]] <- FullTriangle
output[["f"]] <- f
output[["F"]] <- F
output[["f.se"]] <- f.se
output[["F.se"]] <- F.se
output[["sigma"]] <-sigma
output[["Mack.S.E"]] <- FullTriangle.se
class(output) <- c("Mack","list")
return(output)
}

print.Mack <- function(x,...){
n <- ncol(x\$Triangle)
df <- data.frame(Latest=x\$Latest, f=c(NA, x\$f), f.se=c(NA, x\$f.se),
Ultimate=x\$FullTriangle[,n])
df\$IBNR <- df\$Ultimate-df\$Latest
df\$Mack.S.E=x\$Mack.S.E[,n]
df\$CV=df\$Mack.S.E/df\$IBNR
print(df)
}

.errbar <- function (
x, y, yplus, yminus, cap = 0.015,
main = NULL, sub = NULL,
xlab = as.character(substitute(x)),
ylab = if (is.factor(x) || is.character(x)) "" else as.character(substitute(y)),
add = FALSE, lty = 1, type = "p", ylim = NULL, lwd = 1, pch = 16,
errbar.col = par("fg"), Type = rep(1, length(y)), ...)
{
# Based on code by Frank Harrell, Hmisc package, licence: GPL >= 2
if (is.null(ylim))
ylim <- range(y[Type == 1], yplus[Type == 1], yminus[Type ==
1], na.rm = TRUE)
if (is.factor(x) || is.character(x)) {
x <- as.character(x)
n <- length(x)
t1 <- Type == 1
t2 <- Type == 2
n1 <- sum(t1)
n2 <- sum(t2)
omai <- par("mai")
mai <- omai
mai[2] <- max(strwidth(x, "inches")) + 0.25
par(mai = mai)
on.exit(par(mai = omai))
plot(NA, NA, xlab = ylab, ylab = "", xlim = ylim, ylim = c(1,
n + 1), axes = FALSE, ...)
axis(1)
w <- if (any(t2))
n1 + (1:n2) + 1
else numeric(0)
axis(2, at = c(seq.int(length.out = n1), w), labels = c(x[t1],
x[t2]), las = 1, adj = 1)
points(y[t1], seq.int(length.out = n1), pch = pch, type = type,
...)
segments(yplus[t1], seq.int(length.out = n1), yminus[t1],
seq.int(length.out = n1), lwd = lwd, lty = lty, col = errbar.col)
if (any(Type == 2)) {
abline(h = n1 + 1, lty = 2, ...)
offset <- mean(y[t1]) - mean(y[t2])
if (min(yminus[t2]) < 0 & max(yplus[t2]) > 0)
lines(c(0, 0) + offset, c(n1 + 1, par("usr")[4]),
lty = 2, ...)
points(y[t2] + offset, w, pch = pch, type = type,
...)
segments(yminus[t2] + offset, w, yplus[t2] + offset,
w, lwd = lwd, lty = lty, col = errbar.col)
at <- pretty(range(y[t2], yplus[t2], yminus[t2]))
axis(side = 3, at = at + offset, labels = format(round(at,
6)))
}
return(invisible())
}
points(x, y, pch = pch, type = type, ...)
else plot(x, y, ylim = ylim, xlab = xlab, ylab = ylab, pch = pch,
type = type, ...)
xcoord <- par()\$usr[1:2]
smidge <- cap * (xcoord[2] - xcoord[1])/2
segments(x, yminus, x, yplus, lty = lty, lwd = lwd, col = errbar.col)
if (par()\$xlog) {
xstart <- x * 10^(-smidge)
xend <- x * 10^(smidge)
}
else {
xstart <- x - smidge
xend <- x + smidge
}
segments(xstart, yminus, xend, yminus, lwd = lwd, lty = lty,
col = errbar.col)
segments(xstart, yplus, xend, yplus, lwd = lwd, lty = lty,
col = errbar.col)
return(invisible())
}
```
edalmoro/ChainLadderQuantileV1 documentation built on Oct. 1, 2018, 12:23 a.m.