Efficient Implementation of Panjer Recursion

From a Recursion to an Algorithm

Panjer recursion turns a compound aggregate-loss calculation into a sequence of updates: once the probability of zero aggregate loss is known, calculate the probability at grid point 1, then grid point 2, and continue upward. An implementation must keep the frequency model, the severity probabilities, and the grid indexing consistent. A wrong initial value or an update in the wrong order can invalidate every later entry.

Suppose \(N\) is the claim count and claim severities are independent and identically distributed, independent of \(N\), and take nonnegative integer values on a chosen loss grid. Write \(f_j=P(X=j)\), where \(X\) is the severity in grid units, and write \(g_k=P(S=k)\) for \(S\), the sum of the severities of all claims. The aim is to calculate \(g_0,\ldots,g_K\), for a specified maximum grid index \(K\).

Definition: A frequency distribution is in the Panjer \((a,b,0)\) class if its probabilities \(p_n=P(N=n)\) satisfy \(p_n=(a+b/n)p_{n-1}\) for \(n\geq 1\). It is in the \((a,b,1)\) class if this relation holds for \(n\geq 2\), with \(p_0\) and \(p_1\) specified separately. For either class, \(g_0=E[f_0^N]\), or equivalently \(g_0=G_N(f_0)\), where \(G_N\) is the probability generating function of \(N\).

For an \((a,b,0)\) frequency, the update for \(k\geq 1\) is

g_k=\frac{1}{1-af_0}\sum_{j=1}^{k}\left(a+\frac{bj}{k}\right)f_jg_{k-j}.

Only severity indices \(j\leq k\) contribute to \(g_k\): a claim with severity greater than \(k\) cannot be part of an aggregate total of \(k\). If the severity has finite support \(0,\ldots,J\), the sum therefore ends at \(\min(k,J)\). For an \((a,b,1)\) frequency, define \(c=p_1-(a+b)p_0\). The corresponding update has the additional term \(cf_k\) inside the numerator:

g_k=\frac{\displaystyle\sum_{j=1}^{k}\left(a+\frac{bj}{k}\right)f_jg_{k-j}+cf_k}{1-af_0}.

The \((a,b,1)\) formula applies when the stated frequency probabilities obey the defining relation for \(n\geq2\). Do not silently use the \((a,b,0)\) update for that class. In either case, the recursion requires \(1-af_0\ne 0\); if this denominator vanishes, this update cannot be used as written.

Why the Update Uses Earlier Entries

A short coefficient argument explains both the order of calculation and the factors in the formula. For an \((a,b,0)\) frequency, let \(F(z)=\sum_{j\geq0}f_jz^j\) and \(G(z)=\sum_{n\geq0}p_nz^n\) be the severity and frequency probability generating functions. The frequency relation implies

$$ (1-aF(z))G'(F(z))=(a+b)G(F(z))F'(z). $$

The aggregate generating function is \(H(z)=G(F(z))=\sum_{k\geq0}g_kz^k\). Equating the coefficients of \(z^{k-1}\) on both sides gives

k(1-af_0)g_k=\sum_{j=1}^{k}\bigl(ak+bj\bigr)f_jg_{k-j}.

Dividing by \(k(1-af_0)\) gives the stated update. Every \(g_{k-j}\) on its right-hand side has index below \(k\), so the calculation proceeds in ascending order. This is the reason an in-place array can be filled from left to right without requiring future entries.

An Ordered Implementation

For a target grid ending at \(K\), first determine the frequency parameters and the severity probabilities. Then calculate \(g_0\), and only then update the remaining entries in increasing order. For an \((a,b,1)\) model, also supply \(p_0\) and \(p_1\) to calculate its correction term. The following process makes the dependencies explicit.

Fix the grid and inputs.
Choose \(K\). Supply \(f_0,\ldots,f_K\) (or all severity probabilities through the finite support \(J\), if \(J<K\)), as well as the appropriate frequency parameters. The severity probabilities must refer to the same grid units as the aggregate indices.

Initialize the aggregate array.
Set \(g_0=G_N(f_0)=\sum_{n\geq0}p_nf_0^n\). This accounts for positive claim counts whose claims all have zero severity, as well as the possibility of no claims.

Loop over aggregate indices.
For \(k=1,\ldots,K\), sum over \(j=1,\ldots,\min(k,J)\) when severity support ends at \(J\). Use the update for the correct frequency class and divide by \(1-af_0\).

Retain and interpret the requested output.
Store each \(g_k\) at its matching index. The entries give probabilities at grid totals; they do not by themselves describe totals beyond \(K\).

In pseudocode, the core loop for an \((a,b,0)\) frequency is: set \(g[0]\) to \(G_N(f[0])\); for each \(k\) from 1 through \(K\), set a running sum to zero; for each \(j\) from 1 through \(\min(k,J)\), add \((a+bj/k)f[j]g[k-j]\); then set \(g[k]\) to the running sum divided by \(1-af[0]\). The inner loop reads only entries calculated earlier. For an \((a,b,1)\) frequency, add \(cf[k]\) to the running sum before division.

Worked Example: A Binomial Frequency With Zero-Severity Claims

A hypothetical portfolio has \(N\sim\operatorname{Binomial}(2,0.5)\) claims. Each claim has grid severity probabilities \(f_0=0.2,\ f_1=0.5,\ f_2=0.3\), and no other possible severity. Calculate \(g_0,g_1,g_2,g_3\). Use the \((a,b,0)\) recursion and show the ordered updates.

Method and inputs. For a binomial count with two trials and claim probability \(0.5\), the frequency probabilities are \(p_0=0.25,\ p_1=0.5,\ p_2=0.25\). Their ratios give \(a=-1,\ b=3\): \(p_1/p_0=2=a+b\), and \(p_2/p_1=0.5=a+b/2\). The recursion denominator is \(1-af_0=1-(-1)(0.2)=1.2\).

Initialize. Since \(g_0=G_N(f_0)\),

$$ g_0=(0.25+0.5(0.2)+0.25(0.2)^2)=0.25+0.10+0.01=0.36. $$

Calculate \(g_1\). Only \(j=1\) contributes:

g_1=\frac{(a+b)f_1g_0}{1-af_0} =\frac{2(0.5)(0.36)}{1.2} =\frac{0.36}{1.2}=0.30.

Calculate \(g_2\). Both \(j=1\) and \(j=2\) contribute:

g_2=\frac{\left(\frac{a(2)+b(1)}{2}\right)f_1g_1+ \left(\frac{a(2)+b(2)}{2}\right)f_2g_0}{1.2} =\frac{(0.5)(0.5)(0.30)+(2)(0.3)(0.36)}{1.2} =\frac{0.075+0.216}{1.2}=0.2425.

Calculate \(g_3\). The \(j=3\) term is zero because \(f_3=0\):

g_3=\frac{(0)(f_1g_2)+(1)(f_2g_1)+(2)(f_3g_0)}{1.2} =\frac{0+(1)(0.3)(0.30)+0}{1.2} =\frac{0.09}{1.2}=0.075.

The result is \((g_0,g_1,g_2,g_3)=(0.36,0.30,0.2425,0.075)\). As an arithmetic check, the aggregate generating function is \((0.5+0.5F(z))^2=(0.6+0.25z+0.15z^2)^2\); its coefficients through \(z^3\) are the same. The common mistake avoided here is omitting the factor \(1/k\) from the coefficient \(ak+bj\) when using the form of the recursion with \(a+bj/k\).

Worked Example: Geometric Frequency and the Correct Parameters

A hypothetical count has \(P(N=n)=0.6(0.4)^n\), for \(n=0,1,\ldots\). Each claim has severity 0 or 1 grid unit, each with probability \(0.5\). Calculate \(g_0,g_1,g_2\) and identify the frequency parameters.

Method and inputs. The ratio \(p_n/p_{n-1}=0.4\) for every \(n\geq1\), so this is an \((a,b,0)\) frequency with \(a=0.4,\ b=0\). In particular, it is not represented by \(a=0,\ b=0.4\). The severity inputs are \(f_0=f_1=0.5\), with \(f_j=0\) for \(j\geq2\), and the denominator is \(1-af_0=1-(0.4)(0.5)=0.8\).

Initialize. The frequency generating function is \(G_N(s)=0.6/(1-0.4s)\), so

g_0=G_N(0.5)=\frac{0.6}{1-0.4(0.5)} =\frac{0.6}{0.8}=0.75.

Update \(g_1\). Only \(j=1\) contributes:

g_1=\frac{0.4(0.5)g_0}{0.8} =\frac{0.2(0.75)}{0.8} =\frac{0.15}{0.8}=0.1875.

Update \(g_2\). The \(j=2\) severity probability is zero, so only \(j=1\) contributes:

g_2=\frac{0.4(0.5)g_1}{0.8} =\frac{0.2(0.1875)}{0.8} =\frac{0.0375}{0.8}=0.046875.

The answer is \((0.75,0.1875,0.046875)\). The geometric aggregate generating function \(G_N(F(z))=0.6/(1-0.4(0.5+0.5z))=0.75/(1-0.25z)\) has these same first three coefficients. The common mistake avoided is assigning the geometric frequency parameters as \(a=0,\ b=0.4\), which does not satisfy the defining probability ratio.

Worked Example: A Short Severity Support on a Longer Grid

A hypothetical aggregate model has \(N\sim\operatorname{Poisson}(2)\). Each claim has severity 1 grid unit with probability \(0.7\) and 2 grid units with probability \(0.3\). Calculate \(g_0\) through \(g_4\), and count the inner-loop terms used.

Method and inputs. A Poisson frequency with mean 2 has \(a=0,\ b=2\). Here \(f_0=0,\ f_1=0.7,\ f_2=0.3\), and \(f_j=0\) for \(j>2\), so \(J=2\). The initial value is \(g_0=e^{-2}\). Because \(f_0=0\), the denominator is 1.

Update in ascending order. Keeping a common factor of \(e^{-2}\) visible, the first entry is

g_1=2(0.7)g_0=1.4e^{-2}.

At \(k=2\), the contributions are from \(j=1,2\):

g_2=\frac{2}{2}\bigl(0.7g_1+2(0.3)g_0\bigr) =\bigl(0.7(1.4)+0.6\bigr)e^{-2} =1.58e^{-2}.

At \(k=3\), \(j=1,2\) are still the only nonzero severity terms:

g_3=\frac{2}{3}\bigl(0.7g_2+2(0.3)g_1\bigr) =\frac{2}{3}\bigl(0.7(1.58)+0.6(1.4)\bigr)e^{-2} =\frac{973}{750}e^{-2}.

At \(k=4\), the same two severity indices contribute:

g_4=\frac{2}{4}\bigl(0.7g_3+0.6g_2\bigr) =\frac12\left(0.7\left(\frac{973}{750}\right)+0.6(1.58)\right)e^{-2} =\frac{13921}{15000}e^{-2}.

Thus the computed array is \((e^{-2},1.4e^{-2},1.58e^{-2},(973/750)e^{-2},(13921/15000)e^{-2})\). For \(K=4\), the inner-loop lengths are \(1,2,2,2\), for a total of 7 terms. The common mistake avoided is looping over severity indices beyond \(J=2\), where all the inputs are zero and no additional calculation is needed.

Cost, Storage, and Practical Limits

For a severity support ending at \(J\), the number of summand updates needed to calculate through \(K\) is

\sum_{k=1}^{K}\min(k,J).

When \(J\geq K\), this is \(K(K+1)/2\), so the cost grows on the order of \(K^2\). When \(J<K\), it is \(J(J+1)/2+(K-J)J\), which grows on the order of \(KJ\). A finite and relatively short severity support can therefore make a long grid less expensive than the unrestricted quadratic case.

Storing the full output requires \(K+1\) aggregate probabilities. The severity inputs require \(J+1\) entries when the support is \(0,\ldots,J\). Thus the straightforward implementation uses \(O(K+J)\) storage, in addition to the frequency parameters. If only the most recent entries needed by the recursion are required and previous outputs need not be retained, a rolling buffer can reduce the aggregate working storage to at most \(\min(K,J)+1\) entries. That saving is unavailable if the full output array is required.

Keep probabilities and indices aligned throughout. An array element \(f_j\) is the probability of severity \(j\) grid units, not the probability of the \(j\)-th observed claim. Also distinguish the maximum requested aggregate index \(K\) from the severity support endpoint \(J\). Severities above \(K\) cannot contribute to \(g_0,\ldots,g_K\), but the zero-severity probability \(f_0\) still matters to both initialization and the recursion denominator.

Key takeaway: Supply the correct frequency parameters and severity array, initialize with \(g_0=G_N(f_0)\), and calculate each later entry in ascending order. For finite severity support, cap each inner loop at \(\min(k,J)\); estimate cost from the number of terms actually summed, not just the output-array length.

Check Your Understanding

For each problem, identify the required inputs and apply the appropriate ordered update.

A hypothetical portfolio has Poisson claim count with mean 1 and severity probabilities \(f_0=0.25,\ f_1=0.75\). Find \(g_0,g_1,g_2\). (Answer: \(e^{-0.75},\ 0.75e^{-0.75},\ 0.28125e^{-0.75}\))
A hypothetical claim count has \(P(N=n)=0.7(0.3)^n\), and severity probabilities \(f_0=0.2,\ f_1=0.8\). Give \(a,b\), then calculate \(g_0\) and \(g_1\). (Answer: \(a=0.3,\ b=0;\ g_0=35/47\approx0.7447,\ g_1=420/2209\approx0.1901\))
A hypothetical count is \(\operatorname{Binomial}(2,0.25)\), and severity has \(f_0=0.4,\ f_1=0.6\). Find \(a,b\), the denominator, and \(g_0,g_1\). (Answer: \(a=-1/3,\ b=1;\ 1-af_0=17/15;\ g_0=0.7225,\ g_1=0.255\))
A frequency is in the \((a,b,0)\) class with \(a=0.2,\ b=0.5\), and severity support ends at \(J=3\). How many inner-loop summand updates are required to calculate through \(K=10\)? (Answer: \(\sum_{k=1}^{10}\min(k,3)=1+2+3+3+3+3+3+3+3+3=27\) summand updates)
A hypothetical model uses an \((a,b,1)\) frequency with \(p_0=0.4,\ p_1=0.3,\ a=0.2,\ b=0.5\). Severity support ends at \(J=2\), with \(f_2=0.25\), and the recursion denominator is \(0.9\). What correction term must be added to the \(k=2\) numerator before division? (Answer: \(cf_2=[p_1-(a+b)p_0]f_2=[0.3-0.7(0.4)](0.25)=0.005\))

What You'll Learn