5.7. The Pollaczeck-Khinchine Formula
Objectives: Illustrate and compute Pollaczeck-Khinchine formula for probability of eventual ruin for a compound Poisson process.
Audience: Advanced users.
Prerequisites: Advanced risk theory.
Contents:
5.7.1. Helpful References
5.7.2. Classical Risk Theory and the Pollaczeck-Khinchine Formula
The Pollaczeck-Khinchine formula determines the probability of eventual ruin in a portfolio where claims are driven by a compound Poisson process, in terms of starting surplus and the premium rate. Losses are generated by a Poisson process with \(\lambda\) annual expected claims and iid severity \(X\). Losses up to time \(t\) are given by
where \(N(t)\) is Poisson with mean \(\lambda t\). Expected loss per year equals \(\lambda\mathsf{E}[X]\). Premium per year equals \((1+r)\lambda \mathsf{E}[X]\) where \(r\) is the ratio of profit to expected loss. The corresponding expected loss ratio is \(1/(1+r)\). If \(r\le 0\) then eventual ruin is certain, so assume \(r>0\).
Define the integrated severity distribution by
where \(S\) is the survival function of \(X\). \(F_I\) is a thicker tailed distribution than \(F\). Let
denote accumulated surplus to time \(t\) given starting surplus \(u\). \(U\) is called the surplus process. Finally, let
be the probability of eventual ruin.
The Pollaczeck-Khinchine formula says that
where \(F_I^{n*}\) is distribution of the sum of \(n\) independent variables with distribution \(F_I\). Note that
is the probability generating function of a geometric distribution \(M\) with mean \(1/r\) and \(\Pr(M=m)=\frac{r}{1+r}\frac{1}{(1+r)^m}\). Therefore \(\psi(u)=\Pr(Y > u)\) where \(Y\) is an aggregate distribution with frequency \(M\) and severity \(F_I\). A surprising consequence is that the probability of eventual ruin starting with no surplus, \(\psi(0)=1-\Pr(Y=0)=1-\Pr(M=0)=\frac{1}{1+r}\) equals the expected loss ratio!
Embrechts et al. [1997] Section 1.2 shows how to derive the Pollaczeck-Khinchine formula. The key step is to determine the distribution of \(X-(1+r)T\) where \(T\) is the exponential waiting time between claims, and to observe that ruin can occur only at the moment of a claim.
The Pollaczeck-Khinchine formula gives combinations of \(u\) and \(r\) that are consistent with a top-down stability requirement expressed as a target probability of eventual ruin. Overlaying a cost of capital provides a link between \(r\) and \(u\) that determines a minimum viable market size. An example of this method is given below.
Because eventual is the same in days, weeks or years, \(\psi_{X, m}(u)\) is independent of the expected claim count \(\lambda\). In unit of time \(1/\lambda\) all portfolios have an expected claim count of one. Therefore \(\psi^{-1}(p)\) gives a capital requirement (risk measure) that is a function of severity and not frequency, i.e., it is independent of portfolio size. Unlike most risk measures, it does not regard small portfolios as more risky than large ones.
The Cramer-Lundberg formula is an approximation to \(\psi\) that applies for thin tailed severities. It says that
where \(k>0\) is a constant called the adjustment coefficient solving
where \(P=(1+r)\lambda\mathsf{E}[X]\) is the premium. Given a top-down stability requirement, we can work backwards from the Cramer-Lundberg formula to determine a premium.
Exercise. Show that if \(k=-\log(p)/u\) and premium
then the Cramer-Lundberg formula ensures the probability of eventual ruin is \(\le p\). The properties of \(P\) motivate the exponential premium. In turn, the approximation \(P\approx \mathsf{E}[A(1)] + k\mathsf{Var}(A(1))/2\) motivates the variance principle.
Both the Cramer-Lundberg and Pollaczeck-Khinchine formulas assume independent and identically distributed severity and Poisson frequency. These can be reasonable assumptions for the loss process of a small portfolio. The case of a mixed Poisson can be decomposed as a mixture of pure Poisson processes.
5.7.3. FFT Computation
The distribution of \(Y\) can be computed using Fast Fourier transforms in
the same way as any aggregate distribution. Some care is needed when the
margin is very small because the claim count is very large.
Aggregate includes pollaczeck_khinchine() to determine the
integrated distribution \(F_I\) and convolve it with a geometric
frequency.
5.7.4. Using The Pollaczeck-Khinchine Formula I
This section reproduces two examples from the risk vignette for the actuar package.
The first is based on a mean 10 Poisson compound with shape 2 gamma severity. The vignette uses matched-moments discretization and so our numbers do not exactly match, but they are very close. We build the compound, compute some quantiles and tvars, display the density (compare p.8-10). Then using a premium loading of 20% using the expected value premium (p.12), we reproduce the probabilities in Figure 5. Our computation is exact vs. using an approximation.
In [1]: from aggregate import build, qd
In [2]: import matplotlib.pyplot as plt
In [3]: a = build('agg Actuar 10 claims sev gamma 2 poisson'
...: , bs=0.5, log2=8)
...:
In [4]: qd(a)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 10 0.31623 0.31623
Sev 2 1.9999 -7.4963e-05 0.70711 0.71102 1.4142 1.3901
Agg 20 19.999 -7.4963e-05 0.3873 0.38801 0.5164 0.51634
log2 = 8, bandwidth = 1/2, validation: fails sev cv, agg cv.
In [5]: ps = [0.25, .5, .75, .9, .95, .975, .99, .995, .999, 1-1e-14]
In [6]: qd(pd.DataFrame({'p': ps, 'q': a.q(ps)}), float_format=lambda x: f'{x:10.3f}' if x < 1 else f'{x:10.1f}')
p q
0 0.250 14.5
1 0.500 19.5
2 0.750 25.0
3 0.900 30.5
4 0.950 34.0
5 0.975 37.0
6 0.990 41.0
7 0.995 43.5
8 0.999 49.5
9 1.000 115.5
In [7]: fig, axs = plt.subplots(1, 3, figsize=(3 * 3.5, 2.45), constrained_layout=True)
In [8]: ax, ax0, ax1 = axs.flat
In [9]: a.density_df.F.plot(ax=ax, xlim=[0, 60], title='Po-Gamma distribution function');
In [10]: qd(a.density_df.p_total.head(20).reset_index(drop=False), float_format=lambda x: f' {x:<12.5g}' if 0 < x < .1 else f'{x:10.1f}')
loss p_total
0 0.0 5.9175e-05
1 0.5 8.6904e-05
2 1.0 0.00017152
3 1.5 0.00028809
4 2.0 0.00045063
5 2.5 0.00066863
6 3.0 0.0009508
7 3.5 0.0013052
8 4.0 0.001739
9 4.5 0.0022578
10 5.0 0.0028658
11 5.5 0.0035653
12 6.0 0.0043564
13 6.5 0.0052372
14 7.0 0.0062034
15 7.5 0.0072487
16 8.0 0.0083649
17 8.5 0.0095419
18 9.0 0.010768
19 9.5 0.01203
In [11]: qd(a.tvar([.9, .95, .99]))
[35.02023434 38.14387968 44.6199146 ]
In [12]: ruins, find_us, mean, dfi = a.cramer_lundberg(.2)
In [13]: ax0.plot(np.cumsum(dfi), label='integrated')
Out[13]: [<matplotlib.lines.Line2D at 0x7f8089819ed0>]
In [14]: ax0.plot(a.density_df.p_sev.cumsum(), label='severity')
Out[14]: [<matplotlib.lines.Line2D at 0x7f8089819b10>]
In [15]: ax0.set(xlim=[0, 40], title='Severity and integrated severity distributions')
Out[15]: [(0.0, 40.0), Text(0.5, 1.0, 'Severity and integrated severity distributions')]
In [16]: ax0.legend(loc='lower right')
Out[16]: <matplotlib.legend.Legend at 0x7f8087ca3a90>
In [17]: ruins.plot(ax=ax1, xlim=[0, 50],
....: title='Probability of eventual ruin against starting surplus',
....: ylabel='Probability', xlabel='Starting surplus');
....:
The second uses a Pareto severity, where the integrated distribution can be computed exactly. The following code produces the exact values for the probability of eventual default against starting surplus. All the values fall in the range between lower and upper shown on p.19.
In [18]: a = build('agg Actuar2 1 claim sev 4 * pareto 5 - 4 fixed')
In [19]: qd(a)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
X
Freq 1 0
Sev 1 1 -4.2872e-06 1.291 1.2907 4.6476 4.5872
Agg 1 1 -4.2872e-06 1.291 1.2907 4.6476 4.5872
log2 = 16, bandwidth = 1/512, validation: fails sev skew, agg skew.
In [20]: ruins, find_us, mean, dfi = a.cramer_lundberg(.2)
In [21]: ruins.name = 'Prob'
In [22]: bit = ruins.loc[np.arange(0, 51, 5)].to_frame()
In [23]: bit.index = bit.index.astype(int)
In [24]: bit.index.name = 'Initial surplus'
In [25]: qd(bit, ff=lambda x: f'{x:.5f}')
Prob
Initial surplus
0 0.83306
5 0.42670
10 0.23477
15 0.13149
20 0.07439
25 0.04242
30 0.02435
35 0.01406
40 0.00818
45 0.00479
50 0.00282
The actual work, to get the answer as opposed to formatting the result, is only two lines of code in Aggregate (the first and third) vs. 8 in actuar R.
5.7.5. Using The Pollaczeck-Khinchine Formula II
This section illustrates the theory using a lognormal severity with a mean of 50,000 and a CV of 10 (\(\sigma=2.15\)) corresponding to a moderately risky liability line. It compares starting surplus levels for different eventual ruin probabilities assuming a margin \(r=0.1\) with a 1 million and 10 million occurrence limit. It also illustrates simulations of the surplus process in each case with starting surplus calibrated to a 0.05 probability of eventual ruin.
Set up the portfolio.
In [26]: port = build('port PZTest '
....: 'agg Limit1 '
....: '0.1 claims '
....: '1000000 xs 0 '
....: 'sev lognorm 50000 cv 10 '
....: 'poisson'
....: 'agg Limit10 '
....: '0.1 claims '
....: ' 10000000 xs 0 '
....: 'sev lognorm 50000 cv 10 '
....: 'poisson'
....: , bs=500, log2=18, padding=1)
....:
In [27]: qd(port)
E[X] Est E[X] Err E[X] CV(X) Est CV(X) Skew(X) Est Skew(X)
unit X
Limit1 Freq 0.1 3.1623 3.1623
Sev 38064 38060 -0.0001105 3.0931 3.0935 5.9141 5.914
Agg 3806.4 3806 -0.0001105 10.28 10.281 18.845 18.845
Limit10 Freq 0.1 3.1623 3.1623
Sev 47900 47896 -8.7807e-05 5.6958 5.6963 20.969 20.969
Agg 4790 4789.6 -8.7807e-05 18.287 18.289 64.966 64.966
total Freq 0.2 2.2361 2.2361
Sev 42982 42978 -9.7853e-05 4.8898 23.502
Agg 8596.4 8595.6 -9.7853e-05 11.16 11.161 50.728 50.728
log2 = 18, bandwidth = 500, validation: fails sev mean, agg mean.
The left plots show the Pollaczeck-Khinchine formula starting surplus as a function of the eventual ruin probability with margin 0.1 on linear (solid) and log (dashed) scales. The Cramer-Lundberg formula says that the probability of eventual ruin is approximately exponential, which is a straight line on a log scale.
The right column show 500 simulated surplus paths, with \(\times\) indicating ruin scenarios. Capital calibrated to 0.05 eventual ruin probability. Time and volume are symmetric in the model, so volume can be regarded as time for a fixed size portfolio or a varying sized portfolio for a fixed time or a combination. Scale indicates cumulative exposure-years.
In [28]: from aggregate.extensions.pir_figures import fig_9_1
In [29]: fig_9_1(port)
The right hand plots are computed here with only 100 samples, vs. 500 used in the book, and so the approximation is not as accurate.
These simulations show that the probability of eventual ruin is constrained by the buildup of surplus in most scenarios. Defaults occur early in the simulated history. This model could be appropriate for a mutual company—indeed some mutual companies have accumulated substantial amounts of capital. For a stock company, a more realistic approach adds dividends to manage capital.
5.7.6. Market Scale and Viability
Given severity \(X\) and ratio \(r\) of margin to expected loss, the Pollaczeck-Khinchine function \(\psi\) is monotone and hence invertible, allowing us to find \(u_{X,r}(p)=\psi_{X,r}^{-1}(p)\), the starting capital necessary to guarantee probability \(p\) of eventual ruin.
The amount of margin equals \(r\lambda\mathsf{E}[X]\), where \(\lambda\) is the annual expected claim count. Since the expected margin must pay the cost of capital, we get a market viability constraint
where \(\iota\) is the cost of capital. Each element is influenced by different factors:
the hazard and contract design determines \(X\),
the insurance product market determines \(r\),
the capital markets determine \(\iota\), and
a regulator or rating agency determines (or strongly influences) \(p\).
There are two ways to apply this formula.
First, consider a diversifying unit, such as motor liability, where insurers grow by adding new, independent insureds with the same severity. Here, the formula gives a minimum size of market constraint
This function of four variables and \(\lambda\) is:
Increasing and linear in \(\iota\): the market must be larger given more expensive capital.
Decreasing in \(r\): the market can be smaller with a higher margin.
Decreasing in \(p\): the market must be larger to support a stricter capital standard.
Independent of expected severity (because \(\psi\) is homogeneous in \(\mathsf{E}[X]\)) but dependent on the shape of severity (which influences \(\psi\)).
It shows the natural scale using the lognormal example. Size, measured by expected annual claim count, is shown for a range of margins, limits, and stability constraints. If claim frequency is 5%, the table shows that a market with 1M limits is reasonable for all \(p\) and \(r\). For example, the strictest stability constraint \(p=0.01\) and lowest margin rate \(r=0.025\) needs 39,814 claims, or about 800,000 policies, to be viable. With a 10M limit and same \(r\), the market size needs to be >100,000 claims, or about 2.5 million policies, which is less achievable. However, if the margin rate increases to \(r=0.1\), the market size reduces to 9,088 claims or about 180,000 policies.
In [30]: from aggregate.extensions.pir_figures import natural_scale
In [31]: df = natural_scale(port)
index 262144 is out of bounds for axis 0 with size 262144
In [32]: qd(df, float_format=lambda x: f'{x:,.0f}')
Margin 25.000m 50.000m 75.000m 100.000m 150.000m 200.000m 250.000m
Limit p
1.000M 10.000m 39,814 10,182 4,625 2,657 1,230 718 477
50.000m 25,808 6,577 2,978 1,705 784 455 300
100.000m 19,777 5,025 2,269 1,295 593 342 225
250.000m 11,803 2,973 1,331 754 339 193 124
10.000M 10.000m NaN 33,612 15,552 9,088 4,335 2,600 1,765
50.000m NaN 21,571 9,916 5,758 2,714 1,610 1,079
100.000m NaN 16,386 7,489 4,324 2,015 1,182 794
250.000m NaN 9,531 4,280 2,427 1,101 614 382
Note: these numbers differ slightly from the book because of the update parameters used for port.
Second, consider a non-diversifying unit, writing catastrophe exposed business, where insurers grow by covering a greater proportion of each event. Severity becomes market share times an industry severity \(X\), and the number of events is fixed. US hurricane reinsurance is an example. In this case, viability is independent of market share and is controlled by whether the inequality has a solution that is acceptable to both the product market and the capital market. Viability is harder to achieve
with smaller \(\lambda\): rare events are more difficult to insure,
with lower \(p\): higher quality insurance is more expensive,
with higher \(\iota\): more costly capital, and
with lower \(r\) because \(u\) increases quickly.