## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) ## ----pure IBNR example-------------------------------------------------------- library("NetSimR") #set the seed set.seed(0) #Set the reporting delay distribution parameters mu <- 3 sigma <- 1 #Generate reporting delays and preview the distribution x<-rlnorm(1000,mu,sigma) summary(x) hist(x, breaks = 100) ## ----pure IBNR example continued---------------------------------------------- #Generate dates data PoliciesPerDay<-100 PeriodLength<-365*2 InceptionDateStart<-as.Date("2011/1/1") InceptionDate<-seq(InceptionDateStart, by = "day", length.out = PeriodLength) InceptionDate<-rep(InceptionDate,PoliciesPerDay) #assume all policies are yearly DayPolicyDuration<-365 ExpiryDate<-InceptionDate+DayPolicyDuration ValuationDate<-max(InceptionDate) Data<-data.frame(InceptionDate,ExpiryDate) summary(Data) ValuationDate #Generate Claims Counts Data and reporting delays #Assume only 1 claim per policy possible with probability of claiming 10% ClaimFreq<-0.1 Data$ClaimCount<-ifelse(runif(nrow(Data))<ClaimFreq,1,0) #generate accident date assuming uniformly happening throughout the year Data$AccidentDate<-Data$InceptionDate+runif(nrow(Data))*(Data$ExpiryDate-Data$InceptionDate) #generate reporting delays Data$ReportingDelay<-round(Data$ClaimCount*rlnorm(nrow(Data),mu,sigma),1) #generate reporting days Data$NotificationDate<-Data$AccidentDate+Data$ReportingDelay #remove claims reported beyond valuation day Data$ClaimReportedBeforeValuation<-ifelse(Data$NotificationDate>ValuationDate,0,Data$ClaimCount) #Maximum reporting delays Data$MaxReportingDelay<-as.numeric(ValuationDate-Data$AccidentDate)+1 #Earned Duration Data$EarnedDuration<-ifelse(Data$ExpiryDate>ValuationDate,as.numeric(ValuationDate-Data$InceptionDate),as.numeric(Data$ExpiryDate-Data$InceptionDate))/365 #Filter claims only data ModelData<-Data[Data$ClaimReportedBeforeValuation==1,] #model reporting delays with a linear model lm<-lm(log(ModelData$ReportingDelay)~1) #summary(lm) cbind(coefficients(lm),summary(lm)$sigma) #model reporting delays with a truncated model (splitting to monthly periods would be more accurate in terms of maximum reporting delay) library("crch") rTLm<-crch(log(ModelData$ReportingDelay)~1, truncated = TRUE, right = log(ModelData$MaxReportingDelay), dist = "gaussian", link.scale = "identity") #summary(rTLm) coefficients(rTLm) #Generate and predict pure IBNR and UPR Claims Data$PureIBNRClaims<-ifelse(Data$AccidentDate>ValuationDate,0,Data$ClaimCount-Data$ClaimReportedBeforeValuation) Data$UPRClaims<-Data$ClaimCount-Data$ClaimReportedBeforeValuation-Data$PureIBNR PureIBNRPrediction<-PureIBNRLNorm(Data$InceptionDate,Data$ExpiryDate,ValuationDate,coefficients(rTLm)[1],coefficients(rTLm)[2]) Data$PredictedPureIBNRYears<-PureIBNRPrediction$PureIBNRDuration/DayPolicyDuration Data$PredictedUPRYears<-PureIBNRPrediction$UnearnedDuration/DayPolicyDuration Data$NetEarnedDuration<-1-Data$PredictedUPRYears-Data$PredictedPureIBNRYears+0.00001 #Predict Claim Frequency with reduced exposure to adjust for pure IBNR AdjFreqGLM<-glm(Data$ClaimReportedBeforeValuation~1, offset = log(Data$NetEarnedDuration), family = "poisson") #summary(AdjFreqGLM) PredictedFrequency<-exp(AdjFreqGLM$coefficients) #Compare as if not adjusting for pure IBNR UnAdjFreqGLM<-glm(Data$ClaimReportedBeforeValuation~1, offset = log(Data$EarnedDuration+0.00001), family = "poisson") #summary(UnAdjFreqGLM) cbind(AdjustedModel=round(exp(AdjFreqGLM$coefficients),3),UnAdjustedModel=round(exp(UnAdjFreqGLM$coefficients),3),Actual=ClaimFreq) #Compare UPR Claim predictions cbind(Actual=sum(Data$UPRClaims),Predicted=round(sum(Data$PredictedUPRYears)*PredictedFrequency,0)) #Compare Pure IBNR Claim predictions cbind(Actual=sum(Data$PureIBNRClaims),Predicted=round(sum(Data$PredictedPureIBNRYears)*PredictedFrequency,0), Theoretical=327) #Theoretical comes from re-running with the theoretical parameter values #Compare Pure IBNR predictions with Chain Ladder ClaimsAY1RY1<-sum(Data[Data$AccidentDate<as.Date("2012/1/1") & Data$NotificationDate<as.Date("2012/1/1"), ]$ClaimReportedBeforeValuation) ClaimsAY1RY2<-sum(Data[Data$AccidentDate<as.Date("2012/1/1") & Data$NotificationDate<as.Date("2013/1/1"), ]$ClaimReportedBeforeValuation) ClaimsAY2RY1<-sum(Data$ClaimReportedBeforeValuation)-ClaimsAY1RY1-ClaimsAY1RY2 PredictedPureIBNRClaimsChainLadder<-(ClaimsAY1RY2/ClaimsAY1RY1-1)*ClaimsAY2RY1 cbind(Actual=sum(Data$PureIBNRClaims),PredictedGLM=round(sum(Data$PredictedPureIBNRYears)*PredictedFrequency,0),PredictedCL=round(PredictedPureIBNRClaimsChainLadder,0)) #Note: Chain ladder can predict better if periods brake down to months or quarters #Sliced LogNormal-Pareto claim Severity assumption mu=5.6 sigma=1.65 s=10000 alpha=1.5 #Simulate the reserve distribution numberOfSimulations=10000 SimulatedYears <- numeric(length = numberOfSimulations) ExpectedFrequency <- sum(Data$PredictedPureIBNRYears)*PredictedFrequency for (i in 1 :numberOfSimulations){ SimulatedYears[i]=round(sum(qSlicedLNormPareto(runif(rpois(1,ExpectedFrequency)),mu,sigma,s,alpha)),0) } #Visualising pure IBNR reserve distribution head(SimulatedYears) hist(SimulatedYears, breaks = 1000, xlim = c(0,1.5e6)) summary(SimulatedYears) #Note: the strength of the method illustrated is that it can have GLM factors for reporting delay, frequency and severity, and run the simulations per risk