Quantile Regression › Formula

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.