GIRO2015.R

library(lattice)
library(data.table)
library(arm)
library(ChainLadder)
dat <- data.table( # Page D5.17
  originf=factor(rep(0:6, each=7)),  
  devf=factor(rep(0:6, 7)),
  inc.value= c(3511,  3215,  2266,  1712, 1059, 587, 340,
               4001,  3702,  2278, 1180, 956, 629,  NA,
               4355,  3932,  1946, 1522, 1238,  NA,  NA,
               4295,  3455,  2023, 1320,   NA,  NA,  NA,
               4150,  3747,  2320,   NA,   NA,  NA,  NA,
               5102,  4548,    NA,   NA,   NA,  NA,  NA,
               6283,    NA,    NA,   NA,   NA,  NA,  NA))


dat <- dat[order(devf), cum.value:=cumsum(inc.value), by=originf]



xyplot(cum.value + 
         inc.value +
         log(inc.value) ~ devf , groups=originf, data=dat, 
       t="b", par.settings = simpleTheme(pch = 16),
       auto.key = list(space="right",
                       title="Origin\nperiod", cex.title=1,
                       points=FALSE, lines=TRUE, type="b"),
       xlab="Development period", ylab="Amount",
       main="Claims development by origin period",
       layout=c(3,1),
       scales="free")


(model <- lm(log(inc.value) ~ originf + devf, data=na.omit(dat)))


display(model)


par(mfrow=c(1,4))
plot(model)

summary(model)
sig <- summary(model)$sigma

pred <- predict(model, newdata=dat, se.fit = TRUE)

dat <- dat[,
    ':='(y = log(inc.value),
         yhat = pred$fit,
         phat = exp(pred$fit + 0.5*(sig^2+pred$se.fit^2)),
         origin=as.numeric(as.character(originf)),
         dev=as.numeric(as.character(devf))
    )]
dat <- dat[, ':='(cal=origin+dev,
           phat.se = phat * sqrt(exp(pred$residual.scale^2 + pred$se.fit^2)-1)),]

# Reserve
sum(dat[cal > max(origin), phat])

xyplot( phat + I(phat + phat.se) + I(phat-phat.se) + inc.value ~ devf | originf, 
        data=dat, t="b" , lwd=1.5, 
        auto.key = list(points=FALSE, lines=TRUE), as.table=TRUE,
        layout=c(2,4))

confInt <- as.data.frame(
  predict(model, newdata=dat, 
          interval = "confidence"))
names(confInt) <- c('yhat', 'conf.lwr', 'conf.upr')

predInt <- as.data.frame(
  predict(model, newdata=dat,
          interval = "prediction"))
names(predInt) <- c('yhat', 'pred.lwr', 'pred.upr')

dat <- cbind(dat, confInt[,-1], predInt[, -1])

key <- list(
  rep=FALSE, 
  lines=list(col=c("#00526D", "blue", "black"), type=c("p","l", "l"), pch=1),
  text=list(lab=c("Observation","Mean estimate", "95% Prediction interval")),
  rectangles = list(col=adjustcolor("yellow", alpha.f=0.5), border="grey"),
  text=list(lab="95% Confidence interval"))

P <- xyplot( conf.lwr + conf.upr + yhat + 
          y  +  pred.lwr + pred.upr ~ devf | originf,  dat=dat, 
        lwd=1.5, t="l", ylab="log(Claims)", xlab="Dev. Period",
        key=key,as.table=TRUE, layout=c(4,2),
        panel=function(x, y){
          n <- length(x)
          k <- n/3
          upper <- y[(k/2+1):k]
          lower <- y[1:(k/2)]
          x <- x[1:(k/2)]
          panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
                        col = adjustcolor("yellow", alpha.f = 0.5), 
                        border = "grey")
          panel.lines(x, y[(k+1):(k+n/6)], col="blue")
          panel.points(x, y[(k*3/2+1):(k*4/2)], lwd=2, col="#00526D")
          panel.lines(x, y[(k*4/2+1):(k*5/2)], col="black")
          panel.lines(x, y[(k*5/2+1):(k*6/2)], col="black")
        })
P

setwd("G:/Analysis/Outputs/Presentations/GIRO 2015")
png(filename = "ModelIntervals.png",units="px", 
    width=1600, height=1000, res=200)
P
dev.off()



dat <- dat[order(dev), ':='(cum.data=cumsum(inc.value),
                            cum.phat=cumsum(phat)), by=originf]



xyplot( cum.phat + cum.data ~ devf | originf, data=dat, 
       t="b", par.settings = simpleTheme(pch = 16),
       auto.key = list(space="top", columns=2,
                       points=FALSE, lines=TRUE, type="b"),
       xlab="Development period", ylab="Amount",
       main="Claims development by origin period",
       as.table=TRUE,
       layout=c(4,2))

fdm <- model.matrix(~ originf + devf, data=dat[cal > max(origin)])



varcovar <- fdm %*% vcov(model) %*% t(fdm)
(Total.SE <- sqrt(t(dat[cal>max(origin), phat]) %*% (exp(varcovar)-1) %*% 
                    dat[cal>max(origin), phat])) 


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

dat <-  dat[, ':='(d = ifelse(dev < 1, 1, 0),
                   s = ifelse(dev < 1, 0, dev))]        

# Page D5.36
ClaimsVolume <- data.table(
  origin=0:6, 
  volume.index=c(1.43, 1.45, 1.52, 1.35, 1.29, 1.47, 1.91),
  key="origin")
EarningIndex <- data.table(
  cal=0:6, 
  earning.index=c(1.55, 1.41, 1.3, 1.23, 1.13, 1.05, 1),
  key="cal")
setkey(dat, origin, cal)
dat <- dat[ClaimsVolume][EarningIndex]
# Normalise data for volume and earnings
dat <- dat[, logvalue.ind.inf:=log(inc.value/volume.index*earning.index)]
par(mfrow=c(1,1))
with(dat, interaction.plot(dev, origin, logvalue.ind.inf))
points(1+dat$dev, dat$logvalue.ind.inf, pch=16, cex=0.8)

display(model2 <- lm(logvalue.ind.inf ~  d + s, data=dat))

pred <- predict(model2, newdata=dat, se.fit=TRUE)

dat <- dat[ , ':='(
  phat2=exp(pred$fit + 0.5*(pred$se.fit^2 + pred$residual.scale^2))*volume.index/earning.index
)]

sum(dat[ cal > max(origin), phat2])

dat <- dat[order(dev), cum.phat2:=cumsum(phat2), by=originf]


xyplot( cum.phat  + cum.data + cum.phat2 ~ devf | originf, data=dat, 
        t="b", par.settings = simpleTheme(pch = 16),
        auto.key = list(space="top", columns=2,
                        text=c("prediction (full model)","data","prediction (reduced model)"),
                        points=FALSE, lines=TRUE, type="b"),
        xlab="Development period", ylab="Amount",
        main="Claims development by origin period",
        as.table=TRUE,
        layout=c(4,2))









display(model2 <- lm(logvalue.ind.inf ~  d + s, data=dat))


library(brms)
rtools <- "C:\\Rtools\\bin"
gcc <- "C:\\Rtools\\gcc-4.6.3\\bin"
path <- strsplit(Sys.getenv("PATH"), ";")[[1]]
new_path <- c(rtools, gcc, path)
new_path <- new_path[!duplicated(tolower(new_path))]
Sys.setenv(PATH = paste(new_path, collapse = ";"))

bmodel <- brm(log(inc.value) ~ originf + devf + 0, data=na.omit(dat),
                family = "gaussian")

bmodel$model

x <- predict(bmodel, newdata=dat[, c('originf', 'devf')])
x <- x[is.na(dat$inc.value),]

names(x)

sum(exp(x[,1]+0.5*x[,2]^2))

X <- apply(x, 1, function(y) rlnorm(10000, y[1], y[2]))

ultimate <- apply(X, 1, sum)
summary(ultimate)
sd(ultimate)
plot(density(ultimate))

paiddev <- apply(X, 1, cumsum)
dim(paiddev)
paidq <- apply(paiddev, 1, quantile, probs=c(0.5, 0.25, 0.5, 0.75, 0.95))

matplot(t(paidq), t="l")

(model2 <- lm(log(inc.value) ~ originf + d + s, data=na.omit(dat)))
display(model2)
plot(model2)

tail.years <- 6
fdat <- data.frame(
  originf=factor(rep(0:6, 7+tail.years)),
  dev=rep(0:(6+tail.years), each=7)
)
fdat <- with(fdat, 
             data.frame(fdat,
                        d = ifelse(dev < 1, 1, 0),
                        s = ifelse(dev < 1, 0, dev)))        


yhat <- predict(model2, newdata=fdat, se.fit = TRUE)

phat <- exp(yhat$fit + 0.5*(yhat$se.fit^2 + yhat$residual.scale^2))

sum(phat[!is.na(is.na(dat$inc.value))])



library(ChainLadder)
glmReserve(UKMotor, var.power = NULL)

MackChainLadder(UKMotor, est.sigma="Mack")