# Background

## PD models

Financial institutions use Probability of Default (PD) models for various purposes such as client acceptance, provisioning and regulatory capital calculation as required by the Basel accords and the European Capital requirements regulation and directive (CRR/CRD IV).

A PD model is supposed to calculate the probability that a client defaults on its obligations within a one year horizon.

## Backtests

To test whether a model is performing as expected so-called backtests are performed. One such a backtest would be to calculate how likely it is to find the actual number of defaults at or beyond the actual deviation from the expected value (the sum of the client PD values). If this probability turns out to be below a certain threshold the model will be rejected.

For this procedure one would need the CDF of the distribution of the sum of $n$ Bernoulli experiments, each with an individual, potentially unique PD.

# TL;DR: The Question

Does Mathematica have a built-in distribution that describes the sum of a number of Bernoulli draws each with its own probability?

Since many financial institutions divide their portfolios in buckets in which clients have identical PDs, can we optimize the calculation for this situation?

• I would suggest to rename "p-value" to "probability", so that there is no confusion with the ugly concept from frequentist statistics... :) – gwr Nov 18 '16 at 9:39
• @gwr Done, that's indeed better. – Sjoerd C. de Vries Nov 18 '16 at 9:41

This kind of backtest is often performed using approximations, for instance with normal distributions, which may not be always valid. An exact test would be very nice to have.

The probability distribution that is necessary to describe this process is the so-called Poisson Binomial distribution. It describes the sum of independent Bernoulli experiments outcomes (which can be either 0 or 1). Unfortunately, this is not one of the many distributions built into Mathematica.

## Using built-in probability distribution functions

A naive approach would be to calculate the necessary CDF or PDF using TransformedDistribution:

PoissonBinomialDistribution[p_List] :=
TransformedDistribution[
Sum[x[i], {i, Length[p]}],
Table[x[i], {i, Length[p]}] \[Distributed]
(BernoulliDistribution /@ p)
]
]


Trying it for two portfolios, one with all PDs the same and one where the PDs are slightly diffferent:

fixedPDs = ConstantArray[0.05, 12];
randomPDs = RandomVariate[UniformDistribution[{0.04999, 0.05001}], 12];

PDF[PoissonBinomialDistribution[fixedPDs], 6] // AbsoluteTiming
(* {0.0029157, 0.0000106129} *)

PDF[PoissonBinomialDistribution[randomPDs], 6] // AbsoluteTiming
(* {8.98713, 0.0000106155} *)


Interestingly, the timings differ greatly. Mathematica seems to recognize that the first case can be treated as a BinomialDistribution and is fast. The timing of the case with varying PDs (the case we're interested in) is atrocious, though.

The problem is, the number of terms involved in the calculation grows as $\binom{n}{k}$ with $n$ the number of clients and $k$ the number of defaults. In many practical cases this will be prohibitive large. It will take days if not more to produce a result and may take up all available memory.

Here the results for a few (actually rather small) portfolios:

(note the logarithmic scale).

## Fernandez and Williams

A better approach by Fernandez and Williams (2010) is described in the Wikipedia page linked to above.

It is based on the discrete Fourier transform and the PDF can be written simply as:

ClearAll[PoissonBinomialPDFWiki]
PoissonBinomialPDFWiki[p_List, k_Integer] :=
With[{n = Length[p]},
With[{c = Exp[(2*I*Pi)/(n + 1)]},
Re[(1/(n + 1))*Sum[Product[1 + (c^l - 1)*p[[m]], {m, 1, n}]/c^(l*k),
{l, 0, n}]]]]


It is much faster:

PoissonBinomialPDF[fixedPDs, 6] // AbsoluteTiming
(* {0.000563269, 0.0000106129} *)

PoissonBinomialPDF[randomPDs, 6] // AbsoluteTiming
(* {0.000611714, 0.0000106155} *)


## Hong

Even better is the approach by Yili Hong (2011/2013), in On computing the distribution function for the Poisson binomial distribution. Computational Statistics and Data Analysis 59 (2013) 41–51 (journal link, preprint link). As a bonus, it calculates the whole PDF (all possible values) using the Fast Fourier Transform in one sweep. A direct implementation of the algorithm described by Hong would be:

ClearAll[PoissonBinomialPDF]
PoissonBinomialPDF[p_List, k_] :=
Module[{absz, argz, a, b, d, ω = 2 π/(Length[p] + 1), r1,
r2, len = Length[p], f},
absz[j_, l_] := Sqrt[(1 - p[[j]] + p[[j]] Cos[ω l])^2 + (p[[j]] Sin[ω l])^2];
argz[j_, l_] := ArcTan[1 - p[[j]] + p[[j]] Cos[ω l], p[[j]] Sin[ω l]];
d[l_] := Exp[Sum[Log[absz[j, l]], {j, len}]];
a[l_] := d[l] Cos[Sum[argz[j, l], {j, len}]];
b[l_] := d[l] Sin[Sum[argz[j, l], {j, len}]];
r1 = Table[a[i] + b[i] I, {i, Ceiling[len/2]}];
r2 = Table[r1[[len + 1 - i]]\[Conjugate], {i, Ceiling[len/2] + 1, len}];
f = Fourier[Join[{1}, r1, r2]/(len + 1), FourierParameters -> {1, -1}];
If[k < 0 || k > Length[f] - 1, 0, f[[k + 1]]] // Chop
]


Using Mathematica's vectorization capabilities this can be made much faster, by rewriting it as follows:

ClearAll[PoissonBinomialPDF, PoissonBinomialCDF]
PoissonBinomialPDF[p_List] :=
With[{len = Length[p]},
Module[{z, a, b, d, ω, r1, r2},
ω = 2 π/(len + 1);
z = Table[1 - p + p Cos[ω l] + I p Sin[ω l], {l, len}];
d = Times @@@ Abs[z];
a = d Cos[Total[Arg@z, {2}]];
b = d Sin[Total[Arg@z, {2}]];
r1 = a[[;; Ceiling[len/2]]] + b[[;; Ceiling[len/2]]] I;
r2 = Reverse[Take[r1, Floor[len/2]]]\[Conjugate];
Re[Fourier[Join[{1}, r1, r2], FourierParameters -> {-1, -1}]]]
]


For bucketized portfolios only two small changes to the code are necessary:

PoissonBinomialPDF[p_List, count_List] :=
Module[{z, a, b, d, ω, r1, r2, len = Total@count},
ω = 2 π/(len + 1);
z = Table[(1 - p + p Cos[ω l] + I p Sin[ω l])^count, {l, len}];
d = Times @@@ Abs[z];
a = d Cos[Total[Arg[z], {2}]];
b = d Sin[Total[Arg[z], {2}]];
r1 = a[[;; Ceiling[len/2]]] + b[[;; Ceiling[len/2]]] I;
r2 = Reverse[Take[r1, Floor[len/2]]]\[Conjugate];
Re[Fourier[Join[{1}, r1, r2], FourierParameters -> {-1, -1}]]
]


Taking care of a few syntax variations:

PoissonBinomialPDF[p_List, k_Integer] :=
Which[k < 0 || k > Length[p], 0, True, PoissonBinomialPDF[p][[k + 1]]]
PoissonBinomialPDF[p_List, count_List, k_Integer] :=
Which[k < 0 || k > Total[count], 0, True, PoissonBinomialPDF[p, count][[k + 1]]]

PoissonBinomialCDF[p_List] := Accumulate[PoissonBinomialPDF[p]]
PoissonBinomialCDF[p_List, k_Integer] :=
Which[k < 0, 0, k > Length[p], 1, True, PoissonBinomialCDF[p][[k + 1]]]
PoissonBinomialCDF[p_List, count_List] := Accumulate[PoissonBinomialPDF[p, count]]
PoissonBinomialCDF[p_List, count_List, k_Integer] :=
Which[k < 0, 0, k > Total[count], 1, True, PoissonBinomialCDF[p, count][[k + 1]]]


(note that the functions without an integer k in the argument list yield the whole PDF/CDF for values 0..k; be aware that the corresponding indices are 1..k+1)

This makes a huge difference:

(note the log scales both on the x and y axis now).

How does the Wikipedia version compare to Hong?

Clearly, the Hong version is to be preferred.

The bucketized version is indeed much faster than the standard version:

(* create buckets with clients and the same clients in a single, long list *)
fixedPDs = Table[Table[0.001 rep, {rep}], {rep, 100, 500, 100}] // Flatten;
fixedPDsbuck = Table[0.001 rep, {rep, 100, 500, 100}];
counts = Table[rep, {rep, 100, 500, 100}];

PoissonBinomialPDF[fixedPDs, 550] // AbsoluteTiming
(* {8.2356, 0.0221256} *)

PoissonBinomialPDF[fixedPDsbuck, counts, 550] // AbsoluteTiming
(* {0.327213, 0.0221256} *)


For larger portfolios gains are even higher.

## Compiled version of the Hong algorithm implementation

To squeeze out even more performance we can use Compile. In this case we have to be aware that the syntax version of Total used above (Total[Arg[z], {2}]) does not compile. It has to be replaced with Total /@ Arg[z]. Other than that the transformation is straightforward:

ClearAll[PoissonBinomialPDFCompiled1List, PoissonBinomialPDFCompiled2Lists,
PoissonBinomialPDF, PoissonBinomialCDF]
PoissonBinomialPDFCompiled1List =
Compile[{{p, _Real, 1}},
Module[{z, a, b, d, ω, r1, r2, len = Length[p]},
ω = 2 π/(len + 1);
z = Table[1 - p + p Cos[ω l] + I p Sin[ω l], {l, len}];
d = Times @@@ Abs[z];
a = d Cos[Total /@ Arg[z]];
b = d Sin[Total /@ Arg[z]];
r1 = a[[;; Ceiling[len/2]]] + b[[;; Ceiling[len/2]]] I;
r2 = Reverse[Take[r1, Floor[len/2]]]\[Conjugate];
Re[Fourier[Join[{1}, r1, r2], FourierParameters -> {-1, -1}]]
],
CompilationTarget -> "C"];

PoissonBinomialPDFCompiled2Lists =
Compile[{{p, _Real, 1}, {count, _Real, 1}},
Module[{z, a, b, d, ω, r1, r2, len = Total@count},
ω = 2 π/(len + 1);
z = Table[(1 - p + p Cos[ω l] + I p Sin[ω l])^count, {l, len}];
d = Times @@@ Abs[z];
a = d Cos[Total /@ Arg[z]];
b = d Sin[Total /@ Arg[z]];
r1 = a[[;; Ceiling[len/2]]] + b[[;; Ceiling[len/2]]] I;
r2 = Reverse[Take[r1, Floor[len/2]]]\[Conjugate];
Re[Fourier[Join[{1}, r1, r2], FourierParameters -> {-1, -1}]]
],
CompilationTarget -> "C"];


A slight change in the syntax version definitions so that one only has to use the functions names PoissonBinomialPDF and PoissonBinomialCDF:

PoissonBinomialPDF[p_List] := PoissonBinomialPDFCompiled1List[p];
PoissonBinomialPDF[p_List, k_Integer] :=
Which[k < 0 || k > Length[p], 0, True, PoissonBinomialPDFCompiled1List[p][[k + 1]]];
PoissonBinomialPDF[p_List, count_List] := PoissonBinomialPDFCompiled2Lists[p, count];
PoissonBinomialPDF[p_List, count_List, k_Integer] :=
Which[k < 0 || k > Total[count], 0,
True, PoissonBinomialPDFCompiled2Lists[p, count][[k + 1]]
]
PoissonBinomialCDF[p_List] := Accumulate[PoissonBinomialPDFCompiled1List[p]];
PoissonBinomialCDF[p_List, count_List] :=
Accumulate[PoissonBinomialPDFCompiled2Lists[p, count]];
PoissonBinomialCDF[p_List, k_Integer] :=
Which[k < 0, 0, k > Length[p], 1,
True, PoissonBinomialCDF[p][[k + 1]]
]
PoissonBinomialCDF[p_List, count_List, k_Integer] :=
Which[k < 0, 0, k > Total[count], 1,
True, PoissonBinomialCDF[p, count][[k + 1]]
]


Performance comparison:

We usually get an improvement by a factor of about five.

## Performing a backtest

A model backtest would then be performed by calculating the CDF for the given portfolio (using the PDs delivered by the model for each client) and comparing the actually counted number of defaults in a year with, for instance, the 2.5 and 97.5% levels for a 95% confidence interval:

In this example we would reject the model if the number of defaults happened to be less than 6 or greater than 19.

• A computational note: the sequence being FFT'd in PoissonBinomialPDF[] has so-called "conjugate even" symmetry, and thus necessarily has a real transform. There's a more efficient way to get the FFT without needing to generate the reflection r2; I'll see if it's easily adaptable to Mathematica without making the code too long. Meanwhile: why not simplify the evaluation of z to 1 - p + p Exp[I ω Range[len]]? – J. M. will be back soon Jan 3 '17 at 21:32
• @J.M.needshelp. That wouldn't work as p is a list here. Your code would result in a list of non-list elements whereas what is needed here is a list of lists. The Table expression can be written somewhat shorter as Table[1-p+p Exp[I ω l], {l, len}], though. – Sjoerd C. de Vries Jun 17 '18 at 12:13