Quantile Regression

A courier company wants to promise a delivery time. A forecast of the average duration does not tell it what fraction of parcels will arrive late. Quoting a much longer time may reduce late deliveries but make the service less attractive. A conditional 90th percentile addresses the operational question: what duration has at least 90 percent of deliveries at or below it for a route with these predictors? An estimated percentile still needs validation before it becomes a service promise.

Ordinary least squares answers a different question: it fits a linear conditional-mean model when that mean is correctly specified. Fitting the coefficients does not require constant conditional variance. A constant-width prediction band is an additional modeling choice, not a restriction imposed by least squares. Quantile regression instead fits a conditional quantile at a chosen level $\tau$. For a continuous conditional distribution with a unique quantile, the true 0.9 quantile has 90 percent of probability at or below it. Fitted linear quantiles need not achieve that exact conditional fraction in a finite sample.

The mechanism is a change of loss function. Squared error $\rho(u) = u^2$ is minimized by the conditional mean. The pinball loss (an asymmetric absolute-value loss) is minimized by the conditional quantile. Once we make that swap, the familiar linear model machinery (design matrices, gradients, regularization) mostly carries over, but the geometry changes in important ways. The objective becomes piecewise linear and non-differentiable at zero, so we solve a linear program rather than a linear system. The fitted curves are also no longer guaranteed to be ordered across $\tau$.

This chapter is the last stop in our tour of Part V: Linear and Generalized Linear Models. Everything we built for multiple linear regression, the design matrix, the intercept column, the loss-plus-optimizer framing from gradient descent and numerical optimization, reappears here with a different loss. What we gain is a model that speaks directly about risk, spread, and tail behavior, which is exactly the language we will need when we move to time-to-event data in Part VI.

The first advantage is a direct description of different parts of the conditional distribution without choosing a complete distributional family. A mean alone does not determine spread or shape, even for symmetric, light-tailed distributions. Fitting levels such as $\tau = 0.05, 0.25, 0.5, 0.75, 0.95$ estimates several conditional cutoffs; the gaps between them describe spread. The linear functional form remains an assumption at each level. Heteroscedasticity can therefore become a feature we estimate rather than something hidden by a single-width summary.

The second advantage is robustness in the response direction. Because the pinball loss grows linearly rather than quadratically, its response score is bounded at fixed predictors, unlike the least-squares score. This does not bound every finite-sample coefficient change; a sample tail quantile can still be determined by one extreme response. At the median, quantile regression is least absolute deviations regression, a maximum likelihood estimate for the regression coefficients under independent Laplace errors centered at zero with a common fixed positive scale, which we first met in maximum likelihood estimation. At a fixed predictor vector, the estimating score is bounded as the response varies; extreme predictors can still have high leverage. For the pinball loss the residual slopes are $\tau$ on the positive side and $\tau - 1$ on the negative side, with magnitudes $\tau$ and $1 - \tau$ respectively, so no slope ever exceeds $1$ in magnitude.

A useful population property is equivariance under continuous strictly increasing transformations. For such a function $h$, quantiles obey $Q_{h(Y)}(\tau)=h(Q_Y(\tau))$; log and exponential transformations on their respective domains satisfy this condition. If the conditional quantile of $\log Y$ is modeled correctly, exponentiating its estimate targets the corresponding quantile of positive $Y$, without a conditional-mean retransformation correction. This does not say that two restricted linear fits, one on $Y$ and one on $\log Y$, agree exactly after transformation in a finite sample.

Finally, quantile regression is a natural tool for asymmetric decision problems. Inventory levels, staffing, safety margins, and service-level agreements can involve quantile targets when their costs or probability constraints have the appropriate form. A blood-pressure reference chart that flags the top 5 percent of children by age is a quantile curve, not a mean curve. When underage and overage costs are positive and linear in the size of the error, expected pinball loss at the cost-implied $\tau$ matches that decision problem up to a positive scale factor.

The main practical annoyance is quantile crossing. In the population the quantile function is nondecreasing in $\tau$, so $Q(0.4 \mid x) \le Q(0.6 \mid x)$ for every $x$. However, each $\tau$ is estimated by its own optimization, so nothing in the procedure forces the fitted curve at $\tau = 0.4$ to lie below the fitted curve at $\tau = 0.6$. In finite samples the curves can cross; the risk depends on sample size, the design, the predictor domain, and functional-form misspecification. At a fixed predictor value, narrowing a nested pair of target levels cannot increase its population quantile gap, leaving less separation for estimation error to overcome. A crossing fitted value is not just undesirable. It cannot be the ordered quantile function of a conditional distribution and can reverse the endpoints of a prediction band. We will see how to diagnose this, and how rearrangement repairs it, later in the chapter.

The second limitation is that efficiency at the center is worse than least squares when least squares is appropriate. For $n$ independent, identically distributed Gaussian observations with positive variance $\sigma^2$, the sample mean has variance $\sigma^2/n$ while the sample median has asymptotic variance $\pi\sigma^2/(2n)$, so their relative efficiency (the ratio of variances) is $(\sigma^2/n)/(\pi\sigma^2/(2n)) = 2/\pi \approx 0.64$ for the median, meaning we need roughly $1/(2/\pi) = \pi/2 \approx 1.57$ times as much data for the same precision. The same $2/\pi$ asymptotic relative-efficiency factor applies to correctly specified median regression with independent homoscedastic Gaussian errors and a regular design. Tail precision depends on the conditional density and design as well as $\tau$: normal tails have low density and can be difficult to estimate, but uncertainty need not increase monotonically toward the tails for every distribution.

Third, the objective is not smooth, so the fitted coefficients need not change smoothly as the data change; jumps can occur when an optimum switches between non-unique solutions. The general multivariable problem has no least-squares-style matrix formula, although special cases have quantile shortcuts. Ordinary differentiation is insufficient at zero residuals; we use a linear program. Its runtime depends on the sample size, predictor count, matrix representation, conditioning, and solver settings, and the piecewise-linear structure can produce non-unique optima when the data are degenerate (duplicate covariate patterns, for example), which complicates inference. Common regular asymptotic standard-error estimates use resampling or estimates of the conditional density; they do not automatically follow the OLS formula $\sigma^2 (X^\top X)^{-1}$.

Fourth, and most often confused, a fitted quantile band is not a confidence interval. The band $[\hat Q(0.05 \mid x), \hat Q(0.95 \mid x)]$ targets coverage of individual responses $Y$ (here $\hat Q(\tau \mid x)$ is the estimated conditional $\tau$-quantile), a statement about the spread of the data; a confidence interval is a statement about the uncertainty of an estimated parameter such as a regression coefficient. The first can be wide because the data are noisy; the second can be narrow because the sample is large. Both need to be reported, and they answer different questions.

A regression model for the mean answers the question "what is the average response at this $x$?" A regression model for a quantile answers a sharper question: "what cutoff has at most a fraction $\tau$ of responses strictly below it and at least a fraction $\tau$ at or below it at this $x$?" To make that question precise we need three ingredients: a definition of a conditional quantile, a loss function whose minimizer is that quantile, and a finite-sample optimization problem we can solve. The first two are population statements; the third is what turns the idea into an algorithm.

We begin with notation. Let $Y$ be a response variable and let the full predictor vector be $X = (X_0,X_1,\dots,X_p)$, with $X_0=1$ for the intercept. There are $p$ nonconstant predictors and $p+1$ entries in the design vector. The cumulative distribution function of $Y$ is $F_Y(y) = P(Y \le y)$, and the $\tau$-quantile of $Y$, for $\tau \in (0,1)$, is

$$Q_Y(\tau) = \inf\{y : F_Y(y) \ge \tau\}.$$

The infimum handles distributions with flat regions or jumps; where $F_Y$ is continuous and strictly increasing it agrees with the ordinary inverse; otherwise it remains a generalized inverse. Conditioning on the predictors gives the target:

$$Q_{Y \mid X=x}(\tau) = \inf\{y : F_{Y\mid X=x}(y) \ge \tau\},$$

where:

• $F_{Y\mid X=x}$: the conditional cumulative distribution function of $Y$ given $X = x$
• $F_{Y\mid X=x}(y)$: the value of that conditional CDF evaluated at a candidate response value $y$
• $\inf$: the infimum, the greatest lower bound of the set of response values satisfying the inequality
• $Q_{Y \mid X=x}(\tau)$: the conditional $\tau$-quantile of $Y$ at a given predictor value $x$, the smallest response value whose conditional CDF reaches at least $\tau$

For intuition, consider the case where $Y \mid X = x$ is normal with mean $\mu(x)$ and standard deviation $\sigma(x)$. Then

$$Q_{Y \mid X = x}(\tau) = \mu(x) + \sigma(x)\,\Phi^{-1}(\tau),$$

where:

• $\Phi^{-1}(\tau)$: the quantile function of the standard normal distribution, evaluated at probability $\tau$, which is the same $z$-value used for confidence intervals in the z-test
• $\mu(x)$: the conditional mean of $Y$ given $X = x$
• $\sigma(x)$: the conditional standard deviation of $Y$ given $X = x$

At $\tau = 0.5$ we get $\mu(x)$, because $\Phi^{-1}(0.5) = 0$. At $\tau = 0.95$ we get $\mu(x) + 1.645\sigma(x)$. So if we assume normality, the conditional quantile function is just the mean function shifted by a multiple of the conditional standard deviation. Ordinary least squares estimates $\mu(x)$ and, when it assumes constant variance, treats $\sigma$ as a single number. Quantile regression estimates $Q_{Y\mid X=x}(\tau)$ directly, which lets $\mu$ and $\sigma$ both vary with $x$, and lets the shape vary too.

A quantile is a local summary of a distribution: it depends on where probability mass accumulates near the cut, not just on a global average. That local character explains two features that will reappear throughout the chapter. First, different quantiles can respond to the predictors in different ways, so a single slope is not enough. Second, precision depends on the density near the selected cutoff and on the predictor design. Low normal-tail density makes extreme levels difficult in our Gaussian example, but there is no universal ordering of quantile variances. The pinball loss below encodes exactly this local, asymmetric view of the data.

To estimate a quantile we need a loss function minimized at that quantile. In the following expectation calculations assume $\mathbb E[Y^2]<\infty$ for squared error and $\mathbb E|Y|<\infty$ for pinball loss. For the density-based derivatives, assume the distribution has a continuous density near the candidate cutoff. Start with the familiar case. Minimizing the expected squared error over the candidate value $q$,

$$L_{\text{sq}}(q) = \mathbb{E}\left[(Y - q)^2\right],$$

leads to the mean. Differentiating with respect to $q$ (the finite second moment allows the locally integrable derivative to pass through the expectation) gives

$$\frac{d}{dq} L_{\text{sq}}(q) = -2\,\mathbb{E}[Y - q] = -2\left(\mathbb{E}[Y] - q\right),$$

which is zero exactly when $q = \mathbb{E}[Y]$. Setting the derivative to zero and checking that the second derivative, $2 > 0$, confirms a minimum. Squared error is married to the mean.

Now replace squared error with the asymmetric absolute-value loss known as the pinball loss, or check loss. Here $u = y - q$ is the signed residual between an observation $y$ and a candidate value $q$, and $\tau \in (0,1)$ is the target quantile level:

$$\rho_\tau(u) = \begin{cases} \tau\, u, & u \ge 0,\\[2pt] (\tau - 1)\, u, & u < 0. \end{cases}$$

where:

• $u$: the residual, the signed difference between the observation and the candidate value, $u = y - q$
• $\tau$: the target quantile level in $(0,1)$
• $\rho_\tau(u)$: the pinball loss at level $\tau$ for residual $u$
• $q$: the candidate value at which the pinball loss is being evaluated

For $\tau > 0.5$, the positive-residual branch is steeper: its slope magnitude is $\tau$, compared with $1-\tau$ on the negative side. A compact way to write the same function is

$$\rho_\tau(u) = u\left(\tau - \mathbf{1}\{u < 0\}\right),$$

where:

• $u$: the residual, the signed difference between the observation and the candidate value, $u = y - q$
• $\tau$: the target quantile level in $(0,1)$
• $\mathbf{1}\{\cdot\}$: the indicator function, equal to 1 when its argument is true and 0 otherwise
• $q$: the candidate value at which the pinball loss is being evaluated, so that the loss is a function of $u = y - q$ and of $\tau$
• $\rho_\tau(u)$: the pinball loss at level $\tau$ for residual $u$

Check the two branches: if $u \ge 0$ the indicator is zero and we get $\tau u$; if $u < 0$ the indicator is one and we get $(\tau - 1)u$, which is positive because $\tau - 1 < 0$ and $u < 0$. At $\tau = 0.5$ both slopes have magnitude $0.5$, and $\rho_{0.5}(u) = |u|/2$, so the pinball loss is just half the absolute error.

The asymmetry is the whole point. For a high quantile, say $\tau = 0.9$, a positive residual (an observation above the fitted value) is penalized by slope $\tau = 0.9$, while a negative residual is penalized by slope $1 - \tau = 0.1$. For one unrestricted sample cutoff that minimizes this loss, at most 90 percent of observations lie strictly below it and at least 90 percent lie at or below it. For a low quantile the push goes the other way. This is how a single loss function encodes a target probability level. The left branch has slope $\tau - 1$ (magnitude $1 - \tau$) and the right branch has slope $\tau$. These magnitudes are equal only when $\tau = 0.5$. Any other $\tau$ makes the two arms differ in slope, and the loss is asymmetric.

To see that this loss is minimized by the $\tau$-quantile, let $f(y)$ be the probability density function of $Y$ and write the expected loss at a candidate value $q$, where $q$ is the candidate value being scored, $f(y)$ is the density of $Y$, and $\rho_\tau$ is the pinball loss defined above. We obtain

$$L(q) = \mathbb{E}\left[\rho_\tau(Y - q)\right] = \int_{-\infty}^{q} (\tau - 1)(y - q) f(y)\, dy + \int_{q}^{\infty} \tau (y - q) f(y)\, dy.$$

The first integral collects the observations below $q$, where the residual $y - q$ is negative and the loss uses the $(\tau-1)$ branch; the second collects observations above $q$.

where:

• $q$: the candidate value being scored as a quantile
• $f(y)$: the probability density function of $Y$
• $\rho_\tau$: the pinball loss at level $\tau$, defined above
• $L(q)$: the expected pinball loss at the candidate value $q$

Now differentiate with respect to $q$. The finite first moment makes the expected loss finite; the residual-loss derivatives are bounded by 1 in magnitude, which justifies differentiation at continuity points of the CDF. Both integrands vanish at the boundary $y = q$ because each factor contains $(y-q)$, so the Leibniz rule gives no boundary contribution and we only need the derivative of the integrands.

$$\frac{dL}{dq} = \int_{-\infty}^{q} (1 - \tau) f(y)\, dy \; - \; \int_{q}^{\infty} \tau f(y)\, dy.$$

The derivative of $(\tau-1)(y-q)$ with respect to $q$ is $1-\tau$, and that of $\tau(y-q)$ is $-\tau$. At either moving integration endpoint $y=q$, the factor $y-q$ is zero, so each boundary contribution vanishes. Recognizing the remaining integrals as probabilities gives:

where:

• $F(q)$: the cumulative distribution function of $Y$ evaluated at $q$

$$\begin{aligned} \frac{dL}{dq} &= (1 - \tau) F(q) - \tau\left(1 - F(q)\right) \\ &= (1 - \tau) F(q) - \tau + \tau F(q) \\ &= F(q) - \tau. \end{aligned}$$

The two $\tau F(q)$ terms combine: $(1-\tau)F(q) + \tau F(q) = F(q)$, leaving $F(q) - \tau$. Setting this to zero gives $F(q)=\tau$. When the CDF is strictly increasing at this cutoff, it selects the unique $Q_Y(\tau)$; a flat region can give multiple minimizing values. The second derivative is $f(q) \ge 0$, so $L$ is convex and this stationary point is a global minimum. The pinball loss is therefore the loss whose minimizer is the quantile, exactly as squared error is the loss whose minimizer is the mean.

For a finite sample, define $\hat F_n(q)=n^{-1}\sum_i\mathbf 1\{y_i\le q\}$ and let $\hat F_n(q^-)$ be the fraction strictly below $q$. The subgradient interval of the average pinball loss, as a function of $q$, is $[\hat F_n(q^-)-\tau,\hat F_n(q)-\tau]$. A minimum occurs exactly when this interval contains zero, giving the counting rule:

$$\#\{y_i < q\} \le n\tau \le \#\{y_i \le q\}.$$

In words: the number of observations strictly below the estimate is at most $n\tau$, and the number at or below it is at least $n\tau$; dividing by $n$, at most a fraction $\tau$ of the data lie strictly below the estimate and at least a fraction $\tau$ lie at or below it. Ties at exactly $q$ are what allow both inequalities to hold when $n\tau$ is not an integer, so the two sides can disagree by the mass pinned at $q$.

The figure below makes the asymmetry visible. Each curve is a pinball loss for a different $\tau$. The two constant residual slopes, together with the kink at zero, give the quantile counting condition derived above.

Now we restrict attention to linear quantile functions. For a chosen $\tau$, we model

$$Q_{Y \mid X = x}(\tau) = x^\top \beta(\tau) = \beta_0(\tau) + \beta_1(\tau) x_1 + \cdots + \beta_p(\tau) x_p.$$

The convention that $x_0 = 1$ absorbs the intercept $\beta_0(\tau)$, so that $x^\top \beta(\tau) = \sum_{j=0}^{p} \beta_j(\tau) x_j$ with a $(p+1)$-dimensional parameter vector.

Here:

• $\beta(\tau)$: the coefficient vector at quantile level $\tau$, of length $p + 1$, indexed from $j = 0$ to $j = p$
• $\beta_0(\tau)$: the intercept at level $\tau$, multiplying the constant $x_0 = 1$
• $\beta_j(\tau)$: the coefficient on predictor $x_j$ at level $\tau$, for $j = 1, \dots, p$
• $x_j$: the $j$-th nonconstant predictor, for $j = 1, \dots, p$
• $x_1, \dots, x_p$: the $p$ observed predictors, to which the constant $x_0 = 1$ is prepended
• $x \in \mathbb{R}^{p+1}$: the full predictor vector with the prepended constant, so that $x^\top \beta(\tau)$ is the fitted $\tau$-quantile

The notation emphasizes that the coefficients depend on $\tau$: each quantile level gets its own line. The estimate is the minimizer of the average pinball loss over the sample,

$$\hat\beta(\tau) = \arg\min_{\beta \in \mathbb{R}^{p+1}} \frac{1}{n} \sum_{i=1}^{n} \rho_\tau\left(y_i - x_i^\top \beta\right).$$

Compare this to the least-squares objective, which is a smooth quadratic in $\beta$ and, for a full-column-rank design, admits the mathematical solution $(X^\top X)^{-1} X^\top y$ that we derived in ordinary least squares. Here the objective is a sum of piecewise linear functions of $\beta$. It is convex (a sum of convex functions is convex) but it has kinks wherever a residual $y_i - x_i^\top\beta$ is exactly zero, and some optima have zero residuals. Flat sets of optimal coefficients can also contain points with no zero residuals. At a kink the derivative does not exist, so "set the gradient to zero and solve" does not work. We need the machinery of convex optimization.

The standard route is to convert the problem into a linear program, which we can then hand to a mature solver. The trick is to split each residual into its positive and negative parts. Define, for every observation, the positive and negative parts of the residual:

$$u_i = \max\left(y_i - x_i^\top \beta,\, 0\right), \qquad v_i = \max\left(x_i^\top \beta - y_i,\, 0\right),$$

where:

• $y_i$: the observed response for the $i$-th observation
• $x_i^\top \beta$: the fitted quantile at the $i$-th predictor vector
• $u_i$: the positive part of the residual, nonzero only when the observation lies above the fit
• $v_i$: the negative part of the residual, nonzero only when the observation lies below the fit

Both are non-negative, and their difference reconstructs the signed residual: $u_i - v_i = y_i - x_i^\top\beta$. In the linear program we drop the max definitions and treat $u_i \ge 0$ and $v_i \ge 0$ as the unknowns; the optimal solution automatically sets at most one of the pair strictly positive, since subtracting $t=\min(u_i,v_i)>0$ from both preserves their difference and lowers the cost by $[\tau+(1-\tau)]t=t$. Now check that the pinball loss can be written as a weighted sum of them. Note that if $y_i - x_i^\top\beta > 0$ then $u_i = y_i - x_i^\top\beta$ and $v_i = 0$, while if $y_i - x_i^\top\beta < 0$ then $u_i = 0$ and $v_i = x_i^\top\beta - y_i$. If the residual $y_i - x_i^\top\beta$ is positive then $u_i$ equals it, $v_i = 0$, and

$$\tau u_i + (1-\tau) v_i = \tau\left(y_i - x_i^\top\beta\right) = \rho_\tau\left(y_i - x_i^\top \beta\right).$$

If the residual is negative then $u_i = 0$, $v_i = x_i^\top\beta - y_i$, and

$$\tau u_i + (1-\tau) v_i = (1-\tau)\left(x_i^\top\beta - y_i\right) = (\tau - 1)\left(y_i - x_i^\top\beta\right) = \rho_\tau\left(y_i - x_i^\top \beta\right),$$

using $(1-\tau) = -(\tau-1)$ in the middle step. Either way the identity holds, and the identity is linear in the unknowns $\beta, u, v$. That is precisely why the problem becomes a linear program: we minimize a linear function subject to linear constraints and non-negativity.

Written out in full, for a fixed $\tau$ we solve

$$\min_{\beta,\, u,\, v} \; \sum_{i=1}^{n} \left[\tau u_i + (1-\tau) v_i\right] \quad \text{subject to} \quad \begin{cases} x_i^\top \beta + u_i - v_i = y_i, & i = 1, \dots, n,\\ u_i \ge 0, \quad v_i \ge 0, & i = 1, \dots, n, \end{cases}$$

where:

• $\beta$: the coefficient vector, unrestricted in sign
• $u_i, v_i$: the positive and negative parts of the $i$-th residual, both constrained to be non-negative
• $\tau$: the target quantile level
• $x_i$: the $i$-th predictor vector, including the intercept entry
• $y_i$: the $i$-th observed response

To pass this to a solver we stack everything into one vector. Let $X$ be the $n \times (p+1)$ design matrix whose $i$-th row is $x_i^\top$, let $I_n$ denote the $n \times n$ identity matrix, let $\mathbf{1}_n$ be a column of $n$ ones, let $0_{p+1}$ be a column of $p+1$ zeros, and let $z$ be the stacked unknown vector defined below. Collect the unknowns as $z = (\beta, u, v) \in \mathbb{R}^{(p+1) + 2n}$. Then the program reads

$$\min_{z} \; c^\top z$$

subject to the equality constraint and the variable-bound restrictions

$$\begin{bmatrix} X & I_n & -I_n \end{bmatrix} z = y, \qquad z_j \ \text{free for } j = 1,\dots,p+1, \quad z_j \ge 0 \ \text{for } j = p+2,\dots,(p+1)+2n,$$

where:

• $z = (\beta, u, v)$: the stacked vector of unknowns, of length $(p+1) + 2n$
• $c = \left(0_{p+1},\ \tau \mathbf{1}_n,\ (1-\tau)\mathbf{1}_n\right)$: the linear cost vector
• $X$: the $n \times (p+1)$ design matrix
• $I_n$: the $n \times n$ identity matrix
• $y$: the observed response vector

The constraint $\begin{bmatrix} X & I_n & -I_n \end{bmatrix} z = y$ is the residual decomposition $x_i^\top\beta + u_i - v_i = y_i$ in matrix form, and $z_j \ge 0$ enforces non-negativity of the split residuals.

Reading the block structure: the first block $X\beta$ produces the fitted values, adding $u$ and subtracting $v$ reconciles them with the observed $y$, and the cost vector puts weight $\tau$ on positive residuals (observations above the fitted line, that is, under-predictions) and weight $1-\tau$ on negative residuals (observations below the fitted line, that is, over-predictions). The single equation per observation is exactly the residual decomposition we derived. The problem has $2n + p + 1$ unknowns and $n$ equality constraints. That is larger than a least-squares solve. The two identity blocks contain only $2n$ nonzero entries; at fixed predictor count the block matrix becomes sparse as $n$ grows, and a sparse representation lets solvers exploit this structure.

The dual shows how the fitted score limits sensitivity to extreme responses at fixed predictors. The dual of the program above is

$$\max_{a \in \mathbb{R}^n} \; y^\top a \quad \text{subject to} \quad X^\top a = 0, \qquad a_i \in [\tau - 1,\ \tau] \ \text{ for all } i,$$

where:

• $a \in \mathbb{R}^n$: the dual variable vector, one multiplier per observation
• $a_i$: the dual multiplier for the $i$-th observation, confined to the interval $[\tau - 1, \tau]$
• $y$: the observed response vector
• $X$: the $n \times (p+1)$ design matrix

Each observation gets a multiplier $a_i$ that is capped between $\tau - 1$ and $\tau$. At the solution, $a_i = \tau$ for observations with positive residual (points above the fitted line) and $a_i = \tau - 1$ for observations with negative residual (points below the fitted line), with $a_i$ taking any value in $[\tau-1,\tau]$ when the residual is exactly zero. The constraint $X^\top a = 0$ is the first-order condition: the fitted line is positioned so that the signed, capped residuals balance against every column of the design matrix. Every multiplier is confined to the interval $[\tau - 1, \tau]$, whose bound on the magnitude is $\max(\tau, 1 - \tau) \le 1$. Because $a_i$ is bounded independently of $y_i$, an extreme response cannot push its own multiplier beyond this cap. The equilibrium condition $X^\top a = 0$ balances the capped scores against the design columns; it does not guarantee a bound on finite-sample coefficient changes.

For the intercept column, $X_0=\mathbf 1_n$, stationarity gives $\sum_i a_i=0$. Let $r_i=y_i-x_i^\top\beta$, $N_-=\#\{i:r_i<0\}$ and $N_0=\#\{i:r_i=0\}$. At nonzero residuals the dual values are fixed at $\tau$ or $\tau-1$; at zero residuals they may take any value in $[\tau-1,\tau]$. The available zero-residual values balance the intercept precisely when

$$N_-\le n\tau\le N_-+N_0.$$

Here $N_-$ counts negative residuals, $N_0$ counts ties at the fitted value, and $n\tau$ is the target count. This is the same counting inequality derived for a sample quantile. It does not require exactly a fraction $\tau$ of residuals to be negative. The full vector condition also needs valid zero-residual multipliers that satisfy $X^\top a=0$. We will read those multipliers from the LP result rather than treating an arbitrary endpoint of each subgradient interval as an ordinary gradient.

These properties are not decorative. They tell us what quantile regression preserves, what it gives up, and where its finite-sample behavior can differ from the population ideal.

• Monotone equivariance. Population conditional quantiles commute with continuous strictly increasing transformations on the response domain. A selected empirical order statistic has the same property. A general linear quantile-regression estimator need not commute with nonlinear transformations because the transformed function may leave the linear model family.
• Median reduction. At $\tau = 0.5$ the pinball loss is half the absolute error, so the estimator is least absolute deviations regression. For independent Laplace errors centered at zero with a common fixed positive scale, its regression coefficients are maximum likelihood estimates.
• No ordering guarantee across $\tau$. Because each $\tau$ is optimized separately, $\hat\beta(\tau_1)$ and $\hat\beta(\tau_2)$ for $\tau_1 < \tau_2$ are not constrained to produce ordered fitted values at any given $x$. Population quantiles are ordered. Appropriate consistency assumptions can recover that order in the limit, but separate finite-sample fits impose no ordering constraint.
• Asymptotic normality. For independent identically distributed rows, a correctly specified linear conditional quantile, bounded predictors with nonsingular second-moment matrix, and a positive continuous conditional density with adequate local smoothness near that quantile, the usual regular quantile-regression limit applies. Define the density-weighted design moment $D$ and the plain design moment $J$ by

$$D = \mathbb{E}\left[f_{Y\mid X}\left(Q_\tau(x) \mid x\right) x x^\top\right], \qquad J = \mathbb{E}\left[x x^\top\right],$$

where:

• $f_{Y\mid X}(\cdot \mid x)$: the conditional density of $Y$ given $X = x$, evaluated at the quantile $Q_\tau(x)$
• $Q_\tau(x) = x^\top\beta(\tau)$: the true conditional $\tau$-quantile function
• $x$: a predictor vector drawn from the population distribution (the expectation is taken over $X$)
• $D$: the density-weighted design moment; density affects precision through the full sandwich covariance, together with the predictor design and $\tau$
• $J$: the unweighted design moment $\mathbb{E}[x x^\top]$

Then $\sqrt{n}\left(\hat\beta(\tau) - \beta(\tau)\right) \xrightarrow{d} \mathcal{N}\left(0,\ \tau(1-\tau)\, D^{-1} J D^{-1}\right)$ as $n \to \infty$, where $\tau(1-\tau)$ is the variance of the pinball-loss score, $D^{-1} J D^{-1}$ is the sandwich matrix (bread times meat times bread), and $D$ and $J$ are defined immediately above. The density term $D$ records the local information near the quantile. At fixed $\tau$ and predictor distribution, multiplying all conditional densities at the quantile by the same factor $c$ divides the asymptotic covariance by $c^2$; more general precision comparisons depend on the whole density-weighted design. In distributions such as the normal example, low tail density can increase this covariance; the formula does not imply that every distribution has larger variance at more extreme levels.

• Under these regular population conditions, the influence function is bounded as $y$ varies at a fixed predictor vector. It need not be bounded as $x$ varies: high-leverage predictors multiply the capped scores. This population statement does not bound every finite-sample coefficient change.

The formulas are compact, but the logic is easiest to trust when we can compute a few quantiles by hand. We start with a single median, then a slope, and finally the heteroscedastic case that motivates the entire method. The numbers are small enough to check on paper.

Take five response values $y = (2, 3, 5, 8, 10)$ and ask for the median, $\tau = 0.5$. The counting rule says we need a value $q$ with at most $0.5 \times 5 = 2.5$ observations strictly below it and at least $2.5$ at or below it. The value $q = 5$ has two observations below it $(2, 3)$ and three at or below it $(2, 3, 5)$, so it satisfies the rule. Let us confirm by evaluating the average pinball loss at three candidates.

• At $q = 3$ the residuals are $u = (-1, 0, 2, 5, 7)$. With $\tau = 0.5$, $\rho_{0.5}(u) = 0.5|u|$, so the losses are $(0.5, 0, 1.0, 2.5, 3.5)$ and their average is $7.5/5 = 1.5$.
• At $q = 5$ the residuals are $u = (-3, -2, 0, 3, 5)$, losses $(1.5, 1, 0, 1.5, 2.5)$, average $6.5/5 = 1.3$.
• At $q = 6$ the residuals are $u = (-4, -3, -1, 2, 4)$, losses $(2, 1.5, 0.5, 1, 2)$, average $7/5 = 1.4$.

These values compare the same average loss at three candidate medians.

Here the pinball loss at $\tau = 0.5$ is $\rho_{0.5}(u) = 0.5|u|$, and the residual is the signed distance from each observation to the candidate value $q$.

where:

• $u$: the residual $y_i - q$, the signed distance from each observation to the candidate value $q$
• $\rho_{0.5}(u)$: the pinball loss at $\tau = 0.5$, which equals $0.5|u|$
• $q$: the candidate value being evaluated as an estimate of the sample median

The minimum is at $q = 5$, which is the sample median. Moving $q$ up from 5 adds loss on the three points below and removes loss on the two above; because there are more points below than above, that trade is unfavorable. The balance point is exactly where the counts cross, which is the counting rule.

Now consider a one-parameter model with no intercept, $Q(\tau \mid x) = \beta x$, fitted to four observations, where $\beta$ is the single slope to be estimated and $x$ is the predictor:

Because all $x_i > 0$ and the residual in observation $i$ is $y_i - \beta x_i = x_i(r_i - \beta)$, the pinball loss is positively homogeneous of degree one: for $c > 0$, $\rho_\tau(cu) = c\,\rho_\tau(u)$. Writing $y_i - \beta x_i = x_i(r_i - \beta)$ with $x_i > 0$, the objective becomes

$$\sum_{i=1}^{4} \rho_\tau\left(y_i - \beta x_i\right) = \sum_{i=1}^{4} x_i\, \rho_\tau\left(r_i - \beta\right).$$

where:

• $r_i = y_i / x_i$: the ratio of response to predictor for the $i$-th observation
• $x_i$: the $i$-th predictor value, which also serves as the weight
• $\beta$: the single slope parameter being estimated
• $\rho_\tau$: the pinball loss at level $\tau$, which is positively homogeneous of degree one

This is a weighted pinball loss in the single unknown $\beta$, and its minimizer is the weighted $\tau$-quantile of the ratios $r_i$ with weights $x_i$. Order the ratios and accumulate weight: $1.50$ with weight 2, then $2.00$ with weight 1 (cumulative 3), then $3.00$ with weight 3 (cumulative 6), then $3.25$ with weight 4 (cumulative 10). For $\tau = 0.5$ the target cumulative weight is $0.5 \times 10 = 5$. The cumulative weight jumps from $3$ (after $2.00$) to $6$ (after $3.00$), so it crosses the target of $5$ at the ratio $3.00$; by the weighted-median convention we take the smallest ratio whose cumulative weight is at least the target, giving $\hat\beta(0.5) = 3$.

Checking the loss: at $\beta = 3$ the fitted values are $(3, 6, 9, 12)$ and the residuals are $(-1, -3, 0, 1)$, giving an average loss of $0.5(1 + 3 + 0 + 1)/4 = 0.625$. At $\beta = 2.75$ the residuals are $(-0.75, -2.5, 0.75, 2)$, average loss $0.5(0.75 + 2.5 + 0.75 + 2)/4 = 0.75$. At $\beta = 3.25$ the residuals are $(-1.25, -3.5, -0.75, 0)$, average loss $0.5(1.25 + 3.5 + 0.75 + 0)/4 = 0.6875$. Both neighbours are worse, confirming the minimum. Note also that this weighted-median minimizer applies because all slopes share a single parameter $\beta$ multiplying positive weights $x_i$; the objective is a weighted pinball loss in $\beta$ whose minimizer is the weighted $\tau$-quantile of the ratios.

This shortcut is worth remembering. For a single slope, no intercept, and positive predictors, the objective reduces to a weighted $\tau$-quantile of the ratios $r_i=y_i/x_i$ with weights $x_i$. At $\tau=0.5$ this is a weighted median. That is why quantile regression sometimes feels like a ranking operation rather than a fitting operation.

Now the case that matters. Take three illustrative quantile lines with the coefficients below. These are chosen for hand arithmetic, not claimed solver output from the later simulated dataset. Over $0\le x\le10$ their ordering gives a widening nominal 10th-to-90th-percentile band:

The fitted curves are $Q(0.1 \mid x) = 1.0 + 0.6x$, $Q(0.5 \mid x) = 2.0 + 1.1x$, and $Q(0.9 \mid x) = 3.6 + 2.0x$, where:

• $Q(\tau \mid x)$: the fitted conditional $\tau$-quantile at predictor value $x$
• $\beta_0(\tau)$: the intercept of the fitted line at level $\tau$
• $\beta_1(\tau)$: the slope of the fitted line at level $\tau$

Evaluate them at three values of $x$ and record the width of the central 80 percent band, which is defined as $W(x) = Q(0.9 \mid x) - Q(0.1 \mid x)$, where $Q(\tau \mid x)$ is the fitted conditional $\tau$-quantile and $W(x)$ is the resulting band width in response units:

• At $x = 1$: the three values are $1.6$, $3.1$, $5.6$, so $W(1) = 5.6 - 1.6 = 4.0$.
• At $x = 5$: the values are $4.0$, $7.5$, $13.6$, so $W(5) = 13.6 - 4.0 = 9.6$.
• At $x = 10$: the values are $7.0$, $13.0$, $23.6$, so $W(10) = 23.6 - 7.0 = 16.6$.

The width is not constant: because both bounding curves are linear, the width is exactly linear in $x$. Subtracting the two fitted lines $\hat Q(0.9 \mid x) = \beta_0(0.9) + \beta_1(0.9)\,x$ and $\hat Q(0.1 \mid x) = \beta_0(0.1) + \beta_1(0.1)\,x$ term by term,

$$\begin{aligned} W(x) &= \left[\beta_0(0.9) - \beta_0(0.1)\right] + \left[\beta_1(0.9) - \beta_1(0.1)\right] x \\ &= (3.60 - 1.00) + (2.00 - 0.60)\,x \\ &= 2.60 + 1.40\,x, \end{aligned}$$

which reproduces $4.0$, $9.6$, and $16.6$ at $x = 1, 5, 10$ (for example, $W(5) = 2.60 + 1.40 \times 5 = 9.60$). The slope difference $1.40$ is the rate at which the conditional spread grows: each additional unit of $x$ widens the 10-to-90 range by 1.40 units. An additional constant-variance normal residual model would produce a constant-width band. Least squares itself does not force that variance model.

We can also price individual predictions with the pinball loss. Suppose a new observation arrives with $x = 5$ and $y = 13.0$. The $\tau = 0.9$ prediction is $13.6$, so the residual is $u = 13.0 - 13.6 = -0.6$. Because $u < 0$ we use the lower branch: $\rho_{0.9}(-0.6) = (0.9 - 1)(-0.6) = (-0.1)(-0.6) = 0.06$. Overpredicting this response at a high quantile is cheap because the lower branch of $\rho_\tau$ has slope $-(1-\tau) = -0.1$ (in magnitude $0.1$), so a negative residual of size $0.6$ is penalized by only $0.1 \times 0.6 = 0.06$; a true 0.9 quantile of a continuous distribution has 90 percent of its probability below it. Our illustrative coefficients carry no empirical calibration guarantee; the small penalty here follows from the chosen loss. Now take $y = 7.0$ at the same $x$ and score it against the $\tau = 0.1$ line, whose prediction is $4.0$. The residual is $3.0 > 0$, so $\rho_{0.1}(3.0) = 0.1 \times 3.0 = 0.30$. The two residuals, $-0.6$ and $+3.0$, are not of comparable size, so this comparison mainly illustrates the pricing rule: the pinball loss charges $\tau$ per unit on the positive side and $1 - \tau$ per unit on the negative side, with both rates evaluated at the relevant $\tau$.