Negative Binomial in VBA

 Have you ever tried to simulate a negative binomial random variable in a Spreadsheet?

If the answer to that is ‘nope – I’d just use Igloo/Metarisk/Remetrica’ then consider yourself lucky! Unfortunately not every actuary has access to a decent software package, and for those muddling through in Excel, this is not a particularly easy task.

If on the other hand your answer is ‘nope – I’d use Python/R, welcome to the 21st century’. I’d say great, I like using those programs as well, but sometimes for reasons out of your control, things just have to be done in Excel.
​
This is the situation I found myself in recently, and here is my attempt to solve it:

Attempt 0

The first step I took in attempting to solve the problem was of course to Google it, then cross my fingers and hope that someone else has already solved it and this is just going to be a simple copy and paste. Unfortunately when I did search for VBA code to generate a negative binomial random variable, nothing comes up. In fact, nothing comes up when searching for code to simulate a Poisson random variable in VBA. Hopefully if you've found your way here, looking for this exact thing, then you're in luck, just scroll to the bottom and copy and paste my code.

When I Googled it, there were a few solutions that almost solved the problem; there is a really useful Excel add-in called ‘Real statistics’ which I’ve used a few times:
http://www.real-statistics.com/

It's a free excel add-in, and it does have functionality to simulate negative bimonials. If however you need someone to be able to re-run the Spreadsheet, they also will need to have it installed. In that case, you might as well use Python, and then hard code the frequency numbers. Also, I have had issues with it slowing Excel down considerably, so I decided not to use this in this case.

I realised I’d have to come up with something myself, which ideally would meet the following criteria

• Is free (aren’t all the best things)
• Can be bundled into a Spreadsheet
• Runs in a reasonable amount of time

How hard can that be?

Attempt 1

I’d seen a trick before (from Pietro Parodi’s excellent book ‘Pricing in General Insurance’) that a negative binomial can be thought of as a Poisson distribution with a Gamma distribution as the conjugate prior. See the link below for more details:

https://en.wikipedia.org/wiki/Conjugate_prior#Table_of_conjugate_distributions

Since Excel has a built in Gamma inverse, we have simplified the problem to needing to write our own Poisson inverse. We can then easily generate negative binomials using a two step process:

1. Generate the Gamma sample which we will use as our lambda
2. Generate a Poisson sample, where the Poisson dist has mean equal to our gamma sample

Great, so we’ve reduced our problem to just being able to simulate a Poisson in VBA. Unfortunately there’s still no built in Poisson inverse in Excel (or at least the version I have), so we now need a VBA based method to generate this.
​
There is another trick we can use for this - which is also taken from Pietro Parodi - the waiting time for a Poisson dist is an Exponential Dist. And the CDF of an Exponential dist is simple enough that we can just invert it and come up with a formula for generating an Exponential sample. We then set up a loop and add together exponential values, to arrive at Poisson sample. The code for this is give below:

Function Poisson_Inv(Lambda As Double)
    
    s = 0
    N = 0
    Do While s < 1
        u = Rnd()
        s = s - (Application.WorksheetFunction.Ln(u) / Lambda)
        k = k + 1
    Loop
    BH_Poisson_Inv = (k - 1)

End Function

​The VBA code for our negative binomial is therefore:

Function NegBinOld2(b, r)

    Dim Lambda As Double
    Dim N As Long
    
    u = Rnd()

    Lambda = Application.WorksheetFunction.Gamma_Inv(u, r, b)

    N = Poisson_Inv(Lambda)

    NegBinOld2 = N

End Function

Does this do everything we want?

• It’s free (yay!)
• It can be included as a VBA module of a Spreadsheet
• It does run reasonably quickly

 
There are a couple of downside of though:

• Custom VBA functions referenced directly in the Spreadsheet can be a bit buggy when refreshing - it would be good to have more control over this
• Every time we refresh we will get a different result, and there is no way to set a seed to prevent this.

 
This leads us on to Attempt 2

Attempt 2

If we pass the VBA a random uniform sample, then whenever we hit refresh in the Spreadsheet the random sample will refresh, which will force the Negative Binomial to resample. Without this, sometimes the VBA will function will not reload. i.e. we can use the sample to force a refresh whenever we like. Adapting the code gives the following:
​
​

Function NegBinOld(b, r, Rnd1 As Double)

    Dim Lambda As Double
    Dim N As Long
    
    u = Rnd1

    Lambda = Application.WorksheetFunction.Gamma_Inv(u, r, b)

    N = Poisson_Inv(Lambda)

    NegBinOld = N

End Function

So this solves the refresh problem. What about the random seed problem? Even though we now always get the same lambda for a given rand – and personally I quite like to hardcode these in the Spreadsheet once I’m happy with the model, just to speed things up. We still use the VBA rand function to generate the Poisson, this means everytime we refresh, even when passing it the same rand, we will get a different answer and this answer will be non-replicable. This suggests we should somehow use the first random uniform sample to generate all the others in a deterministic (but still pseudo-random) way.

Attempt 3

The way I implemented this was to the set the seed in VBA to be equal to the uniform random we are passing the function, and then using the VBA random number generator (which works deterministically for a given seed) after that. This gives the following code:
​

Function NegBin(b, r, Rnd1 As Double)
    
    Rnd (-1)
    Randomize (Rnd1)
    
    Dim Lambda As Double
    Dim N As Long
       
    u = Rnd()

    Lambda = Application.WorksheetFunction.Gamma_Inv(u, r, b)

    N = Poisson_Inv(Lambda)

    NegBin = N

End Function

So we seem to have everything we want – a free, quick, solution that can be bundled in a Spreadsheet, which allows other people to rerun without installing any software, and we’ve also eliminated the forced refresh issue. What more could we want?

The only slight issue with the last version of the negative binomial is that our parameters are still specified in terms of ‘b’ and ‘r’. Now what exactly are ‘b’ and ‘r’ and how do we relate them to our sample data? I’m not quite sure.... The next trick is shamelessly taken from a conversation I had with Guy Carp’s chief Actuary about their implementation of severity distributions in MetaRisk.

Attempt 4

​Why can't we reparameterise the distribution using parameters that we find useful, instead of feeling bound by using the standard statistics textbook definition (or even more specifically the list given in the appendix to ‘Loss Models – from data to decisions’, which seems to be somewhat of an industry standard), why can't we redefine all the parameters from all common actuarial distributions using a systematic approach for parameters?

Let's imagine a framework where no matter which specific severity distribution you are looking at, the first parameter contains information about the mean (even better if it is literally scaled to the mean in some way), the second contains information about the shape or volatility, the third contains information about the tail weight, and so on. This makes fitting distributions easier, it makes comparing the goodness of fit of different distributions easier, and it make sense checking our fit much easier.

I took this idea, and tied this in neatly to a method of moments parameterisation, whereby the first value is simply the mean of the distribution, and the second is the variance over the mean. This gives us our final version:
​

Function NegBin(Mean, VarOMean, Rnd1 As Double)
    
    Rnd (-1)
    Randomize (Rnd1)
    
    Dim Lambda As Double
    Dim N As Long
    
    b = VarOMean - 1
    r = Mean / b
    
    u = Rnd()

    Lambda = Application.WorksheetFunction.Gamma_Inv(u, r, b)

    N = Poisson_Inv(Lambda)

    NegBin = N

End Function

Function Poisson_Inv(Lambda As Double)
    
    s = 0
    N = 0
    Do While s < 1
        u = Rnd()
        s = s - (Application.WorksheetFunction.Ln(u) / Lambda)
        k = k + 1
    Loop
    BH_Poisson_Inv = (k - 1)

End Function