Combining risks II: zero-inflated distributions

Zeroes mess things up

Imagine I'm modeling risk: I have a set of low-probability risks, each of which has a distribution of possible losses, and each of which is independent. It's possible that none, some, or all of these risks happen. What I want is to know the overall risk distribution, in other words, what's the probability of losing $x? If each of these risks was certain to occur, it would be easy to combine them, but the presence of zeroes really messes things up as we'll see.

Let me make the problem clearer by laying out a simple case where we're combining the losses from two events:

• for an event \(F\), if it happens, there's a probability distribution for the range of losses \(f\), the probability of the event happening is 5% (\(P(F)\)),
• for an event \(G\), if it happens, there's a probability distribution for the range of losses \(g\), the probability of the event happening is 6% (\(P(G)\)),

We have four possibilities:

1. None of the events occur. 
2. Event \(F\) occurs and event \(G\) does not occur.
3. Event \(G\) occurs and event \(F\) does not occur.
4. Event \(F\) and \(G\) occur.

How do we model the overall loss distribution? Let's start with the simple part.

Combining probabilities

Working out the probabilities is standard and simple probability theory: 

1. None of the events occur: \(1 - P(F))(1-PG)))\)
2. Event F occurs and event G does not occur: \(P(F)(1-P(G))\)
3. Event G occurs and event F does not occur: \(P(G)(1-P(F))\)
4. Event F and G occur: \(P(F)P(G)\)

The sum of these four probabilities will be 1.

All of this doesn't tell us anything about the overall distribution of losses; how likely is it I will lose $x?

If there were no zeroes

If there were no zeros, combining the loss probabilities would be simple. In this case, we would have to assume the events did occur (so, \(P(F) = P(G) = 1\)). The overall loss distribution \(h\) is given by:

\[ h = f * g = \mathcal{F}^{-1}(\mathcal{F}(f)\mathcal{F}(g)) \]

in other words, it's the convolution of the distributions \(f\) and \(g\), where \(\mathcal{F}\) is the Fourier Transform.

Model zeroes with the dirac delta function

The problem comes down to how we represent zero losses as a continuous function so we can use the same math we use to handle continuous probability distributions. Fortunately, quantum physics has solved that problem for us with the dirac delta function. This isn't really a function, but I won't dig into that too much. 

Here's its definition.

\[ \delta(x) = \begin{cases} +\infty, & x = 0 \\ 0, & x \neq 0 \end{cases}\]

Here are some properties we're going to use later.

\[ \int_{-\infty}^{\infty} \delta(x)\,dx = 1 \]

\[ \int_{-\infty}^{\infty} f(x)\,\delta(x)\,dx = f(0) \]

\[ \int_{-\infty}^{\infty} f(x)\,\delta(x-a)\,dx = f(a) \]

\[ \delta(x) = \lim_{\epsilon \to 0^+} \frac{1}{\epsilon\sqrt{\pi}} e^{-x^2/\epsilon^2} \]

Zeroes and probability distributions

I started this blog with an example:

• for an event \(F\), if it happens, there's a probability distribution for the range of losses \(f\), the probability of the event happening is 5% (\(P(f)\)),

Here's how we can represent this using the dirac delta function:

\[ h(x) = P(F)\,f(x) + (1-P(F))\,\delta(x) \]

In words, this says we multiply the dirac delta function by the probability of it happening (\((1 - P(f))\)) and \(f\) by the probability of it happening (\(P(F)\)). We can see why this is so from the equation below.

\[ \int_{-\infty}^{\infty} h(x)\,dx = P(F)\int_{-\infty}^{\infty} f(x)\,dx + (1-P(F))\int_{\infty}^{\infty}\delta(x)\,dx = P(F) + (1-P(F)) = 1\]

As an aside if the dirac delta function is at \(a\) instead of zero, then:

\[ h(x) = P(F)\,f(x) + (1-P(F))\,\delta(x-a)\]

What this means is we have a way of representing the zero loss case and the loss case together. This is called a zero-inflated probability distribution, If we were to plot it as a probability distribution, we'd see a spike at zero and the probability distribution \(f\).

Zero inflated probability distributions and overall risk

Two risks

To go back to my starting example, we have two risk, \(h\) and \(j\) which are zero inflated probability distributions. We know that:

\[ h(x) = P(F)\,f(x) + (1-P(F))\,\delta(x) \]

\[ j(x) = P(G)\,g(x) + (1-P(G))\,\delta(x) \]  

To get the combined probability distribution we use convolution and Fourier Transforms:

\[ k = h * j = \mathcal{F}^{-1}(\mathcal{F}(h)\mathcal{F}(j)) \]

which is:

\[\begin{equation} \mathcal{F}\{h\}\,\mathcal{F}\{j\} = pq\,F G + p(1-q)\,F + (1-p)q\,G + (1-p)(1-q). \end{equation}\]

We can define these symbols:

\[F(\nu) = \mathcal{F}\{f\}(\nu)\]

and:

\[G(\nu) = \mathcal{F}\{g\}(\nu)]\]

Using these relationships:

\[\mathcal{F}^{-1}\{FG\} = f * g\]

\[\mathcal{F}^{-1}\{F\} = f\]

\[\mathcal{F}^{-1}\{G\} = g\]

\[\mathcal{F}^{-1}\{1\} = \delta\]

\[\mathcal{F}\{\delta\} = 1\]

We have:

\[\begin{equation} \begin{aligned} k(x) ={}& P(f)P(g)\,(f*g)(x) + P(f)\bigl(1-P(g)\bigr)\,f(x) \\ &+ \bigl(1-P(f)\bigr)P(g)\,g(x) + \bigl(1-P(f)\bigr)\bigl(1-P(g)\bigr)\,\delta(x). \end{aligned} \end{equation}\]

This equation tells us how to calculate the overall probability distribution, and so calculate the risk of losing \($x\) (or more formally, the risk of losses in the range \($x\) to \($y\)).

Notice there are contributions from zero events happening, \(F\) only happening, \(G\) alone happening, and \(F\) and \(G\) happening.

Multiple risks

The formula for two events is annoying, but you can work with it. What happens if there are multiple risks, defined like this?

\[h_i(x) = P(f_i)\,f_i(x) + \bigl(1-P(f_i)\bigr)\,\delta(x)\]

In this case, the overall probability distribution is given by:

\[ k = \mathcal{F}^{-1}\!\left(\prod_i \Bigl[ P(f_i)\,\mathcal{F}\{f_i\}(\nu) + \bigl(1-P(f_i)\bigr) \Bigr]\right) \]

That's much harder to calculate, especially with what I'm going to tell you next.

Fourier pairs and real-world difficulties

In some cases, the Fourier Transform of some distributions gives another distribution, or at least a nice closed formula. Sadly, most of the time, that isn't true. Here's what that means in practice; we start with some nice distributions, we zero inflate them, and then we have a calculation involving their Fourier Transform and the inverse of the combined formula, neither of which is a closed form. So what can we do? The answer is to use discrete approximations.

Using the discrete approximation

Using the discrete approximation involves digitizing the data, i.e., taking samples at regular discrete intervals. To do this, we have to use the same step size and the same starting point for all the equations. Once we've digitized the data, we can use the discrete Fourier Transform to get the Fourier Transform, do the products, then use the inverse discrete Fourier Transform to get to the overall loss distribution. This turns out to be difficult for all kinds of reasons. I'm just going to go over a few of the problems here. Bear in mind, we have finite computers with finite memory, and in most cases, we want the overall distribution calculation to be done quickly.

How do we represent distributions that go off to infinity? Many continuous distributions, including the normal, range from \(-\infty\) to \(\infty\). We can't digitize an infinite range, so we have to decide where to stop and start. A good rule of thumb is to stop when the value of the probability function is below some threshold number, e.g., \(10^{-8}\). What that threshold is depends on the calculations you're performing. Bear in mind that this truncation will change the overall calculations. You have to choose your threshold so that the effect is unimportant.

What grid size (step size) do we choose? If we choose too narrow a grid, computation time goes up a lot. If we choose too wide a grid, we may well miss peaks in the data.

Where do we start the grid from? If we start from the wrong place, we can miss edges of sharp functions. In practice, we start from zero to get the zero-inflated part of the distribution (which also involves another careful calculation).

The list of difficulties goes on. None of them are showstoppers, but what it does mean is that you need a level of expertise to perform these calculations correctly. This is not a place for beginners.

Central limit theorem to the rescue?

The central limit theorem (CLT) may well be helpful here. To put it simply, the CLT says that when you take samples from a distribution and sum them, the distribution of the sum follows a normal distribution. In most undergraduate texts, the CLT is formulated in a restrictive way; they say the samples have to be from the same distribution. That's not actually true; the samples have to be from distributions with certain properties.

In principle then, if we have a large enough number of risks, then the overall risk will follow a normal distribution, and we can figure out the properties of the normal distribution from the properties of the underlying  risk distributions. This would collapse the calculation time by orders of magnitude. To figure out when the CLT applies, we can use the Berry-Esseen theorem.

An alternative to all this

There is a way to get to the overall probability distribution without running into the closed form problem, and that's to use Monte Carlo analysis.

This problem is a particularly simple form of Monte Carlo. At its base, all you do is choose random numbers from the distributions and multiply them. In reality, things are a little more complicated, but still pretty straightforward. In practice, you might use quasi Monte-Carlo to increase the efficiency. 

The problem with Monte Carlo is accuracy. How many samples do you need to take to get an accurate representation of the overall probability distribution? It's likely over 32,000 for any moderately complicated problem, and maybe more. Quasi Monte Carlo gives you a lot of benefits, but ties you to powers of 2 for the sample size.

If an analytic solution is possible, the analytic solution will always be more accurate than Monte Carlo, and it may even be faster to calculate.

Discussion

I have to say, this whole process leaves me unsatisfied. 

I love the fact that we can write down an equation combining multiple risks, but I hate the fact that the digitization process is so awkward; if we could push digitization after the initial Fourier Transforms, that would be so much easier. Simplification comes down to finding the closed-form formula for the Fourier Transforms of distributions, but only a few distributions have closed form formula. 

I really dislike the arbitrary nature of selecting a grid size.

I like the fact that the CLT should kick in and simplify the process if you have enough risks, but I'm uncomfortable because of what this implies for computation time. It means that computation time increases with every risk, up to the point where the CLT kicks in, then computation time falls off a cliff. That feels a little arbitrary to me.

I dislike choosing Monte Carlo because of the effort needed to do it and the risk of an inaccurate result. It feels to me like if you have the equations for the underlying distributions, you really ought to be able to combine them analytically.

Conclusions

You can combine zero-inflated distributions to get overall distributions, but the process is complex and requires a strong knowledge of calculus and advanced probability theory.

AI statement

This blog was entirely human written with no help from AI. I formatted the formula with the help of AI.

Combining risks I: summing random variables

What is a random variable?

Let's start with two problems.

1. You're an insurance company that provides car and house insurance. You have a forecast of a range of claims for cars and a forecast of a range of claims for houses. How do you work out the overall range of claims' forecast?
2. You are a manufacturing company producing widgets in a multi-stage process. For each stage, you know the range of defect rates. How do you combine them to get an overall range of defect rates?

We can model each type of insurance or each manufacturing stage as a random variable. To simplify, a random variable is the result of a random process that gives values that follow a probability distribution. There's a much more precise definition in probability and statistics, but that's way beyond where I want to take this blog post. For now, what you need to know is that random variables are common, important, and have distinct mathematical properties.

As with probability distributions, there are two types of random variables: continuous and discrete. I'm going to focus on continuous random variables here. 

(Gemini)

Random variables in pictures

Let's imagine we're modeling car insurance claims. What's the next claim value going to be? We don't know exactly what it is (remember, the random part), but we do know it follows a probability distribution and very often, we know what the probability distribution is. Here's a normal probability distribution I'm going to use for illustrative purposes.

If our insurance claims follows this normal probability distribution, we don't know what the next claim will be, but we can assign it probabilities. For example, there's a 25% probability it will be between $3,026 and $4,371.

Some notation and definitions

A continuous random variable \(X\) is made up of a set of measurements \(\{x_0, x_1, \ldots\}\) that follow a probability distribution \(f(x)\) and a cumulative distribution function \(F\). Here are some things that must be true:

1. It must follow a probability distribution function (PDF) defined as \(P(a \le X \le b) = \int_a^b f(x), dx\). The PDF must have these properties:
   1. \(f(x) \ge 0\) for all \(x\)
   2. The total area "under the curve" for the PDF is 1, or: \(\int_{-\infty}^{\infty} f(x), dx = 1\)
2. There must be a cumulative distribution function (CDF) defined as \(F(x) = P(X \le x)\)

The expectation value is the mean value and is denoted by \(E(X)\). The variance is denoted by \(\text{Var}(X)\).

Summing random variables

Let's go back to the insurance company selling car and house insurance. Let's assume for now that the chart on the left is the probability distribution for the value of car claims in one year and the chart on the right is the value of house claims in a year. We'll also assume car and house claims are independent. How do we get the total claims distribution?

The problem is, any car claim value can combine with any house claim value. The depth of the problem's a bit easier to understand if we switch to the discrete domain. Let's say there are only three possible claim values for the house (\(f\)) and claims insurance (\(g\)) If every house insurance possibility can occur with any car insurance claim, then we can write the sum of the random variables as:

\[(f*g)[z] = \sum_{m=0}^{2} f[m]\,g[z-m]\]

That's fine in the discrete case, but there are an infinite number of values in the continuous case; we would have to multiply an infinite number of numbers by an infinite number of numbers, Or would we? It turns out, there's a short-cut that transforms convolution into multiplication.

Let's start with some basics. If we represent the sum as \(S = X_0 + X_1 + ...\) then:

\[ S_n = \sum_{i=1}^{n} X_i \]

by definition. By some math I won't go into, the mean (expectation value) is always the given by the sum of the expectation values:

\[E[S_n] = \sum_{i=1}^{n} E[X_i]\]

the formula for the sum of the variances is a little more complicated. If all the \(X\) variables are independent, then:

\[\text{Var}(S_n) = \sum_{i=1}^{n} \text{Var}(X_i)\]

if they aren't, then:

\[\text{Var}(S_n) = \sum_{i=1}^n \text{Var}(X_i) + 2\sum_{i<j} \text{Cov}(X_i, X_j)\]

So far, this is standard. Let's say each \(X\) has a minimum value of \(min_i\) and a maximum of \(max_i\) then:

\[\sum_{i=1}^{n} \min_i \ \le\ S_n \ \le\ \sum_{i=1}^{n} \max_i\]

which is a fancy way of saying the minimum is the sum of the minimum values and the maximum is the sum of the maximum values.

We also know that:

\[\int_{-\infty}^{\infty} f_{S_n}(z)\, dz = 1\]

The bottom line is, we have some information about the sum probability distribution, but we don't have the actual probability distribution of the sum. To get the probability distribution of the sum, we're going to turn to convolution.

Convolution to get the combined probability distribution

In probability theory, there are a couple of quantities defined for probability distributions. You can use these quantities to reason about the behavior of probability distributions and do math with them. There are two quantities that concern us: the moment generating function and the characteristic function. 

• The moment generating function encodes the mean, the variance, the skewness etc. (which are called moments). It's a convenient packaging of a distribution that you can use to get the moments and do math with distributions. The problem is, it doesn't always exist for all probability functions. 
• The characteristic function completely characterizes the probability distribution; you can get moments from it and do math with probability distributions, but it is harder to work with than the moment generating functions. By contrast with the moment generating function, it exists for all probability functions. 

We're going to be focused on the characteristic function here. This is how it's defined: \[\varphi_X(t) = E\left[e^{itX}\right] = \int_{-\infty}^{\infty} e^{itx} f(x)\, dx\]

unless you're deeply into math, you probably won't recognize the equation on the right, so let me tell you what it is; it's the Fourier Transform of the probability distribution. To complete the definitions, you can get \(f\) from \(\varphi\) via the inverse Fourier Transform:

\[f_X(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-itx}\, \varphi_X(t)\, dt\]

I came across Fourier Transforms in Physics. They were originally developed to model heat flow and have expanded into signal processing. They're all about decomposing a signal into a series of overlapping waves. So why have they popped up in probability theory? The answer is subtle, but it relates to this term \(e^{itX}\) in the definition of the characteristic function. We can expand it to \(e^{itX} = cos(tX) + i·sin(tX)\) which looks like waves - which in turn provides the links to Fourier analysis. This argument is a bit hand-wavy, but the it would take pages to get into the detail.

Let's turn to adding two random variables. To make it easier, I'll go back to the insurance example of a company offering car and house insurance. How do we add these two distributions to get overall losses?

The process of summing two random variables is called convolution. Here's how it's defined. The symbol for convolution is \(*\).

\[f_Z(z) = f_X(z) * f_Y(z) = (f_X * f_Y)(z) = \int_{-\infty}^{\infty} f_X(x)\, f_Y(z - x)\, dx\]

Once again, we're back in the world of Fourier Transforms. If we use the symbol \(\mathcal{F}\) to represent a Fourier Transform and \(\mathcal{F}^{-1}\) to represent an inverse Fourier Transform, then:

\[f_Z = f_X * f_Y = \mathcal{F}^{-1}\Big\{\mathcal{F}\{f_X\} \cdot \mathcal{F}\{f_Y\}\Big\}\]

What does going into the world of Fourier Transforms achieve? Two things:

1. It changes the calculation from convolution into multiplication of two functions; a much easier operation.
2. It means we have access to the huge amount of Fourier Transform literature and methods. When it comes to doing this on a computer, we can use the vast amount of work others have done to create Fourier Transform libraries. 

Of course, we can sum more than two random variables. Here are the relevant equations:

\[\varphi_{X_i}(t) = \int_{-\infty}^{\infty} e^{itx}\, f_{X_i}(x)\, dx\]

\[\varphi_{S_n}(t) = \prod_{i=1}^{n} \varphi_{X_i}(t)\]

\[f_{S_n}(z) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-itz}\, \varphi_{S_n}(t)\, dt\]

or more succinctly:

\[f_{S_n} = \mathcal{F}^{-1}\left\{ \prod_{i=1}^{n} \mathcal{F}\{f_{X_i}\} \right\}\]

Here's what you do in words and not symbols:

1. Get the probability distributions for each random variable.
2. Calculate the Fourier Transform for each random variable.
3. Multiply the Fourier Transforms.
4. Do an inverse Fourier Transform on the result.

The final step gives us the sum of the random variables.

Problems in the real world

This sounds easy, but it isn't. Although you can find the Fourier Transform for many probability distributions (search for characteristic function), multiplying them can be tricky. Even worse, transforming the product using an inverse Fourier Transform can be very, very hard. 

I'm going to take one example, let's imagine we want to sum two random variables, both with a beta distribution. The Fourier Transforms of a beta function looks like this:

\[\varphi_X(t) = {}_1F_1(\alpha;\ \alpha+\beta;\ it)\]

where:

\[{}_1F_1(a; b; z) = \sum_{k=0}^{\infty} \frac{a^{(k)}}{b^{(k)}} \cdot \frac{z^k}{k!}\]

Now, imagine multiplying two of these beasts and transforming them with an inverse Fourier Transform!

In fact, for many of these types of calculations, no closed-form solution exists.

How you solve this problem is a little complex, but first, let's take a detour and talk about some easier calculations.

Some short cuts exist

Normal distributions are easy. The convolution of a normal with another normal is a normal! The mean and variance of the summed normal distribution is:

\[\mu_{S_n} = E[S_n] = \sum_{i=1}^{n} \mu_i = \sum_{i=1}^{n} E[X_i]\]

\[\sigma_n = \sum_{i=1}^{n} \sigma_i^2\]

Uniform distributions are a bit harder. There are three cases:

1. All the distributions are exactly the same (i.e. iid or independent identically distributed) and are Uniform(0,1), then the sum follow an Irwin-Hall distribution.
   
   \[f_{S_n}(z) = \frac{1}{(n-1)!} \sum_{k=0}^{\lfloor z \rfloor} (-1)^k \binom{n}{k} (z-k)^{n-1}, \qquad 0 \le z \le n\]
   
   
2. All the random variables follow the same distribution, but it isn't Unform(0,1), instead is it's Uniform(a, b).
   
   \[ f_{S_n}(z) = \frac{1}{(b-a)(n-1)!} \sum_{k=0}^{\left\lfloor \frac{z-na}{b-a} \right\rfloor} (-1)^k \binom{n}{k} \left(\frac{z-na}{b-a} - k\right)^{n-1}, \qquad na \le z \le nb\]
   
   
3. The random variables follow different Uniform distributions. In this case, there is no convenient closed-form expression for the sum. 

Here are some of the other shortcuts.

\begin{array}{|l|l|l|} \hline \textbf{Distribution} & \textbf{Condition} & \textbf{Sum is} \\ \hline \text{Normal} & X_i \sim N(\mu_i, \sigma_i^2) & N\left(\sum \mu_i,\ \sum \sigma_i^2\right) \\ \hline \text{Poisson} & X_i \sim \text{Pois}(\lambda_i) & \text{Pois}\left(\sum \lambda_i\right) \\ \hline \text{Binomial} & X_i \sim \text{Bin}(n_i, p),\ \text{same } p & \text{Bin}\left(\sum n_i,\ p\right) \\ \hline \text{Exponential (same rate)} & X_i \sim \text{Exp}(\lambda),\ \text{i.i.d.} & \text{Gamma}(n, \lambda) \\ \hline \text{Gamma (same rate)} & X_i \sim \text{Gamma}(k_i, \lambda) & \text{Gamma}\left(\sum k_i,\ \lambda\right) \\ \hline \text{Chi-squared} & X_i \sim \chi^2_{k_i} & \chi^2_{\sum k_i} \\ \hline \text{Negative Binomial (same } p\text{)} & X_i \sim \text{NB}(r_i, p) & \text{NB}\left(\sum r_i,\ p\right) \\ \hline \text{Cauchy} & X_i \sim \text{Cauchy}(x_{0,i}, \gamma_i) & \text{Cauchy}\left(\sum x_{0,i},\ \sum \gamma_i\right) \\ \hline \end{array}

What happens if you're trying to sum random variables with different distributions or your sum isn't one of the ones listed here? You have two choices: slog through the math yourself or step into the digital world.

Summing random variables in practice

Often, there's no closed form expression for the sum, so you have to do thing a little differently. Let's assume you're trying to sum several random variables with very different probability distributions. What can you do?

The first step is to go as far as you can with the standard formula I gave above. For example, if you have multiple normally distributed random variables, sum them. 

In many cases, you can figure out the characteristic function for the distributions; these are the Fourier Transforms. You can now write out the Fourier Transform of your sum as a series of characteristic function multiplications. Now, evaluate that expression on a grid. For example, if every characteristic function is a function of x, choose x values from a minimum to a maximum with a step size. You may need to use hundreds of thousands or even millions of x values. Finally, use a digital inverse FFT (Fast Fourier Transform) to transform back and get the distribution of the sum.

This process sometimes needs a great amount of computing power. If the power you need is greater than what's available there is another approach. This involves Monte Carlo simulation, or better still, Quasi Monte Carlo simulation. I'm not going to go into that here, it's way too involved.

The overall process of summing random variables is obviously more involved than what I've described, but you get the gist.

Zero-inflation is a problem

So far, I've only handled the case where the probabilities are 1, i.e., there will be car insurance claims or house claims. What happens if that isn't the case? What happens if there's a probability of no losses whatsoever? This is a more complicated situation I'll deal with in another blog post.

Summary

To sum two or more random variables, we use convolution, and Fourier Transforms are very much part of that process. In some cases, there are nice closed-form expressions for the sum, but most of the time, there are not. In practice, we need to digitize the data and step into the world of the discrete Fourier Transform.

AI statement

All the writing was don by me (a human). I used AI to generate the formulae and the charts.