Neumann Series on the Recursive Moments of Copula-Dependent Aggregate Discounted Claims

We study the recursive moments of aggregate discounted claims, where the dependence between the inter-claim time and the subsequent claim size is considered. Using the general expression for them−th order moment proposed by Léveillé and Garrido (Scand. Actuar. J. 2001, 2, 98–110), which takes the form of the Volterra integral equation (VIE), we used the method of successive approximation to derive the Neumann series of the recursive moments. We then compute the first two moments of aggregate discounted claims, i.e., its mean and variance, based on the Neumann series expression, where the dependence structure is captured by a Farlie–Gumbel–Morgenstern (FGM) copula, a Gaussian copula and a Gumbel copula with exponential marginal distributions. Insurance premium calculations with their figures are also illustrated.


Introduction
As the occurrence of catastrophe events becomes more frequent, the assumption of independence between event occurrence and claim severity is no longer sufficient in insurance risk modeling, given its impact on pricing and reserving, capital allocation solvency, as well as regulatory systems.The February 2009, Victorian bushfire in Australia (10,200 insurance claims amounting to approximately AUD 1.2 billion), the February 2011, Christchurch earthquake (USD 13 billion insured economic losses), the 2011 Great Eastern Japanese earthquake (loss amounting to as much as USD 40 billion), as well as the 2012 Hurricane Sandy (an expected loss of USD 25 billion) are the examples of this effect (see [1,2]).
In dealing with the dependency between the inter-claim arrivals and claim sizes, various approaches have been proposed in previous studies that can be noticed in [3][4][5][6][7][8][9][10][11], as well as the references therein.Regardless of the model used, we notice that previous research focused on either examining the expression of the moments of the aggregate discounted claims, Z(t), as can be seen in [6,[11][12][13][14], or by finding the related ruin measures and the ruin probability expressions, just like in [3][4][5]10,15].
Assuming the Poisson claim arrival process with claim sizes following mixed exponential distributions, [7] obtained the explicit expressions of the actuarial net premiums and the variances of the discounted aggregate claims from the Laplace transform of the distribution of the shot noise process, which was derived using the martingale approach.The first two moments of the aggregate discounted claims were obtained in [9] assuming the dependency between the claim sizes and the rates of claim occurrence affected by a Markovian environment, called the circumstance process.A delayed renewal process was also explored in [12][13][14][15], as well as [11] to accommodate the epochs between claim arrival and the observation of the risk process.
The asymptotical behaviour of a conditional tail probability dependence structure of claim sizes given the inter-claim arrival time was studied in [5,8].Assuming that the conditional tail of claim size given the inter-claim time satisfies a certain condition for a bounded inter-claim time and a really huge claim size, [5] obtained the asymptotic tail probabilities of the discounted aggregate claims.Three copulas were indicated as fulfilling this assumption, which are the Farlie-Gumbel-Morgenstern (FGM) and the Frank and Ali-Mikhail-Haq (AMH) copulas, and the Weibull claim size was paired with exponential inter-claim arrival time in their numerical example.On the other hand, [8] explored the analytical properties related to the same dependence structure described by the survival copulas, such as their local and global uniformity.
Conditioning on the first arrival and using a renewal theory argument, [12] derived a useful expression for the m−th recursive moment, whereby the inter-claim arrival time and the claim severity are assumed to be independent.The same conditioning argument was then applied in [6,10], assuming the FGM copula and then solved using the Laplace transform approach.More recently, [11] also adopted the same technique to derive the recursive moments of a Sparre Andersen risk process assuming a fairly general dependence structure between the inter-claim time and subsequent claim size variables, providing a simplified moments expression for assuming Erlang weights.Four types of copula were showcased in their examples of joint distribution between the said variables, which are the polynomial copula, the Bernstein copula, the generalized FGM copula, the extended FGM copula (references for these copulas are available in Section 3 of [11]).
The recursive moment equation resulting from the technique used in [6,12] takes the form of a Volterra integral equation of the second kind, which is widely used in the fields of mathematical physics, such as the electromagnetic and viscoelasticity fields, to represent the dynamics of materials that contain memory (refer to [16][17][18][19][20]).We are interested in using the same technique and then extend the recursive moments obtained in [6], so that it can be applied to any continuous bivariate distribution to accommodate the dependency between the two variables.To do so, we solve the recursive expression of the moments using the Neumann series obtained via the Picard method of successive approximations, upon which a selection of bivariate distributions can be applied, including bivariate copula.
This article is structured as follows.Section 2 will introduce the general framework of the continuous time renewal risk model together with its recursive moments with exponentially distributed inter-claim time and general claim size distribution.The dependency between the claim size and inter-claim time are then specified using a bivariate copula.For that purpose, we consider three copulas, which are the FGM copula, the Gaussian copula, which is a type of elliptical copula, and the Gumbel copula, an Archimedean type of copula, which is a natural candidate to represent an extreme value copula that caters for the one-sided dependence structure (see [21]).
In Section 3, we introduce the Volterra integral equation, which will be solved using the successive approximations method, leading to the Neumann series expression of the recursive moments, which is the main result of this paper.The Neumann series expression of the recursive moments allows the flexibility to capture various dependence structures provided by copula probability density functions (pdf).
Section 4 starts with the comparison between the value of moments obtained by our Neumann series expression assuming the FGM copula and the closed form solution by [6].We then present the numerical analysis, showing the value of moments across the dependence parameter for each copula considered, assuming an exponentially distributed claim size.The illustration and comparison of moments, as well as premium values under the standard deviation principle are also included in this section.Section 5 concludes the article.

Model Setup
We consider a continuous time renewal risk model as in [6], whereby Z = {Z(t)} t≥0 with: In this model, N = {N (t)} t≥0 is a homogeneous Poisson process and X i is a non-negative random variable (r.v.) representing the claim amount occurring at time T i for i = 1, 2, . . ., N (t).The instantaneous rate of net interest, δ, is assumed to be deterministic.We also define the inter-claim time variable r.v.W j as: The variables, X j and W j , are assumed to be continuous.In this study, we relax the independent assumption between the inter-claim time, W j , and the claim size, X j , and we let {(X j ,W j )} j∈N to form a sequence of independent and identically distributed (i.i.d) random vectors, whose components are dependent.

Recursive Moments of Aggregate Discounted Claims
Conditioning on the arrival of the first claim as in [6,10,12], and knowing that E(X m |W = s) = ∞ 0 x m f X|W =s (x)dx for m ≥ 1, we have the general form of the m−th moments of aggregate discounted claim as the following: where f X,W (x, s) is the bivariate probability density function (pdf) of the pair, X j and W j .
In this study, the joint pdf is described via a copula, C θ (u, v), whose pdf is given by with dependence parameter θ (see, e.g., [22,23] for a general review on copulas).The bivariate pdf of (X, W ) at (x, s) can be represented as: where f (•) and F (•) are the marginal pdf and cdf for r.v.'s X and W .

Copula Used
We are interested to calculate the first, second and mth moment of aggregate discounted claims under three copulas: the FGM copula, the Gaussian copula and the Gumbel copula.Their pdfs are given by: The FGM copula is used in this study due to its simplicity and analytical tractability.It is also used to verify our numerical results in Section 4.1 with [6].The well-known elliptical family member, the Gaussian copula, is chosen as, to the best of our knowledge, the effect of elliptical copula in terms of the dependence between claim size and inter-claim time have not been explored extensively.The Gumbel copula is also chosen, since it could be adopted by an insurance company, who assumes risks with extreme magnitude, having the tendency to occur together, as pointed out by De Matteis in [21].Many standard statistical texts offer illustrations of copula scatter plots with various dependence structure, for which we refer to [22][23][24].

Linear Integral Equations
The most general form of linear integral equation (IE) is given by: where Ψ(T ) is the solution to the IE that we need to obtain, g(T ) is a given function and K(T, s) is the kernel for the IE.Equation ( 8) can be a homogeneous/non-homogenous, Volterra/Fredholm IE of the 1st/the 2nd kind, for which readers are referred to the conditions given in Section 2.1 of [18].Linear IE can be solved either numerically using methods, such as the Runge-Kutta and collocation methods (see, e.g., [25,26]), or solved explicitly, such as by obtaining its Neumann series via the Picard method of successive approximations or using the Laplace transform method.

Volterra IE of the 2nd Kind
If we have g(T ) = 0, h(T ) = 1, and b(T ) = T , Equation ( 8) becomes: which is a non-homogeneous Volterra integral equation of the second kind.The Volterra IE is widely used in the areas of viscoelasticity and electromagnetic to compute the dynamics of materials that "contain" memory, other than being useful in renewal theory and demography (see, e.g., [16,27], as well as the references therein for a more rigorous treatment on Volterra integral equations).We easily notice that the moments provided by Equations ( 2) and ( 3) take the form of Equation ( 9) and attempt to derive the explicit solution of the recursive expressions using Neumann series in the next subsection.
A unique and continuous solution, Ψ(T ), is obtainable if we have a combination of a continuous kernel, K(T, s), in the region a ≤ s ≤ T ≤ b(T ) with a function, g(T ), that is continuous in the region a ≤ T ≤ b(T ), even though it is not a requirement for the kernel function, K(T, s), to be continuous (see page 1 of [28] and page 5 of [20]).For the case of a discontinuous kernel function, we need to check if K(T, s) fulfills the three regularity conditions set on page 3 of [27], and, hence is an L 2 -function.
In the case of the first and second moments, µ Z (T ) and µ Z (T ), the function, g(T ), is represented by the following equations, respectively: where F W (s) = 1 − e −βs is the inter-claim time cdf.As X and W are continuous r.v.'s, and by corollary 2.2.6 of [23] on copula continuity, g(T ) is the continuous function for s ∈ [0, T ] and x ∈ [0, ∞], since it is the sum and product of continuous functions.The kernel function is also continuous, as it is an exponential function given by: Additionally, it is a bounded function in the square Π = {(T, s) : a ≤ T ≤ b(T ), a ≤ s ≤ T }.

Neumann Series
In this section, we will find the Neumann series of the Volterra IE assuming the exponentially distributed inter-claim arrival time and a general claim size with continuous pdf.To do so, we start with a proposition from Chapter 3 of [27], which used the Picard method of successive approximation.

Proposition 1. Neumann Series for a Volterra IE of the 2nd Kind
For the Volterra IE of the 2nd kind, as in Equation (9), where g(T ) and K(T, s) are L 2 -functions, its Neumann series is given by: where Γ(T, s; λ) = ∞ n=1 λ n−1 K n (T, s) is the unique resolvent kernel and K n (T, s) is the n-th-iterated kernel function satisfying the recurrence formula: with K 1 (T, s) = K(T, s).
In order to prove our theorem, it is necessary to find the resolvent kernel, which is obtained in the following lemma.Lemma 2. Consider the kernel function given by Equation (12).For m = 1, 2, . .., its resolvent kernel is therefore given by: Γ(T, s; λ) = e −mδ(T −s) (15) Proof.Using Equation ( 14), we obtain K 2 (T, s), K 3 (T, s), for even n, the resolvent kernel is then obtained by summing up K m (T, s) as follows: Now, we can obtain the expression for the first and second moment, which is the main result of this article.
Theorem 3. The explicit solution of the first two moments are given by: where Proof.Applying Proposition 1 and Lemma 2 to Equation ( 3) with m = 1, 2, the results follow.
Section 4 will numerically illustrate the computation of the first and second moment under three copulas, assuming that the claim sizes are exponentially distributed, i.e., X ∼ exp(α).We do not proceed to obtain the closed form solution of the Neumann series expression for the higher moments, as it is tedious and time consuming.However, they are obtainable using the results provided in this section.

Numerical Illustration
We now present numerical illustration of the Neumann series expression for the first two moments.We start our discussion by presenting the scatter plots of each copula in Figures 1-3, where the marginals are exponential distribution, which is in line with the assumptions used in the numerical computations of the moments in this section.All computations were done using Mathematica.

Numerical Accuracy of Neumann Series Expression for Moments
Recall that Equation ( 16) has at most triple integration involved, while Equation ( 17) has up to sextuple numerical integration.This implies that the computation of Equation ( 17) is expected to be close to the solution by [6], due to numerical approximation error, and the values would vary according to selected software packages.To evaluate the performance of the main results, we compare the numerical values returned by our Neumann series (under the column Neumann of Table 1), using the FGM copula, with the numerical values given by the closed form solution in [6] (under the column BCLM of Table 1).
Our calculations showed that the Neumann series expression for the first moment gives the same value as the closed form solution presented in [6].On the other hand, the second moment gives a slightly different value at θ = −0.9 and 0.9, when the r.v.'s, X and W , are highly dependent.After a close scrutiny of the programming messages, we noticed that this is caused by numerical approximation errors of the quadruple, quintuple and sextuple integrations that are not present in the calculation of the first moment.To improve the accuracy of the Neumann series expression for higher order moments, the reader can use other software packages or use Monte Carlo simulation.

Moments of the Aggregate Discounted Claims
Setting δ = 0.04, α = 0.01 and β = 1 for the case of exponential claim inter-arrival time and exponential claim sizes, respectively, we show the values of moments of the aggregate discounted claims for each copula used in this study.We present the values of the first and second moments of the compound distribution, i.e., µ Z (5) and µ (2) Z (5), as well as the variance under each copula in Tables 2-4.The term "spread", which is defined as the difference between the values returned by θ at both ends, i.e., θ −0.95 − θ 0.95 for FGM and Gaussian, and θ 1 − θ 100 for Gumbel, are also shown in Tables 2 and 3. Our calculations showed that all copula exhibit decreasing values as θ increases, in line with [11].Intuitively, a negative dependence structure represented by the pair of short inter-claim waiting time (or frequent claim occurrence within a given time period) with huge claim size will only prompt the insurer to charge a higher premium as opposed to the positive dependence structure.
As we have expected, the values of moments do not vary much across θ when calculated under the FGM copula, as opposed to the Gaussian and Gumbel copulas.Being an extreme copula, values of the first moment calculated under the Gumbel copula also showed the widest spread of the first moment.
Table 5 shows the values of the first moment as a function of α and β, respectively, for which we use the Gaussian copula at θ = −0.9.It shows that increasing the inter-claim waiting time parameter, β, results in increasing the mean value of the aggregate discounted claims, and vice versa in the case of the claim size parameter.Given an average value of inter-claim arrival time, β, the mean of aggregate discounted claims gets lower as we have a lower average claim size, given by 1 α .On the other hand, given an average value of the claim size, the mean of the aggregate discounted claims gets bigger as the inter-claim arrival time gets shorter, which implies more frequent claim occurrences.This scenario is illustrated in Figure 4 for θ = −0.9(left hand side of the diagram), as well as θ = 0 (right hand side of the diagram).Table 5.Values of µ Z (5) under the Gaussian copula at θ = −0.9.

Premium Calculation under FGM, Gaussian and Gumbel copulas
We now compute the loaded premium related to the risk of an insurance portfolio represented by Z(T ), where the dependence structure is captured by a copula.For that purpose, the first two moments will be useful in the premium calculation based on the expected value principle, the variance principle, as well as the standard deviation (SD) premium principle, as the following: Table 6 exhibits the loaded premium according to the SD principle under the three copulas considered, with κ = 0.1, while Figures 5 and 6 illustrate the range of premiums under the copulas studied according to the SD premium principle.

Conclusions
In this paper, we utilized copulas to capture the dependence structure between the inter-claim arrival time and claim sizes in classical actuarial risk theory.To do so, we represented the expression for the m−th order moment proposed in [6,12] in the form of the Volterra integral equation (VIE) of the second kind, which is widely used in renewal theory, demographics, electromagnetism and viscoelasticity.We derived the Neumann series expression for this recursive equation using the Picard method of successive approximations, based on which we computed the first two moments of the aggregate discounted claims.For the dependence structure between the inter-claim arrival time and claim sizes, we used a Farlie-Gumbel-Morgenstern copula, a Gaussian copula and a Gumbel copula with exponential marginal distributions.We showed the values of moments of the aggregate discounted claims, as well as the loaded premium for each copula used in this study.
It would be of interest to derive the expression for Equations ( 2) and (3) using other joint pdfs between X and W .Other copulas with different claim size distributions for X may be considered in the proposed approach, which we leave for further research.We can also consider the Monte Carlo simulation, as well as other numerical methods to solve the VIE (such as Runge-Kutta and the collocation methods), as the next objective of further research to deal with the computation of higher moments.

Figure 4 .
Figure 4. Sensitivity of the first moment under the Gaussian copula at θ = 0 and θ = −0.9 with respect to claim size and inter-claim time averages.

Figure 5 .
Figure 5.The loaded premium under FGM and Gaussian copulas based on the SD premium principle.

Figure 6 .
Figure 6.The loaded premium under the Gumbel copula based on the SD premium principle.

Table 6 .
Loaded premium according to the SD principle under various copulas.