Application of a Vine Copula for Multi-Line Insurance Reserving

Abstract

This article introduces a novel use of the vine copula which captures dependence among multi-line claim triangles, especially when an insurance portfolio consists of more than two lines of business. First, we suggest a way to choose an optimal joint loss development model for multiple lines of business that considers marginal distribution, vine copula structure, and choice of family for each pair of copulas. The performance of the model is also demonstrated with Bayesian model diagnostics and out-of-sample validation measures. Finally, we provide an implication of the dependence modeling, which allows a company to analyze and establish the risk capital for whole portfolio.

1. Introduction

There have been various approaches to determine loss reserve for a single line of business with statistical models. According to Mack (1993), the chain-ladder method, which has been used as a rule of thumb for the determination of reserve, can be interpreted as a nonparameteric stochastic model. After that, many univariate stochastic reserving models have been developed and used for both determination of point estimates of reserve and risk management. For detailed discussions on univariate stochastic reserving models, see England and Verrall (2002).

However, it is also certain that most of the insurance companies do not run only a single line of business so one needs to consider possible dependence among the lines of business in reserve modeling. Consequently, stochastic reserving models need to be extended to multivariate frameworks. In this regard, some actuarial literature focused on the extension of the chain-ladder method to multivariate cases, such as Braun (2004), Schmidt (2006), and Merz and Wüthrich (2008) which are followed by Shi et al. (2012), which used a bivariate normal distribution to model multi-line reserves. Besides the dependence among multiple lines of business, there have been some works on the dependence between paid claim triangles and incurred claim triangles such as Zhang (2010) and Merz and Wüthrich (2010).

Note that the methods mentioned above are less flexible since it is not allowed to disentangle marginal distributions and multivariate association structure. Therefore, use of a copula was introduced to the multivariate reserving problem since this allows us to consider marginal distribution and association structure in a separate way. Based on this idea, Shi and Frees (2011) proposed the use of a bivariate Gaussian copula to model dependence among multi-line reserves where the marginal distributions are chosen as Gaussian and gamma, respectively. Peters et al. (2014) considered dependence between paid-incurred claim triangles via bivariate copula. Further, Abdallah et al. (2015) used the idea of the hierarchical archimedian copula to capture calendar year effect and multi-line dependence simultaneously.

Although there have been some works on the multivariate reserving in actuarial literature, applications of all the aforementioned works are restricted only to possible dependence between two lines of business. Indeed, it is natural that a property and casualty insurance company runs more than two lines of business, so we need to incorporate such a high-dimensional dependence structure in our reserve modeling.

In that regard, we apply the idea of a vine copula in this article, which uses bivariate copulas as its building blocks and connects them with vine structure to describe the high-dimensional association in a flexible way. The vine copula has been widely used in the actuarial and financial literature. Loaiza Maya et al. (2015) investigated dependency among the exchange rates of Latin American countries using the vine copula. Reboredo and Ugolini (2016), Arreola Hernandez et al. (2017), and Trucíos et al. (2020) used the vine copula to assess systemic risks due to possible dependence among economic subjects. Surprisingly, usage of the vine copula in property and casualty insurance literature does exist but is scarce. For example, Shi and Yang (2018) used a vine copula to capture the serial dependence of claim amounts and derive the experience ratemaking factor.

This paper has been organized as follows. In Section 2, basic concept of the vine copula is introduced and the model selection procedure to be implemented is specified. In Section 3, we describe the data used for our empirical analysis. In Section 4, we go through the model selection procedure to determine marginal distributions and the vine copula structure to be used in our analysis. In Section 5, we discuss the implications of estimation results from the perspective of enterprise risk management using the predictive distribution of unpaid claims. Section 6 addresses practical issues for implementing the proposed methodology. Finally, we conclude this article in Section 7 by providing some future directions of research.

2. Proposed Methodology

Suppose an insurance company owns a portfolio which consists of N multiple lines of business. By assuming balanced observations in multiple triangles, one can write the multivariate cumulative paid claims as ${\mathbf{Y}}_{ij}=(Y_{ij}^{\left(1\right)},Y_{ij}^{\left(2\right)},\dots ,Y_{ij}^{\left(N\right)})$ where $n=1,\dots ,N$ indicates a claim triangle from the $n^{th}$ line of business, $i=1,\dots ,I$ means the $i^{th}$ accident years, and $j=1,\dots ,J$ denotes $j^{th}$ development lag. In general, it is of interest to predict the cumulative paid claim for the next year given information up to the current year, which can be written as $\mathbb{E}\left[Y_{i,j+1}^{\left(n\right)}|Y_{i,j}^{\left(n\right)}\right]$. Therefore, instead of working with ${\mathbf{Y}}_{ij}$, one can directly model the incremental development of claims, or age-to-age factors as follows:

$${\mathbf{D}}_{ij}=(D_{ij}^{\left(1\right)},D_{ij}^{\left(2\right)},\dots ,D_{ij}^{\left(N\right)}),$$

where $D_{i,j+1}^{\left(n\right)}:=Y_{i,j+1}^{\left(n\right)}/Y_{i,j}^{\left(n\right)}$ so that $\mathbb{E}\left[Y_{i,j+1}^{\left(n\right)}|Y_{i,j}^{\left(n\right)}\right]=Y_{i,j}^{\left(n\right)}\mathbb{E}\left[D_{i,j+1}^{\left(n\right)}\right]$. Since each of $D^{\left(n\right)}={\left(D_{ij}^{\left(n\right)}\right)}_{i=1,\dots ,I,j=1,\dots ,J}$ are observed from the same business line, it is natural that we model the marginal distribution of $D^{\left(n\right)}$ and the dependence structure among $D^{\left(n\right)}$ via copulas to jointly model $D^{\left(n\right)},n=1,2,\dots ,N$.

According to Sklar (1959), if all marginal distribution functions are continuous, then there is a unique function $C:[0,1]\to [0,1]$ such that

$$H(x_1,x_2,\dots ,x_m)=C(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_m\left(x_m\right)),$$

where $F_i$ denotes marginal distribution function of $X_i$ and H denotes joint distribution function of $(X_1,X_2,\dots ,X_m)$. Therefore, use of copulas allows us to capture the association among the joint response random variables, which may follow different marginal distributions. In that sense, by letting $\theta =({\theta }^{\left(1\right)},\dots ,{\theta }^{\left(N\right)})$, where ${\theta }^{\left(n\right)}$ is the parameter vector for the marginal distribution of $D^{\left(n\right)}$ and $\varphi$ is the copula parameter, the likelihood for joint distribution is given as follows:

$$\begin{array}{cc}\hfill \ell (\theta ,\varphi |\mathbf{d})& =\sum _{i=1}^I\sum _{j=2}^{J+1-I}logh\left(d_{ij}^{\left(1\right)},\dots ,d_{ij}^{\left(N\right)}|\theta ,\varphi \right)\hfill \\ & =\sum _{i=1}^I\sum _{j=2}^{J+1-I}\left[\sum _{n=1}^{N}log{f}^{\left(n\right)}\left(d_{ij}^{\left(n\right)}|{\theta }^{\left(n\right)}\right)+log{c}_{\varphi }\left(F^{\left(1\right)}\left(d_{ij}^{\left(1\right)}\right|{\theta }^{\left(1\right)}),\dots ,F^{\left(N\right)}\left(d_{ij}^{\left(N\right)}\right|{\theta }^{\left(N\right)})\right)\right].\hfill \end{array}$$

Here h means the joint density of $(D^{\left(1\right)},\dots ,D^{\left(N\right)})$ and we assume that $d_{ij}^{\left(n\right)}(i=1,\dots ,I,j=2,\dots ,J)$ are independent for fixed n.

Depending on the sign of dependence, we may use different families of copulas. For example, in order to capture positive dependence, one can suggest the following bivariate copulas:

$$\begin{array}{cc}\hfill \mathrm{Gumbel}:C^{\mathrm{G}}(u,v)& =exp\left[-{\left({(-log\left(u\right))}^{\varphi }+{(-log\left(v\right))}^{\varphi }\right)}^{1/\varphi }\right],\varphi \in [1,\infty )\hfill \\ \hfill \mathrm{Clayton}:C^{\mathrm{C}}(u,v)& ={\left\{u^{-\varphi }+v^{-\varphi }-1\right\}}^{-1/\varphi },\varphi \in (0,\infty ).\hfill \end{array}$$

Although those two copulas can only capture positive dependence, one can easily rotate those in order to capture reversed tail behavior and negative dependence as follows:

$$\begin{array}{cc}\hfill {90}^{\circ }\mathrm{rotated}\mathrm{Gumbel}:C^{\mathrm{G}90}(u,v)& =v-C^G(1-u,v),\hfill \\ \hfill {90}^{\circ }\mathrm{rotated}\mathrm{Clayton}:C^{\mathrm{C}90}(u,v)& =v-C^C(1-u,v),\hfill \\ \hfill \mathrm{Survival}\mathrm{Gumbel}:C^{\mathrm{SG}}(u,v)& =u+v-1+C^G(1-u,1-v),\hfill \\ \hfill \mathrm{Survival}\mathrm{Clayton}:C^{\mathrm{SC}}(u,v)& =u+v-1+C^C(1-u,1-v),\hfill \\ \hfill {270}^{\circ }\mathrm{rotated}\mathrm{Gumbel}:C^{\mathrm{G}270}(u,v)& =u-C^G(u,1-v),\hfill \\ \hfill {270}^{\circ }\mathrm{rotated}\mathrm{Clayton}:C^{\mathrm{C}270}(u,v)& =u-C^C(u,1-v).\hfill \end{array}$$

Further, the Gaussian copula and Frank copula are also prevalent choices which can capture both positive and negative dependence:

$$\begin{array}{cc}\hfill \mathrm{Gaussian}:C^{\mathrm{N}}(u,v)& ={\mathsf{\Phi }}_2({\mathsf{\Phi }}^{-1}\left(u\right),{\mathsf{\Phi }}^{-1}\left(v\right);\varphi ),\varphi \in [-1,1],\hfill \\ \hfill \mathrm{Frank}:C^{\mathrm{F}}(u,v)& =-\frac{1}{\varphi }log\left[1+\frac{(exp(-\varphi u)-1)(exp(-\varphi v)-1)}{exp(-\varphi )-1}\right],{\displaystyle \varphi \in \mathbb{R}\setminus \left\{0\right\}},\hfill \end{array}$$

where ${\mathsf{\Phi }}_2(·,·)$ stands for the cumulative distribution function of a bivariate, standard, normal, random variable with correlation $\varphi$ and $\mathsf{\Phi }(·)$ stands for the cumulative distribution function of a univariate, standard, normal random variable.

Given families of copulas, allow us to consider various facets of possible dependence among the lines of business, including the upper and lower tail dependence properties. For example, the Clayton copula can capture a positive association and has lower tail dependence but no upper tail dependence. The Gumbel copula can capture positive associations and has upper tail dependence but no lower tail dependence. Further, both Frank and Gaussian copulas are symmetric and able to capture positive and negative dependence, but they have no tail dependencies. Note that upper tail dependence of original copulas corresponds to lower tail dependence of survival copulas, and vice versa. Further, if an original copula can capture a positive association, then the corresponding ${90}^{\circ }$ or ${270}^{\circ }$ rotated copula can capture a negative association; for example, ${90}^{\circ }$ rotated Clayton and ${270}^{\circ }$ rotated Gumbel copulas can be used to capture negative associations. We refer the readers to Nelsen (1999), Embrechts et al. (2001), and Hua and Joe (2011) for a more detailed explanation. The variability of tail behaviors is potentially important considering the data used in our empirical analysis. Here we capture the dependencies among the incremental development factors of different lines of business so that it is natural that the incremental development factors will be quite large in the initial years and become smaller in the mature years. In this regard, upper tail dependence corresponds to the associations among the loss development factors in the initial years, whereas lower tail dependence corresponds to the associations among the loss development factors in the mature years.

If $N=2$, then the dependence can be captured by a bivariate copula, which has been explored in previous works such as Shi and Frees (2011). However, we may not preclude the possibility that $N>2$, and one can suggest the use of vine copula when we have more than two lines of business to capture the multiple dependence. The concept of vine copula was introduced by Aas et al. (2009) who extended the research conducted by Bedford and Cooke (2001) and Bedford and Cooke (2002) to show how multivariate data, which exhibit complex patterns of dependence in the tails, can be modeled using a cascade of pair-copulas, acting on two variables at a time. Later Czado (2010) and Joe and Kurowicka (2011) developed this idea further and provided us with a definite road map for constructing vine copulas to represent complex multivariate data. The core concept behind construction of vine copulas is that a d-dimensional density $f(x_1,x_2,\dots ,x_d)$ can be represented as a product of pair copula densities and marginal densities. For example, let us consider $d=3$ dimension. A possible way to represent the joint density $f(x_1,x_2,x_3)$ is as follows:

$$f(x_1,x_2,x_3)=f_{3|12}\left(x_3\right|x_1,x_2)f_{2|1}\left(x_2\right|x_1)f_1\left(x_1\right).$$

Additionally, since

$$\begin{array}{cc}\hfill f_{2|1}\left(x_2\right|x_1)& =c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))f_2\left(x_2\right),\hfill \\ \hfill f_{3|12}\left(x_3\right|x_1,x_2)& =c_{13|2}(F_{1|2}\left(x_1\right|x_2),F_{3|2}\left(x_3\right|x_2)|x_2)f_{3|2}\left(x_3\right|x_2),\hfill \\ \hfill f_{3|2}\left(x_3\right|x_2)& =c_{23}(F_2\left(x_2\right),F_3\left(x_3\right))f_3\left(x_3\right),\hfill \end{array}$$

where $c_{ij}(F_i\left(x_i\right),F_j\left(x_j\right))$ is the joint copula density of random variables $X_i$ and $X_j$. One can write out the joint density of $X_1,X_2,X_3$ as follows:

$$\begin{array}{cc}\hfill f(x_1,x_2,x_3)& =f_1\left(x_1\right)f_2\left(x_2\right)f_3\left(x_3\right)\left(\mathrm{marginals}\right)\hfill \\ & \times c_{12}(F_1\left(x_1\right),F_2\left(x_2\right))c_{23}(F_2\left(x_2\right),F_3\left(x_3\right))(\mathrm{unconditional}\mathrm{pair})\hfill \\ & \times c_{13|2}(F_{1|2}\left(x_1\right|x_2),F_{3|2}\left(x_3\right|x_2)|x_2)(\mathrm{conditional}\mathrm{pair}).\hfill \end{array}$$

Note that in general, $c_{13|2}$ is affected by $x_2$, not only through $F_{1|2}\left(x_1\right|x_2)$ and $F_{3|2}\left(x_3\right|x_2)$, but also directly through $x_2$. However, such a general form of conditional coplula density is cumbersome to be used in statistical inference. For this reason, it is tempting to use a simplified form of conditional copula density so that $c_{13|2}$ is affected by $x_2$ only through $F_{1|2}\left(x_1\right|x_2)$ and $F_{3|2}\left(x_3\right|x_2)$. Haff et al. (2010) showed that we may write each conditional pair copula density in a simplified form for elliptical distributions with a positive definite scale matrix. Further, it was also shown that use of the simplified form could be a good approximation even in a case wherein such an assumption could not be fully validated. Hence, we will use a simplified form of conditional copula density in this article hereafter.

The sequential dependence structure via pair copulas is described with a regular vine, or R-vine. Let $\mathcal{V}=(T_1,\dots ,T_{d-1})$ be a set of trees where each tree $T_j=(N_j,E_j)$ is connected and $N_1=\{1,\dots ,d\}$. Here $N_j$ and $E_j$ mean the nodes and edges of $T_j$, respectively. If $N_j=E_{j-1}$ for $j=2,\dots ,d$, $\mathcal{V}$ is called a vine. Further, if

$$|e\cap e^{\prime }|=1\mathrm{for}\{e,e^{\prime }\}\in E_j\mathrm{and}j=2,\dots ,d-1,$$

(1)

then $\mathcal{V}$ is called a R-vine. There are two important subclasses of R-vine, C-vine and D-vine. If there is a node $n\in N_j$ so that the degree of that node is $d-j$ in an R-vine, it is called a C-vine. D-vine is an R-vine wherein the degree of every node is at most 2.

The above decomposition of the joint density in terms of vine structure is not unique. For example, when $d=3$, we have the following three available R-vine structures as in Figure 1. Note that in the case of 3-dimensional R-vine structures, C-vine and D-vine are identical, which is no more true for $d>3$.

Note that the number of regular vine structures on d-dimensional variables grows super-exponentially so that there are $d!·2^{\left(\genfrac{}{}{0pt}{}{d-2}{2}\right)}$ regular vine structures on d variables, as shown in Joe and Kurowicka (2011). Therefore, it is already meaningless to display every possible vine structure if $d\ge 5$.

Therefore, for the calibration of vine copula model, we need to carefully determine the following components based on observed data:

• Marginal distribution for each of $D^{\left(n\right)}$;
• Optimal d-dimensional vine structure;
• Copula family for each of pair copula in the selected vine structure.

For statistical inference, one may consider either the frequentist or the Bayesian approach and subsequently choose the optimal structure via a model selection procedure. However, it can be quite computationally demanding to choose an optimal model in both approaches. For example, as mentioned in Table 1 of Gruber and Czado (2015), the size of search space of vine copula model for $d=5$ is already (1.3559 × 10${}^{11}$) when we consider seven families of copula for each of the pair copulas, which is definitely computationally infeasible. In this regard, Dissmann et al. (2013) proposed a top-down approach to choose the vine structure from a frequentist view. Gruber and Czado (2018) proposed sequential vine copula model selection from a Bayesian perspective. Schamberger et al. (2017) also considered a Bayesian approach to model dependence structure with a factor copula in order to handle the dimensionality issue. Recently, Kreuzer and Czado (2019) showed that the pair copula family can be chosen together with their parameters using the Hamiltonian MCMC.

The Bayesian approach can be beneficial in reserve modeling since it can incorporate the uncertainty of parameter estimation and reflect it to compute the predictive distribution of future unpaid claims, as proposed in Shi et al. (2012). Therefore, we propose the following two-step approach, which can be considered as a compromise of Dissmann et al. (2013) and Gruber and Czado (2018):

• Estimate the parameters with Bayesian inference and choose marginal distribution based on Bayesian model selection criteria, such as deviance information criterion (DIC) and the logarithm of the pseudomarginal likelihood (LPML). The DIC for each marginal distribution of $n^{th}$ line, proposed by Spiegelhalter et al. (2002), is defined as

  $$\begin{array}{c}\hfill {\mathrm{DIC}}^{\left(n\right)}=-4\int \ell \left({\theta }^{\left(n\right)}\right|{\mathbf{d}}^{\left(n\right)})\pi \left({\theta }^{\left(n\right)}\right|{\mathbf{d}}^{\left(n\right)})d\theta +2\ell \left({\tilde{\theta }}^{\left(n\right)}\right|{\mathbf{d}}^{\left(n\right)}),\end{array}$$

  where ${\tilde{\theta }}^{\left(n\right)}=\mathbb{E}\left[{\theta }^{\left(n\right)}|{\mathbf{d}}^{\left(n\right)}\right]$ and ${\mathbf{d}}^{\left(n\right)}=\{d_{i,2}^{\left(n\right)},\dots ,d_{i,J-i+1}^{\left(n\right)}|i=1,\dots I\}$ are the observed values of incremental loss developments from $n^{th}$ line of business. Since this integral is hardly available in a closed form, one can estimate DIC as follows:

  $${\widehat{\mathrm{DIC}}}^{\left(n\right)}=-4\sum _{s=1}^S\ell \left({\theta }_{\left[s\right]}^{\left(n\right)}\right|{\mathbf{d}}^{\left(n\right)})+2\ell \left(\frac{1}{S}\sum _{s=1}^S{\theta }_{\left[s\right]}^{\left(n\right)}|{\mathbf{d}}^{\left(n\right)}\right),$$

  where ${\left\{{\theta }_{\left[s\right]}^{\left(n\right)}\right\}}_{s=1}^S$ are MCMC samples generated from the posterior density. Note that we prefer models with smaller DIC values.
  LPML is calculated based on the conditional predictive ordinate (CPO), which was proposed by Gelfand et al. (1992) and Geisser (2017). The CPO for $d_{i,j}^{\left(n\right)}$ is defined as follows:

  $${\mathrm{CPO}}_{i,j}^{\left(n\right)}=\int f\left(d_{i,j}^{\left(n\right)}\right|\theta )\pi \left({\theta }^{\left(n\right)}\right|{\mathbf{d}}_{(-i,j)}^{\left(n\right)})d{\theta }^{\left(n\right)},$$

  where $d_{(-i,j)}^{\left(n\right)}={\mathbf{d}}^{\left(n\right)}-d_{i,j}^{\left(n\right)}$. Since CPO is usually not readily available in closed form, one can estimate ${\mathrm{CPO}}_{i,j}^{\left(n\right)}$ as follows, according to Gelfand and Dey (1994):

  $${\widehat{\mathrm{CPO}}}_{i,j}^{\left(n\right)}=S{\left(\sum _{s=1}^S\frac{1}{f\left(d_{(i,j)}^{\left(n\right)}\right|{\theta }_{\left[s\right]})}\right)}^{-1}.$$
  Finally, according to Ibrahim et al. (2014), ${\mathrm{CPO}}_{i,j}^{\left(n\right)}$ can be summarized as LPML as follows:

  $${\mathrm{LPML}}^{\left(n\right)}=\sum _{i=1}^I\sum _{j=2}^{J+1-I}log\left({\widehat{\mathrm{CPO}}}_{i,j}^{\left(n\right)}\right).$$
  We prefer models with larger LPML values.
• Based on the fitted marginal model and corresponding posterior means of marginal parameters ${\widehat{\theta }}^{\left(n\right)}$ for $n=1,\dots ,N$, generate probability integral transfrom (PIT) ${\widehat{U}}_{ij}^{\left(n\right)}=F^{\left(n\right)}\left(d_{ij}^{\left(n\right)}|{\widehat{\theta }}^{\left(n\right)}\right)$ and optimize the following:

  $$\ell \left(\varphi |\widehat{U}\right)=\sum _{i=1}^I\sum _{j=2}^{J+1-I}\left[log{c}_{\varphi }\left({\widehat{U}}_{ij}^{\left(1\right)},\dots ,{\widehat{U}}_{ij}^{\left(N\right)}\right)\right],$$

  (2)

  which ends up with the optimal vine copula structure and copula family for each pair copula using the top-down algorithm of Dissmann et al. (2013). The algorithm starts with the choice of first tree $T_1$ and subsequently chooses $T_j$ for $j=2,\dots ,d-1$ while we keep the constraint of having R-vine structure in (1). The optimal tree $T_j$ is chosen by solving the following optimization problem:

  $$\underset{E_j,N_j}{argmin}\sum _{e\in E_j}{w}_e,$$

  where $w_e$ is pairwise BIC of an edge e in our implementation.
  In our search for the best family for each of pair copula, we use the aforementioned bivariate Gaussian, Frank, Clayton, Gumbel, and their rotated copulas. Note that this approach also utilizes the inference by margin (IFM) method proposed by Joe and Xu (1996), since it decomposes estimation of marginal distribution and copula structure separately for model selection a with lesser computational burden. It should be noted that such computational convenience is obtained at the expense of potential estimation bias, as shown in Louzada and Ferreira (2016).

3. Data Description

For the empirical analysis, we used a publicly available dataset from ACE Limited 2013 Global Loss Triangles, which consists of aggregated claim developments for insurance operations in North America (4 lines of business), insurance operations overseas (3 lines of business), and global reinsurance operations (2 lines of business). A description and the corresponding index for each line of business are given in

Table 1. Description of the lines of business.

Index	Description
1	North American Workers’ Compensation
2	North American General Liability
3	North American Other Casualty
4	North American Non-Casualty
5	Overseas General Casualty
6	Overseas General Non-Casualty
7	Overseas General Personal Accident
8	Global Reinsurance Property
9	Global Reinsurance Non-Property

. We incorporate a 9-dimensional vine copula structure to capture possible dependence among all lines of business. Indeed, the current COVID-19 situation could be a clear example of showing the inappropriateness of ignoring potential dependencies among different countries. Note that the dataset and code for data analysis are attached as Supplementary Materials.

Although we consider all possible lines of business simultaneously here, an actuary who applies vine copulas to model dependence among the lines of business should be careful, since the complexity of the vine copula structure increases super-exponentially as the dimensionality grows. Therefore, it is a role of an actuary to apply his/her industrial experience or knowledge for the dependence modeling. Further, it is also possible to observe more than 10 lines of business according to the classification in Schedule P of the National Association of Insurance Commissioners (NAIC) annual report, including but not limited to:

• Homeowners/farmowners;
• Private passanger auto liability/medical;
• Commercial auto/truck liability/medical;
• Worker’s compensation;
• Special liability;
• Other liability;
• Fidelity/surety.

Therefore, it is also possible to merge some lines of business into one category in order to not only enhance homogeneity and credibility in the analysis but also avoid too much complexity of the vine structure due to the large number of lines of business.

The loss development triangles, which are training sets in terms of predictive analytics, can be expressed as follows:

$${\mathcal{D}}_{1:I}=\{Y_{ij}^{\left(n\right)}:1\le i\le I\mathrm{and}1\le j\le min(I,I+1-i),n=1,\dots ,9\}.$$

(3)

Our task is to predict the cumulative paid claims for the next years, which are described as follows:

$${\mathcal{D}}_{I+k}=\{Y_{ij}^{\left(n\right)}:1+k\le i\le I\mathrm{and}j=I+1+k-i,n=1,\dots ,9\}.$$

(4)

For enterprise risk management (ERM) purposes, insurance companies are usually interested in the distribution of aggregate unpaid claims after one year of loss development, which is determined by incremental paid losses defined as follows:

$$L^{\left(n\right)}=\sum _{i=2}^I\left(Y_{i,I+2-i}^{\left(n\right)}-Y_{i,I+1-i}^{\left(n\right)}\right)=\sum _{i=2}^I{D}_{i,I+2-i}^{\left(n\right)}·Y_{i,I+1-i}^{\left(n\right)},L=\sum _{n=1}^N{L}^{\left(n\right)}.$$

(5)

This formulation validates us to model $D_{i,j}^{\left(n\right)}$ rather than $Y_{i,j}^{\left(n\right)}$ because one can see that

$$\mathbb{E}\left[L|{\mathcal{D}}_{1:I}\right]=\sum _{n=1}^N\sum _{i=2}^I\mathbb{E}\left[D_{i,I+2-i}^{\left(n\right)}|{\mathcal{D}}_{1:I}\right]·Y_{i,I+1-i}^{\left(n\right)}.$$

A simple approach for the analysis of L is the silo approach, which means we just aggregate all lines of business so that the paid losses of the same development lag and accident year are merged and modeled altogether. In this case, one can write L as follows:

$$L=\sum _{i=2}^I\left(Y_{i,I+2-i}^{•}-Y_{i,I+1-i}^{•}\right)=\sum _{i=2}^I{D}_{i,I+2-i}^{•}·Y_{i,I+1-i}^{•},$$

where $Y_{i,j}^{•}={\sum }_{n=1}^N{Y}_{i,j}^{\left(n\right)}$ and $D_{i,j}^{•}=Y_{i,j+1}^{•}/Y_{i,j}^{•}$. However, such an approach ignores the heterogeneity of the lines of business and could be problematic. Suppose there are two lines of business, personal and commercial auto insurances, and reported losses of these lines are given in

Table 2. An illustrative example of heterogeneity among business lines. 

	Paid Losses			Incremental Development
Year	Personal	Commerical	Aggregate	Personal	Commercial	Aggregate
2	5,000,000	800,000	5,800,000	-	-	-
3	5,200,000	1,200,000	6,400,000	3.92%	40.55%	9.84%
4	5,300,000	1,500,000	6,800,000	1.90%	22.31%	6.06%

. Note that the incremental development is calculated as $log{Y}_{j+1}/Y_j$. One can see that due to the difference in volume, the volatility of loss development for the commercial auto line is wiped out if we analyze the aggregated data. Therefore, the vine copula model can be considered as a flexible generalization of two extreme models, silo and independent approaches, so that one can consider possible dependence and heterogeneity among the lines of business simultaneously.

4. Model Selection and Parameter Estimation

For marginal distribution, we apply the idea of the cross-classified model for each line of business as in Shi and Frees (2011) and Taylor and McGuire (2016), which regresses $Y_{ij}^{\left(n\right)}$ to linear predictor ${\gamma }^{\left(n\right)}+{\alpha }_i^{\left(n\right)}+{\delta }_j^{\left(n\right)}$ as follows:

$$\begin{array}{cc}\hfill \mathrm{Lognormal}:& \mathbb{E}\left[log{Y}_{ij}^{\left(n\right)}\right]={\gamma }^{\left(n\right)}+{\alpha }_i^{\left(n\right)}+{\delta }_j^{\left(n\right)},\hfill \\ \hfill \mathrm{Gamma}:& \mathbb{E}\left[Y_{ij}^{\left(n\right)}\right]=exp({\gamma }^{\left(n\right)}+{\alpha }_i^{\left(n\right)}+{\delta }_j^{\left(n\right)}),\hfill \end{array}$$

(6)

where ${\alpha }_i^{\left(n\right)}$ corresponds to the effect for the $i^{th}$ accident year and ${\delta }_j^{\left(n\right)}$ corresponds to the cumulative effect up to $j^{th}$ development lag, for $n^{th}$ line of business, respectively.

Recall that it is in our interests to predict $\mathbb{E}\left[Y_{i,j+1}^{\left(n\right)}|Y_{i,j}^{\left(n\right)}\right]=Y_{i,j}^{\left(n\right)}\mathbb{E}\left[D_{i,j+1}^{\left(n\right)}\right]$. Therefore, one can rewrite (6) in the following way:

$$\begin{array}{cc}\hfill \mathrm{Lognormal}:& \mathbb{E}\left[log{D}_{i,j+1}^{\left(n\right)}\right]={\eta }_{j+1}^{\left(n\right)},\hfill \\ \hfill \mathrm{Gamma}:& \mathbb{E}\left[D_{i,j+1}^{\left(n\right)}\right]=exp\left({\eta }_{j+1}^{\left(n\right)}\right),\hfill \end{array}$$

(7)

where ${\eta }_{j+1}^{\left(n\right)}={\delta }_{j+1}^{\left(n\right)}-{\delta }_j^{\left(n\right)}$ and $D_{i,j+1}^{\left(n\right)}=Y_{i,j+1}^{\left(n\right)}/Y_{i,j}^{\left(n\right)}$. Therefore, ${\eta }_{j+1}^{\left(n\right)}$ can be interpreted as the incremental loss development factor from $j^{th}$ lag to $j+1^{th}$ lag. In this regard, we utilize two candidate distributions in (7) and compare the model selection diagnostics to choose a suitable distribution for each of marginal models.

Further, it is natural to expect that for a fixed accident year, paid claim amounts gradually increase until they are fully mature and developed, while the magnitude of development gets smaller as development lag increases, which means the incremental loss development factor is expected to be greater than 1 but monotonically decreasing until it converges to 1. In terms of cross-classified model, such a statement is equivalent to the following mathematical condition:

$${\eta }_2^{\left(n\right)}\ge {\eta }_3^{\left(n\right)}\ge \dots {\eta }_{L\left(n\right)}^{\left(n\right)}=0,$$

(8)

for a large enough integer $L\left(n\right)$, or there are ${\zeta }_t^{\left(n\right)}\in [0,\infty )$ for $t=2,\dots ,L\left(n\right)$ and $n=1,\dots ,9$ such that ${\eta }_j^{\left(n\right)}={\sum }_{t=j}^{L\left(n\right)}{\zeta }_t^{\left(n\right)}$. Therefore, we may suggest the following four models as the candidates for marginal distribution of each line of business:

• LNU: The cross-classfied lognormal model in (7);
• LNC: The cross-classfied lognormal model in (7) with (8);
• GamU: The cross-classfied gamma model in (7);
• GamC: The cross-classfied gamma model in (7) with (8).

All of the four models are fitted via RStan because of its flexibility. For example, one can incorporate the constraint in (8) by forcing a lower limit of ${\zeta }_t^{\left(n\right)}$ as 0 in RStan and using a diffuse uniform prior on a positive real number for ${\zeta }_t^{\left(n\right)}$, $n=1,\dots ,9$, and $t=1,\dots ,J$. For the unconstrained models, a diffuse uniform prior on the positive real number is directly used for ${\eta }_{j+1}^{\left(n\right)}$, $n=1,\dots ,9$, and $j=2,\dots ,J$. For each marginal model, four chains with 1000 iterates are used; the first 500 iterates are discarded for burn-in, which usually requires computation time for MCMC sampling of less than a second.

Table 3. Estimated regression coefficients for marginal models.

	Unconstrained Lognormal Model
	T1			T2			T3			T4			T5			T6			T7			T8			T9
	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%
${\eta }_2$	0.79	0.73	0.85	1.26	1.19	1.33	0.66	0.62	0.70	0.36	0.32	0.40	0.72	0.70	0.74	0.81	0.79	0.83	0.53	0.53	0.54	1.45	1.28	1.61	1.80	1.74	1.86
${\eta }_3$	0.32	0.26	0.38	0.61	0.54	0.69	0.22	0.18	0.26	0.05	0.02	0.09	0.27	0.24	0.29	0.19	0.18	0.21	0.09	0.08	0.09	0.28	0.10	0.47	0.68	0.62	0.74
${\eta }_4$	0.19	0.13	0.26	0.43	0.35	0.51	0.13	0.09	0.17	0.03	0.00	0.07	0.16	0.14	0.19	0.07	0.05	0.10	0.03	0.03	0.04	0.14	0.02	0.31	0.38	0.31	0.44
${\eta }_5$	0.14	0.07	0.21	0.24	0.16	0.33	0.08	0.03	0.13	0.03	0.00	0.07	0.11	0.08	0.14	0.03	0.01	0.05	0.01	0.01	0.02	0.12	0.01	0.28	0.26	0.19	0.33
${\eta }_6$	0.10	0.02	0.18	0.21	0.12	0.31	0.06	0.01	0.11	0.03	0.00	0.07	0.07	0.04	0.10	0.02	0.00	0.04	0.01	0.00	0.02	0.12	0.01	0.30	0.17	0.10	0.25
${\eta }_7$	0.08	0.01	0.16	0.10	0.02	0.21	0.04	0.00	0.09	0.03	0.00	0.08	0.05	0.01	0.08	0.02	0.00	0.04	0.01	0.00	0.02	0.14	0.01	0.34	0.12	0.03	0.20
${\eta }_8$	0.08	0.01	0.18	0.10	0.01	0.21	0.04	0.00	0.10	0.04	0.00	0.10	0.04	0.01	0.08	0.02	0.00	0.05	0.01	0.00	0.02	0.17	0.01	0.43	0.08	0.01	0.17
${\eta }_9$	0.10	0.01	0.24	0.11	0.01	0.26	0.06	0.01	0.13	0.06	0.00	0.14	0.04	0.00	0.09	0.03	0.00	0.06	0.01	0.00	0.03	0.23	0.02	0.57	0.10	0.01	0.22
	Unconstrained Gamma Model
	T1			T2			T3			T4			T5			T6			T7			T8			T9
	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%
${\eta }_2$	0.81	0.76	0.86	1.27	1.21	1.34	0.67	0.63	0.70	0.37	0.33	0.41	0.72	0.70	0.74	0.81	0.79	0.83	0.53	0.53	0.54	1.60	1.43	1.76	1.81	1.76	1.87
${\eta }_3$	0.32	0.27	0.38	0.62	0.55	0.69	0.22	0.18	0.26	0.05	0.01	0.10	0.27	0.24	0.29	0.19	0.17	0.21	0.09	0.08	0.09	0.29	0.11	0.47	0.68	0.63	0.74
${\eta }_4$	0.19	0.13	0.25	0.44	0.36	0.51	0.13	0.09	0.17	0.03	0.00	0.07	0.16	0.14	0.19	0.08	0.05	0.10	0.03	0.03	0.04	0.14	0.02	0.31	0.38	0.32	0.44
${\eta }_5$	0.14	0.07	0.20	0.25	0.17	0.33	0.08	0.03	0.13	0.03	0.00	0.07	0.11	0.09	0.14	0.03	0.01	0.05	0.01	0.01	0.02	0.12	0.01	0.29	0.26	0.20	0.33
${\eta }_6$	0.09	0.02	0.17	0.21	0.13	0.30	0.06	0.02	0.11	0.03	0.00	0.07	0.07	0.04	0.10	0.02	0.00	0.04	0.01	0.00	0.02	0.13	0.01	0.31	0.17	0.10	0.25
${\eta }_7$	0.07	0.01	0.15	0.10	0.02	0.20	0.03	0.00	0.08	0.03	0.00	0.09	0.04	0.01	0.08	0.02	0.00	0.04	0.01	0.00	0.02	0.14	0.01	0.35	0.12	0.04	0.20
${\eta }_8$	0.08	0.01	0.17	0.10	0.01	0.20	0.04	0.00	0.10	0.04	0.00	0.11	0.04	0.01	0.08	0.02	0.00	0.05	0.01	0.00	0.02	0.18	0.01	0.45	0.08	0.01	0.17
${\eta }_9$	0.10	0.01	0.22	0.10	0.00	0.26	0.05	0.00	0.13	0.06	0.01	0.15	0.03	0.00	0.08	0.02	0.00	0.06	0.01	0.00	0.02	0.25	0.02	0.64	0.09	0.01	0.22
	Constrained Lognormal Model
	T1			T2			T3			T4			T5			T6			T7			T8			T9
	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%
${\eta }_2$	0.79	0.74	0.85	1.26	1.19	1.32	0.66	0.63	0.70	0.36	0.32	0.40	0.72	0.70	0.74	0.81	0.79	0.83	0.53	0.53	0.54	1.45	1.29	1.61	1.80	1.75	1.85
${\eta }_3$	0.33	0.27	0.38	0.61	0.54	0.69	0.21	0.18	0.25	0.07	0.04	0.11	0.27	0.24	0.29	0.19	0.17	0.21	0.09	0.08	0.09	0.33	0.20	0.48	0.68	0.62	0.73
${\eta }_4$	0.20	0.15	0.26	0.43	0.36	0.50	0.14	0.10	0.17	0.05	0.03	0.08	0.16	0.14	0.19	0.08	0.05	0.10	0.03	0.03	0.04	0.21	0.12	0.33	0.38	0.32	0.44
${\eta }_5$	0.15	0.10	0.20	0.27	0.20	0.34	0.09	0.06	0.12	0.04	0.02	0.06	0.11	0.09	0.14	0.04	0.02	0.06	0.02	0.01	0.02	0.16	0.08	0.25	0.26	0.20	0.32
${\eta }_6$	0.11	0.07	0.16	0.20	0.13	0.27	0.06	0.04	0.09	0.03	0.01	0.05	0.08	0.05	0.10	0.02	0.01	0.04	0.01	0.01	0.02	0.12	0.05	0.20	0.18	0.13	0.24
${\eta }_7$	0.08	0.03	0.12	0.13	0.06	0.20	0.04	0.02	0.07	0.02	0.01	0.04	0.05	0.03	0.07	0.02	0.01	0.03	0.01	0.00	0.01	0.09	0.03	0.16	0.12	0.07	0.18
${\eta }_8$	0.05	0.01	0.10	0.08	0.02	0.15	0.03	0.01	0.05	0.01	0.00	0.03	0.03	0.01	0.06	0.01	0.00	0.02	0.00	0.00	0.01	0.06	0.01	0.12	0.08	0.03	0.13
${\eta }_9$	0.03	0.00	0.07	0.04	0.00	0.10	0.01	0.00	0.03	0.01	0.00	0.02	0.02	0.00	0.04	0.01	0.00	0.01	0.00	0.00	0.01	0.03	0.00	0.08	0.04	0.00	0.09
	T1			T2			T3			T4			T5			T6			T7			T8			T9
	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%	Mean	5%	95%
${\eta }_2$	0.81	0.76	0.86	1.27	1.21	1.33	0.67	0.63	0.70	0.37	0.33	0.40	0.72	0.70	0.74	0.81	0.79	0.83	0.53	0.53	0.54	1.60	1.44	1.76	1.81	1.76	1.86
${\eta }_3$	0.33	0.27	0.38	0.62	0.55	0.69	0.22	0.18	0.25	0.07	0.04	0.10	0.27	0.24	0.29	0.19	0.18	0.21	0.09	0.08	0.09	0.33	0.21	0.48	0.68	0.63	0.74
${\eta }_4$	0.20	0.15	0.25	0.44	0.37	0.51	0.13	0.10	0.17	0.05	0.03	0.07	0.16	0.14	0.19	0.08	0.06	0.10	0.03	0.03	0.04	0.22	0.12	0.32	0.38	0.32	0.44
${\eta }_5$	0.15	0.10	0.19	0.27	0.20	0.34	0.09	0.06	0.12	0.04	0.02	0.06	0.11	0.09	0.14	0.04	0.02	0.06	0.02	0.01	0.02	0.16	0.08	0.25	0.26	0.21	0.32
${\eta }_6$	0.11	0.06	0.15	0.20	0.14	0.27	0.06	0.03	0.09	0.03	0.01	0.05	0.08	0.05	0.10	0.02	0.01	0.04	0.01	0.01	0.02	0.12	0.05	0.20	0.18	0.13	0.24
${\eta }_7$	0.08	0.04	0.12	0.13	0.06	0.19	0.04	0.02	0.07	0.02	0.01	0.04	0.05	0.03	0.07	0.02	0.01	0.03	0.01	0.00	0.01	0.09	0.03	0.16	0.12	0.07	0.18
${\eta }_8$	0.05	0.02	0.09	0.08	0.02	0.15	0.03	0.01	0.05	0.01	0.00	0.03	0.03	0.01	0.06	0.01	0.00	0.02	0.00	0.00	0.01	0.06	0.01	0.12	0.08	0.03	0.13
${\eta }_9$	0.03	0.00	0.06	0.04	0.00	0.10	0.01	0.00	0.04	0.01	0.00	0.02	0.02	0.00	0.04	0.01	0.00	0.01	0.00	0.00	0.01	0.03	0.00	0.08	0.04	0.00	0.09

provides a summary of the estimates of parameters in all models for the marginal components, which are summarized with posterior mean and upper/lower bounds of $90\%$ of the Bayesian credible interval for each parameter. One can see that ${\widehat{\eta }}_9^{\left(n\right)}$ tends to be greater than or equal to ${\widehat{\eta }}_8^{\left(n\right)}$ in the unconstrained models, whereas ${\widehat{\eta }}_8^{\left(n\right)}>{\widehat{\eta }}_9^{\left(n\right)}$ for all $n=1,\dots ,9$ in the constrained models, which is a more natural pattern in loss development analysis.

When posterior MCMC samples of parameters are obtained, it is quite important to make sure that these samples converge to (proper) posterior distributions. In this regard, we use $\widehat{R}$ statistics proposed by Gelman and Rubin (1992) and traceplots which enable us to judge the convergence of MCMC sampling with the generated samples. Basically, $\widehat{R}$ statistics tell us that the samples from distinct chains are well-mixed so that the chains are considered converged if $\widehat{R}\simeq 1$.

Table A1. $\widehat{R}$ values for marginal models.

	Unconstrained Lognormal Model
	T1	T2	T3	T4	T5	T6	T7	T8	T9
${\eta }_2$	1.000	1.001	1.000	1.000	1.001	1.002	1.002	1.000	1.003
${\eta }_3$	1.000	0.999	1.000	1.002	1.000	0.999	1.000	1.003	1.002
${\eta }_4$	1.000	1.001	0.999	1.000	1.002	0.999	1.001	1.000	1.002
${\eta }_5$	1.001	1.001	1.002	1.000	1.003	0.999	1.001	1.001	1.000
${\eta }_6$	1.006	1.001	0.999	1.003	0.999	1.000	1.002	1.000	1.000
${\eta }_7$	1.000	1.001	1.002	1.000	1.000	1.003	0.999	1.000	1.001
${\eta }_8$	0.999	1.002	1.000	1.006	0.999	1.001	1.000	1.000	1.000
${\eta }_9$	1.000	0.999	0.999	1.001	1.002	1.002	1.002	1.002	1.001
	Unconstrained Gamma Model
	T1	T2	T3	T4	T5	T6	T7	T8	T9
${\eta }_2$	0.999	1.000	1.000	1.000	1.006	1.000	0.999	1.000	0.999
${\eta }_3$	1.000	1.000	0.999	1.002	1.001	0.999	1.004	1.000	1.003
${\eta }_4$	1.002	1.001	1.004	1.000	1.001	1.000	1.004	0.999	1.001
${\eta }_5$	1.004	1.000	1.002	1.002	0.999	1.001	1.000	1.003	1.002
${\eta }_6$	1.003	1.002	1.002	1.002	1.000	0.999	0.999	1.000	0.999
${\eta }_7$	1.000	1.004	1.001	1.001	1.002	0.999	1.001	1.000	1.000
${\eta }_8$	1.003	1.000	1.003	1.001	1.000	1.003	1.002	1.003	1.000
${\eta }_9$	1.000	1.002	1.000	1.002	0.999	1.004	1.011	1.000	1.001
	Constrained Lognormal Model
	T1	T2	T3	T4	T5	T6	T7	T8	T9
${\eta }_2$	0.999	0.999	0.999	1.000	1.001	0.999	1.000	1.000	0.999
${\eta }_3$	1.002	1.002	1.001	1.001	1.001	1.002	1.003	1.000	1.002
${\eta }_4$	0.998	1.001	1.002	1.000	1.001	0.999	1.000	0.998	1.002
${\eta }_5$	1.003	1.000	0.999	0.999	0.999	1.001	1.001	1.000	1.003
${\eta }_6$	1.002	0.999	0.999	0.999	1.002	1.000	0.999	0.999	1.003
${\eta }_7$	0.999	1.001	1.000	1.000	0.999	0.999	0.999	0.999	1.000
${\eta }_8$	1.001	1.002	1.000	1.000	1.001	0.999	1.000	1.001	0.998
${\eta }_9$	1.001	1.000	1.000	1.001	1.000	0.999	1.001	1.000	0.999
	Constrained Gamma Model
	T1	T2	T3	T4	T5	T6	T7	T8	T9
${\eta }_2$	0.999	1.003	0.999	1.002	1.000	1.002	0.999	0.999	1.000
${\eta }_3$	0.999	1.000	0.999	1.004	1.003	1.001	0.999	1.000	0.999
${\eta }_4$	1.000	0.999	1.003	1.001	1.004	1.002	0.999	0.999	1.004
${\eta }_5$	0.999	1.001	0.999	1.001	1.001	1.002	0.999	0.999	1.000
${\eta }_6$	0.999	1.001	1.000	1.000	1.000	0.998	0.999	1.000	1.000
${\eta }_7$	1.000	1.001	0.999	1.000	0.999	0.999	1.000	1.001	1.000
${\eta }_8$	1.000	0.999	1.000	1.000	0.998	1.000	1.000	1.001	0.998
${\eta }_9$	1.001	1.000	1.001	1.000	1.000	1.001	1.000	1.002	1.000

in Appendix A shows that $\widehat{R}\simeq 1$ for all parameters, which tells us that the chains are well-mixed. One can also see that MCMC chains are well-mixed by looking at selected traceplots of maginal models, as shown in Figure A1. It indicates different MCMC chains end up with similar empirical distributions, which is a necessary condition of the convergence of the MCMC algorithm to the correct posterior distribution.

After the convergence of Bayesian models is assessed, we compare the goodness-of-fit of Bayesian models using DIC and LPML. Based on the DIC and LPML values for fitted marginal models in

Table 4. DIC and LPML for marginal components.

		T1	T2	T3	T4	T5	T6	T7	T8	T9
DIC	LNU	−57.1	−46.5	−90.3	−84.2	−125.1	−140.5	−224	16.8	−61.8
	GamU	−38.6	−7.7	−74.4	−75.3	−105.9	−123.6	−214.1	47.6	−14.9
	LNC	−62.9	−50.5	−95.9	−88.4	−129.2	−144.1	−228.5	11.9	−66.7
	GamC	−44.2	−11.8	−79.2	−79.6	−110.5	−127.3	−218	42.8	−19.2
LPML	LNU	−326.2	−129.3	−117.2	−291.7	−83.2	−127.3	−43.2	−241.8	−155.2
	GamU	−379.2	−168.6	−144.3	−364.8	−103.1	−162	−62.1	−342.7	−215.2
	LNC	−368.1	−132.7	−123.4	−300.6	−82.6	−118.1	−43.4	−242	−144.8
	GamC	−419.6	−176.4	−145.5	−392.5	−100.1	−161.3	−57.5	−350.5	−219.5

, it turns out that lognormal distribution is favored compared to gamma distribution. Further, DIC values of constrained lognormal distribution are always less than ones of unconstrained lognormal distribution, while LPML values have no definitive patterns between constrained and unconstrained lognormal distributions. Therefore, we suggest to use constrained lognormal distribution as the marginal distribution of every line of business, based on the model selection diagnostics and interpretability of the regression coefficients while unconstrained lognormal distributions are considered as a benchmark.

Once a marginal distribution for each triangle is specified, the dependence structure can be modeled with a 9-dimensional vine copula. Since marginal distributions are fixed and calibrated, it is possible to optimize (2). For the optimization routine, we use the algorithm of Dissmann et al. (2013) which is readily available in an R package VineCopula, as a function RVineStructureSelect, which explores both the optimal vine structure and the choice of family for each of the pair copulas sequentially. For a detailed explanation of such an implementation, see Chapter 8 of Czado (2019).

In our search for the vine structure, we also consider the sparsity of the vine copula. Since we need $\left(\genfrac{}{}{0pt}{}{d}{2}\right)$ pair copulas to construct a d-dimensional vine structure, the number of required pair copulas quadratically increases as the dimension increases, which adds both complexity and computational burden to the optimization scheme. Therefore, in this article, we also incorporate the idea of a sparse copula so that a pair copula is considered to be an independent copula unless the estimated association Kendall’s Tau is significantly different from 0 for the pair copula. More specifically, we use the following test statistic T based on the estimated value of Kendall’s Tau, $\widehat{\tau }$ as proposed in Genest and Favre (2007):

$$T=\sqrt{\frac{9m(m-1)}{2(2m+5)}}\times \left|\widehat{\tau }\right|,$$

where T asymptotically follows a standard normal distribution if m is large enough. For more approaches of incorporating sparsity in a high-dimensional vine copula, see Gruber and Czado (2018) and Nagler et al. (2019).

According to the optimization with the data, the optimal vine structure is given

Table 5. Specification of vine structure.

Tree	Edge	Family	$\mathit{\varphi }$	$\mathit{\tau }$
1	9,8	Gumbel	1.87	0.46
1	9,6	${90}^{\circ }$ rotated Clayton	−1.20	−0.37
1	8,1	Frank	9.40	0.65
1	1,7	${270}^{\circ }$ rotated Clayton	−0.20	−0.09
1	7,2	${90}^{\circ }$ rotated Clayton	−0.95	−0.32
1	6,4	Frank	5.87	0.50
1	9,5	${270}^{\circ }$ rotated Gumbel	−1.24	−0.20
1	4,3	Survival Gumbel	1.61	0.38
2	9,4;6	${90}^{\circ }$ rotated Gumbel	−1.35	−0.26
3	9,7;8,1	Survival Clayton	0.52	0.21

which only displays non-independent pair copulas. In the table, $\varphi$ and $\tau$ mean estimated copula parameter and Kendall’s Tau for each pair copula, respectively. Some empirical relationships are observed by the estimated copula parameters, which seem intuitive considering the nature of the lines of business. For example, it might be quite natural to have strongly positive dependence among the development patterns of line 4 (North American non-casualty) and line 6 (overseas general non-casualty) due to inherent characteristics of the same claim type. Dependence between line 8 (global reinsurance property) and line 9 (global reinsurance non-property) is captured by the Gumbel copula, which indicates there is positive association between claim settlements in the initial years (that corresponds to upper tail part) for a global reinsurance business unit. Such a positive tail dependence might have originated from an internal decision on the claim processing in the same business unit. Note that a company has limited resources to deal with ligitations related to the claim adjustments so that claim adjustments can be delayed in a line by focusing on expediation in claim adjustments in another line, which might explain negative associations between lines 6 and 9, lines 2 and 7, and lines 5 and 9. One can see that the estimated vine structure exhibits enough of a degree of sparsity because only 10 pair copulas are non-independent among $\left(\genfrac{}{}{0pt}{}{9}{2}\right)=36$ pair copulas in the vine structure so that most of the dependencies arise from tree 1, as shown in Figure 2. There are nine lines of business that correspond to the nodes of tree 1 so that we have nine nodes and eight edges as an R-vine structure which captures the first layer of dependence among the lines of business. Since the dependence of a pair of copulas is quite weak except for those in tree 1, we only display the R-vine structure for tree 1.

We also visualize the dependence structure via normalized contour plots in Figure 3, which displays the association between two probabilty integral transforms in a normalized scale. For example, if a pair copula $c(u_1,u_2)$ is given, then the corresponding normalized contour plot gives us a contour plot with the following density g:

$$g({\mathrm{\Phi }}^{-1}\left(u_1\right),{\mathrm{\Phi }}^{-1}\left(u_2\right))=c(u_1,u_2)\varphi \left({\mathrm{\Phi }}^{-1}\left(u_1\right)\right)\varphi \left({\mathrm{\Phi }}^{-1}\left(u_2\right)\right).$$

Normalized contour plots for trees 1, 2, and 3 in the chosen vine structure are provided in Figure 3. One can see that some pair copulas (for example, $c_{4,3}$ and $c_{8,1}$) demonstrate strong tail dependence or skewed shape of the contours, which cannot be captured with naive use of eillptical copulas such as Gaussians or $t-$copulas. Therefore, Figure 3 substantiates that the choice of pair copula families are enough to describe the dependence structure observed in the data. In the figure, the first row corresponds to normalized contour plots for tree 3 where one can see concentric circles as contour plots, which indicates independence of corresponding pair copulas. On the other hand, one can see many skewed contour plots for tree 3 from the third row of the figure, which indicates varying dependence of corresponding pair copulas.

5. Validation and Actuarial Implication

Once a predictive model is calibrated based on the upper loss triangles ${\mathcal{D}}_{1:9}$, which is defined in (3), it needs to be validated based on available ${\mathcal{D}}_{9+k}$, which is defined in (4). For this purpose, ${\mathcal{D}}_{10}=\{Y_{ij}^{\left(n\right)}:2\le i\le 9\mathrm{and}j=11-i,n=1,\dots ,9\}$, the latest diagonals are used as a validation set.

Under the proposed lognormal model for marginal components, one can show that

$$\mathbb{E}\left[Y_{i,j+1}^{\left(n\right)}|Y_{i,j}^{\left(n\right)}\right]=\mathbb{E}\left[D_{i,j+1}^{\left(n\right)}\right]Y_{i,j}^{\left(n\right)}=exp\left({\eta }_{i,j+1}^{\left(n\right)}+\frac{1}{2}{\sigma }^{2\left(n\right)}\right)Y_{i,j}^{\left(n\right)},$$

so that a point estimate of $Y_{i,11-i}^{\left(n\right)}$ is given as ${\widehat{Y}}_{i,11-i}^{\left(n\right)}=Y_{i,10-i}^{\left(n\right)}\times exp\left({\widehat{\eta }}_{11-i}^{\left(n\right)}+\frac{1}{2}{\widehat{\sigma }}^{2\left(n\right)}\right)$, where ${\widehat{\eta }}^{\left(n\right)}$ and ${\widehat{\sigma }}^{\left(n\right)}$ are obtained as posterior means from the MC samples of parameters. Since it is of interest to predict the future unpaid claim L as defined in (5), here we apply three models to project the unpaid claims; the independent model which assumes independence among all lines of business, the silo model which assumes a perfect positive association among all lines of business, and the copula model which utilizes the vine copula structure specified and optimized in Section 4.

Due to linearity of expectation, it is natural to expect that point estimates of reserves under the copula model would be more or less the same as the estimates under independent model. However, insurance companies are often interested in not only the predictive mean, but also the confidence interval of the estimated mean and risk measures for enterprise risk management. In this regard, the proposed vine copula model allows us to consider the impacts of associations among the lines of business that might lead to more accurate estimation of risk measures. For example, if there are positive associations among the lines of business, then the independent model will underestimate Value at Risk (VaR) or Conditional tail expectation (CTE) since a worse scenario can affect all lines of business adversely.

Table 6. Summary of unpaid claims prediction with the fitted models.

	Unconstrained		Constrained
	Copula	Silo	Copula	Silo	Actual
AY = 2005	529,426	219,222	148,767	81,657	90,662
AY = 2006	264,689	167,787	185,653	131,065	109,420
AY = 2007	299,024	196,099	319,881	227,435	209,649
AY = 2008	454,503	404,983	461,885	414,982	338,359
AY = 2009	405,808	386,094	447,697	409,172	304,544
AY = 2010	651,228	610,943	703,166	615,770	450,724
AY = 2011	942,324	1,071,867	986,753	1,070,549	762,328
AY = 2012	2,182,202	2,649,858	2,181,708	2,648,604	1,837,070
Total	5,729,205	5,706,854	5,435,511	5,599,236	4,102,756

summarizes the point estimates of unpaid claims for the subsequent calendar year based on the copula model and the silo model. Note that point estimates of unpaid claims from the independent model should be the same as ones from the copula model, as long as we model each line of business with the same marginal distribution. Further, although we primarily use constrained models for the marginal components, here we also provide the point estimates of unpaid claims from the models where we do not impose any constraints on the parameters.

From the table, it is shown that copula model with constraints turns out to be the best model in terms of prediction of aggregate unpaid claims for the subsequent year. We observe that under unconstrained models, the predicted values of unpaid claims can be severely exaggerated, especially for mature years (in our case, AY = 2005). That is because the value of ${\eta }_9^{\left(n\right)}$ can be overestimated in the unconstrained models, as observed in Table 3. Therefore, it implies that naive implementation of unconstrained model for multi-line reserving problem may end up with poor prediction.

In order to evaluate the prediction performance, here we use validation measures such as root mean squared error (RMSE) and mean absolute error (MAE) defined as follows:

$$\begin{array}{c}\hfill \mathrm{RMSE}=:\sqrt{\frac{1}{8}\sum _{i=2}^9{({\widehat{Y}}_{i,11-i}^{•}-Y_{i,11-i}^{•})}^2},\mathrm{MAE}=:\frac{1}{8}\sum _{i=2}^9|{\widehat{Y}}_{i,11-i}^{•}-Y_{i,11-i}^{•}|,\end{array}$$

where ${\widehat{Y}}_{i,11-i}^{•}$ is the estimate of $Y_{i,11-i}^{•}$ under each method. We prefer the models with smaller values of RMSE and/or MAE. According to

Table 7. Summary of validation performance measures of the fitted models.

	Unconstrained		Constrained
	Copula	Silo	Copula	Silo
RMSE	234,540	318,849	190,381	315,934
MAE	203,306	203,900	166,594	189,311

, the copula model with constraints still turns out to be the best model in terms of RMSE and MAE, which evaluate the performance of predicting the point estimates.

However, it is of less concern to get point estimates of L, the aggregate unpaid claims. Actuaries (and companies) want to know a range of estimates, or the predictive distribution of the aggregate unpaid claims to assess the quality of decision making on reserve and retain appropriate amounts of risk capital for the whole company. Thus, here we provide a way to simulate predictive samples of $L^{\left(n\right)}$ and describe the predictive distribution of $L^{\left(n\right)}$ for $n=1,\dots ,9$ and subsequently $L={\sum }_{n=1}^9{L}^{\left(n\right)}$ with the following steps, which are similar to those of Gao (2018):

• Generate a 9-dimensional uniform random vector $(u_{ij}^{r:\left(1\right)},\dots ,u_{ij}^{r:\left(9\right)})$ based on the specified copula structure in a sequential way. For example, one can first generate $(u_{ij}^{r:\left(1\right)},u_{ij}^{r:\left(2\right)})$ from $C_{12}$. After that, $u_{ij}^{r:\left(3\right)}$ is generated subsequently due to the following identity:

  $$\mathcal{U}[0,1]\stackrel{d}{=}{w}_3=C_{3|12}\left(u_3\right|u_1,u_2)=h_{3|1;2}\left(C_{3|2}\left(u_3\right|u_2),C_{1|2}\left(u_1\right|u_2)\right),$$

  where $h_{3|1;2}(w,v)=\frac{\partial }{\partial w}{C}_{13;2}(w,v)$. One can continue such implementation to obtain a random sample of d-dimensional uniform vector, which is readily available with RVineSim function in R package VineCopula. (For details, see Chapter 6 of Czado (2019).)
• Based on the lognormal assumption of marginal components, $(D_{i,j}^{r:\left(1\right)},\dots ,D_{i,j}^{r:\left(9\right)})$ are generated as follows; $D_{i,j}^{r:\left(n\right)}=exp\left({\widehat{\eta }}_j^{r:\left(n\right)}+{\mathrm{\Phi }}^{-1}\left(u_{ij}^{r:\left(n\right)}\right){\widehat{\sigma }}^{r:\left(n\right)}\right)$ where ${\widehat{\eta }}_j^{r:\left(n\right)}$ and ${\widehat{\sigma }}^{r:\left(n\right)}$ are $r^{th}$ MC samples of ${\eta }_j^{\left(n\right)}$ and ${\sigma }^{\left(n\right)}$ from the marginal model for the $n^{th}$ line of business, respectively.
• Repeat steps 1 to 2 to get $L^{r:\left(n\right)}$, the MC samples of $L^{\left(n\right)}$ for $r=1,\dots ,R$ where

  $$L^{r:\left(n\right)}=\sum _{i=2}^9{\widehat{Y}}_{i,11-i}^{r:\left(n\right)}-Y_{i,10-i}^{\left(n\right)}=\sum _{i=2}^9\left(D_{i,11-i}^{r:\left(n\right)}-1\right)Y_{i,10-i}^{\left(n\right)}\mathrm{and}{L}^{r:•}=\sum _{n=1}^9{L}^{r:\left(n\right)}.$$

Since $\{Y_{i,10-i}^{\left(n\right)}:2\le i\le 9,n=1,\dots ,9\}\subset {\mathcal{D}}_{1:9}$, only $\{D_{i,11-i}^{\left(n\right)}:2\le i\le 9,n=1,\dots ,9\}$ needs to be simulated to describe the distribution of L. In the simulation scheme, possible uncertainties of estimated parameter values are already considered since the marginal parameters are estimated with the Bayesian approach.

Figure 4 shows the predictive kernel densities of L based on the Monte Carlo samples generated with the aforementioned approach from the constrained model. In the silo approach, all nine triangles are aggregated into one triangle based on the development lags and accident years, so that loss development patterns of all business lines are implicitly assumed to be identical, which may not be true in some cases. One can see that the predictive distribution of L under the copula model is more jagged than one under the independent model, but the copula model is still closer to the independent model than the silo model. It agrees with the proposed vine copula structure described in Figure 2, which shows a sufficient degree of sparsity.

The discrepancy among predictive distributions based on copula, independent, and silo models also can be quantified via Hellinger distance. Hellinger distance, originally proposed by Hellinger (1909), is used to quantify the "distance" between two distributions. For example, if there are two functions f and g which are probability density functions of two distributions F and G, the square of Hellinger distance is calculated as follows:

$$H^2(F,G)=1-\int \sqrt{f\left(x\right)g\left(x\right)}dx.$$

One can see that $0\le H(F,G)\le 1$ is always satisfied with Hellinger distance so that it is well calibrated and interpreted, where $H(F,G)=0$ means F and G are identical almost surely and $H(F,G)=1$ means F and G are totally different. Due to this property, Hellinger distance has been widely used in the statistics literature, such as Beran (1977), and estimation of $H(F,G)$ using generated samples from F and G can be easily implemented via R packages such as statip. Here, ${\widehat{H}}^2\left(L\right|\mathrm{Copula},L\left|\mathrm{Independent}\right)$, the square of the Hellinger distance between the kernel density of L under copula model with constraints (illustrated with blue solid line in Figure 4) and the density under the independent model with constraints (illustrated with black dotted line in Figure 4) is about $4.33\%$ while ${\widehat{H}}^2\left(L\right|\mathrm{Copula},L\left|\mathrm{Silo}\right)=13.54\%$. Therefore, one can see that the proposed vine copula structure can capture weak dependencies among the lines of business.

Note that Hellinger distance is a special case of the functional Bregman divergence proposed in Goh and Dey (2014) and Jeong (2020), which is given as follows:

$$D_{\psi }(F,G)=\int \psi \left(\frac{f\left(x\right)}{g\left(x\right)}\right)g\left(x\right)dx,$$

where $\psi (·)$ is a strictly convex and differentiable function with positive support and $\psi \left(1\right)=0$. One can see that $D_{\psi }(F,G)=H^2(F,G)$ if $\psi \left(z\right)=1-\sqrt{z}$, and $D_{\psi }(F,G)$ is Kullback–Leibler divergence if $\psi \left(z\right)=-logz$.

Based on the predictive distributions of L and $L^{\left(n\right)}$ for $n=1,\dots ,9$, one can calculate the risk margins for the unpaid claims as in

Table 8. Estimated risk margins for the unpaid claims.

	VaR			CTE
	90%	95%	99%	90%	95%	99%
Marginal T1	384,383	420,016	494,059	432,782	464,901	521,869
Marginal T2	1,403,464	1,520,387	1,720,617	1,560,576	1,657,167	1,893,831
Marginal T3	655,532	694,812	766,766	707,338	740,161	805,153
Marginal T4	1,657,500	1,802,081	2,126,876	1,858,735	1,993,388	2,313,396
Marginal T5	777,885	807,011	869,457	818,302	844,170	895,520
Marginal T6	691,142	714,401	757,403	721,538	741,557	791,119
Marginal T7	405,697	413,424	428,243	416,023	423,117	439,641
Marginal T8	586,402	696,629	1,009,717	757,769	878,869	1,199,343
Marginal T9	556,892	591,841	658,353	603,369	633,144	703,385
Marginal Total	7,118,897	7,660,602	8,831,491	7,876,431	8,376,473	9,563,257
Independent Aggregate	6,163,307	6,363,592	6,755,125	6,428,484	6,608,753	6,978,484
Copula Aggregate	6,145,848	6,376,847	6,774,463	6,435,712	6,610,684	6,995,674
Silo Aggregate	6,434,983	6,692,370	7,309,860	6,815,535	7,076,775	7,586,896

. Both VaR and CTE are estimated in an empirical way—in other words, based on the generated samples of the predictive distribution of L or $L^{\left(n\right)}$. Usually, it is expected that a company can get a benefit from the diversification of risks for multiple lines of business so that the risk capital for aggregate reserve is less than the sum of the risk capitals for the reserve of each line, which is known as subadditivity of risk measure. Although subadditivity is not guaranteed for VaR, it is guaranteed for CTE and we can quantify the diversification effects by taking the differences of risk capitals by subtracting the risk capital for aggregate reserve from the sum of the risk capitals for the reserve of each line.

Table 9. Risk reduction from diversification.

	VaR			CTE
	90%	95%	99%	90%	95%	99%
Independent	955,590	1,297,010	2,076,366	1,447,947	1,767,720	2,584,773
Copula	973,049	1,283,755	2,057,028	1,440,719	1,765,789	2,567,583
Silo	683,914	968,232	1,521,631	1,060,896	1,299,698	1,976,361

summarizes measured diversification effects for given models. We can observe that the measured diversification effect of the copula model is much greater than that of the silo model because of weak dependence captured in the vine structure, although it is still less than that of the independent model, which assumes perfect independence among the lines of business.

6. Practical Issues for Implementation

Note that a loss reserving dataset from a single company was utilized in this article because a primary insurer will be more interested in analysis of dependence among the lines of business in the same company for effective enterprise risk management, rather than whole industry. However, this methodology can be applied for analysis of reserving trend in a dataset of several companies as well. For example, analyzing dependencies among the loss development profiles of several companies can be of interest for the regulatory authorities and reinsurance companies.

One can be also concerned about the predictive ability of these models, which can be distorted by catastrophic claims. In this regard, insurance companies often decompose the loss development profiles into two layers, primary and excess layers, to minimize the distortion impact from the presence of catastrophic claims (Dew and Hedges 1998), and our proposed methodology can be applied to analyze dependencies among the loss development for the primary layer. Note that such distinction of layers is usually internal and may not be disclosed to the public, which also applies to the dataset used in this article.

7. Conclusions

In this article, we explored and introduced a novel approach which considered possible dependence among the multiple lines of business via vine copula. Use of a vine copula is very flexible and it allows us to model high-dimensional dependence among the business lines, not only bivariate dependence. In the case of traditional copula models, one needs to calibrate many copula models with different families, and go through model validation procedure to choose the best one among the candidates. On the other hand, our proposed vine copula structure enables us to explore the optimal copula structure for the given data in a unified manner.

For model selection, a stepwise approach can be applied to choose the marginal distributions, copula structure, and family for each pair of copulas subsequently. Our empirical analysis on a synthetic insurance portfolio consisting of nine lines of business showed weak associations among the lines.

Further, it was also shown that the naive implementation of a cross-classified model may result in a counterintuitive pattern of loss development, so one might consider a constrained cross-classified model to mitigate this issue. In our work, constraints on the development lag parameters are naturally incorporated via prior elicitation in a Bayesian framework. We expect that a more thorough discussion of the constrained cross-classified model in terms of variable selection could be one of the future directions of this work.