2.1 Unbalanced feature of reserving data

As described in section 1, loss reserving data typically exhibit unbalanced nature. We consider for illustration a real data set from a Canadian insurer collected from 2003 to 2012 (denoted by years 1–10). This data set is used for illustration in Côté et al. (Reference Côté, Genest and Abdallah2016) and is provided in Tables C.1 and C.2 in Appendix C. The two lines of business (LOBs) used for illustration are Bodily Injury line and Accident Benefits (excluding Disability Income).

Figure 1 provides heat maps of incremental loss ratios on the left and a plot of incremental loss ratios for accident year 2003 from the two LOBs on the right. For a given accident year, the loss ratio increment for development year j is defined as the ratio of incremental claim payments in that development year to the earned premium for the accident year. Within a single loss triangle, one can observe a quite significant variation in claim observations across development years for any particular accident year. As shown in the heat map in Figure 1, the claim activity for the Bodily Injury line is low in development year 0, then peaks in the next few years and dies out after the peak. This typical pattern is also shown in the plot of loss ratio on the right-hand side of Figure 1 for accident year 2003. For the Accident Benefits line, the claim activity is the highest in development year 0 or 1, then drops quickly as we approach later development year. This pattern is also shown in the plot of loss ratios for accident year 2003 of the Accident Benefit line. The plot of ratios on the right-hand side of Figure 1 also indicates the difference in development patterns for two different business lines. We can say that the Accident Benefits line is shorter-tailed than the Bodily Injury line. A variation can be observed across claim observations that come from the same accident year and the same development year, simply due to different claim development patterns across these business lines. This is in addition to the variation between loss values in different development lags and from different loss triangles, such as cells in development year 1 from the Bodily Injury line and cells in development year 10 from the Accident Benefit line.

Figure 1 Loss ratios from Bodily Injury line and Accident Benefits line (from a Canadian insurer). (a) Heat maps of Bodily Injury line (top) and Accident Benefit line (bottom). (b) Plot of loss ratios for accident year 2003.

Overall, Figure 1 shows a large variation across claim observations in a loss reserving data set. Within a single loss triangle, there is variation across development years due to the development pattern of claims over time. Different claim development patterns can also result in variation between observations across loss triangles. Typically, one often does not expect dependence between lines that have different tail lengths, for example, an Auto Property Damage line is often independent of a Workers Compensation line. However, lines with different tail lengths can still have some association. One of such examples is a portfolio of Accident Benefits line and the Bodily Injury line in the above illustration.

With the variations between claim observations within and across triangles, we refer to loss reserving data as unbalanced data. This data feature creates a number of challenges in applying common shock models to reserving data, which will be discussed in the remainder of this section. For generality and completeness, the focus is placed on the unbalanced feature of data consisting of multiple LOBs.

2.2 General common shock framework

Consider N loss triangles of claim cells $Y_{i,j}^{(n)}$ . The notation $Y_{i,j}^{(n)}$ can represent incremental claim payments or counts. We have the indices $i\, (i=1,...,I)$ representing the accident period, $j\,(j=0,...,J)$ representing the development period and $n\,(n=1,...,N)$ representing the business line. It also follows that the claims $Y_{i,j}^{(n)}$ belong to the calendar period $t=i+j-1,\,(t=1,...,T)$ .

Let $\mathcal{S}^{(n)} = \lbrace\mathcal{S}_s^{(n)};\,s=1,...,S\rbrace$ be a partition of the set of all claims $Y_{i,j}^{(n)}$ from business line n. Also assume that all partitions are the same for different lines n for simplicity. Denote by $\boldsymbol{\pi}_{(i,j)}={s}$ a unique mapping of claim $Y_{i,j}^{(n)}$ to a set $\mathcal{S}_s^{(n)}$ in the partition. For example, the partition $\mathcal{S}^{(n)} = \lbrace\mathcal{S}_s^{(n)};\,s=1,...,I\rbrace$ , where $\mathcal{S}_s^{(n)} = \lbrace Y_{s,j}^{(n)};\, j=1,...,J\rbrace$ represents a partition of claims by accident year. The selection of the partition $\mathcal{S}^{(n)}$ is very flexible and can be specified for different types of dependence.

Many multivariate models with different types of dependence can be generalised by the common shock framework in Avanzi et al. (Reference Avanzi, Taylor and Wong2018) with

(1) \begin{equation} Y_{i,j}^{(n)} = \kappa_{i,j}^{(n)}W_{\boldsymbol{\pi}_{(i,j)}} + \lambda_{i,j}^{(n)}U_{\boldsymbol{\pi}_{(i,j)}}^{(n)} + Z_{i,j}^{(n)} \end{equation}

where $\boldsymbol{\pi}_{(i,j)}=s$ denote the unique mapping of the claim $Y_{i,j}^{(n)}$ to the corresponding subset $\mathcal{S}_s^{(n)}$ in the partition $\mathcal{S}^{(n)}$ , and $W_{\boldsymbol{\pi}_{(i,j)}},\,U_{\boldsymbol{\pi}_{(i,j)}}^{(n)}$ , and $Z_{i,j}^{(n)}$ are independent stochastic variates. The common shock $W_{\boldsymbol{\pi}_{(i,j)}}$ introduces dependence across all business lines $n=,1...,N$ on claims that belong to the subsets $\mathcal{S}_s^{(n)}$ . For example, accident year dependence across lines can be captured using the partition where $\mathcal{S}_s^{(n)} = \lbrace Y_{s,j}^{(n)};\, j=1,...,J\rbrace$ . The other common shock $U_{\boldsymbol{\pi}_{(i,j)}}^{(n)}$ introduces dependence across claims within the set $\mathcal{S}_s^{(n)}$ of business line n only, such as development year dependence with the partition set specification $\mathcal{S}_s^{(n)} =\lbrace Y_{i,s}^{(n)};\, i=1,...,I\rbrace$ . Overall, the flexibility of choice of the subsets $\mathcal{S}_s^{(n)}$ allows different dependence structures to be captured. Lastly, the idiosyncratic component, which is unique to the claim $Y_{i,j}^{(n)}$ , is denoted by $Z_{i,j}^{(n)}$ . Scaling factors, denoted by $\kappa_{i,j}^{(n)},\,{\lambda}_{i,j}^{(n)}$ , control the extent to which the set-wide common shock contributes to individual members of the set. In this section, we have wished to preserve the link to the general notation of Avanzi et al. (Reference Avanzi, Taylor and Wong2018) through the use of the notation $\boldsymbol{\pi}_{(i,j)}$ . This notation will be simplified in section 2.3 for specific examples.

2.3 Balancing common shock proportions in loss reserving data

As a result of the unbalanced feature of reserving data, a common shock model can create problems in the absence of careful modelling.

Consider a special case of equation (1) for dependence within a business line (i.e. $W_{\boldsymbol{\pi}_{(i,j)}}=0$ ). Further specify accident period dependence (i.e. $\boldsymbol{\pi}_{(i,j)}^{(n)} = p$ for the mapping of subsets in the partition where $\mathcal{S}_s^{(n)} = \lbrace Y_{s,j}^{(n)};\, j=1,...,J\rbrace$ ). This allows us to simplify $U_{\boldsymbol{\pi}_{(i,j)}}^{(n)}={X_{i}}^{(n)}$ . Hence, the general framework is reduced to

(2) \begin{align} Y_{i,j}^{(n)} = \lambda_{i,j}^{(n)}{X_{i}^{(n)}} + Z_{i,j}^{(n)} \end{align}

Consequently, the proportionate contribution of the common shock to the expected value of the total observation is

(3) \begin{align} \frac{\lambda_{i,j}^{(n)}\text{E}\left[{X_{i}^{(n)}}\right]}{ \lambda_{i,j}^{(n)}\text{E}\left[{X_{i}^{(n)}}\right] + \text{E}\left[Z_{i,j}^{(n)}\right]} \end{align}

If the scaling factor is removed, that is, $\lambda_{i,j}^{(n)}=1$ , this proportion has an inverse relationship with the mean of the idiosyncratic component $\text{E}\left[Z_{i,j}^{(n)}\right]$ . As a result, in a set of loss cells in a triangle that are dependent and share a common shock, the cells with large values have a smaller proportion of common shock contribution and vice versa. This is because claims within the same accident period, or within the same calendar period belong to different development periods. As explained in section 2.1, their values can vary significantly due to the variation in claim activity across development periods. This issue can also be observed in the case of calendar period dependence (i.e. $\boldsymbol{\pi}_{(i,j)}^{(n)} = s$ for the mapping of subsets in the partition where $\mathcal{S}_s^{(n)} = \lbrace Y_{i,s-i+1}^{(n)};\, i=1,...,J\rbrace$ ).

A similar issue is encountered for a portfolio of dependent business lines with differing tail lengths, such as the two business lines Bodily Injury and Accident Benefits in the illustration in Section 2.1. We consider a special case of equation (1) that allows for dependence between business lines only (i.e. ${U_{\boldsymbol{\pi}_{i,j}}^{(n)}=0}$ ). Further specify cellwise dependence (i.e. partition mapping where $\mathcal{S}_{i,j}^{(n)} = \lbrace Y_{i,j}^{(n)}\rbrace$ ). This allows us to simplify $W_{\boldsymbol{\pi}_{(i,j)}}={V_{i,j}}$ . The contribution of the common shock to the total expected observation is then given by

(4) \begin{align}\frac{\kappa_{i,j}^{(n)}\text{E}\left[V_{i,j}\right]}{ \kappa_{i,j}^{(n)}\text{E}\left[V_{i,j}\right] + \text{E}\left[Z_{i,j}^{(n)}\right]}\end{align}

If the scaling factor is removed, that is, $\kappa_{i,j}^{(n)}$ , this proportion also has an inverse relationship with the mean of the idiosyncratic component $\text{E}\left[Z_{i,j}^{(n)}\right]$ . As explained in section 2.1, values of claims in a portfolio of multiple triangles can vary in two main ways: across development years within a loss triangle, and across loss triangles. As a result, the proportion of common shock varies within and across loss triangles, wherein loss cells with larger values have smaller common shock contributions. In the case of pairwise dependence considered above, the disproportion is typically a result of varying tail lengths across business lines. However, it is worth noting that unbalanced common shock proportions can also be typically observed for accident year dependence, or calendar dependence across business lines from the same cause.

We consider the case of accident year dependence across the two triangles illustrated in section 2.1 (i.e. partition mapping $\mathcal{S}_s = \lbrace Y_{s,j}^{(n)};\, j=1,...,J; n=1,...,N\rbrace$ , and we can simplify $W_{\boldsymbol{\pi}_{(i,j)}}={V_{i}}$ ). For illustration, the mean of the common shock $E\left[{V_{i}}\right]$ is set to 5% of the loss ratios in the first development year of each accident year in the Bodily Injury line. The contributions of common shock are shown in Figure 2 assuming no scaling terms. With accident year dependence across business lines, claims within the same accident year share the same common shock. These include claims from different development years within and across loss triangles. Because their values vary due to different claim activities within and across lines, their common shock proportions also vary. Specifically, common shock proportions are significantly smaller in areas with high claim activity, and larger in areas with low claim activity, as shown in Figure 2.

Figure 2 Heat maps of common shock contributions in Bodily Injury line (top) and Accident Benefit line (bottom) without using scaling terms.

In general, quite significant variations in common shock proportions can be observed within and across segments in the absence of careful modelling as a result of the unbalanced nature of loss reserving data. One may wish to confine the relation of the common shock to total observations over the entire range of the triangles.

2.4 Maintaining model parsimony

The most straightforward solution to the balancing common shock proportions within and across triangles is to have cell-specific scaling factors $\kappa_{i,j}^{(n)},\,\lambda_{i,j}^{(n)}$ to adjust the common shock effects for each total observation $Y_{i,j}^{(n)}$ . However, this implies that 2IJN new parameters are required for the entire range of triangles of observed data and outstanding claims to be predicted. Given that the variation in claim observations typically occurs across development periods, one may simplify the scaling factors to be column-specific $\kappa_{i,j}^{(n)}=\kappa_{j}^{(n)},\,\lambda_{i,j}^{(n)}=\lambda_{j}^{(n)}$ . However, this still results in 2JN new parameters.

Loss triangle data typically have a small sample size. While the presence of scaling factors can mitigate the impact of the unbalanced nature of reserving data, it also adds many more parameters to the model. If scaling factors are not chosen carefully, it may result in overfitting and the number of parameters to be estimated can even exceed the number of observations.

2.5 Maintaining distributional tractability

On some occasions, parameters $\lambda_{i,j}^{(n)}$ and $\kappa_{i,j}^{(n)}$ need to be specified such that the total observation $Y_{i,j}^{(n)}$ follows a specific distribution (Avanzi et al., Reference Avanzi, Taylor and Wong2018). This is referred to as distributional tractability, or closure under the taking of marginals, which is considered in Joe (Reference Joe1997, Chapter 4) to be one of the four desirable properties of a multivariate model (see also section 1).

Consider as an example the common shock Tweedie framework in Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b). This framework is developed for cellwise dependence across business lines (i.e. $\mathcal{S}_{i,j}^{(n)} = \lbrace Y_{i,j}^{(n)}\rbrace$ , $U_{i,j}^{(n)}=0$ ). Fitting this into the general common shock structure in equation (1) and simplifying $W_{\boldsymbol{\pi}_{(i,j)}}={V_{i,j}}$ , we have

(5) \begin{align} Y_{i,j}^{(n)} & = \kappa_{i,j}^{(n)}{V_{i,j}}+Z_{i,j}^{(n)} \end{align}

where the two components ${V_{i,j}},\,Z_{i,j}^{(n)}$ are assumed to be independent and have Tweedie distributions

(6) \begin{align} {V_{i,j}} & \sim \text{Tweedie}_p(\alpha,\beta) \end{align}

(7) \begin{align} Z_{i,j}^{(n)} & \sim \text{Tweedie}_p(\eta_i^{(n)}\nu_j^{(n)}, \gamma^{(n)}) \\[8pt]\nonumber \end{align}

Parameter p is the power parameter which specifies a member of the Tweedie family, for example, $p=1$ corresponds to a Poisson distribution. The representation of Tweedie distributions used is the reproductive representation (Jørgensen, Reference Jørgensen1997, Chapter 4). This representation specifies a Tweedie random variable using a location (or mean) parameter and a dispersion parameter. In the above model specification, parameters $\alpha$ and $\eta_i^{(n)}\nu_j^{(n)}$ are the location parameters, and parameters $\beta$ and $\gamma^{(n)}$ are the dispersion parameters. The reproductive representation has a distinctive property, wherein the weighted average of independent Tweedie variables with the same power parameter p and the same location parameter is also a Tweedie variable with the same power and location parameters. The weighting factors are determined using dispersion parameters of the component variables in the weighted average.

It then follows that the mean and variance of the two components ${V_{i,j}},\,Z_{i,j}^{(n)}$ are

(8) \begin{align} \text{E}\left[{V_{i,j}}\right] &= \alpha, & \qquad \text{Var}\left[{V_{i,j}}\right] &= \beta\alpha^p\end{align}

(9) \begin{align} \text{E}\left[Z_{i,j}^{(n)}\right] &= \eta_i^{(n)}\nu_j^{(n)}, & \qquad \text{Var}\left[Z_{i,j}^{(n)}\right] &= \gamma^{(n)}(\eta_i^{(n)}\nu_j^{(n)})^p \\[8pt]\nonumber \end{align}

As stated in Remark 2.2 of Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b), the most simple parametrisation is used for the common shock component ${V_{i,j}}$ with parameters $\alpha$ and $\beta$ .

As mentioned earlier in this section, it can be desirable to maintain distributional tractability, or closure under the taking of marginals for ease of interpretation. It follows from the form of closure under addition of the Tweedie family of distributions, as proven in Jørgensen (1997, Chapter 3), that a specific choice of $\kappa_{i,j}^{(n)}$ is required to ensure that $Y_{i,j}^{(n)}$ also has a Tweedie distribution. This choice is

(10) \begin{align} \kappa_{i,j}^{(n)} = \left(\dfrac{\alpha}{ \eta_i^{(n)}\nu_j^{(n)} }\right)^{1-p}\dfrac{\gamma^{(n)}}{\beta} \end{align}

The mean expression is given by

(11) \begin{equation} \text{E}\left[Y_{i,j}^{(n)}\right] = \left(\dfrac{\alpha}{\eta_i^{(n)}\nu_j^{(n)}}\right)^{2-p}\dfrac{\gamma^{(n)}}{\beta} \eta_i^{(n)}\nu_j^{(n)} +\eta_i^{(n)}\nu_j^{(n)} \end{equation}

where the first term in the summation is the contribution from the common shock and the second term is the contribution from the idiosyncratic component. The expected contribution of the common shock to the total expected observation is

(12) \begin{equation} \dfrac{\left(\dfrac{\alpha}{\eta_i^{(n)}\nu_j^{(n)}}\right)^{2-p}\dfrac{\gamma^{(n)}}{\beta} }{\left(\dfrac{\alpha}{\eta_i^{(n)}\nu_j^{(n)}}\right)^{2-p}\dfrac{\gamma^{(n)}}{\beta} +1} \end{equation}

The following observation can be made on the effect of the power parameter p:

1. – If $p<2$ : The above ratio increases as $\nu_j^{(n)}$ decreases. As a result, the proportion of common shock is understated in early development periods and overstated in late development periods (Avanzi et al., Reference Avanzi, Taylor and Wong2018). In a portfolio of segments with varying tail lengths, the larger the discrepancy between the tail lengths (i.e. between $\nu_j^{(n)}$ and ${\nu_j^{(m)}}$ ), the larger the variation in the common shock contributions. The behaviour of the above ratio has been examined with respect to development factor $\nu_j^{(n)}$ in particular because variation within and across lines of business is mainly driven by the development pattern of claims as explained in section 2.1. As a result, one would expect the development factors to vary the most.

2. – If $p>2$ : The opposite observation is made for the relationship between the above ratio and $\nu_j^{(n)}$ (i.e. the above ratio decreases as $\nu_j^{(n)}$ decreases).

3. – If $p=2$ : In this special case, the common shock contribution is simplified to

   (13) \begin{equation} \dfrac{\dfrac{\gamma^{(n)}}{\beta}}{\dfrac{\gamma^{(n)}}{\beta}+1} \end{equation}
   which is now independent of accident and development periods. Consequently, the common shock contributes proportionately to the total observations over the entire range of the triangles. It is also worth emphasising that specifying $p=2$ gives the multivariate gamma case of the multivariate Tweedie framework.

The above analyses and examples show that the choices of scaling factors $\kappa_{i,j}^{(n)}$ and $\lambda_{i,j}^{(n)}$ are subject to many constraints. To accurately capture the dependence structure, these parameters are required to balance the common shock proportions within all claim observations over the entire range of the triangles. However, this can result in overfitting, which can be a critical issue in loss reserving due to small sample size data. Furthermore, the specification of these parameters may need to be restricted in some cases for the purpose of preserving distributional tractability. It is then the aim of this paper to find a solution that compromises between these conflicting issues with a specific application on the common shock Tweedie approach in Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b).

3.1 Theoretical framework

Claims are first standardised using a common unit of exposure such as the number of claims, or the total amount of premium collected, to ensure consistency across accident periods and business lines. Recall the specification of the common shock Tweedie model in Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b) described in section 2.5,

(14) \begin{align} {Y_{i,j}^{(n)}} & {= \kappa_{i,j}^{(n)}V_{i,j}+Z_{i,j}^{(n)}} \end{align}

where

(15) \begin{align}{V_{i,j}} & {\sim \text{Tweedie}_p(\alpha,\beta) } \end{align}

(16) \begin{align}{Z_{i,j}^{(n)}} & {\sim \text{Tweedie}_p(\eta_i^{(n)}\nu_j^{(n)}, \gamma^{(n)}) }\end{align}

(17) $$\kappa _{i,j}^{(n)}{\rm{ }} = {\rm{ }}{\left( {{\alpha \over {\eta _i^{(n)}\nu _j^{(n)}}}} \right)^{1 - p}}{{{\gamma ^{(n)}}} \over \beta }r{\rm{ }}$$

Recall that $\alpha$ and $\eta_i^{(n)}\nu_j^{(n)}$ are location (mean) parameters, and $\beta$ and $\gamma^{(n)}$ are dispersion parameters of ${V_{i,j}}$ and $Z_{i,j}^{(n)}$ , respectively.

As shown in equation (10), the common shock scaling factor has to be specified in the above form that involves parameters of the common shock $V_{i,j}$ and the idiosyncratic component $Z_{i,j}^{(n)}$ . However, due to the unbalanced feature of reserving data with $\nu_j^{(n)}$ varying across development lag j and business line n, the common shock contributes disproportionately to the total observation $Y_{i,j}^{(n)}$ . It is also desirable to maintain model parsimony.

Given the above considerations, we can replace the non-cell-specific parameter $\alpha$ in the scaling factor with column-specific parameter

(18) \begin{align} {{\alpha_j}} &{={\tilde{c}} \left(\prod_{n}\text{E}\left[Z_{i,j}^{(n)}\right]\right)^{\frac{1}{N}} = {\tilde{c}}\left(\prod_{n}{\eta}_i^{(n)}{\nu}_j^{(n)}\right)^{\frac{1}{N}}} \end{align}

(19) \begin{align} &{\approx c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}} \\[8pt]\nonumber \end{align}

The parameter $\alpha_j$ is also the location parameter of the common shock ${V_{i,j}}$ . As a result, we approximately have

(20) \begin{align} {V_{i,j}} & {\sim \text{Tweedie}_p({\alpha_j},\beta) = { \text{Tweedie}_p\left(c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}},\beta\right)} } \end{align}

Essentially, the common shock parameter $\alpha_j$ is proportional to the geometric average of idiosyncratic components of claims which share the same common shock component. In this case, these are claims in the same accident period and development period as the framework is used to capture cellwise dependence. This geometric average can then be simplified by removing accident period factors because we can reasonably expect limited variation across accident periods as a result of claims standardisation, assuming no significant changes occur across accident periods.

The above specification of scaling factor aims to balance the impact of unbalanced feature in loss reserving data which are mainly introduced by variations in development factors $\nu_j^{(n)}$ . Using this specification, the common shock proportion is given by

(21) \begin{equation} {\dfrac{\left(\dfrac{c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}{\eta_i^{(n)}\nu_j^{(n)}}\right)^{2-p}\dfrac{\gamma^{(n)}}{\beta} }{\left(\dfrac{c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}{\eta_i^{(n)}\nu_j^{(n)}}\right)^{2-p}\dfrac{\gamma^{(n)}}{\beta} +1} } \end{equation}

This does not provide a complete balance of common shock proportions because the effect of $\nu_j^{(n)}$ is reduced by a factor $\sqrt[N]{\nu_j^{(n)}}$ . However, it still provides quite a significant improvement over the original framework. This will be demonstrated in the simulation illustration in section 4. This specification can also preserve distributional tractability of the framework. In addition, model parsimony is retained as the total number of parameters in the framework is unchanged. This can be considered an effective solution given the three constraints discussed in section 2.

In addition to the above treatment for unbalanced data, we also introduce a treatment for negative claims

(22) \begin{align} Y_{i,j}^{(n)} + \xi^{(n)} & = \left(\dfrac{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}{ \eta_i^{(n)}\nu_j^{(n)} }\right)^{1-p}\dfrac{\gamma^{(n)}}{\beta}{V_{i,j}}+Z_{i,j}^{(n)} \end{align}

where a translation factor is used and defined such that

(23) \begin{align} \xi^{(n)} & = \begin{cases} 0 & \text{if } {\min\lbrace Y_{i,j}^{(n)},\,\forall i,\,j \rbrace\ge 0},\\ \ge -\min\lbrace Y_{i,j}^{(n)}\rbrace & \text{if } {\min\lbrace Y_{i,j}^{(n)},\,\forall i,\,j\rbrace < 0} \\ \end{cases} \end{align}

The translation is only needed for a loss triangle if it contains at least one negative value and it must be large enough to offset the smallest negative value. It is worth emphasising that in this case, while its lower bound is deterministic, the actual value of $\xi^{(n)}$ still has to be estimated. The generalisation of this treatment to the general common shock framework in Avanzi et al. (Reference Avanzi, Taylor and Wong2018) is straightforward.

Following from the above specification, the marginal density is then given by

(24) \begin{align} Y_{i,j}^{(n)} + \xi^{(n)} & \sim \text{Tweedie}_p\left( \eta_i^{(n)}\nu_j^{(n)} \left[\left(\dfrac{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}{ \eta_i^{(n)}\nu_j^{(n)} }\right)^{2-p} \dfrac{\gamma^{(n)}}{\beta}+1\right],\right.\nonumber\\ &\qquad \left.\gamma^{(n)}\left[\left(\dfrac{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}{ \eta_i^{(n)}\nu_j^{(n)} }\right)^{2-p} \dfrac{\gamma^{(n)}}{\beta}+1\right]^{1-p}\right) \end{align}

where the first parameter is the location parameter and also the mean of $Y_{i,j}^{(n)} + \xi^{(n)}$ . The second parameter is the dispersion parameter. It follows that the vector of translated claims in the same position across all triangles

(25) \begin{align} {{}_{\xi}\boldsymbol{Y}_{i,j} = \begin{pmatrix} {Y}_{i,j}^{(1)} + {\xi}^{(1)}\\ {Y}_{i,j}^{(2)} + {\xi}^{(2)}\\ \vdots\\ {Y}_{i,j}^{(N)} + {\xi}^{(N)}\\ \end{pmatrix}} \end{align}

has a multivariate Tweedie distribution with the multivariate density

(26) \begin{align} & f_{{}_{\xi}\boldsymbol{Y}_{i,j}}\left(y_{i,j}^{(1)}+\xi^{(1)},...,y_{i,j}^{(N)}+\xi^{(N)}\right) =\nonumber\\ &\qquad\int_{0}^{A_{i,j}} f_{{V_{i,j}}}(w_{i,j})\prod_{n=1}^{N}f_{Z_{i,j}^{(n)}}\left(y_{i,j}^{(n)} + \xi^{(n)}-\left(\dfrac{{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}}{ \eta_i^{(n)}\nu_j^{(n)} }\right)^{1-p}\dfrac{\gamma^{(n)}}{\beta}w_{i,j}\right)dw_{i,j} \end{align}

where

(27) \begin{equation} A_{i,j} = {\min \left(\left(\dfrac{ \eta_i^{(1)}\nu_j^{(1)} }{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}\right)^{1-p}\dfrac{\beta}{\gamma^{(1)}}(y_{i,j}^{(1)}+\xi^{(1)}),...,\left(\dfrac{ \eta_i^{(N)}\nu_j^{(N)} }{ {c\sqrt[N]{\nu_j^{(1)}...\nu_j^{(N)}}}}\right)^{1-p}\dfrac{\beta}{\gamma^{(N)}}(y_{i,j}^{(N)}+\xi^{(N)})\right)} \end{equation}

and where $f(.)$ is the Tweedie density in reproductive form (see also Jørgensen, 1997, Chapter 4).

3.2 Model estimation with Bayesian inference

Bayesian inference is used for model estimation. Bayesian estimation has gained its popularity in the loss reserving literature due to rapid computing advancements and Markov Chain Monte Carlo (MCMC) methods that allow the calculation of intractable posterior densities to be performed significantly faster (Avanzi et al. Reference Avanzi, Taylor, Vu and Wong2016b; Verrall et al., Reference Verrall, Hössjer and Björkwall2012). In addition, the incorporation of prior densities in the calculation of posterior densities is a natural way to allow for parameter error in modelling (Shi et al., Reference Shi, Basu and Meyers2012; England et al., Reference England, Verrall and Wüthrich2012). Another aim of using a Bayesian setup is to also estimate the power parameter p and translation parameter $\xi^{(n)}$ with allowance of parameter uncertainty. This is to formalise the estimation of these parameters as they are often estimated heuristically in practice. It is worth emphasising that the Bayesian structure is not integral to our model but serves as a device for estimation.

A two-step procedure is used for estimation, similar to that in Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b). The first stage is the estimation of all parameters except c and $\beta$ of the common shock ${{V_{i,j}}}$ . This stage, however, gives the estimate of a ratio of these parameters denoted as

(28) \begin{align} \delta = \frac{c^{2-p}}{\beta} \end{align}

as can observed from equation (20). This is followed by the multivariate stage that estimates c and $\beta$ conditional on estimates of other parameters from the first stage. The motivation for this procedure comes from properties of the common shock Tweedie framework. Claim observations in the same position across triangles in this framework follow a multivariate Tweedie distribution, and each observation itself also has a marginal Tweedie distribution. In addition, the multivariate density has an integral calculation, as shown in equation (22). This can prolong the estimation of the posterior density, making the tuning and convergence of MCMC much more difficult.

A Bayesian setup requires the specification of the likelihood functions, prior densities and, if posterior densities are not in closed form, computational algorithms used to approximate them. The likelihood functions follow from equation (20) for the first stage and equation (22) for the second stage.

Prior densities are then specified. Prior densities can be chosen to be informative or uninformative. Uninformative priors assign equal possibilities to all values in the feasible set of parameter values, whereas informative priors convey some prior preference for certain values of the parameters. However, the use of informative priors can significantly improve the convergence rate, especially when the parameter dimension is large (Congdon, Reference Congdon2010). Parameter estimates from univariate Tweedie model (Alai & Wüthrich, Reference Alai and Wüthrich2009) can assist in the specification of informative prior densities for parameters $\eta_{i}^{(n)}$ , $\nu_j^{(n)}$ and $\gamma^{(n)}$ . A preliminary analysis of the dependence structure can help select informative prior densities for the common shock parameters c and $\beta$ . Regarding the prior densities for p and $\xi^{(n)}$ , some constraints need to be taken into account. In particular, p is not defined in (0,1), and $\xi^{(n)}$ has a lower bound as per its specification in equation (19).

Putting together the likelihood and prior specifications, the posterior density in the first stage is given by

(29) \begin{equation} f_{\boldsymbol{\Omega}|\boldsymbol{Y^U}}(\boldsymbol{\Omega}|\boldsymbol{Y^U}) \propto \left(\prod_{i,j,n}f_{Y_{i,j}^{(n)}+ {\xi^{(n)}}}\left({y_{i,j}^{(n)}}+ {\xi^{(n)}}|\boldsymbol{\Omega}\right)\right)f_{p}(p)f_{\boldsymbol{\xi}}(\boldsymbol{\xi})f_{\delta}(\delta)f_{\boldsymbol{ {\eta}}}(\boldsymbol{ {\eta}})f_{\boldsymbol{ {\nu}}}(\boldsymbol{ {\nu}})f_{\boldsymbol{\gamma}}(\boldsymbol{\gamma}) \end{equation}

where

\begin{gather*} \boldsymbol{\Omega} = \begin{pmatrix} p\\ \boldsymbol{\xi}\\ \delta\\ \boldsymbol{ {\eta}} \\ \boldsymbol{ {\nu}}\\ \boldsymbol{\gamma}\\ \end{pmatrix}, \,\boldsymbol{\xi} = \begin{pmatrix} \xi^{(1)}\\ \xi^{(2)}\\ \vdots\\ \xi^{(N)} \end{pmatrix},\,\boldsymbol{\eta_i} = \begin{pmatrix} \eta_i^{(1)}\\ \eta_i^{(2)}\\ \vdots\\ \eta_i^{(N)} \end{pmatrix},\, {\boldsymbol{\eta} = \begin{pmatrix} \boldsymbol{\eta_2}\\ \boldsymbol{\eta_3}\\ \vdots\\ \boldsymbol{\eta_I} \end{pmatrix}},\, \boldsymbol{\nu_j} = \begin{pmatrix} \nu_j^{(1)}\\ \nu_j^{(2)}\\ \vdots\\ \nu_j^{(N)} \end{pmatrix},\, {\boldsymbol{\nu} = \begin{pmatrix} \boldsymbol{\nu_1}\\ \boldsymbol{\nu_2}\\ \vdots\\ \boldsymbol{\nu_J} \end{pmatrix}},\,\boldsymbol{\gamma} = \begin{pmatrix} \gamma^{(1)}\\ \gamma^{(2)}\\ \vdots\\ \gamma^{(N)} \end{pmatrix} \end{gather*}

and where $\textbf{\textit{Y}}^{\boldsymbol{U}}$ is a vector of claim observations in the upper claim triangles.

From the model structure in equation (18), we have that all claims $Y_{i,j}^{(n)}$ are independent conditional on common shock. Hence, the joint likelihood can be written as a product of two separate parts: a product of the densities of claims conditional on common shock, and the density of the common shock. In the first stage of the estimation procedure, the likelihood obtained is the first part of the joint likelihood. As also mentioned earlier, this stage provides the estimates of mean parameters $\boldsymbol{\nu},\,\boldsymbol{\eta}$ and dispersion parameters $\boldsymbol{\gamma}$ of the idiosyncratic variables $Z_{i,j}^{(n)}$ , translation parameters $\boldsymbol{\xi}$ and power parameter p. This stage also provides the estimate of $\delta$ which is a function of parameters c and $\beta$ of the common shock $V_{i,j}$ .

In the second estimation step, we work with the joint likelihood directly since common shock components are not observed. In this step, the estimation of c and $\beta$ is carried out conditioning on estimates of other parameters in the first step, including $\delta$ which is a function of c and $\beta$ . The multivariate Tweedie density of $\textbf{\textit{Y}}_{i,j}$ is used to obtain the likelihood in this estimation. The posterior density in this step is given by

(30) \begin{equation} f_{c|\boldsymbol{Y^U,\Omega}}({c}|\boldsymbol{Y^U,\Omega}) \propto \left(\prod_{i,j}f_{\boldsymbol{\boldsymbol{{}_\xi Y}_{i,j}}}\left(\boldsymbol{{}_\xi y}_{i,j}|c,\boldsymbol{\Omega}\right)\right)f_{c}(c) \end{equation}

The posterior densities in both stages are not in recognisable forms, hence MCMC algorithms are required for the estimation. The MCMC algorithm used is Metropolis-Hastings, which is a popular class of MCMC algorithms when the posterior distribution is not in a recognisable form. Random walk Metropolis-Hastings algorithms are used for marginal estimation and multivariate estimation. Proposal densities are chosen (tuned) so that acceptance probabilities are within desirable ranges. The tuning process can be done manually using classical Metropolis-Hastings algorithms. Alternatively, it can be done automatically in adaptive Metropolis-Hastings algorithms using coerced acceptance rates (Haario et al., Reference Haario, Saksman and Tamminen2001; Vihola, Reference Vihola2012).

4.1 An illustration with unbalanced data and negative claims

A data set consisting of two business lines, one of which has a negative claim observation, is simulated. The two loss triangles are represented in Tables A.1 and A.2 in Appendix A. These two triangles consist of simulated claim observations. Each observation in the triangles is drawn from the multivariate Tweedie model for unbalanced data represented in section 3. For simplicity, these observations are assumed to have been adjusted for changes in exposure across accident years.

The marginal fitting is first performed. Parameters are first transformed using the log transformation, and uniform prior densities are used. 200,000 simulations are run and 100,000 simulations are discarded as the burn-in period. The sample chain is thinned by accepting every 5th iteration to reduce the serial dependence between iterations. MCMC paths of some parameters are given in Figure A.1. A similar procedure is performed for the multivariate estimate. The estimates of c and $\beta$ are obtained from this step. Parameter estimates are provided in Table A.3 in Appendix A.

To evaluate the Bayesian inference used for estimation, we compare the true parameter values with 90% confidence intervals obtained from the posterior distributions of these parameters. The results the true values always lie within the corresponding 90% confidence intervals. This indicates the accuracy of the estimation procedure.

We have calibrated the model on the same simulated data set using sub-triangles of dimension $5\times 5$ to assess the robustness of the proposed calibration method. The results show that true values also fall within the 90% confidence intervals of the estimates from this calibration. However, confidence intervals are generally larger than those from the calibration that uses full-size triangles. This is expected due to higher uncertainty in the estimates coming from smaller sample size.

To evaluate the bias in the resulting reserve predictions, forecasts of outstanding claims using our model are compared with the true forecasts as well as forecasts from a multivariate chain ladder model. The true forecasts are calculated as the expected value of outstanding claims using true parameter values. The multivariate chain ladder model used is the model in Pröhl & Schmidt (Reference Pröhl and Schmidt2005). The results are shown in Table 1.

Table 1. Comparison of outstanding claims forecasts

It can be observed from Table 1 that the true forecasts fall within the 90% confidence intervals of the balanced multivariate Tweedie model forecasts. The forecasts from our model are also closer to the true forecasts than those from the multivariate chain ladder model. However, it is worth noting that the simulated data were generated from the multivariate Tweedie model in this illustration.

We further assess bias in the resulting dependence structure by comparing the true cellwise Pearson correlation coefficients for the outstanding claims and the cellwise Pearson correlation coefficients calculated using the parameter estimates. Residual ratios, defined as ratios of estimated Pearson correlation coefficients to true Pearson correlation coefficients, are provided in Table 2. The ratios are close to 1, indicating that the cellwise dependence in the data is well captured.

Table 2. Residual ratios of estimated Pearson correlation coefficients to true Pearson correlation coefficients

We acknowledge that the use of the two-step Bayesian inference does not provide the full picture due to the dependence between the estimated parameters and the reserve being a non-linear function in terms of these parameters. However, this calibration approach was selected due to a number of advantages as mentioned in section 3. These include overcoming the difficulties in dealing with Tweedie densities which are not in tractable form and enhancing computational speed. The above analyses indicate that

1. – the calibration method can capture the dependence structure well;

2. – the resulting reserve predictions show no apparent bias and they are in line with the chain ladder predictions.

Therefore, even though we may not get the full picture, the above results give us confidence that this would not have a material impact on the performance of the calibration.

4.2 A comparison of performances of multivariate Tweedie models on unbalanced data

A natural question arises regarding the performance of the multivariate Tweedie approach for unbalanced data compared to the original multivariate Tweedie approach introduced in Avanzi et al. (Reference Avanzi, Taylor, Vu and Wong2016b). To be able to assess their performances more accurately, this comparison is performed on a simulated data set whose underlying model is known. True common shock contributions are also known, and these serve as the benchmark for the comparison.

To not put any particular framework at a disadvantage, the synthetic data used for this illustration are simulated from a mixture of models. We deliberately select a (extreme) data set to which neither of the frameworks is properly adapted. In particular, two loss triangles of 10 development periods and 10 accident periods are generated such that the dependence is strong in the first four development periods, and not as strong in the last six periods. The common shock components are generated with column-specific mean parameters $\alpha_j = c_j\sqrt{\nu_j^{(1)}\nu_j^{(2)}}$ with $c_j=0.5$ for $1\leq j\leq4$ , and $c_j=0.02$ for $5\leq j\leq 10$ . The second business line is also simulated to be longer-tailed than the first. Similar to the previous illustration, each observation in the two triangles is drawn from the multivariate Tweedie model for unbalanced data represented in section 3. These observations are assumed to have been standardised for accident year effect for simplicity. The two loss triangles are presented in Tables B.1 and B.2 in Appendix B.

Heat maps of ratios of fitted common shock proportions to true proportions are given in Figure 3 for triangle 1. Fitted values are calculated using posterior median of parameters, and true values are calculated using true parameter values. The modified Tweedie model provides a very good fit for the first four development periods. The goodness of fit is considerably less satisfactory in the later development periods when the true common shock proportion drops. The discrepancy is more significant for the first business line which has shorter tail development. The original common shock Tweedie model provides a poor goodness of fit overall, especially in early development periods. The proportions of common shock are underestimated in early development periods and overestimated in later periods. Even though not reported here, similar results are also observed in heat maps of ratios of fitted common shock proportions to true proportions for triangle 2.

Figure 3 Heat maps of ratios of fitted common shock proportions to true proportions for triangle 1 (top: Tweedie framework modified for unbalanced data, bottom: original common shock Tweedie framework).

Overall the modified Tweedie framework does not fully eliminate the issues of unbalanced data across development periods. However, there is a reduction. The fitting is quite good in early development periods but is less satisfactory in later periods. The use of the geometric average of column factors across multiple triangles may contribute to this performance as the geometric average may not be close to some individual column factors if the development patterns are too different. However, it is worth emphasising that the example used has quite an extreme variation in common shock proportions across development periods, and one should not expect such radical variation in practice. In addition, the poor performance also arises from the discrepancy between the modified model with a constant scaling term c and the true model generating the data (with column-specific scaling term $c_j$ ). With this specification, it is not surprising that the earlier (large) development periods dominate the estimation of c. We do not expect good results because of model misspecification, but we can arrive at two main conclusions: the modified framework outperforms the original framework; and the common shock proportions are mis-estimated in the higher development periods, where amounts are small and do not contribute significantly to total liability. It is also worth noting that the poor estimation of common shock proportion does not affect mean forecasts, only dependency between the triangles, and then only where the magnitudes of the forecasts are small.

5.1 Preliminary analysis

A preliminary analysis is performed to assess the suitability of this data set. This includes the assessment of the tails, as well as the dependence structure.

5.1.1 Analysis of the tails

From the plots of loss ratios provided earlier in Figure 1, it can be observed that the Bodily Injury line has longer claims development than the Accident Benefits line. Tail lengths of the two business lines are also assessed using age-to-age development factors

(31) \begin{equation} {f_j^{(n)} = \dfrac{\sum\limits_{i=1}^{I-j}Y_{i,j+1}^{(n)}}{\sum\limits_{i=1}^{I-j}Y_{i,j}^{(n)}}}\end{equation}

Results are given in Table 3. It can be observed that the development factors of the Bodily Injury dominate those of the Accident Benefits line for all development periods, except in the final year. However, this blip may be a false signal due to the truncation of data at the last development period and only one single observation is made in this final year. Hence, the Bodily Injury line is convincingly longer-tailed than the Accident Benefits line.

Table 3. Claims development factors for each development period

5.1.2 Explanatory dependence analysis

A heuristic dependence analysis is performed by fitting to each line a Tweedie Generalised Linear Model (GLM) with a log-link and the chain ladder mean structure

(32) \begin{equation} {a_i^{(n)} + b_j^{(n)}}\end{equation}

This is to remove fixed accident period and development period effects. Correlations between GLM Pearson residuals of the two lines are given in Table 4. The dependence between residuals is strong and significant after allowing for fixed accident period and development period effects.

To examine whether this strong correlation comes from calendar year effects that can impact both lines simultaneously, we also perform another GLM analysis with an additional fixed calendar year effect in the mean structure

(33) \begin{equation} {a_i^{(n)} + b_j^{(n)} + h_t^{(n)}}\end{equation}

Correlations between GLM Pearson residuals of the two lines are then given in Table 5. The correlation coefficients have been reduced, however, not very significantly.

Table 4. Correlation coefficients between cellwise GLM residuals and their corresponding p-values

Table 5. Correlation coefficients between cellwise GLM residuals and their corresponding p-values after removing fixed calendar year effects

Heat maps of residual ratios are given in Figure C.1 in Appendix C. Residual ratios are defined as ratios of observed values to GLM fitted values with the mean structure specified in equation (29). There are some common cellwise patterns that are quite obvious from the heat maps, for example, low payments in development year 7 compensated by accelerated payments in years 8–9 in the first accident year, payment dips in accident year 4 and development year 2. These show that this data set is suitable for illustration of the model.

Results from the preliminary analysis show that this data set is suitable to be used for illustration of the model.

5.2 Estimation and goodness-of-fit assessment

Bayesian inference is used for estimation. The marginal fitting is first performed. 400,000 simulations are run, and 300,000 simulations are discarded as the burn-in period. The sample chain is thinned by accepting every $5^{\text{th}}$ iteration to reduce the serial dependence between iterations. The multivariate fitting is then performed with 90,000 simulations, and the first 30,000 are discarded as the burn-in period. The chain is then thinned by selecting every $3^{\text{th}}$ iteration. Summary statistics are then computed on these posterior samples. The results are given in Tables C.3 and C.4 of Appendix C.

Marginal and multivariate goodness of fits are assessed. Marginal goodness of fit is assessed using QQ plots of residuals in Figure 4. The plot shows that the fit is quite off in the right tail of the Bodily Injury line and slightly off in both tails of the Accident Benefit line. The goodness of fit in other regions, however, is reasonable. This may be a result of the restriction of using the same power parameter p for both lines. However, the multivariate Tweedie framework still provides marginal flexibility with flexible choices of p. For comparison, similar QQ plots are performed for a common shock normal model in Figure 5. It can be observed that the Tweedie marginals provide a much better fit than the normal marginals (with power parameter $p=0$ ).

Figure 4 QQ plots of residuals from common shock Tweedie model ( $p=1.829$ ).

Figure 5 QQ plots of residuals from common shock normal model.

Multivariate goodness of fit is assessed by comparing the empirical bivariate marginals of real data observations and of back-fitted values. These are obtained using the empirical cumulative distribution functions of claim observations from each triangle. Because of the use of a Bayesian inference, various sets of back-fitted data can be generated. A path is randomly chosen for illustration. Scatter plots of these empirical bivariate marginals are presented in Figure 6. It can be observed that the model can capture the general positive dependence structure in the data.

Figure 6 Plots of empirical bivariate marginals for observed values and back-fitted values.

To look for any trace of dependence not captured by the model, we examine the residuals from model fitting. These residuals are obtained as the differences between observations and fitted values, where the latter are calculated using posterior estimates. The Pearson correlation coefficient of these residuals reduces to 0.1204 (p-value 0.3812). This is much weaker than the correlation coefficient of 0.3416 of GLM Pearson residuals in section 5.1.2 and is also insignificant. The insignificant correlation indicates that our model has explained away most of the dependence in the data.

5.3 Common shock proportions

Predictive distributions of outstanding claim observations in the lower triangles can be calculated using predictive Bayesian inference. Using parameter estimates, the contributions of common shock within each cell in the two triangles are calculated and given in Tables 6 and 7. It can be observed that there is only a very mild variation in the common shock proportions within and across triangles. We can relate this result to the challenges coming from applying a common shock model to loss reserving data which has an unbalanced nature discussed in section 2. It shows that the proposed approach has provided a balance of common shock proportions across all loss cells within and across loss triangles.

Table 6. Proportions of common shock to the expected total observations calculated using parameter estimates – Bodily Injury

Table 7. Proportions of common shock to the expected total observations calculated using parameter estimates – Accident Benefits

5.4 Outstanding claims forecast

To obtain the distributions of the outstanding claims, posterior samples of parameters from the Bayesian inference are used to project claims in lower triangles. This projection utilises the specification in equations (7), (16) and (18). This gives a set of samples of future claims in the lower triangles. Using this set, summary statistics of the total outstanding claims distributions are given in Table 8 and kernel densities of outstanding claims are given in Figure 7. Summary statistics provided include the posterior mean, standard deviation, $\text{VaR}_{75\%}$ and $\text{VaR}_{95\%}$ of the distribution of total outstanding claims for each line, as well as for both lines.

Table 8. Summary statistics of outstanding claims distributions

The empirical bivariate marginals of total reserves are shown in Figure 8. For illustration purpose, we show the scatter plot of total reserves from 1,000 posterior samples. The plot shows a mild positive dependence structure in the total outstanding claims across two lines. This is accompanied by a Pearson correlation of 0.0855 (p-value $<$ 2.2e-16). There is no clear evidence of a concentration in the tail regions of the dependence. There can be diversification across claims within a single loss triangle. Hence, the dependence on the aggregate reserves from each line is mild and tail dependence may not be apparent.

The two business lines do not have a comonotonic dependence structure, and this allows the insurer to gain some diversification benefits when they set their risk margins. Using the specification of a risk margin under APRA’s Prudential Standards GPS 340, we have the following definition of risk margin and diversification benefit

(34) \begin{align}&\text{Risk margin}_{\chi\%}[Y] = \max\left\lbrace \text{VaR}_{\chi\%}[Y] - \text{E}[Y];\,\dfrac{1}{2}\text{SD}[Y]\right\rbrace\end{align}

(35) \begin{align}& \text{Diversification benefit}\nonumber\\&\quad = \dfrac{\left(\text{Risk margin}_{\chi\%}[Y_1]+\text{Risk margin}_{\chi\%}[Y_2]\right) - \text{Risk margin}_{\chi\%}[Y_1+Y_2]}{\text{Risk margin}_{\chi\%}[Y_1]+\text{Risk margin}_{\chi\%}[Y_2]} \times 100\% \\[8pt]\nonumber\end{align}

Risk margin $_{75\%}$ and Risk margin $_{95\%}$ , as well as associated diversification benefits, are provided in Table 9. It can then be observed that diversification benefits can be gained as a result of allowing for (non-comonotonic) dependence across business lines.

Table 9. Risk margin and diversification benefits statistics

Figure 7 Kernel densities of predictive distributions of total outstanding claims in each line of business and in the aggregate portfolio.

Figure 8 Plot of empirical bivariate marginals of total reserves (using 1,000 posterior samples).