Online Prediction of Extreme Conditional Quantiles via B-Spline Interpolation

Conditional quantile estimation is a statistical technique that focuses on modeling and analyzing conditional quantiles of a response variable concerning one or more predictor variables. The techniques involved in conditional quantile estimation include quantile regression and various non-parametric approaches. The motivation behind conditional quantile estimation is typically different from that of traditional mean-based modeling. While traditional regression techniques, like linear regression, focus on estimating the conditional mean of the response variable, conditional quantile methods seek to estimate quantiles, which offer a more comprehensive view of the conditional distribution. This approach is particularly useful when dealing with non-normally distributed data, outliers, or when the interest is in understanding the variability at different points in the distribution.

We consider here the setting where the response $Y\in\mathbb{R}$ depends on covariates $\mathbf{X}\in\mathbb{R}^{p}$. Denote the conditional distribution of the variable $Y$ as $F_{Y}(\cdot|\mathbf{X})$ and marginal distribution of $Y$ as $F_{Y}$. Suppose observations $(y_{i},\mathbf{x}_{i}),i=1,2,\ldots,N$, where $y_{i}$ is generated from $F_{Y}(\cdot|\mathbf{x}_{i})$ given $\mathbf{X}=\mathbf{x}_{i}$. The conditional quantile of $Y$ given $\mathbf{X}=\mathbf{x}$ at the quantile level $\tau\in(0,1)$ is defined as $Q_{Y}(\tau|\mathbf{x})=F_{Y}^{-1}(\tau|\mathbf{x})$, where $F_{Y}^{-1}(\cdot|\mathbf{x})$ is the generalized inverse of the conditional distribution function of $Y$. The conditional quantile $Q_{Y}(\tau|\mathbf{x})$ can be calculated as the solution to the following optimization problem

where $\rho_{\tau}(u)=(\tau-I(u\leq 0))u$ is the so-called check function (Koenker and Bassett Jr, 1978). Quantile regression is one of the most widely used techniques for conditional quantile estimation. It extends traditional regression by estimating the conditional quantiles directly. We follow the standard parametric assumption that the true conditional quantile has a linear form as

where $(\alpha(\tau),\boldsymbol{\beta}^{T}(\tau))$ is a parameter vector associated with the level $\tau$ and $(\alpha_{0}(\tau),\boldsymbol{\beta}_{0}^{T}(\tau))$ denotes the true parameter. With observations $(y_{i},\mathbf{x}_{i}),\,i=1,2,\ldots,N$, the parameter can be estimated by solving an optimization problem corresponding to the empirical counterpart of (2.1):

The conditional quantile estimation of $Y$ given $\mathbf{X}=\mathbf{x}$ is thus calculated as

The empirical linear estimation of the conditional quantile works effectively when the quantile level $\tau=\tau_{0}$ is moderately high. However, as the quantile level $\tau=\tau_{N}$ approaches $1$ and $N(1-\tau_{N})\rightarrow c,c\in[0,\infty)$, the expected number of observations exceeding $Q_{Y}(\tau_{N}|\mathbf{x})$ becomes insufficient. In such cases, empirical linear estimation of the conditional quantile is no longer feasible, necessitating extrapolation beyond observed data. Extrapolation is the process of estimating or predicting values beyond the range of the available data. It is a common practice in various fields, including finance, environmental science, insurance, and reliability analysis, where the primary interest lies in understanding and managing extreme events, such as financial market crashes, natural disasters, rare disease occurrences, or equipment failures. To address this need for extrapolation, we turn to the GPD, which offers an excellent approximation of the tail distribution of the conditional distribution of $Y$ given $\mathbf{X}=\mathbf{x}$. More specifically, given a large threshold $u$, $\mathbf{X}=\mathbf{x}$, and $z>0$, the exceedance $Y-u\mid Y>u$ given $\mathbf{X}=\mathbf{x}$ satisfies

where $\tau_{u}(\mathbf{x})=F_{Y}(u|\mathbf{x})\in(0,1)$ and we require $1+\gamma_{0}(\mathbf{x})z/\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))>0$ for $\gamma_{0}(\mathbf{x})\neq 0$. Here $\gamma_{0}(\mathbf{x})$ is the extreme value index, and $\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))$ is a scale related to the quantile level of the threshold $u$. When $\gamma_{0}(\mathbf{x})<0$, there is a finite right endpoint $u^{*}=u-\frac{\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))}{\gamma_{0}(\mathbf{x})}$ in the support of the distribution of $Y$, that is, $F_{Y}(y|\mathbf{X})=1$ for all $y\geq u^{*}$. When $\gamma_{0}(\mathbf{x})=0$, the exceedance $Y-u\mid Y>u$ given $\mathbf{X}=\mathbf{x}$ has an exponential distribution with mean $\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))$. When $\gamma_{0}(\mathbf{x})>0$, the exceedance $Y-u\mid Y>u$ given $\mathbf{X}=\mathbf{x}$ has a heavy tail with up to $\frac{1}{\gamma_{0}(\mathbf{x})}$-th finite moments. Note that for any higher threshold $u^{\prime}>u$, the exceedance $Y-u\mid Y>u$ given $\mathbf{X}=\mathbf{x}$ again follows the generalized Pareto distribution with the same shape parameter $\gamma_{0}(\mathbf{x})$ but a different scale parameter.

Once the threshold $u$ is set, the next step involves estimating the parameters of the GPD. In the absence of covariates, a standard way of estimating the GPD parameters is the maximum likelihood method, which provides asymptotically normal estimators in the unconditional case with $\gamma>-1/2$ (Drees et al, 2004). To estimate the covariate-dependent variables $\gamma_{0}(\mathbf{x})$ and $\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))$, we also use the maximum likelihood method. Specifically, given observations $(y_{i},\mathbf{x}_{i}),\,i=1,2,\ldots,N$ and the threshold $u$, we define the exceedances of the observations above the threshold as

where $(a)_{+}=\max\{0,a\}$. In other words, $z_{i}=0$ whenever the value $y_{i}$ is below the threshold. Denoe the log-likelihood function as

and the loss function is calculated only on positive $z_{i}$ as

Then, to estimate $\gamma_{0}(\mathbf{x})$ and $\sigma_{0}(\mathbf{x};\tau_{u}(\mathbf{x}))$, we can apply the censored maximum likelihood method via solving the following optimization problem

Inverting the distribution function in (2.5), we have an approximation of the quantile for probability level $\tau_{N}>\tau_{0}$ by taking $u=Q_{Y}(\tau_{0}|\mathbf{x})$ and $z=Q_{Y}(\tau_{N}|\mathbf{x})-Q_{Y}(\tau_{0}|\mathbf{x})$. As a result, we have that

where for simplicity, we define

We denote the solution of the problem (2.8) as $\widehat{\gamma}(\mathbf{x})$ and $\widehat{\sigma}(\mathbf{x};\tau_{0})$. Therefore, using plug-in estimators, the prediction of the extreme conditional quantile $Q_{Y}(\tau_{N}|\mathbf{x})$ is

It is important to note that the estimations in (2.3) and (2.8) are usually implemented individually, such as the stages two and three of the algorithm in [12]. However, the combination of (2.3) and (2.8) is essentially a bilevel programming problem when the threshold $u$ in (2.8) equals the conditional quantile $Q(\tau|\mathbf{x})$ determined by (2.2) and (2.3). In the following, we apply algorithms of bilevel programming to obtain the joint estimation in (2.3) and (2.8).

In this section, we formulate the extreme quantile estimation via GPD approximation as a bilevel programming problem. Bilevel programming, a class of optimization problems characterized by their hierarchical structure, has found wide-ranging applications across diverse domains including economics, engineering, transportation, and environmental management (Sinha et al, 2017). These problems feature a unique interplay between two optimization tasks: an upper-level decision-maker seeks to optimize a given objective while factoring in the response of a lower-level decision-maker, who, in turn, optimizes their objective in response to the upper-level decisions. To illustrate this concept in our context, we combine (2.3) and (2.8) and obtain

where $\mathbf{x}$ and the prefixed level $\tau_{0}$ is given. The conditional quantile regression and the maximum likelihood function are denoted as the lower- and upper-level functions of (2.12), respectively.

The hierarchical structure of bilevel programming introduces unique challenges that pose difficulties for existing solution methods. For instance, penalty function methods address bilevel optimization by solving a sequence of unconstrained optimization problems. These unconstrained problems are derived by introducing a penalty term that quantifies the extent of constraint violation. The penalty term typically relies on a parameter, taking a zero value for feasible solutions and a positive value (indicating minimization) for infeasible ones. However, many penalty function methods require the upper-level function to be convex, which conflicts with the non-convex nature of the log-likelihood function (2.7). This discrepancy arises from the fact that the Hessian matrix $\frac{\partial l(\gamma,\sigma|z)}{\partial\gamma\partial\sigma}$ is not semi-positive definite. Consequently, various penalty methods, such as those proposed by Aiyoshi and Shimizu (1981) and Lv et al (2007), tend to be ineffective in addressing the problem. Additionally, nested evolutionary algorithms have gained popularity for dealing with bilevel problems (Sinha et al, 2014). These algorithms involve solving a lower-level optimization problem for each upper-level member, which are typically employed in two main ways. The first approach combines an evolutionary algorithm at the upper level with classical algorithms at the lower level (Mathieu et al, 1994). In contrast, the second approach leverages evolutionary algorithms at both the upper and lower levels (Angelo et al, 2013). While these nested strategies can be effective, they face a significant computational burden, rendering them impractical for large-scale bilevel problems.

To handle problem (2.12) effectively, we first reformulate the lower-level function as

The above problem can be transformed into a linear program as

where $\lambda_{i}^{+},\lambda_{i}^{-},i=1,2,\ldots,N$ are Lagrange multiplicators. With (2.13), the bilevel programming (2.12) can be rewritten as

where we let $z_{i}=(y_{i}-u)_{+}$ and $u=\alpha+\boldsymbol{\beta}^{T}\mathbf{x}$. Then, we replace the constraints in problem (2.14) with its Karush-Kuhn-Tucker (KKT) conditions. These conditions manifest as Lagrangian and complementarity constraints, simplifying the bilevel optimization problem into a single-level constrained optimization problem. Define the augmented Lagrangian (AL) function as

where $t_{i}^{+},t_{i}^{-},i=1,2,\ldots,N,$ are the Lagrange multiplicators. The derivatives of $L_{N}^{(\text{AL})}$ with respect to the parameters $\alpha,\boldsymbol{\beta},\lambda^{+}_{i},\lambda^{-}_{i}$ are given as follows:

Therefore, the problem (2.14) is reformulated as

where $\Theta_{\mathbf{x},\tau_{0}}=(\gamma,\sigma,\alpha,\boldsymbol{\beta},\lambda^{+}_{i},\lambda^{-}_{i},t^{+}_{i},t^{-}_{i},u_{i})$ denotes the variable set for a prefixed $\mathbf{x}$ and level $\tau_{0}$. Problem (2.16) is a special case of mathematical programs with equilibrium constraints (MPEC), in which the constraints include the complementarity conditions

MPEC plays an important role in various fields, such as the Stackelberg game in economic sciences, and problem (2.16) could be efficiently solved via interior point algorithm or active set algorithm. For more details on these algorithms, we refer readers to the work of [16]. In the following, we leave out the solutions of the Lagrange multiplicators. and denote the solutions of problem (2.16) as $\Theta_{\mathbf{x},\tau_{0}}=(\widetilde{\alpha}(\tau_{0}),\widetilde{\boldsymbol{\beta}}(\tau_{0}),\widetilde{\gamma}(\mathbf{x}),\widetilde{\sigma}(\mathbf{x};\tau_{0}))$ for simplicity.

Note that the solution to (2.16) is derived given a specific covariate $\mathbf{x}$. This leads to estimations $\widetilde{\gamma}(\mathbf{x})$ and $\widetilde{\sigma}(\mathbf{x};\tau_{0})$ that are inherently covariate-dependent. To illustrate this, recall that $\widetilde{\gamma}(\mathbf{x})$ and $\widetilde{\sigma}(\mathbf{x};\tau_{0})$ are utilized in modeling the exceedances over a threshold. We estimate this threshold through (2.4), based on a prefixed $\tau_{0}$. Consequently, the estimated threshold becomes intricately linked to the covariates, $\mathbf{x}$. This, in turn, results in varying exceedances for different values of $\mathbf{x}$, ultimately leading to distinct GPD models and covariate-dependent estimations, namely $\widetilde{\gamma}(\mathbf{x})$ and $\widetilde{\sigma}(\mathbf{x};\tau_{0})$.

In practice, we often encounter situations where we need to deal with an infinite number of future observations, in addition to the finite historical observations at hand. A typical example is the realm of online streaming data, continually generated from various sources (Gaber et al, 2005). Online streaming data finds applications across diverse industries, including finance, transportation, and e-commerce. Notably, streaming data volumes can be vast, occasionally approaching infinity. When it comes to the task of conditional quantile estimation with these infinite observations from online streaming data, applying (2.16) entails solving the problem an infinite number of times. This introduces a significant computational complexity that poses a substantial challenge.

To alleviate the computational burden associated with solving (2.16) for possible infinite observations, we employ B-spline interpolation. To clarify, let us denote the finite offline observations as $(y_{i}^{(\text{off})},\mathbf{x}_{i}^{(\text{off})})$ for $i=1,\ldots,N^{(\text{off})}$ and the future online observations as $\mathbf{x}_{i}^{(\text{on})}$ for $i=1,\ldots,N^{(\text{on})}$ where we allow $N^{(\text{on})}$ to tend to infinite. It is important to note that we assume the covariates for the two observation sets are distinct. B-spline interpolation empowers us to solve (2.16) only $N^{(\text{off})}$ times, obtaining estimations for any number of data points, even as $N^{(\text{on})}$ approaches infinity. In more specific terms, assume that we have obtained $\Theta_{\mathbf{x}_{i}^{(\text{off})},\tau_{0}}=(\widetilde{\alpha}(\tau_{0}),\widetilde{\boldsymbol{\beta}}(\tau_{0}),\widetilde{\gamma}(\mathbf{x}_{i}^{(\text{off})}),\widetilde{\sigma}(\mathbf{x}_{i}^{(\text{off})};\tau_{0}))$ for $i=1,\ldots,N^{(\text{off})}$, for given offline observations by solving (2.16). Notably, the estimations $\widetilde{\alpha}(\tau_{0})$ and $\widetilde{\boldsymbol{\beta}}(\tau_{0})$ are solely linked to the choice of $\tau_{0}$, which is assumed to be the same for both offline and online observations in our paper. As a result, for $\mathbf{x}_{i}^{(\text{on})}$, the estimations of $\alpha$ and $\boldsymbol{\beta}$ remain consistent with those for $\mathbf{x}_{i}^{(\text{off})}$. Therefore, we only need to apply B-spline interpolation to estimate the extreme value index $\gamma(\mathbf{x})$ and the scale $\sigma(\mathbf{x},\tau_{0})$ for $\mathbf{x}=\mathbf{x}_{i}^{(\text{on})}$.

B-splines, which stand for “Basis splines”, hold a pivotal role in the realm of numerical analysis and computer-aided design. They provide a powerful method for representing and approximating complex curves and surfaces, as well as for solving various mathematical and engineering problems. B-splines are piecewise polynomials, and the positions where the pieces meet are known as knots, which are defined as a non-decreasing sequence of real numbers

where $M$ is the knot number. Provided that $M\geq d+2$, we can define the B-splines of degree $d\geq 0$ over the knots $\xi$ using a recursive formulation. Specifically, the $j$-th B-spline $B_{j,d,\boldsymbol{\xi}}(x):\mathbb{R}\rightarrow\mathbb{R}$ can be calculated by

where $\xi_{j}\leq\xi_{j+1}\leq\cdots\leq\xi_{j+d+1}$ are $d+2$ real numbers taken from the knot sequence $\xi$. When $\xi_{j}=\xi_{j+d+1}$, $B_{j,d,\boldsymbol{\xi}}(x)$ is identically zero. The foundation of B-splines begins with the initial B-spline of order zero, represented as:

In simpler terms, a B-spline operates over $M$ knots, strictly adhering to their non-decreasing order. B-splines contribute meaningfully only within the range defined by the first and last knots, remaining zero elsewhere. It is worth noting that while two B-splines can share some of their knots, two B-splines defined over exactly the same set of knots are identical. In essence, a B-spline is uniquely characterized by its knot sequence.

B-splines provide a powerful method for approximating complex functions. In scenarios where the covariants are multi-dimensional, we estimate the covariate-dependent variables $\gamma(\mathbf{x})$ and the scale $\sigma(\mathbf{x},\tau_{0})$ as

Here, $\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)},\ldots,\boldsymbol{\xi}^{(p)}$ represent the $p$ sets of knots, with each set containing $M$ knots. It is assumed that $\mathbf{x}(i)$ falls within the range defined by the boundary knots $\xi^{(i)}_{1}$ and $\xi^{(i)}_{M}$. The choice of $D$, representing the maximum degree of approximation, is crucial. While higher degrees tend to yield more accurate estimates, they also come with increased computational complexity. Here, we opt for $D=3$, a well-known choice commonly referred to as cubic spline interpolation. To determine the coefficients $t_{\gamma,j},j=1,\ldots,J$, for $\mathbf{s}=[s_{1},s_{2},\ldots,s_{N}]^{\top}\in\mathbb{R}^{N}$, we introduce the B-spline matrix $\mathbf{B}_{d,\boldsymbol{\xi}}^{\mathbf{s}}\in\mathbb{R}^{N\times J}$ as

We further define the feature vector $\mathbf{v}_{j}$ as follows

where $j=1,2,\ldots,p$ and $\mathbf{x}_{i}^{{}^{(\text{off})}}(j)$ represents the $j$-th entry of observation $\mathbf{x}_{i}^{{}^{(\text{off})}}$. For pairs $(\widetilde{\gamma}(\mathbf{x}_{i}^{(\text{off})}),\mathbf{x}_{i}^{{}^{(\text{off})}}),i=1,2,\ldots,N^{(\text{off})}$, we we can determine the coefficients $t_{\gamma,j},j=1,\ldots,J$ by solving the following matrix equation

with the boundary knots satisfy

for $i=1,\ldots,p$. The solution to this equation is denoted as $\widetilde{t}_{\gamma,j},j=1,\ldots,J$. Then for new observation $\mathbf{x}^{{(\text{on})}}$, we can estimate the corresponding extreme value index by

when $\boldsymbol{\xi}^{(i)}_{1}\leq\mathbf{x}^{(\text{on})}(i)\leq\boldsymbol{\xi}^{(i)}_{M}$. For cases where $\mathbf{x}(i)$ falls beyond the range of boundary knots $\xi^{(i)}_{1}$ and $\xi^{(i)}_{M}$, the estimation of the extreme value index can be expressed as follows

where $\widetilde{\mathbf{x}}_{\gamma}(i)$ is determined by minimizing the absolute difference as

In summary, we estimate $\gamma(\mathbf{x})$ as

Similarly, the scale parameter $\sigma(\mathbf{x},\tau_{0})$ can be estimated using B-spline interpolation as:

In this expression, the coefficients $\widetilde{t}_{\sigma,j}$ for $j=1,2,\ldots,J$ are determined by solving the matrix equation given by:

and $\widehat{\mathbf{x}}_{\sigma}(i)$ satisfying

In this way, the estimation of extreme conditional quantile for observation $\mathbf{x}^{(\text{on})}$ is calculated as

We summarize our method in Algorithm 1.

It is important to emphasize that B-spline interpolation, while highly effective for approximating scale and shape parameters, does introduce a degree of bias into the estimation process. This trade-off is essential as it significantly reduces computational time, making it a valuable tool for online observations. For this specific context, there is no need to engage in solving the censored maximum likelihood problem required for GPD distribution fitting when dealing with online data. Instead, we capitalize on historical offline data to perform the approximation using B-spline interpolation. It is noteworthy that the solutions derived from equations like (2.17) and (2.20) exclusively rely on the historical dataset. Consequently, we can obtain these solutions well in advance of the arrival of new online data. This proactive approach implies that once the online data is observed, we can directly substitute the results into (2.21) to efficiently obtain the parameter estimations, thus drastically reducing the computational burden associated with online data processing.