#This code was obtained from https://www.kaggle.com/code/floser/glm-neural-nets-and-xgboost-for-insurance-pricing, with some adjustment for the purposes of this assignment.

dat <- read.csv("/content/MotorcycleData.csv",header=TRUE)

dat$Gender <- factor(dat$Gender)
dat$ClaimCount <- pmin(dat$ClaimCount,4)
dat$ClaimCost <- pmin(dat$ClaimCost,100000000)
dat$Exposure <- pmin(dat$Exposure,1)

# inspect data set
dim(dat)  # count rows and colums
head(dat) # show first few policies
tail(dat) # show last few policies

#library(tidyverse) # work quietly to use %>% and make code tidy:
suppressMessages(library(tidyverse))

#-	GeographicZone: we choose a continuous feature component for the GeographicZone code and we therefore map {A,...,F} → {1,...,6};
#-	EngineClass: we build 6 categorical classes by merging vehicle power groups bigger and equal to 9 (totally 6 labels);
#-	VehicleAge: we build 3 categorical classes [0,1), [1,10], (10,∞);
#-	Age: we build 7 categorical classes [18,21), [21,26), [26,31), [31,41), [41,51), [51,71), [71, ∞);
#-	BonusClass: continuous variable;
#-	Gender: binary feature component;
#-	PolicyDuration: Duration is chosen as continuous feature component;

dat1 <- dat
dat1$GeographicZoneGLM <- as.integer(as.factor(dat1$GeographicZone))
dat1$EngineClassGLM <- as.factor(pmin(dat1$EngineClass,9))
VehicleAgeGLM <- cbind(c(0:110), c(1, rep(2,10), rep(3,100)))
dat1$VehicleAgeGLM <- as.factor(VehicleAgeGLM[dat1$VehicleAge+1,2])
dat1[,"VehicleAgeGLM"] <-relevel(dat1[,"VehicleAgeGLM"], ref="2")
AgeGLM <- cbind(c(18:100), c(rep(1,21-18), rep(2,26-21), rep(3,31-26), rep(4,41-31), rep(5,51-41), rep(6,71-51), rep(7,101-71)))
dat1$AgeGLM <- as.factor(AgeGLM[dat1$Age-17,2])
dat1[,"AgeGLM"] <-relevel(dat1[,"AgeGLM"], ref="5")
dat1$BonusClassGLM <- as.integer(pmin(dat1$BonusClass, 150))
dat1$PolicyDurationGLM <- as.numeric(dat1$PolicyDuration)
dat1$Severity <- dat1$ClaimCost/dat1$ClaimCount
head(dat1) # show first few policies and generated features

# Modified function PDX: Weighted Poisson Deviance
PDW <- function(pred, obs, ex = rep(1, length(obs))) {
  200 * sum( ex * ( pred - obs  + log((obs / pred) ^ (obs )))) / sum(ex)
}

# Function PD2: Print Poisson Deviance learn/test
PDW2 <- function(txt, l.c, l.x, l.e, t.c, t.x, t.e) {
    sprintf("%s, Learn/Test: %.2f%% / %.2f%%", txt, PDW(l.c, l.x, l.e), PDW(t.c, t.x, t.e))
}

# Function CF2: Print claim frequency
CF2 <- function(txt, l.c, l.x, t.c, t.x) {
    sprintf("%s: %.2f%% / %.2f%%", txt, sum(l.c)/sum(l.x)*100, sum(t.c)/sum(t.x)*100)
}

# Function Benchmark.GLM2: Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
Benchmark.GLM2 <- function(txt, pred) {
  index <- ((PDW(pred, test$ClaimCount, test$Exposure) - PDW(test$fit.cf, test$ClaimCount, test$Exposure)) / (PDW(test$fitGLM2, test$ClaimCount, test$Exposure) - PDW(test$fit.cf, test$ClaimCount, test$Exposure))) * 100
  sprintf("GLM2-Improvement-Index (PD test) of %s: %.1f%%", txt, index)
}

#In this case, the learning sample randomly takes 80% and the test sample 20% of the data sets.

# for later use we create five 20%-subsamples ("folds") and take the last fold as the holdout data set
k <- 5
set.seed(42)
fold <- sample(1:k, nrow(dat1), replace = TRUE)
dat1$fold <- fold
learn <- dat1[dat1$fold != 5,]    # 80%
test <- dat1[ dat1$fold == 5,]    # 20%
CF2("Claim Frequency (Actual) Learn/Test", learn$ClaimCount,learn$Exposure, test$ClaimCount,test$Exposure)
learnsev <- learn[learn$ClaimCount>0,]
testsev <- test[test$ClaimCount>0,]
dat1sev <- dat1[dat1$ClaimCount>0,]
head(learnsev, n = 30)
head(testsev, n = 30)

#The model "INT" consists solely of the intercept. It should give an impression of how much or little of the overall variation can be explained by classical and machine learning models.

# Model INT "predictions"
cf <- sum(learn$ClaimCount)/sum(learn$Exposure) # claim frequency
learn$fit.cf <- cf*learn$Exposure
test$fit.cf <- cf*test$Exposure
# Print Poisson Deviance (weighted)
PDW2("Poisson Deviance INT", learn$fit.cf,learn$ClaimCount,learn$Exposure, test$fit.cf,test$ClaimCount,test$Exposure)




# Print claim frequency actual vs predicted (Intercept model)
CF2("Claim Frequency INT, Test-Sample, Actual/Predicted", test$ClaimCount,test$Exposure, test$fit.cf,test$Exposure)

# Basic GLM
f.glm2 <- ClaimCount ~ Age + Gender + GeographicZoneGLM + EngineClassGLM + VehicleAge + BonusClassGLM + PolicyDurationGLM
d.glm2 <- glm(f.glm2, data=learn, offset=log(Exposure), family=poisson())

summary(d.glm2) # print fit statistics and parameter list

# make predictions
learn$fitGLM2 <- fitted(d.glm2)
test$fitGLM2 <- predict(d.glm2, newdata=test, type="response")
dat$fitGLM2 <- predict(d.glm2, newdata=dat1, type="response") # for the CANN

# Print claim frequency actual vs predicted
CF2("Claim Frequency GLM2, Test-Sample, Actual/Predicted", test$ClaimCount,test$Exposure, test$fitGLM2,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance GLM2", learn$fitGLM2,learn$ClaimCount,learn$Exposure, test$fitGLM2,test$ClaimCount,test$Exposure)



#LASSO
# data preparation
install.packages("glmnet")
suppressMessages(library(glmnet))

f.lasso <- as.formula(ClaimCount ~ (AgeGLM + Gender + EngineClassGLM + VehicleAgeGLM + BonusClassGLM + PolicyDurationGLM + GeographicZoneGLM)^2)
x <- sparse.model.matrix(f.lasso, dat1)[, -1]

learn.glmnet <- x[dat1$fold != 5,]    # 80%
test.glmnet <- x[dat1$fold == 5,]    # 20%

# cross-validation
# claim number and offset
{t1 <- proc.time()
   cv.lasso <- cv.glmnet(x = learn.glmnet, y=learn$ClaimCount, alpha = 1, offset = log(learn$Exposure), family = "poisson", nfolds = 5)
   (proc.time()-t1)}

plot(cv.lasso)

cv.lasso$lambda.min

# Show (a few) non-zero coefficients
df_lasso_coef <- as.data.frame(as.matrix(coef(cv.lasso, s = "lambda.min")))
df_lasso_coef_filtered <- subset(df_lasso_coef, lambda.min != 0)
dim(df_lasso_coef_filtered)
head(df_lasso_coef_filtered)

test$fitLASSO <-  predict(cv.lasso, newx = test.glmnet, newoffset = log(test$Exposure), type = "response", s = "lambda.min")[,1]
learn$fitLASSO <- predict(cv.lasso, newx = learn.glmnet, newoffset = log(learn$Exposure), type = "response" , s = "lambda.min")[,1]

# Print claim frequency actual vs predicted
CF2('Claim Frequency LASSO, Test-Sample, Actual/Predicted', test$ClaimCount,test$Exposure, test$fitLASSO,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance LASSO",
learn$fitLASSO,learn$ClaimCount,learn$Exposure,
test$fitLASSO,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
Benchmark.GLM2("LASSO", test$fitLASSO)


#Compare predictions between models
head(learn, n = 30)



#Ridge (Poisson) Regression
{t1 <- proc.time()
cvRIDGE <- cv.glmnet(x=learn.glmnet, y=learn$ClaimCount/learn$Exposure, alpha = 0, weights = learn$Exposure, family = "poisson", nfolds = 5)
 (proc.time()-t1)}

plot(cvRIDGE)

cvRIDGE$lambda.min

test$fitRIDGE <-  exp(predict(cvRIDGE, newx = test.glmnet, s = "lambda.min"))[,1] * test$Exposure
learn$fitRIDGE <- exp(predict(cvRIDGE, newx = learn.glmnet, s = "lambda.min"))[,1] * learn$Exposure

# Print claim frequency actual vs predicted
CF2('Claim Frequency RIDGE, Test-Sample, Actual/Predicted', test$ClaimCount,test$Exposure, test$fitRIDGE,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance RIDGE", learn$fitRIDGE,learn$ClaimCount,learn$Exposure, test$fitRIDGE,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
Benchmark.GLM2("RIDGE", test$fitRIDGE)


#GAM
suppressMessages(library(mgcv))

# bam(): Generalized additive models for very large datasets, see https://cran.r-project.org/web/packages/mgcv/mgcv.pdf
#        "The advantage of bam is much lower memory footprint than gam, but it can also be much faster, for large datasets."

{t1 <- proc.time()
  d.gam3 <- bam(ClaimCount ~ s(VehicleAge) + s(Age) + BonusClassGLM + EngineClassGLM + Gender + GeographicZoneGLM + PolicyDurationGLM + offset(log(Exposure)),
                ,data=learn, scale=-1, family=poisson)
(proc.time()-t1)}
summary(d.gam3) # print fit statistics


# make predictions
learn$fitGAM3 <- fitted(d.gam3)
test$fitGAM3 <- predict(d.gam3, newdata=test, type="response")
# Print claim frequency actual vs predicted
CF2("Claim Frequency GLM2, Test-Sample, Actual/Predicted", test$ClaimCount,test$Exposure, test$fitGAM3,test$Exposure)

# Print Poisson Deviance
PDW2('Poisson Deviance GAM3', learn$fitGAM3,learn$ClaimCount,learn$Exposure, test$fitGAM3,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
Benchmark.GLM2("GAM3", test$fitGAM3)


#Compare predictions between models
head(learn, n = 30)





#XG Boost
#Data Preparation

num_features <- c('BonusClass','Age','EngineClass','VehicleAge','PolicyDuration')
cat_features <- c('GeographicZone','Gender')
features <- c(num_features, cat_features)

# Feature pre-processing for XGBoost and lightGBM

df_feat <- dat[,c(features)]
df_feat$GeographicZone <- as.integer(as.factor((df_feat$GeographicZone)))
print(paste("Features before one-hot-encoding:",dim(df_feat)[2]))

# one-hot encoding for VehBrand
df_feat <- as.data.frame(model.matrix( ~ 0 +. ,data = df_feat))
print(paste("Features after one-hot-encoding: ",dim(df_feat)[2]))

# split data into learn and test
X_learn <- df_feat[dat1$fold != 5,]    # 80%
X_test  <- df_feat[dat1$fold == 5,]    # 20%


# create monotonic constraints vector: increasing: +1, decreasing: -1, unconstraint: 0
mtc <- c(rep( 0, length(X_learn)))
head(X_learn,2) # check feature names
mtc[1] <- 1 # BonusClass increasing
print(mtc)


install.packages("xgboost")
suppressMessages(library(xgboost))
# construct xgb.DMatrices (these are the internal data structures used by XGBoost during training)
dtrain <- xgb.DMatrix(data = data.matrix(X_learn), label = learn$ClaimCount/learn$Exposure, weight=learn$Exposure)
# make label wrong (shouldn't change anything): change to fold-number
dtest <- xgb.DMatrix(data = data.matrix(X_test), label = test$fold/test$Exposure, weight=test$Exposure)




#XGBoost without constraints
# train model (using tuned hyperparameters, see corresponding Python notebook)
{t1 <- proc.time()
  fit.xgb0 <- xgb.train(data=dtrain, objective='count:poisson',
                       nrounds = 500, max_depth = 5, eta = 0.05, tree_method = "hist")
(proc.time()-t1)}

## make predictions
learn$fitXGB0 <- predict(fit.xgb0 , newdata = dtrain) * learn$Exposure
test$fitXGB0 <- predict(fit.xgb0, newdata = dtest) * test$Exposure

# Print claim frequency actual vs predicted
CF2('Claim Frequency XGB, Test-Sample, Actual/Predicted', test$ClaimCount,test$Exposure, test$fitXGB0,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance XGB", learn$fitXGB0,learn$ClaimCount,learn$Exposure, test$fitXGB0,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
test$fit.cf <- test$Exposure * sum(learn$ClaimCount)/sum(learn$Exposure) # (recalculate INT-Model)
Benchmark.GLM2("XGB unconstraint", test$fitXGB0)


# XGBoost with monotonically increasing condition with BonusClass to account for a tariff structure (as an example)

# train model (using tuned hyperparameters, see corresponding Python notebook)
{t1 <- proc.time()
  fit.xgb <- xgb.train(data=dtrain, objective='count:poisson', monotone_constraints = mtc,
                       nrounds = 500, max_depth = 5, eta = 0.05, tree_method = "hist")
(proc.time()-t1)}

## make predictions
learn$fitXGB <- predict(fit.xgb , newdata = dtrain) * learn$Exposure
test$fitXGB <- predict(fit.xgb, newdata = dtest) * test$Exposure

# Print claim frequency actual vs predicted
CF2('Claim Frequency XGB, Test-Sample, Actual/Predicted', test$ClaimCount,test$Exposure, test$fitXGB,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance XGB", learn$fitXGB,learn$ClaimCount,learn$Exposure, test$fitXGB,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
test$fit.cf <- test$Exposure * sum(learn$ClaimCount)/sum(learn$Exposure) # (recalculate INT-Model)
Benchmark.GLM2("XGB", test$fitXGB)



#Light GBM is a fast and powerful gradient tree boosting machine

install.packages("R6", repos = "http://cran.r-project.org")
install.packages("lightgbm")
suppressMessages(library(lightgbm))

# LGB matrices (frequency-weighted-version)
# https://towardsdatascience.com/how-to-handle-the-exposure-offset-with-boosted-trees-fd09cc946837
llearn <- lgb.Dataset(data = data.matrix(X_learn), label = learn$ClaimCount/learn$Exposure, weight=learn$Exposure)
ltest <- lgb.Dataset(data = data.matrix(X_test), label = test$ClaimCount/test$Exposure, weight=test$Exposure)


params <- list(
  objective = "poisson"
  , metric = "poisson"
  , learning_rate = 0.1
  , monotone_constraints = mtc
  , num_leaves=31
  , n_estimators=100
)

# train lgb-model (using default hyperparameters, see corresponding Python notebook)
valids <- list(test = ltest)
{t2 <- proc.time()
LGB <- lgb.train(
  params = params
  , data = llearn
  , valids = valids
)
(proc.time()-t2)}

## make predictions
learn$fitLGB <- predict(LGB, data.matrix(X_learn)) * learn$Exposure
test$fitLGB <- predict(LGB, data.matrix(X_test)) * test$Exposure

# Print claim frequency actual vs predicted
CF2('Claim Frequency lightGBM, Test-Sample, Actual/Predicted', test$ClaimCount,test$Exposure, test$fitLGB,test$Exposure)

# Print Poisson Deviance
PDW2("Poisson Deviance lightGBM", learn$fitLGB,learn$ClaimCount,learn$Exposure, test$fitLGB,test$ClaimCount,test$Exposure)

# Improvement in Poisson Deviance on test set compared to GLM2-INT-Improvement
test$fit.cf <- test$Exposure * sum(learn$ClaimCount)/sum(learn$Exposure) # (recalculate INT-Model)
Benchmark.GLM2("LGB", test$fitLGB)


#Explaining lightGBMs predictions with TreeSHAP
#SHAP (SHapley Additive exPlanations) is an explainable AI method that helps interpret the output of any machine learning model. It fairly calculates the contribution of each feature to the final prediction. SHAP values range from negative to positive, with higher values indicating a greater impact on the outcome. It is based on the concept of Shapley values, a cooperative game theory that assigns a value to each player's contribution to a group's outcome. SHAP can be used to explain the model's decision for a single prediction instance or to provide a global overview of the model's behavior.

#For more information and an implementation in R, we recommend the tutorial "SHAP for Actuaries: Explain any Model" of M. Mayer, D. Meier, and M.V. Wüthrich (2023), https://github.com/actuarial-data-science/Tutorials/tree/master/14%20-%20SHAP.

# SHAP
install.packages(c("shapviz","kernelshap", "lightgbm"))
library(shapviz)
library(kernelshap)

# Generate explainer object with TreeSHAP to explain the LightGBM model
{t2 <- proc.time()
 shap_lgb_test <- shapviz(LGB, X_pred = data.matrix(X_test))
(proc.time()-t2)}

#SHAP importance: summary plot
sv_importance(shap_lgb_test, kind = "bee")

# SHAP feature importance
sv_importance(shap_lgb_test, show_numbers = TRUE)

# SHAP dependence plots (note the different y-scales)
options(repr.plot.width=20, repr.plot.height = 15)

theme_set(theme_gray(base_size = 7))

sv_dependence(shap_lgb_test, colnames(X_test), jitter_width = 0.05, alpha = 0.2) &
  ylim(-0.8, 1.2)

options(repr.plot.width=12, repr.plot.height = 7)

head(test,4)

# Waterfall plot of first observation (log-scale)
sv_waterfall(shap_lgb_test, row_id = 1)

# Waterfall plot of sixth observation
sv_waterfall(shap_lgb_test, row_id = 6)

# Waterfall plot of seventh observation
sv_waterfall(shap_lgb_test, row_id = 7)

# Waterfall plot of eighth observation
sv_waterfall(shap_lgb_test, row_id = 8)


#Compare predictions between models
head(learn, n = 40)




# Basic GLM for severity Age continuous
f.glmsev <- Severity ~ Age + Gender + GeographicZoneGLM + EngineClassGLM + VehicleAge + BonusClassGLM + PolicyDurationGLM
d.glmsev <- glm(f.glmsev, data=learnsev, family=Gamma())
summary(d.glmsev) # print fit statistics and parameter list


# Basic GLM for severity Age categorical
f.glmsev2 <- Severity ~ AgeGLM + Gender + GeographicZoneGLM + EngineClassGLM + VehicleAge + BonusClassGLM + PolicyDurationGLM
d.glmsev2 <- glm(f.glmsev2, data=learnsev, family=Gamma())
summary(d.glmsev2) # print fit statistics and parameter list

# make predictions
learnsev$fitGLMSev <- fitted(d.glmsev)
testsev$fitGLMSev <- predict(d.glmsev, newdata=testsev, type="response")
dat1sev$fitGLMSev <- predict(d.glmsev, newdata=dat1sev, type="response")
dat1sev$fitGLMSev2 <- predict(d.glmsev2, newdata=dat1sev, type="response")

#Output severity dataset with predictions
write.csv(dat1sev, "/content/Dat1sev.csv", row.names = FALSE)


