Robust estimation and inference for expected shortfall regression with many regressors

Abstract

Expected shortfall (ES), also known as superquantile or conditional value-at-risk, is an important measure in risk analysis and stochastic optimisation and has applications beyond these fields. In finance, it refers to the conditional expected return of an asset given that the return is below some quantile of its distribution. In this paper, we consider a joint regression framework recently proposed to model the quantile and ES of a response variable simultaneously, given a set of covariates. The current state-of-the-art approach to this problem involves minimising a non-differentiable and non-convex joint loss function, which poses numerical challenges and limits its applicability to large-scale data. Motivated by the idea of using Neyman-orthogonal scores to reduce sensitivity to nuisance parameters, we propose a statistically robust and computationally efficient two-step procedure for fitting joint quantile and ES regression models that can handle highly skewed and heavy-tailed data. We establish explicit non-asymptotic bounds on estimation and Gaussian approximation errors that lay the foundation for statistical inference, even with increasing covariate dimensions. Finally, through numerical experiments and two data applications, we demonstrate that our approach well balances robustness, statistical, and numerical efficiencies for expected shortfall regression.

expected shortfall, heavy-tailed distribution, Huber loss, Neyman orthogonality, quantile regression

1 Introduction

Expected shortfall (ES), also known as superquantile or conditional value-at-risk (VaR), has been recognised as an important risk measure with versatile applications in finance (Acerbi & Tasche, 2002; Rockafellar & Uryasev, 2002), management science (Ben-Tal & Teboulle, 1986; Du & Escanciano, 2017), operations research (Rockafellar et al., 2014; Rockafellar & Uryasev, 2000), and clinical studies (He et al., 2010). For example, in finance, expected shortfall refers to the expected return of an asset or investment portfolio conditional on the return being below a lower quantile of its distribution, namely its VaR. In their Fundamental Review of the Trading Book (Basel Committee, 2016, 2019), the Basel Committee on Banking Supervision confirmed the replacement of VaR with ES as the standard risk measure in banking and insurance.

Let Y be a real-valued random variable with finite first-order absolute moment, $\mathbb{E}|Y|<\mathrm{\infty }$⁠, and let $F_Y$ be its cumulative distribution function (CDF). For any $\alpha \in (0,1)$⁠, the quantile and ES at level α are defined as

$$Q_{\alpha }(Y)=F_Y^{-1}(\alpha )=inf\{y\in \mathbb{R}:F_Y(y)\ge \alpha \}\text{and}{\mathrm{ES}}_{\alpha }(Y)=\mathbb{E}\{Y|Y\le Q_{\alpha }(Y)\},$$

(1)

respectively. If $F_Y$ is continuous, the α-level ES is equivalently given by (see, e.g. Lemma 2.16 of McNeil et al., 2015)

$${\mathrm{ES}}_{\alpha }(Y)=\frac{1}{\alpha }{\int }_0^{\alpha }{Q}_u(Y)\mathrm{d}u.$$

(2)

For instance, in socio-economics applications, Y is the income and ${\mathrm{ES}}_{\alpha }(Y)$ can be interpreted as the average income for the sub-population whose income falls below the α-quantile of the entire population. We refer the reader to Chapter 6 of Shapiro et al. (2014) and Rockafellar and Royset (2014) for a thorough discussion of ES and its mathematical properties.

With the increasing focus on ES and its desired properties as a risk measure, it is natural to examine the impact of a p-dimensional explanatory vector X, on the tail behaviour of Y through ES. One motivating example is the Job Training Partnership Act (JTPA), a large publicly funded training programme that provides training for adults with barriers to employment and out-of-school youths. The goal is to examine whether the training programme improves future income for adults with low-income earnings (Bloom et al., 1997), for which quantile regression (QR)-based approaches have been proposed (Abadie et al., 2002; Chernozhukov & Hansen, 2008). For example, the 0.05-quantile of the post-programme income refers to the highest income earning of those who have the 5% lowest income among the entire population, while the 0.05-ES concerns the average income earning within this sub-population and therefore is more scientifically relevant in the JTPA study.

Compared to the substantial body of literature on QR, extant works on ES estimation and inference in the presence of covariates are scarce. We refer the reader to Scaillet (2004), Cai and Wang (2008), Kato (2012), Linton and Xiao (2013), and Martins-Filho et al. (2018) for non-parametric conditional ES estimation, and more recently to Dimitriadis and Bayer (2019), Patton et al. (2019), and Barendse (2020) in the context of (semi-)parametric models. As suggested in Patton et al. (2019), this is partly because regulatory interest in ES as a risk measure is only recent, and also due to the fact that this measure is not elicitable (Gneiting, 2011). Let $\mathcal{P}$ be a class of distributions on ${\mathbb{R}}^d$⁠. We say that a statistical functional $\theta :\mathcal{P}\to D$ with $D\subseteq {\mathbb{R}}^p$ (⁠$p\ge 1$⁠) is elicitable relative to the class $\mathcal{P}$ if there exists a loss function $\rho :{\mathbb{R}}^d\times {\mathbb{R}}^p\to \mathbb{R}$ such that $\theta (F)={\mathrm{argmin}}_{\theta \in D}{\mathbb{E}}_{Z\sim F}\rho (Z,\theta )$ for any $F\in \mathcal{P}$⁠. Here, ${\mathbb{E}}_{Z\sim F}$ means that the expectation is taken with respect to the random variable Z that follows the distribution F. For example, the mean is elicitable using the $L_2$-loss, and the median is elicitable using the $L_1$-loss. Although the ES is not elicitable on its own, it is jointly elicitable with the quantile using a class of strictly consistent joint loss functions (Fissler & Ziegel, 2016). Based on this result, Dimitriadis and Bayer (2019) and Patton et al. (2019) proposed a joint regression model for the conditional α-level quantile and ES of Y, given the covariates $X\in {\mathbb{R}}^p$⁠. In this work, we focus on (conditional) linear joint quantile-ES models:

$$Q_{\alpha }(Y|X)=X^{\mathrm{T}}{\beta }^{*}\text{and}{\mathrm{ES}}_{\alpha }(Y|X)=X^{\mathrm{T}}{\theta }^{*}.$$

(3)

Equivalently, we have $\epsilon =Y-X^{\mathrm{T}}{\beta }^{*}$ and $\xi =Y-X^{\mathrm{T}}{\theta }^{*}$⁠, where the conditional α-quantile of $\epsilon$ and the conditional α-level expected shortfall of $\epsilon$⁠, given $X\in {\mathbb{R}}^p$⁠, are zero. More generally, one may allow the quantile and the ES models to depend on different covariate vectors $X_q$ and $X_e$⁠, respectively. In this case, the conditional α-quantile and α-ES of $\epsilon$ and $\xi$⁠, respectively, given $X=(X_q^{\mathrm{T}},X_e^{\mathrm{T}}{)}^{\mathrm{T}}$⁠, are assumed to be zero.

To jointly estimate ${\beta }^{*}$ and ${\theta }^{*}$⁠, Dimitriadis and Bayer (2019) and Patton et al. (2019) considered an M-estimator, defined as the global minimum of any member of a class of strictly consistent joint loss functions over some compact set (Fissler & Ziegel, 2016). The joint loss function, which will be specified in equation (5), is non-differentiable and non-convex. Dimitriadis and Bayer (2019) employed the derivative-free Nelder–Mead algorithm to minimise the resulting non-convex loss, which is a heuristic search method that may converge to non-stationary points. From a statistical perspective, they further established consistency and asymptotic normality for the global minima. However, from a computational perspective, finding the global minimum of a non-convex function is generally intractable: approximating the global minimum of a k-times continuously differentiable function $f:{\mathbb{R}}^p\to \mathbb{R}$ to $ϵ$-accuracy requires at least as many as $(1/ϵ{)}^{p/k}$ evaluations (ignoring problem-dependent constants) of the function and its first k-derivatives (Nemirovski & Yudin, 1983). The lack of differentiability makes this problem even more challenging numerically. To mitigate the computational effort, Barendse (2020) proposed a two-step procedure by first estimating the quantile parameters via standard QR, followed by least squares regression with generated response variables. Although computationally efficient, the ensuing estimator is sensitive to heavy-tailed error distributions due to the use of $L_2$-loss for fitting possibly highly skewed data in the second step; see Section 3.1 for a rigorous statement.

In this paper, we propose a new two-stage method for robust estimation and inference under a joint quantile and expected shortfall regression model (3), with a particular focus on the latter. Compared to existing approaches, our proposed method is robust against heavy-tailed errors without compromising statistical efficiency under light-tailed distributions. Computationally, our method can be implemented via fast and scalable gradient-based algorithms. The main contributions of this work are summarised as follows:

• Our method builds upon a recent approach to joint quantile and expected shortfall regression via a two-step procedure (Barendse, 2020). However, a general non-asymptotic theory for this approach has yet to be established. To fill this gap, we establish a finite-sample theoretical framework for the two-step ES estimator when the dimension of the model, p, increases with the number of observations, n. Specifically, we provide explicit upper bounds, as a function of $(n,p)$⁠, on the estimation error (under $L_2$-risk) and (uniform) Gaussian approximation errors; see Online Supplementary Section A. We also construct asymptotically valid (entrywise) confidence intervals for the ES parameters. The main computational effort of this two-step procedure is the QR fit in stage one. Therefore, we recommend using the convolution-smoothed QR method (Fernandes et al., 2021), which can be solved using fast first-order algorithms that are scalable to very large-scale problems (He et al., 2023). Our non-asymptotic theory allows the dimension p to grow with the sample size, which paves the way for analysing series/projection estimators under joint non-parametric quantile-ES models (Belloni et al., 2019) and penalised estimators under high-dimensional sparse models.

• The standard two-step estimator is a least squares estimator (LSE) with generated response variables. As a result, it is sensitive to the tails of the distribution of Y. We propose a robust ES regression method that applies adaptive Huber regression (Zhou et al., 2018) in the second step to address this issue. The resulting estimator achieves sub-Gaussian deviation bounds even when the (conditional) distribution of $Y|X$ only has Pareto-like tails. To achieve a trade-off between bias and robustness, we propose using a diverging robustification parameter $\tau =\tau (n,p)>0$⁠. In practice, we choose this hyper-parameter using a data-driven mechanism (L. Wang et al., 2021), guided by the non-asymptotic results presented in Section 4 and inspired by the censored equation approach introduced in Hahn et al. (1990). We have also developed efficient algorithms to compute standard and robust two-step ES estimators under additional constraints. These constraints ensure that the fitted ES does not exceed the fitted quantile at each observation. We refer the reader to the Online Supplementary Section D for more details.

• We conduct thorough numerical comparisons between the two-step estimator and the proposed robust variant with the joint M-estimator of Dimitriadis and Bayer (2019) on large synthetic data sets generated from a location-scale model, with both light- and heavy-tailed error distributions. To compute the joint M-estimator, we use the R package esreg, which is available at https://cran.r-project.org/package=esreg. To implement the proposed robust two-step procedure, we use a combination of R packages, quantreg or conquer and adaHuber. Our results show that the proposed robust ES regression approach achieves satisfying statistical performance, a higher degree of robustness against heavy-tailed error distributions, and superior computational efficiency and stability. We also demonstrate the effectiveness of our approach through numerical experiments and two real data examples.

In this work, the term ‘robustness’ specifically refers to the robustness against heavy-tailed distributions, as revealed by non-asymptotic deviation analysis dating back to Catoni (2012). In Catoni (2012)’s study of univariate mean estimation, it was found that while the sample mean has the optimal minimax mean squared error among all mean estimators, its deviation is worse for non-Gaussian samples than for Gaussian ones. Moreover, the worst-case deviation is sub-optimal when the sampling distribution has heavy tails. To be more specific, let $X_1,\dots ,X_n$ be independent copies of X with mean μ and variance ${\sigma }^2>0$⁠. Applying Chebyshev’s inequality to the empirical mean ${\overline{X}}_n=(1/n)\sum _{i=1}^n{X}_i$ yields a polynomial-type deviation bound: $|{\overline{X}}_n-\mu |\le \sigma \sqrt{1/(n\delta )}$ holds with probability at least $1-\delta$ for any $\delta \in (0,1)$⁠. Furthermore, if X is sub-Gaussian, meaning that $\mathbb{E}{e}^{\lambda (X-\mu )}\le e^{{\sigma }^2{\lambda }^2/2}$ for all $\lambda \in \mathbb{R}$⁠, then ${\overline{X}}_n$ can be referred to as a sub-Gaussian estimator, as it satisfies with probability $1-\delta$ that $|{\overline{X}}_n-\mu |\le \sigma \sqrt{2\log (2/\delta )/n}$⁠. In order to obtain sub-Gaussian deviations under a condition of bounded second moments, Fan et al. (2017) considered the Huber mean estimator ${\hat{\mu }}_{\tau }={\mathrm{argmin}}_{\theta \in \mathbb{R}}\sum _{i=1}^n{\ell }_{\tau }(X_i-\theta )$⁠, which is closely related to the method proposed by Catoni (2012). Here ${\ell }_{\tau }(\cdot )$ denotes the Huber loss; see definition (18). Theorem 5 in Fan et al. (2017) establishes that for any $v\ge \sigma$⁠, ${\hat{\mu }}_{\tau }$ with $\tau =v\sqrt{n/\log (1/\delta )}$ and $\delta \in (0,1)$ satisfies $|{\hat{\mu }}_{\tau }-\mu |\le 4v\sqrt{\log (1/\delta )/n}$ with probability at least $1-2\delta$ as long as $n\ge 8\log (1/\delta )$⁠. While Fan et al. (2017) does not explicitly state this, the divergence of τ in this context is also intended to strike a balance between bias and robustness. In comparison to (univariate) mean estimation, the problem of regression with growing dimensions and generated response variables present new technical challenges and requires more nuanced analysis. Nevertheless, the underlying phenomenon is quite similar.

1.1 Notation

For any two vectors $u=(u_1,\dots ,u_k{)}^{\mathrm{T}}$ and $v=(v_1,\dots ,v_k{)}^{\mathrm{T}}\in {\mathbb{R}}^k$⁠, we define their inner product as $u^{\mathrm{T}}v=⟨u,v⟩=\sum _{j=1}^k{u}_j{v}_j$⁠. We use $\Vert \cdot {\Vert }_p$$(1\le p\le \mathrm{\infty })$ to denote the ${\ell }_p$-norm in ${\mathbb{R}}^k$⁠: $\Vert u{\Vert }_p=(\sum _{i=1}^k|u_i{|}^p{)}^{1/p}$ for $p\ge 1$ and $\Vert u{\Vert }_{\mathrm{\infty }}=\underset{1\le i\le k}{max}|u_i|$⁠. Let ${\mathbb{S}}^{k-1}=\{u\in {\mathbb{R}}^k:\Vert u{\Vert }_2=1\}$ be the unit sphere in ${\mathbb{R}}^k$ under ${\ell }_2$-norm. Given a positive semi-definite matrix $A\in {\mathbb{R}}^{k\times k}$ and $u\in {\mathbb{R}}^k$⁠, let $\Vert u{\Vert }_A:=\Vert A^{1/2}u{\Vert }_2$⁠. Given an event/subset $\mathcal{A}$⁠, $\mathbb{1}(\mathcal{A})$ or ${\mathbb{1}}_{\mathcal{A}}$ denotes the zero-one indicator function for $\mathcal{A}$⁠. For two real numbers a and b, we write $a\wedge b=min\{a,b\}$ and $a\vee b=max\{a,b\}$⁠. For two sequences $\{a_n{\}}_{n\ge 1}$ and $\{b_n{\}}_{n\ge 1}$ of non-negative numbers, we write $a_n\lesssim b_n$ if $a_n\le C{b}_n$ for some constant $C>0$ independent of n, $a_n\gtrsim b_n$ if $b_n\lesssim a_n$⁠, and $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $a_n\gtrsim b_n$⁠.

2 Preliminaries and background

2.1 The joint regression framework

Assume we observe a sequence of data vectors $\{(Y_i,X_i){\}}_{i=1}^n$⁠, where $Y_i\in \mathbb{R}$ is the response variable, and $X_i\in {\mathbb{R}}^p$ is a p-dimensional vector of explanatory variables (covariates). For some fixed probability level $\alpha \in (0,1)$⁠, denote the conditional α-level quantile and ES of $Y_i$ given the covariates $X_i$ as $Q_{\alpha }(Y_i|X_i)$ and ${\mathrm{ES}}_{\alpha }(Y_i|X_i)$⁠, respectively. For the latter, we adhere to the definition ${\mathrm{ES}}_{\alpha }(Y_i|X_i)=\mathbb{E}\{Y_i|Y_i\le Q_{\alpha }(Y_i|X_i),X_i\}$⁠.

We consider the joint regression framework introduced in Dimitriadis and Bayer (2019) for modelling the conditional quantile and expected shortfall. For some probability level $\alpha \in (0,1)$⁠, assume that

$$Q_{\alpha }(Y_i|X_i)=X_i^{\mathrm{T}}{\beta }^{*},{\mathrm{ES}}_{\alpha }(Y_i|X_i)=X_i^{\mathrm{T}}{\theta }^{*},$$

(4)

where ${\beta }^{*},{\theta }^{*}\in {\mathbb{R}}^p$ are the unknown true underlying parameters for quantile and ES, respectively. Fissler and Ziegel (2016) explained that quantile and ES are jointly elicitable and proposed a class of strictly consistent joint loss functions for quantile and ES estimation. Let $G_1$ be an increasing and integrable function, and let ${\mathcal{G}}_2$ be a three times continuously differentiable function such that both ${\mathcal{G}}_2$ and its derivative $G_2={\mathcal{G}}_2^{\prime }$ are strictly positive. The proposed joint loss function in Fissler and Ziegel (2016) takes the form

$$\begin{array}{rl}s(\beta ,\theta ;Y,X)& =\{\alpha -\mathbb{1}(Y\le X^{\mathrm{T}}\beta )\}\{G_1(Y)-G_1(X^{\mathrm{T}}\beta )\}\\ & +\frac{G_2(X^{\mathrm{T}}\theta )}{\alpha }\left\{{\displaystyle \underset{=:s_0(\beta ,\theta ;Y,X)}{{\displaystyle \underbrace{\alpha X^{\mathrm{T}}(\theta -\beta )-(Y-X^{\mathrm{T}}\beta )\mathbb{1}(Y\le X^{\mathrm{T}}\beta )}}}}\right\}-{\mathcal{G}}_2(X^{\mathrm{T}}\theta ).\end{array}$$

(5)

This general form also includes the joint loss function proposed by Acerbi and Székely (2014) by taking $G_1(x)=-(W/2)x^2$ for some $W\in \mathbb{R}$ and ${\mathcal{G}}_2(x)=\alpha x^2/2$⁠.

In the regression framework with a fixed number of covariates, Dimitriadis and Bayer (2019) established the consistency and asymptotic normality of the M-estimator $({\tilde{\beta }}^{\mathrm{T}},{\tilde{\theta }}^{\mathrm{T}}{)}^{\mathrm{T}}$⁠, defined as

$$\left(\begin{array}{c}\tilde{\beta }\\ \tilde{\theta }\end{array}\right)\in {\displaystyle \underset{\beta ,\theta \in \mathrm{\Theta }}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}s(\beta ,\theta ;Y_i,X_i),$$

(6)

where $\mathrm{\Theta }\subseteq {\mathbb{R}}^p$ is the parameter space, assumed to be compact, convex, and has non-empty interior. The main challenge of the aforementioned approach is that the objective function in equation (6) is non-differentiable and non-convex for any feasible choice of the functions $G_1$ and $G_2$ (Fissler & Ziegel, 2016). Note from definition (1) that the expected shortfall depends on the quantile, not vice versa. The estimation and inference of ${\theta }^{*}$ is thus the main challenge. It is, however, infeasible to estimate a single regression model for ES through M-estimation, that is, by minimising some strictly consistent loss function (Dimitriadis & Bayer, 2019).

In the joint regression framework, if the main goal is to estimate and forecast ES, then ${\beta }^{*}$ can be naturally viewed as a nuisance parameter. Motivated by the idea of using Neyman-orthogonal scores to reduce sensitivity with respect to nuisance parameters (Barendse, 2020; Chernozhukov et al., 2018; Neyman, 1979) proposed a two-stage procedure that bypasses non-convex optimisation problems. In the first stage, an estimate $\hat{\beta }$ of ${\beta }^{*}$ is obtained via standard QR. The second step employs an orthogonal score with fitted thresholding quantiles to estimate ${\theta }^{*}$⁠. The key observation is as follows. Define the function

$$\begin{array}{rl}{\psi }_0(\beta ,\theta ;X)& =\mathbb{E}\{s_0(\beta ,\theta ;Y,X)|X\}\\ & =\alpha X^{\mathrm{T}}\theta -\mathbb{P}(Y\le X^{\mathrm{T}}\beta |X)\mathbb{E}(Y|Y\le X^{\mathrm{T}}\beta ,X)+\{\mathbb{P}(Y\le X^{\mathrm{T}}\beta |X)-\alpha \}{X}^{\mathrm{T}}\beta ,\end{array}$$

(7)

where $s_0$ is given in equation (5). Under model (4), we have ${\psi }_0({\beta }^{*},{\theta }^{*};X)=0$ almost surely over X. Let $F_{Y|X}$ be the conditional distribution function of Y given X. Provided that $F_{Y|X}$ is continuously differentiable, taking the gradient with respect to β on both sides of the above equality yields

$${\mathrm{\partial }}_{\beta }{\psi }_0(\beta ,\theta ;X)=\{F_{Y|X}(X^{\mathrm{T}}\beta )-\alpha \}X,\text{ for any}\beta ,\theta \in {\mathbb{R}}^p.$$

We hence refer to the following property:

$${{\mathrm{\partial }}_{\beta }{\psi }_0(\beta ,\theta ;X)|}_{\beta ={\beta }^{*}}=\{F_{Y|X}(X^{\mathrm{T}}{\beta }^{*})-\alpha \}X=0$$

(8)

as Neyman orthogonality.

2.2 Two-step ES estimation via Neyman-orthogonal score

We start with a detailed overview of the two-step approach proposed by Barendse (2020) using the Neyman-orthogonal score (7) under the joint model (4). In Section 3.1, we will develop a non-asymptotic (finite-sample) theory for the two-step ES estimator, $\hat{\theta }$⁠, under the regime in which p is allowed to increase with the sample size n. We further develop asymptotic normality results for individual coordinates, or more generally linear projections, of $\hat{\theta }$⁠, in the increasing-dimension regime ‘$p^2/n=o(1)$’. Our non-asymptotic results and techniques pave the way for analysing high-dimensional sparse quantile-ES models.

The first step involves computing the standard QR estimator of ${\beta }^{*}$⁠:

$$\hat{\beta }\in {\displaystyle \underset{\beta \in {\mathbb{R}}^p}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\rho }_{\alpha }(Y_i-X_i^{\mathrm{T}}\beta ),$$

(9)

where ${\rho }_{\alpha }(u)=\{\alpha -\mathbb{1}(u<0)\}u$ is the check function (Koenker & Bassett, 1978). The second step is motivated by the orthogonal score $s_0$ in equation (5). Specifically, let $\hat{\mathcal{L}}(\beta ,\theta )=(1/n)\sum _{i=1}^n{s}_i^2(\beta ,\theta )$ be the joint empirical loss with

$$s_i(\beta ,\theta ):=s_0(\beta ,\theta ;Y_i,X_i)=\alpha X_i^{\mathrm{T}}\theta -\mathbb{1}(Y_i\le X_i^{\mathrm{T}}\beta )Y_i+\{\mathbb{1}(Y_i\le X_i^{\mathrm{T}}\beta )-\alpha \}{X}_i^{\mathrm{T}}\beta .$$

(10)

Given $\hat{\beta }$ obtained from the first step, the ES estimator $\hat{\theta }$ of ${\theta }^{*}$ is computed as

$$\hat{\theta }\in {\displaystyle \underset{\theta \in {\mathbb{R}}^p}{argmin}}\hat{\mathcal{L}}(\hat{\beta },\theta ).$$

(11)

For any β fixed, the function $\theta \mapsto \hat{\mathcal{L}}(\beta ,\theta )$ is convex with gradient and Hessian given by

$${\mathrm{\partial }}_{\theta }\hat{\mathcal{L}}(\beta ,\theta )=\frac{2\alpha }{n}{\displaystyle \sum _{i=1}^n}{s}_i(\beta ,\theta )X_i\text{and}{\mathrm{\partial }}_{\theta }^2\hat{\mathcal{L}}(\beta ,\theta )=\frac{2{\alpha }^2}{n}{\displaystyle \sum _{i=1}^n}{X}_i{X}_i^{\mathrm{T}},$$

respectively. By the first-order condition, the ES regression estimator $\hat{\theta }$ satisfies the moment condition ${\mathrm{\partial }}_{\theta }\hat{\mathcal{L}}(\hat{\beta },\hat{\theta })=0$⁠, and has a closed-form expression

$$\hat{\theta }=\hat{\beta }+{\left({\displaystyle \sum _{i=1}^n}{X}_i{X}_i^{\mathrm{T}}\right)}^{-1}\frac{1}{\alpha }{\displaystyle \sum _{i=1}^n}(Y_i-X_i^{\mathrm{T}}\hat{\beta })X_i\mathbb{1}(Y_i\le X_i^{\mathrm{T}}\hat{\beta }),$$

(12)

provided that $\mathbb{X}=(X_1,\dots ,X_n{)}^{\mathrm{T}}\in {\mathbb{R}}^{n\times p}$ is full-rank.

 

Remark 1

When p is large, we suggest using the convolution-smoothed QR (conquer) estimator (Fernandes et al., 2021; He et al., 2023) in the first step, which can be computed by fast and scalable gradient-based algorithms. Given a smoothing parameter/bandwidth $h>0$⁠, the conquer estimator ${\hat{\beta }}_h$ minimises the convolution-smoothed loss function $\beta \mapsto \sum _{i=1}^n{\rho }_{\alpha ,h}(Y_i-X_i^{\mathrm{T}}\beta )$ with ${\rho }_{\alpha ,h}(u)=({\rho }_{\alpha }*K_h)(u)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\rho }_{\alpha }(u)K_h(v-u)\mathrm{d}v$⁠, where $K_h(u):=(1/h)K(u/h)$ for some symmetric, non-negative kernel function K, and * is the convolution operator. We refer to Fernandes et al. (2021) and He et al. (2023) for more details, including both asymptotic and finite-sample properties of ${\hat{\beta }}_h$ when p is fixed and growing as well as the bandwidth selection.

Define $p\times p$ matrices $\mathrm{\Sigma }=\mathbb{E}(X{X}^{\mathrm{T}})$ and $\mathrm{\Omega }=\mathbb{E}({\omega }^{2}X{X}^{\mathrm{T}})$ with $\omega :=(Y-X^{\mathrm{T}}{\beta }^{*})\mathbb{1}(Y\le X^{\mathrm{T}}{\beta }^{*})+\alpha X^{\mathrm{T}}({\beta }^{*}-{\theta }^{*})$ satisfying $\mathbb{E}(\omega |X)=0$ under model (4). Provided that $p=p_n$ satisfies $p^2/n\to 0$⁠, we will show in the Online Supplementary Theorem A.3 that ${\hat{\theta }}_j$ is asymptotically normal:

$$\frac{\alpha \sqrt{n}({\hat{\theta }}_j-{\theta }_j^{*})}{\sqrt{({\mathrm{\Sigma }}^{-1}\mathrm{\Omega }{\mathrm{\Sigma }}^{-1}{)}_{jj}}}\mathrm{d}\to \mathcal{N}(0,1)\text{as}n,p\to \mathrm{\infty }.$$

As a direct implication, an asymptotically valid entrywise confidence interval for ${\theta }^{*}$ can be constructed as follows. Recall that $(\hat{\beta },\hat{\theta })$ is the joint quantile-ES regression estimators given in equations (9) and (11), respectively. Define the estimated ‘residuals’ as

$${\hat{\epsilon }}_i=Y_i-X_i^{\mathrm{T}}\hat{\beta }\text{and}{\hat{\omega }}_i={\hat{\epsilon }}_i\wedge 0+\alpha X_i^{\mathrm{T}}(\hat{\beta }-\hat{\theta }).$$

(13)

We then use the sample analogue of Σ and a plug-in estimator of Ω, namely, $\hat{\mathrm{\Sigma }}=\sum _{i=1}^n{X}_i{X}_i^{\mathrm{T}}/n$ and $\hat{\mathrm{\Omega }}=\sum _{i=1}^n{\hat{\omega }}_i^2{X}_i{X}_i^{\mathrm{T}}/n$⁠. Consequently, we construct (approximate) 95% confidence interval for each coefficient as

$$[{\hat{\theta }}_j-\frac{1.96}{\alpha \sqrt{n}}({\hat{\mathrm{\Sigma }}}^{-1}\hat{\mathrm{\Omega }}{\hat{\mathrm{\Sigma }}}^{-1}{)}_{jj}^{1/2},{\hat{\theta }}_j+\frac{1.96}{\alpha \sqrt{n}}({\hat{\mathrm{\Sigma }}}^{-1}\hat{\mathrm{\Omega }}{\hat{\mathrm{\Sigma }}}^{-1}{)}_{jj}^{1/2}],j=1,\dots ,p.$$

(14)

3 Robust expected shortfall regression

3.1 Motivation

The two-step estimator $\hat{\theta }$ given in equation (12) is essentially an LSE with generated response variables. While the two-step procedure is computationally efficient and enjoys nice asymptotic properties, due to the use of the least squares type loss, it is sensitive to outliers or heavy-tailed data that is ubiquitous in various areas such as climate, insurance claims, and genomics data. In particular, heavy-tailedness has become a well-known stylised fact of financial returns and stock-level predictor variables (Cont, 2001). Since the expected shortfall is a quantity that describes the tail behaviour of a distribution, it is important to construct an estimator that is robust to the power-law or Pareto-like tails.

To motivate the need for a robust ES estimator, we start with the non-regression setting in which $X_i\equiv 1$⁠. The two-step ES estimator (12) can then be simplified as

$${\widehat{\mathrm{ES}}}_{\alpha }=\frac{1}{\alpha n}{\displaystyle \sum _{i=1}^n}{Y}_i\mathbb{1}\{Y_i\le {\hat{Q}}_{\alpha }\}+{\hat{Q}}_{\alpha }\{1-\hat{F}({\hat{Q}}_{\alpha })/\alpha \},$$

(15)

where $\hat{F}$ is the empirical CDF of Y and ${\hat{Q}}_{\alpha }={\hat{F}}^{-1}(\alpha )$ is the sample quantile. The estimator ${\widehat{\mathrm{ES}}}_{\alpha }$ (15) coincides with the ES estimate (4) in Bassett et al. (2004), although the latter is motivated differently by the following property:

$${\mathrm{ES}}_{\alpha }(Y)=\mathbb{E}(Y)-\frac{1}{\alpha }{\displaystyle \underset{\beta \in \mathbb{R}}{min}}\mathbb{E}{\rho }_{\alpha }(Y-\beta ).$$

Since $|\hat{F}({\hat{Q}}_{\alpha })-\alpha |\le 1/n$⁠, up to higher order terms, ${\widehat{\mathrm{ES}}}_{\alpha }$ equals $(\alpha n{)}^{-1}\sum _{i=1}^n{Y}_i\mathbb{1}\{Y_i\le {\hat{Q}}_{\alpha }\}$ which, by the consistency of sample quantiles, is first-order equivalent to the ‘oracle’ ES estimator $\widehat{{\mathrm{ES}}_{\alpha }^{\mathrm{ora}}}:=(\alpha n{)}^{-1}\sum _{i=1}^n{Y}_i\mathbb{1}\{Y_i\le Q_{\alpha }(Y)\}$⁠.

Since the truncated variable $Y_i\mathbb{1}\{Y_i\le Q_{\alpha }(Y)\}$ can be highly left-skewed with heavy tails, the corresponding empirical mean is sensitive to the (left) tails of the distribution of Y, and hence lacks robustness against heavy-tailed data. Specifically, let $X_1,\dots ,X_n$ be i.i.d. random variables with mean μ and variance ${\sigma }^2>0$⁠. When $X_i$ is sub-Gaussian (i.e. $\mathbb{E}(e^{\lambda X_i})\le e^{{\lambda }^2{\sigma }^2/2}$ for any $\lambda \in \mathbb{R}$⁠), it follows from the Chernoff bound (Chernoff, 1952) that

$$\mathbb{P}\{|{\overline{X}}_n-\mu |\ge \sigma \sqrt{2\log (2/\delta )/n}\}\le \delta ,\mathrm{valid}\mathrm{for}\mathrm{any}\delta \in (0,1).$$

(16)

In other words, the sample mean ${\overline{X}}_n=(1/n)\sum _{i=1}^n{X}_i$ satisfies the sub-Gaussian deviation bound. On the other hand, the following proposition provides a lower bound for the deviations of the empirical mean $(1/n)\sum _{i=1}^n{Y}_i\mathbb{1}\{Y_i\le Q_{\alpha }(Y)\}$ when the distribution of Y is the least favourable among all heavy-tailed distributions with mean zero and variance ${\sigma }^2$⁠.

 

Proposition 1

For any value of the standard deviation $\sigma >0$ and any probability level $\delta \in (0,e^{-1}]$⁠, there exists some distribution with mean zero and variance ${\sigma }^2$ such that for any $\alpha \in (0,1)$⁠, the i.i.d. sample $\{Y_i{\}}_{i=1}^n$ of size n drawn from it satisfies

$$\mathbb{P}[\frac{1}{n}{\displaystyle \sum _{i=1}^n}{Y}_i\mathbb{1}(Y_i\le Q_{\alpha })-\mathbb{E}\{Y\mathbb{1}(Y\le Q_{\alpha })\}\le -\sigma \sqrt{\frac{1}{\delta n}}\cdot \frac{1-e\delta }{\sqrt{2e}}]\ge \delta ,$$

(17)

as long as $n\ge e\delta /\alpha$⁠, where $Q_{\alpha }=Q_{\alpha }(Y)$ is the αth quantile of Y.

Together, the upper and lower bounds (16) and (17) show that the worst case deviations of the empirical mean are sub-optimal when the underlying distribution is heavy-tailed (as opposed to having Gaussian-like thin tails). If Y follows a heavy-tailed distribution, such as the t- or Pareto distributions, then the left-truncated variables $Z_i:=Y_i\mathbb{1}\{Y_i\le Q_{\alpha }(Y)\}$ have not only heavy but also asymmetric tails. In this case, the empirical mean $(\alpha n{)}^{-1}\sum _{i=1}^n{Z}_i$ can be a sub-optimal estimator of ${\mathrm{ES}}_{\alpha }(Y)$⁠.

3.2 Robust estimation and inference via adaptive Huber regression

To robustify the ES regression estimator (12) in the presence of skewed heavy-tailed observations, we utilise the idea of adaptive Huber regression in Zhou et al. (2018). For some $\tau >0$⁠, the Huber loss (Huber, 1973) takes the form

$$\begin{array}{r}{\ell }_{\tau }(u)=\{\begin{array}{cc}{u}^2/2& \text{if}|u|\le \tau ,\\ \tau |u|-{\tau }^2/2& \text{if}|u|>\tau .\end{array}\end{array}$$

(18)

We propose a robust/Huberised ES regression estimator defined as

$${\hat{\theta }}_{\tau }\in {\displaystyle \underset{\theta \in {\mathbb{R}}^p}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\ell }_{\tau }(s_i(\hat{\beta },\theta )),$$

(19)

where $s_i(\hat{\beta },\theta )$ is as defined in equation (10), and $\tau >0$ is a robustification parameter that should be calibrated adaptively from data.

To see this, we consider the oracle Huber ES estimator defined as

$${\hat{\theta }}_{\tau }^{\mathrm{ora}}\in {\displaystyle \underset{\theta \in {\mathbb{R}}^p}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\ell }_{\tau }(s_i({\beta }^{*},\theta ))={\displaystyle \underset{\theta \in {\mathbb{R}}^p}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\ell }_{\tau }(Z_i-\alpha X_i^{\mathrm{T}}\theta ),$$

(20)

where $Z_i=(Y_i-X_i^{\mathrm{T}}{\beta }^{*})\mathbb{1}(Y_i\le X_i^{\mathrm{T}}{\beta }^{*})+\alpha X_i^{\mathrm{T}}{\beta }^{*}$⁠. For any $\tau >0$⁠, ${\hat{\theta }}_{\tau }^{\mathrm{ora}}$ is an M-estimator of its population counterpart

$${\theta }_{\tau }^{*}={\displaystyle \underset{\theta \in {\mathbb{R}}^p}{argmin}}\mathbb{E}\{{\ell }_{\tau }(Z_i-\alpha X_i^{\mathrm{T}}\theta )\}.$$

Let ${\psi }_{\tau }(t)={\ell }_{\tau }^{\prime }(t)=\mathrm{sign}(t)min(|t|,\tau )$ be the derivative of the Huber loss. By the convexity of the Huber loss, ${\theta }_{\tau }^{*}$ must satisfy the first-order condition $\mathbb{E}\{{\psi }_{\tau }(Z_i-\alpha X_i^{\mathrm{T}}{\theta }_{\tau }^{*})X_i\}=0$⁠. On the other hand, define the ES deviations ${\omega }_i=Z_i-\alpha X_i^{\mathrm{T}}{\theta }^{*}$⁠, satisfying $\mathbb{E}({\omega }_i|X_i)=0$ and $\mathbb{E}({\omega }_i)=0$⁠. Since the conditional distribution of ${\omega }_i$ given $X_i$ is asymmetric, in general we have $\mathbb{E}\{{\psi }_{\tau }(Z_i-\alpha X_i^{\mathrm{T}}{\theta }^{*})X_i\}=\mathbb{E}\{{\psi }_{\tau }({\omega }_i)X_i\}\ne 0$⁠, which in turn implies that ${\theta }_{\tau }^{*}\ne {\theta }^{*}$⁠. We thus refer to their difference under the ${\ell }_2$-norm, $\Vert {\theta }_{\tau }^{*}-{\theta }^{*}{\Vert }_2$⁠, as the robustification bias. Proposition 2 provides an upper bound for the robustification bias, which depends on τ and some moment parameters. In particular, τ needs to diverge for the robustification bias to diminish.

 

Proposition 2

Assume that $\epsilon :=Y-X^{\mathrm{T}}{\beta }^{*}$ satisfies ${\mathrm{var}}_X(\epsilon \wedge 0)\le {\overline{\sigma }}^2$ almost surely for some constant $\overline{\sigma }>0$⁠, and that ${\kappa }_4=\underset{u\in {\mathbb{S}}^{p-1}}{sup}\mathbb{E}⟨u,{\mathrm{\Sigma }}^{-1/2}X{⟩}^4<\mathrm{\infty }$⁠, where $\mathrm{\Sigma }=\mathbb{E}(X{X}^{\mathrm{T}})$ is positive definite. Then, for any $\tau \ge 2{\kappa }_4^{1/4}\overline{\sigma }$⁠, we have $\Vert {\theta }_{\tau }^{*}-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\le 2{\overline{\sigma }}^2/(\alpha \tau )$⁠.

In Section 4, we investigate the finite-sample properties of the robust ES estimator ${\hat{\theta }}_{\tau }$ obtained via equations (9) and (19): our results include a deviation inequality for $\Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}$ (Theorem 1), the Bahadur representation (Theorem 2), and a Berry–Esseen bound for linear projections of ${\hat{\theta }}_{\tau }$ and ${\hat{\theta }}_{\tau }^{\mathrm{ora}}$ (Theorem 3). With a properly chosen τ that is of order $\tau \asymp \overline{\sigma }\sqrt{n/p}$⁠, we will show that $\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\lesssim \overline{\sigma }\sqrt{p/n}$ with high probability. Moreover, for any deterministic vector $a\in {\mathbb{R}}^p$⁠, the standardised statistic $\alpha \sqrt{n}⟨a,{\hat{\theta }}_{\tau }-{\theta }^{*}⟩/{\varrho }_a$ converges in distribution to $\mathcal{N}(0,1)$⁠, where ${\varrho }_a^2=a^{\mathrm{T}}{\mathrm{\Sigma }}^{-1}\mathrm{\Omega }{\mathrm{\Sigma }}^{-1}a$ and $\omega =(Y-X^{\mathrm{T}}{\beta }^{*})\mathbb{1}(Y\le X^{\mathrm{T}}{\beta }^{*})+\alpha X^{\mathrm{T}}({\beta }^{*}-{\theta }^{*})$⁠. Our theoretical analysis reveals two attractive properties of the adaptive Huberised ES estimator ${\hat{\theta }}_{\tau }$⁠: (i) the non-asymptotic deviation upper bounds for ${\hat{\theta }}_{\tau }$ are much smaller in order than those for $\hat{\theta }$ at any given confidence level and (ii) the asymptotic relative efficiency of ${\hat{\theta }}_{\tau }$ to $\hat{\theta }$ is one. Moreover, Theorem 3 shows that the two-step robust estimator (with estimated conditional quantiles) is asymptotically equivalent to the oracle Huberised estimator (20) (assuming ${\beta }^{*}$ were known). This further justifies the usefulness of the Neyman-orthogonal score, which makes the QR estimation error first-order negligible.

Consistent estimators of Σ and $\mathrm{\Omega }=\mathbb{E}({\omega }^{2}X{X}^{\mathrm{T}})$ are useful for statistical inference. Given the pair of quantile-ES regression estimators $(\hat{\beta },{\hat{\theta }}_{\tau })$⁠, with a slight abuse of notation we use ${\hat{\epsilon }}_i$ and ${\hat{\omega }}_i$ to denote the fitted QR and ES residuals as in equation (13) except with $\hat{\theta }$ replaced by ${\hat{\theta }}_{\tau }$⁠. As discussed in Section 2.2, a naive estimate of Ω is $\hat{\mathrm{\Omega }}=(1/n)\sum _{i=1}^n{\hat{w}}_i^2{X}_i{X}_i^{\mathrm{T}}$⁠. In the presence of heavy-tailed errors ${\epsilon }_i$⁠, even the ‘oracle’ estimate $\tilde{\mathrm{\Omega }}=(1/n)\sum _{i=1}^n{w}_i^2{X}_i{X}_i^{\mathrm{T}}$ performs poorly and tends to overestimate. Motivated by Huber regression, we further propose a simple truncated estimator of Ω given by

$${\hat{\mathrm{\Omega }}}_{\gamma }=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\psi }_{\gamma }^2({\hat{\omega }}_i)X_i{X}_i^{\mathrm{T}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}min\{|{\hat{\omega }}_i|,\gamma {\}}^2{X}_i{X}_i^{\mathrm{T}},$$

(21)

where $\gamma =\gamma (n,p)>0$ is a second robustification parameter. Consequently, we construct approximate 95% robust confidence intervals for ${\theta }_j^{*}$’s as

$$[{\hat{\theta }}_{\tau ,j}-\frac{1.96}{\alpha \sqrt{n}}({\hat{\mathrm{\Sigma }}}^{-1}{\hat{\mathrm{\Omega }}}_{\gamma }{\hat{\mathrm{\Sigma }}}^{-1}{)}_{jj}^{1/2},{\hat{\theta }}_{\tau ,j}+\frac{1.96}{\alpha \sqrt{n}}({\hat{\mathrm{\Sigma }}}^{-1}{\hat{\mathrm{\Omega }}}_{\gamma }{\hat{\mathrm{\Sigma }}}^{-1}{)}_{jj}^{1/2}],j=1,\dots ,p.$$

(22)

The convergence rate of ${\hat{\mathrm{\Omega }}}_{\gamma }$ with a suitably chosen γ will be discussed in Section 4.

As previously discussed, the robustification parameter τ plays a crucial role in achieving a balance between bias and robustness against heavy-tailed error distributions. This balance is necessary because of the asymmetric nature of the ES residual $\omega =\epsilon \wedge 0+\alpha X^{\mathrm{T}}({\beta }^{*}-{\theta }^{*})$ with $\epsilon =Y-X^{\mathrm{T}}{\beta }^{*}$⁠. Assuming that the (conditional) variance of ${\epsilon }_{-}=\epsilon \wedge 0$ is bounded, i.e. ${\mathrm{var}}_X({\epsilon }_{-})\le {\overline{\sigma }}^2$ (almost surely) for some $\overline{\sigma }>0$⁠, Theorem 1 suggests that to achieve a tight deviation bound at the $1-\delta$ confidence level for any given $\delta \in (0,1)$⁠, the robustification parameter $\tau =\tau (n,p)$ should be of order $\overline{\sigma }\sqrt{n/(p+\log {\delta }^{-1})}$⁠. In practice, the scale of $\overline{\sigma }$ is typically unknown. A useful heuristic is to substitute it with the sample standard deviation of the negative QR residuals $\{{\hat{\epsilon }}_{i,-}=min(Y_i-X_i^{\mathrm{T}}\hat{\beta },0){\}}_{i=1}^n$⁠, which we denote by $\hat{\sigma }$⁠. Here, $\hat{\beta }$ refers to the first-stage QR estimator. Using $\hat{\tau }=\hat{\sigma }\sqrt{n/(p+\log {\delta }^{-1})}$ as a data-driven proxy for τ, the resulting estimator is also location and scale equivariant.

In the following, we present a refined data-driven approach for selecting τ that consistently outperforms the previously mentioned rule of thumb in the numerical experiments conducted in Section 6. This approach is adapted from the method proposed in L. Wang et al. (2021) and draws inspiration from the censored equation approach originally introduced by Hahn et al. (1990) as a proof technique for deriving robust weak convergence theory for self-normalised sums. Note that for each $\tau >0$⁠, the Huber ES estimator ${\hat{\theta }}_{\tau }$ can be defined equivalently as the solution to the estimating equation $\sum _{i=1}^n{\psi }_{\tau }({\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta )X_i=0$⁠, $\theta \in {\mathbb{R}}^d$⁠, where ${\hat{Z}}_i={\hat{\epsilon }}_{i,-}+\alpha X_i^{\mathrm{T}}\hat{\beta }$ are the generated response variables, and $\hat{\beta }$ denotes the initial QR estimator. Since the optimal choice of τ is proportional to the noise scale, we propose to estimate ${\theta }^{*}$ and the unknown noise scale simultaneously by solving the following system of equations for $(\theta ,s)\in {\mathbb{R}}^p\times (0,\mathrm{\infty })$⁠:

$$\begin{array}{rl}{g}_1(\theta ,s)& :=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\psi }_k\left(\frac{{\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta }{s}\right)X_i=0,\\ g_2(\theta ,s)& :=\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\psi }_k^2\left(\frac{{\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta }{s}\right)-1=0,\end{array}$$

where $k=k(n,p)=\sqrt{n/(p+\log n)}$⁠. Since the Huber loss is convex, the first (vector) equation in θ with s fixed can be solved using either the iteratively reweighted least squares algorithm or the Broyden–Fletcher–Goldfarb–Shanno algorithm. For the second equation with θ fixed, it can be shown that the function $s\mapsto g_2(\theta ,s)$ is non-increasing, as demonstrated by its derivative

$$\frac{\mathrm{\partial }}{\mathrm{\partial }s}{g}_2(\theta ,s)=-\frac{2}{n{s}^3}{\displaystyle \sum _{i=1}^n}({\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta {)}^2\mathbb{1}\{|{\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta |\le sk\}\le 0.$$

Proposition 3 in L. Wang et al. (2021) further guarantees that the equation $g_2(\theta ,s)=0$ with θ fixed has a unique solution, provided that $\sum _{i=1}^n\mathbb{1}(|{\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta |>0)>n/k^2=p+\log n$⁠. Based on these observations, we propose the following alternating algorithm, which begins at iteration 0 with an initial estimate ${\theta }^0=\hat{\theta }$⁠, the two-step LSE given in equation (11), or equivalently equation (12). At each iteration $t=1,2,\dots$⁠, the procedure involves two steps:

• Compute the ES ‘residuals’ ${\omega }_i^t={\hat{Z}}_i-\alpha X_i^{\mathrm{T}}{\theta }^{t-1}$ using the previous estimate ${\theta }^{t-1}$⁠. Let $s^t$ be the solution to the equation $(1/n)\sum _{i=1}^n(|{\omega }_i^t/s|\wedge k{)}^2=1$⁠, $s>0$⁠.

• Compute the updated estimate ${\theta }^t\in {\mathrm{argmin}}_{\theta \in {\mathbb{R}}^p}\sum _{i=1}^n{\ell }_{{\tau }^t}({\hat{Z}}_i-\alpha X_i^{\mathrm{T}}\theta )$⁠, where ${\tau }^t=s^{t}k$⁠.

Given a prespecified tolerance $ϵ>0$ (e.g. $ϵ={10}^{-5}$⁠), the algorithm will terminate at the tth iteration if $max\{\Vert g_1({\theta }^t,s^t){\Vert }_2,|g_2({\theta }^t,s^t)|\}\le ϵ$⁠, or if the maximum number of iterations is reached. Our numerical experiments in Section 6 show that this algorithm generally achieves convergence after only a small number of iterations. Intuitively, we attribute the algorithm’s fast convergence to the observation that ${\hat{\theta }}_{\tau }$ changes gradually as τ varies. This graduate change cause the residuals to behave similarly over a range of τ values. We refer the reader to the Online Supplementary Section B for a detailed elaboration on the motivations behind our proposed data-driven method.

4 Statistical theory

This section presents non-asymptotic high probability bounds for the error $\Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_2$ of the Huberised two-step ES estimator ${\hat{\theta }}_{\tau }$⁠, as defined in equation (19). Additionally, we establish a non-asymptotic Bahadur representation for ${\hat{\theta }}_{\tau }$⁠, which is a crucial step towards obtaining a Berry–Esseen-type bound for Gaussian approximation. Throughout this section, we write $X=(x_1,\dots ,x_p{)}^{\mathrm{T}}\in {\mathbb{R}}^p$ with $x_1\equiv 1$⁠. Without loss of generality, we assume that the random predictors $x_2,\dots ,x_p$ have zero means, that is, ${\mu }_j=\mathbb{E}(x_j)=0$ for $j=2,\dots ,p$⁠. This makes the later sub-Gaussian assumption more reasonable; see Condition 2 below. Otherwise, we set $Z=(1,z_2,\dots ,z_p{)}^{\mathrm{T}}=(1,x_2-{\mu }_2,\dots ,x_p-{\mu }_p{)}^{\mathrm{T}}$⁠. With this notation, the joint model (4) becomes $Q_{\alpha }(Y|Z)={\beta }_0^{†}+\sum _{j=2}^p{z}_j{\beta }_j^{*}$ and ${\mathrm{ES}}_{\alpha }(Y|Z)={\theta }_0^{†}+\sum _{j=2}^p{z}_j{\theta }_j^{*}$⁠, where ${\beta }_1^{†}={\mu }^{\mathrm{T}}{\beta }^{*}$ and ${\theta }_1^{†}={\mu }^{\mathrm{T}}{\theta }^{*}$ with $\mu =(1,{\mu }_2,\dots ,{\mu }_p{)}^{\mathrm{T}}$⁠. The sub-Gaussian assumption can then be imposed on Z, and our analysis naturally applies to $\{(Y_i,Z_i){\}}_{i=1}^n$⁠.

In the context of a joint (linear) quantile and ES regression model, we initiate by establishing a high probability bound, explicitly dependent on n and p, for the QR estimator $\hat{\beta }$ (9). To this end, we impose the following conditions on the covariates and the conditional distribution of Y given X.

 

Condition 1

The conditional density function of $\epsilon :=Y-X^{\mathrm{T}}{\beta }^{*}$ given X, denoted by $f_{\epsilon |X}$⁠, exists and is continuous on its support. Moreover, there exist constants $\underset{\_}{f},l_0>0$ such that $f_{\epsilon |X}(0)\ge \underset{\_}{f}$ and $|f_{\epsilon |X}(t)-f_{\epsilon |X}(0)|\le l_0|t|$ for all $t\in \mathbb{R}$ almost surely (over X).

 

Condition 2

The random covariate vector $X\in {\mathbb{R}}^p$ is sub-Gaussian, that is, there exists some (dimension-free) constant ${\upsilon }_1\ge 1$ such that $\mathbb{P}(|u^{\mathrm{T}}W|\ge {\upsilon }_{1}t)\le 2e^{-t^2/2}$ for all $t\ge 0$ and $u\in {\mathbb{S}}^{p-1}$⁠, where $W={\mathrm{\Sigma }}^{-1/2}X$ and $\mathrm{\Sigma }=\mathbb{E}(X{X}^{\mathrm{T}})$ is positive definite. Let ${\kappa }_l=\underset{u\in {\mathbb{S}}^{p-1}}{sup}\mathbb{E}|u^{\mathrm{T}}W{|}^l$ for $l\ge 1$⁠.

Condition 1 imposes regularity conditions on the random error distributions, accommodating heteroskedastic error distributions and not requiring the existence of any moment. Condition 2 is used to guarantee that population and empirical quantities (e.g. the objective or gradient function or the gradient function) are uniformly close to each other in a compact region. It can be replaced by a boundedness assumption, which will lead to similar results. For example, $X=(x_1,\dots ,x_p{)}^{\mathrm{T}}$ is compactly supported with either $\Vert X{\Vert }_{\mathrm{\infty }}\le C_X$ or $\Vert {\mathrm{\Sigma }}^{-1/2}X{\Vert }_2\le B_X$⁠, where $C_X$ is an absolute constant and $B_X$ is usually proportional to $\sqrt{p}$⁠.

 

Proposition 3

Under Conditions 1 and 2, the QR estimator $\hat{\beta }$ given in equation (9) satisfies, for any $t\ge 0$⁠, that $\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\le C_1{\underset{\_}{f}}^{-1}\sqrt{(p+t)/n}$ holds with probability at least $1-e^{-t}$ as long as $n\ge C_2{l}_0^2{\underset{\_}{f}}^{-4}(p+t)$⁠, where $C_1,C_2>0$ are constants depending only on ${\upsilon }_1$⁠.

While QR has been extensively studied since the seminal work of Koenker and Bassett (1978), there remains a paucity of literature that addresses its finite-sample properties, particularly in terms of high probability bounds. Proposition 3 revisits Theorem 2.1 originally presented in Pan and Zhou (2021). For the sake of completeness, we provide a self-contained and simplified proof in the Online Supplementary Section G.9. Shifting our focus to ES regression, which involves conditional expectations, we additionally impose the following moment condition on the random error $\epsilon$⁠.

 

Condition 3

The conditional CDF $F_{\epsilon |X}$ of $\epsilon$ given X is continuously differentiable and satisfies $|F_{\epsilon |X}(t)-F_{\epsilon |X}(0)|\le \overline{f}|t|$ for all $t\in \mathbb{R}$⁠. Moreover, the negative part of $\epsilon$⁠, denoted by ${\epsilon }_{-}=\epsilon \wedge 0$⁠, satisfies ${\mathrm{var}}_X({\epsilon }_{-})\le {\overline{\sigma }}^2$ almost surely (over X), where ${\mathrm{var}}_X$ denotes the conditional variance given X.

Condition 3 asserts that the conditional variance of the negative part of the QR residual $\epsilon =Y-X^{\mathrm{T}}{\beta }^{*}$ is bounded. In our theoretical analysis, we assume $\overline{\sigma }$ to be a constant for convenience. More generally, one can assume a form of $\epsilon =\sigma (X)\eta$⁠, where $\sigma :{\mathbb{R}}^p\to (0,\mathrm{\infty })$ is a positive function on ${\mathbb{R}}^p$ (not necessarily bounded), and η is independent of X satisfying $\mathrm{var}(\eta \mathbb{1}(\eta \le 0))\le {\overline{\sigma }}^2$⁠. In this case, an additional moment assumption on $\sigma (X)$⁠, such as boundedness $\mathbb{E}\{\sigma {(X)}^4\}$⁠, would suffice.

Our next result establishes high probability bounds for the estimation error of ES regression, conditioning on the event that $\hat{\beta }$ falls within a local neighbourhood of ${\beta }^{*}$⁠.

 

Theorem 1

Assume Conditions 2 and 3 hold. For any $t>0$⁠, let $r_0>0$ be such that $r_0\lesssim \overline{\sigma }$ and $\overline{f}{r}_0^2\lesssim \overline{\sigma }\sqrt{(p+t)/n}$⁠. Then, the two-step robust α-ES (⁠$0<\alpha \le 1/2$⁠) estimator ${\hat{\theta }}_{\tau }$ with $\tau =c_0\overline{\sigma }\sqrt{n/(p+t)}$ (for any $c_0\ge 1$⁠) satisfies that, with probability at least $1-3e^{-t}$ conditioned on the event $\{\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\le r_0\}$⁠,

$$\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\le C_1\overline{\sigma }\sqrt{\frac{p+t}{n}}+C_2(\sqrt{\frac{p+t}{n}}{r}_0+\overline{f}{r}_0^2)$$

(23)

provided that the sample size obeys $n\ge C_3(p+t)$⁠, where $C_1>0$ is a constant depending on $({\upsilon }_1,c_0)$ and $C_2,C_3>0$ depend only on ${\upsilon }_1$⁠.

For any $\delta \in (0,1)$⁠, the robust estimator ${\hat{\theta }}_{\tau }$ with $\tau \asymp \overline{\sigma }\sqrt{n/(p+\log (1/\delta ))}$ satisfies with probability at least $1-\delta$ conditioned on $\{\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\le r_0\}$ that

$$\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\lesssim \overline{\sigma }\sqrt{\frac{p+\log (1/\delta )}{n}}+\sqrt{\frac{p+\log (1/\delta )}{n}}{r}_0+\overline{f}{r}_0^2.$$

The above bound is proportional to $\log (1/\delta )$⁠, in contrast to the bound for the two-step LSE, which is proportional to $1/\delta$⁠, as demonstrated in the Online Supplementary Theorem A.1 in Section A. This observation suggests that the Huberised estimator is much more robust to heavy tails from a non-asymptotic perspective, compared to the two-step LSE. Specifically, in cases where the error variables only have finite variance, the worst-case deviations of $\hat{\theta }$ are considerably larger than those of ${\hat{\theta }}_{\tau }$⁠.

Bias-robustness trade-off

 

Remark 2

The choice of τ stated in Theorem 1 is a reflection of the bias-robustness trade-off. As discussed in Section 3.2, the robust estimator ${\hat{\theta }}_{\tau }$ can be viewed as an M-estimator of ${\theta }_{\tau }^{*}={\mathrm{argmin}}_{\theta }\mathbb{E}\{{\ell }_{\tau }(Z_i-\alpha X_i^{\mathrm{T}}\theta )\}$⁠, which differs from the true ES regression coefficient ${\theta }^{*}$ due to the asymmetry of ES ‘residuals’ ${\omega }_i=Z_i-\alpha X_i^{\mathrm{T}}{\theta }^{*}$⁠. Consider the decomposition

$$\Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\le {\displaystyle \underset{\mathrm{robustification}\mathrm{bias}}{{\displaystyle \underbrace{\Vert {\hat{\theta }}_{\tau }-{\theta }_{\tau }^{*}{\Vert }_{\mathrm{\Sigma }}}}}}+{\displaystyle \underset{\mathrm{robust}\mathrm{estimation}\mathrm{error}}{{\displaystyle \underbrace{\Vert {\theta }_{\tau }^{*}-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}}}}}.$$

As long as $\tau \gtrsim \overline{\sigma }$ under Condition 3, Proposition 2 ensures that $\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }_{\tau }^{*}{\Vert }_{\mathrm{\Sigma }}\le 2{\overline{\sigma }}^2/\tau$⁠. Examining the proof of Theorem 1, we see that

$$\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\lesssim \overline{\sigma }\sqrt{\frac{p+t}{n}}+\tau \frac{p+t}{n}+\frac{{\overline{\sigma }}^2}{\tau }+r_0(\sqrt{\frac{p+t}{n}}+\frac{\overline{\sigma }}{\tau })+r_0^2$$

with high probability conditioned on the event $\{\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\le r_0\}$⁠. We therefore select $\tau \asymp \overline{\sigma }\sqrt{n/(p+t)}$ in order to minimise the upper bound as a function of τ.

A uniform bound over τ

 

Remark 3

Recall from Proposition 3 that with probability at least $1-n^{-1}$⁠, $\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\lesssim {\underset{\_}{f}}^{-1}\sqrt{(p+\log n)/n}$ as long as $n\gtrsim p+\log n$⁠. Complementing the proof of Theorem 1 with a discretisation argument, we can obtain a more general result that holds for a range of τ values. Specifically, let $\overline{\tau }\ge \underset{\_}{\tau }>0$ be such that $\overline{\sigma }\lesssim \underset{\_}{\tau }\le \overline{\tau }\lesssim \overline{\sigma }\sqrt{n/(p+\log n)}$⁠. Then, with probability at least $1-C{n}^{-1}$ for some absolute constant $C\ge 1$⁠,

$${\displaystyle \underset{\tau \in [{\displaystyle \underset{\_}{\tau }},\overline{\tau }]}{sup}}\alpha \Vert {\hat{\theta }}_{\tau }-{\theta }^{*}{\Vert }_{\mathrm{\Sigma }}\lesssim \overline{\sigma }\sqrt{\frac{p+\log n}{n}}+\frac{{\overline{\sigma }}^2}{{\displaystyle \underset{\_}{\tau }}}+max\{\overline{f},1/{\displaystyle \underset{\_}{\tau }}\}\frac{p+\log n}{{{\displaystyle \underset{\_}{f}}}^{2}n},$$

(24)

as long as $n\gtrsim p+\log n$⁠. The proof of the uniform upper bound in equation (24) is provided in the Online Supplementary Section G.9. As ensured by this uniform bound, a data-driven choice of τ within the aforementioned range can be used.

If, in addition to Condition 3, some higher order moment of ${\epsilon }_{-}$ is bounded, namely, ${\mathbb{E}}_X\{|{\epsilon }_{-}-{\mathbb{E}}_X({\epsilon }_{-}){|}^k\}\le {\alpha }_k$ almost surely (over X) for some $k>2$⁠, the second term on the right-hand side of equation (24) will become ${\alpha }_k{\underset{\_}{\tau }}^{1-k}$⁠. In order to attain tight (finite-sample) concentration bounds, the robustification parameter $\tau =\tau (n,p)$ should not exceed $\sqrt{n/(p+\log n)}$ in magnitude. Conversely, τ should demonstrate sufficiently rapid growth in order for the bias term, controlled by ${\overline{\sigma }}^2{\underset{\_}{\tau }}^{-1}$ or ${\alpha }_k{\underset{\_}{\tau }}^{1-k}$ (in case higher order moments of ${\epsilon }_{-}$ are bounded), to decay at a comparable rate to the stochastic error.

Unlike the two-step LSE $\hat{\theta }$⁠, the robust counterpart ${\hat{\theta }}_{\tau }$ does not possess a closed-form expression. As a pivotal step in deriving Gaussian approximation results, the following theorem furnishes a non-asymptotic Bahadur representation for ${\hat{\theta }}_{\tau }$⁠, complete with explicit error bounds depending on $(n,p)$ and the first-stage QR estimation error.

 

Theorem 2

Assume the same conditions as in Theorem 1. For any $t>0$⁠, the α-ES estimator ${\hat{\theta }}_{\tau }$ with $\tau \asymp \overline{\sigma }\sqrt{n/(p+t)}$ satisfies that, with probability at least $1-6e^{-t}$ conditioned on $\{\Vert \hat{\beta }-{\beta }^{*}{\Vert }_{\mathrm{\Sigma }}\le r_0\}$⁠,

$${\Vert \alpha {\mathrm{\Sigma }}^{1/2}({\hat{\theta }}_{\tau }-{\theta }^{*})-\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\psi }_{\tau }({\omega }_i){\mathrm{\Sigma }}^{-1/2}{X}_i\Vert }_2\lesssim \overline{\sigma }\frac{p+t}{n}+\overline{f}{r}_0^2+r_0\sqrt{\frac{p\log n+t}{n}}$$

(25)

as long as $n\gtrsim p+t$⁠, where ${\omega }_i={\epsilon }_i\wedge 0+\alpha X_i^{\mathrm{T}}({\beta }^{*}-{\theta }^{*})$⁠.

Lastly, we present the following Gaussian approximation result that bounds the Kolmogorov distance between the distribution of the standardised statistic $\alpha \sqrt{n}{a}^{\mathrm{T}}({\hat{\theta }}_{\tau }-{\theta }^{*})/{\varrho }_a$ and the standard normal distribution, uniformly over all deterministic vectors $a\in {\mathbb{R}}^p$⁠, where ${\varrho }_a^2=a^{\mathrm{T}}{\mathrm{\Sigma }}^{-1}\mathrm{\Omega }{\mathrm{\Sigma }}^{-1}a$⁠. A similar conclusion applies to the oracle robust estimate ${\hat{\theta }}_{\tau }^{\mathrm{ora}}$ (20). The following theorem shows that the two-step robust estimator obtained via equations (9) and (19) is asymptotically equivalent to the oracle Huberised estimator (20), assuming ${\beta }^{*}$ is known.

 

Theorem 3

In addition to Conditions 1–3, assume that there exist constants $\overline{\sigma },{\alpha }_3>0$ such that

$${\mathrm{var}}_X({\epsilon }_{-})\ge {{\displaystyle \underset{\_}{\sigma }}}^2\text{ and }{\mathbb{E}}_X\{|{\epsilon }_{-}-{\mathbb{E}}_X({\epsilon }_{-}){|}^3\}\le {\alpha }_3\text{almost surely over}X.$$

(26)

Then, the robust α-level $(\alpha \in (0,1/2])$ ES estimator ${\hat{\theta }}_{\tau }$ with $\tau \asymp \overline{\sigma }\sqrt{n/(p+\log n)}$ satisfies

$$\begin{array}{rl}& {\displaystyle \underset{a\in {\mathbb{R}}^p,t\in \mathbb{R}}{sup}}|\mathbb{P}(\alpha \sqrt{n}⟨a,{\hat{\theta }}_{\tau }-{\theta }^{*}⟩/{\varrho }_a\le t)-\mathrm{\Phi }(t)|\\ & \lesssim \frac{{\alpha }_3}{{{\displaystyle \underset{\_}{\sigma }}}^3}\sqrt{\frac{p+\log n}{n}}+(\overline{f}/{{\displaystyle \underset{\_}{f}}}^2\vee {\alpha }_3^{1/3})\frac{p\sqrt{\log n}+\sqrt{p}\log n}{{\displaystyle \underset{\_}{\sigma }}\sqrt{n}}.\end{array}$$

(27)

Moreover, the oracle Huberised ES estimator ${\hat{\theta }}_{\tau }^{\mathrm{ora}}$ (20) with the same τ satisfies

$${\displaystyle \underset{a\in {\mathbb{R}}^p,t\in \mathbb{R}}{sup}}|\mathbb{P}(\alpha \sqrt{n}⟨a,{\hat{\theta }}_{\tau }^{\mathrm{ora}}-{\theta }^{*}⟩/{\varrho }_a\le t)-\mathrm{\Phi }(t)|\lesssim \frac{{\alpha }_3}{{{\displaystyle \underset{\_}{\sigma }}}^3}\sqrt{\frac{p+\log n}{n}}.$$

(28)

The above Gaussian approximation result lays the theoretical foundation for the statistical inference problems of testing the linear hypothesis $H_0:a^{\mathrm{T}}{\theta }^{*}=c_0$ vs. $H_1:a^{\mathrm{T}}{\theta }^{*}\ne c_0$ and constructing confidence intervals for $a^{\mathrm{T}}{\theta }^{*}$⁠, where $a\in {\mathbb{R}}^p$ and $c_0\in \mathbb{R}$ are predetermined. Given the joint quantile and ES regression estimates $(\hat{\beta },{\hat{\theta }}_{\tau })$⁠, let ${\hat{\mathrm{\Omega }}}_{\gamma }$ be the truncated estimator of $\mathrm{\Omega }=\mathbb{E}({\omega }^{2}X{X}^{\mathrm{T}})$ defined in equation (21) with $\gamma =\gamma (n,p)>0$ denoting a second robustification parameter. Then, we consider the robust test statistic $T_a=\alpha \sqrt{n}(a^{\mathrm{T}}{\hat{\theta }}_{\tau }-c_0)/{\hat{\varrho }}_{a,\gamma }$ for testing $H_0:a^{\mathrm{T}}{\theta }^{*}=c_0$⁠, and the (approximate) 100$(1-c)$% confidence interval $a^{\mathrm{T}}{\hat{\theta }}_{\tau }\pm z_{c/2}{\hat{\varrho }}_{a,\gamma }/(\alpha \sqrt{n})$ for $a^{\mathrm{T}}{\theta }^{*}$⁠, where ${\hat{\varrho }}_{a,\gamma }^2:=a^{\mathrm{T}}{\hat{\mathrm{\Sigma }}}^{-1}{\hat{\mathrm{\Omega }}}_{\gamma }{\hat{\mathrm{\Sigma }}}^{-1}a$ is a robust variance estimator and $z_{c/2}$ is the upper $(c/2)$-percentile of $\mathcal{N}(0,1)$⁠. The consistency of ${\hat{\varrho }}_{a,\gamma }^2$ with a properly chosen γ is investigated in the Online Supplementary Section C.

5 Non-parametric expected shortfall regression

In this section, we consider non-parametric models for joint quantile and expected shortfall regression. For a predetermined quantile level $\alpha \in (0,1)$⁠, the goal is to estimate the unknown (conditional) quantile and expected shortfall functions $f_q^{*}(x)=Q_{\alpha }(Y|X=x)$ and $f_e^{*}(x)={\mathrm{ES}}_{\alpha }(Y|X=x)$⁠, with an emphasis on the latter. By equation (1), $f_q^{*}$ and $f_e^{*}$ can be identified as

$$f_q^{*}={\displaystyle \underset{f_q}{argmin}}\mathbb{E}{\rho }_{\alpha }(Y-f_q(X))\text{and}{f}_e^{*}={\displaystyle \underset{f_e}{argmin}}\mathbb{E}\{Y-f_e(X){\}}^2{\mathbb{1}}_{\{Y\le f_q^{*}(X)\}}.$$

Motivated by the two-step procedure developed under joint linear models, in the following we propose a non-parametric ES estimator using the series regression method (Andrews, 1991; Eubank & Spiegelman, 1990; Newey, 1997). Such a non-parametric estimate is carried out by regressing the dependent variable on an asymptotically growing number of approximating functions of the covariates, and therefore is closely related to the estimator define in equation (11) under the so-called many regressors model (Belloni et al., 2019), that is, the dimension $p=p_n$ is allowed to grow with n. The idea of series estimation is to first approximate $f_q^{*}$ and $f_e^{*}$ by their ‘projections’ on the linear spans of $m_1$ and $m_2$ series/basis functions, respectively, and then fit the coefficients using the observed data. Specifically, we approximate functions $f_q^{*}$ and $f_e^{*}$ by linear forms $U{(x)}^{\mathrm{T}}\beta$ and $V{(x)}^{\mathrm{T}}\theta$⁠, where

$$U(x)=(u_1(x),\dots ,u_{m_1}(x){)}^{\mathrm{T}}\text{and}V(x)=(v_1(x),\dots ,v_{m_2}(x){)}^{\mathrm{T}}$$

are two vectors of series approximating functions of dimensions $m_1$ and $m_2$⁠. Here both $m_1$ and $m_2$ may increase with n. We thus define the vectors of quantile and ES series approximation coefficients as

$${\beta }^{*}\in {\displaystyle \underset{\beta \in {\mathbb{R}}^{m_1}}{argmin}}\mathbb{E}{\rho }_{\alpha }(Y-U{(X)}^{\mathrm{T}}\beta )\text{and}{\theta }^{*}\in {\displaystyle \underset{\theta \in {\mathbb{R}}^{m_2}}{argmin}}\mathbb{E}\{Y-V{(X)}^{\mathrm{T}}\theta {\}}^2{\mathbb{1}}_{\{Y\le f_q^{*}(X)\}}.$$

(29)

Given independent observations $(Y_i,X_i)$⁠, $1\le i\le n$ from $(Y,X)\in \mathbb{R}\times \mathcal{X}$ with $\mathcal{X}$ denoting a compact subset of ${\mathbb{R}}^p$⁠, we write $U_i=U(X_i)\in {\mathbb{R}}^{m_1}$ and $V_i=V_i(X_i)\in {\mathbb{R}}^{m_2}$⁠. Extending the two-step approach described in Section 2.2, we first define the (conditional) quantile series estimator of $f_q^{*}(x)=Q_{\alpha }(Y|X=x)$ (Belloni et al., 2019):

$${\hat{f}}_q(x)=U{(x)}^{\mathrm{T}}\hat{\beta },x\in \mathcal{X},\text{where}\hat{\beta }={\hat{\beta }}_{m_1}\in {\displaystyle \underset{\beta \in {\mathbb{R}}^{m_1}}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\rho }_{\alpha }(Y_i-U_i^{\mathrm{T}}\beta ).$$

(30)

With generated response variables ${\hat{Z}}_i=\alpha {\hat{f}}_q(X_i)+\{Y_i-{\hat{f}}_q(X_i)\}{\mathbb{1}}_{\{Y_i\le {\hat{f}}_q(X_i)\}}$⁠, the second-stage ES series estimator is given by

$${\hat{f}}_e(x)=V{(x)}^{\mathrm{T}}\hat{\theta },x\in \mathcal{X},\text{where}\hat{\theta }={\hat{\theta }}_{m_2}\in {\displaystyle \underset{\theta \in {\mathbb{R}}^{m_2}}{argmin}}\frac{1}{n}{\displaystyle \sum _{i=1}^n}({\hat{Z}}_i-\alpha V_i^{\mathrm{T}}\theta {)}^2.$$

(31)

Commonly used series functions with good approximation properties include B-splines, polynomials, Fourier series and compactly supported wavelets. We refer to Newey (1997) and Chen (2007) for a detailed description of these series functions. In the context of QR, Chen (2007) established the consistency and rate of convergence at a single quantile index. More recently, Belloni et al. (2019) developed a large sample theory for the quantile series coefficient process, including convergence rate and uniform strong approximations. The choice of the parameter $m_1$⁠, also known as the order of the series estimator, is crucial for establishing the balance between bias and variance.

Note that the quantile series estimator ${\hat{f}}_q$ in equation (30) has been well studied by Belloni et al. (2019). Because the number of regressors increases with the sample size, conventional central limit theorems are no longer applicable to capture the joint asymptotic normality of the regression coefficients. The growing dimensionality is the primary source of technical complications. Our theoretical analysis under the joint linear model (4), which leads to novel non-asymptotic high probability bounds, can be used as a starting point for studying the two-step non-parametric ES series estimator ${\hat{f}}_e$ defined in equation (31). Of particular interest is to develop a uniform inference procedure for the conditional ES function $f_e^{*}$ and its theory. That is, at a given confidence level $1-\gamma$⁠, we aim to construct a pair of functional estimates $[{\hat{f}}_e^L,{\hat{f}}_e^U]$ from $\{(Y_i,X_i){\}}_{i=1}^n$ such that

$$\mathbb{P}\{{\hat{f}}_e^L(x)\le f_e^{*}(x)\le {\hat{f}}_e^U(x)\text{ for all}x\in \mathcal{X}\}\to 1-\gamma ,\text{as}n\to \mathrm{\infty }.$$

Since a significant amount of additional work is still needed, including explicit characterisations of the ES series approximation error and the impact of first-stage non-parametric QR estimation error, we leave a rigorous theoretical investigation of ${\hat{f}}_e$ to future work. Although we have only focussed on series methods, there are other non-parametric techniques that offer superior empirical and theoretical performance. Among those, deep neural networks have stood out as a promising tool for non-parametric estimation, from least squares, logistic to QR (Farrell et al., 2021; Schmidt-Hieber, 2020; Shen et al., 2021). It is practically useful to construct deep learning implementations of two-step estimators and statistically important to deliver valid inferences on finite-dimensional parameters following first-step estimation (of both quantile and ES functions) using deep learning. A detailed investigation of these problems is beyond the present scope but of future interest.

6 Numerical studies and real data examples

6.1 Monte Carlo experiments

In this section, we assess the numerical performance of the proposed method for fitting expected shortfall regression. For its R implementation, we first obtain a QR estimate via the quantreg library, and in step two use the adaHuber library to solve (19) with the robustification parameter selected adaptively as described in Section 3.2.

We compare the proposed two-step adaptive Huber ES estimator (2S-AH) to several competitors: (i) the joint regression estimate (joint) via FZ loss minimisation, implemented via the R library esreg with the default option; (ii) the two-step LSE (12) (2S-LS); and (iii) the oracle two-step ‘estimator’ (2S-oracle). Recall that the two-step procedure first obtains a QR estimator $\hat{\beta }$ via either standard (Koenker & Bassett, 1978) or smoothed QR regression (He et al., 2023), and subsequently computes the ES estimator based on fitted quantile thresholds $\{X_i^{\mathrm{T}}\hat{\beta }{\}}_{i=1}^n$⁠. The oracle method refers to the two-step ES estimate based on the true quantile thresholds $\{X_i^{\mathrm{T}}{\beta }^{*}{\}}_{i=1}^n$⁠.

In our simulation studies, we first generate ${\gamma }^{*}=({\gamma }_1^{*},\dots ,{\gamma }_p^{*}{)}^{\mathrm{T}}$ and ${\eta }^{*}=({\eta }_1^{*},\dots ,{\eta }_p^{*}{)}^{\mathrm{T}}$ independently, where ${\gamma }_j^{*}$s are independent Rademacher random variables and ${\eta }_j^{*}{\sim }_{\mathrm{i}.\mathrm{i}.\mathrm{d}.}0.5\cdot \mathrm{Bernoulli}(1/2)$⁠. Data are then generated from the heteroscedastic model

$$Y_i=X_i^{\mathrm{T}}{\gamma }^{*}+X_i^{\mathrm{T}}{\eta }^{*}\cdot {\epsilon }_i,$$

(32)

where $X_i=(X_{i1},\dots ,X_{ip}{)}^{\mathrm{T}}$ with $X_{ij}{\sim }_{\mathrm{i}.\mathrm{i}.\mathrm{d}.}\mathrm{Unif}(0,1.5)$⁠, and the random noise ${\epsilon }_i$ follows one of the following two distributions: (i) standard normal distribution and (ii) t-distribution with $\nu >2$ degrees of freedom (⁠$t_{\nu }$⁠). Given ${\gamma }^{*}$ and ${\eta }^{*}$⁠, the true quantile and expected shortfall regression coefficients are ${\beta }^{*}={\gamma }^{*}+Q_{\alpha }(\epsilon )\cdot {\eta }^{*}$ and ${\theta }^{*}={\gamma }^{*}+{\mathrm{ES}}_{\alpha }(\epsilon )\cdot {\eta }^{*}$⁠, where $Q_{\alpha }(\epsilon )$ and ${\mathrm{ES}}_{\alpha }(\epsilon )$ are the α-level quantile and expected shortfall of $\epsilon$⁠, respectively.

We first set the dimension $p=20$ and sample size $n=\lceil 50p/\alpha \rceil$⁠, where the quantile level α takes values in $\{0.05,0.1,0.2\}$⁠. Simulation results on the relative ${\ell }_2$-error $\Vert \hat{\theta }-{\theta }^{*}{\Vert }_2/\Vert {\theta }^{*}{\Vert }_2$⁠, averaged over 200 replications, are reported in Tables 1 and 2 under the $\mathcal{N}(0,1)$ and $t_{2.5}$ noise model, respectively. All four methods have very similar performance across different quantile levels in the normal model, while in the presence of heavy-tailed errors, the proposed robust method achieves consistently more favourable performance. This demonstrates that the use of adaptive Huber regression (in stage two) gains robustness against heavy-tailed errors without compromising statistical efficiency when the error distribution is light-tailed.

Table 1.

Mean relative ${\ell }_2$-error $\Vert \hat{\theta }-{\theta }^{*}{\Vert }_2/\Vert {\theta }^{*}{\Vert }_2$ (and standard error), averaged over 200 replications, when ${\epsilon }_i\sim \mathcal{N}(0,1)$⁠, $p=20$⁠, $n=\lceil 50p/\alpha \rceil$⁠, and $\alpha =\{0.05,0.1,0.2\}$

	$\mathcal{N}(0,1)$ noise
Method	$\alpha =0.05$	$\alpha =0.1$	$\alpha =0.2$
2S-AH	0.130 (0.003)	0.150 (0.003)	0.171 (0.004)
2S-LS	0.130 (0.003)	0.150 (0.003)	0.171 (0.004)
joint	0.130 (0.003)	0.151 (0.003)	0.177 (0.004)
2S-oracle	0.129 (0.003)	0.149 (0.003)	0.171 (0.004)

Table 2.

Mean relative ${\ell }_2$-error $\Vert \hat{\theta }-{\theta }^{*}{\Vert }_2/\Vert {\theta }^{*}{\Vert }_2$ (and standard error), averaged over 200 replications, when ${\epsilon }_i\sim t_{2.5}$⁠, $p=20$⁠, $n=\lceil 50p/\alpha \rceil$ and $\alpha =\{0.05,0.1,0.2\}$

	$t_{2.5}$ noise
Method	$\alpha =0.05$	$\alpha =0.1$	$\alpha =0.2$
2S-AH	0.484 (0.008)	0.470 (0.009)	0.429 (0.008)
2S-LS	0.612 (0.013)	0.606 (0.016)	0.532 (0.013)
joint	0.581 (0.012)	0.567 (0.014)	0.511 (0.013)
2S-oracle	0.612 (0.013)	0.607 (0.016)	0.532 (0.013)

In a more extreme setting where $\alpha =0.01$⁠, Figure 1 shows the boxplots of squared ${\ell }_2$-errors for three ES estimates (2S-LS, 2S-AH, and joint) under the normal and $t_3$ models. Although the 2S-LS estimator is easy-to-compute, it is more sensitive to heavy-tailedness than the joint estimator obtained via FZ loss minimisation. We further compare the proposed method with the joint regression approach in terms of computational efficiency. The computational time in seconds averaged over 50 independent replications, for the two methods with growing $(n,p)$ subject to $n=\lceil 50p/\alpha \rceil$ (⁠$\alpha \in \{0.05,0.1,0.2\}$⁠) are reported in Figure 2. These numerical results show evidence that our R implementation of the robust two-step method can be faster than the esreg library for the joint regression approach by several orders of magnitude.

Figure 1.

Boxplots of squared total ${\ell }_2$-errors (including the intercept when its true value is 2), based on 500 replications, for three ES regression estimators (2S-LS, 2S-AH, and joint) at quantile level $\alpha =0.01$⁠. The mean squared errors of these three estimators are 0.1219, 0.0983, and 0.1119 in the normal model, and 1.7401, 0.8017, and 1.2542 in the $t_3$ model.

Figure 2.

Average elapsed time (in seconds) over 50 replications for the proposed method implemented by a combination of quantreg and adaHuber and the joint regression approach implemented by esreg under $\mathcal{N}(0,1)$ and $t_{2.5}$ error models when $\alpha \in \{0.05,0.1,0.2\}$⁠. The sample size is set to be $n=\lceil 50p/\alpha \rceil$⁠. The solid and dashed lines correspond to the proposed method and the joint regression approach, respectively.

To shed some light on the drastic difference in numerical efficiency between the two methods, note that the joint regression approach (Dimitriadis & Bayer, 2019) relies on the Nelder–Mead simplex method, which is sensitive to the starting values and not guaranteed to converge to a local minimum. The convergence of the Nelder–Mead method is already very slow for large-scale problems because it is a direct search method based on function comparison. And due to its sensitivity to starting values, Dimitriadis and Bayer (2019) proposed to re-optimise the model (several times) with the perturbed parameter estimates as new starting values. This explains, to some extent, the fast increase in the runtime of esreg as both n and p grow. The function in quantreg that fits linear QR is coded in fortran, and thus is very fast in larger problems. The computation of adaptive Huber regression is based on the Barzilai-Borwein gradient descent method (Barzilai & Borwein, 1988), implemented via RcppArmadillo in adaHuber.

Next, we construct entrywise (approximate) 95% confidence intervals (CIs) for the expected shortfall regression parameter ${\theta }^{*}$⁠. The CI for the two-step estimator is based on equation (14) (non-robust) and equation (22) (robust), and we use the default option in the esreg package to implement Dimitriadis and Bayer (2019)’s method. To evaluate the accuracy and reliability of the CIs, we compute the empirical coverage probability and interval width based on 500 independent replications, then averaged over the p slope coefficients. Results for $p=20$ and $n=\lceil 50p/\alpha \rceil$ (⁠$\alpha \in \{0.05,0.1,0.2\}$⁠) are reported in Tables 3 and 4.

Table 3.

Empirical coverage probability and mean width (based on 500 replications) of 95% confidence intervals averaged over $p=20$ variables when $n=\lceil 50p/\alpha \rceil$⁠, $\alpha =\{0.05,0.1,0.2\}$⁠, and ${\epsilon }_i\sim \mathcal{N}(0,1)$

$\mathcal{N}(0,1)$	$\alpha =0.05$		$\alpha =0.1$		$\alpha =0.2$
Method	Coverage	Width	Coverage	Width	Coverage	Width
2S-AH	0.950	0.595	0.949	0.660	0.948	0.744
joint	0.946	0.584	0.944	0.651	0.942	0.740
2S-LS	0.950	0.595	0.949	0.661	0.948	0.745

Table 4.

Empirical coverage probability and mean width (based on 500 replications) of 95% confidence intervals averaged over $p=20$ variables when $n=\lceil 50p/\alpha \rceil$⁠, $\alpha =\{0.05,0.1,0.2\}$⁠, and ${\epsilon }_i\sim t_{2.5}$

$t_{2.5}$	$\alpha =0.05$		$\alpha =0.1$		$\alpha =0.2$
Method	Coverage	Width	Coverage	Width	Coverage	Width
2S-AH	0.947	3.633	0.946	2.790	0.948	2.243
joint	0.959	5.771	0.959	3.571	0.954	2.872
2S-LS	0.952	4.521	0.950	3.397	0.953	2.687

Once again, all three methods perform similarly under normal errors, while the robust approach gives the narrowest CIs while maintaining the desired coverage level under $t_{2.5}$ errors. Together, the results in Tables 2 and 4 demonstrate the robustness of the proposed method, as indicated by the theoretical investigations in Section 4.

6.2 Data application I: health disparity

Iron deficiency is one of the most common nutritional deficiency worldwide and is one of the leading cause of anaemia (Camaschella, 2015). Being able to detect iron deficiency is essential in medical care for patients with inflammation, infection, or chronic disease. It is also important in preventive care since iron deficiency tends to present signs of a more serious illness such as gastrointestinal malignancy (Rockey & Cello, 1993). One measure of iron deficiency that has proven to be useful is the soluble transferrin receptor (sTRP), a carrier protein for transferrin (Mast et al., 1998). A high value of sTRP indicates iron deficiency.

The scientific goal here is to assess whether there is any disparity in sTRP levels among four different ethnic groups: Asian, Black, Mexican American, and White. To this end, we analyse a data set obtained from the National Health and Nutrition Examination Survey from 2017 to 2020 (pre-covid). In this data set, the response variable sTRP was measured for female participants who range in age from 20 to 49 years. The covariates of interest are three dummy variables that correspond to Asian, Mexican American, and Black, using White as the baseline. We adjust for demographic variables such as age, education level, and healthy diet throughout our analysis. For simplicity, we remove all participants with missing values on the covariates and the final data set consists of $n=1,689$ observations and $p=7$ covariates.

As an exploratory analysis, in Figure 3 we plot the quantile curves of sTRP measurements at levels from 50% to 99% for each of the four different ethnic groups. In this data set, the sTRP values range from 1.24 to 35.1 mg/L. We note that the normal range for females is between 1.9 and 4.4 mg/L (Kratovil et al., 2007), and values that are much higher than 4.4 mg/L indicate severe iron deficiency. We see from Figure 3 that the majority of the population have sTRP levels within the normal range. However, there are large disparities between Black and the other three ethnic groups, reflected in higher quantiles of the marginal distributions of sTRP.

Figure 3.

The soluble transferrin receptor levels (mg/L) vs. quantile levels (ranging from 0.5 to 0.99) for the female population in four different ethnic groups: Asian, Black, Mexican American, and White. The orange horizontal dashed line indicates the upper bound of the normal range (1.9–4.4 mg/L) for transferrin receptors among females.

To quantify the statistical significance of the aforementioned disparity, we fit robust expected shortfall regression at $\alpha =0.75$ (upper tail), with the robustification parameter tuned by the procedure described in Section 3.2. This is equivalent to fitting the proposed 2S-AH method at level $1-\alpha$ (see Section 3) after flipping the signs of both the response and the covariates. We also implement the standard QR at level α.

Table 5 reports the estimated coefficients and the associated 95% confidence intervals for the three indicator covariates on the ethnic groups Asian, Mexican American, and Black, using White as a baseline. We see that both the quantile and robust expected shortfall regression methods are able to detect a health disparity between Black and White. Specifically, the estimated robust ES regression coefficient and 95% CI (in the parenthesis) is 3.03 (1.88, 4.19) vs. its QR counterparts’ 0.86 (0.37, 1.35). With the use of QR (at level 0.75), we do not observe a statistically significant health disparity between Asian and White. In contrast, 2S-AH detects health disparity between Asian and White with an estimated coefficient 2.34 (0.59, 4.09). We also see that the QR detects health disparity between Mexican American and White, but the effect size is close to zero. In summary, ES regression complements QR, and can be more effective, as a tool to detect health disparity especially when it only occurs in the tail of the conditional distribution.

Table 5.

The estimated regression coefficients (and 95% confidence intervals) for three dummy variables: Asian, Black, and Mexican American, using White as a baseline

	Asian	Black	Mexican American
QR	0.31 (⁠$-$0.02, 0.64)	0.86 (0.37, 1.35)	$-$0.22 (⁠$-$0.42, -0.01)
ES regression (2S-AH)	2.34 (0.59, 4.09)	3.03 (1.88, 4.19)	0.13 (⁠$-$0.76, 1.03)

Note. Results of the upper-tail robust ES regression method 2S-AH and standard QR at quantile level $\alpha =0.75$ are reported.

6.3 Data application II: JTPA

We consider the JTPA study, a publicly funded training programme that provides training for adults with the goal of improving their earnings. Specifically, we focus on the Title II sub-programme of the JPTA study that is mainly offered to adults with barriers to employment and out-of-school youths. This data set was previously analysed in Bloom et al. (1997). It consists of 30 months of accumulated earnings for 6,102 females and 5,102 males, with 16 covariates that are related to the demographics of the individuals such as age, race, and the indicator variable that indicates whether the individual received JPTA training. After removing individuals with zero income, there are 4,576 males and 5,296 females. Our goal is to assess the effect of JPTA training on participants’ earnings with an emphasis on the low-income population that is employed, for both male and female sub-groups.

To this end, we fit an expected shortfall regression model using the proposed robust method with $\alpha =\{0.05,0.1,0.2\}$⁠. The robustification parameter τ is selected automatically via the procedure described in Section 3.2. Specifically, we regress the 30-month accumulated earnings on the JPTA training to assess the effect of JPTA training on low-income individuals, adjusting for whether individuals completed high school, race, Hispanic/non-Hispanic, marital status, working less than 13 weeks in the past year, and age. We report the estimated regression coefficient for the binary variable JPTA training and its associated 95% confidence intervals. The results are summarised in Table 6.

Table 6.

The estimated regression coefficient of the binary predictor JPTA training (and its 95% confidence interval) for the proposed robust method and the standard QR at quantile level $\alpha \in \{0.05,0.1,0.2\}$

Male sub-group	$\alpha =0.05$	$\alpha =0.1$	$\alpha =0.2$
QR	465 (255, 675)	882 (603, 1,161)	2031 (1,431, 2,603)
ES regression (2S-AH)	283 (149, 418)	552 (333, 771)	1,093 (641, 1,546)
Female sub-group	$\alpha =0.05$	$\alpha =0.1$	$\alpha =0.2$
QR	202 (76, 328)	480 (307, 653)	1,086 (719, 1,452)
ES regression (2S-AH)	123 (41, 205)	300 (146, 453)	672 (385, 958)

Note. Results are rounded to the closest integer.

From Table 6, we see that 95% confidence intervals for the robust method do not contain zero for all $\alpha \in \{0.05,0.1,0.2\}$⁠. This indicates that the JPTA training is statistically effective to improve earnings for the low-income population. Specifically, for the male sub-population, the estimated ES-effects of JPTA training are 283, 552, and 1,093 dollars at levels $0.05$⁠, $0.1$⁠, and $0.2$⁠, respectively. To further assess whether the estimated effects are scientifically meaningful, we compute the average 30-month accumulated earnings below the quantile levels $0.05,0.1$⁠, and $0.2$ for the male sub-group, which are 214, 566, and 1,496, respectively. We find that the JPTA training doubles the average income for individuals with income below the quantile levels 0.05 and 0.1, and becomes less effective for individuals with higher income. Similar findings are also observed for the female sub-group.

7 Conclusion and discussions

This paper considers expected shortfall regression under a joint quantile and ES model recently proposed in Dimitriadis and Bayer (2019) and Patton et al. (2019). The existing approach is based on a joint M-estimator, defined as the global minimum of any member of a class of strictly consistent joint loss functions (Fissler & Ziegel, 2016) over some compact set. Since the loss function is non-differentiable and non-convex, the computation of such a joint M-estimator is intrinsically difficult especially when the dimensionality is large. To circumvent the aforementioned challenge, Barendse (2020) proposed a two-step procedure for estimating the joint quantile and ES model based on Neyman orthogonalisation: the first step involves fitting the QR, and the second step employs the Neyman-orthogonal scores to estimate the ES parameters. Due to the use of $L_2$-loss in the second step, the resulting estimator is sensitive to heavy-tailed error distributions.

To address the robustness and computation concerns simultaneously, we propose a robust two-step method that applies adaptive Huber regression (Zhou et al., 2018) in the second step. The key is the use of a diverging robustification parameter for bias-robustness trade-off, tuned by a convenient data-driven mechanism. The proposed method can be efficiently implemented by a combination of R packages quantreg/conquer and adaHuber. The Python code that implements both our proposed methods and the existing non-convex optimisation-based methods (Dimitriadis & Bayer, 2019; Peng & Wang, 2022) is now publicly available at https://github.com/WenxinZhou/conquer. We establish a finite-sample theoretical framework for this two-step method, including deviation bound, Bahadur representation and (uniform) Gaussian approximations, in which the dimension of the model, p, may depend on and increase with the sample size, n. Robust confidence intervals/sets are also constructed. Numerical experiments further demonstrate that the proposed robust ES regression approach achieves satisfying statistical performance, high degree of robustness (against heavy-tailed data) and superior computational efficiency and stability. Through two data applications on health disparity and the JTPA study, we illustrate that ES regression complements QR as a useful tool to explore heterogeneous covariate effects on the average tail behaviour of the outcome.

Although we restrict attention to (joint) linear models in this work, our non-asymptotic theory and the underpinning techniques pave the way for analysing (i) series/projection estimators under joint non-parametric quantile-ES models and (ii) penalised estimators under high-dimensional sparse quantile-ES models. We leave these extensions in future research. One limitation in our data analysis for the JTPA study is that we do not account for potential selection bias. Specifically, as pointed out by Abadie et al. (2002), out of all subjects that were assigned to participate in the training programme, only approximately 60% of them (compliers) actually committed to the training programme. These individuals may simply have higher motivation in improving their earnings, and thus, the training status is likely positively correlated with potential income earnings. Generalising the proposed method to estimate the complier expected shortfall treatment effect, using an instrumental variable approach previously considered in Abadie et al. (2002), is another direction for future research.

The ES regression methods considered in this paper are suited for a fixed quantile level $\alpha \in (0,1)$⁠, independent of the sample size. For extreme quantiles satisfying $\alpha ={\alpha }_n\to$ 0 or 1 as $n\to \mathrm{\infty }$⁠, both the FZ loss minimisation method (see equations (5) and (6)) and two-step procedures perform poorly because observations become scarce at that level, i.e. $\alpha n$ is not large enough. In fact, if dimension p is fixed, Online Supplementary Theorem A.1 and Theorem 1 imply that the two-step ES regression estimates, robust, and non-robust, are consistent if ${\alpha }_n^{2}n\to \mathrm{\infty }$ as $n\to \mathrm{\infty }$⁠. In the case where ${\alpha }_n^{2}n=O(1)$⁠, these methods are no longer useful and one may need to resort to extreme value theory (de Haan & Ferreira, 2006; H. J. Wang et al., 2012), which provides the statistical tools for a feasible extrapolation into the tail of the variable of interest. A more detailed discussion on modelling the extremes is deferred to the Online Supplementary Section E.

Acknowledgments

We sincerely thank the Joint Editor, Qiwei Yao, as well as the Associate Editor and the three reviewers for their numerous constructive comments and valuable suggestions, which have significantly helped us improve the quality of this work.

Funding

X.H. is supported by NSF Grants DMS-1914496 and DMS-1951980. K.M.T. is supported by NSF DMS-2113346 and NSF CAREER DMS-2238428. W.-X.Z. acknowledges the support from NSF DMS-2113409.

Data availability

The two data sets used in Sections 6.2 and 6.3 are publicly available at https://www.cdc.gov/nchs/nhanes/index.htm and https://economics.mit.edu/people/faculty/josh-angrist/angrist-data-archive, respectively.

Supplementary material

Supplementary material is available online at Journal of the Royal Statistical Society: Series B.

Author notes

Conflict of interest: None declared.