R – How to Setup Panjer’s Recursion Correctly in Actuarial Science

actuarial-sciencebinningrself-study

I have a table of $k=(0,1,2,3,4,5,6)$ and $number=(40544,8082,1205,145,20,3,1)$

I need to fit data by a Compound Poisson-Gamma distribution and then make a discretization and compare results with observations.

My attemp is:

The compound Poisson-Gamma distribution can be realized as Negative Binomial distibution:

library(actuar)
library(fitdistrplus)

number <- c(rep(0, 40544), rep(1, 8082), rep(2, 1205), rep(3, 145), rep(4, 20), rep(5, 3), rep(6, 1))

fitnb <- fitdist(number, "nbinom")
shape = fitnb$estimate[1] # size = 2.080522         
lambda = fitnb$estimate[2] # mu = 0.2205776

As the result I know the lambda and shape parameters for the discretisation:

n = 6 # we have values from 0 to 6
fx <- discretize(pgamma(x, shape = shape, 1), 
                 method = "unbiased", from = 0, to = n, step = 1.0, 
                 lev = levgamma(x, shape = shape, 1))
sum(fx) # 0.9804088 < 1

Following the documentation and answer I appied the aggregateDist() function with changed parameters lambda and convolve:

The recursion will fail to start if the expected number of claims is
too large. One may divide the appropriate parameter of the frequency
distribution by $2^n$ and convolve the resulting distribution $n$ =
convolve times.

Fsc <- aggregateDist("recursive", 
                     model.freq = "poisson", 
                     model.sev = fx, 
                     lambda = lambda/(2^n), 
                     convolve = n, x.scale = 1)

But I seen again the

Warning message:
In panjer(fx = model.sev, dist = dist, p0 = p0, x.scale = x.scale,  :
  maximum number of recursions reached before the probability distribution was complete

I was surprised on the output:

summary(Fsc)
#    Aggregate Claim Amount Empirical CDF:
#            Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
#    0.000000e+00 0.000000e+00 0.000000e+00 4.261483e-01 0.000000e+00 3.200000e+04

Finally, I have tried to compare obtained results (black) with original data (red):

plot(Fsc, do.points = FALSE, verticals = TRUE, xlim = c(0, n))
femp <- c(40544, 8082, 1205, 145, 20, 3, 1)/50000
sum(femp) # 1
plot(stepfun(0:6, diffinv(femp)), pch = 19, col = "red", add = TRUE)

Edit. After the jbowman's comment I have tried to use the "negative binomial":

# https://rdrr.io/cran/actuar/src/R/panjer.R
# r <- size     # size     
# p <- 1/(1+mu) # prob

Fsc_ng <- aggregateDist("recursive", 
                         model.freq = "negative binomial", 
                         model.sev = fx, 
                         size = shape, prob = 1/(1+lambda), echo = TRUE, maxit = 10)

# x Pr[S = x]   Cumulative probability
# 0 0.68415174  0.68415174
# 1 0.083957895 0.76810964
# 2 0.079291596 0.84740123
# 3 0.055756862 0.9031581
# 4 0.036292288 0.93945038
# 5 0.022994652 0.96244504
# 6 0.01271653  0.97516157
# 7 0.0066424692    0.98180404
# 8 0.0040147675    0.9858188
# 9 0.0023141156    0.98813292
# 10    0.0013060718    0.98943899

Question. How to setup the parameters of aggregateDist() function?

Best Answer

There is no need to discretize a gamma distribution; the negative binomial distribution is already discrete. Once you have the parameter estimates for the negative binomial, you are done, except for the plotting:

size <- 2.080522
mu <- 0.2205776

est_cdf <- pnbinom(0:6, mu=mu, size=size)
emp_cdf <- cumsum(c(40544, 8082, 1205, 145, 20, 3, 1)) / 50000

plot(emp_cdf ~ c(0:6), type = 'b', pch=17, xlab="X values", ylab="CDF")
lines(est_cdf ~ c(0:6), type = 'b', lty=2, pch=16, col=2)
legend("topleft", col=c(1,2), pch=c(17,16),
       legend = c("Empirical CDF",
                  "NB Estimated CDF")
)

which generates:

showing an almost perfect fit, with the red lines and points on top of the black ones. For a numerical comparison of the negative binomial fit (est_cdf) with the data (emp_cdf), we get:

est_cdf
[1] 0.8108675 0.9725817 0.9964581 0.9995712 0.9999502 0.9999944 0.9999994
emp_cdf
[1] 0.81088 0.97252 0.99662 0.99952 0.99992 0.99998 1.00000

Related Solutions

Solved – Negative Binomial “Process”

Several stochastic processes lead to marginal counts having a Negative Binomial (NB) distribution and can therefore be called NB processes. Among them, the NB Lévy Process is of special interest since increments (counts) over non-overlapping time intervals are independent, a property shared with the Poisson Process, a Gamma process and the Wiener Process. The count $N_t$ on an interval of length $t$ has the NB distribution $$ N_t \sim \textrm{NB}(r,\,p), \quad r = \gamma t $$ so the process depends on the two parameters $\gamma >0$ (with the dimension of an inverse time) and the probability $p$ ($0 < p < 1$). The expectation is proportional to the interval length, and so is its variance $$ \mathbb{E}(N_t) = \gamma t \, (1-p)/p \qquad \textrm{Var}(N_t) = \gamma t \, (1-p)/p^2. $$ The variance is greater than the mean (overdispersion), and the index of dispersion $\textrm{Var}(N_t)/\mathbb{E}(N_t) = 1/p$ does not depend on $t$. When $p$ is close to $1$ and $\gamma (1-p)$ is close to $\lambda >0$, the process behaves like a Poisson Process with rate $\lambda$. An explanation for overdispersion is that several events can happen at the same time, so a small interval can contain more than one event.

It is easy to fit such a process by Maximum Likelihood when the intervals have different lengths. In this case we face a NB regression with a link function differing from the default link in NB GLMs. A special likelihood maximisation is useful.

The article by T.J. Kozubowski and K. Podgorski provide theoretical results as well as an illustration.

Curiously enough, this process does not seem to be frequently used as such by statisticians.

Solved – Quantiles of a compound gamma/negative binomial distribution

As a practical answer to the real questions you're addressing, such high quantiles will generally be quite sensitive to issues with model choice (especially such things as whether you model the right censoring and how heavy the tails are in the components).

But in any case - especially when dealing with high quantiles where ordinary simulation becomes impractical - that has great value in its own right; it's an interesting question from both theoretical and practical standpoints.

A couple of other approaches to this problem are (1) using the Fast Fourier Transform and (2) direct numerical integration.

One useful reference on this topic is Luo and Shevchenko (2009)$^{[1]}$.

In it they develop an adaptive direct numerical integration approach that's faster than simulation and competitive with FFT.

The more traditional approach in actuarial work was been (3) Panjer recursion, which can be found in numerous texts. Embrechts and Frei (2009)$^{[2]}$ discuss and compare Panjer recursion and FFT. (Note that both of these techniques involve discretization of the continuous distribution.)

On the other hand, doing a very unsophisticated version of simulation, and with no effort whatever made to be efficient, generating from a compound gamma-negative binomial isn't particularly onerous. This is timing on my kids' little laptop:

system.time(replicate(100000,sum(rgamma(MASS:::rnegbin(1,4,2),5,.1))))
   user  system elapsed 
   2.82    0.00    2.84

I think 2.8 seconds to generate 100$^{\,}$000 simulations of a compound distribution on a slow little laptop really isn't bad. With some effort to be efficient (of which one might suggest many possibilities), I imagine that could be made a good deal faster.

Here's the ecdf for $10^6$ simulations (which took about 29 seconds):

$\hspace{1cm}$ enter image description here

We see the characteristic discrete jump at zero you expect to see with a compound distribution.

[While it should be easy to make simulation a lot faster, all three alternatives mentioned here - if carried out sensibly - should be a lot faster still.]

You should note that the actuar package supports computation with compound distributions, and offers several methods for calculation with them.

See, for example, this vignette which discusses this facility.

[Of possibly some further passing interest, note that there is an R package for the Poisson-lognormal distribution -- poilog; if you need that distribution at some point it may be useful.]

Added in edit:

A potential quick approximation where the gamma shape parameter isn't changing -

In the gamma case, because a convolution of gammas with constant shape parameter is another gamma, you could write down the distribution of $Y|N=n$, and then evaluate the cdf and the density at a large number of grid-values at each $n$, then simply accumulate the sum directly (rather as one would for a KDE). The direct calculation only yields a lower bound to the true quantile, but if the negative binomial is not heavy tailed it should be quite rapid.

References:

[1]: Luo, X. and Shevchenko, P.V. (2009),
"Computing Tails of Compound Distributions Using Direct Numerical Integration,"
Journal of Computational Finance, 13 (2), 73-111.
[arXiv preprint available here]

[2]: Embrechts, P., and Frei, M. (2009),
"Panjer recursion versus FFT for compound distributions,"
Mathematical Methods of Operations Research, 69:3 (July) pp 497-508.
[seems to be a pre-publication version here]

Best Answer

Related Solutions

Solved – Negative Binomial “Process”

Solved – Quantiles of a compound gamma/negative binomial distribution

Related Question