Multivariate Methods

Objectives: Multivariate distributions from shared mixing variables; 2 dimensional FFT techniques.

Audience: Users interested in directly estimating multivariate distributions.

Prerequisites: DecL, general use of aggregate, probability.

See also:

Contents:

• Helpful References

• Two Shortcomings

• Negative Multinomial Distribution and Related Frequency Distributions

• Correlated Aggregate Distributions

Helpful References

Two Shortcomings

aggregate methods are fundamentally univariate. Only one loss is tracked through each calculation. As a result, there are several programs it is impossible to model directly:

• The net of an occurrence cat program with an aggregate limit (need to track the net and cession to know when the aggregate limit is exhausted).

• The total cession from a specific and agg program (need to track two ceded results).

These are dual problems. It is possible to model them using two dimensional FFTs.

Negative Multinomial Distribution and Related Frequency Distributions

When we consider mixed Poisson distributions we often regard \(G\) as carrying inter-risk correlation, or more evocatively “contagion”, information about weather, the state of the economy and inflation, gas prices etc. Hence if we have two related frequency variables \(N_1\) and \(N_2\) we should expect to use the same \(G\) and produce a bivariate mixed Poisson where, conditional on \(G=g\), \(N_i\) has a Poisson distribution with mean \(n_i g\) and \(N_1\) and \(N_2\) are conditionally independent. The MGF of such a distribution will be

\[\begin{split}M(\zeta_1,\zeta_2) &= \mathsf{E}(e^{\zeta_1N_1+\zeta_2N_2}) \notag \\ &= \mathsf{E}(\mathsf{E}(e^{\zeta_1N_1+\zeta_2N_2}|G ) ) \notag \\ &= \mathsf{E}_G(\mathsf{E}(e^{\zeta_1N_1}|G ) \mathsf{E}(e^{\zeta_2N_2}|G ) ) \notag \\ &= \mathsf{E}_G(\exp(G(n_1(e^{\zeta_1}-1)+n_2(e^{\zeta_2}-1)))) \notag \\ &= M_G(n_1(e^{\zeta_1}-1) + n_2(e^{\zeta_2}-1)).\end{split}\]

For example, if \(G\) is a gamma random variable with MGF

\[M_G(\zeta) = (1-\beta \zeta)^{-k}\]

(mean \(k\beta\), variance \(k\beta^2\)) we get a bivariate frequency distribution with MGF

\[\begin{split}M(\zeta_1,\zeta_2) &= [1- \beta(n_1(e^{\zeta_1}-1) + n_2(e^{\zeta_2}-1))]^{-k} \notag \\ &=[1+\beta \sum_i n_i -\beta \sum_i n_ie^{\zeta_i}]^{-k} \notag \\ &=(Q -\sum_i P_ie^{\zeta_i})^{-k}\end{split}\]

where \(P_i=\beta n_i\), \(P=\sum_i P_i\) and \(Q=1+P\). Equation ([nmnmgf]) is the moment generating function for a negative multinomial distribution, as defined in Johnson, Kotz and Kemp. The negative multinomial distribution has positively correlated marginals as expected given its construction with a common contagion \(G\).

The form of the moment generating function for negative multinomial distribution can be generalized allowing us to construct multivariate frequency distributions \((N_1,\dots,N_t)\) where

1. Each \(N_i\) is a negative binomial.

2. The sum \(N_1+\cdots + N_t\) under the multivariate distribution is also negative binomial. (In general, the sum of independent negative binomials will not be negative binomial.)

3. The \(N_i\) are correlated.

We will call such multivariate frequencies, with common mixing distributions, \(G\)-mixed multivariate Poisson distributions.

Evolution of Claims Over Time

Here is an application of the NMN distribution. If \(A\) is an aggregate distribution representing ultimate losses we may want to determine a decomposition \(A=\sum_t D_t\) into a sum of losses paid at time \(t\) for \(t=1,\dots,T\).

If \(A=X_1+\cdots +X_N\) has a compound Poisson distribution then such a decomposition is easy to arrange. Let \(\pi_t\) be the expected proportion of ultimate losses paid at time \(t\), so \(\sum_{t=1}^{t=T} \pi_t=1\). By definition we mean

\[\mathsf{E}(D_t) = \pi_t \mathsf{E}(A).\label{meandt}\]

(Equation ([meandt]) is a different assumption to

\[\mathsf{E}(D_{t})=\pi_t\mathsf{E}(A|\text{information available at $t-1$})=\pi_t A_{t-1},\]

which is closer to the problem actually faced by the reserving actuary. Our \(\pi_t\)’s are prior estimates assumed known at time 0. These types of differences have interesting implications for actuarial methods and they are explored further in Mack.) Now we seek a decomposition

\[A=D_1+D_2+\cdots+D_T\label{decomp}\]

but we know only ([meandt]). The simplest approach is to assume that severity \(X\) is independent of time and that \(\pi_t n\) of the total \(n\) claims are paid at time \(t\). If we further assume that the number of claims paid at time \(t\) is also Poisson, then the moment generating function of \(D_1+\cdots+D_T\) is given by

\[\begin{split}M_{D_1+\cdots+D_T}(\zeta) &= \prod_t \exp( \pi_t n(M_X(\zeta)-1)) \notag \\ & = \exp(n(\sum_t \pi_t M_X(\zeta)-1)) \notag \\ &= \exp(n(M_X(\zeta)-1)) \notag \\ & = M_A(\zeta).\end{split}\]

Thus we have a very simple decomposition for ([decomp]): the individual \(D_t\) are independent compound Poisson variables with expected claim count \(\pi_t n\) and severity distribution \(X\).

Moving one step further, it is often observed in practice that average severity increases with \(t\) so the assumption that \(X\) is fixed for all \(t\) is unrealistic. It may be better to assume that losses which close at time \(t\) are samples of a random variable \(X_t\). As above, we assume that the expected number of such losses is \(\pi_t' n\) where \(n\) is the expected ultimate number of claims, and \(\pi_t'\) adjusts the original \(\pi_t\) for the difference in average severity \(\mathsf{E}(X)\) vs. \(\mathsf{E}(X_t)\). Now

\[\begin{split}M_{D_1+\cdots+D_T}(\zeta) &= \prod_t \exp( \pi_t' n(M_{X_t}(\zeta)-1)) \notag \\ & = \exp(n(\sum_t \pi_t' M_{X_t}(\zeta)-1)) \notag \\ &= \exp(n(M_{X'}(\zeta)-1)) \notag \\ & = M_A(\zeta)\end{split}\]

where \(X'\) is a mixture of the \(X_t\) with weights \(\pi_t'\). Equation ([CPDecomp2]) is a standard result in actuarial science, see Bowers et al.

If we try to replicate the compound Poisson argument using a negative binomial distribution for \(N\) we will clearly fail. However if \(X\) is defined as a mixture of \(X_t\) with weights \(\pi_t\), as before, then we can write

\[M_{D_1,\dots,D_T}(\zeta_1,\dots,\zeta_T) = (Q-\sum_t P\pi_t M_{X_t}(\zeta_t))^{-k}\label{nmn1}\]

and so

\[M_A(\zeta) =M_{D_1,\dots,D_T}(\zeta,\dots,\zeta) = (Q-\sum_t P_t M_{X_t}(\zeta))^{-k}=(Q-PM_X(\zeta))^{-k} \label{nmn2}\]

where \(P_t:=\pi_t P\). Equation ([nmn1]) is the MGF for a negative multinomial distribution, as defined in the previous section and Johnson, Kotz and Kemp. As we have seen the negative multinomial distribution has positively correlated marginals, in line with our prior notions of liability dynamics. It therefore provides a good model for the decomposition of ultimate losses into losses paid each period.

Related Multivariate Frequency Distributions

We can use the same trick with other mixing distributions than the gamma. The Poisson inverse Gaussian (PIG) distribution is an inverse Gaussian mixture of Poissons, just as the negative binomial distribution is a gamma mixture. The MGF is

\[M(\zeta) = \exp(-\tau (\sqrt{1+\beta(e^{\zeta}-1)}-1)). \label{pig-pgf}\]

The mean is \(\tau\beta\) and the variance is \(\tau\beta(1+\beta)\). We can define a multivariate PIG (MPIG) by

\[M(\zeta_1,\dots,\zeta_T) = \exp(-\tau (\sqrt{1+\sum\beta_i(e^{\zeta_i}-1)}-1)). \label{mpig-pgf}\]

Sichel’s distribution is an generalized-inverse Gaussian mixture of Poissons. The MGF is

\[M(\zeta) = \frac{K_{\gamma}(\omega\sqrt{1-2\beta(e^{\zeta}-1)})}{ K_{\gamma}(\omega)(1-2\beta(e^{\zeta}-1))^{\gamma/2}}.\]

The mean and variance are given in Johnson, Kotz and Kemp [page 456]. Clearly we can apply the same techniques to get another multivariate frequency distribution.

The Poisson-Pascal distribution is a Poisson-stopped sum of negative binomials. It has moment generating function

\[M(\zeta) = \exp(\theta ((1-P(e^{\zeta}-1))^{-k}-1))\]

and so will also yield another multivariate family. The mean and variance are given by

\[\begin{split}\mu = \theta kP \\ \mu_2 = \theta kP(Q+kP).\end{split}\]

Excess count interpretation of \(G\)-mixed multivariate Poisson distributions

The reader has probably realized that a \(G\)-mixed multivariate Poisson seems closely related to a single \(G\)-mixed Poisson and a series of indicator variables, combining results from the previous sub-sections with Section 1.6. Let \(N\) be \(G\)-mixed Poisson with parameter \(n\) and \(\mathsf{var}(G)=c\). Let \((N_1,N_2)\) be \(G\)-mixed bivariate Poisson with parameters \(n_1\) and \(n_2\) and the same \(G\), so the MGF of \((N_1,N_2)\) is

\[M_1(\zeta_1,\zeta_2)=M_G(n_1(e^{\zeta_1}-1)+n_2(e^{\zeta_2}-1)).\]

Finally let \((I,J)\) be a bivariate distribution supported on \(\{0,1\}\times\{0,1\}\) with

\[\begin{split}\mathsf{Pr}(I=0,J=0) &= p_{00} \\ \mathsf{Pr}(I=1,J=0) &= p_{10} \\ \mathsf{Pr}(I=0,J=1) &= p_{01} \\ \mathsf{Pr}(I=1,J=1) &= p_{11}\end{split}\]

and \(\sum p_{ij}=1\).

We can define a new bivariate distribution from \((I,J)\) and \(N\) as

\[(M_1,M_2)=(I_1,J_1)+\cdots +(I_N,J_N).\]

The MGF of \((M_1,M_2)\) is

\[M_2(\zeta_1,\zeta_2)=M_G(n(p_{11}e^{\zeta_1+\zeta_2} + p_{10}e^{\zeta_1} + p_{01}e^{\zeta_2} + p_{00}).\]

Thus, if \(p_{11}=0\) we see the single-frequency sum of the bivariate \((M_1,M_2)\) is actually a \(G\)-mixed bivariate Poisson. If \(p_{00}=0\) then \(n=n_1+n_2\), otherwise \((1-p_{00})n=n_1+n_2\) and there are some extraneous “zero” claims. However, if \(p_{11}\not=0\) then the single frequency sum is not a \(G\)-mixed bivariate Poisson.

Here is an interesting interpretation and application of \((I,J)\). We can regard \(I\) as an indicator of whether a claim has been reported at time \(t\) and \(J\) and indicator of whether the claim is closed. Then

\[\begin{split}\mathsf{Pr}(I=0,J=0) &=\text{meaningless} \\ \mathsf{Pr}(I=1,J=0) &=\text{reported claim which closes without payment} \\ \mathsf{Pr}(I=0,J=1) &=\text{claim not yet reported which closes with payment} \\ \mathsf{Pr}(I=1,J=1) &=\text{claim reported and closed with payment}.\end{split}\]

Combining with a distribution \(N\) of ultimate claims we can use ([singleFreq]) to produce \((M_1,M_2)=(I_1+\cdots+I_N,J_1+\cdots+J_N)\)—a bivariate distribution of (claims reported at time \(t\), ultimate number of claims)! Note the value \((0,0)\) is a meaningless annoyance (it scales \(n\)) and we assume \(p_{00}=0\). The three other parameters can easily be estimated using standard actuarial methods.

Given such a bivariate and a known number of claims reported we can produce a posterior distribution of ultimate claims. Furthermore, in all these techniques we can extend the simple count indicators \((I,J)\) to be the distribution of case incurred losses and ultimate losses. Then we would get a bivariate distribution of case incurred to date and ultimate losses. I believe there is a lot of useful information that could be wrought from these methods and that they deserve further study. They naturually give confidence intervals on reserve ranges, for example.

We end with a numerical example illustrating the theory we have developed and showing another possible application. Rather than interpreting \(p_{ij}\) as reported and ultimate claims we could interpret them as claims from line A and line B, where there is some expectation these claim would be correlated. For example A could be auto liability and B workers compensation for a trucking insured. Let \(c=0.02\) be the common contagion and \(n=250\). Then let

\[\begin{split}\mathsf{Pr}(I=0,J=0) &= 0 \\ \mathsf{Pr}(I=1,J=0) &= 0.45 \\ \mathsf{Pr}(I=0,J=1) &= 0.05 \\ \mathsf{Pr}(I=1,J=1) &= 0.50.\end{split}\]

We interpret \(I\) as indicating a workers compensation claim and \(J\) as indicating an auto liability claim. The distribution says that when there is an auto liability claim (\(J=1\)) there is almost always an injury to the driver, resulting in a workers compensation claim (\(I=1\)). However, there are many situations where the driver is injured but there is no liability claim—such as back injuries. Overall we expect \(250(0.45+0.50)=237.5\) workers compensation claims and \(250(0.05+0.5)=137.5\) auto liability claims and 250 occurrences.

We will consider the single-frequency bivariate distribution and the negative multinomial. We have seen that the negative multinomial distribution will be slightly different because \(p_{11}\not=0\). The appropriate parameters are \(n_1=250(p_{10}+p_{11})=237.5\) and \(n_1=250(p_{01}+p_{11})=137.5\). Figure 1.1 shows the negative multinomial bivariate (top plot) and the single-frequency bivariate aggregate of \((I,J)\) (bottom plot). Because of the correlation between \(I\) and \(J\), \(p_{11}=0.5\), the lower plot shows more correlation in aggregates and the conditional distributions have less dispersion. Figure 1.2 shows the two marginal distributions, which are negative binomial \(c=0.02\) and mean 237.5 and 137.5 respectively, the sum of these two variables assuming they are independent (labelled “independent sum”), the sum assuming the negative multinomial joint distribution (“NMN Sum”) which is identical to a negative binomial with \(c=0.02\) and mean \(350=237.5+137.5\), the total number of claims from both lines, and finally, the sum with dependent \((I,J)\) (“bivariate sum”). The last sum is not the same as the negative binomial sum; it has a different MGF.

Figure 1.2 also shows the difference between the sum of two independent negative binomials with means \(n_1\) and \(n_2\) and contagion \(c\) and a negative binomial with mean \(n_1+n_2\) and contagion \(c\). The difference is clearly very material in the tails and is an object lesson to modelers who subdivide their book into homogeneous parts but then add up those parts assuming independence. Such an approach is wrong and must be avoided.

As the contagion \(c\) increases the effects of \(G\)-mixing dominate and the difference between the two bivariate distributions decreases, and conversely as \(c\) decreases to zero the effect is magnified. The value \(c=0.02\) was selected to balance these two effects.

Comparison of negative multinomial (top) and single frequency bivariate claim count (bottom) bivariate distributions.

Comparison of negative multinomial and single frequency bivariate claim count marginal and total distributions.

Correlated Aggregate Distributions

Here we extend some of the ideas in Section 1.7.3 from plain frequency distributions to aggregate distributions. Begin with bivariate aggregate distributions. There are two different situations which commonly arise. First we could model a bivariate severity distribution and a univariate count distribution:

\[(A,B)=(X_1,Y_1)+\cdots+(X_N, Y_N).\]

Equation ([modelone]) arises naturally as the distribution of losses and allocated expense, ceded and retained losses, reported and ultimate claims, and in many other situations. Secondly we could model

\[(A,B)=(X_1+\cdots+X_M, Y_1+\cdots+Y_N)\]

where \(X_i\) and \(Y_j\) are independent severities and \((M,N)\) is a bivariate frequency distribution. ([modeltwo]) could be used to model losses in a clash policy.

We will use the following notation. \(A=X_1+\cdots+X_M\) and \(B=Y_1+\cdots+Y_N\) are two aggregate distributions, with \(X_i\) iid and \(Y_j\) iid, but neither \(X\) and \(Y\) nor \(M\) and \(N\) necessarily independent. Let \(\mathsf{E}(X)=x\) and \(\mathsf{E}(Y)=y\), \(\mathsf{var}(X)=v_x\) and \(\mathsf{var}(Y)=v_y\). Let \(\mathsf{E}(M)=m\), \(\mathsf{E}(N)=n\), \(c\) be the contagion of \(M\) and \(d\) that of \(N\). Hence \(\mathsf{var}(M)=m(1+cm)\) and \(\mathsf{var}(N)=n(1+dn)\).

Will now calculate the correlation coefficient between \(A\) and \(B\) in four situations.

Correlated Severities, Single Frequency

Assume that the bivariate severity distribution \((X,Y)\) has moment generating function \(M_{(X,Y)}(\zeta,\tau)\). Also assume that the claim count distribution \(N\) is a \(G\)-mixed Poisson. Then, just as for univariate aggregate distributions, the MGF of the bivariate aggregate \((A,B)\) is

\[M_{(A,B)}(\zeta,\tau)= M_G( n(M_{(X,Y)}(\zeta,\tau)-1)).\label{bivMGF}\]

Therefore, since \(\mathsf{E}(G)=1\) and \(\mathsf{E}(G^2)=1+c\),

\[\begin{split}E(AB) &= \frac{\partial^2 M_{(A,B)}}{\partial\zeta\partial\tau}\Big\vert_{(0,0)} \notag \\ &= M_G''(0)n^2\frac{\partial M_{(X,Y)}}{\partial\zeta} \frac{\partial M_{(X,Y)}}{\partial\zeta} + M_G'(0)n \frac{\partial^2 M_{(X,Y)}}{\partial\zeta\partial\tau} \notag \\ &=(1+c)n^2xy + n\mathsf{E}(XY) \notag \\ &=(1+c)n^2xy + n\mathsf{cov}(X,Y) + nxy.\end{split}\]

The value of \(\mathsf{cov}(X,Y)\) will depend on the particular bivarate severity distribution.

For example, suppose that \(Z\) represents ground up losses, \(X\) represents a retention to \(a\) and \(Y\) losses excess of \(a\) (per ground up claim), so \(Z=X+Y\). Then \((X,Y)\) is a bivariate severity distribution. Since \(Y\) is zero when \(Z\le a\) we have \(\mathsf{cov}(X,Y)=(a-x)y\).

Bivariate Frequency

The second method for generating correlated aggregate distributions is to use a bivariate frequency distribution. So, suppose \((M,N)\) has a \(G\)-mixed bivariate Poisson distribution. The variance of \(A\) is given by Equation ([varAgg]). To compute the covariance of \(A\) and \(B\) write the bivariate MGF of \((A,B)\) as

\[M_{(A,B)}(\zeta,\eta)=M(\zeta,\eta)=M_G(m(M_X(\zeta)-1) +n(M_Y(\eta)-1))=M_G(\psi(\zeta,\eta))\]

where the last equality defines \(\psi\). Then, evaluating at the partial derivatives at zero, we get

\[\begin{split}\mathsf{E}(AB) &= \frac{\partial^2 M}{\partial\zeta\partial\eta} \notag \\ &= \frac{\partial^2 M_G}{\partial t^2} \frac{\partial \psi}{\partial\zeta} \frac{\partial \psi}{\partial\eta} + \frac{\partial M_G}{\partial t} \frac{\partial^2 \psi}{\partial\zeta\partial\eta} \notag \\ &= (1+c)mxny.\end{split}\]

Hence

\[\mathsf{cov}(A,B)=\mathsf{E}(AB)-\mathsf{E}(A)\mathsf{E}(B)=cmnxy.\]

Parameter Uncertainty

It is common for actuaries to work with point estimates as though they are certain. In reality there is a range around any point estimate. We now work through one possible implication of such parameter uncertainty. We will model \(\mathsf{E}[A]=R\) and \(\mathsf{E}[B]=S\) with \(R\) and \(S\) correlated random variables, and \(A\) and \(B\) conditionally independent given \(R\) and \(S\). We will assume for simplicity that the severities \(X\) and \(Y\) are fixed and that the uncertainty all comes from claim counts. The reader can extend the model to varying severities as an exercise. \(R\) and \(S\) pick up uncertainty in items like the trend factor, tail factors and other economic variables, as well as the natural correlation induced through actuarial methods such as the Bornheutter-Ferguson.

Suppose \(\mathsf{E}[R]=r\), \(\mathsf{E}[S]=s\), \(\mathsf{var}(R)=v_r\), \(\mathsf{var}(S)=v_s\) and let \(\rho\) be the correlation coefficient between \(R\) and \(S\).

By ([varAgg]) the conditional distribution of \(A \mid R\) is a mixed compound Poisson distribution with expected claim count \(R/x\) and contagion \(c\). Therefore the conditional variance is

\[\begin{split}\mathsf{var}(A \mid R) &= \mathsf{E}[M \mid R]\mathsf{var}(X)+\mathsf{var}(M \mid R)\mathsf{E}[X]^2 \\ &= R/x v_x + R/x(1+cR/x) x^2 \\ &= xR(1+ v_x/x^2) + cR^2,\end{split}\]

and the unconditional variance of \(A\) is

\[\begin{split}\mathsf{var}(A) &= \mathsf{E}[\mathsf{var}(A \mid R)] + \mathsf{var}(\mathsf{E}[A \mid R]) \\ &= \mathsf{E}[xR(v_x/x^2+1)+cR^2] + \mathsf{var}(R) \\ &= xr(v_x/x^2+1)+c(v_r+r^2) + v_r.\end{split}\]

Next, because \(A\) and \(B\) are conditionally independent given \(R\) and \(S\),

\[\begin{split}\mathsf{cov}(A,B) &= \mathsf{E}[\mathsf{cov}(A,B \mid R,S)] + \mathsf{cov}(\mathsf{E}[A \mid R], \mathsf{E}[B \mid S]) \\ &= \mathsf{cov}(R, S).\label{simpleCov}\end{split}\]

Note Equation ([simpleCov]) is only true if we assume \(A\not=B\).

Parameter Uncertainty and Bivariate Frequency

Finally, suppose \(\mathsf{E}[A]=R\), \(\mathsf{E}[B]=S\) with \(R\) and \(S\) correlated parameters and conditional on \((R,S)\) suppose that \((M,N)\) has a \(G\)-mixed bivariate Poisson distribution. By ([covMNM]) \(\mathsf{cov}(A,B \mid R,S)=cRS\). The unconditional variances are as given in ([varA]). The covariance term is

\[\begin{split}\mathsf{cov}(A,B) &= \mathsf{E}[\mathsf{cov}(A,B \mid R,S)] + \mathsf{cov}(\mathsf{E}[A \mid R], \mathsf{E}[B \mid S]) \\ &= c\mathsf{E}[RS] + \mathsf{cov}(R,S) \\ &= (1+c)\mathsf{cov}(R,S) + crs \\ &= \rho \sqrt{v_rv_s}(1+c)+crs.\end{split}\]