x <- c(6500,16600,14500,27800,9900,6000,27400,20600,3500,1150,10280,12150,1563162,3500,5500,6000,16150,4500,20000,51540,12500,141811,784760,43000,15600,26840,3000,115100,15500,12700,49200,6300,37500,4000,4600,8500,3000,7700,14600,78400,59700,20600,5000,3500,50900,32000,28850,8000,7400,59750,75900,3000,66900,104500,50300,38750,17500,66400,17670,3949,11150,9653,5000,3000,3000,14050,3000,25000,4000,4500,4750,55000,2000,5500,4750,5000,4500,3000,5000,280000,9500,4500,4500,30000,1500,8000,4500,4800,3000,8000,2500,30000,12000,3500,9600,16000,4000,5500,4000,5000,6000,1200,1000,6200,6300,5500,5000,3500,3500,1250,1300,4600,4000,3600,5000,17300,6000,8000,6000,5500,13500,2500,5500,4500,5900,7000,13000,6700,6800,5000,6200,6300,6100,3600,3300,3500,2900,3000,3150,3250,9900,8700,7300,9600,4000,4500,4100,3350,3950,4250,3850,3950,3500,4000,4550,3000,13300,13900,30000,4550,14650,14400,3500,13550,13900,3900,11000,12000,28000,9000,22000,29000,5500,3000,2500,208009850,18000,27500,20000,28000,22500,20500,27600,28500,18000,18500,18900,11700,3500,3900,20800,21000,5000,8600,10500,10600,9000,9500,12000,12500,11800,7500,9750,9000,8000,9800,5600,9700,75503500,6000,11550,10500,5000,12100,7000,7000,8000,8500,9500,5000,4600,4550,3700,2500,2900,5500,5850,3550,2650,4550,6500,7500,6850,8000,10250,6000,7750,6950,6900,7800,7250,25000,28000,19150,32000,28000,27000,19000,19550,16350,20150,20800,12000,25500,26800,29550,24350,13000,25000,49500,44000,27000,25800,9000,10000,8550,8000,9500,8100,11000,11000,7350,7550,4050,14000,14500,13000,15000,19500,20000,25000,30000,21000,39500,43500,21500,8900,18500,18000,19500,8700,19500,6300,12000,15000,18000,40000,36000,25000,39000,37000,30000,33000,37000,25000,18000,20000,21000,19800,35700,60000,19850,25890,52000,15000, 17000,19850,8000,16500,8500,11000,10000,12000,13500,18000,6500,6850,25000,22000,32550,19000,11000,23000,21550,22000,23550,30550, 35000,9800,23000,26800,23550,23850,15000,26700,26850,32000,70000,73000,71000,29000,35500,76550,78000,50000,46000,47500,57500,35000,36000,58500,57000,50100,39000,60550,48100,56550,3000,4200,3500,2700,11000,12500,12000,8790,14600,6700,10000,19300,12000,12500,20000,25990,27340,20000,15000,64000,25000,29870,32340,39000,32000,19000,15600,14000,28500,34000 )
summary(x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1000 5500 12000 738568 25000 208009850
X <- x[x<max(x)] ##Remove largest value erroneous
hist(X,breaks=c(seq(1,1e5,1000),Inf),xlim=c(0,100000))
library(tweedie)
qqTweedie <- function(xi,p,mu,phi) {
names <- c("Poisson","Gamma","Inverse Gaussian","Positive Stable")
plot(qtweedie(p,xi,mu,phi),quantile(X,probs=p),
main=paste0("Power = ",xi," (",names[xi],")"))
qqline(X,prob=c(0.25,0.75), col="blue", lty=2,
distribution=function(p) qtweedie(p,xi,mu,phi))
}
p <- seq(0.02,0.98,length=100)
par(mfrow=c(2,2))
lapply(c(1:4),qqTweedie,p=p,mu=1,phi=1)
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
## Below we can see the Gamma and the Inverse Gaussian distributions explain the data up to claims of ~40,000. The Gamma distribution underestimates the frequency of larger claims, while the Inverse Gaussian distribution overestimates their frequency. So let's try power=2.5.
par(mfrow=c(1,1))
xi <- 2.5
plot(qtweedie(p,xi,1,1),quantile(X,probs=p),main=paste0("Power = ",xi))
qqline(X,prob=c(0.25,0.75), col="blue", lty=2,
distribution=function(p) qtweedie(p,xi,1,1))
##So claims data seems to follow a tweedie distribution with power=2.5.
##Estimate mu and phi, given power=2.5. This is a non-linear optimization problem in 2 dimensions, so we use package nloptr. It turns out that convergence depends on having starting parameters relatively close the the optimal values, so there is a fair amount of trial and error to get nlopt(...) to converge.
library(nloptr)
F <- function(params){ # Note: xi, Q, and p are defined external to F
mu <- params[1]
phi <- params[2]
return(sum(Q - qtweedie(p,xi,mu,phi))^2)
}
xi <- 2.5
Q <- quantile(X,p)
opt <- nloptr(x0=c(mu=1e4,phi=.01), eval_f=F, ub=c(5e4,.1), lb = c(1,0),
opts = list(algorithm="NLOPT_LN_COBYLA",maxeval=1e3,print_level=1))
## iteration: 1
## f(x) = 669610961370.851562
## iteration: 2
## f(x) = 18466377580.265007
## iteration: 3
## f(x) = 41423048609.784180
## iteration: 4
## f(x) = 289044517444.089600
## iteration: 5
## f(x) = 40528335758.356323
## iteration: 6
## f(x) = 77907814515.691422
## iteration: 7
## f(x) = 25827259096.741714
## iteration: 8
## f(x) = 2647843852.953338
## iteration: 9
## f(x) = 1102070017.938444
## iteration: 10
## f(x) = 959253309.772546
## iteration: 11
## f(x) = 10708481380.087471
## iteration: 12
## f(x) = 3241876263.611691
## iteration: 13
## f(x) = 2076747094.759489
## iteration: 14
## f(x) = 105223880.551220
## iteration: 15
## f(x) = 110202637.786527
## iteration: 16
## f(x) = 13736407.511191
## iteration: 17
## f(x) = 88679.840348
## iteration: 18
## f(x) = 80230239.862883
## iteration: 19
## f(x) = 106066866.552006
## iteration: 20
## f(x) = 27523407.562641
## iteration: 21
## f(x) = 1187009.681452
## iteration: 22
## f(x) = 5424908.978762
## iteration: 23
## f(x) = 2176134.012825
## iteration: 24
## f(x) = 600365.346888
## iteration: 25
## f(x) = 134722.833537
## iteration: 26
## f(x) = 119371.921899
## iteration: 27
## f(x) = 139288.218873
## iteration: 28
## f(x) = 17223.750400
## iteration: 29
## f(x) = 996.645319
## iteration: 30
## f(x) = 109.978934
## iteration: 31
## f(x) = 30381.880766
## iteration: 32
## f(x) = 3297.522489
## iteration: 33
## f(x) = 86.773409
## iteration: 34
## f(x) = 5940.051104
## iteration: 35
## f(x) = 251.790211
## iteration: 36
## f(x) = 5145.746167
opt$solution
## [1] 1.892641e+04 9.279728e-03
mu <- opt$solution[1]
phi <- opt$solution[2]
par(mfrow=c(1,1))
hist(X,breaks=c(seq(1,1e5,1000),Inf),xlim=c(0,1e5))
x <- seq(1,1e5,1e3)
lines(x,dtweedie(x,xi,mu,phi),col="red")
