The Discrete Fourier Transform and Cyclical Overflow

Introduction

A common tool of casualty actuaries is the collective-risk model S = X₁ + . . . + X_N, according to which aggregate loss S is composed of a random number N of independent, identically distributed claims X. Although there are several techniques for deriving probability distributions for S, the discrete Fourier transform (DFT) is arguably the most elegant and powerful when severity X is both discrete and finite. However, it suffers from two drawbacks. First, its use of complex numbers may daunt some actuaries; certainly it does not easily program into a spreadsheet.^[1] And, second, the finite support of S raises the issue of overflow, especially when the support of claim count N is unlimited. Since one can set an arbitrarily large common support for X and S, this issue is not a practical problem;^[2] however, how the true distribution of S overflows the DFT-derived distribution is an interesting topic. Moreover, understanding the overflow issue may enable and encourage others to employ the DFT. It is to this end that we will start in Section 2 with the necessary complex algebra. From there, in Section 3 we will define abstractly the DFT as an invertible transformation (a finite matrix of complex numbers). Bringing in probability, we will in Section 4 derive the DFT of S = X₁ + . . . + X_N in terms of N and X. Finally, in Section 5 we will show what happens to overflow in the DFT, relating it back to the complex algebra of Section 2.

The n^th roots of unity

Euler’s formula e^iθ = cos θ + i sin θ forms the basis for the n^th roots of unity.^[3] The magnitude of z = e^iθ is |z| = cos² θ + sin² θ = 1. Hence, e^iθ is a point on the complex unit circle, and the function z(θ) = e^iθ can be imagined as a point moving counterclockwise around this circle as θ increases. Plugging θ = 2π into Euler’s formula, we have e^{i • 2π} = cos 2π + i sin 2π = 1. Since the circumference of a circle measures 2π radians, the exponential function is periodic along the imaginary axis.

Now consider the equation zⁿ = 1, where n is a positive integer. For any integer k,

\[ \left(e^{i \frac{2 \pi k}{n}}\right)^{n}=e^{i \frac{2 \pi k}{n} n}=e^{i \cdot 2 \pi k}=\left(e^{i \cdot 2 \pi}\right)^{k}=(1)^{k}=1 \]

So numbers of the form \(\omega_{k}=e^{i \frac{2 \pi k}{n}}\) are n^th roots of unity. Furthermore, since \(\omega_{n}=e^{i \frac{2 \pi n}{n}}=e^{i \cdot 2 \pi}=1,\) within this formally infinite set only n elements are distinct, viz., \(\omega_{k}=e^{i \frac{2 \pi k}{n}}\) for k ∈ {0, 1, . . . , n − 1}. Multiplication of roots of unity involves modular arithmetic,^[4] since

\[ \begin{aligned} \omega_{k} \omega_{l}=e^{i \frac{2 \pi k}{n}} e^{i \frac{2 \pi l}{n}}=e^{i \frac{2 \pi(k+1)}{n}}= & \omega_{k+l}=\omega_{\bmod (k+l, n)} \\ & \Rightarrow k+1 \equiv \bmod (k+l, n) . \end{aligned} \]

Because an n^th degree polynomial has at most n distinct roots,^[5] there can be no other n^th roots of unity.

The operation of multiplication on G, the set of the n^th roots of unity, constitutes an algebraic group (Clark 1984, 17f). In addition to closure (if a, b ∈ G, then a · b ∈ G), multiplication on G has the three defining properties:

1. ∀a, b, c ∈ G (a · b) · c = a · (b · c)

2. ∃e ∈ G : ∀a ∈ G a · e = e · a = a

3. ∀a ∈ G ∃b ∈ G : a · b = b · a = e

The properties are (1) associativity, (2) identity existence, and (3) inverse existence. That multiplication is associative hardly needs to be mentioned. The identity element is ω₀ = 1. Abstractly, the identity element must be unique, for if e₁ and e₂ satisfy the identity property, then e₁ = e₁ · e₂ = e₂. Similarly, the inverse of an element must be unique, for, if both a · b₁ = b₁ · a = e and a · b₂ = b₂ · a = e, then

\[ b_{1}=b_{1} \cdot e=b_{1} \cdot\left(a \cdot b_{2}\right)=\left(b_{1} \cdot a\right) \cdot b_{2}=e \cdot b_{2}=b_{2} . \]

For group \(G\), the inverse of \(\omega_k\) is \(\omega_k^{-1}=\omega_{n-k}=\) \(\omega_{\text {mad }(n-k) n}\); the third part of the equation applies to the identity element: \(\omega_0^{-1}=\omega_0\). But, again from Euler’s formula:

\[ \begin{aligned} \omega_{k}^{-1}=e^{-i \frac{2 \pi k}{n}} & =\cos \left(-\frac{2 \pi k}{n}\right)+i \sin \left(-\frac{2 \pi k}{n}\right) \\ & =\cos \left(\frac{2 \pi k}{n}\right)-i \sin \left(\frac{2 \pi k}{n}\right)=\bar{\omega}_{k} . \end{aligned} \]

This means that the inverse of a root of unity is its complex conjugate.

Of all the \(n^{\text {mh }}\) roots of unity, \(\omega_1=e^{\frac{2 \pi}{3}}\) is the principal, from which the others can be derived by exponentiation \(\omega_k=\omega_1^k\). We will often drop the subscript from the principal root and designate all the roots as \(\omega_k=\omega_1^k=\omega^k\). The value of \(n\) will be understood; however, were it forgotten, it could always be reclaimed as the smallest positive integer for which \(\omega^r=1\).

Finally, as to the sum of the n^th roots of unity:

\[ (\omega-1) \sum_{k=0}^{n-1} \omega^{k}=\sum_{k=1}^{n} \omega^{k}-\sum_{k=0}^{n-1} \omega^{k}=\omega^{n}-\omega^{0}=1-1=0 . \]

Either ω = 1 or \(\sum_{k=0}^{n-1} \omega^{k}=0\) But the principal root is unity if and only if n = 1. Hence

\[ \sum_{k=0}^{n-1} \omega^{k}=\left\{\begin{array}{ll} 1 & \text { if } n=1 \\ 0 & \text { if } n>1 \end{array} .\right. \]

Intuitively, the formula for n > 1 is obvious from the equal dispersion of the roots about the circle; the origin is the center of their mass.^[6]

The discrete Fourier transform

The DFT is simply an n × n matrix Ω involving the n^th roots of unity ω^k. The jk^th element of Ω is ω^{j • k}. As in the previous section, ω is the principal n^th root of unity, and the indexing starts at zero. Define the two n × 1 vectors:

\[ \mathbf{n}=\left[\begin{array}{c} 0 \\ \vdots \\ n-1 \end{array}\right] \text { and } \mathbf{p}=\left[\begin{array}{c} p_{0} \\ \vdots \\ p_{n-1} \end{array}\right] \]

The DFT of \(\mathbf{p}\) is \(\boldsymbol{\Omega} \mathbf{p}\). Understanding the exponentiation as elementwise, we can express this as \(\boldsymbol{\Omega}=\omega^{\mathrm{nn}}\), which matrix is symmetric. But it is unnecessary to perform the exponentiation on \(\omega\). Due to the cyclicality of the multiplication, the elements of \(\mathbf{n n}^{\prime}\) can be reduced modulo \(n\) to the set \(\{0,1, \ldots, n-1\}\), at which point they are just as much \(n^{\text {th }}\)-root subscripts as they are exponents. In sum, \(\boldsymbol{\Omega}=\omega^{m n^{\prime}}=\omega^{m a d\left(i n n^{\prime}: m\right)}\).

To reclaim vector \(\mathbf{p}\), one needs the inverse DFT, or \(\boldsymbol{\Omega}^{-1}\), whose existence we will now show. Define the conjugate of matrix \(\boldsymbol{\Omega}\), or \(\overline{\boldsymbol{\Omega}}\), as the matrix of the conjugated elements of \(\boldsymbol{\Omega}\). According to Section 2 , the conjugate of a root of unity is the inverse of the root. So the \(j k^{\text {th }}\) element of \(\overline{\boldsymbol{\Omega}}\) is \(\omega^{-j k}\). Since \(\boldsymbol{\Omega}\) is symmetric, so too is \(\overline{\boldsymbol{\Omega}}\). Now the \(j k^{\text {th }}\) element of \(\overline{\boldsymbol{\Omega}} \boldsymbol{\Omega}\) is the dot product of the \(j^{\text {th }}\) row of \(\bar{\Omega}\) and the \(k^{\text {th }}\) column of \(\boldsymbol{\Omega}\). Hence:

\[ [\overline{\mathbf{\Omega}} \boldsymbol{\Omega}]_{j k}=\sum_{l=0}^{n-1} \bar{\omega}_{j l} \omega_{t k}=\sum_{l=0}^{n-1} \omega_{j l}^{-1} \omega_{l k}=\sum_{l=0}^{n-1} \omega^{-j l} \omega^{k}=\sum_{l=0}^{n-1}\left(\omega^{k-j}\right)^{l} . \]

The formula for the geometric series gives us

\[ \begin{aligned} {[\overline{\mathbf{\Omega}} \boldsymbol{\Omega}]_{j k} } & =\sum_{l=0}^{n-1}\left(\omega^{k-j}\right)^{l} \\ & =\left\{\begin{array}{ll} n & \text { if } \omega^{k-j}=1 \\ \frac{\left(\omega^{k-j}\right)^{n}-1}{\left(\omega^{k-j}\right)-1} & \text { if } \omega^{k-j} \neq 1 \end{array}=\left\{\begin{array}{ll} n & \text { if } \omega^{k-j}=1 \\ 0 & \text { if } \omega^{k-j} \neq 1 \end{array} .\right.\right. \end{aligned} \]

But \(\omega^{k-j}=1\) if and only if \(\bmod (k-j, n)=0\). Since \(j, k \in\{0,1, \ldots, n-1\}, \bmod (k-j, n)=0\) if and only if \(j=k\). So, finally,

\[ [\overline{\boldsymbol{\Omega}} \boldsymbol{\Omega}]_{j k}=\left\{\begin{array}{ll} n & \text { if } j=k \\ 0 & \text { if } j \neq k \end{array} \Rightarrow \overline{\boldsymbol{\Omega}} \boldsymbol{\Omega}=n \mathbf{I}_{n} .\right. \]

And since both \(\boldsymbol{\Omega}\) and \(\overline{\boldsymbol{\Omega}}\) are symmetric, \(\boldsymbol{\Omega} \overline{\boldsymbol{\Omega}}=\) \(\boldsymbol{\Omega}^{\prime} \overline{\boldsymbol{\Omega}}^{\prime}=(\overline{\boldsymbol{\Omega}} \boldsymbol{\Omega})^{\prime}=\left(n \mathbf{I}_n\right)^{\prime}=n \mathbf{I}_n\). Therefore, \(\boldsymbol{\Omega}\) is invertible: \(\boldsymbol{\Omega}^{-1}=\overline{\boldsymbol{\Omega}} / n .{ }\)^[7]

The discrete Fourier transform as a characteristic function

The relevance of the previous two sections will become clear in connection with the characteristic function, which is a moment generating function (mgf) with an imaginary argument. The characteristic function of any real-valued random variable \(X\), or \(\varphi_X\), is \(\varphi_x(t)=E\left[e^{i x}\right]\). Because \(\left|E\left[e^{i x}\right]\right|=\left|\int_{x=-\infty}^{\infty} e^{i x} d F_{F}(x)\right| \leq\) \(\int_{x=-\infty}^{\infty}\left|e^{i x \mid}\right| d F_{X}(x)=\int_{x=-\infty}^{\infty} 1 \cdot d F_{F}(x)=1,\) the characteristic function, unlike the mgf, always converges. However, since here we are concerned with discrete and finite random variables, whose mgfs do converge, this is of no advantage. Its advantage here lies in the simplicity and accuracy of its inversion.

First, we give a probability interpretation to vector \(\mathbf{p}\) of Section 3. The \(k^{\text {th }}\) element of \(\mathbf{p}_x\) will be the probability for random variable \(X\) to equal \(k\) : \(\operatorname{Pr} o b[X=k]=p_k\). So now the elements of \(\mathbf{p}_x\) must be non-negative real numbers whose sum is unity. Second, we define \(t_k=\frac{2 \pi k}{n}\) for \(k \in\{0,1, \ldots, n-1\}\). These \(n\) values of \(t\) are the radian measures on the unit circle of the \(n^{\text {th }}\) roots of unity, so \(\omega_k=e^{\frac{3+z}{7}}=e^a\). Accordingly,

\[ \varphi_{x}\left(t_{j}\right)=E\left[e^{i t_{j} X}\right]=\sum_{k=0}^{n-1} e^{i t_{j} k} p_{k}=\sum_{k=0}^{n-1}\left(\omega_{j}\right)^{k} p_{k}=\sum_{k=0}^{n-1} \omega^{j k} p_{k} . \]

The last expression is the \(j^{\text {th }}\) element of the matrix product \(\boldsymbol{\Omega}_{\mathbf{p}_X}\). In fact, defining \(\mathbf{t}\) as the \(n \times 1\) vector whose \(j^{\text {th }}\) element is \(t_j\), and applying \(\varphi_x\) elementwise, we have \(\varphi_X(\mathbf{t})=\boldsymbol{\Omega} \mathbf{p}_x\). Therefore, the DFT is a discrete form of the characteristic function, from which one can reclaim the probabilities as \(\mathbf{p}_x=\boldsymbol{\Omega}^{-1} \varphi_X(\mathbf{t})=\) \(\frac{\overline{\boldsymbol{\Omega}}}{n} \varphi_x(\mathbf{t})\).

Even so, it seems that we are just playing games, jumping from probabilities to the characteristic function and back again. Must not one know the probabilities before knowing the characteristic function? The answer is “Not necessarily.” The collective-risk model is a case in point, for if S = X₁ + . . . + X_N:

\[ \varphi_{s}(t)=E\left[e^{i s}\right]=E\left[e^{i\left(X_{1}+\cdots+x_{N}\right)}\right]=E\left[\prod_{k=-}^{N} e^{i x_{k}}\right] . \]

Separate N from the X variables by conditional expectation:

\[ E\left[\prod_{k=1}^{N} e^{i x_{k}}\right]=\underset{N}{E}\left[E\left[\prod_{k=1}^{N} e^{i i x_{k}} \mid N\right]\right]=\underset{N}{E}\left[\prod_{k=1}^{N} E\left[e^{i i x_{k}}\right] \mid N\right] . \]

Finally, simplify due to the identical distribution of the X variables:

\[ \begin{aligned} \underset{N}{E}\left[\prod_{k=1}^{N} E\left[e^{i X_{k}}\right] \mid N\right] & =\underset{N}{E}\left[\prod_{k=1}^{N} \varphi_{X}(t) \mid N\right] \\ & =\underset{N}{E}\left[\varphi_{x}(t)^{N} \mid N\right]=E\left[\varphi_{x}(t)^{N}\right] \end{aligned} \]

Therefore, \(\varphi_s(t)=E\left[\varphi_x(t)^N\right]\), in agreement with Heckman and Meyers (1983, 34), Klugman, Panjer, and Willmot (1998, 316), and Wang (1998, 864). The point is made: one can obtain the characteristic function of \(S\) without knowledge of its probabilities, and by DFT inversion (i.e., \(\mathbf{p}_s=\frac{\bar{\Omega}}{n} \varphi_s\) ) work back to the probabilities.

If \(N\) is a constant, or \(\operatorname{Prob}[N=k]=1\), then \(S=X_1+\cdots+X_k\) and \(\varphi_s(t)=\varphi_x(t)^k\). Therefore, \(\mathbf{p}_s=\frac{\overline{\boldsymbol{\Omega}}}{n} \varphi_s=\frac{\overline{\boldsymbol{\Omega}}}{n}\left(\boldsymbol{\Omega}_{\mathbf{p}_x}\right)^k\). This is the formula for the \(k^{\text {th }}\) convolution of \(X\), which holds for any non-negative integer \(k\).^[8] We will check it for \(k=0\). In this case \(\left(\boldsymbol{\Omega} \mathbf{p}_x\right)^k=1_n\), an \(n \times 1\) vector of ones, and \(\mathbf{p}_s=\frac{\overline{\boldsymbol{\Omega}}}{n} 1_n\), which contains the row averages of \(\overline{\boldsymbol{\Omega}}\). Its \(j^{\text {th }}\) element is

\[ \begin{aligned} {\left[\mathbf{p}_{s}\right]_{j} } & =\frac{1}{n} \sum_{l=0}^{n-1} \bar{\omega}_{j l}=\frac{1}{n} \sum_{l=0}^{n-1} \omega^{-j l}=\frac{1}{n} \sum_{l=0}^{n-1}\left(\omega^{-j}\right)^{l} \\ & =\frac{1}{n}\left\{\begin{array}{ll} n & \text { if } j=0 \\ \frac{\left(\omega^{-j}\right)^{n}-1}{\omega^{-j}-1} & \text { if } j>0 \end{array}\right. \\ & =\left\{\begin{array}{ll} 1 & \text { if } j=0 \\ 0 & \text { if } j>0 \end{array}\right. \end{aligned} \]

So Prob[S = 0] =1, as it must.

Again, if \(\operatorname{Prob}[X=1]=1, S=1_1+\cdots+1_N=N\). In this case \(\varphi_X(t)=E\left[e^{i x}\right]=e^{i t}\) and

\[ \varphi_{s}(t)=E\left[\varphi_{X}(t)^{N}\right]=E\left[\left(e^{i t}\right)^{N}\right]=E\left[e^{i N}\right]=\varphi_{N}(t) . \]

This means that \(S\) is distributed as \(N\); in the realm of the DFT, \(\mathbf{p}_s=\mathbf{p}_N\).

In general, \(E\left[z^N\right]\) is called the probability generating function (Klugman, Panjer, and Willmot 1998, 201) and denoted as \(P_N(z)\). Therefore, \(\varphi_s(t)=E\left[\varphi_X(t)^N\right]=\) \(\mathrm{P}_N\left(\varphi_x(t)\right)\). Klugman’s Appendix B contains the probability generating functions of many claim-count distributions, and our purpose here does not require us to duplicate his and others’ work. So we turn at last to the overflow problem.

The discrete Fourier transform and overflow

Example 4.11 of Klugman, Panjer, and Willmot (1998, 319f) demonstrates the use of the DFT to obtain the distribution of \(S=X_1+\cdots+X_N\), where \(N\) is Poisson-distributed with a mean of 3 , and \(\mathbf{p}_x=\left[\begin{array}{llllll}0 & 0.5 & 0.4 & 0.1 & 0 & \ldots\end{array}\right]^{\prime}\). The demonstration is performed twice, once over the support \(\left[\begin{array}{llll}0 & 1 & \ldots & 7\end{array}\right]^{\prime}\) and once over the support \(\left.\begin{array}{llll}0 & 1 & \ldots & 4,095\end{array}\right]^{\prime}\). For the first support, \(n=8\); for the second, \(n=4,096\). The purpose of the twofold demonstration is to show the importance of choosing a value for \(n\) large enough for the overflow to be negligible. Having replicated his calculations, we present the results in Table 1.

Klugman remarks:

. . . all the entries for n = 4,096 are accurate to the five decimal places presented. On the other hand, with n = 8, the FFT gives values that are clearly distorted. If any generalization can be made it is that more of the extra probability has been added to the smaller values of S. [p. 319]

In other words, the values for n = 8 are distorted by the overflow, but the lower the value of s, the greater the error. This is true enough, but more generalization can be made.

The probability for S to exceed 7, or to overflow the procedure for n = 8, is about 17.7 percent. The remarkable point is that the DFT preserves all the probability, regardless of n. Why does it not simply cut off accurately at s = 7? Instead, it compresses all the probability into the eight values. The explanation will appear from the following example.

Starting with a Bernoulli(0.5)-distributed X, we derived its next four convolutions and display them in Table 2. The probabilities of C(n) are binomial(n, 0.5).

Next, we calculated those probabilities according to the DFT with n = 4. Accordingly,

\[ \begin{array}{c} \mathbf{n}=\left[\begin{array}{l} 0 \\ 1 \\ 2 \\ 3 \end{array}\right] \mathbf{p}_{x}=\left[\begin{array}{l} 0.5 \\ 0.5 \\ 0.0 \\ 0.0 \end{array}\right] \\ \boldsymbol{\Omega}=\left[\begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{array}\right] \quad \mathbf{\Omega}^{-1}=\frac{1}{4}\left[\begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{array}\right] \end{array} \]

Hence, \(\mathbf{p}_{C(n)}=\boldsymbol{\Omega}^{-1}\left(\boldsymbol{\Omega}_{\mathbf{p}_x}\right)^n\). Table 3 shows the DFT result.

As expected, C(2) and C(3) agree. But overflow begins with C(4), and the form of the overflow is cyclical. What should have been in the k = 4 row added to the k = 0 row; what should have been in the k = 5 row added to the k = 1 row.

The DFT cannot ignore overflow; rather it must recycle it. Instead of calculating probabilities for S = X₁+ . . . + X_N, it calculates them for mod(S, n). Its arithmetic is as modular as that of the G group of Section 2:

\[ \omega_{k} \omega_{l}=\omega_{k+l}=\omega_{\bmod (k+l, n)} \Rightarrow k+l \equiv \bmod (k+l, n) . \]

Alternatively, one may imagine the support to be infinite, over the set of natural numbers \(\{0,1,2, \ldots\}\). However, there is no least root of unity; one must settle for some finite \(n\) and its principal root. So to match the roots of unity with the natural numbers, one must cycle endlessly through \(\omega^0 \quad \omega^1 \quad \ldots \quad \omega^{n-1}\). The roots of unity are like ID variables according to which the probabilities may be summarized, as in the summation:

\[ \sum_{k=0}^{\infty} \omega^{k} p_{k}=\sum_{k=0}^{\infty} \omega^{\bmod (k, n)} p_{k}=\sum_{k=0}^{n-1}\left(\omega^{k} \sum_{\bmod (l, n)=k} p_{l}\right)=\sum_{k=0}^{n-1}\left(\omega^{k} p_{k}^{*}\right) . \]

Summarizing the n = 4,096 probabilities of Table 1 according to mod(s, 8) will produce the n = 8 probabilities.

Conclusion

Many actuaries have neither the temperament nor the fortitude to penetrate the thicket of Σ operators found in many treatments of the DFT. Hopefully, they will find this treatment gentler, starting as it does with the n^th roots of unity. The addition of a dash of group theory, about which much more could be said,^[9] is suggestive of the modular arithmetic underlying the DFT and of the behavior of any overflow. The magic^[10] of the DFT is not dispelled by its being opened to examination. But to those who understand the preceding sections of this paper, any impediments to working DFT magic will pertain only to computer hardware and software.^[11]