R codes for Tweedie Compound Poisson Gamma

Question

I found the following R codes to fit a Tweedie Compound Poisson Gamma distribution. I have to fit it to my 399 claim amounts. I have seen the following R codes ptweedie.series(q, power, mu, phi) and dtweedie.series(y, power, mu, phi). However I fail to understand the codes fully and after importing my data into R, how to proceed? Thanks in advance.

Answer 1

First a note: importing your dataset from the comments above yielded 398 claims, not 399. One of these was 4 orders of magnitude larger than the median claim. So I suspect a typo. In the analysis that follows I excluded that sample, leaving 397.

A quick look at the Wikipedia entry for Tweedie Distributions reveals that this is actually a family of exponential distributions distinguished by the power parameter (xi in the R documentation). Power=1 yields the Poisson distribution, power=2 yields the Gamma distribution, power=3 yields the inverse Gaussian distribution, and so on. The Tweedie distributions are also defined for non-integer power. The parameter mu is the mean, and phi is a dispersion parameter, related to variance.

So the basic question, as I understand it, is which combination of power, mu, and phi yield a distribution which best fits your claims data?

One way to assess whether a distribution fits a sample is the Q-Q plot. This plots quantiles of your sample vs. quantiles of the test distribution. If the sample is distributed as the test distribution, the Q-Q plot should be a straight line. In R code (and with X as your vector of samples):

summary(X)        # NOTE: max/median > 1e4 !!!
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 1.00e+03 5.50e+03 1.20e+04 5.47e+05 2.50e+04 2.08e+08 
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)

Both the Gamma and the Inverse Gaussian distributions explain your 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 your claims data seems to follow a tweedie distribution with power=2.5. The next step is to 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))
opt$solution
# [1] 1.884839e+04 9.735325e-03

Finally, we confirm that the solution does indeed fit the data well.

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")