1 Introduction

This is the third notebook of a series of three that outlines and elaborates upon code used to replicate the central scenario in the paper of Maximilien Baudry “NON-PARAMETRIC INDIVIDUAL CLAIM RESERVING IN INSURANCE”

https://www.institutdesactuaires.com/global/gene/link.php?doc_id=11747&fg=1
http://chaire-dami.fr/files/2016/10/Reserving-article.pdf

In this notebook we step through the process to apply machine learning techniques in order to create reserve estimates following the techniques set out in sections 3 and 4 of Baudry’s paper.

2 Preparations before modelling

The reserving data structures built in this notebook are from a simulated phone insurance dataset. The creation of that simulated data and machine learning data structures for reserving have been set out in the first and second notebooks of this series.

Rather than repeat that code in this notebook the key routines and functions have been put into the r script Baudry_functions_v2.R enabling us to quickly load them for use in this notebook.

# Importing packages

library(knitr)
library(rmdformats)
library(data.table)
library(magrittr)
library(lubridate)
library(ggplot2)
library(cowplot)
library(repr)
library(kableExtra)
library(formattable)

library(xgboost)
library(tidymodels)
library(SHAPforxgboost)

source("./Baudry_functions_v2.R")

3 Simulate data and prepare database

We can now quickly repeat the process from Notebooks 1 and 2 using the functions in the script.

3.1 Simulate policy and claim data

We start with a simulated phone insurance policy and claim dataset. I am calling the function from Notebook 1 of this series to create the dataset. Using a fixed seed will ensure you get a reproducible simulated dataset.

dt_PhoneData <- simulate_central_scenario(1234)

3.2 Join policy and claim data

We then join the simulated phone insurance policy and claim dataset. I have taken the code from Notebook 2 and wrapped it in to a function to perform the join and subsequent data tidying.

Note that if you are working with your own company data you would likely need to amend this function to work with the specifics of your data.

dt_polclaim <- join_policy_claim(dt_PhoneData,
                                 date_pol_start = "date_UW",
                                 date_pol_end = "date_lapse",
                                 date_occur = "date_occur")

3.3 Timeslice policy claim

We then time slice the joined policy and claim dataset. Again, I have taken the code from Notebook 2 and wrapped it into a function to perform the time slicing, given a set of dates that define the time slice periods.

Again, note that this function would need to be re-written to work with your own company data and more importantly it currently would not work correctly with partial payment claims.

lst_Date_slice <- floor_date(seq(as.Date("2016/1/1"), as.Date("2019/06/30"), by = 30), unit= "second") 
dt_polclaim <- time_slice_polclaim(dt_polclaim, lst_Date_slice)

3.4 Select valuation date

Before we create the RBNS and IBNER datasets we must pick a valuation point from the available timeslices. I’m selecting the 10th point for illustration, which is the 10th 30 day period from 01/01/2016, ie a valuation date of 27/09/2016.

i <- valuation <- 10
t_i <- lst_Date_slice[i] 
delta <- min(i, length(lst_Date_slice) - i + 1)

3.5 Creating RBNS dataset

We start with the RBNS dataset and by defining the list of features that will be used in the RBNS model dataset.

As you follow this example you may observe that we are not using the historic claim payment information as an explanatory feature in the RBNS reserve prediction model. Baudry’s paper and approach does allow for the inclusion of such information so it may seem somewhat strange not to include it in this worked example.

I believe the reason that it is not included is a result of the way data has been simulated, specifically that simulated claims are settled with a single payment. That assumption combined with the fact that this method is being applied to individual claim level data prevents the previous time period’s payment information being used by a machine learning model. That is because the previous period’s payment amount can only ever take the value zero for a claim requiring an RBNS reserve. (Claims where the previous paid is non-zero have been paid and settled ie no longer require an RBNS reserve). Consequently, in this central simulated example, prior payment amount is uninformative in the machine learning process and therefore cannot be included as a feature in the RBNS dataset.

The final point to note is that the code shared here deals only with the simplified payment process from Baudry’s central scenario, ie claims are settled with a single payment. Real world data with partial payments and prior payments as explanatory features would require material changes to be made to the code we have shared. Such changes are left for interested readers to make.

#define modelVars
RBNS_model_vars <- c("clm_number",
                     "pol_number",
                     "j",
                     "k",
                     "date_pol_start",
                     "date_occur",
                     "date_report",
                     "date_pay",
                     "Cover",
                     "claim_type",
                     "Brand",
                     "Model",
                     "Price",
                     "target"
    )


# Create a combined TRAIN dataset for k = 1 and all valid j delay values
dt_RBNS_train <- RBNS_Train(dt_polclaim, t_i, i, k = 1, lst_Date_slice, RBNS_model_vars)

# Create a combined TEST dataset for k = 1 and all valid j delay values
dt_RBNS_test <- RBNS_Test(dt_polclaim, t_i, delta, k = 1, lst_Date_slice, RBNS_model_vars)

The train and test datasets are then joined into a single dataset and a small amount of tidying is done to make them ready for use.

# Add a flag to determine which rows are from the trainset and which from the test set
dt_RBNS_train[, flgTrain := 1]
dt_RBNS_test[, flgTrain := 0]

# combine into a single RBNS dataset   
dt_All_RBNS <- rbind(dt_RBNS_train, dt_RBNS_test)
#write.csv(dt_All_RBNS,"dt_All_RBNS.csv", row.names = F)

The important aspects of the tidying relate to creating useable delay metrics from the numerous dates and converting some character features such as cover and claim type into factors.

# order and create some delay fields
setkey(dt_All_RBNS, clm_number, k, j)
    
dt_All_RBNS[, Count := .N , by =clm_number]

#create delay measures and convert time measures from intervals of seconds to intervals of days
dt_All_RBNS[, ':='(
  delay_uw_occ = ifelse(year(date_occur) == 2199,
                        -1,
                        ceiling((as.numeric(date_occur) - as.numeric(date_pol_start)) / (24 * 60 * 60))
                        ),
  delay_occ_rep = ifelse(year(date_occur) == 2199,
                         -1,
                         ceiling((as.numeric(date_report) - as.numeric(date_occur)) / (24 * 60 * 60))
                         ),
  delay_uw_val = ceiling((as.numeric(t_i) - as.numeric(date_pol_start)) / (24 * 60 * 60)),
  delay_rep_pay = ceiling((as.numeric(date_pay) - as.numeric(date_report)) / (24 * 60 * 60)),
  
  date_uw = ceiling(as.numeric(date_pol_start) / (24 *  60 * 60)),
  Cover = as.factor(Cover),
  claim_type = as.factor(claim_type)
  )]

3.6 Creating IBNR dataset

We can then create the IBNR dataset following a similar procedure starting by defining the list of features that will be used in the IBNR model. Then we call the functions that convert the timesliced policy and claim data into the IBNR train and test datasets.

#define IBNR modelVars
IBNR_model_vars <- c("clm_number",
                     "pol_number",
                     "j",
                     "k",
                     "exposure",
                     "date_pol_start",
                     "date_occur",
                     "date_report",
                     "date_pay",
                     "Cover",
                     "Brand",
                     "Model",
                     "Price",
                     "target")
    
# Create a combined TRAIN dataset for k = 1 and all valid j delay values
lst_IBNR_train <- IBNR_Train(dt_polclaim, t_i, i, k = 1,lst_Date_slice, IBNR_model_vars)

# Create a combined TEST dataset for k = 1 and all valid j delay values
dt_IBNR_test <- IBNR_Test(dt_polclaim, t_i, delta, k = 1,lst_Date_slice, IBNR_model_vars)

The train and test datasets are then joined into a single dataset and a small amount of tidying is done to make them ready for use.

lst_IBNR_train$Freq[, flgTrain := 1]
lst_IBNR_train$Loss[, flgTrain := 2]
dt_IBNR_test[, flgTrain := 0]

dt_All_IBNR <- rbind(lst_IBNR_train$Freq,
                     lst_IBNR_train$Loss,
                     dt_IBNR_test)
#write.csv(dt_All_IBNR,"dt_All_IBNR.csv", row.names = F)

The important aspects of the tidying relate to creating useable delay metrics from the numerous dates and converting some character features such as cover and claim type into factors.

# order and create some delay fields
setkey(dt_All_IBNR, clm_number, k, j)
    
dt_All_IBNR[, Count := .N , by =clm_number]
dt_All_IBNR[,':='( delay_uw_occ = ifelse(year(date_occur) == 2199,
                                        -1,
                                        ceiling((as.numeric(date_occur) - as.numeric(date_pol_start))
                                                  /(24*60*60))
                                          ),
                   delay_occ_rep = ifelse(year(date_occur) == 2199,
                                          -1,
                                          ceiling((as.numeric(date_report) - as.numeric(date_occur))
                                                  /(24*60*60))
                                          ),
                   delay_rep_pay = ifelse(year(date_occur) == 2199,
                                          -1,
                                          ceiling((as.numeric(date_pay) - as.numeric(date_report))
                                                  /(24*60*60))
                                          ),
                   delay_uw_val = ceiling((as.numeric(t_i) - as.numeric(date_pol_start))/(24*60*60)),
                   date_uw = ceiling(as.numeric(date_pol_start)/(24*60*60)),
                   Cover = as.factor(Cover))]

That finishes the data preparation. We can now move onto the process of building models of RBNS and IBNR and using them to predict the current reserve requirements.

4 Model build

I will use the xgboost R machine learning package to build the models. In Baudry’s original paper he used the extraTrees package; a variant upon Random Forests.

Xgboost has been selected as it is the most popular implementation of the gradient boosting machine (GBM) algorithm. GBMs have been known since around 2015 to be the most accurate algorithmic technique for regression and classification tasks.

The code and modeling process for each reserve type is very similar so they are shown in the tabbed sections below rather that repeating in-line.

4.1 RBNS model

Starting with the RBNS reserves we will first build a model upon the training data using a cross validated approach to enable selection of key xgboost hyperparameters.

Then using optimal parameters we retrain on the complete training dataset. This trained model is then scored against the test data in order to create predictions from which the RBNS reserve can be calculated.

4.1.1 Creating xgboost dataset

Xgboost requires it’s data to be in a specific format; a special matrix form called a DMatrix, for which there is a function xgb.DMatrix. Not all variable types can be passed to xgb.DMatrix in particular categorical or nominal variables such as Brand have to be converted from text to numeric values.

For this example I’ve chosen to use the parsnip package and create a recipe that converts nominal predictor values into a one-hot-encoded form.

RBNS_predictors <- c("j",
                     "k",
                     "Cover",
                     "claim_type",
                     "Brand",
                     "Model",
                     "Price",
                     #"date_uw",
                     #"delay_uw_occ",
                     "delay_occ_rep")

rowList_RBNS <- list(train=dt_All_RBNS[, which(flgTrain==1)],
                test=dt_All_RBNS[, which(flgTrain==0)],
                all = 1:nrow(dt_All_RBNS))


RBNS_rec <- recipe( ~ ., data = dt_All_RBNS[, RBNS_predictors, with = FALSE]) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  prep()


df.RBNS_train <- bake(RBNS_rec, new_data = dt_All_RBNS[rowList_RBNS$train,] )
df.RBNS_test <- bake(RBNS_rec, new_data = dt_All_RBNS[rowList_RBNS$test,] )
df.RBNS_all <- bake(RBNS_rec, new_data = dt_All_RBNS )


xgb.RBNS_DMat.train <- xgb.DMatrix(data = as.matrix(df.RBNS_train),
                              label = dt_All_RBNS[rowList_RBNS$train, target])

xgb.RBNS_DMat.test <- xgb.DMatrix(data = as.matrix(df.RBNS_test),
                              label = dt_All_RBNS[rowList_RBNS$test, target])

xgb.RBNS_DMat.all <- xgb.DMatrix(data = as.matrix(df.RBNS_all),
                             label = dt_All_RBNS[, target])

4.1.2 Fit initial model using cross validation

Having prepared the data for xgboost I can now fit an initial model. I’ve used a simple set of hyper-parameters and used cross validation to select the optimal number of boosted trees (nrounds) for these hyper-parameter selections by calling the xgb.cv function with early stopping.

I have used the reg:tweedie objective function based upon inspection of the target variable.

summary(dt_All_RBNS[rowList_RBNS$train, target])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00   76.44    0.00  913.00

hist(dt_All_RBNS[rowList_RBNS$train, target])

In a real world example a process of hyper-parameter tuning would be undertaken in order to select more optimal xgboost hyper-parameters.

param <- list(
  objective = "reg:tweedie",
  max_depth = 2L,            # tree-depth
  subsample = 0.75,          # randomly sample rows before fitting each tree
  colsample_bytree = 0.8,    # randomly sample columns before fitting each tree
  min.child.weight = 10,     # minimum weight per leaf
  eta = 0.1                  # Learning rate
)


# Train model with cross validation
set.seed(1984) # for repeatability

xgb_RBNS_CV <- xgb.cv(
  params                 = param,
  data                   = xgb.RBNS_DMat.train,
  nrounds                = 500,        # Maximum number of trees to build
  nfold = 5,
  early_stopping_rounds  = 10L,        # Stops algorithm early if performance has not improved in 
  print_every_n          = 10L,        # How often to print to console
  prediction             = TRUE        # Keeps the predictions
)

## [1]  train-tweedie-nloglik@1.5:197.331256+2.882719   test-tweedie-nloglik@1.5:197.332748+11.169201 
## Multiple eval metrics are present. Will use test_tweedie_nloglik@1.5 for early stopping.
## Will train until test_tweedie_nloglik@1.5 hasn't improved in 10 rounds.
## 
## [11] train-tweedie-nloglik@1.5:76.005264+0.635454    test-tweedie-nloglik@1.5:76.023500+4.728197 
## [21] train-tweedie-nloglik@1.5:34.399348+0.339311    test-tweedie-nloglik@1.5:34.444048+2.235487 
## [31] train-tweedie-nloglik@1.5:22.642844+0.121015    test-tweedie-nloglik@1.5:22.659766+1.318212 
## [41] train-tweedie-nloglik@1.5:20.103612+0.194808    test-tweedie-nloglik@1.5:20.130332+1.028970 
## [51] train-tweedie-nloglik@1.5:19.614061+0.213816    test-tweedie-nloglik@1.5:19.648131+0.982517 
## [61] train-tweedie-nloglik@1.5:19.488429+0.224332    test-tweedie-nloglik@1.5:19.531806+0.963364 
## [71] train-tweedie-nloglik@1.5:19.438496+0.229845    test-tweedie-nloglik@1.5:19.484217+0.961241 
## [81] train-tweedie-nloglik@1.5:19.411388+0.236892    test-tweedie-nloglik@1.5:19.464310+0.958345 
## [91] train-tweedie-nloglik@1.5:19.392717+0.238659    test-tweedie-nloglik@1.5:19.457403+0.962381 
## [101]    train-tweedie-nloglik@1.5:19.378791+0.238025    test-tweedie-nloglik@1.5:19.453831+0.959511 
## [111]    train-tweedie-nloglik@1.5:19.367772+0.238858    test-tweedie-nloglik@1.5:19.450139+0.962756 
## [121]    train-tweedie-nloglik@1.5:19.356093+0.239343    test-tweedie-nloglik@1.5:19.458886+0.961838 
## Stopping. Best iteration:
## [113]    train-tweedie-nloglik@1.5:19.364408+0.238899    test-tweedie-nloglik@1.5:19.449566+0.961121

Having fit the model we store the out-of-fold predictions.

dt_All_RBNS[rowList_RBNS$train, preds_oof := xgb_RBNS_CV$pred]

4.1.3 Fit final model on all training data

Having fit the model using 5 fold cross validation we observe the optimum number of fitting rounds to be 113.

We can then use this to train a final model on all the data.

xgb_RBNS_Fit <- xgb.train(
  params                 = param,
  data                   = xgb.RBNS_DMat.train,
  nrounds                = xgb_RBNS_CV$best_iteration,
# base_score             = 1,
  watchlist              = list(train=xgb.RBNS_DMat.train, test=xgb.RBNS_DMat.test) ,
  print_every_n          = 10
)

## [1]  train-tweedie-nloglik@1.5:197.150055    test-tweedie-nloglik@1.5:82.789551 
## [11] train-tweedie-nloglik@1.5:75.865128 test-tweedie-nloglik@1.5:32.100452 
## [21] train-tweedie-nloglik@1.5:34.546490 test-tweedie-nloglik@1.5:14.854555 
## [31] train-tweedie-nloglik@1.5:22.767673 test-tweedie-nloglik@1.5:9.884780 
## [41] train-tweedie-nloglik@1.5:20.138584 test-tweedie-nloglik@1.5:8.672418 
## [51] train-tweedie-nloglik@1.5:19.639982 test-tweedie-nloglik@1.5:8.421473 
## [61] train-tweedie-nloglik@1.5:19.497528 test-tweedie-nloglik@1.5:8.333148 
## [71] train-tweedie-nloglik@1.5:19.442001 test-tweedie-nloglik@1.5:8.295179 
## [81] train-tweedie-nloglik@1.5:19.415724 test-tweedie-nloglik@1.5:8.274471 
## [91] train-tweedie-nloglik@1.5:19.400574 test-tweedie-nloglik@1.5:8.264827 
## [101]    train-tweedie-nloglik@1.5:19.389402 test-tweedie-nloglik@1.5:8.261550 
## [111]    train-tweedie-nloglik@1.5:19.382881 test-tweedie-nloglik@1.5:8.261103 
## [113]    train-tweedie-nloglik@1.5:19.381304 test-tweedie-nloglik@1.5:8.260632

dt_All_RBNS[, preds_full := predict(xgb_RBNS_Fit,xgb.RBNS_DMat.all)]

4.1.4 Inspect model fit

Having fitted the full model we can then inspect the model fit. The traditional way of inspecting global model feature importance is to use the gains chart.

#default feature importance by gain
featImp_RBNS <- xgb.importance(xgb_RBNS_Fit, feature_names = colnames(xgb.RBNS_DMat.train))
xgb.plot.importance(featImp_RBNS, main="Feature Importance - RBNS")

An increasingly popular and more robust approach is to use SHAP values https://github.com/slundberg/shap. The SHAP equivalent of the feature importance chart is shown below.

# Return the SHAP values and ranked features by mean|SHAP|
shap_values <- shap.values(xgb_model = xgb_RBNS_Fit, X_train = as.matrix(df.RBNS_train))

# Prepare the long-format data:
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train =  as.matrix(df.RBNS_train))

# **SHAP summary plot**
shap.plot.summary(shap_long, dilute = nrow(df.RBNS_train)/10000)

A second useful chart is the partial dependence plot. This shows how the values of a predictive (input) feature influence the predicted (output) value, while holding the values of all other predictive (input) features constant. It is also known as the marginal effect.

Here we show the partial dependence plots for the top 4 SHAP features.

fig_list <- lapply(names(shap_values$mean_shap_score)[1:4], 
                   shap.plot.dependence,
                   data_long = shap_long,
                   dilute = nrow(shap_long)/ 10000)

gridExtra::grid.arrange(grobs = fig_list, ncol = 2)

The feature importance and partial dependency plots provide quick insight into the model.

We see that claim development period, j, is the most important feature and that the RBNS reserve is smaller for larger values of j.
The second most important feature is claimtype where breakage claims are associated with lower RBNS reserves. This is as expected from the data generating process where breakage claims come from a Beta distribution with a lower mean.
The third most important feature is phone price where there is linear increase in RBNS reserve with increasing phone price. This is also as expected from the data generating process.
The fourth feature is phone brand which again follows expectations from the data generating process.

SHAP values can also be used to show the components of a single prediction. In the following plot we show the top 4 components for each row of the data and zoom in at row 500.

# choose to show top 4 features by setting `top_n = 4`, set 6 clustering groups.  
plot_data <- shap.prep.stack.data(shap_contrib = shap_values$shap_score,
                                  data_percent = 10000/nrow(shap_long),
                                  top_n = 4,
                                  n_groups = 6)
  
# choose to zoom in at location 500, set y-axis limit using `y_parent_limit`  
# it is also possible to set y-axis limit for zoom-in part alone using `y_zoomin_limit`  
shap.plot.force_plot(plot_data, zoom_in_location = 500, y_parent_limit = c(-1,1))

4.1.5 Summarising RBNS reserves

By comparing the model predictions to the simulated claims run-off we can get a feel for the accuracy of the machine learning approach. Below I aggregate the RBNS predictions by claim occurrence month and compare them to the known simulated claim run-off.

dt_All_RBNS [, date_occur_YYYYMM := as.character(year(date_occur) + month(date_occur)/100 )]

dt_RBNS_summary <- dt_All_RBNS[rowList_RBNS$test,.(preds = sum(preds_full), target = sum(target)), keyby = date_occur_YYYYMM]

date_occur_YYYYMM	preds	target	Diff	Diff_pcnt
2016.04	1,444	911	533	58.5%
2016.05	5,056	5,188	-132	-2.5%
2016.06	12,569	14,257	-1,688	-11.8%
2016.07	56,427	53,815	2,612	4.9%
2016.08	325,836	326,631	-795	-0.2%
2016.09	622,619	622,179	440	0.1%
Total	1,023,951	1,022,981	970	0.1%

You can now jump back up to the beginning of the modeling section and select the IBNR frequency modeling tab.

4.2 IBNR Frequency model

IBNR reserves are built by multiplying the outpts of a claim frequency and claim severity model. The general process of building the model follows that of the RBNS reserves.

4.2.1 Creating xgboost dataset

Again we convert the data into a specific matrix form called a DMatrix. However before doing so we aggregate the data for efficiency of model fit times. There could be a loss of accuracy in aggregating data, so in practice this is something you would want to experiment with.

The other point to note is the values of flgTrain used to identify the training and test rows in our dataset. Recall from Notebook 2, in the IBNR dataset training rows for the frequency model have a flgtrain value of 1 whereas the Severity training rows have a value of 2.

IBNR_predictors <- c("j",
                     "k",
                     "Cover",
                     "Brand",
                     "Model",
                     "Price",
                     "date_uw")

# aggregate the data ... does this lead to loss of variance and accuracy?
dt_All_IBNR [, date_pol_start_YYYYMM := as.character(year(date_pol_start) + month(date_pol_start)/100 )]

dt_All_IBNR_F <- dt_All_IBNR[, .(exposure = sum(exposure),
                                   target_cost = sum(target),
                                   target_count = sum(target>0)),
                               by= c(IBNR_predictors, "date_pol_start_YYYYMM", "flgTrain")]

dt_All_IBNR_F <- dt_All_IBNR_F[exposure>0]


# setup train and test rows
rowList_IBNR_F <- list(train=dt_All_IBNR_F[, which(flgTrain==1)],
                     test=dt_All_IBNR_F[, which(flgTrain==0)],
                     all = dt_All_IBNR_F[, which(flgTrain!=2)])

# setup data for xgboost
IBNR_rec <- recipe( ~ ., data = dt_All_IBNR_F[, IBNR_predictors, with = FALSE]) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  prep()

df.IBNR_F_train <- bake(IBNR_rec, new_data = dt_All_IBNR_F[rowList_IBNR_F$train,] )
df.IBNR_F_test <- bake(IBNR_rec, new_data = dt_All_IBNR_F[rowList_IBNR_F$test,] )
df.IBNR_F_all <- bake(IBNR_rec, new_data = dt_All_IBNR_F[rowList_IBNR_F$all,] )

xgb.IBNR_F_DMat.train <- xgb.DMatrix(data = as.matrix(df.IBNR_F_train),
                              weight = dt_All_IBNR_F[rowList_IBNR_F$train, exposure],
                              label = dt_All_IBNR_F[rowList_IBNR_F$train, target_count])

xgb.IBNR_F_DMat.test <- xgb.DMatrix(data = as.matrix(df.IBNR_F_test),
                             weight = dt_All_IBNR_F[rowList_IBNR_F$test, exposure],
                             label = dt_All_IBNR_F[rowList_IBNR_F$test, target_count])

xgb.IBNR_F_DMat.all <- xgb.DMatrix(data = as.matrix(df.IBNR_F_all),
                            weight = dt_All_IBNR_F[rowList_IBNR_F$all, exposure],
                            label = dt_All_IBNR_F[rowList_IBNR_F$all, target_count])

4.2.2 Fit initial model using cross validation

I have used the count:poisson objective function based upon inspection of the target variable.

hist(dt_All_IBNR_F[rowList_IBNR_F$train, target_count])

param <- list(
  objective = "count:poisson",
  max_depth = 2L,           # tree-depth
  subsample = 0.7,          # randomly sample rows before fitting each tree
  colsample_bytree = 0.8,   # randomly sample columns before fitting each tree
  min.child.weight = 10,    # minimum weight per leaf
  eta = 0.1               # Learning rate
  #monotone_constraints = monotone_Vec # Monotonicity constraints
)


# Train model with cross validation
set.seed(1984) # for repeatability

xgb_IBNR_F_CV <- xgb.cv(
 params                 = param,
 data                   = xgb.IBNR_F_DMat.train,
 nrounds                = 2000,        # Maximum number of trees to build
 nfold = 5,

 early_stopping_rounds  = 50L,        # Stops algorithm early if performance has not improved in n rounds
 print_every_n          = 50L,        # How often to print to console
 #base_score             = 0.001,       # Model starting point
   prediction             = TRUE        # Keeps the predictions
)

## [1]  train-poisson-nloglik:0.509656+0.000600 test-poisson-nloglik:0.509657+0.002292 
## Multiple eval metrics are present. Will use test_poisson_nloglik for early stopping.
## Will train until test_poisson_nloglik hasn't improved in 50 rounds.
## 
## [51] train-poisson-nloglik:0.162205+0.001478 test-poisson-nloglik:0.162599+0.004882 
## [101]    train-poisson-nloglik:0.136502+0.001222 test-poisson-nloglik:0.137285+0.005015 
## [151]    train-poisson-nloglik:0.132466+0.001164 test-poisson-nloglik:0.133936+0.004689 
## [201]    train-poisson-nloglik:0.130754+0.001130 test-poisson-nloglik:0.132831+0.004558 
## [251]    train-poisson-nloglik:0.129548+0.001087 test-poisson-nloglik:0.132461+0.004548 
## [301]    train-poisson-nloglik:0.128583+0.001078 test-poisson-nloglik:0.132294+0.004463 
## [351]    train-poisson-nloglik:0.127722+0.001084 test-poisson-nloglik:0.132064+0.004437 
## [401]    train-poisson-nloglik:0.126975+0.001102 test-poisson-nloglik:0.131920+0.004426 
## [451]    train-poisson-nloglik:0.126253+0.001077 test-poisson-nloglik:0.131680+0.004414 
## [501]    train-poisson-nloglik:0.125648+0.001070 test-poisson-nloglik:0.131608+0.004397 
## [551]    train-poisson-nloglik:0.125077+0.001016 test-poisson-nloglik:0.131505+0.004303 
## [601]    train-poisson-nloglik:0.124579+0.000978 test-poisson-nloglik:0.131449+0.004381 
## Stopping. Best iteration:
## [579]    train-poisson-nloglik:0.124774+0.001003 test-poisson-nloglik:0.131380+0.004327

Having fit the model we store the out-of-fold predictions for both the claim frequency and the claim counts.

 dt_All_IBNR_F[rowList_IBNR_F$train, preds_oof_IBNR_F := xgb_IBNR_F_CV$pred]
 dt_All_IBNR_F[rowList_IBNR_F$train, preds_oof_IBNR_Nos := exposure * preds_oof_IBNR_F]

4.2.3 Fit final model on all training data

Having fit the model using 5 fold cross validation we observe the optimum number of fitting rounds to be 579.

We can then use this to train a final model on all the data.

xgb_IBNR_F_Fit <- xgb.train(
   params                 = param,
   data                   = xgb.IBNR_F_DMat.train,
   nrounds                = xgb_IBNR_F_CV$best_iteration,
# base_score             = 1,
   watchlist              = list(train=xgb.IBNR_F_DMat.train, test=xgb.IBNR_F_DMat.test) ,
   print_every_n          = 50
 )

## [1]  train-poisson-nloglik:0.509601  test-poisson-nloglik:0.497376 
## [51] train-poisson-nloglik:0.163437  test-poisson-nloglik:0.122574 
## [101]    train-poisson-nloglik:0.136768  test-poisson-nloglik:0.088578 
## [151]    train-poisson-nloglik:0.132694  test-poisson-nloglik:0.084188 
## [201]    train-poisson-nloglik:0.131227  test-poisson-nloglik:0.083316 
## [251]    train-poisson-nloglik:0.130050  test-poisson-nloglik:0.082979 
## [301]    train-poisson-nloglik:0.129073  test-poisson-nloglik:0.083072 
## [351]    train-poisson-nloglik:0.128206  test-poisson-nloglik:0.083234 
## [401]    train-poisson-nloglik:0.127501  test-poisson-nloglik:0.083556 
## [451]    train-poisson-nloglik:0.126870  test-poisson-nloglik:0.083565 
## [501]    train-poisson-nloglik:0.126326  test-poisson-nloglik:0.083704 
## [551]    train-poisson-nloglik:0.125827  test-poisson-nloglik:0.083802 
## [579]    train-poisson-nloglik:0.125509  test-poisson-nloglik:0.083817

 dt_All_IBNR_F[rowList_IBNR_F$all, preds_full_IBNR_Nos := predict(xgb_IBNR_F_Fit,xgb.IBNR_F_DMat.all)]

4.2.4 Inspect model fit

Having fitted the full model we can then inspect the model fit. The traditional way of inspecting global model feature importance is to use the gains chart.

#default feature importance by gain
featImp_IBNR_F <- xgb.importance(xgb_IBNR_F_Fit, feature_names = colnames(xgb.IBNR_F_DMat.train))
xgb.plot.importance(featImp_IBNR_F, main="Feature Importance - IBNR Frequency")

An increasingly popular and more robust approach is to use SHAP values https://github.com/slundberg/shap. The SHAP equivalent of the feature importance chart is shown below.

# Return the SHAP values and ranked features by mean|SHAP|
shap_values <- shap.values(xgb_model = xgb_IBNR_F_Fit, X_train = as.matrix(df.IBNR_F_train))

# Prepare the long-format data:
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train =  as.matrix(df.IBNR_F_train))

# **SHAP summary plot**
shap.plot.summary(shap_long, dilute = nrow(df.IBNR_F_train)/10000)

Here we show the partial dependence plots for the top 4 SHAP features.

fig_list <- lapply(names(shap_values$mean_shap_score)[1:4], 
                   shap.plot.dependence,
                   data_long = shap_long,
                   dilute = nrow(shap_long)/ 10000)

gridExtra::grid.arrange(grobs = fig_list, ncol = 2)

The feature importance and partial dependency plots provide quick insight into the model.

We see that claim development period, j, is the most important feature and that the IBNR frequency is smaller for larger values of j.
The second most important feature is phone price with the IBNR claim count increasing as price increases. Although this is not a direct feature in the data generating process there is a link with higher theft frequencies and higher phone prices.
The third most important feature is phone model with IBNR frequency increasing with model type. This is expected from the data generating process as theft claim frequency increases with model type.
The fourth feature is phone cover and it should be no suprise to see that Breakage only, the most basic cover level is associated with lower IBNR claim counts.

SHAP values can also be used to show the components of a single predction. In the following plot we show the top 4 components for each row of the data and zoom in at row 500.

# choose to show top 4 features by setting `top_n = 4`, set 6 clustering groups.  
plot_data <- shap.prep.stack.data(shap_contrib = shap_values$shap_score,
                                  data_percent = 10000/nrow(shap_long),
                                  top_n = 4,
                                  n_groups = 6)
  
# choose to zoom in at location 500, set y-axis limit using `y_parent_limit`  
# it is also possible to set y-axis limit for zoom-in part alone using `y_zoomin_limit`  
shap.plot.force_plot(plot_data, zoom_in_location = 500, y_parent_limit = c(-1,1))

4.2.5 Summarising IBNR claim counts

By comparing the model predictions to the simulated claims run-off we can get a feel for the accuracy of the machine learning approach. Below I aggregate the IBNR claim count predictions by claim occurrence month and compare them to the known simulated claim run-off.

dt_All_IBNR_F_summary <- dt_All_IBNR_F[rowList_IBNR_F$test,.(preds = sum(preds_full_IBNR_Nos), target = sum(target_count)), keyby = date_pol_start_YYYYMM]

date_pol_start_YYYYMM	preds	target	Diff	Diff_pcnt
2016.01	206	168	38	22.6%
2016.02	178	176	2	1.4%
2016.03	184	160	24	14.8%
2016.04	174	164	10	6.2%
2016.05	155	195	-40	-20.7%
2016.06	149	179	-30	-16.7%
2016.07	89	170	-81	-47.9%
2016.08	66	162	-96	-59.2%
2016.09	66	59	7	12.3%
Total	1,267	1,433	-166	-11.6%

You can now jump back up to the beginning of the modeling section and select the IBNR severity modeling tab.

4.3 IBNR Severity model

IBNR reserves are built by multiplying the outputs of a claim frequency and claim severity model. The general process of building the model follows that of the RBNS reserves.

4.3.1 Creating xgboost dataset

Again we convert the data into a specific matrix form called a DMatrix.

The other point to note is the values of flgTrain used to identify the training and test rows in our dataset. Recall from Notebook 2, in the IBNR dataset training rows for the frequency model have a flgtrain value of 1 whereas the Severity training roes have a value of 2.

IBNR_predictors <- c("j",
                     "k",
                     "Cover",
                     "Brand",
                     "Model",
                     "Price",
                     "date_uw")

dt_All_IBNR [, date_pol_start_YYYYMM := as.character(year(date_pol_start) + month(date_pol_start)/100 )]

dt_All_IBNR_S <- dt_All_IBNR[, .(exposure = sum(target>0),
                                 target_cost = sum(target)),
                            by= c(IBNR_predictors, "date_pol_start_YYYYMM", "flgTrain")]

dt_All_IBNR_S <- dt_All_IBNR_S[exposure>0]

# setup train and test rows
rowList_IBNR_S <- list(train=dt_All_IBNR_S[, which(flgTrain==2)],
                     test=dt_All_IBNR_S[, which(flgTrain==0)],
                     all = dt_All_IBNR_S[, which(flgTrain!=1)])

# setup data for xgboost

IBNR_rec <- recipe( ~ ., data = dt_All_IBNR_S[, IBNR_predictors, with = FALSE]) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  prep()

df.IBNR_S_train <- bake(IBNR_rec, new_data = dt_All_IBNR_S[rowList_IBNR_S$train,] )
df.IBNR_S_test <- bake(IBNR_rec, new_data = dt_All_IBNR_S[rowList_IBNR_S$test,] )
df.IBNR_S_all <- bake(IBNR_rec, new_data = dt_All_IBNR_S[rowList_IBNR_S$all,] )



xgb.IBNR_S_DMat.train <- xgb.DMatrix(data = as.matrix(df.IBNR_S_train),
                              weight = dt_All_IBNR_S[rowList_IBNR_S$train, exposure],
                              label = dt_All_IBNR_S[rowList_IBNR_S$train, target_cost])

xgb.IBNR_S_DMat.test <- xgb.DMatrix(data = as.matrix(df.IBNR_S_test),
                             weight = dt_All_IBNR_S[rowList_IBNR_S$test, exposure],
                             label = dt_All_IBNR_S[rowList_IBNR_S$test, target_cost])

xgb.IBNR_S_DMat.all <- xgb.DMatrix(data = as.matrix(df.IBNR_S_all),
                            weight = dt_All_IBNR_S[rowList_IBNR_S$all, exposure],
                            label = dt_All_IBNR_S[rowList_IBNR_S$all, target_cost])

4.3.2 Fit initial model using cross validation

I have used the reg:gamma objective function based upon inspection of the target variable.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     3.0   125.5   255.0   353.9   486.5  2011.0

param <- list(
  objective = "reg:gamma",
  max_depth = 2L,           # tree-depth
  subsample = 0.7,          # randomly sample rows before fitting each tree
  colsample_bytree = 0.8,   # randomly sample columns before fitting each tree
  min.child.weight = 10,    # minimum weight per leaf
  eta = 0.1               # Learning rate
  #monotone_constraints = monotone_Vec # Monotonicity constraints
)

# Train model with cross validation
set.seed(1984) # for repeatability

xgb_IBNR_S_CV <- xgb.cv(
 params                 = param,
 data                   = xgb.IBNR_S_DMat.train,
 nrounds                = 2000,        # Maximum number of trees to build
 nfold = 5,

 early_stopping_rounds  = 50L,        # Stops algorithm early if performance has not improved in n rounds
 print_every_n          = 50L,        # How often to print to console
 #base_score             = 0.001,       # Model starting point
   prediction             = TRUE        # Keeps the predictions
)

## [1]  train-gamma-nloglik:725.805481+6.987979 test-gamma-nloglik:725.496777+28.009135 
## Multiple eval metrics are present. Will use test_gamma_nloglik for early stopping.
## Will train until test_gamma_nloglik hasn't improved in 50 rounds.
## 
## [51] train-gamma-nloglik:10.027304+0.047412  test-gamma-nloglik:10.030309+0.237108 
## [101]    train-gamma-nloglik:6.820591+0.008904   test-gamma-nloglik:6.828172+0.038464 
## [151]    train-gamma-nloglik:6.810284+0.008580   test-gamma-nloglik:6.826305+0.036504 
## Stopping. Best iteration:
## [141]    train-gamma-nloglik:6.811574+0.008407   test-gamma-nloglik:6.825954+0.036434

Having fitted an initial model the out-of-fold predictions are stored.

 dt_All_IBNR_S[rowList_IBNR_S$train, preds_oof_IBNR_S := xgb_IBNR_S_CV$pred]
 dt_All_IBNR_S[rowList_IBNR_S$train, preds_oof_IBNR_Cost := exposure * preds_oof_IBNR_S]

4.3.3 Fit final model on all training data

Having fit the model using 5 fold cross validation we observe the optimum number of fitting rounds to be 141.

We can then us this to train a final model on all the data.

xgb_IBNR_S_Fit <- xgb.train(
   params                 = param,
   data                   = xgb.IBNR_S_DMat.train,
   nrounds                = xgb_IBNR_S_CV$best_iteration,
# base_score             = 1,
   watchlist              = list(train=xgb.IBNR_F_DMat.train, test=xgb.IBNR_F_DMat.test) ,
   print_every_n          = 50
 )

## [1]  train-gamma-nloglik:NaN test-gamma-nloglik:NaN 
## [51] train-gamma-nloglik:NaN test-gamma-nloglik:NaN 
## [101]    train-gamma-nloglik:NaN test-gamma-nloglik:NaN 
## [141]    train-gamma-nloglik:NaN test-gamma-nloglik:NaN

Having trained the model the predictions are stored.

 dt_All_IBNR_S[rowList_IBNR_S$all, preds_full_IBNR_Cost := predict(xgb_IBNR_S_Fit,xgb.IBNR_S_DMat.all)]

4.3.4 Inspect model fit

Having fitted the full model we can then inspect the model fit. The traditional way of inspecting global model feature importance is to use the gains chart.

#default feature importance by gain
featImp_IBNR_S <- xgb.importance(xgb_IBNR_S_Fit, feature_names = colnames(xgb.IBNR_S_DMat.train))
xgb.plot.importance(featImp_IBNR_S, main="Feature Importance - IBNR Severity")

An increasingly popular and more robust approach is to use SHAP values https://github.com/slundberg/shap. The SHAP equivalent of the feature importance chart is shown below.

# Return the SHAP values and ranked features by mean|SHAP|
shap_values <- shap.values(xgb_model = xgb_IBNR_S_Fit, X_train = as.matrix(df.IBNR_S_train))

# Prepare the long-format data:
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train =  as.matrix(df.IBNR_S_train))

# **SHAP summary plot**
shap.plot.summary(shap_long, dilute = max(nrow(df.IBNR_S_train),10000)/10000)

Here we show the partial dependence plots for the top 4 SHAP features.

fig_list <- lapply(names(shap_values$mean_shap_score)[1:4], 
                   shap.plot.dependence,
                   data_long = shap_long,
                   dilute = nrow(shap_long)/ 10000)

gridExtra::grid.arrange(grobs = fig_list, ncol = 2)

The feature importance and partial dependency plots provide quick insight into the model.

We see that phone price is the most important feature and has a linear relationship with IBNR severity as we would expect from the data generating process.
The second and third most important features relate to phone cover. We see that cover BOT is associated with higher IBNR costs whereas cover B is associated with lower.
The fourth feature is phone brand and follows the pattern we would expect from the data simulation process.

SHAP values can also be used to show the components of a single prediction. In the following plot we show the top 4 components for each row of the data and zoom in at row 500.

# choose to show top 4 features by setting `top_n = 4`, set 6 clustering groups.  
plot_data <- shap.prep.stack.data(shap_contrib = shap_values$shap_score,
                                  data_percent = 10000/max(nrow(shap_long),10000),
                                  top_n = 4,
                                  n_groups = 6)

## The SHAP values of the Rest 5 features were summed into variable 'rest_variables'.

# choose to zoom in at location 500, set y-axis limit using `y_parent_limit`  
# it is also possible to set y-axis limit for zoom-in part alone using `y_zoomin_limit`  
shap.plot.force_plot(plot_data, zoom_in_location = 500, y_parent_limit = c(-1,1))

## Data has N = 1111 | zoom in length is 111 at location 500.

4.3.5 Summarising IBNR claim costs

By comparing the model predictions to the simulated claims run-off we can get a feel for the accuracy of the machine learning approach. Below I aggregate the IBNR claim cost predictions by claim occurrence month and compare them to the known simulated claim run-off.

dt_All_IBNR_S_summary <- dt_All_IBNR_S[rowList_IBNR_S$test,.(preds = sum(preds_full_IBNR_Cost), target = sum(target_cost)), keyby = date_pol_start_YYYYMM]

date_pol_start_YYYYMM	preds	target	Diff	Diff_pcnt
2016.01	57,279	52,782	4,497	8.5%
2016.02	54,144	55,615	-1,471	-2.6%
2016.03	51,057	47,200	3,857	8.2%
2016.04	53,099	56,975	-3,876	-6.8%
2016.05	63,855	64,833	-978	-1.5%
2016.06	52,763	52,989	-226	-0.4%
2016.07	51,072	50,276	796	1.6%
2016.08	51,240	53,502	-2,262	-4.2%
2016.09	18,392	18,989	-597	-3.1%
Total	452,900	453,161	-261	-0.1%

5 Summary

In this third notebook we have stepped through the process to apply machine learning techniques in order to create reserve estimates following the techniques set out in sections 3 and 4 of Baudry’s paper.

Specifically we have shown how to apply the xgboost machine learning algorithm and illustrated how it can, with little human input, create reasonable reserve estimates as shown below.

Reserve	Prediction	Ground_Truth	Diff	Diff_pcnt
RBNS	1,023,951	1,022,981	970	0.1%
IBNR	452,900	453,161	-261	-0.1%
TOTAL	1,476,851	1,476,142	709	0.0%

In addition we have seen how we can create explanations for the reserve predictions using feature importance and partial dependence plots. We have also used SHAP values which can be used to explain both global features and individual predictions.

In so doing, it is hoped that this series of notebooks encourage further exploration of Baudry’s paper and the wider field of machine learning in reserving. The code supporting this notebook along with the supporting files and folders available in a zip file here.

Download and extract the zip file to a local directory and then open the R project file Baudry_3.rproj in your local R software installation. In the root of the project folder you will see two files;

Notebook_3_ApplyMachineLearningReserving_v1.Rmd - which is the source code used to recreate this notebook
Notebook_3_ApplyMachineLearningReserving_v1.R - the equivalent code provided as an R script

Please note that, depending upon your R installation , you may have to install R libraries before you can run the code provided. R will warn you if you have missing dependencies and you can then install them from CRAN.

NON-PARAMETRIC INDIVIDUAL CLAIM RESERVING IN INSURANCE

Notebook 3 of 3: Applying machine learning for reserving

NON-PARAMETRIC INDIVIDUAL CLAIM RESERVING IN INSURANCE

1 Introduction

2 Preparations before modelling

3 Simulate data and prepare database

3.1 Simulate policy and claim data

3.2 Join policy and claim data

3.3 Timeslice policy claim

3.4 Select valuation date

3.5 Creating RBNS dataset

3.6 Creating IBNR dataset

4 Model build

4.1 RBNS model

4.1.1 Creating xgboost dataset

4.1.2 Fit initial model using cross validation

4.1.3 Fit final model on all training data

4.1.4 Inspect model fit

4.1.5 Summarising RBNS reserves

4.2 IBNR Frequency model

4.2.1 Creating xgboost dataset

4.2.2 Fit initial model using cross validation

4.2.3 Fit final model on all training data

4.2.4 Inspect model fit

4.2.5 Summarising IBNR claim counts

4.3 IBNR Severity model

4.3.1 Creating xgboost dataset

4.3.2 Fit initial model using cross validation

4.3.3 Fit final model on all training data

4.3.4 Inspect model fit

4.3.5 Summarising IBNR claim costs

5 Summary