On the Evaluation of the Distribution of a General Multivariate Collective Model: Recursions versus Fast Fourier Transform

: With the purpose of introducing dependence between different types of claims, multivariate collective models have recently gained a lot of attention. However, when it comes to the evaluation of the corresponding compound distribution, the problems increase with the dimensionality of the model. In this paper, we consider a multivariate collective model that generalizes a model already studied from the point of view of recursive and FFT evaluation of its distribution, and we extend the same study to the general model. With the intention to see which method works better for this general model, we compare the recursive method with the FFT technique, and emphasize the advantages and drawbacks of each one, based on numerical examples.


Introduction
Recently, Robe-Voinea and Vernic (2016aVernic ( , 2016bVernic ( , 2017Vernic ( , 2018 and Vernic (2018) studied the recursive and Fast Fourier Transform (FFT)-based evaluation of the distribution of the following multivariate collective model: which may arise in different contexts (see, e.g., the discussion in Section 14.1 of Reference Sundt and Vernic (2009)), from which we mention the case where a policyholder has m types of policies, such as auto, home, business, etc., that can be simultaneously affected by some claim events, such as floods, storms or earthquakes. More precisely, in this case, ∑ N j l=0 U jl denotes the aggregate claims affecting solely the policy of type j, while N 0 denotes the random variable (r.v.) number of claims simultaneously affecting all m types of policies, with L jk denoting the size of the kth such claim corresponding to the policy of type j. The assumptions under which this model was considered are: Each set of claim sizes (U jl ) l≥1 are non-negative, independent and identically distributed (i.i.d.) r.v.s, 1 ≤ j ≤ m; they are also independent of the claim numbers and of the other claim sizes, (L 1k , ..., L mk ) included; the random vectors (L 1k , ..., L mk ) k≥1 are non-negative i.i.d. as the generic random vector L = (L 1 , ..., L m ), and independent of the claim numbers, while the components of L, however, are dependent; by convention, U j0 = L j0 = 0, 1 ≤ j ≤ m.
Note that the above model assumes that a claim event affects either a single type of insurance line or all the insurance lines at once; there is no middle way, i.e., an event cannot affect only, say, lines 1 and 2, without causing claims in the other lines.
To overcome this drawback, in this paper we consider the more general multivariate collective model: where • The m-variate claim size random vectors X (i 1 ...i k ) i i≥1 are i.i.d. as the generic m-variate random vector X (i 1 ...i k ) , whose jth univariate component X (i 1 ...i k ) j = 0 if j / ∈ {i 1 , ..., i k } , meaning that X (i 1 ...i k ) results from those claim events simultaneously affecting solely the lines {i 1 , ..., i k }; these events are counted by the r.v. N i 1 ...i k . Moreover, the X The components of the random vector number of claims N = N i 1 ...i k 1≤k≤m;1≤i 1 <...<i k ≤m are dependent r.v.s, in total (maximum) ν = 2 m − 1.
We adopt the actuarial terminology in which the distribution of S is called "compound" and the distribution of N is called "counting".
To evaluate the distribution of this model, we shall consider that all the claim distributions are of the discrete type (e.g., they have been previously discretized; this is a usual assumption for collective models). We start the next section by presenting the exact formula of the p.f. of S based on convolutions, which, unfortunately, is unpractical. Therefore, we also aim at developing recursions for the evaluation of this distribution, an approach that requires the introduction of supplementary assumptions under which it is possible to obtain recursive formulas; examples of such recursions are given in Section 2.1. Apart from the restrictive assumptions, another important drawback of recursions is that they become very time consuming when the dimensionality m of the model increases (see the numerical examples in Section 2.3). To overcome these drawbacks, in Section 2.2 we propose the use of the Fast Fourier Transform (FFT) technique, which can be applied whenever we know the form of the characteristic function of S and which is very efficient when we want to evaluate the distribution's tail. However, this remarkably fast method is an approximate one, and we must pay a special attention to its specific errors; this aspect is illustrated by the numerical examples discussed in Section 2.3.
For simplicity, let us introduce more notation: We denote by f S the p.f. of S, by g and ϕ the probability generating function (pgf) and the characteristic function (cf), respectively, of a r.v., which will be indexed with the r.v.'s name. Also, n, t, x, y are vectors whose corresponding dimension results from the context, 0 is the zero-vector, while the difference x − y is componentwise.
x i we denote the sum of the components of the vector x and by f * n the n-fold convolution of f . To shorten the formulas, we rewrite the sum ∑ m k=1 ∑ 1≤i 1 <...<i k ≤m as ∑ 1≤k≤m 1≤i 1 <...<i k ≤m .

Evaluation of the Compound Distribution
We start by presenting the exact formula of the p.f. of S based on convolutions. This formula is so complex that, in general, it cannot be directly applied to find the distribution of S. Proposition 1. The p.f. of the multivariate collective model (2) is given by where n = (n i 1 ...i k ) 1≤k≤m;1≤i 1 <...<i k ≤m , n i 1 ...i k ∈ N.
Proof. We have which immediately yields the result.
We shall also need the pgf and the cf of S. Proposition 2. The pgf and cf of the general multivariate collective model (2) are, respectively, given by Proof. We prove only the pgf formula (the one for the cf follows along the same lines). Considering the independence assumptions of the model, we have hence the formula (3).

Recursive Evaluation
Due to the difficulty of directly applying the exact formula from Proposition 1, we present in the following examples of alternative recursive formulas for obtaining the p.f. of S under some supplementary assumptions. These assumptions are chosen such that the multivariate compound distribution of S can be rewritten as a compound distribution with a univariate counting distribution, for which we can apply the already existing recursions.
Proof. Due to the independence of the random vectors X (i 1 ...i k ) , we have that g X + = ∏ 1≤k≤m;1≤i 1 <...<i k ≤m g X (i 1 ...i k ) ; therefore, we can rewrite the pgf (5) as meaning that in this case, the distribution of Model (2) is also a compound distribution, with a univariate Poisson counting distribution. More precisely, S can also be rewritten as Regarding model (8), withÑ satisfying Panjer's recursion (see Panjer (1981)) with parameters a, b ∈ R, i.e., from Reference Sundt (1999) (see, also, formulas (15.4) and, respectively, (15.5) in Sundt and Vernic (2009)) it holds that Since in our caseÑ ∼ Po (λ + ) , we have a = 0 and b = λ + . Based on this, we insert Equation (9) into Equation(10) and obtain for x l ≥ 1, We know that X .., i k } , hence, concerning the argument y of f i 1 ...i k (y) , we can take the components y i k+1 = ... = y i m = 0. Therefore, if l / ∈ {i 1 , ..., i k } , clearly y l = 0 in the argument y of f i 1 ...i k (y) , which yields the first stated formula. The second formula results in a similar way by inserting Equation (9) into Equation (11), while the starting value is immediate from f S (0) = g S (0) and from the above form of g S . This completes the proof.

A1
The p.f. of the total number of claims N tot = ∑ 1≤k≤m;1≤i 1 <...<i k ≤m N i 1 ...i k satisfies Panjer's recursion for a, b ∈ R. A2 Given N tot = n, the conditional distribution of the random vector number of claims N is assumed to be multinomial Mnom(n; p) with parameters n ∈ N and p = ( Under these assumptions, the pgf, the cf of S and two alternative recursive formulas are presented in the following. Proposition 4. Under the assumptions (A1 and A2), the pgf and cf of the general multivariate collective model (2) become, respectively, Proof. To obtain the pgf formula, we recall that the pgf of the multinomial distribution Mnom(n; p) is (see, e.g., Johnson et al. (1997) Inserting this into Equation (3) easily yields Equation (12). Equation (13) follows in a similar way, which completes the proof.
Proposition 5. Under the assumptions (A1 and A2) of Model (2), with starting value while for x + > 0, and y = (y 1 , ..., y m ) is such that y i 1 , ..., y i m is a permutation of its components.
Proof. Considering the assumptions (A1 and A2), we rewrite Model (2) as where C 0 = 0, while the random vectors C 1 , C 2 , ... are i.i.d. as the m-variate random vector C with the p.f.
We use again Equations (10) and (11). By inserting Equation (16) into Equation (10), the stated formula of the constant K is easily obtained and, for x l ≥ 1, Using reasoning similar with the one used in the proof of Proposition 3, we obtain Equation (14). Similarly, Equations (11) and (16) lead to Equation (15). This completes the proof. Particular case: m = 3. Let us now have a look at a recursive formula in the trivariate case, where the general Model (2) is S = (S 1 , S 2 , S 3 ) with For example, Equation (15) becomes where

Case 3 Assumptions
Another assumption under which recursive formulas already exist is the univariate mixed Poisson counting distribution. To this purpose, we assume that, given that a positive univariate r.v. Θ takes the value θ, the r.v.s Then, the pgf of S given Θ = θ becomes, from Equation (3): where λ + = ∑ 1≤k≤m;1≤i 1 <...<i k ≤m λ i 1 ...i k . This is the pgf of a compound distribution with univariate Poisson Po (θλ + ) counting distribution and multivariate claims distribution having p.f. h = ∑ 1≤k≤m;1≤i 1 <...<i k ≤m λ i 1 ...i k λ + f i 1 ...i k ; hence, the conditional distribution of S, given Θ = θ, can be evaluated based on Equations (10) and (11), with a = 0 and b = θλ + . To find the unconditional distribution of S, we use the technique described in Chapter 20 of Sundt and Vernic (2009). Therefore, with U denoting the distribution function of Θ, we introduce the auxiliary functions v i (x) = ∞ 0 θ i f S|Θ=θ (x)dU (θ) , i = 0, 1, 2, ..., and note that f S = v 0 . Multiplying Equations (10) and (11) by θ i and integrating yields the following two recursions for with starting value v i (0) = ∞ 0 θ i e θλ + (h(0)−1) dU (θ). Therefore, the algorithm for evaluating f S (y) for all 0 ≤ y ≤ x is more complex and implies the backward evaluation of all v i (y) , 0 ≤ y ≤ x, i = 1, ..., x + (here backward means by decreasing i, see, e.g., the algorithm in Section 20.4.1 in Reference Sundt and Vernic (2009)). Being very time consuming, we don't insist on this algorithm. However, we note that the recursions can be refined under the assumption that the mixing distribution U is of the continuous type, with the density denoted by u satisfying the condition This is also called Willmot's mixing distribution, see Reference Willmot (1993).

Remark 1.
In view of the FFT, we also display the formula of the cf of S given Θ = θ, β+λ + , and it follows that we can use Equations (10) and (11) to obtain direct recursions for f S , i.e., with starting value f S (x) = β β + λ + −λ + h(0) δ . Moreover, regarding the cf, we easily obtain

Fast Fourier Transform Evaluation
The recursive method is an exact one, but, as already mentioned in the introduction, it has some important drawbacks: It can be applied only on some particular models and it becomes quite slow with the increasing of the dimensionality of S. A much faster and less restrictive way to evaluate the p.f. of S is provided by the Fast Fourier Transform method, which is an approximate technique used to strongly reduce the computing time, especially when evaluating the distribution's tail. As an advantage, this method can be applied to any model as long as its cf (4) (on which it is based) has a closed form, even if there is no recursive formula available. Therefore, the FFT technique received special consideration in the actuarial literature (see, e.g., References Bühlmann (1984), Embrechts et al. (1993), Jin and Ren (2014) or Robe-Voinea and Vernic (2018)). It consists of an algorithm that computes the discrete Fourier transform of a multivariate function, as well as its inverse, extremely fast. Let f (x) denote an m-variate function defined on the integer support m × j=1 0, 1, ..., r j − 1 ; then its discrete Fourier transform,f , and, respectively, the inverse mapping, can defined by (definition consistent with the functions fftn and ifftn in Matlab) In general, the FFT method requires that the values r j are powers of two for all j. For the multivariate model (2), this algorithm becomes: (2) Step 1. After setting the truncation point for each claim size random vector X (i 1 ...i k ) at the same r = (r 1 , ..., r m ), the corresponding truncated claim size distribution is obtained as f (i 1 ...i k ) = f i 1 ...i k (j) 0≤j≤r−1 ; if necessary, the resulting f (i 1 ...i k ) will be filled with zeros (e.g., to constraint the r j s to be powers of two).

FFT Algorithm for model
Step 2. Apply the m-dimensional FFT to each f (i 1 ...i k ) , which results in the multidimensional tablẽ f (i 1 ...i k ) .
Step 3. Use Equation (4) in the general case to obtain the discrete cf Step 4. Apply the multidimensional IFFT toφ S to obtain the p.f. of S.
Usually, to find the optimal r j s, one gradually increases them until the differences between the actual solutions and the previous ones are under a certain threshold (e.g., we increase r j as 32, 64, 128, 256 etc.). However, when dealing with heavy tailed claim size distributions, the results of this method can be strongly affected by a specific error caused by the discrete Fourier transform, which consists of placing under the truncation point the compound probability mass which is in fact above this point. This so-called "aliasing error" (AE) can be significantly reduced by applying to the claim size distributions an exponential change of measure, hence, forcing the tails of these distributions to decrease at an exponential rate; this transformation is known under the name of "exponential tilting" (for more details on this transformation see, e.g., Reference Grübel and Hermesmeier (1999)).

Particular cases:
Under the particular assumptions considered in the previous section to allow for a recursive evaluation, one should use the following formulas at Step 3 of the above algorithm: -When N ∼ MPo(λ;λ),φ S is given by Equation (6); -Under the Case 2 assumptions (A1 and A2),φ S is given by Equation (13); -Under the Case 3 mixed Poisson assumption,φ S is given by Equation (18).

Numerical Illustration
In this section, we consider a particular trivariate model (2) with for which we implemented both the recursive formulas and the FFT algorithm, under different assumptions.
The expected value of X (1) and the variances of X (2) , X (12) , X (13) , X (123) do not exist, hence we can see the effect of the exponential tilting in the heavy-tailed case. To discretize these distributions, we used the method of rounding considering the span h = 1 (good enough for illustration, but not optimal, see the discussion in Reference Robe-Voinea and Vernic (2018)).
Concerning the FFT method, as discussed in Section 2.2, we increased the truncation point r = (r, r, r) (we took r 1 = r 2 = r 3 = r for simplicity) from 16 till 128 (unfortunately, r = 256 generated an "out of memory" warning), and noticed that r = 128 yielded enough accurate results (for our data) compared to the exact method (see Tables 1, 3 and 5). Moreover, we also varied the tilting parameter θ = (θ 1 = θ 2 = θ 3 ) and noticed that an increasing of θ improves the results till θ = 7/r, while a larger value like θ = 9/r doesn't significantly improve the results (see Table 4 in Example 2).
As expected, there is an important difference between the computing times requested by the two methods. This difference increases with the increasing of the truncation point r and becomes really huge for r = 32 in Example 1 and for r = 128 in Examples 2 and 3. Therefore, we decided to compare the resulting p.f.s only up to a certain right endpoint denoted by x M = (x M − 1, x M − 1, x M − 1) , even if the support of the FFT was much larger. Note that the discretization time was not taken into account in the displayed computing times since discretization is needed by both methods (the total discretization time up to r = 128 was about 160 s).
To emphasize the differences between the FFT and the recursive results, we used the cumulative distribution function (cdf), the AE and the maximum absolute error evaluated between the exact p.f. and the FFT one; these last two are defined, respectively, by We shall now present three examples based on the three particular cases considered in Section 2.1. From these examples, we also note that in cdf terms, F FFT S (x M ) ≥ F S (x M ), an inequality caused by the AE that places compound mass below the truncation point.
Example 1. We assume that N ∼ MPo(λ;λ), where λ = 1, λ 1 = 3, λ 2 = 3, λ 3 = 2, λ 12 = 2, λ 13 = 1.7, λ 123 = 1.5; since for this particular model, the recursive method (we implemented Equation (7)) implies the evaluation of the p.f. f X + (i.e., multivariate convolutions), the corresponding computing time increases tremendously with r. Therefore, starting with r = 32, we took only x M = 20, which needed about 30 minutes only for the convolution part. However, the FFT was ready in only a few s even for r = 128, see Table 1, where we also display a comparison of the accuracy of the two methods. This example clearly emphasizes the speed discrepancy between the two methods and the important advantage of the FFT speed. Example 2. We now assume that N tot follows a Poisson distribution Po (λ = 5) , for which we recall that a = 0, b = λ and g N tot = e λ(t−1) . Numerically, we took the multinomial parameters p 1 = 0.3, p 2 = 0.2, p 3 = 0.2, p 12 = 0.15, p 13 = 0.1, p 123 = 0.05. We implemented the recursive Equation (17) and performed it up to the maximum x M = 70 in about 35 min. The speed difference between the two methods can be seen in Table 2, where we displayed the relative computing times Rec/FFT (for r = 128, FFT took about 8 s). The accuracy comparison of the two methods is presented in Table 3 and the effect of changing the tilting parameters in Table 4, both supporting the above conclusions regarding the choices of r and θ.  Example 3. This example is related to Case 3, i.e., N follows a mixed Poisson distribution and, for simplicity, we let Θ ∼ Ga (δ, β) . Therefore, we implemented recursion (19) and the FFT based on Equation (20). The values of the parameters are: δ = β = 2, λ 1 = λ 2 = 2.5, λ 3 = λ 12 = 2, λ 13 = 1.7, λ 123 = 1.5. The comparison between the two methods is presented in Table 5, from where we note once again that a value of r = 128 is sufficient to obtain good enough results by FFT (at least for these data). Concerning the computing times, the values were similar with the ones obtained in Example 2, see Table 2. 1.0650 × 10 −4 8.3639 × 10 −5 3.5553 × 10 −5 6.4368 × 10 −6 Max.err 1.6674 × 10 −7 3.2871 × 10 −8 6.8457 × 10 −9 1.5338 × 10 −9

Conclusions
In this paper, we proposed a general multivariate collective model that allows for dependence between the r.v.s number of claims, and, moreover, between the different r.v.s claim sizes. Since the evaluation of the resulting compound distribution is not straightforward, we discussed two types of techniques to deal with it: The recursive method that was presented in Section 2.1 and the FFT algorithm that was described in Section 2.2. Unfortunately, even if the recursive method has the advantage of being exact, it has two main drawbacks compared with the FFT method: First, recursions are available under some restrictive assumptions and second, they become very slow with the increasing of the dimensionality of the model. On the other hand, the main drawback of the FFT method consists in its specific errors, especially the aliasing error. However, the FFT technique is so fast compared with the exact recursions, that it is quite worthwhile to use it, especially when values from the tail of the compound distribution are needed (nevertheless, it is important to pay attention when choosing optimal values for the truncation points and for the tilting parameters). Another advantage of the FFT is that specific functions are already implemented in existing software, even for higher dimensions, with, eventually, the disadvantage of memory limitation.
To conclude, we would recommend the following approach: If recursive formulas are available for the considered model, they should be used to evaluate the compound distribution until some reasonable (in computing time terms) upper limit is reached, and then the FFT method should be applied for a more extended domain; to validate the accuracy of the FFT results, they should be compared with the ones obtained by the recursive method.
Funding: This research received no external funding.