# install.packages("tidyr")
# install.packages("readxl")
# install.packages("moments")
# install.packages("evd")
# install.packages("actuar")
# install.packages("EnvStats")
# install.packages("melt")
# install.packages("evd")
# install.packages("EnvStats")
library(EnvStats)
##
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
##
## predict, predict.lm
library(tidyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(moments)
##
## Attaching package: 'moments'
## The following objects are masked from 'package:EnvStats':
##
## kurtosis, skewness
library(evd)
library(actuar)
##
## Attaching package: 'actuar'
## The following objects are masked from 'package:evd':
##
## dgumbel, pgumbel, qgumbel, rgumbel
## The following objects are masked from 'package:EnvStats':
##
## dpareto, ppareto, qpareto, rpareto
## The following objects are masked from 'package:stats':
##
## sd, var
## The following object is masked from 'package:grDevices':
##
## cm
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:EnvStats':
##
## boxcox
library(stats)
library(fitdistrplus)
## Loading required package: survival
library(melt)
rm(list = ls())
library(readxl)
path <- "/Users/julie/Desktop/concordia/Y3-winter2024/actu459/Data.xlsx"
Data <- read_excel(path)
View(Data)
#This makes the Exp^2 column consistent, and removes unneeded columns.
####################### Graphs ################
# Initialize Data
#We use average claim amounts since it is more representative.
Data$Avg.Per.Claim <- Data$Exp.Amt/Data$`S.Ad. No.`
ExpAmt<-sort(as.matrix(Data[,1]))
AvgAmt<-sort(as.matrix(Data[,4]))
log_Avg <- log(AvgAmt)
inv_Avg <- 1/AvgAmt
# Histogram
par(mfrow = c(1,2))
hist(ExpAmt, main = "Histogram of Expense Amount", xlab = "Expense Amount")
hist(AvgAmt, main = "Histogram of Average Amount", xlab = "Average Amount")
hist(log_Avg, main = " Histogram of Log-Transformed Average Amount", xlab = "Log-Transformed Average Amount")
hist(inv_Avg, main = "Histogram of Inverse-Transformed Average Amount", xlab = "Inverse-Transformed Average Amount")
# Boxplot
boxplot(ExpAmt, main = "Boxplot of Expense Amount")
boxplot(AvgAmt, main = "Boxplot of Average Amount")
boxplot(log_Avg, main = "Boxplot of Log-Transformed Average Amount")
boxplot(inv_Avg, main = "Boxplot of Inverse-Transformed Average Amount")
# Empirical CDF
par(mfrow = c(1,2))
ECDF_Exp <- ecdf(ExpAmt)
plot(ECDF_Exp, main = "ECDF of Expense Amount")
ECDF_Avg <- ecdf(AvgAmt)
plot(ECDF_Avg, main = "ECDF of Average Amount")
ECDF_log_Avg <- ecdf(log_Avg)
plot(ECDF_log_Avg, main = "ECDF of Log-Transformed Average Amount")
ECDF_inv_Avg <- ecdf(inv_Avg)
plot(ECDF_inv_Avg, main = "ECDF of Inverse-Transformed Average Amount")
# MRL
mrlplot(ExpAmt, main = "MRL of Expense Amount" )
mrlplot(AvgAmt, main = "MRL of Average Amount")
mrlplot(log_Avg, main = "MRL of Log-Transformed Average Amount")
mrlplot(inv_Avg, main = "MRL of Inverse-Transformed Average Amount")
# LEV Functions
LEV_exp_amt <- elev(Data$Exp.Amt)
plot(LEV_exp_amt, main = "LEV of Expense Amount")
LEV_Avg <- elev(AvgAmt)
plot(LEV_Avg, main = "LEV of Average Amount")
LEV_log_avg <- elev(log_Avg)
plot(LEV_log_avg, main = "LEV of Log-Transformed Average Amount")
LEV_inv_avg <- elev(inv_Avg)
plot(LEV_inv_avg, main = "LEV of Inverse-Transformed Average Amount")
LEV_exp_amt <- elev(Data$Exp.Amt)
plot(LEV_exp_amt)
# Normal QQ plot
qqnorm(ExpAmt, main = "Normal Q-Q Plot of expense Amount")
qqnorm(AvgAmt, main = "Normal Q-Q Plot of Average Amount")
qqnorm(log_Avg, main = "Normal Q-Q Plot of Log-Transformed Average Amount")
qqline(log_Avg, col = "blue")
qqnorm(inv_Avg, main = "Normal Q-Q Plot of Inverse-Transformed Amount")
qqline(inv_Avg, col = "blue")
################## Estimation ############################
N <- 614
#iid assumption verification
par(mfrow=c(1,1))
acf(Data$Avg.Per.Claim)
#Gives how many sd away from mean is largest observation.
(max(Data$Avg.Per.Claim)-mean(Data$Avg.Per.Claim))/sd(Data$Avg.Per.Claim)
## [1] 15.19699
#Do we jedd other ACFs?
acf(Data$Exp.Amt)
acf(Data$Exp2)
acf(Data$Avg.Per.Claim)
#Percentiles
p <- c(seq(0,1,0.25),0.95,0.99)
Q <- quantile(Data$Avg.Per.Claim,p)
log(Q)
## 0% 25% 50% 75% 100% 95% 99%
## 1.704748 2.843595 3.158453 3.671063 10.003559 4.798863 7.209141
#Add coefficient of variation figure.
Avg.Claim.Out<-Data$Avg.Per.Claim[!Data$Avg.Per.Claim %in% boxplot.stats(Data$Avg.Per.Claim)$out]
length(Avg.Claim.Out)
## [1] 548
#Number of outliers:
Nb.Out = 614 - length(Avg.Claim.Out)
#Percentage of outliers
Percent.Out = Nb.Out/N
## For censoring part:
#First, locate the smallest outlier:
CensPoint <- quantile(Data$Avg.Per.Claim[Data$Avg.Per.Claim %in% boxplot.stats(Data$Avg.Per.Claim)$out],0.5)
#72.90 ; use this as censoring point:
right <- ifelse(Data$Avg.Per.Claim<=CensPoint,Data$Avg.Per.Claim,CensPoint)
left <- ifelse(Data$Avg.Per.Claim<CensPoint,Data$Avg.Per.Claim,CensPoint)
censdata <- data.frame(left,right)
#### FITTING A BETA DISTRIBUTION ####
#Cullen and Frey : guides choice of distribution...
descdist(Data$Avg.Per.Claim, discrete = FALSE) #including outliers.
## summary statistics
## ------
## min: 5.5 max: 22105
## median: 23.53416
## mean: 173.527
## estimated sd: 1443.146
## estimated skewness: 13.13156
## estimated kurtosis: 187.4369
#keeping 95% of dataset:
Avg.93 <- Data$Avg.Per.Claim[Data$Avg.Per.Claim<quantile(Data$Avg.Per.Claim, 0.93)]
descdist(Avg.93, discret=FALSE)
## summary statistics
## ------
## min: 5.5 max: 103.1658
## median: 22.58293
## mean: 28.7032
## estimated sd: 18.15888
## estimated skewness: 1.693746
## estimated kurtosis: 5.657583
descdist(Avg.Claim.Out)
## summary statistics
## ------
## min: 5.5 max: 72.56949
## median: 22.15154
## mean: 26.3284
## estimated sd: 14.15446
## estimated skewness: 1.357088
## estimated kurtosis: 4.363216
BetaInitial <- fitdist(1/Data$Avg.Per.Claim,"beta", method='mle')
plot(BetaInitial,las=1)
#P-P plot looks pretty good as well as ecdf vs fitted cdf. Beta shows potential.
gofstat(BetaInitial)
## Goodness-of-fit statistics
## 1-mle-beta
## Kolmogorov-Smirnov statistic 0.0921204
## Cramer-von Mises statistic 1.5949580
## Anderson-Darling statistic 8.8699446
##
## Goodness-of-fit criteria
## 1-mle-beta
## Akaike's Information Criterion -2761.847
## Bayesian Information Criterion -2753.007
vcov(BetaInitial)
## shape1 shape2
## shape1 0.0107190 0.2320991
## shape2 0.2320991 6.5881527
solve(vcov(BetaInitial))
## shape1 shape2
## shape1 393.35995 -13.857980
## shape2 -13.85798 0.640001
#confreg(BetaInitial,level=0.95)
#Already not too bad of a fit on data that has not been transformed...
#Let's try keeping vlaues up to 95th percentile
Beta.93 <- fitdist(1/Avg.93,"beta",method="mle")
plot(Beta.93)
gofstat(Beta.93)
## Goodness-of-fit statistics
## 1-mle-beta
## Kolmogorov-Smirnov statistic 0.03328026
## Cramer-von Mises statistic 0.12132030
## Anderson-Darling statistic 0.88447058
##
## Goodness-of-fit criteria
## 1-mle-beta
## Akaike's Information Criterion -2748.344
## Bayesian Information Criterion -2739.649
vcov(Beta.93)
## shape1 shape2
## shape1 0.04361505 0.883117
## shape2 0.88311699 20.656072
solve(vcov(Beta.93))
## shape1 shape2
## shape1 170.683127 -7.2972815
## shape2 -7.297282 0.3603954
#If we remove all outliers:
Beta.Out <- fitdist(1/Avg.Claim.Out,"beta",method="mle")
plot(Beta.Out)
#Way better test-stats...
gofstat(Beta.Out)
## Goodness-of-fit statistics
## 1-mle-beta
## Kolmogorov-Smirnov statistic 0.02174861
## Cramer-von Mises statistic 0.05312875
## Anderson-Darling statistic 0.55313762
##
## Goodness-of-fit criteria
## 1-mle-beta
## Akaike's Information Criterion -2676.457
## Bayesian Information Criterion -2667.844
vcov(Beta.Out)
## shape1 shape2
## shape1 0.0634575 1.244377
## shape2 1.2443769 27.591884
solve(vcov(Beta.Out))
## shape1 shape2
## shape1 136.297856 -6.1469492
## shape2 -6.146949 0.3134662
#Censored beta:
Beta.Cens <- fitdistcens(1/censdata,"beta")
plot(Beta.Cens)
#goftest(Beta.Cens$censdata[,1],x=1/Data$Avg.Per.Claim,distribution="beta","ks")
#Comparing parameters:
par0 <- BetaInitial$estimate
par1 <- Beta.93$estimate
par2 <- Beta.Out$estimate
par3 <- Beta.Cens$estimate
x <- seq(0,1,length.out=512)
y0 <- dbeta(x,par0[1],par0[2])
y1 <- dbeta(x,par1[1],par1[2])
y2 <- dbeta(x,par2[1],par2[2])
y3 <- dbeta(x,par3[1],par3[2])
#KS stat for beta using y3 as defined below::
max(y3-density(1/censdata[,1])$y) #does not work
## [1] 17.22192
#Variance covariance
vcov(Beta.Cens)
## shape1 shape2
## shape1 0.02198489 0.4729226
## shape2 0.47292260 12.3398111
#Fisher information
solve(vcov(Beta.Cens))
## shape1 shape2
## shape1 259.057084 -9.928349
## shape2 -9.928349 0.461542
plot(density(1/Avg.93),
main="Four Proposed Beta distributions",
type="h",col="lightgray",
xlim=c(0,0.2),
ylim=c(0,20),
xlab="",
ylab="y")
par(new=TRUE)
plot(x,y0,axes=FALSE,type="l", col="red",
xlim=c(0,.20),
ylim=c(0,20),
ylab="")
par(new=TRUE)
plot(x,y1,axes=FALSE,type="l", col="green",
xlim=c(0,0.2),
ylim=c(0,20),
ylab="")
par(new=TRUE)
plot(x,y2,axes=FALSE,type="l", col="blue",
xlim=c(0,0.2),
ylim=c(0,20),ylab="")
par(new=TRUE)
plot(x,y3,axes=FALSE,type="l", col="orange",
xlim=c(0,0.2),
ylim=c(0,20),ylab="")
grid()
c <- 0.95
confint(BetaInitial,level=c)
## 2.5 % 97.5 %
## shape1 1.753044 2.158885
## shape2 37.818040 47.879475
confint(Beta.93,level=c)
## 2.5 % 97.5 %
## shape1 3.273351 4.091997
## shape2 66.156747 83.972409
confint(Beta.Out,level=c)
## 2.5 % 97.5 %
## shape1 3.831968 4.819428
## shape2 75.027989 95.618578
confint(Beta.Cens,level=c)
## 2.5 % 97.5 %
## shape1 2.456661 3.03788
## shape2 52.709520 66.47947
#### FITTING A LOGNORMAL
#Fit on raw data.
LnormFit <- fitdist(Data$Avg.Per.Claim, method="mle","lnorm")
plot(LnormFit,las=1)
gofstat(LnormFit)
## Goodness-of-fit statistics
## 1-mle-lnorm
## Kolmogorov-Smirnov statistic 0.1511205
## Cramer-von Mises statistic 5.3359125
## Anderson-Darling statistic 31.0964268
##
## Goodness-of-fit criteria
## 1-mle-lnorm
## Akaike's Information Criterion 5896.082
## Bayesian Information Criterion 5904.922
#No need for GOF to see it is not appropriate...
#Covariance matrix + Fisher Info
vcov(LnormFit)
## meanlog sdlog
## meanlog 0.001571764 0.0000000000
## sdlog 0.000000000 0.0007858747
#Fisher information is the inverse of the above matrix:
solve(vcov(LnormFit))
## meanlog sdlog
## meanlog 636.2278 0.000
## sdlog 0.0000 1272.467
#With data up to 95 percentile:
Lnorm.93 <- fitdist(Avg.93, method="mle","lnorm")
plot(Lnorm.93,las=1)
gofstat(Lnorm.93)
## Goodness-of-fit statistics
## 1-mle-lnorm
## Kolmogorov-Smirnov statistic 0.07064631
## Cramer-von Mises statistic 0.71271113
## Anderson-Darling statistic 4.31037386
##
## Goodness-of-fit criteria
## 1-mle-lnorm
## Akaike's Information Criterion 4581.808
## Bayesian Information Criterion 4590.503
#Covariance matrix + Fisher Info
vcov(Lnorm.93)
## meanlog sdlog
## meanlog 0.0005187194 0.0000000000
## sdlog 0.0000000000 0.0002593518
#Fisher information is the inverse of the above matrix:
solve(vcov(Lnorm.93))
## meanlog sdlog
## meanlog 1927.824 0.000
## sdlog 0.000 3855.766
#Same as beta, now fit without any outliers:
Lnorm.Out <- fitdist(Avg.Claim.Out, method="mle","lnorm")
plot(Lnorm.Out,las=1)
gofstat(Lnorm.Out)
## Goodness-of-fit statistics
## 1-mle-lnorm
## Kolmogorov-Smirnov statistic 0.05241185
## Cramer-von Mises statistic 0.34180086
## Anderson-Darling statistic 2.23548474
##
## Goodness-of-fit criteria
## 1-mle-lnorm
## Akaike's Information Criterion 4226.394
## Bayesian Information Criterion 4235.007
#Covariance matrix + Fisher Info
vcov(Lnorm.Out)
## meanlog sdlog
## meanlog 0.0004391117 0.0000000000
## sdlog 0.0000000000 0.0002195476
#Fisher information is the inverse of the above matrix:
solve(vcov(Lnorm.Out))
## meanlog sdlog
## meanlog 2277.325 0.00
## sdlog 0.000 4554.82
#Censored Lnorm at 50p outliers
Lnorm.Cens <- fitdistcens(censdata,"lnorm")
Est<-Lnorm.Cens$estimate
plot(Lnorm.Cens)
Lnorm.Cens$aic
## [1] 5293.267
#KS stat for censored:
max(y3-density(censdata[,1])$y)
## [1] 17.91717
#Covariance matrix + Fisher Info
vcov(Lnorm.Cens)
## meanlog sdlog
## meanlog 0.0007050687 0.000000000
## sdlog 0.0000000000 0.000352527
#Fisher information is the inverse of the above matrix:
solve(vcov(Lnorm.Cens))
## meanlog sdlog
## meanlog 1418.302 0.000
## sdlog 0.000 2836.662
c <- 0.95
confint(LnormFit,level=c)
## 2.5 % 97.5 %
## meanlog 3.3192513 3.474659
## sdlog 0.9274317 1.037321
confint(Lnorm.93,level=c)
## 2.5 % 97.5 %
## meanlog 3.1533900 3.2426680
## sdlog 0.5126682 0.5757963
confint(Lnorm.Cens,level=c)
## 2.5 % 97.5 %
## meanlog 3.2548490 3.3589353
## sdlog 0.6211609 0.6947603
confint(Lnorm.Out,level=c)
## 2.5 % 97.5 %
## meanlog 3.1047805 3.1869225
## sdlog 0.4615027 0.5195848
#Comparing parameters:
par0 <- LnormFit$estimate
par1 <- Lnorm.93$estimate
par2 <- Lnorm.Out$estimate
par3 <- Lnorm.Cens$estimate
#All 3 lnorm models models on same graph
x <- seq(0,150,length.out=512)
y0 <- dlnorm(x,meanlog=par0[1],sdlog=par0[2])
y1 <- dlnorm(x,meanlog=par1[1],sdlog=par1[2])
y2 <- dlnorm(x,meanlog=par2[1],sdlog=par2[2])
y3 <- dlnorm(x,meanlog=par3[1],sdlog=par3[2])
#KS stat for
max(y3-density(censdata[,1])$y)
## [1] 0.02013393
plot(density(Avg.93),
main="Lognormal - Censored Data",
type="h",col="lightgray",
xlim=c(0,130),
ylim=c(0,0.05),
xlab="",
ylab="y")
par(new=TRUE)
plot(x,y0,axes=FALSE,type="l", col="red",
xlim=c(0,130),
ylim=c(0,0.05),
ylab="")
par(new=TRUE)
plot(x,y1,axes=FALSE,type="l", col="green",
xlim=c(0,130),
ylim=c(0,0.05),
ylab="")
par(new=TRUE)
plot(x,y2,axes=FALSE,type="l", col="blue",
xlim=c(0,130),
ylim=c(0,0.05),ylab="")
par(new=TRUE)
plot(x,y3,axes=FALSE,type="l", col="orange",
xlim=c(0,130),
ylim=c(0,0.05),ylab="")
grid()
#Green and orange lines are essentially the same. So we keep the model that removes the less outliers.
#### FITTING AN INVERSE GAMMA
#Fit on raw data.
InvGammaFit <- fitdist(Data$Avg.Per.Claim, method="mle","invgamma")
plot(InvGammaFit,las=1)
gofstat(InvGammaFit)
## Goodness-of-fit statistics
## 1-mle-invgamma
## Kolmogorov-Smirnov statistic 0.09530957
## Cramer-von Mises statistic 1.72230201
## Anderson-Darling statistic 9.59005987
##
## Goodness-of-fit criteria
## 1-mle-invgamma
## Akaike's Information Criterion 5590.212
## Bayesian Information Criterion 5599.052
#Covariance matrix + Fisher Info
vcov(InvGammaFit)
## shape scale
## shape 0.01134174 0.2588032
## scale 0.25880320 7.6089949
#Fisher information is the inverse of the above matrix:
solve(vcov(InvGammaFit))
## shape scale
## shape 393.83689 -13.3954946
## scale -13.39549 0.5870417
#With data up to 95 percentile:
InvGamma.93 <- fitdist(Avg.93, method="mle","invgamma")
plot(InvGamma.93,las=1)
gofstat(InvGamma.93)
## Goodness-of-fit statistics
## 1-mle-invgamma
## Kolmogorov-Smirnov statistic 0.03606245
## Cramer-von Mises statistic 0.14608910
## Anderson-Darling statistic 1.02460964
##
## Goodness-of-fit criteria
## 1-mle-invgamma
## Akaike's Information Criterion 4556.263
## Bayesian Information Criterion 4564.958
#Covariance matrix + Fisher Info
vcov(InvGamma.93)
## shape scale
## shape 0.04781335 1.022359
## scale 1.02235887 24.942577
#Fisher information is the inverse of the above matrix:
solve(vcov(InvGamma.93))
## shape scale
## shape 169.25087 -6.9373397
## scale -6.93734 0.3244432
#Same as beta, now fit without any outliers:
InvGamma.Out <- fitdist(Avg.Claim.Out, method="mle","invgamma")
plot(InvGamma.Out,las=1)
gofstat(InvGamma.Out)
## Goodness-of-fit statistics
## 1-mle-invgamma
## Kolmogorov-Smirnov statistic 0.02085832
## Cramer-von Mises statistic 0.05651847
## Anderson-Darling statistic 0.57009807
##
## Goodness-of-fit criteria
## 1-mle-invgamma
## Akaike's Information Criterion 4218.222
## Bayesian Information Criterion 4226.834
#Covariance matrix + Fisher Info
vcov(InvGamma.Out)
## shape scale
## shape 0.07001513 1.451382
## scale 1.45138187 33.643844
#Fisher information is the inverse of the above matrix:
solve(vcov(InvGamma.Out))
## shape scale
## shape 135.078486 -5.8272314
## scale -5.827231 0.2811075
#Censored Inv Gamma at 50p outliers
InvGamma.Cens <- fitdistcens(censdata,"invgamma")
Est<-InvGamma.Cens$estimate
plot(InvGamma.Cens)
#KS stat for censored:
-max(y3-density(1/censdata[,1])$y)
## [1] -4.105395e-05
#Covariance matrix + Fisher Info
vcov(InvGamma.Cens)
## shape scale
## shape 0.02373708 0.5382852
## scale 0.53828522 14.5933317
#Fisher information is the inverse of the above matrix:
solve(vcov(InvGamma.Cens))
## shape scale
## shape 257.595713 -9.5015976
## scale -9.501598 0.4189975
c <- 0.95
confint(InvGammaFit,level=c)
## 2.5 % 97.5 %
## shape 1.799991 2.217454
## scale 40.429862 51.242751
confint(InvGamma.93,level=c)
## 2.5 % 97.5 %
## shape 3.420791 4.277933
## scale 72.519648 92.096766
confint(InvGamma.Cens,level=c)
## 2.5 % 97.5 %
## shape 2.547649 3.151586
## scale 57.133411 72.108014
confint(InvGamma.Out,level=c)
## 2.5 % 97.5 %
## shape 4.017964 5.055192
## scale 82.672790 105.409671
par(mfrow=c(1,1))
IGammaInitial <- fitdist(Data$Avg.Per.Claim,"invgamma",method="mle")
IGamma.93 <- fitdist(Avg.93,"invgamma",method="mle")
IGamma.Out <- fitdist(Avg.Claim.Out,"invgamma",method="mle")
IGamma.Cens <- fitdistcens(censdata, "invgamma")
#Comparing parameters:
par0 <- IGammaInitial$estimate
par1 <- IGamma.93$estimate
par2 <- IGamma.Out$estimate
par3 <- IGamma.Cens$estimate
#All 3 lnorm models models on same graph
x <- seq(0,150,length.out=1000)
y0 <- dinvgamma(x,par0[1],1/par0[2])
y1 <- dinvgamma(x,par1[1],1/par1[2])
y2 <- dinvgamma(x,par2[1],1/par2[2])
y3 <- dinvgamma(x,par3[1],1/par3[2])
plot(density(Avg.93),
main="Four Proposed Inverse Gamma distributions",
type="h",col="lightgray",
xlim=c(0,130),
ylim=c(0,0.05),
xlab="",
ylab="y")
par(new=TRUE)
plot(x,y0,axes=FALSE,type="l", col="red",
xlim=c(0,130),
ylim=c(0,0.05),
ylab="")
par(new=TRUE)
plot(x,y1,axes=FALSE,type="l", col="green",
xlim=c(0,130),
ylim=c(0,0.05),
ylab="")
par(new=TRUE)
plot(x,y2,axes=FALSE,type="l", col="blue",
xlim=c(0,130),
ylim=c(0,0.05),ylab="")
par(new=TRUE)
plot(x,y3,axes=FALSE,type="l", col="orange",
xlim=c(0,130),
ylim=c(0,0.05),ylab="")
grid()
#Green and orange lines are essentially the same. So we keep the model that removes the less outliers.