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

Abstract

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.

1. Introduction

Recently, Robe-Voinea and Vernic (2016a, 2016b, 2017, 2018) and Vernic (2018) studied the recursive and Fast Fourier Transform (FFT)-based evaluation of the distribution of the following multivariate collective model:

$$(S_1,\dots ,S_m)=\left(\sum _{l=0}^{N_1}{U}_{1l}+\sum _{k=0}^{N_0}{L}_{1k},\dots ,\sum _{l=0}^{N_m}{U}_{ml}+\sum _{k=0}^{N_0}{L}_{mk}\right),m\ge 2,$$

(1)

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, ${\sum }_{l=0}^{N_j}{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 ${\left(U_{jl}\right)}_{l\ge 1}$ are non-negative, independent and identically distributed (i.i.d.) r.v.s, $1\le j\le m;$ they are also independent of the claim numbers and of the other claim sizes, $(L_{1k},\dots ,L_{mk})$ included; the random vectors ${(L_{1k},\dots ,L_{mk})}_{k\ge 1}$ are non-negative i.i.d. as the generic random vector $\mathbf{L}=(L_1,\dots ,L_m),$ and independent of the claim numbers, while the components of $\mathbf{L}$, however, are dependent; by convention, $U_{j0}=L_{j0}=0,1\le j\le 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:

$$\mathbf{S}=(S_1,\dots ,S_m)=\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}\sum _{i=0}^{N_{i_1\dots i_k}}{\mathbf{X}}_i^{\left(i_1\dots i_k\right)},$$

(2)

where

• The m-variate claim size random vectors ${\left({\mathbf{X}}_i^{\left(i_1\dots i_k\right)}\right)}_{i\ge 1}$ are i.i.d. as the generic m-variate random vector ${\mathbf{X}}^{\left(i_1\dots i_k\right)},$ whose jth univariate component $X_j^{\left(i_1\dots i_k\right)}=0$ if $j\notin \left\{i_1,\dots ,i_k\right\},$ meaning that ${\mathbf{X}}^{\left(i_1\dots i_k\right)}$ results from those claim events simultaneously affecting solely the lines $\left\{i_1,\dots ,i_k\right\}$; these events are counted by the r.v. $N_{i_1\dots i_k}$. Moreover, the ${\mathbf{X}}_i^{\left(i_1\dots i_k\right)}$s are also independent of the other claim size random vectors (i.e., of each ${\left({\mathbf{X}}_i^{\left(j_1\dots j_t\right)}\right)}_{i\ge 1}$, where $\left\{i_1,\dots ,i_k\right\}\ne \left\{j_1,\dots ,j_t\right\}$) and of the claim numbers. We let $X_{i,j}^{\left(i_1\dots i_k\right)}$ denote the jth univariate component of ${\mathbf{X}}_i^{\left(i_1\dots i_k\right)},$ $f_{i_1\dots i_k}$ the probability function (p.f.) of ${\mathbf{X}}^{\left(i_1\dots i_k\right)}$ (in the discrete case) and, by convention, ${\mathbf{X}}_0^{\left(i_1\dots i_k\right)}=0$.
• The components of the random vector number of claims $\mathbf{N}={\left(N_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}$ are dependent r.v.s, in total (maximum) $\nu =2^m-1.$

We adopt the actuarial terminology in which the distribution of $\mathbf{S}$ is called “compound” and the distribution of $\mathbf{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 $\mathbf{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 $\mathbf{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_{\mathbf{S}}$ the p.f. of $\mathbf{S}$, by g and $\phi$ 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, $\mathbf{n},\mathbf{t},\mathbf{x},\mathbf{y}$ are vectors whose corresponding dimension results from the context, $\mathbf{0}$ is the zero-vector, while the difference $\mathbf{x}-\mathbf{y}$ is componentwise. By $x_{+}={\displaystyle \sum _{i=1}^m}{x}_i$ we denote the sum of the components of the vector $\mathbf{x}$ and by $f^{\ast n}$ the n-fold convolution of f. To shorten the formulas, we rewrite the sum ${\sum }_{k=1}^m{\sum }_{1\le i_1<\dots <i_k\le m}$ as ${\displaystyle \sum _{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}$.

2. Evaluation of the Compound Distribution

We start by presenting the exact formula of the p.f. of $\mathbf{S}$ based on convolutions. This formula is so complex that, in general, it cannot be directly applied to find the distribution of $\mathbf{S}$.

Proposition 1.

The p.f. of the multivariate collective model (2) is given by

$$f_{\mathbf{S}}\left(\mathbf{x}\right)=\sum _{\mathbf{n}}Pr\left(\mathbf{N}=\mathbf{n}\right)\sum _{\begin{array}{ccc}\mathbf{x}=& \sum & {\mathbf{x}}_{i_1\dots i_k}\\ & 1\le k\le m& \\ & 1\le i_1<\dots <i_k\le m& \end{array}}\left(\prod _{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}{f}_{i_1\dots i_k}^{\ast n_{i_1\dots i_k}}\left({\mathbf{x}}_{i_1\dots i_k}\right)\right),$$

where $\mathbf{n}={\left(n_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m},n_{i_1\dots i_k}\in \mathbb{N}$.

Proof. 

We have

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& Pr\left(\mathbf{S}=\mathbf{x}\right)=\sum _{\mathbf{n}}Pr\left(\mathbf{N}=\mathbf{n}\right)Pr\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}\sum _{i=0}^{n_{i_1\dots i_k}}{\mathbf{X}}_i^{\left(i_1\dots i_k\right)}=\mathbf{x}\right)\hfill \\ & =& \sum _{\mathbf{n}}Pr\left(\mathbf{N}=\mathbf{n}\right)\sum _{\begin{array}{ccc}\mathbf{x}=& \sum & {\mathbf{x}}_{i_1\dots i_k}\\ & 1\le k\le m& \\ & 1\le i_1<\dots <i_k\le m& \end{array}}Pr\left(\bigcap _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}\left(\sum _{i=0}^{n_{i_1\dots i_k}}{\mathbf{X}}_i^{\left(i_1\dots i_k\right)}={\mathbf{x}}_{i_1\dots i_k}\right)\right),\hfill \end{array}$$

which immediately yields the result. ☐

We shall also need the pgf and the cf of $\mathbf{S}$.

Proposition 2.

The pgf and cf of the general multivariate collective model (2) are, respectively, given by

$$\begin{array}{ccc}\hfill g_{\mathbf{S}}\left(\mathbf{t}\right)& =& g_{\mathbf{N}}\left({\left(g_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}\right),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\phi }_{\mathbf{S}}\left(\mathbf{t}\right)& =& g_{\mathbf{N}}\left({\left({\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}\right).\hfill \end{array}$$

(4)

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

$$\begin{array}{ccc}\hfill g_{\mathbf{S}}\left(\mathbf{t}\right)& =& \mathbb{E}\left[\prod _{j=1}^m{t}_j^{S_j}\right]=\mathbb{E}\left[\prod _{j=1}^m{t}_j^{{\sum }_{k=1}^m{\sum }_{1\le i_1<\dots <i_k\le m}{\sum }_{i=0}^{N_{i_1\dots i_k}}{X}_{i,j}^{\left(i_1\dots i_k\right)}}\right]\hfill \\ & =& \mathbb{E}\left[\prod _{k=1}^m\prod _{1\le i_1<\dots <i_k\le m}\mathbb{E}\left[\left.\prod _{j=1}^m{t}_j^{{\sum }_{i=0}^{N_{i_1\dots i_k}}{X}_{i,j}^{\left(i_1\dots i_k\right)}}\right|\mathbf{N}\right]\right]\hfill \\ & =& \mathbb{E}\left[\prod _{k=1}^m\prod _{1\le i_1<\dots <i_k\le m}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}^{N_{i_1\dots i_k}}\left(\mathbf{t}\right)\right],\hfill \end{array}$$

hence the formula (3). ☐

2.1. 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 $\mathbf{S}$ under some supplementary assumptions. These assumptions are chosen such that the multivariate compound distribution of $\mathbf{S}$ can be rewritten as a compound distribution with a univariate counting distribution, for which we can apply the already existing recursions.

2.1.1. Case 1 Assumptions

As in Reference Robe-Voinea and Vernic (2017), we assume that $\mathbf{N}$ follows the multivariate Poisson distribution $MPo(\lambda ;\tilde{\lambda })$ with parameters $\lambda >0$ and $\tilde{\lambda }={\left({\lambda }_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}>\mathbf{0},$ having the pgf (see, e.g., Johnson et al. (1997))

$$g_{\mathbf{N}}\left(\mathbf{t}\right)=exp\left\{\lambda \left(\prod _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{t}_{i_1\dots i_k}-1\right)+\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{\lambda }_{i_1\dots i_k}\left(t_{i_1\dots i_k}-1\right)\right\}.$$

As a consequence, Proposition 2 easily yields the following pgf and cf

$$\begin{array}{ccc}\hfill g_{\mathbf{S}}\left(\mathbf{t}\right)& =& exp\left\{\lambda \left(\prod _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)+\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{\lambda }_{i_1\dots i_k}\left(g_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)\right\},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\phi }_{\mathbf{S}}\left(\mathbf{t}\right)& =& exp\left\{\lambda \left(\prod _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)+\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{\lambda }_{i_1\dots i_k}\left({\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)\right\}.\hfill \end{array}$$

(6)

Also, two recursive formulas for evaluating the distribution of $\mathbf{S}$ are obtained in the following proposition, where we denote by $f_{{\mathbf{X}}_{+}}$ the p.f. of the sum r.v. ${\sum }_{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}{\mathbf{X}}^{\left(i_1\dots i_k\right)}$.

Proposition 3.

Under the assumption that $\mathbf{N}\sim MPo(\lambda ;\tilde{\lambda })$ it holds that

$$\begin{array}{cc}{f}_{\mathbf{S}}\left(\mathbf{x}\right)& =\sum _{k=1}^m\sum _{\begin{array}{c}1\le i_1<\dots <i_k\le m\\ l\in \left\{i_1,\dots ,i_k\right\}\end{array}}{\lambda }_{i_1\dots i_k}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},\dots ,y_{i_k}\right)\le \left(x_{i_1},\dots ,x_{i_k}\right)\\ y_{i_{k+1}}=\dots =y_{i_m}=0\end{array}}\frac{y_l}{x_l}{f}_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\\ & \hspace{1em}+\lambda \sum _{0<\mathbf{y}\le \mathbf{x}}\frac{y_l}{x_l}{f}_{{\mathbf{X}}_{+}}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),x_l\ge 1,x_j\ge 0,\forall j\ne l,\end{array}$$

and

$$\begin{array}{cc}{f}_{\mathbf{S}}\left(\mathbf{x}\right)& =\frac{1}{x_{+}}\left[\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{\lambda }_{i_1\dots i_k}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},\dots ,y_{i_k}\right)\le \left(x_{i_1},\dots ,x_{i_k}\right)\\ y_{i_{k+1}}=\dots =y_{i_m}=0\end{array}}{y}_{+}{f}_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right.\\ & \left.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\lambda \sum _{0<\mathbf{y}\le \mathbf{x}}{y}_{+}{f}_{{\mathbf{X}}_{+}}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right],x_{+}\ge 1,\end{array}$$

(7)

with starting value

$$f_{\mathbf{S}}\left(\mathbf{0}\right)=exp\left\{\lambda \prod _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{f}_{i_1\dots i_k}\left(\mathbf{0}\right)+\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{\lambda }_{i_1\dots i_k}{f}_{i_1\dots i_k}\left(\mathbf{0}\right)-{\lambda }_{+}\right\},$$

where ${\lambda }_{+}=\lambda +{\sum }_{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}{\lambda }_{i_1\dots i_k}$. In the above formulas, $\mathbf{y}=\left(y_1,\dots ,y_m\right)$ is such that $\left(y_{i_1},\dots ,y_{i_m}\right)$ is a permutation of its components.

Proof. 

Due to the independence of the random vectors ${\mathbf{X}}^{\left(i_1\dots i_k\right)},$ we have that $g_{{\mathbf{X}}_{+}}={\displaystyle \prod _{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}};$ therefore, we can rewrite the pgf (5) as

$$\begin{array}{ccc}\hfill g_{\mathbf{S}}\left(\mathbf{t}\right)& =& exp\left\{{\lambda }_{+}\left(\frac{\lambda }{{\lambda }_{+}}\left(g_{{\mathbf{X}}_{+}}\left(\mathbf{t}\right)-1\right)+\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}\left(g_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)\right)\right\}\hfill \\ & =& exp\left\{{\lambda }_{+}\left(\frac{\lambda }{{\lambda }_{+}}{g}_{{\mathbf{X}}_{+}}\left(\mathbf{t}\right)+\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)\right\},\hfill \end{array}$$

meaning that in this case, the distribution of Model (2) is also a compound distribution, with a univariate Poisson counting distribution. More precisely, $\mathbf{S}$ can also be rewritten as

$$\mathbf{S}=\sum _{k=0}^{\tilde{N}}{\mathbf{C}}_k,$$

(8)

where $\tilde{N}\sim Po\left({\lambda }_{+}\right)$, ${\mathbf{C}}_0=\mathbf{0},$ while the random vectors ${\mathbf{C}}_1,{\mathbf{C}}_2,\dots$ are i.i.d. as the m-variate random vector $\mathbf{C}$ having the mixture p.f.

$$f_{\mathbf{C}}\left(\mathbf{x}\right)=\frac{\lambda }{{\lambda }_{+}}{f}_{{\mathbf{X}}_{+}}\left(\mathbf{x}\right)+\sum _{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}{f}_{i_1\dots i_k}\left(\mathbf{x}\right).$$

(9)

Regarding model (8), with $\tilde{N}$ satisfying Panjer’s recursion (see Panjer (1981)) with parameters $a,b\in \mathbb{R}$, i.e.,

$$Pr\left(\tilde{N}=n\right)=\left(a+\frac{b}{n}\right)Pr\left(\tilde{N}=n-1\right),\forall n\ge 1,$$

from Reference Sundt (1999) (see, also, formulas (15.4) and, respectively, (15.5) in Sundt and Vernic (2009)) it holds that

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& \frac{1}{1-a{f}_{\mathbf{C}}\left(\mathbf{0}\right)}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(a+b\frac{y_l}{x_l}\right)f_{\mathbf{C}}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),x_l\ge 1,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& \frac{1}{1-a{f}_{\mathbf{C}}\left(\mathbf{0}\right)}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(a+b\frac{y_{+}}{x_{+}}\right)f_{\mathbf{C}}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),\mathbf{x}>\mathbf{0}.\hfill \end{array}$$

(11)

Since in our case $\tilde{N}\sim Po\left({\lambda }_{+}\right),$ we have $a=0$ and $b={\lambda }_{+}$. Based on this, we insert Equation (9) into Equation(10) and obtain for $x_l\ge 1,$

$$f_{\mathbf{S}}\left(\mathbf{x}\right)={\lambda }_{+}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\frac{y_l}{x_l}\left(\frac{\lambda }{{\lambda }_{+}}{f}_{{\mathbf{X}}_{+}}\left(\mathbf{y}\right)+\sum _{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}{f}_{i_1\dots i_k}\left(\mathbf{y}\right)\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}).$$

We know that $X_j^{\left(i_1\dots i_k\right)}=0$ if $j\notin \left\{i_1,\dots ,i_k\right\},$ hence, concerning the argument $\mathbf{y}$ of $f_{i_1\dots i_k}\left(\mathbf{y}\right),$ we can take the components $y_{i_{k+1}}=\dots =y_{i_m}=0$. Therefore, if $l\notin \left\{i_1,\dots ,i_k\right\},$ clearly $y_l=0$ in the argument $\mathbf{y}$ of $f_{i_1\dots i_k}\left(\mathbf{y}\right),$ 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_{\mathbf{S}}\left(\mathbf{0}\right)=g_{\mathbf{S}}\left(\mathbf{0}\right)$ and from the above form of $g_{\mathbf{S}}$. This completes the proof. ☐

2.1.2. Case 2 Assumptions

Similarly to Robe-Voinea and Vernic (2016a, 2016b), the supplementary assumptions are now:

A1 

The p.f. of the total number of claims $N_{tot}={\sum }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{N}_{i_1\dots i_k}$ satisfies Panjer’s recursion for $a,b\in \mathbb{R}$.

A2 

Given $N_{tot}=n,$ the conditional distribution of the random vector number of claims $\mathbf{N}$ is assumed to be multinomial $Mnom(n;\mathbf{p})$ with parameters $n\in \mathbb{N}$ and $\mathbf{p}={\left(p_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m},$ where $p_{i_1\dots i_k}\in (0,1)$ such that ${\sum }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}=1.$ Therefore, with $\mathbf{n}={\left(n_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}$ and $n={\sum }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{n}_{i_1\dots i_k},$

$$Pr(\mathbf{N}=\mathbf{n}|N_{tot}=n)={\displaystyle \frac{n!}{{\prod }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{n}_{i_1\dots i_k}!}}\prod _{1\le k\le m;1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}^{n_{i_1\dots i_k}}.$$

Under these assumptions, the pgf, the cf of $\mathbf{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,

$$\begin{array}{ccc}\hfill g_{\mathbf{S}}\left(\mathbf{t}\right)& =& g_{N_{tot}}\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)\right),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\phi }_{\mathbf{S}}\left(\mathbf{t}\right)& =& g_{N_{tot}}\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)\right).\hfill \end{array}$$

(13)

Proof. 

To obtain the pgf formula, we recall that the pgf of the multinomial distribution $Mnom(n;\mathbf{p})$ is (see, e.g., Johnson et al. (1997)) $g\left(\mathbf{t}\right)={\left({\sum }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{t}_{i_1\dots i_k}\right)}^n$, so that the pgf of $\mathbf{N}$ becomes for $\mathbf{t}={\left(t_{i_1\dots i_k}\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m},$

$$\begin{array}{ccc}\hfill g_{\mathbf{N}}\left(\mathbf{t}\right)& =& \mathbb{E}\left[\mathbb{E}\left[\left.\prod _{k=1}^m\prod _{1\le i_1<\dots <i_k\le m}{t}_{i_1\dots i_k}^{N_{i_1\dots i_k}}\right|N_{tot}\right]\right]=\mathbb{E}\left[{\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{t}_{i_1\dots i_k}\right)}^{N_{tot}}\right]\hfill \\ & =& g_{N_{tot}}\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{t}_{i_1\dots i_k}\right).\hfill \end{array}$$

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

$$f_{\mathbf{S}}\left(\mathbf{0}\right)=g_{\mathbf{S}}\left(\mathbf{0}\right)=g_{N_{tot}}\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{f}_{i_1\dots i_k}\left(\mathbf{0}\right)\right),$$

the following recursive formula holds for $x_l\ge 1,x_j\ge 0,\forall j\ne l,$

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& K\sum _{k=1}^m\left\{\sum _{\begin{array}{c}1\le i_1<\dots <i_k\le m\\ l\in \left\{i_1,\dots ,i_k\right\}\end{array}}{p}_{i_1\dots i_k}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},\dots ,y_{i_k}\right)\le \left(x_{i_1},\dots ,x_{i_k}\right)\\ y_{i_{k+1}}=\dots =y_{i_m}=0\end{array}}\left(a+b\frac{y_l}{x_l}\right)f_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right.\hfill \\ & & \left.+a\sum _{\begin{array}{c}1\le i_1<\dots <i_k\le m\\ l\notin \left\{i_1,\dots ,i_k\right\}\end{array}}{p}_{i_1\dots i_k}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},\dots ,y_{i_k}\right)\le \left(x_{i_1},\dots ,x_{i_k}\right)\\ y_{i_{k+1}}=\dots =y_{i_m}=0\end{array}}{f}_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right\},\hfill \end{array}$$

(14)

while for $x_{+}>0,$

$$f_{\mathbf{S}}\left(\mathbf{x}\right)=K\sum _{{}_{\begin{array}{c}1\le k\le m\\ 1\le i_1<\dots <i_k\le m\end{array}}}{p}_{i_1\dots i_k}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},\dots ,y_{i_k}\right)\le \left(x_{i_1},\dots ,x_{i_k}\right)\\ y_{i_{k+1}}=\dots =y_{i_m}=0\end{array}}\left(a+b\frac{y_{+}}{x_{+}}\right)f_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),$$

(15)

where $K={\left[1-a{\sum }_{1\le k\le m;1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{f}_{i_1\dots i_k}\left(\mathbf{0}\right)\right]}^{-1}$ and $\mathbf{y}=\left(y_1,\dots ,y_m\right)$ is such that $\left(y_{i_1},\dots ,y_{i_m}\right)$ is a permutation of its components.

Proof. 

Considering the assumptions (A1 and A2), we rewrite Model (2) as

$$\mathbf{S}=\sum _{k=0}^{N_{tot}}{\mathbf{C}}_k,$$

where ${\mathbf{C}}_0=\mathbf{0},$ while the random vectors ${\mathbf{C}}_1,{\mathbf{C}}_2,\dots$ are i.i.d. as the m-variate random vector $\mathbf{C}$ with the p.f.

$$f_{\mathbf{C}}\left(\mathbf{y}\right)=\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{f}_{i_1\dots i_k}\left(\mathbf{y}\right).$$

(16)

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\ge 1,$

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& K\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(a+b\frac{y_l}{x_l}\right)\left(\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}{f}_{i_1\dots i_k}\left(\mathbf{y}\right)\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\hfill \\ & =& K\sum _{k=1}^m\sum _{1\le i_1<\dots <i_k\le m}{p}_{i_1\dots i_k}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(a+b\frac{y_l}{x_l}\right)f_{i_1\dots i_k}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}).\hfill \end{array}$$

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 $\mathbf{S}=(S_1,S_2,S_3)$ with

$$\begin{array}{ccc}\hfill S_1& =& \sum _{i=0}^{N_1}{X}_{i,1}^{\left(1\right)}+\sum _{i=0}^{N_{12}}{X}_{i,1}^{\left(12\right)}+\sum _{i=0}^{N_{13}}{X}_{i,1}^{\left(13\right)}+\sum _{i=0}^{N_{123}}{X}_{i,1}^{\left(123\right)},\hfill \\ \hfill S_2& =& \sum _{i=0}^{N_2}{X}_{i,2}^{\left(2\right)}+\sum _{i=0}^{N_{12}}{X}_{i,2}^{\left(12\right)}+\sum _{i=0}^{N_{23}}{X}_{i,2}^{\left(23\right)}+\sum _{i=0}^{N_{123}}{X}_{i,2}^{\left(123\right)},\hfill \\ \hfill S_3& =& \sum _{i=0}^{N_3}{X}_{i,3}^{\left(3\right)}+\sum _{i=0}^{N_{13}}{X}_{i,3}^{\left(13\right)}+\sum _{i=0}^{N_{23}}{X}_{i,3}^{\left(23\right)}+\sum _{i=0}^{N_{123}}{X}_{i,3}^{\left(123\right)}.\hfill \end{array}$$

For example, Equation (15) becomes

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& K\left\{\sum _{i=1}^3{p}_i\sum _{\begin{array}{c}1\le y_i\le x_i\\ y_2=y_3=0\end{array}}\left(a+b\frac{y_i}{x_{+}}\right)f_i\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right.\hfill \\ & & +\sum _{1\le i_1<i_2\le 3}{p}_{i_1{i}_2}\sum _{\begin{array}{c}\mathbf{0}<\left(y_{i_1},y_{i_2}\right)\le \left(x_{i_1},x_{i_2}\right)\\ y_{i_3}=0\end{array}}\left(a+b\frac{y_{i_1}+y_{i_2}}{x_{+}}\right)f_{i_1{i}_2}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\hfill \\ & & \left.+p_{123}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(a+b\frac{y_{+}}{x_{+}}\right)f_{123}\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y})\right\},\hfill \end{array}$$

(17)

where $K={\left[1-a\left({\sum }_{i=1}^3{p}_i{f}_i\left(\mathbf{0}\right)+{\sum }_{1\le i_1<i_2\le 3}{p}_{i_1{i}_2}{f}_{i_1{i}_2}\left(\mathbf{0}\right)+p_{123}{f}_{123}\left(\mathbf{0}\right)\right)\right]}^{-1}$.

2.1.3. 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. $\Theta$ takes the value $\theta ,$ the r.v.s $N_{i_1\dots i_k}$ are all i.i.d. Poisson distributed such that $N_{i_1\dots i_k\left|\Theta =\theta \right.}\sim Po\left(\theta {\lambda }_{i_1\dots i_k}\right),{\lambda }_{i_1\dots i_k}>0,1\le k\le m,1\le i_1<\dots <i_k\le m.$ Then, the pgf of $\mathbf{S}$ given $\Theta =\theta$ becomes, from Equation (3):

$$\begin{array}{ccc}\hfill g_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{t}\right)& =& g_{\mathbf{N}\left|\Theta =\theta \right.}\left({\left(g_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}\right)\hfill \\ & =& exp\left\{\theta {\lambda }_{+}\left(\sum _{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}{g}_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-1\right)\right\},\hfill \end{array}$$

where ${\lambda }_{+}={\sum }_{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}{\lambda }_{i_1\dots i_k}$. This is the pgf of a compound distribution with univariate Poisson $Po\left(\theta {\lambda }_{+}\right)$ counting distribution and multivariate claims distribution having p.f. $h={\sum }_{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}\frac{{\lambda }_{i_1\dots i_k}}{{\lambda }_{+}}{f}_{i_1\dots i_k};$ hence, the conditional distribution of $\mathbf{S}$, given $\Theta =\theta ,$ can be evaluated based on Equations (10) and (11), with $a=0$ and $b=\theta {\lambda }_{+}$. To find the unconditional distribution of $\mathbf{S},$ we use the technique described in Chapter 20 of Sundt and Vernic (2009). Therefore, with U denoting the distribution function of $\Theta$, we introduce the auxiliary functions

$$v_i\left(\mathbf{x}\right)={\int }_0^{\infty }{\theta }^i{f}_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{x}\right)dU\left(\theta \right),i=0,1,2,\dots ,$$

and note that $f_{\mathbf{S}}=v_0.$ Multiplying Equations (10) and (11) by ${\theta }^i$ and integrating yields the following two recursions for $v_i$

$$\begin{array}{ccc}\hfill v_i\left(\mathbf{x}\right)& =& \frac{{\lambda }_{+}}{x_l}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}{y}_{l}h\left(\mathbf{y}\right)v_{i+1}\left(\mathbf{x}-\mathbf{y}\right),x_l\ge 1,\hfill \\ \hfill v_i\left(\mathbf{x}\right)& =& \frac{{\lambda }_{+}}{x_{+}}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}{y}_{+}h\left(\mathbf{y}\right)v_{i+1}\left(\mathbf{x}-\mathbf{y}\right),x_{+}>0,\hfill \end{array}$$

with starting value $v_i\left(\mathbf{0}\right)={\int }_0^{\infty }{\theta }^i{e}^{\theta {\lambda }_{+}\left(h\left(\mathbf{0}\right)-1\right)}dU\left(\theta \right)$. Therefore, the algorithm for evaluating $f_{\mathbf{S}}\left(\mathbf{y}\right)$ for all $\mathbf{0}\le \mathbf{y}\le \mathbf{x}$ is more complex and implies the backward evaluation of all $v_i\left(\mathbf{y}\right),\mathbf{0}\le \mathbf{y}\le \mathbf{x},i=1,\dots ,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

$$\frac{d}{d\theta }lnu\left(\theta \right)=\frac{{\sum }_{i=0}^k{\eta }_i{\theta }^i}{{\sum }_{i=0}^k{\chi }_i{\theta }^i}.$$

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 $\mathbf{S}$ given $\Theta =\theta ,$

$${\phi }_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{t}\right)=exp\left\{\theta \left(\sum _{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}{\lambda }_{i_1\dots i_k}{\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)-{\lambda }_{+}\right)\right\},$$

where

$${\phi }_{\mathbf{S}}\left(\mathbf{t}\right)={\int }_0^{\infty }{\phi }_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{t}\right)dU\left(\theta \right).$$

(18)

Particular case: Simpler recursions are obtained when $\Theta$ is gamma $Ga\left(\delta ,\beta \right)$ distributed, with $\delta ,\beta >0$. In this case, the univariate mixed Poisson $Po\left(\theta {\lambda }_{+}\right)$ distribution becomes a Negative Binomial distribution $NB\left(\delta ,\frac{\beta }{\beta +{\lambda }_{+}}\right),$ which satisfies Panjer’s recursion with $a=\frac{{\lambda }_{+}}{\beta +{\lambda }_{+}}$ and $b=\left(\delta -1\right)\frac{{\lambda }_{+}}{\beta +{\lambda }_{+}}$. Since

$$f_{\mathbf{S}}\left(\mathbf{x}\right)={\int }_0^{\infty }{f}_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{x}\right)dU\left(\theta \right)=\sum _n{h}^{\ast n}\left(\mathbf{x}\right){\int }_0^{\infty }Pr\left({\tilde{N}}_{\theta }=n\right)dU\left(\theta \right)=\sum _{n}Pr\left(\tilde{N}=n\right)h^{\ast n}\left(\mathbf{x}\right),$$

where ${\tilde{N}}_{\theta }\sim Po\left(\theta {\lambda }_{+}\right),$ hence $\tilde{N}\sim NB\left(\delta ,\frac{\beta }{\beta +{\lambda }_{+}}\right),$ and it follows that we can use Equations (10) and (11) to obtain direct recursions for $f_{\mathbf{S}},$ i.e.,

$$\begin{array}{ccc}\hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& \frac{{\lambda }_{+}}{\beta +{\lambda }_{+}-{\lambda }_{+}h\left(\mathbf{0}\right)}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(1+\left(\delta -1\right)\frac{y_l}{x_l}\right)h\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),x_l\ge 1,\hfill \\ \hfill f_{\mathbf{S}}\left(\mathbf{x}\right)& =& \frac{{\lambda }_{+}}{\beta +{\lambda }_{+}-{\lambda }_{+}h\left(\mathbf{0}\right)}\sum _{\mathbf{0}<\mathbf{y}\le \mathbf{x}}\left(1+\left(\delta -1\right)\frac{y_{+}}{x_{+}}\right)h\left(\mathbf{y}\right)f_{\mathbf{S}}(\mathbf{x}-\mathbf{y}),\mathbf{x}>\mathbf{0},\hfill \end{array}$$

(19)

with starting value $f_{\mathbf{S}}\left(\mathbf{x}\right)={\left(\frac{\beta }{\beta +{\lambda }_{+}-{\lambda }_{+}h\left(\mathbf{0}\right)}\right)}^{\delta }.$ Moreover, regarding the cf, we easily obtain

$$\begin{array}{ccc}\hfill {\phi }_{\mathbf{S}}\left(\mathbf{t}\right)& =& \frac{{\beta }^{\delta }}{\Gamma \left(\delta \right)}{\int }_0^{\infty }{\phi }_{\mathbf{S}\left|\Theta =\theta \right.}\left(\mathbf{t}\right){\theta }^{\delta -1}{e}^{-\beta \theta }d\theta \hfill \\ & =& {\left(\frac{\beta }{\beta -{\sum }_{{}_{1\le k\le m;1\le i_1<\dots <i_k\le m}}{\lambda }_{i_1\dots i_k}{\phi }_{{\mathbf{X}}^{\left(i_1\dots i_k\right)}}\left(\mathbf{t}\right)+{\lambda }_{+}}\right)}^{\delta }.\hfill \end{array}$$

(20)

2.2. 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 $\mathbf{S}$. A much faster and less restrictive way to evaluate the p.f. of $\mathbf{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\left(\mathbf{x}\right)$ denote an m-variate function defined on the integer support $\underset{j=1}{\stackrel{m}{\times }}\left\{0,1,\dots ,r_j-1\right\}$; then its discrete Fourier transform, $\tilde{f}$, and, respectively, the inverse mapping, can defined by (definition consistent with the functions fftn and ifftn in Matlab)

$$\begin{array}{ccc}\hfill \tilde{f}\left(\mathbf{c}\right)& =& \sum _{x_1=0}^{r_1-1}\dots \sum _{x_m=0}^{r_m-1}f\left(\mathbf{x}\right)exp\left\{-2\pi i\sum _{j=1}^m\frac{x_j{c}_j}{r_j}\right\},c_j=0,\dots ,r_j-1,1\le j\le m,\hfill \\ \hfill f\left(\mathbf{x}\right)& =& \frac{1}{{\displaystyle \prod _{j=1}^m}{r}_j}\sum _{c_1=0}^{r_1-1}\dots \sum _{c_m=0}^{r_m-1}\tilde{f}\left(\mathbf{c}\right)exp\left\{2\pi i\sum _{j=1}^m\frac{x_j{c}_j}{r_j}\right\},x_j=0,\dots ,r_j-1,1\le j\le m.\hfill \end{array}$$

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:

FFT Algorithm for model (2)

Step 1. After setting the truncation point for each claim size random vector ${\mathbf{X}}^{\left(i_1\dots i_k\right)}$ at the same $\mathbf{r}=(r_1,\dots ,r_m)$, the corresponding truncated claim size distribution is obtained as ${\mathbf{f}}^{\left(i_1\dots i_k\right)}={\left[f_{i_1\dots i_k}\left(\mathbf{j}\right)\right]}_{\mathbf{0}\le \mathbf{j}\le \mathbf{r}-\mathbf{1}}$; if necessary, the resulting ${\mathbf{f}}^{\left(i_1\dots i_k\right)}$ will be filled with zeros (e.g., to constraint the $r_j$s to be powers of two).

Step 2. Apply the m-dimensional FFT to each ${\mathbf{f}}^{\left(i_1\dots i_k\right)}$, which results in the multidimensional table ${\tilde{\mathbf{f}}}^{\left(i_1\dots i_k\right)}.$

Step 3. Use Equation (4) in the general case to obtain the discrete cf ${\tilde{\phi }}_{\mathbf{S}}\left(\mathbf{j}\right)=g_{\mathbf{N}}\left({\left({\tilde{\mathbf{f}}}^{\left(i_1\dots i_k\right)}\left(\mathbf{j}\right)\right)}_{1\le k\le m;1\le i_1<\dots <i_k\le m}\right),\mathbf{0}\le \mathbf{j}\le \mathbf{r}-\mathbf{1}$.

Step 4. Apply the multidimensional IFFT to ${\tilde{\phi }}_{\mathbf{S}}$ to obtain the p.f. of $\mathbf{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 $\mathbf{N}\sim MPo(\lambda ;\tilde{\lambda }),{\tilde{\phi }}_{\mathbf{S}}$ is given by Equation (6);

-

Under the Case 2 assumptions (A1 and A2), ${\tilde{\phi }}_{\mathbf{S}}$ is given by Equation (13);

-

Under the Case 3 mixed Poisson assumption, ${\tilde{\phi }}_{\mathbf{S}}$ is given by Equation (18).

2.3. Numerical Illustration

In this section, we consider a particular trivariate model (2) with

$$\begin{array}{ccc}\hfill S_1& =& \sum _{i=0}^{N_1}{X}_{i,1}^{\left(1\right)}+\sum _{i=0}^{N_{12}}{X}_{i,1}^{\left(12\right)}+\sum _{i=0}^{N_{13}}{X}_{i,1}^{\left(13\right)}+\sum _{i=0}^{N_{123}}{X}_{i,1}^{\left(123\right)},\hfill \\ \hfill S_2& =& \sum _{i=0}^{N_2}{X}_{i,2}^{\left(2\right)}+\sum _{i=0}^{N_{12}}{X}_{i,2}^{\left(12\right)}+\sum _{i=0}^{N_{123}}{X}_{i,2}^{\left(123\right)},\hfill \\ \hfill S_3& =& \sum _{i=0}^{N_3}{X}_{i,3}^{\left(3\right)}+\sum _{i=0}^{N_{13}}{X}_{i,3}^{\left(13\right)}+\sum _{i=0}^{N_{123}}{X}_{i,3}^{\left(123\right)},\hfill \end{array}$$

for which we implemented both the recursive formulas and the FFT algorithm, under different assumptions.

As claim size distributions, we considered only type II Pareto distributions with the purpose to emphasize the effect of the exponential tilting on the FFT technique. We recall that the decumulative distribution (or survival) function of the m-variate type II Pareto distribution $P{a}_{m}II\left(\alpha ,{\left({\sigma }_i\right)}_{i=1,\dots ,m}\right),\alpha ,{\sigma }_i>0,i=1,\dots ,m,$ is given by

$$\overline{F}\left(\mathbf{x}\right)={\left(1+\sum _{i=1}^m\frac{x_i}{{\sigma }_i}\right)}^{-\alpha },x_i>0,i=1,\dots ,m.$$

The expected value of each marginal exists only if $\alpha >1$, while the variance exists only when $\alpha >2.$ We took (mainly from the numerical Example 4 in Reference Robe-Voinea and Vernic (2018))

$$\begin{array}{ccc}\hfill X^{\left(1\right)}& \sim & P{a}_{1}II\left(1,1\right),X^{\left(2\right)}\sim P{a}_{1}II\left(2,2\right),X^{\left(3\right)}\sim P{a}_{1}II\left(3,1\right),\hfill \\ \hfill {\mathbf{X}}^{\left(12\right)}& \sim & P{a}_{2}II\left(1.5,1,2\right),{\mathbf{X}}^{\left(13\right)}\sim P{a}_{2}II\left(2,1,1\right),{\mathbf{X}}^{\left(123\right)}\sim P{a}_{3}II\left(1.5;2,2,2\right).\hfill \end{array}$$

The expected value of $X^{\left(1\right)}$ and the variances of $X^{\left(2\right)},{\mathbf{X}}^{\left(12\right)},{\mathbf{X}}^{\left(13\right)},{\mathbf{X}}^{\left(123\right)}$ 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 $\mathbf{r}=\left(r,r,r\right)$ (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 $\theta =\left({\theta }_1={\theta }_2={\theta }_3\right)$ and noticed that an increasing of $\theta$ improves the results till $\theta =7/r,$ while a larger value like $\theta =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 $\mathbf{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 ${\mathbf{x}}_M=\left(x_M-1,x_M-1,x_M-1\right),$ 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

$$\begin{array}{ccc}\hfill AE& =& \sum _{\mathbf{0}\le \mathbf{x}\le {\mathbf{x}}_M}\left|f_{\mathbf{S}}\left(\mathbf{x}\right)-f_{\mathbf{S}}^{FFT}\left(\mathbf{x}\right)\right|,\hfill \\ \hfill Max.err& =& \underset{\mathbf{0}\le \mathbf{x}\le {\mathbf{x}}_M}{max}\left|f_{\mathbf{S}}\left(\mathbf{x}\right)-f_{\mathbf{S}}^{FFT}\left(\mathbf{x}\right)\right|.\hfill \end{array}$$

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_{\mathbf{S}}^{FFT}\left({\mathbf{x}}_M{)\ge F}_{\mathbf{S}}\right({\mathbf{x}}_M),$ an inequality caused by the AE that places compound mass below the truncation point.

Example 1.

We assume that $\mathbf{N}\sim MPo(\lambda ;\tilde{\lambda }),$ where $\lambda =1,{\lambda }_1=3,{\lambda }_2=3,{\lambda }_3=2,{\lambda }_{12}=2,{\lambda }_{13}=1.7,{\lambda }_{123}=1.5$; since for this particular model, the recursive method (we implemented Equation (7)) implies the evaluation of the p.f. $f_{{\mathbf{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. Example 1: Comparing recursive and FFT methods for $h=1,\theta =7/r$ and various $r,x_M$.

	${r}={{x}}_{{M}}=16$	${r}=32,{{x}}_{{M}}=20$	${r}=64,{{x}}_{{M}}=20$	${r}=128,{{x}}_{{M}}=20$
Rec. $F_{\mathbf{S}}\left({\mathbf{x}}_M\right)$	0.219737	0.312845	0.312845	0.312845
FFT $F_{\mathbf{S}}^{FFT}\left({\mathbf{x}}_M\right)$	0.219884	0.312909	0.312855	0.312847
FFT time up to $\mathbf{r}$	0.016 s	0.124 s	0.952 s	9.484 s
$AE$	1.4743 × ${10}^{-4}$	6.3771 × ${10}^{-5}$	1.0251 × ${10}^{-5}$	1.7488 × ${10}^{-6}$
$Max.err$	8.8571 × ${10}^{-8}$	2.0810 × ${10}^{-8}$	3.5244 × ${10}^{-9}$	6.9685 × ${10}^{-10}$

, 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\left(\lambda =5\right),$ for which we recall that $a=0,b=\lambda$ and $g_{N_{tot}}=e^{\lambda \left(t-1\right)}.$ 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. Example 2: Relative performances of the two methods when varying r ($\theta =7/r$).

	${r}=16$	${r}=32$	${r}=64$
Rec/FFT in time	12	130	781

, 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. Example 2: Comparing recursive and FFT methods for $h=1,\theta =7/r$ and various $r,x_M$.

	${r}={{x}}_{{M}}=16$	${r}={{x}}_{{M}}=32$	${r}={{x}}_{{M}}=64$	${r}=128,{{x}}_{{M}}=70$
Rec. $F_{\mathbf{S}}\left({\mathbf{x}}_M\right)$	0.80035	0.91543	0.96436	0.96804
FFT $F_{\mathbf{S}}^{FFT}\left({\mathbf{x}}_M\right)$	0.80039	0.91544	0.96436	0.96804
$AE$	4.3580 × ${10}^{-5}$	1.4294 × ${10}^{-5}$	3.8798 × ${10}^{-6}$	9.0937 × ${10}^{-7}$
$Max.err$	7.9393 × ${10}^{-7}$	1.8276 × ${10}^{-7}$	4.1642 × ${10}^{-8}$	9.6893 × ${10}^{-9}$

and the effect of changing the tilting parameters in

Table 4. Example 2: Comparing recursive and FFT methods for $h=1,r=64$ and various $\theta$.

	Recursion	FFT no tilt.	$\begin{array}{c}\mathbf{FFT}\mathbf{tilt}.\\ \mathit{\theta }=5/\mathit{r}\end{array}$	$\begin{array}{c}\mathbf{FFT}\mathbf{tilt}.\\ \mathit{\theta }=7/\mathit{r}\end{array}$	$\begin{array}{c}\mathbf{FFT}\mathbf{tilt}.\\ \mathit{\theta }=9/\mathit{r}\end{array}$
$F_{\mathbf{S}}\left({\mathbf{x}}_M\right)$	0.96436	0.96863	0.96439	0.96436	0.96436
$AE$		4.2693 × ${10}^{-3}$	2.8485 × ${10}^{-5}$	3.8798 × ${10}^{-6}$	5.4867 × ${10}^{-6}$
$Max.err$		4.5823 × ${10}^{-5}$	3.0770 × ${10}^{-7}$	4.1642 × ${10}^{-8}$	2.7813 × ${10}^{-8}$

, both supporting the above conclusions regarding the choices of r and $\theta .$

Example 3.

This example is related to Case 3, i.e., $\mathbf{N}$ follows a mixed Poisson distribution and, for simplicity, we let $\Theta \sim Ga\left(\delta ,\beta \right).$ Therefore, we implemented recursion (19) and the FFT based on Equation (20). The values of the parameters are: $\delta =\beta =2,{\lambda }_1={\lambda }_2=2.5,{\lambda }_3={\lambda }_{12}=2,{\lambda }_{13}=1.7,{\lambda }_{123}=1.5.$ The comparison between the two methods is presented in

Table 5. Example 3: Comparing recursive and FFT methods for $h=1,\theta =7/r$ and various $r,x_M$.

	${r}={{x}}_{{M}}=16$	${r}={{x}}_{{M}}=32$	${r}={{x}}_{{M}}=64$	${r}=128,{{x}}_{{M}}=70$
Rec. $F_{\mathbf{S}}\left({\mathbf{x}}_M\right)$	0.49044	0.72191	0.88701	0.90087
FFT $F_{\mathbf{S}}^{FFT}\left({\mathbf{x}}_M\right)$	0.49055	0.72200	0.88705	0.90088
$AE$	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}$

, 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.

3. 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.