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.
99,233 characters · 19 sections · 78 citation commands
Deep Quantile and Deep Composite Model Regression
In actuarial modeling, we typically use regression models to estimate expected values of claim sizes exploiting systematic effects of features or covariates. These estimated expected values are used for insurance pricing, claims forecasting and claims reserving. The most commonly used regression model is the generalized linear model (GLM) with the exponential dispersion family (EDF) as the underlying distributional model. Parameter estimation within this EDF-GLM framework is performed by maximum likelihood estimation (MLE). In this setup, MLE is equivalent to minimizing the corresponding deviance loss function, the latter being a strictly consistent scoring function for mean. Strictly consistent scoring functions facilitate M-estimation and forecast evaluation; further details are provided below. For actuarial problems, there are three main issues with this approach:
Our paper tackles the above points; we first discuss the related literature and then we are going to describe our proposal. To deal with the problem of the lack of any simple off-the-shelf distribution that fits the entire range of the data, i.e., the body and tail of the data, one often either uses a composite (splicing) model or a mixture model. These two approaches can be fitted with the expectation-maximization (EM) algorithm. Composite models are studied in Cooray--Ananda CoorayAnanda, Scollnik Scollnik, Pigeon--Denuit PigeonDenuit, Gr\"un--Miljkovic Grun and Parodi Parodi. These papers do not consider a regression structure and covariates. The only actuarial publications that we are aware of considering composite regression models are Gan--Valdez GanValdez and Laudag\'e et al. Laudage. Both proposals choose very specific distributional assumptions above and below a given threshold (splicing point) to obtain analytical tractability. In our proposal, we do not fix the threshold itself, but a quantile level instead. This allows us for direct model fitting below and above that quantile, and, moreover, we can have a flexible regression structure for both the body and the tail of the data. For mixture models, not considered in this paper, we refer to Fung et al. Fung2021 and the literature therein. Both approaches, composite and mixture models, typically use the EM algorithm for model fitting. Our approach is based on the (simpler) stochastic gradient descent (SGD) algorithm.
Focusing on point (2) from above, there are several ways in extending (or modifying) a GLM. These include, e.g., generalized additive models (GAMs), regression trees and tree boosting, or feed-forward neural (FN) network regression models; we refer to Hastie et al. HTF. Here, we focus on FN network regression models as a straightforward and powerful extension of GLMs. Moreover, we see that an FN network architecture can easily be designed so that it can simultaneously fulfill different estimation tasks. This is a crucial property that we need in our derivations, as we jointly estimate different regression models for the body and tail of the data.
Concerning point (3) from above, quantile regression has gained quite some popularity in the machine learning community to quantify uncertainty, see Meinshausen Meinshausen and Takeuchi et al. Takeuchi. Quantile regression has been introduced by Koenker--Bassett KoenkerBassett, and it is widely used in statistical modeling these days; for a recent monograph we refer to Uribe--Guillen UribeGuillen. Our proposal extends the approach of Richman Richman3 who combines joint estimation of a quantile and the mean within an FN network architecture, see Listing 8 in Richman Richman3.
We consider two different, though closely related problems in this paper. On the one hand, we extend the work of Richman Richman3 by proposing an FN network architecture that allows for consistent multiple quantile regression respecting monotonicity of the quantiles at different levels. On the other hand, this FN network architecture is at the basis of the extension of the joint quantile and ES regression considered in Guillen et al. Guillen2021. The main building blocks in our estimation procedure are strictly consistent scoring functions for our target functionals. These are functions of the estimated model prediction and the response variable which are minimised in expectation by the correctly specified model. Strictly consistent scoring functions are at the core of consistent M-estimation (where M stands for minimisation) in a regression setup, see Dimitriadis et al. DimitriadisFisslerZiegel2020. If a target functional admits a strictly consistent scoring function, it is called elicitable. The elicitability of the mean and the quantile thus allow for mean and quantile regression, using the squared loss (or a Bregman divergence) or the pinball loss (or a generalized piecewise linear loss), respectively. In model selection and, more generally, forecast ranking, Gneiting Gneiting and Gneiting--Raftery GneitingRaftery advocate for the usage of strictly consistent scoring functions since they honour honest and truthful forecasts.
In a first step, we perform multiple quantile regression using the sum of pinball losses as strictly consistent scoring functions. Since quantiles are monotone in their probability levels, we modify the FN network architecture of Richman Richman3 such that we can guarantee this monotonicity.
In a second step, we extend the work of Guillen et al. Guillen2021 in several ways. In risk management, besides knowing quantiles, one is also interested in estimating ES, which reflects the average outcome beyond the quantile. Since ES is not elicitable Gneiting, stand-alone ES regression is generally not possible. Fissler--Ziegel FisslerZiegel2016 showed that the pair of ES together with the quantile at the same probability level is elicitable. This paved the way to joint quantile--ES regression, see Dimitriadis--Bayer DimitriadisBayer and Guillen et al. Guillen2021 for (G)LM approaches. We first extend the result of Fissler--Ziegel FisslerZiegel2016, showing the elicitability of the composite triplet consisting of a quantile together with the lower and upper ES at the same probability level, and characterizing the entire class of strictly consistent scoring functions for it. Hence, we are in the position to estimate a full composite regression model in a one-step procedure. This is in contrast to the two-step estimation approach of Barendse Barendse2020, first estimating the quantile and then in a second step the ES below and above the quantile. Such a two-step estimation approach becomes problematic in the presence of common parameters in the three components of the regression model. Second, we motivate the particular choice of the strictly consistent scoring function in a data-driven manner, striving for efficient parameter estimation. To this end, we use optimality results known for mean regression which are akin to the optimality results derived in Dimitriadis et al. DimitriadisFisslerZiegel2020 for the pair of the quantile and ES. The third novelty we provide is the FN network architecture that respects the necessary monotonicity property of the composite triplet in an estimation context. This architecture can directly be fitted using the stochastic gradient descent (SGD) algorithm. This fitting approach does neither require a two-step fitting approach nor the EM algorithm as it is used in related problems. Moreover, our fitting turns out to be robust, and we do not encounter the stability issues as in the two-step approaches and the EM algorithm.
{\bf Organization of this manuscript.} In Section (ref) we review the concepts of strictly consistent scoring functions used in estimation and forecast evaluation. We recall known characterization results of strictly consistent scores for the mean and for the quantile. We extend the results of Fissler--Ziegel FisslerZiegel2016 introducing the class of strictly consistent scoring functions for the composite triplet of the quantile, lower ES and upper ES at the same probability level. In Section (ref) we discuss how these quantities can be estimated within a neural network regression framework, and we also discuss asymptotically efficient choices of scoring functions. In Section (ref) we give a real data example which demonstrates the suitability of our proposal. Section (ref) concludes and provides an outlook to related open problems.
We start from the decision-theoretic approach of forecast evaluation developed by Gneiting Gneiting and Gneiting--Raftery GneitingRaftery, also used for M-estimation in Dimitriadis et al. DimitriadisFisslerZiegel2020. This provides us with suitable choices of scoring functions for model fitting, model selection, and forecast evaluation. Denote by ${\cal F}$ the class of considered distribution functions $F \in {\cal F}$. Let ${\mathbb Y} \subseteq {\mathbb R}$ be an interval, a halfline or the whole real line ${\mathbb R}$, such that the support of each $F\in {\cal F}$ is contained in ${\mathbb Y}$. Choose an action space ${\mathbb A} \subseteq {\mathbb R}$ from which we select actions $a \in {\mathbb A}$ to estimate a statistics $A(F) \in {\mathbb A}$ of $F \in {\cal F}$. That is, we have a functional
that we try to estimate. Commonly used functionals are the mean functional $A(F)={\mathbb E}_F[Y]$ for $Y\sim F \in {\cal F}$ and quantiles. Given a probability level $\tau \in (0,1)$, the $\tau$-quantile of $F \in {\cal F}$ is given by the (left-continuous) generalized inverse
To receive an intuitive understanding we usually speak about a functional (ref) that attains a single value $A(F)$ in the action space ${\mathbb A}$. Often this functional is obtained by an optimization (M-estimator), or by finding the roots of a given function (Z-estimator). Since, in general, we do not want to assume uniqueness of such a solution, we should understand this functional as a set-valued map
where ${\cal P}({{\mathbb A}})$ is the power set of ${\mathbb A}$. E.g., the set-valued $\tau$-quantile of a distribution $F$ is given by
where $F(t-) = \lim_{x\uparrow t}F(x)$. This defines a closed interval and its lower endpoint corresponds to the left-continuous generalized inverse $F^{-1}(\tau)$ given in (ref). To keep notation simple, we identify $\widehat{a}$ and $A(F)$ if $A(F)=\{\widehat{a}\}$ is a singleton.
In order to evaluate the accuracy of actions $a$ for the statistics $A(F)$ (for unknown $F$) we consider a scoring function (also called loss function)
This describes the loss of an action $a \in {\mathbb A}$ if a realization $y $ of $Y \sim F$ materializes. To incentivize truthful forecasts, Gneiting Gneiting advocates that this scoring function $L$ should be strictly consistent for the functional $F\mapsto A(F)$ of interest.
Gneiting Gneiting shows that strictly consistent scoring functions are linked to proper scoring rules, which serve to evaluate probabilistic forecasts (i.e. actions $a$ taking the form of probability distributions) and which are discussed in detail in Gneiting--Raftery GneitingRaftery. On the estimation side, we also need (strict) consistency of scoring functions to obtain consistency of M-estimators, see Chapter 5 in Van der Vaart VanderVaart2 and Dimitriadis et al. DimitriadisFisslerZiegel2020. This raises the question of which functionals $A: {\cal F} \to {\mathbb A}$ admit strictly consistent scoring functions $L$, i.e., are elicitable.
The general question of elicitability is studied and discussed in detail in Gneiting Gneiting. For instance, the mean functional $A(F)={\mathbb E}_F[Y]$ and the $\tau$-quantile $q_\tau$ defined in (ref) are elicitable; see Theorems (ref) and (ref), below. We start with useful properties of scoring functions $L$ that are going to be used in the next theorems. Assume ${\mathbb Y}={\mathbb A} \subseteq {\mathbb R}$ is an interval with non-empty interior. We set:
These regularity conditions can often be weakened. Especially the positivity in (S0) can be relaxed. However, it is a nice property to have in optimizations as it gives us a natural lower bound to the minimal score that can be achieved.
The following result goes back to Savage Savage.
The maps $L_{\phi}$ defined in (ref) are called {\it Bregman divergences} and, basically, the previous theorem says that all strictly consistent scoring functions for the mean functional are Bregman divergences. The most prominent Bregman divergence arises by setting $\phi(y) = y^2$, yielding the squared loss $L(y;a) = (y-a)^2$. On the other hand, the variance functional $\mathbb V_F(Y)$ is not elicitable (Gneiting Gneiting and Osband Osband1985), i.e., there does not exist any strictly consistent scoring function for estimating the variance according to a minimization (ref).
The following elicitability results for quantiles (ref) originates from Thomson Thomson and Saerens Saerens.
Members of the class (ref) are called generalized piecewise linear losses, see Gneiting Gneiting2011b. Basically, it is this theorem that allows us to consider quantile regression, introduced by Koenker--Bassett KoenkerBassett, as it tells us that quantiles can be estimated from strictly consistent scoring functions. This is going to be outlined below. There is still the freedom of the choice of the strictly increasing function $g$. The most commonly used choice is the identity function $g(y)=y$. This then provides us with the so-called {\it pinball loss} for given probability level $\tau \in (0,1)$
In view of Theorem (ref), the choice of pinball loss requires that ${\cal F}$ contains only distributions with a finite first moment.
We introduce scoring functions closely related to the pinball loss (ref) where the positivity assumption (S0) has been relaxed:
for $y,a\in{\mathbb R}$ and for $\tau \in (0,1)$. Note that $S^+_\tau(y;a) = S^-_\tau(y;a) + y$, moreover, it holds for the pinball loss
An immediate consequence from Theorem (ref) is the following corollary.
Observe that the three functions $S^-_\tau$ $S^+_\tau$ and $L_\tau$ only differ in terms that do not depend on $a$, and therefore the argument of these minimizations is the same. The pinball loss $L_\tau(y;a)$ has the advantage to satisfy condition (S0).
From the previous section we know that we can estimate $\tau$-quantiles (ref) by minimizing (ref) under the strictly consistent pinball loss (ref). In actuarial science we are often interested in also considering the lower ES and the upper ES
of a random variable $Y \sim F$. Since $\mathrm{ES}_\tau^-$ and $\mathrm{ES}_\tau^+$ are law-determined, we can interpret them as functionals (ref). The monotonicity of the generalized inverse $p\mapsto F^{-1}(p)$ immediately yields
The ES, considered as functionals, turn out not to be elicitable on families of distributions $\mathcal F$ which offer the necessary flexibility in most modeling situations; see Gneiting Gneiting and Weber Weber2006. Thus, there is no strictly consistent scoring function suitable for M-estimation of ES. Fissler--Ziegel FisslerZiegel2016 have proved that ${\rm ES}^-_\tau(Y)$ is {\it jointly} elicitable with the $\tau$-quantile $F^{-1}(\tau)$, the latter also being called Value-at-Risk (VaR) in risk management. We are going to extend this result in Theorem (ref), establishing the elicitability of the composite triplet $({\mathrm{ES}}^-_\tau, q_\tau, {\mathrm{ES}}^+_\tau)$. We first need the following lemma.
Proof. This follows directly along the lines of Lemmas 2.3 and 3.3 in Embrechts--Wang EmbrechtsWang2015. {\scriptsize $\Box$}
The intuitive interpretation is that the score in (ref) is a combination of a Bregman divergence (second line) and a generalized piecewise linear loss (first line in combination with $S_{\tau}^-$ and $S_{\tau}^+$). This structure will be exploited in the proof.
Proof of Theorem (ref). First, fix some $v\in{\mathbb Y}$. The map \[ (y;e^-,e^+)\mapsto L(y;e^-,v,e^+) = \left\langle \nabla \Phi(e^-, e^+),
\right\rangle -\Phi(e^-, e^+) + b_{v}(y), \] where the remainder $b_v(y)$ does not depend on $(e^-, e^+)$, is a Bregman divergence. Hence, if $\Phi$ is (strictly) convex, it is (strictly) ${\mathcal F}$-consistent for the functional \[ F\mapsto \big(-\tfrac{1}{\tau} \mathbb E_F\big[ S^-_\tau(Y;v) ], \tfrac{1}{1-\tau} \mathbb E_F\big[ S^+_\tau(Y;v) ]\big). \] Second, for fixed $(e^-,e^+)\in {\mathbb Y}^2$, the map \[ (y;v)\mapsto L(y;e^-,v,e^+) = (\mathds{1}_{\{y\le v\}} - \tau)G_{e^-, e^+}(v) - \mathds{1}_{\{y\le v\}} G_{e^-, e^+}(y) + b_{e^-,e^+}(y), \] where the remainder $b_{e^-,e^+}(y)$ does not depend on $v$, is a generalized piecewise linear loss (ref) not necessarily satisfying the positivity (S0). Hence, if $G_{e^-, e^+}$ is (strictly) increasing, it is (strictly) ${\mathcal F}$-consistent for $q_\tau$.
Combining these two observations yields the (strict) ${\mathcal F}$-consistency of $L$.
Finally, $L(y;y,y,y)=0$ can be verified by a direct computation, and the non-negativity of $L$ follows from its consistency and the fact that for a random variable $Y$ which deterministically equals the constant $y\in{\mathbb R}$, it holds that $\mathrm{ES}_\tau^-(Y) = \mathrm{ES}_\tau^+(Y) = y$ and $q_\tau(Y) = \{y\}$. {\scriptsize $\Box$}
Theorem (ref) extends the result by Frongillo--Kash FrongilloKash2020 asserting that an elicitable functional (in our case the $\tau$-quantile) is jointly elicitable with finitely many associated Bayes risks (here corresponding to $\mathrm{ES}_\tau^-$ and $\mathrm{ES}_\tau^+$).
Theorem (ref), below, shows that the scores of the form (ref) are basically the only ${\mathcal F}$-consistent scores for $(\mathrm{ES}_\tau^-, q_\tau, \mathrm{ES}_\tau^+)$. The argument -- which can be found in its entirety in the Supplementary Section (ref) -- uses Osband's principle, which originates from Kent Osband's seminal thesis Osband1985; see Gneiting Gneiting for an intuitive exposition and Fissler--Ziegel FisslerZiegel2016 for a precise technical formulation. It exploits first- and second-order conditions stemming from the optimization problem of (strict) consistency (ref). Hence, we need to impose smoothness conditions on the expected score, which play the role of condition (S2), but are slightly weaker. This is reflected in Assumption (ref), and the fact that we work within a subclass of ${\mathcal F}\subseteq {\mathcal F}_{\rm cont}^\tau$, the class of distribution functions on ${\mathbb R}$ which are continuously differentiable and whose $\tau$-quantiles are singletons. The second kind of conditions are richness conditions on the underlying class of distributions ${\mathcal F}$. On the one hand, the first- and second-order conditions only yield local assertions about the expected scores. Hence, ${\mathcal F}$ needs to be rich enough such that the functional maps surjectively to the action domain considered. Recalling monotonicity condition (ref), this means that we can only provide conditions on action spaces contained in $\{(a_1,a_2,a_3)\in{\mathbb R}^3\colon a_1\le a_2\le a_3\}$. Second, the richness conditions ensure that the functional `varies sufficiently' such that it can be distinguished from any other functional.\footnote{E.g., on the class of symmetric distributions with positive densities, the mean and the median coincide and cannot be distinguished. Such phenomena need to be excluded.} This is reflected in Assumptions (ref) and (ref). Finally, since Osband's principle actually characterizes the gradient of the expected score, one needs to be in the position to integrate this gradient (Assumption (ref)) and to approximate the pointwise values of the score with expectations (Assumption (ref)).
Supplementary Section (ref) provides all technical details. In particular, the proof of Theorem (ref) can be found in Subsection (ref).
There remain the choices of the functions $g$ and $\Phi$ in (ref). Especially, the choice of the strictly convex function $\Phi$ such that $G_{e^-, e^+}$ in (ref) is strictly increasing is not obvious. We discuss the following particularly convenient three choices for given $\tau \in (0,1)$
where $\phi, \phi_+, \phi_-\colon {\mathbb Y} \to {\mathbb R}$ are strictly convex and the (sub-)gradients satisfy $\phi'_+<0$, $\phi'_->0$. Moreover, if we choose $g\colon {\mathbb Y} \to {\mathbb R}$ to be an increasing function (not necessarily strictly increasing), we obtain for $G_{e^-, e^+}$, defined in (ref),
Since $\phi'_+<0$ and $\phi'_->0$, $G_{e^-, e^+}$ is strictly increasing, even if $g$ is constant. This yields the following three types of scores: For $\Phi$ in (ref) we get
Similarly, for $\Phi$ in (ref) we receive
and finally for (ref)
Scores of the form (ref) can directly be interpreted to be the sum of scores for the pair $(q_\tau, \mathrm{ES}_\tau^-)$ as deduced in Fissler--Ziegel FisslerZiegel2016 and for the pair $(q_\tau, \mathrm{ES}_\tau^+)$ as derived in Nolde--Ziegel NoldeZiegel2017. On the other hand, (ref) can be deduced from the sum of a scoring function for $(q_\tau, \mathrm{ES}_\tau^+)$ and for the mean, using the so-called revelation principle. This principle originates from Osband's thesis Osband1985 and it has been made rigorous in Gneiting Gneiting. It asserts that any bijection of an elicitable functional is elicitable and it makes the corresponding strictly consistent scoring functions explicit. In the case of (ref), this bijection is \[ (\mathrm{ES}_\tau^-,q_\tau, \mathrm{ES}_\tau^+) \mapsto (q_\tau, \mathrm{ES}_\tau^+, \tau \mathrm{ES}_\tau^- + (1-\tau)\mathrm{ES}_\tau^+) = (q_\tau, \mathrm{ES}_\tau^+, {\mathbb E}). \] For the scores in (ref), the corresponding bijection reads similar.
The next section shows how to use strictly consistent scoring functions for parameter estimation via M-estimation in a regression context. Here, the (strict) ${\mathcal F}$-consistency of the corresponding score is crucial to obtain the consistency of the M-estimator, i.e., that the M-estimator converges in probability to the true value as the sample size goes to infinity; see Dimitriades et al. DimitriadisFisslerZiegel2020. Knowing that the estimator is consistent, it is of interest to have fast convergence, i.e., an efficient estimator. Under asymptotic normality of the estimator, a more efficient estimator has a strictly smaller asymptotic covariance matrix.\footnote{As usual, we mean this with respect to the Loewner order. That is, a covariance matrix $A$ is strictly smaller than a covariance matrix $B$ of the same dimension if $A\neq B$ and if $B-A$ is positive semi-definit.} In the absence of any explanatory information (that is, when estimating an intercept only model), the choice of the strictly consistent scoring function is immaterial since the intercept estimators under different strictly consistent scores will always coincide on finite samples and correspond to the functional of the empirical distribution function. In a more interesting regression scenario, however, using feature information, the question of efficiency becomes important, this is discussed in our setup in Section (ref), below.
The previous sections have discussed estimation theory of functionals of (unkown) distribution functions $F$. We now lift this framework to a regression context where random variables $Y$ are supported by covariates (feature information). Assume that the random variable $Y$ is established with feature information $\boldsymbol{X} \in {\cal X} \subset \{1\}\times {\mathbb R}^q$, where ${\cal X}$ is the feature space of all potential explanatory variables $\boldsymbol{X}$. We assume for a datum $(Y,\boldsymbol{X})$ that the conditional distribution of $Y$, given $\boldsymbol{X}$, is described by a distribution function $F_{Y|\boldsymbol{X} = \boldsymbol{x}}$, $\boldsymbol{x}\in{\cal X}$, or in short $F_{Y|\boldsymbol{x}}$, and the conditional $\tau$-quantile of this random variable is given by the left-continuous generalized inverse
Note that we label distributions with subscripts $Y|\boldsymbol{x}$, now, to clearly indicate that we are considering the conditional distribution of $Y$, given feature information $\boldsymbol{X} = \boldsymbol{x} \in {\cal X}$. This gives us the existence of a regression function $Q_\tau :{\cal X} \to {\mathbb R}$ such that
Quantile regression tries to determine this regression function $Q_\tau$ from a given function class ${\cal Q}$ based on i.i.d. observations $(Y_i,\boldsymbol{X}_i)$, $1\le i \le n$. The classical approach of Koenker--Bassett KoenkerBassett makes a GLM assumption by postulating the existence of a strictly monotone and smooth link function $h:{\mathbb R} \to {\mathbb R}$ such that
with regression parameter $\boldsymbol{\beta} \in {\mathbb R}^{q+1}$ and $\langle \cdot, \cdot \rangle$ denoting the scalar product in the Euclidean space ${\mathbb R}^{q+1}$. In this case the function class ${\cal Q}$ is parametrized by a regression parameter $\boldsymbol{\beta} \in {\mathbb R}^{q+1}$, that we try to optimally determine from i.i.d. data $(Y_i,\boldsymbol{X}_i)$, $1\le i \le n$. The choice of the link $h$ acts as a hyper-parameter, that is not part of the optimization process.
We can choose a strictly consistent scoring function $L_\tau$ (ref) for the $\tau$-quantile $F_{Y|\boldsymbol{x}}^{-1}(\tau)$ and estimate the regression parameter $\boldsymbol{\beta} \in {\mathbb R}^{q+1}$ by
subject to existence. Typically, we do not know the true distribution function and, therefore, cannot explicitly evaluate the right-hand side of the above optimization. Using an empirical version based on i.i.d. data $(Y_i,\boldsymbol{X}_i)$, $1\le i \le n$, motivates M-estimator
In practice, the GLM structure (ref) is often too restrictive. This emphasizes to use a more flexible function class ${\cal Q}$. FN networks provide the building blocks for such a more flexible function class. This motivates {\it deep quantile regression} which has been introduced to the actuarial literature by Richman Richman3. We replace (ref) by
where $\boldsymbol{z}^{(d:1)}: {\cal X} \to \{1\}\times {\mathbb R}^{r_d}$ is an FN network of depth $d \in {\mathbb N}$, regression parameter $\boldsymbol{\beta} \in {\mathbb R}^{r_d+1}$ and link $h$. The FN network $\boldsymbol{z}^{(d:1)}$ is a composition of $d$ FN layers
with FN layers $\boldsymbol{z}^{(m)}$ for $1\le m \le d$ involving further parameters $\boldsymbol{w}^{(m)}_j \in {\mathbb R}^{r_{m-1}+1}$, and with $r_{m-1}+1$ describing the input dimension to FN layer $\boldsymbol{z}^{(m)}$; for a detailed exhibition of FN networks we refer to Section 7.2 in W\"uthrich--Merz WM2021. Altogether this FN network approach (ref) is parametrized by
Thus, we fix a depth $d \in {\mathbb N}$, FN layer dimensions $r_1,\ldots, r_d \in {\mathbb N}$, the activation functions in the FN layers and link function $h$ (as hyper-parameters), then the function class ${\cal Q}$ is parametrized by $\boldsymbol{\vartheta}$, and the “optimal” member for i.i.d. data $(Y_i,\boldsymbol{X}_i)$, $1\le i \le n$, is found by
On purpose “optimal” has been written in quotation marks. On a finite sample of size $n$ the solution to (ref) will likely in-sample overfit to the learning data ${\cal L}=(Y_i,\boldsymbol{X}_i)_{1\le i \le n}$ because already small networks are fairly flexible. Therefore, this model is usually fit with a SGD algorithm that explores an {\it early stopping rule}, i.e., that selects an estimate $\widehat{\boldsymbol{\vartheta}}_\tau$ that describes the systematic effects in the data ${\cal L}$ and not the noisy part. This fitting approach is state-of-the-art, and it is described in detail in Section 7.2.3 of W\"uthrich--Merz WM2021. Therefore, we will not repeat it here.
Often we do not want to estimate quantiles for only one probability level, but we would like to study quantiles at different levels. For illustrative purposes we consider two probability levels $0<\tau_1 < \tau_2 < 1$, and a generalization to more than two probability levels is straightforward. A naive way is to individually estimate regression functions $\boldsymbol{x} \mapsto Q_{\tau_l}(\boldsymbol{x})=F_{Y|\boldsymbol{x}}^{-1}(\tau_l)$ for $l=1,2$ using (ref) and (ref). We call this a naive approach because these individual estimations may violate the monotonicity property of quantiles, i.e., for all $\boldsymbol{x}$ we require
To simplify this outline, we assume that the random variable $Y$ is positive, a.s., which implies that the generalized inverse $\tau \mapsto F_{Y|\boldsymbol{x}}^{-1}(\tau) > 0$ has a positive range. This motivates the choice of a link function $h$ with positive support ${\mathbb R}_+$. For enforcing the monotonicity of quantiles for different levels (ref), we propose to jointly model these quantiles. For our first proposal we choose two link functions $h$ and $h_+$ being both positively supported. In analogy to (ref), this motivates {\it joint deep quantile regression} for probability levels $\tau_1 <\tau_2$
for a network parameter $\boldsymbol{\vartheta}=(\boldsymbol{w}^{(1)}_1, \ldots, \boldsymbol{w}^{(d)}_{r_d}, \boldsymbol{\beta}_{\tau_1}, \boldsymbol{\beta}_{\tau_2})^\top$. Thus, we choose a common deep FN network $\boldsymbol{z}^{(d:1)}$ that is {\it shared} by both quantiles, and the bigger quantile $Q_{\tau_2}(\boldsymbol{x})$ is modeled by a positive difference $h_+^{-1}\langle \boldsymbol{\beta}_{\tau_2}, \boldsymbol{z}^{(d:1)}(\boldsymbol{x})\rangle \ge 0$ to the smaller quantile $Q_{\tau_1}(\boldsymbol{x})$. We call (ref) an {\it additive approach} with base level $Q_{\tau_1}(\boldsymbol{x})$, and in network jargon we say that the FN network learns a common representation $\boldsymbol{z}_i = \boldsymbol{z}^{(d:1)}(\boldsymbol{x}_i)$ of features $\boldsymbol{x}_i$, $1\le i \le n$, which is then used in the two GLMs
Alternatively, for positive random variables $Y$, a.s., we can choose the upper quantile $Q_{\tau_2}(\boldsymbol{x})$ as base level, and to ensure positivity we can multiplicatively decrease this upper quantile. For this we choose the sigmoid function for $h_\sigma^{-1}=(1+\exp\{-x\})^{-1} \in (0,1)$ which motivates the {\it multiplicative approach} for probability levels $\tau_1 <\tau_2$
Also in this case monotonicity is guaranteed because by assumption $h_\sigma^{-1}\langle \boldsymbol{\beta}_{\tau_1}, \boldsymbol{z}^{(d:1)}(\boldsymbol{x}) \rangle \in (0,1)$.
Since the learned representations $\boldsymbol{z}_i = \boldsymbol{z}^{(d:1)}(\boldsymbol{x}_i)$ need to fit both quantiles simultaneously we need to learn these representations jointly. This motivates the optimization problem under regression assumption (ref) or (ref), respectively, and up to over-fitting (see discussion after (ref))
where the weights $\eta_1,\eta_2 >0$ are chosen such that both quantiles contribute roughly equally to the total score. This is used to ensure that the quality of the estimate of both quantiles is roughly similar.
A deep composite regression model now only requires little changes compared to the deep multiple quantile regression. Again, we assume that $Y>0$, a.s. We then aim at estimating the composite triplet
where we highlight the conditional structure of $Y$, given $\boldsymbol{X}=\boldsymbol{x}$. Thus, in addition to the quantile regression of the previous sections, we aim at estimating regression functions $\boldsymbol{x} \mapsto E_\tau^s(\boldsymbol{x})={\rm ES}^s_\tau(Y|\boldsymbol{x})$ for $s \in \{-,+\}$. Thanks to Theorem (ref), we know that we can jointly estimate the triplet $(E^-_\tau(\boldsymbol{x}), ~Q_{\tau}(\boldsymbol{x}),~E^+_\tau(\boldsymbol{x}))$ using the strictly consistent scoring function (ref). In analogy to (ref), we choose positively supported link functions $h$ and $h_+$. This motivates {\it deep composite model regression}
for network parameter $\boldsymbol{\vartheta}=(\boldsymbol{w}^{(1)}_1, \ldots, \boldsymbol{w}^{(d)}_{r_d}, \boldsymbol{\beta}_{1}, \boldsymbol{\beta}_{2}, \boldsymbol{\beta}_{3})^\top$. Again, we choose a common deep FN network $\boldsymbol{z}^{(d:1)}$ that is shared by the $\tau$-quantile and the lower and upper ES. Alternatively to (ref), we could also use a multiplicative approach similar to (ref). We obtain the following M-estimator
where $L$ is a strictly consistent score given in (ref).
As already discussed at the end of Subsection (ref), in a regression context, the choice of the strictly consistent scoring function influences the (asymptotic) variance of the M-estimator. Therefore, we would like to provide some guidance on how to choose a score of the form (ref) in a data driven manner. For simplicity, we will focus on scores of the form (ref) and (ref), mainly discussing the choice of $\phi$, $\phi_-$ and $\phi_+$. Moreover, motivated by modeling claim sizes, we assume that $Y>0$, a.s., such that ${\mathbb Y} = (0,\infty)$.
To this end, first recall a classical efficiency result in the context of mean regression; see Newey--McFadden NeweyMcFadden1994. If $\mu(\boldsymbol{x}) = {\mathbb E}[Y|\boldsymbol{X}=\boldsymbol{x}]$ is the estimable regression function and $\sigma(\boldsymbol{x})^2 = \mathbb V(Y|\boldsymbol{X}=\boldsymbol{x})$ is the conditional variance, the most efficient Bregman score (ref) should satisfy
for some $c>0$ and for all $\boldsymbol{x}\in {\cal X}$. Since the third line of (ref) is basically a Bregman score for the mean, we use condition (ref) to come up with a choice for $\phi$ in (ref). For $\phi_+$ in (ref) and (ref), we suggest to exploit a similar relation for the truncated variance (recalling that $\mathrm{ES}_\tau^+$ is also a truncated mean). In particular, Theorem 4.3 in Dimitriadis et al. DimitriadisFisslerZiegel2020 suggests
for some $c_+>0$ and for all $\boldsymbol{x}\in {\cal X}$, where $\sigma_\tau^+(\boldsymbol{x})^2=\mathbb V(Y| Y> F^{-1}_{Y|\boldsymbol{x}}(\tau), \ \boldsymbol{X}=\boldsymbol{x})$. Similarly,
for some $c_->0$ and for all $\boldsymbol{x}\in {\cal X}$, where $\sigma_\tau^-(\boldsymbol{x})^2=\mathbb V(Y| Y\le F^{-1}_{Y|\boldsymbol{x}}(\tau), \ \boldsymbol{X}=\boldsymbol{x})$.
To render this approach feasible, we suggest the following. Fit a pre-estimate for the mean $\widehat \mu(\boldsymbol{x})$, using a strictly consistent score for mean estimation, and an FN network for $\boldsymbol{x} \mapsto \widehat \mu(\boldsymbol{x})$. This allows us to study the squared Pearson's residuals, resulting in a non-parametric regression problem, for $1\le i \le n$,
where the error terms $u_i$ should be centered ${\mathbb E}[u_i|\boldsymbol{X}_i]=0$. The goal is to solve (ref) for $c>0$ and $\phi''$. We suggest to simplify this problem and to turn it into a parametric problem by considering the following one-dimensional parametric family for $\phi$:
where $y>0$. This provides us with Bregman divergences, see (ref),
For $b\not\in(1,2)$, these are exactly the deviance losses within Tweedie's family Tweedie for power variance parameters $p=2-b$, see Example 4.11 in W\"uthrich--Merz WM2021. The case $b=2$ corresponds to the Gaussian distribution, $b=1$ to the Poisson distribution, $b=0$ to the gamma distribution and $b=-1$ to the inverse Gaussian distribution; for $b \in (1,2)$ there are no Tweedie's distributions, see Theorem 2 in J{\o}rgensen Jorgensen2. Thus, analyzing relation (ref) will motivate the specific choice of $b$ and of $\phi_b$, respectively, such that
and it allows us to select $c>0$.
For $\phi_+$ and $\phi_-$, we use a similar, though slightly more complicated, approach. First, we come up with a pre-estimate for the conditional quantile function $\widehat Q_\tau(\boldsymbol{x})$, along the lines of Subsection (ref). Using this pre-estimate $\widehat Q_\tau(\boldsymbol{x})$, we split our sample into two distinct sets, introducing the index sets $\mathcal I_- = \{i\in \{1,\ldots, n\}| Y_i\le \widehat Q_\tau(\boldsymbol{X}_i)\}$ and $\mathcal I_+ = \{i\in \{1,\ldots, n\}| Y_i > \widehat Q_\tau(\boldsymbol{X}_i)\}$. Fit conditional mean models $\widehat E_\tau^-(\boldsymbol{x})$ and $\widehat E_\tau^+(\boldsymbol{x})$ to the two data sets $\mathcal I_-$ and $\mathcal I_+$. Then, we can analyze
Again, one needs to estimate the parameters $c_-$ and $c_+$ as well as the functions $\phi''_-$ and $\phi''_+$, and, for simplicity, we again suggest to use a member of the parametric family (ref), so that we only need to select $b_-$, $c_-$, and $b_+$, $c_+$, respectively. Then, one obtains $\phi_-= \phi_{b_-}$ and $\phi_{+} = \phi_{b_+}$. However, one should invoke the restrictions that $\phi'_{+}<0$ and $\phi'_{-}>0$. Therefore, we have the restriction $b_->1$ for $\phi_{b_-}$ and $b_+<1$ for $\phi_{b_+}$.
We dispense with a discussion of the optimal choice of $g$. Equally so, we ignore conditions of the type of equation (4.19) in Dimitriadis et al. DimitriadisFisslerZiegel2020, stipulating that for $G_{e^-, e^+}$ in (ref) it holds that $G'_{e^-, e^+}\big(Q_\tau(\boldsymbol{x})\big) = c_\tau f_{Y|\boldsymbol{x}}\big(Q_\tau(\boldsymbol{x})\big)$ for some $c_\tau >0$ and for all $\boldsymbol{x}\in\mathcal X$. Here, $f_{Y|\boldsymbol{x}}$ is the conditional density of $Y$ given $\boldsymbol{X}=\boldsymbol{x}$. Unless we use strong parametric assumptions on $f_{Y|\boldsymbol{x}}$ to come up with a reasonable pre-estimate, we would need to resort to kernel density estimation methods here, which amounts to a considerable computational complexity. We think that in the given setting, the precise estimation of the splicing point corresponding to the estimation of the conditional quantile is less important than the precise estimation of the lower and upper conditional ES, yielding a good estimate of the overall expected claim size. Hence, we suggest to either set $g$ to be constant 0 or to use a multiple of the classical pinball loss, arising from $g(y) = c_\tau y$, $c_\tau>0$.
We conclude that the overall claim size (pure risk premium) can be calculated by
The beauty of this approach now is that we can fit different models above and below the splicing point $F_{Y|\boldsymbol{x}}^{-1}(\tau)$, accounting for different properties in the main body and the tail of the data. Thus, we can have different distributions above and below the splicing point, reflected in using different scoring functions above and below $F_{Y|\boldsymbol{x}}^{-1}(\tau)$ (implied by $\phi_-$ and $\phi_+$), resulting in different regression functions potentially using covariates $\boldsymbol{x} \in {\cal X}$ in different ways, i.e., we can have different regression functions in the tail and the main body of the data.
We present a real data example with claim amounts describing medical expenses in compulsory Swiss accident insurance. In total we have 267,992 claims with positive claim amounts, $Y_i>0$, and these claim amounts range from 1 to 691,066 CHF. Figure (ref) shows the empirical density and the log-log plot of these claim amounts. The empirical density is unimodal, and from the log-log plot we conclude that the tail is moderately heavy-tailed, but it is not regularly varying (which would correspond to an asymptotic straight line in the log-log plot).
These claim amounts are supported by 7 features. We have 3 categorical features for the `labor sector' of the injured, the `injury type' and the `injured body part', 1 binary feature telling whether the injury is a `work or leisure' accident, and 3 continuous features corresponding to the `age' of the injured, the `reporting delay' of the claim and the `accident quarter' (capturing seasonality in claims because leisure activities differ in summer and winter times). A preliminary analysis shows that all these features have predictive power, i.e., they are explaining systematic effects in the claim amounts.
For our analysis, we partition the entire data into learning data ${\cal L}=(Y_i,\boldsymbol{x}_i)_{1\le i \le n}$ that is used for model fitting, and test data ${\cal T}=(Y^\dagger_t, \boldsymbol{x}^\dagger_t)_{1\le t \le T}$ which we (only) use for an out-of-sample analysis. We do this partition stratified w.r.t. the claim amounts and in a ratio of $9:1$. This results in learning data ${\cal L}$ of size $n=241,193$ and in test data ${\cal T}$ of size $T=26,799$. We hold on to the same partition in all examples studied. For network fitting we further partition the learning data ${\cal L}=(Y_i,\boldsymbol{x}_i)_{1\le i \le n}$ into training data ${\cal U}$ and validation data ${\cal V}$. Thus, ${\cal U} \cup {\cal V}$ and ${\cal T}$ are disjoint (and assumed to be independent) so that we can perform a proper out-of-sample forecast evaluation. For a detailed discussion of such a partition of the data for SGD fitting we refer to Section 7.2.3 and, in particular, to Figure 7.7 in W\"uthrich--Merz WM2021.
We apply the framework of Section (ref) to perform a deep multiple quantile regression. We choose 3 probability levels $0<\tau_1< \tau_2 < \tau_3<1$ that we simultaneously estimate; the specific choices considered are $(\tau_1,\tau_2,\tau_3) = (10\%, 50\%, 90\%)$.
We first discuss pre-processing of feature components before choosing the deep FN network architecture $\boldsymbol{z}^{(d:1)}$. For the 3 categorical variables we use embedding layers of dimension 2, i.e., they are treated by a contextualized embedding. Embedding layers are explained in Section 7.4.1 of W\"uthrich--Merz WM2021. We use the {\sf R} library {\tt keras} for our implementation, and these embedding layers are encoded on lines 3-13 of Listing (ref) in the supplementary material. The binary variable is encoded by 0-1 and the continuous variables are pre-processed by the MinMaxScaler to ensure that they live on the same scale, see formula (7.30) in W\"uthrich--Merz WM2021 for the MinMaxScaler. Based on this feature encoding we use an FN network of depth $d=3$ having input dimension $r_0=3\cdot 2 + 1 + 3=10$ (for the categorical, binary and continuous feature components). This 10-dimensional variable enters the network on line 15 of Listing (ref). For the deep FN network architecture $\boldsymbol{z}^{(d:1)}$ we choose depth $d=3$ with $(r_1,r_2,r_3)=(20,15,10)$ hidden neurons in the hidden layers, and the hyperbolic tangent activation function $\Psi$. This is encoded on lines 16-18 of Listing (ref) and gives us learned representations $\boldsymbol{z}_i =\boldsymbol{z}^{(d:1)}(\boldsymbol{x}_i) \in {\mathbb R}^{r_d+1}$ of dimension $r_d+1=11$ of features $\boldsymbol{x}_i$. We remark that for insurance data of sample size of roughly 100,000 to 500,000 and with 10 to 20 feature components we have had good experiences by an FN network architecture of this complexity, confirmed by the various examples in W\"uthrich--Merz WM2021. Therefore, we hold on to this choice.
Next we implement an additive structure (ref) for deep multiple quantile regression
This requires that we (re-)use learned representations $\boldsymbol{z}_i =\boldsymbol{z}^{(d:1)}(\boldsymbol{x}_i)$ in the last hidden layer three times and the sequence of quantiles should be monotonically increasing in $\tau_j$. We choose as link functions $h$ and $h_+$ the log-link which provides us with the exponential function for their inverses, and we then model these quantiles recursively to preserve monotonicity. Lines 20-26 of Listing (ref) give the corresponding {\sf R} code, and line 28 outputs these ordered quantiles. The multiplicative approach is rather similar, and the corresponding changes in the {\sf R} code are shown in Listing (ref) in the supplementary material.
This network architecture has $r=834$ network parameters that need to be fitted to the learning data ${\cal L}$. We use pinball losses (ref) for the probability levels $(\tau_1,\tau_2,\tau_3) = (10\%, 50\%, 90\%)$; Listing (ref) in the supplementary material shows the corresponding {\sf R} code. These pinball losses then enter the compilation of the model on line 8 of Listing (ref). Moreover, we use the {\tt nadam} version of SGD which usually has a good fitting performance. We fit the two architectures (additive and multiplicative) to our learning data ${\cal L}$, using 80% of the learning data ${\cal L}$ as training data ${\cal U}$ and 20% as validation data ${\cal V}$ to explore the early stopping rule to prevent from over-fitting. To reduce the randomness of SGD fitting we repeat this procedure for 20 different starting points of the algorithm, and we calculate the average predictor over these 20 runs. The results are given in Table (ref).
Table (ref) shows the average out-of-sample pinball losses $L_{\tau_j}$ on the test data ${\cal T}$ of the two approaches, for $\tau_j \in \{10\%, 50\%, 90\%\}$. We note that the figures of the two approaches are very similar, and we cannot give a clear preference to one of the two approaches. We further explore these results.
Table (ref) shows the out-of-sample empirical coverage ratios on test data ${\cal T}$ of the event that the observation $Y_t^\dagger$ is smaller or equal to the estimated quantile, that is, we evaluate
where $Q_{\tau_j}(\boldsymbol{x}^\dagger_t)$ is the estimated quantile for level $\tau_j \in \{10\%, 50\%, 90\%\}$ (using either the additive or the multiplicative approach) and evaluated in the features $\boldsymbol{x}^\dagger_t$ of the out-of-sample observations $Y_t^\dagger$, $1\le t \le T$. Table (ref) verifies that the fitted deep multiple quantile regression finds the quantiles (on portfolio level) very well because $\widehat{\tau}_j \approx \tau_j$; we emphasize that these are out-of-sample figures.
Figure (ref) shows the estimated quantiles $Q_{\tau_j}(\boldsymbol{x}^\dagger_t)$ for individual features $\boldsymbol{x}^\dagger_t$ at probability levels $\tau_j \in \{10\%, 50\%, 90\%\}$ (blue, black, blue). The individual observations are ordered w.r.t. the estimated median in black. We observe that this ordering does not imply monotonicity for the other quantiles, especially for bigger claim potentials. This indicates heteroskedasticity in our data. The red dots show the corresponding (out-of-sample) observed claim amounts $Y^\dagger_t$.
We have analyzed these quantiles also in different granularity, for instance, we have analyzed the empirical coverage ratios on feature levels. Also there the figures look good, but we refrain from giving more plots because we would like to focus on the deep composite model regression proposed in Section (ref).
To fit a deep composite regression model we need a first preliminary step to explore relationship (ref). This motivates the choices of $\phi$, $\phi_+$ and $\phi_-$ in (ref)-(ref); for $g$ we use the identity function $g(y)=y$, giving us the pinball loss. For this preliminary step, we choose exactly the same network architecture (ref) as for deep quantile regression, the only change is that we replace the pinball loss by a strictly consistent scoring function for the mean functional, see Theorem (ref). As Bregman divergence we choose the gamma deviance loss, which corresponds to choice $b=0$ in (ref). We remark that the gamma model is the most popular model for insurance claim size modeling, and it often provides a first reasonable choice for a regression model; this is also the case for this preliminary step.
We fit three networks of depth $d=3$ with $(r_1,r_2,r_3)=(20,15,10)$ neurons and exponential output activation $h^{-1}$ to (a) all learning data ${\cal L}$, (b) the learning data having observations $Y_i > Q_{90\%}(\boldsymbol{x}_i)$, and (c) the learning data having observations $Y_i \le Q_{90\%}(\boldsymbol{x}_i)$, where $Q_{90\%}(\boldsymbol{x}_i)$ is the estimated $\tau=90\%$ quantile from the previous example of Section (ref). This allows us to analyze (ref) providing a choice for $\phi$, (ref) giving a choice for $\phi_+$, and (ref) giving a choice for $\phi_-$ for a deep composite model regression at probability level $\tau=90\%$.
Figure (ref) shows a spline fit and a linear regression to the estimated squared Pearson's residuals $(Y_i - \widehat{\mu}(\boldsymbol{x}_i))^2$, where $\widehat{\mu}(\cdot)$ is the estimated mean functional and where both axis are on the log-scale. The left-hand side shows all claims, the middle shows the situation where we only fit the network to the claims above the estimated quantile $Q_{90\%}(\boldsymbol{x}_i)$, and the right-hand side only considers the claims below that quantile.
Table (ref) gives the linear regression estimates for $c$ and $b$ in (ref) for the three cases. We observe that in all three cases we receive $b<1$ which results in derivatives $\phi'_b(y)=y^{b-1}/(b-1)<0$. This implies that we can only work under scoring function (ref), because $\phi_-$ requires $b>1$, e.g., the square loss function with $b=2$ would work for $\phi_-$ but this will not provide optimal convergence rates according to (ref).
We are now ready to fit the deep composite regression model for the claims $Y$ w.r.t. the probability level $\tau=90\%$. To this end, we implement the regression function (ref) and as deep FN network we use the same architecture as in the additive deep multiple quantile regression approach. This regression model is then fitted under the strictly consistent scoring function (ref) for the joint $\tau$-quantile and lower and upper ES (using the parameters of Table (ref)).
We fit this network under scoring function (ref), which is encoded in Listing (ref), using the {\tt nadam} version of SGD. Figure (ref) shows the out-of-sample estimated lower ES, $E^-_{90\%}(\boldsymbol{x}_t^\dagger)$, and upper ES, $E^+_{90\%}(\boldsymbol{x}_t^\dagger)$, against the estimated quantiles, $Q_{90\%}(\boldsymbol{x}_t^\dagger)$, of 2,000 randomly selected instances $\boldsymbol{x}_t^\dagger$, and the cyan lines present spline fits to all out-of-sample instances. Basically, the gaps between the cyan lines and the diagonal orange line describe the differences between the ES and the quantile. This gap is of constant size between the lower ES and the quantile above 7 (on the log-scale) which means that the lower ES is a fixed ratio of the quantile. The structure of the upper ES relative to the quantile is more complicated as the gap is becoming smaller with bigger values for the 90% quantile.
Based on these estimates we can now determine the expected value of $Y$, given $\boldsymbol{x}$,
We compare these estimated means to the ones obtained from the (plain vanilla) deep gamma model used in the preparatory Section (ref); note that the minimization of the gamma deviance loss with $b=0$ in (ref) is equivalent to the maximization of the log-likelihood function of the gamma distribution.
Figure (ref) (lhs) compares the deep gamma model to the fitted deep composite model. For small estimated means $\mu(\boldsymbol{x})$ the deep gamma and the deep composite model are rather similar, however, for large estimated means the deep gamma model provides clearly smaller estimates, see gray dotted lines in Figure (ref) (lhs). It seems that the deep gamma model under-estimates claims with large expected payments. To verify this, and to check the robustness of the results, we fit a second deep composite model to the data. To this end, we choose a different strictly consistent scoring function, namely, we choose (ref) which separates small and large claims in an additive way through $\phi_-$ and $\phi_+$. For the main body of the claims (below the 90%-quantile) we choose the square loss function for $\phi_-$ (which corresponds to choice $b=2$ in (ref)), and for the claims above the probability level we choose $b=0$ for $\phi_+$ in (ref), which corresponds to the gamma deviance loss. We call this second model `deep composite model 2'.
Figure (ref) (rhs) compares the two fitted deep composite models, using different score functions $L$, which implicitly implies different distributional assumptions in an MLE context. We observe that under both score functions we receive rather similar results, which verifies the robustness of the fittings. There are differences because SGD fitting with early stopping typically finds different `good' solutions, and this inherent fluctuation in SGD can only be reduced by averaging (blending) over model calibrations.
Analogously to checking the empirical out-of-sample coverage ratios in Table (ref), which validated the calibration of our deep quantile models, we now check calibration for the composite triplet consisting of the quantile, the lower ES and the upper ES. For the calibration of $Q_{90\%}$ we again use the empirical coverage defined in (ref), which should be close to $90\%$. Unfortunately, the lower and the upper ES not only fail to be elicitable, but also turn out not to be identifiable; see Supplement (ref) for a brief discussion of identifiability. Hence, we cannot check the calibration of the models $E_{90\%}^-$ and $E_{90\%}^+$ standalone, but we can only evaluate the goodness-of-fit of these models jointly with the corresponding quantile model $Q_{90\%}$. The corresponding joint strict identification functions are given as the first and third component in (ref). The empirical out-of-sample identification functions are
Values of $\widehat v_-$ and $\widehat v_+$ close to 0 indicate well calibrated models for $(E_{90\%}^-, Q_{90\%})$ and $(E_{90\%}^+, Q_{90\%})$, respectively.
Table (ref) reports the empirical coverage ratios $\widehat \tau$, for $\tau=90\%$, along with the empirical identification functions $\widehat v_-$ and $\widehat v_+$ for the two deep composite models and for the deep gamma model; having a gamma distributional assumption from the deep gamma approach of Section (ref) we can also calculate the composite triplet under this gamma model assumption. In Section (ref) we have fitted the means $\widehat{\mu}(\boldsymbol{x}^\dagger_t)$ under the gamma assumption. Moreover, using these means we can estimate a dispersion parameter. If we use the deviance dispersion estimate we receive in this gamma model a dispersion parameter of $1/\gamma=1.83$. This allows us to calculate the corresponding quantiles and ES in the gamma model. The lower and upper ES in the gamma model are obtained by
where $Y\sim \Gamma_{\gamma, \mu}$ denotes the gamma distribution with mean $\mu>0$ and shape parameter $\gamma>0$.
From Table (ref) we conclude that the deep composite regression models meet the right coverage ratio of 90% very well, whereas the deep gamma model results in a too high coverage ratio which indicates that the gamma distributional assumption does not match the tail of the observed data. This carries over to the identification functions of the lower and upper ES. The first deep composite model seems to be slightly better calibrated than the second one, which is indicated by values of $\widehat v_-$ and $\widehat v_+$ closer to zero for the first model. The deep gamma model exhibits a poor calibration reflected in high absolute values of $\widehat v_-$ and $\widehat v_+$.
Figure (ref) gives marginal estimations of the composite triplet of lower ES, the 90% quantile and the upper ES for the features `labor sector', `injury type' and `injured body part'; the $x$-axis is ordered w.r.t. the average 90% quantiles. We observe monotonicity between the average lower ES and the $90\%$ quantiles (black and blue dots), but this monotonicity gets lost w.r.t. the upper ES (red dots). This indicates that we have different regression functions for the main body and the tail of the data. That is, we have developed a very flexible (deep regression) model where features $\boldsymbol{x} \in {\cal X}$ impact predictions differently in the body and the tail of the distributions. This concludes the example.
We present deep quantile and deep composite model regression. For the composite model, we use a conditional quantile as splicing point, not an absolute threshold. This provides flexibility e.g. in the presence of heteroskedasticity. We utilize a network architecture which respects the natural ordering of quantiles at different probability levels as well as the natural ordering of the composite triplet consisting the the lower expected shortfall, the quantile and the upper expected shortfall at the same probability level. While strictly consistent scoring functions for tuples of quantiles are available, e.g. in the form of sums of pinball losses, and thus M-estimation can be performed for deep quantile models, we first derive and introduce the class of strictly consistent scoring functions for the composite triplet. Addressing the specific choice of the strictly consistent score for the composite triplet, we discuss data-driven choices which potentially increase the efficiency in estimation. The suitability of our methods is illustrated on a real data example with claim amounts describing medical expenses in compulsory Swiss accident insurance.
A relevant extension of our methods is a generalized composite model using two (or even more) splicing points. This allows one to model the influence of features differently for small claim sizes, large claim sizes and the body of the data. This extension seems particularly beneficial since there is empirical evidence that the (conditional) distribution of claim sizes is different in nature for small, claim sizes, the body, and for large claim sizes. Hence, this extension has the potential to increase accuracy of overall claim modeling, which is relevant in insurance pricing. The positive results on the elicitability of the range value at risk (RVaR) -- or interquantile expectation -- together with the two corresponding quantiles in Fissler--Ziegel FisslerZiegel2021 suggest the existence of strictly consistent scoring functions for the “extended composite quintuple” consisting of two quantiles (say, the 10% and 90% quantiles) together with the lower ES, upper ES and the range value at risk. However, the specific data-driven choice of the score motivated by efficiency considerations is far from being clear since the class of scoring functions for RVaR and two quantiles is less flexible than the one involving ES and a quantile. E.g. it is shown in FisslerZiegel2021 that there are no positively homogeneous strictly consistent scores for the former case. Therefore, we defer a detailed study of this extension to future research.
There is no competing interest.
{
}