Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
95,515 characters · 17 sections · 36 citation commands
Bivariate Distribution Regression with Application to Insurance Data
Data with bivariate discrete and continuous outcomes are often encountered in various areas, including economics, insurance risk analysis, and production management. Characterizing the dependencies among the outcomes and their joint distributional features is a crucial yet challenging task due to their complexities, especially when conditioning on observed variables. In insurance data analysis, companies maintain a record of the number of claims and the average claim amount from their policyholders of non-life insurance. The study of their joint distribution, conditional on policyholders' attributes, plays an important role in insurance companies' decisions and access risks. Note that there are a plethora of financial/insurance network studies tang2022insurance following the collapse and near-failure of the insurance giant American International Group in 2008. Our proposed method analyses the basic risk of the insurance business that is the foundation of these networks.
This paper proposes a semiparametric estimation method for the conditional joint distribution of bivariate outcomes, using the distribution regression (DR) approach. williams1972analysis introduce the DR approach to analyze ordered categorical outcomes by using multiple binary regressions. foresi1995conditional first extend the DR approach to characterize a conditional distribution, and various studies hall1999methods, chernozhukov2013inference, rothe2013misspecification propose the DR approach in different contexts.
In this paper, we extend the existing research to the conditional joint distribution of a pair of discrete and continuous outcomes. We first apply the factorization formulation of the bivariate joint distribution and then use the DR method to separately estimate two conditional distributions: the distribution of the discrete outcome conditional on covariates and the distribution of the continuous outcome conditional on the covariates and the discrete outcome. Incorporating the discrete outcome as an additional covariate allows us to characterize the dependency between the two outcomes conditional on the covariates in a simple yet flexible manner. Moreover, the combination of the two regression results can uncover the joint conditional distribution and its characteristics. The joint conditional distribution given a set of covariates can provide a robust statistical basis for computing and optimizing conditional risk measures such as the Value-at-Risk and Expected Shortfall noyan2013optimization.
Our approach addresses several issues which may be considered outstanding in multivariate modeling. First, the estimation method in this paper is semiparametric, in that a collection of binary outcome regressions are used to characterize a joint conditional distribution through the factorization formulation, instead of imposing global parametric restrictions as of the vast existing literature olkin1961multivariate,cox1992response,gueorguieva2001correlated. Hence, it is useful when researchers know little about the underlying distributions and their parametric form.
Next, our method can be applied to discrete, continuous, or mixed distribution outcomes and flexibly accommodate their dependencies. One of the popular strategies to construct the multivariate distributions conditional on covariates is the copula regression models shi2018pair. The flexibility of copula models mainly lies in the possibility of specifying the marginals and the dependency among outcome variables separately. The theoretical foundation for the application of copulas is Sklar's theorem, which guarantees the uniqueness of the copula function for continuous outcome variables. In the presence of discrete or mixed distribution outcomes, however, the copula function is not unique. To address this issue, yang2020nonparametric study nonparametric estimation of copulas for discrete outcomes. For multivariate mixed outcomes, YangLu2020NCEf proposes a nonparametric estimator of copulas with the marginal specification based on standard parametric mixed distribution. This specification implies the proportion of zeros in the mixed variables plays a key role in the finite sample performance of the estimator. Our proposed method can serve as an alternative approach in these circumstances.
Third, our approach is easy to implement and computationally fast among the class of semiparametric and nonparametric methods, even when the number of covariates is moderately large. In fact, one can use standard statistical software to implement our approach by fitting only two sets of parametric binary regressions locally over the space of outcome variables. In those local models, parameters can be considered as “pseudo-parameters", which extract local information of distributions of interest and can be estimated at the parametric rate under certain regularity conditions White1982. We provide the limiting distribution of our estimator, while the limit process depends on unknown nuisance parameters. To circumvent the issue of nonpivotal limit processes, we consider the exchangeable bootstrap praestgaard1993Exchangeably and show its validity, extending the result of chernozhukov2013inference.
Our paper complements and extends the recent studies on the multivariate extensions of the DR method. meier2020multivariate proposes a method to estimate the joint conditional distribution function by directly applying the DR approach with an indicator function over a multi-dimensional grid. As we discussed in Section (ref), this direct application could face a practical issue even when the number of girds is moderately large. In contrast, our factorization approach resolves this issue by estimating univariate conditional distribution sequentially. Also, klein2022multivariate introduce a DR-type approach that imposes a global structure on the conditional joint distribution function, whereas our approach is semiparametric.
We conduct extensive simulation studies to examine the finite-sample performance of our proposed DR approach under various data generating processes. In particular, we consider two popular parametric models in the insurance literature, the hierarchical model garrido2016generalized and the copula model czado2012mixed, and a non-standard distribution constructed through transforming a bivariate Gaussian density. The simulation results show that the proposed method performs consistently well across all these setups, whereas the existing parametric approaches perform well only when the model is correctly specified. These results underscore the importance of our semiparametric approach in the context of finite samples.
For empirical application, we analyze a French insurance portfolio. For each policyholder, their characteristics are collected together with their past claim experience, which includes the discrete number of claims made (claim frequency) and the average cost per claim (severity). Traditional parametric approaches have been widely used in the insurance literature czado2012mixed, and the increased complexity of insurance data has been driving the development of nonparametric methods. The average cost per claim follows a mixed distribution: a probability mass at zero corresponding to no claims and an otherwise positive claim from a skewed and long-tailed distribution. Hence, a naive specification of parametric claim distributions is often unsatisfactory. The data is very large, consisting of more than 400,000 observations, yet most policyholders did not report any claims. In this case, the copula approach is less robust in capturing the joint distribution YangLu2020NCEf. More importantly, the average severity exhibits a clear multi-modality in the dataset and high skewness when excluding the zero-count observations. This raw data feature suggests that the existing popular parametric approaches are insufficient in this case. The proposed DR approach demonstrates superior performance in both in-sample and out-of-sample results against the existing parametric hierarchical and copula models.
The rest of the paper is organized as follows. In Section 2, we present the DR approach for modeling the bivariate discrete and continuous outcomes. Section 3 sets out the asymptotic properties of our DR method. We provide simulation results in Section 4. In Section 5, we compare our proposed method against existing approaches in application to study a real insurance data set. We conclude this paper in Section 6. The proof of the main results is given in Appendix.
In this section, we illustrate how the DR approach characterizes a univariate conditional distribution. In the following, we let $W$ be an outcome variable with support $\mathcal{W}\in\mathbb{R}$, which can be discrete, continuous, or mixed, and let $X$ be a $d_x\times1$ vector of covariates with support $\mathcal{X}\in \mathbb{R}^d_x$. The DR approach models the conditional CDF of $W$ given $X=x$ by fitting a parametric linear-index model targeting an arbitrary location of the outcome. More specifically, letting $\Lambda: \mathbb{R} \to [0,1]$ be a known link function, we model the conditional distribution function as, for $(x, w) \in \mathcal{X}{\times}\mathcal{W}$.
where $P: \mathcal{X}\rightarrow \mathbb{R}^{d}$ is a known transformation of the conditioning variables, $\alpha(w)\in\mathbb{R}^{d}$ is a vector of unknown parameters specific to the location $w$. The useful link functions include logit, probit, log-log, etc.
Suppose that the data consist of a random sample $\{(X_{i}, W_{i})\}_{i=1}^{n}$ from the distribution of $(X, W)$ with the sample size of $n$. We can estimate model ((ref)) as binary choice models for the outcomes $\1\{W \le w\}$ under the maximum likelihood framework, where $\1\{\cdot\}$ is the indicator function:
and then we can estimate the conditional distribution function by
Applying the above modeling and estimation procedures on a sequence of locations over the outcome support, the collection of estimation results can characterize the whole conditional distribution.
Our interests lie in the conditional distributional features of bivariate outcomes consisting of continuous and discrete random variables. Practitioners are equipped to grasp the complete picture emerging from the bivariate dependence structure and the influence of specific covariates on various aspects of the variables given the joint conditional distribution. This section outlines the construction of the joint conditional distribution and its associated estimation procedures.
In what follows, we denote by $Y$ a continuous random variable with the support $\mathcal{Y} \subset \mathbb{R}$ and $Z$ a discrete random variable with the finite support $\mathcal{Z} \subset \mathbb{R}$. Let $X$ be a $d_{x} \times 1$ vector of covariates with its support $\mathcal{X} \subset \mathbb{R}^{d_{x}}$. We define $F_{Y|X,Z}$ and $F_{Z|X}$ as the conditional distributions of $Y$ given $X, Z$ and $Z$ given $X$, respectively. Then, we can write the joint distribution function of $(Y, Z)$ conditional on $X$, using the factorization formulation, as follows: for $(x, y, z) \in \mathcal{X}{\times}\mathcal{Y} {\times} \mathcal{Z}$,
We estimate the conditional distributions $F_{Y|X,Z}$ and $F_{Z|X}$ by applying the DR method separately and then obtain the joint conditional distribution as in equation ((ref)). The DR approach fits a parametric linear-index model targeting an arbitrary location of the outcome. The collection of estimation results over outcome locations can characterize the conditional distribution. More specifically, letting $\Lambda: \mathbb{R} \to [0,1]$ be a known link function, we model the conditional distribution function as, for $(x, y, z) \in \mathcal{X}{\times}\mathcal{Y} {\times} \mathcal{Z}$.
where $P_1: \mathcal{X}\times\mathcal{Z}\rightarrow \mathbb{R}^{d_1}$ and $P_2: \mathcal{X}\rightarrow\mathbb{R}^{d_2}$ are two transformations, $\alpha(y)\in\mathbb{R}^{d_1}$ and $\beta(z)\in\mathbb{R}^{d_2}$ are two vectors of unknown parameters. Those unknown parameters are specific to the points of interest, $y$ or $z$, which can be regarded as pseudo-parameters to characterize the conditional distribution at those points, as discussed in the following subsection. By setting the link function as the normal or logistic distribution function, we can consider the models as probit or logit models, respectively. For each outcome, while one can select a different link function $\Lambda(\cdot)$, we use the same notation for simplicity.
There are several advantages of using the DR method to estimate the conditional joint distribution in ((ref)). First, since DR is a local parametric regression, it is easy to implement and computationally fast, even when the number of covariates is moderately large. Second, it characterizes the conditional distributions by collecting regression results over the supports $\mathcal{Y}$ and $\mathcal{Z}$. Thus, the proposed method naturally encapsulates the dependence without global parametric assumptions such as the parametric copula structure. Third, the transformation $P_i$ allows for a flexible enough effect of covariates. For sufficiently rich transformation, one can approximate the conditional distribution function arbitrarily well without extra concern about the choice of the link function. Lastly, the outcome variable of interest can be discrete, continuous, or mixed distributions. This extends some existing works of multidimensional distributional regression klein2022multivariate generally focusing on continuous distributions.
\paragraph{Example} As an illustration, we consider automobile insurance. In an actuarial study, researchers can often observe the number of claims $Z$, the average severity $Y$ and some covariates $X$ for individual policyholders. Insurance companies face claim losses from each individual policyholder and the fixed overhead cost of each claim, denoted by $k >0$. Then, the aggregate claim amount and the total cost of a policyholder are expressed as \[ S := Y \cdot Z \ \ \ \mathrm{and} \ \ \ C :=Y \cdot Z +k \cdot Z. \] The conditional distribution functions of the aggregate claim $F_{S|X}(s|x)$ and of the total cost $F_{C|X}(c|x)$ can be written as
and
For the purpose of risk management, we can consider a risk measure as a transformation of the distribution function to a scalar value. For instance, we can consider a Value-at-Risk (VaR) measure conditional on policyholders' attributes $x \in \mathcal{X}$, given by
This measure is used to estimate the amount of total cost given policyholders' information at a tail event taking place with probability $(1-\tau)$. Similarly, we can consider a Expected Shortfall (ES) given by
This measure is used to evaluate the expected loss on a portfolio in the worst $100\tau\%$ of cases.
Suppose that the data consist of a random sample $\{(X_{i}, Y_{i}, Z_{i})\}_{i=1}^{n}$ from the distribution of $(X, Y, Z)$ with the sample size of $n$. We can consider ((ref)) as models that account for the probability of the events $\{Y \le y\}$ and $\{Z \le z\}$ conditional on the covariates. Thus, we can estimate the models as binary choice models for the outcomes $\1\{Y \le y\}$ and $\1\{Z \le z\}$ under the maximum likelihood framework. More specifically, the estimators are defined as the maximizers of the log-likelihood functions,
where
Using the maximum likelihood estimators, we can estimate the conditional distributions,
for any $(x, y, z) \in \mathcal{X}{\times}\mathcal{Y}{\times}\mathcal{Z}$.
In practice, for the discrete variable $Z$ with support $\mathcal{Z}:=\{z^{(1)}, z^{(2)},\ldots, z^{(L)}\}$, we estimate $\widehat{F}_{Z|X}(z^{(l)}|x)$ for all $z^{(l)}\in\mathcal{Z}, l=1,2,\ldots,L$ to construct the estimator $\widehat{F}_{Z|X}$, as illustrated in Figure (ref)(a). And for the continuous variable $Y$, one can estimate $\widehat{F}_{Y|X, Z}(y^{(j)}|x, z)$ for sufficiently many discrete points $y^{(j)}\in\mathcal{Y}, j=1,2,\ldots,K$ to construct the estimator $\widehat{F}_{Y|X, Z}$, as illustrated in Figure (ref)(b). Computationally, we are estimating $K+L$ local binary regressions in total. Our formulation brings computational benefits from the alternative formulation of DR in meier2020multivariate building the estimator over a grid which requires $KL$ local optimizations.
One important property that characterizes $\widehat{F}_{Y|X, Z}(y|x, z)$ is monotonicity, i.e., conditional distribution functions are non-decreasing by definition. Yet, the estimated distribution functions in ((ref)) do not necessarily satisfy monotonicity in finite samples. We can monotonize the conditional distribution estimators using the rearrangement method proposed by chernozhukov2009improving. The rearranging procedure could yield finite-sample improvement chetverikov2018econometrics and allow for a straightforward application of the functional delta method when we transform the estimated distributions by Hadamard differentiable maps.
The introduced estimation procedures can be applied directly to study any variables on a finite support. However, if the variable of interest has an infinite support, such as the interval $[0, \infty)$ or $(-\infty, \infty)$, the DR approach can be augmented with extreme value theory to provide reliable estimates and insights. The extreme value theory offers methodology for studying the tail behavior of the variable, which allows for extrapolation beyond the range of the available data embrechts2013modelling. Specifically, to adapt the DR approach for variables with infinite support, we first apply the standard DR approach to the finite support of the variable, which is extracted from the dataset, to construct the conditional distribution. Then, to extrapolate the extreme tail, we fit a generalized extreme value distribution based on several conditional CDF values on the tail of the estimated conditional distribution by the method of moments. Finally, we can obtain the conditional distribution on the whole support by combining the conditional distribution on the finite support with the fitted extreme value distribution on the extreme tails.
In this section, we first provide the functional central limit theorems for the estimators of the conditional distribution functions and their transformations. Then, we introduce the exchangeable bootstrap for our estimators and establish its validity for practical inference. Detailed proof of all theoretical results is provided in Appendix (ref). In what follows, let $\|\cdot\|$ be the Euclidean norm for vectors and we denote by $\ell^{\infty}(T)$ the collection of all bounded functions defined on set $T$.
As the population counterpart of the log likelihoods, we define $\ell_{y}(\cdot):= \mathbb{E}[\ell_{i, y}(\cdot)]$ and $\ell_{z}(\cdot) := \mathbb{E}[\ell_{i, z}(\cdot)]$. Then, the true parameters $\alpha_{0}(y)$ and $\beta_{0}(z)$ are defined as the solution to the following maximization problems,
We denote the second derivative of the population log likelihood evaluated at the true parameters by $H_{0, y} := \nabla^2 \ell_{y}\big(\alpha_{0}(y)\big), $ and $H_{0, z} := \nabla^2 \ell_{z} \big ( \beta_{0}(z) \big)$. We group the true parameters into a vector as well as the estimators, by defining
and let $\Theta$ denote the parameter space.\footnote{ The parameter space can be defined for each $(y, z) \in \mathcal{Y}{\times}\mathcal{Z}$, while we suppress the dependency for notational simplicity. } Also, for $\theta:=(\alpha', \beta')\in \Theta$, we introduce a vector of the first derivatives of functions in ((ref)) as
where $R(u):=\lambda(u)/\big\{\Lambda(u)[1-\Lambda(u)]\big\}$.
To obtain the asymptotic results, the following assumptions are imposed.
Assumptions:
Assumption A1 is imposed in the research conducted by chernozhukov2013inference. The bounded supports of covariates ensures that $\mathbb{E}\|X_{i}\|^2 < \infty$. It's worth noting that the compactness assumption on the supports is essential for the uniform valid statistical inference over the entire state space, but it is not necessary for estimation purposes. Assumption A2 ensures that standard optimization procedures based on derivatives can easily obtain the maximum likelihood estimators. A similar condition is assumed in the research by chernozhukov2013inference and both the logit and probit links satisfy this condition. Assumption A3 guarantees the existence of the true parameters. Even when the model in ((ref)) is miss-specified, we can consider the true parameters as pseudo-parameters satisfying the first-order conditions, $\nabla \ell_{y} \big ( \alpha_{0}(y) \big) =0$ and $\nabla \ell_{z} \big ( \beta_{0}(z) \big) = 0$, under assumptions A2 and A3, and thus the estimators of the parameters can be interpreted under the quasi-likelihood framework for each $y \in \mathcal{Y}$ and $z \in \mathcal{Z}$ Huber1967, White1982. Assumption A4 is required to ensure that the information matrices are invertible over the supports. Assumption A5 is required to obtain the limit process of our estimators over the supports for statistical inference.
Under the assumptions above, the proposition below provides the limit process of the estimators $ \hat{\theta}(y, z) $ over $\mathcal{Y} \times \mathcal{Z}$.
In the proof of Proposition (ref), we use the convexity property in Assumption A2 to obtain the limiting processes over the support $\mathcal{Y}$, following the argument used for quantile regression pollard1991asymptotics, kato2009asymptotics. The result in Proposition (ref) shows that the covariance function exhibits the sandwich form as the covariance matrix is obtained under the quasi-likelihood framework.
The distribution function estimators, $\widehat{F}_{Y|X, Z}$ and $\widehat{F}_{Z|X}$, are a transformation of the estimator $\hat{\theta}(\cdot)$ as in ((ref)). Let $\mathbb{D}:= \ell^{\infty}(\mathcal{Y})^{d_1} {\times} \ell^{\infty}(\mathcal{Z})^{d_2} $ and define $ \phi( \hat{\theta}) = \big ( \widehat{F}_{Y|X, Z}, \widehat{F}_{Z|X} \big ) $, where the map $\phi: \mathbb{D}_{\phi}\subset \mathbb{D}\mapsto \mathbb{S}_{\phi}$, given by
It can be shown that the map $\phi$ is Hadamard differentiable at $\theta \in \mathbb{D}_{\phi}$ tangentially to $\mathbb{D}$ with the Hadamard derivative $(a, b) \mapsto \phi_{\theta_{0}(\cdot)}'(a, b)$ is given by
The theorem below shows the joint asymptotic distribution of the distribution function estimators, applying the functional delta method with the Hadamard derivative in the above display. Furthermore, we can easily derive the asymptotic distribution of the estimator of distributional characteristics, such as the Value-at-Risk conditional on covariates and distributional features of $Y$ and $Z$ after some transformation if the distributional characteristics are obtained through Hadamard differentiable maps.
The limiting processes presented in the above proposition and theorem depend on unknown nuisance parameters and may complicate inference in finite samples. The subsequent subsection introduces the bootstrap scheme and reveals its validity.
To deal with the issue of nonpivotal limit processes, we consider a resampling method called the exchangeable bootstrap praestgaard1993Exchangeably, van1996weak. This resampling scheme consistently estimates limit laws of relevant empirical distributions and thus, using the functional delta method, consistently estimates the limit process of the estimator.
For the resampling scheme, we introduce a vector of random weights $(W_{1}, \dots, W_{n})$. To establish the validity of the bootstrap, we assume that the random weights satisfy the following conditions.
Assumption B. Let $(W_{1}, \dots, W_{n})$ be $n$ scalar, nonnegative random variables, which are identically distributed, independent of the original sample, and satisfy the following conditions: for some $\epsilon>0$,
As van1996weak explain, this resampling scheme encompasses a variety of bootstrap methods, such as the empirical bootstrap, subsampling, wild bootstrap and so on. These conditions are employed by chernozhukov2013inference for the inference of counterfactual distributions.
Given the random weights, we obtain the bootstrap estimator
by maximizing the log likelihoods:
Then, we can obtain the bootstrap counterparts of the conditional distribution estimators:
as well as their transformation $ \nu \big( \widehat{F}_{Y|X, Z}^{\ast}, \widehat{F}_{Z|X}^{\ast} \big) $.
For the validity of the resampling method explained above, we need to introduce the notion of conditional weak convergence in probability, following van1996weak. For some normed space $\mathbb{Q}$, let $BL_{1}(\mathbb{Q})$ denote the space of all Lipschitz continuous functions from $\mathbb{Q}$ to $[-1,1]$. Given the original sample $\{(X_{i}, Y_{i}, Z_{i})\}_{i=1}^{n}$, consider a random element $B_{n}^{\ast}:= g(\{(X_{i}, Y_{i}, Z_{i})\}_{i=1}^{n}, \{W_{i}\}_{i=1}^{n})$ as a function of the original sample and the random weight vector generating the bootstrap draw. The bootstrap law of $B_{n}^{\ast}$ is said to consistently estimate the law of some tight random element $B$ or $B_{n} \rightsquigarrow^{p} B$ if
where $\mathbb{E}_{n}$ is the expectation with respect to $\{W_{i}\}_{i=1}^{n}$ conditional on the original sample.
In the theorem provided below, we first show that the exchangeable bootstrap provides a method to consistently estimate the limit process of a pair of conditional distributions. Additionally, we show that the limit process of the Hadamard differentiable transform can be estimated using the functional delta method.
In practice, we monotonize the bootstrap counterparts of the conditional distribution estimators, using the rearrangement method proposed by chernozhukov2009improving.
This section presents Monte Carlo simulation results to reveal the finite-sample properties of the proposed method, compared with the existing methods.
Let $F_{Z}$ and $F_{Y}$ be two parametric distributions left to be specified for $Z$ and $Y$, respectively, and let $g_{Z}$ and $g_{Y}$ be two proper link functions for modeling the conditional means of $Z$ and $Y$ under GLM framework, respectively. We compare our method with the following two popular parametric models: the hierarchical model of garrido2016generalized
and the Gaussian copula regression model of czado2012mixed
where $C_{\eta}:[0,1]^{2} \to [0,1]$ is the Gaussian copula with the correlation parameter $\eta$.
We consider the three different data generating processes (DGPs) below for obtaining the samples of $\big\{(X_{i}, Y_{i}, Z_{i})\big\}_{i=1}^{n}$ with the sample size $n=2,000$. We set $\alpha=(0.5,1,1,1)'$ for all DGPs, while we consider two cases of $\beta$ for each DGP in order to consider that the probability of $Z$ taking zero differs across the cases, as follows:
We select parameter values to ensure that $\Pr(Z=0)$ is roughly between 0.20-0.25 and 0.65-0.70 under Case 1 and Case 2, respectively. In all DGPs, we assume that $Y$ equals 0 if $Z$ is 0, thus the outcome $Y$ follows a mixed distribution. We use regressors $X=\big(1, X^{(1)},X^{(2)},X^{(3)}\big)'$ with $X^{(j)}$ randomly generated from the standard uniform distribution for $j=1, 2, 3$.
For all estimation models, we first model the two univariate conditional distributions of $Z$ and $Y|Y>0$, given $P(Y=0|X)\equiv P(Z=0|X)$, the mixed distribution of $Y$ and the joint distribution of $(Y,Z)$ can be obtained. For comparisons under DGPs 1 and 2, the true marginal distribution families are assumed for $F_Z$ and $F_Y$ of the two parametric models. Thus, the hierarchical model in ((ref)) and the copula model in ((ref)) are correctly specified under DGPs 1 and 2, respectively. On the other hand, the DR is miss-specified in all DGPs. Taking miss-specification into account, we consider a transformation that includes pairwise products of regressors in addition to the original regressors for all estimation models. For the DR approach, the logit link function is applied and the discretization points of the support of $Y$ are chosen as empirical quantiles of $\{Y_{i}\}_{i=1}^{n}$ for probabilities 1%, 2%, $\ldots$, 100%.
We compare the performance of our DR method with the two competing models on estimating the conditional mean, conditional standard deviation, 95% ES, and 0.95th conditional quantiles of $C=Y\cdot Z+Z$ given $X=x$, which are denoted by $\mathbb{E}(C|x)$, $Std(C|x)$, $ES_{0.95}(C|x)$ and $Q_{0.95}(C|x)$, respectively. For the DR method, we obtain the conditional CDF $F_{C|X}(c|x)$ across a series of discrete points over the support of $C$, which allows us to estimate these measures numerically. Specifically, upon obtaining the sequence of conditional CDF values, we could employ the inverse transformation method to generate a $10,000$ random samples, and approximate the desired measures by their sample counterparts. We consider values of regressors $x = (1, x_{1}, 0.5, 0.5)'$ with $x_1$ taking 0.25, 0.50, or 0.75 for comparison. Under each DGP, we present the estimated errors measured by bias and mean square of errors (MSE) based on 1,000 times Monte Carlo simulations.
\paragraph{Poisson-GB2 Hierarchical DGP}
For this DGP, we assume that $F_Z$ and $F_Y$ follow the Poisson and GB2 distributions, respectively, for both parametric models. Accordingly, the log and identical link functions are adopted for $g_Z$ and $g_Y$, respectively. Thus, under this DGP, the hierarchical model is correctly specified, and the Gaussian copula regression uses the correctly specified margins but the misspecified dependence structure.
All of the simulation results are presented in Table (ref). As expected, the hierarchical model has the best estimation performance, while the DR approach can always provide a comparative performance in most cases. In particular, the DR approach performs better than the hierarchical model on estimating the ES. Both of them consistently outperform the Gaussian copula regression model for all quantities.
\paragraph{Poisson-Gamma Gaussian Copula DGP}
For this DGP, we assume that $F_Z$ and $F_Y$ follow the Poisson and Gamma distributions, respectively, for both parametric models, and the log link function is adopted for both the $g_Z$ and $g_Y$. Therefore, the Gaussian copula model is correctly specified for this simulation exercise, and the hierarchical model is correctly specified in the margin but with a different dependence structure. The comparison results are provided in Table (ref). The hierarchical model gives the best estimate of the mean value, but the DR Method performs comparatively well. The correctly specified Gaussian copula model performs better than the hierarchical model for the standard deviation, ES, and quantile in case 1, while the hierarchical model has better performance in most scenarios in case 2. The DR approach reveals a great advantage in estimating the quantile, and it is reasonably comparable with the better model for other quantities.
\paragraph{Truncated Bivariate Normal DGP}
The distributions of $Z$ and $Y$ in this DGP are not the generally adopted distribution families; it is thus difficult for the parametric models to specify the appropriate distributions directly. Based on the sample information, we specify the Poisson and Gamma distributions as the marginal distributions for both parametric models. As shown in Table (ref), the hierarchical model fails to estimate all quantities properly in case 1, but it can give desirable results in case 2, especially for the mean. The copula model provides the best estimation of the mean and standard deviation in case 1, while the DR approach greatly outperforms the two parametric models on estimating the quantile and ES in most scenarios.
Overall, the simulation results reveal an advantage of our semiparametric approach on estimating the higher order moments, quantiles, and ES. The parametric models always estimate the mean properly in that they are mean-based regression models, while the DR approach is reasonably comparable to the correctly specified model. In practice, the exact distributional characteristics are never known exactly, and our approach provides an estimation procedure for flexibly modeling the joint distribution of multiple random variables conditional on some covariates.
In Appendix (ref), we present additional simulation results. First, the 95th conditional quantile, mean, and standard deviation of $Y$ under each DGP are compared. Besides, another DGP based on the hierarchical model with Negative Binomial and Log Normal distributions are considered. In this section, all comparisons are conducted by looking at three different covariate values. In the appendix, we provide additional comparison results by looking at $1000$ randomly generated covariates from the uniform distributions, which is considered as a cohort. For each DGP, the 95% VaR and ES of $C$ for this cohort are explored in both cases.
To test the empirical application of our method, we analyze motor third-party liability policies from an unknown French insurance company. In non-life insurance, the Collective Risk Model (CRM) has become one of the most crucial decision-making models. With recent developments, different dependent structures between the frequency and severities are accommodated to extend the traditional CRM garrido2016generalized, czado2012mixed. The hierarchical and copula strategies for mixed bivariate modeling are widely applied by considering the claim frequency as the discrete outcome and the individual or the average claim severity as the continuous outcome. As discussed previously, the proposed DR approach can also be applied to study the insurance data without worrying about the model's specifications. For the real data analysis, we compare our method with two popular dependent CRMs: the hierarchical model and the Gaussian copula regression model.
The French Motor Third-Part Liability data we used are publicly available (R-Package CASdatasets). The data comprises 413,169 observations; each consists of a set of characteristics associated with the policyholder and their past claim experience. Specifically, the datasets consist of the number of claims (frequency), individual claim amounts, and several rating factors (listed in Table (ref)) for each policyholder observed mostly in one year. In this application, the discrete and continuous outcomes of interest are the claim frequency $Z$ and the average severity $Y$. The average severity for each policyholder is created by taking an average of the claim amounts over the number of claims made. In addition, we obtain the corresponding aggregate claim amounts $S$ of each policyholder for analysis. We create dummies for the non-ordered categorical variables for model construction, e.g., the car brands and administrative regions, and we have a total of $20$ covariates. In our analysis, all rating factors are considered covariates $X$ when modeling the claim frequency and average severity in different estimation models.
A couple of conclusions can be drawn from a preliminary statistical analysis of the dataset. First, the data is heavily concentrated on zero, with 397,779 (or 96.28%) policyholders not having made any claims at all. For the policyholders who made claims, the distribution of the average severity is very skewed (refer to panel (a) in Figure (ref)), that the mean of the distribution exceeds the 75% quantile and the median is close to its $75\%$ quantile. This phenomenon is well documented in the literature YangLu2020NCEf. In addition, the bimodal shape implies that grasping its behavior fully would be difficult with traditional parametric approaches.
For the hierarchical and Gaussian copula models, we use a GLM with Poisson distribution and log link function for modeling the conditional distribution of the claim frequency $Z$. Given that the empirical distribution of the average severity $Y$ is right-skewed and long-tailed, a GLM with GB2 distribution and identical link function is applied for modeling its conditional distribution in both parametric models. In the proposed DR approach, we use the logit link function and identical transformation for both the frequency and the average severity. The support of $Z$ is $\{0,1,\ldots,4\}$ in this dataset, and the discretization points for estimating the distribution of $Y$ are chosen as the $0.1\%, 0.2\%,\ldots,100\%$ quantiles of positive $\{Y_i\}_{i=1}^{n_1}$.
To run the model comparison, we randomly separate the total dataset into the training set with $n_{1}=$ 300,000 policyholders and the validation set with $n_{2}=$ 113,169 policyholders. In the training set, there are only 12,139 policyholders who had made claims.
We first look at the claim frequency $Z$. Table (ref) reports the observed frequency and the fitted frequency using the estimated models and the chi-square statistics. The fitted frequency is calculated as $n_1\widehat{P}(Z=z)$, where the estimated probabilities $\widehat{P}(Z=z)$ can be obtained from the $\widehat{F}_{Z|X=x_i}$ as follows, $$ \widehat{P}(Z=z)=\frac{1}{n_1}\sum _{i=1}^{n_1}\left(\widehat{F}_{Z|X}(z|x_i)-\widehat{F}_{Z|X}(z-1|x_i)\right). $$ The much smaller test statistic of the proposed method suggests a better performance than the global assumption of the Poisson distribution. More importantly, unlike the parametric approach here, the DR approach can capture the upper tail of the distribution that helps insurers oversee and manage their claim-handling expenses.
Note here $Y$ has a mixed distribution. For the average severity $Y$, we set it at zero where the policyholder did not make a claim, that is, $P(Y=0|X)\equiv P(Z=0|X)$. Thus, for all estimation models, we directly model the univariate conditional distributions of $Z$ and $Y|Y>0$ so as to characterize the conditional joint distribution of $(Y,Z)$. We look at the empirical distribution of positive $Y$ and the estimated distributions, $n_1^{-1}\sum_{i=1}^{n_1}\widehat{F}_{Y|X=x_i,Z>0}$, where $$\widehat{F}_{Y|X,Z>0}(y|x):=\frac{\int_{\mathcal{Z}\setminus\{0\}}\widehat{F}_{Y|X,Z}(y|x,z)d\widehat{F}_{Z|X}(z|x)}{1-\widehat{F}_{Z|X}(0|x)}, $$ by the DR, the copula, and the hierarchical models. In Figure (ref), we plot the histograms created based on samples generated according to the estimated distributions. All of the distributions are extremely long-tailed, with the maximal sample values of the true data, by the DR, hierarchical, and copula models are 210,837, 222,466, 142,466, and 762,466, respectively. So, we truncate the distribution up to 15,000 in all histograms to reduce the visual distraction from the tail. The results show that the DR method fits the true distribution much better, while the hierarchical and copula models completely fail to capture the distribution mode.
In addition, we shall investigate the in-sample and out-of-sample performance on estimating the distribution of the aggregative claim amount $S=Z\cdot Y$, which is a quantity of great interest for any insurer. The conditional distribution of $S$ is provided in the Example above and its unconditional distribution is obtained by averaging over the covariate. To demonstrate the efficiency of our estimated distribution, we shall present the estimated out-of-sample CDF via 300 bootstrap samples, randomly selected through permutations with replacements.
We plot the estimated unconditional distribution of positive $S$, $n_1^{-1}\sum_{i=1}^{n_1}\widehat{F}_{S|X=x_i,Z>0}(s|x_i)$ over its empirical counterpart, the mean of $\1 \{S_{i}\leq s|S_{i}>0\}$ in Figure (ref). First, in panel (a), we plot the point estimate of the unconditional CDF of positive $S$ constructed using the three methods against the empirical in-sample CDF. Based on these results, both the hierarchical and the copula models underestimate the CDF function in the lower tail of the distribution and vice versa in the upper tail. For insurance companies, these models are often used to forecast, i.e., provide guidelines for pricing and risk management of future policyholders. Hence, out-of-sample performances are important. In panels (b)-(d) of Figure (ref), we present the out-of-sample distribution forecast of positive $S$. The results are quite similar to that of the in-sample ones. The DR approach stands as the preferred approach, with only a slight overestimate of the CDF around the median. The results clearly demonstrate the superior performance of our method across the entire distribution and the robustness to the empirical distribution of $S$.
For risk management purposes, we illustrate the use of our method in the analyses of risk factors on the out-of-sample dataset. The driver's age and the gas type are often critical risk factors in insurance ratemaking and risk management, so in this section, we focus our analysis on four different policyholder cohorts separated by these two factors.
First, we explore the dependence between the claim amount and frequency. Based on samples generated via the estimated joint conditional distribution by the DR approach, we compute the Pearson, Kendal's tau, and Spearman's rho correlation coefficients for the four policyholder cohorts, given in Table (ref). All of the correlation coefficients show that there is no significant relationship between these two variables, while this result is misleading in that the probability of a policyholder incurring only one claim is around 0.95 for all cohorts. Furthermore, in Figure (ref), we visualize the joint distributions for each policyholder cohort using boxplots constructed based on the samples. There is clear evidence that the extreme average severity and frequency are negatively associated, which is consistent with our belief that drivers who incur only one claim are more likely to make extremely large claim amounts, while drivers who file several claims are typically involved in minor accidents. Further comparisons among different cohorts show that drivers over 30 are more likely to be involved in severe accidents than those using petrol cars, while there is no significant difference for drivers using petrol and diesel cars.
We assume that the total cost associated with each policyholder is given as $C=Y\cdot Z+200\cdot Z$, where $200$ is assumed for the fixed claim handling expense when a claim arises. As referred to previously, given the estimated conditional distribution, we can assess the change in the VaR and ES across different cohorts of policyholders based on the quantities $VaR_{\tau}(C|x)$ and $ES_{\tau}(C|x)$. Such analysis allows analysts to identify the cost leaders in the portfolio and make adequate risk management adjustments. We demonstrate the out-of-sample risk measures $VaR_{\tau}(C|x)$ and $ES_{\tau}(C|x)$ for $\tau=0.98, 0.99$ computed using our proposed method and that of the competing parametric approaches in Tables (ref) and (ref), respectively. In each table, the first two rows corresponding to “unconditional" are results obtained based on all of the out-of-sample policyholders, and in the remaining rows, we present the results of four different policyholder cohorts.
We present the point estimates and 95% confidence intervals of VaR and ES constructed through 300 bootstrap samples for all scenarios in both parametric models. As the support of the claim severity distribution is typically assumed to be infinite in insurance, the proposed approach extrapolates the extreme tail of the conditional distribution by fitting a generalized extreme value distribution, as discussed in Section (ref). The extreme value theory offers a channel for extrapolation outside the range of the available data, as demonstrated in the column labeled "DR(EVD)." However, using bootstrap inference towards the boundary of the support raises challenging theoretical problems. Therefore, in this case, only the point estimates for each scenario are provided, and inference is left for future study. On the other hand, the DR approach proposed on compact support captures the tail behavior of the distribution comparable to the DR with extreme value extrapolation. Additionally, the compact support argument provides the convenience of standard statistical inference via bootstrap samples.
First, for $VaR_{\tau}(C|x)$, the results reveal that the proposed approach provides much more accurate point estimates than the hierarchical and copula models. We reiterate that the point estimates obtained by extrapolating the conditional distribution using the extreme value distribution are consistent with the results obtained by using the DR approach alone. Additionally, the DR approach without extrapolation provides narrow confidence intervals that accurately capture the true out-of-sample results in all scenarios. Both parametric models, especially the copula model, tend to underestimate the 98% VaR and overestimate the 99% VaR for all cohorts. Empirically, we observe the younger cohort's higher risk profile, and the use of diesel intensifies it. The estimated VaR across the three approaches confirms this increase in risk profile. For $ES_{\tau}(C|x)$, our findings from Table (ref) show that in all cases, the proposed approach provides more accurate point estimates that are closer to the empirical results than the hierarchical model, while the copula model completely fails to estimate the ES properly. The ES is a risk measure that is more sensitive to extreme values compared to VaR. Therefore, the DR approach augmented with the extreme value distribution provides more precise point estimates compared to the DR approach without extrapolation. Furthermore, our results confirm that the risk profile is higher for younger cohorts and the use of diesel intensifies the risk.
This paper proposes a semiparametric method based on the distribution regression approach. The advantage of using this proposed method is three-fold. First, we avoid imposing too rigid parametric assumptions, which makes our approach robust for analyzing real data. Secondly, the covariates are incorporated to influence the whole distribution instead of only affecting the distributions' location parameters. Finally, by including the discrete outcome as a covariate in the conditional distribution of the continuous outcome, our model captures intricate dependence structures between the two outcomes. While our analysis in this paper focuses on bivariate modeling of one discrete and one continuous outcome, the method can be easily extended to allow any random variables, including mixed distributions. The simulation examples under different scenarios demonstrate the robustness of our method. The empirical study shows that our method can extract interesting features of a motor insurance portfolio critical to pricing and risk management in real data applications.
\setstretch{1}