Asymptotic Linearization of Quantile Regression

Why Linearize Quantile Regression?

Quantile regression estimates how a chosen conditional quantile changes with covariates. Its estimator minimizes a piecewise-linear check-loss objective, so it does not have the same smooth first-order condition as least squares. To understand its large-sample behavior, we study how the objective changes when the coefficient vector moves a small distance from its population value. Under suitable conditions, that local change is approximately quadratic. The resulting approximation expresses estimation error as an average of observable score contributions.

This linear representation is useful because it separates the estimator's first-order behavior into two parts: a score determined by whether each outcome lies above or below the target conditional quantile, and a matrix describing the local curvature of the objective. The conditional density at the quantile matters: a high density makes the quantile more sharply identified, while a low density makes small changes in probability correspond to larger changes in the outcome scale.

Setup and Regularity Conditions

For independent observations $(Y_i,X_i)$, let $Y_i$ be a scalar outcome and let $X_i$ be a $d$-dimensional regressor vector, usually including an intercept. Fix a quantile level $\tau$ strictly between 0 and 1. Suppose the conditional $\tau$-quantile is linear:

Q_\tau(Y_i\mid X_i)=X_i^\top\beta_\tau.

Write $U_i=Y_i-X_i^\top\beta_\tau$ for the true-model residual. The check loss is $\rho_\tau(u)=u\{\tau-\mathbf{1}(u<0)\}$, and the sample estimator minimizes the average check loss over coefficient vectors. The event $U_i=0$ is exactly an outcome at the fitted population quantile.

Definition: The conditional-density-weighted design matrix for fixed-quantile regression is $$ J_\tau=E\!\left[f_{Y\mid X}(X^\top\beta_\tau\mid X)XX^\top\right], $$ where $f_{Y\mid X}(\cdot\mid X)$ is the conditional density of $Y$ given $X$, evaluated at the target conditional quantile. The score contribution is $X_i\psi_\tau(U_i)$, where $\psi_\tau(u)=\tau-\mathbf{1}(u<0)$.

A standard asymptotic linearization requires conditions that make the target locally identified and the objective's local curvature well behaved. In particular, assume the observations are independent and identically distributed; the conditional quantile is correctly specified by the linear model; and there is no conditional atom at the target, so $P(U_i<0\mid X_i)=\tau$. Assume the conditional density exists and is continuous in a neighborhood of the target, with enough local regularity to approximate probabilities by density times interval length. Assume the regressor moments needed for the expansion are finite, and that $J_\tau$ is finite and nonsingular. Also assume the estimator is consistent and the usual uniform control of the local objective holds. These conditions are standard sufficient conditions, not a claim that every possible set of assumptions must be identical.

The no-atom and positive-curvature conditions matter especially after considering quantile identification with atoms and ties. At an atom, the probability can jump at the target rather than change smoothly with a small movement in the coefficient. Then a density-based quadratic expansion need not apply. Nonsingularity of $J_\tau$ also rules out directions in which the regressors provide no local information about the coefficient vector.

Local Expansion of the Check-Loss Objective

Define the sample objective as $Q_n(\beta)=\sum_{i=1}^n\rho_\tau(Y_i-X_i^\top\beta)$. To examine a local perturbation, move from $\beta_\tau$ to $\beta_\tau+h/\sqrt n$, where $h$ is a fixed $d$-vector. The scale $1/\sqrt n$ is the relevant scale for a regular estimator. The check-loss identity below separates a linear score term from a remainder:

\rho_\tau(u-v)-\rho_\tau(u) =-v\psi_\tau(u) +\int_0^v\{\mathbf{1}(u\leq s)-\mathbf{1}(u\leq 0)\}\,ds \quad (u\ne 0).

Here $v=X_i^\top h/\sqrt n$. The first term contributes the empirical score. For the integral term, conditional density regularity implies that, for small $v$, its conditional expectation given $X_i$ is approximately $\tfrac12 f_{Y\mid X}(X_i^\top\beta_\tau\mid X_i)v^2$. Summing across observations therefore produces the quadratic curvature involving $J_\tau$.

Center at the population coefficient.
Set $U_i=Y_i-X_i^\top\beta_\tau$, so the local residual after perturbation is $U_i-X_i^\top h/\sqrt n$.

Separate score and curvature.
Apply the check-loss identity. The score terms add to $-h^\top n^{-1/2}\sum_i X_i\psi_\tau(U_i)$; the expected remainder adds to $\tfrac12h^\top J_\tau h$.

Control the approximation.
Under the stated regularity conditions, the centered remainder is negligible uniformly over bounded $h$, giving the local quadratic approximation.

Minimize the quadratic.
The unique minimizer of the limiting quadratic is $J_\tau^{-1}n^{-1/2}\sum_iX_i\psi_\tau(U_i)$, which yields the linear representation.

More precisely, the local objective difference has the expansion

Q_n\!\left(\beta_\tau+\frac{h}{\sqrt n}\right)-Q_n(\beta_\tau) = -h^\top\frac{1}{\sqrt n}\sum_{i=1}^nX_i\psi_\tau(U_i) +\frac12h^\top J_\tau h+o_p(1).

The sign of the score term follows from the fact that increasing a coefficient changes the residual from $U_i$ to $U_i-X_i^\top h/\sqrt n$. Minimizing the displayed quadratic gives the following main result.

Theorem: Under the stated fixed-dimensional regularity conditions, a consistent quantile regression estimator satisfies $$ \sqrt n(\hat\beta_\tau-\beta_\tau) = J_\tau^{-1}\frac{1}{\sqrt n}\sum_{i=1}^nX_i\psi_\tau(U_i)+o_p(1). $$ Equivalently, its first-order influence contribution for observation $i$ is $J_\tau^{-1}X_i\psi_\tau(U_i)$. The equality is an asymptotic approximation: it does not assert that the estimator is exactly a linear function of the data in a finite sample.

Worked Examples

Worked Example: An Intercept-Only Median

Suppose a hypothetical sample has $n=8$ outcomes $7, 9, 9.6, 10.2, 10.8, 11.1, 12,$ and $14$, measured in outcome units. Consider the intercept-only median model, with population median $\beta_{0.5}=10$. Assume the population density at 10 is $0.20$ per outcome unit and is smooth and positive near 10. These assumptions give $J_{0.5}=f_Y(10)=0.20$, which is nonsingular for this one-parameter model.

Three observations are below 10 and five are above it, with no observation equal to 10. Thus $\psi_{0.5}(U_i)=0.5-1=-0.5$ for each of the three lower observations and $0.5$ for each of the five higher observations. The score sum is

\sum_{i=1}^8\psi_{0.5}(U_i) =3(-0.5)+5(0.5)=1.

The first-order approximation to the estimator's deviation from the population median is therefore

\hat\beta_{0.5}-\beta_{0.5} \ \approx\ \frac{1}{n}J_{0.5}^{-1}\sum_{i=1}^n\psi_{0.5}(U_i) = \frac{1}{8}\left(\frac{1}{0.20}\right)(1) =0.625

in outcome units. This is a linearization, not an exact finite-sample calculation. In fact, the conventional sample median for these eight values is $(10.2+10.8)/2=10.5$, a deviation of $0.5$ from 10. The difference between $0.625$ and $0.5$ is not a contradiction: the approximation describes first-order large-sample behavior, and this sample is small.

Worked Example: A Binary Regressor and Matrix Curvature

Consider a hypothetical median regression with $X=(1,Z)^\top$, where $Z$ is binary and half the population has $Z=0$ and half has $Z=1$. Suppose the conditional density at the median is $0.20$ for $Z=0$ and $0.10$ for $Z=1$, in inverse outcome units. Both groups have a continuous, positive density near their conditional median. Then

J_{0.5} =E[f_{Y\mid X}(X^\top\beta_{0.5}\mid X)XX^\top] = \begin{pmatrix} 0.15&0.05\\ 0.05&0.05 \end{pmatrix}.

For example, the upper-left entry is $(0.5)(0.20)+(0.5)(0.10)=0.15$, and the lower-right entry is $(0.5)(0)+(0.5)(0.10)=0.05$. Its determinant is $0.15(0.05)-(0.05)^2=0.005$, so this matrix is nonsingular.

Now suppose a hypothetical sample has 50 observations in each group. Among the $Z=0$ observations, 24 residuals relative to the true conditional median are below zero; among the $Z=1$ observations, 26 are below zero. There are no residual ties at zero. At the median, the group score sums are $25-24=1$ for $Z=0$ and $25-26=-1$ for $Z=1$. Hence the full score vector is

\sum_{i=1}^{100}X_i\psi_{0.5}(U_i) = \begin{pmatrix}1+(-1)\\-1\end{pmatrix} = \begin{pmatrix}0\\-1\end{pmatrix}.

The inverse curvature matrix is

J_{0.5}^{-1} = \begin{pmatrix} 10&-10\\ -10&30 \end{pmatrix}.

Thus the coefficient deviation is approximated by

\hat\beta_{0.5}-\beta_{0.5} \ \approx\ \frac{1}{100} \begin{pmatrix} 10&-10\\ -10&30 \end{pmatrix} \begin{pmatrix}0\\-1\end{pmatrix} = \begin{pmatrix}0.10\\-0.30\end{pmatrix}.

The approximation is an intercept increase of $0.10$ outcome units and a $Z=1$ versus $Z=0$ median difference decrease of $0.30$ outcome units. The sign and scale depend on both the score imbalance and the density-weighted design, not on the score counts alone.

Worked Example: Checking Local Curvature in a Direction

For a hypothetical regression with $X=(1,Z)^\top$, suppose $P(Z=1)=0.25$, and the conditional error density at the target quantile is $0.2542$ per outcome unit in both groups. The density is assumed continuous and positive near the target, and the two distinct regressor values give a full-rank design. Since $E(XX^\top)=\begin{pmatrix}1&0.25\\0.25&0.25\end{pmatrix}$, the curvature matrix is

J_\tau = 0.2542 \begin{pmatrix} 1&0.25\\ 0.25&0.25 \end{pmatrix} = \begin{pmatrix} 0.2542&0.06355\\ 0.06355&0.06355 \end{pmatrix}.

Take the coefficient direction $h=(0.2,-0.1)^\top$. The quadratic curvature in this direction is

h^\top J_\tau h =(0.2)^2(0.2542) +2(0.2)(-0.1)(0.06355) +(-0.1)^2(0.06355)

=0.010168-0.002542+0.0006355 =0.0082615 \approx0.008262.

This value is positive, as expected for a nonzero direction under a positive-definite curvature matrix. In the local objective expansion, the quadratic contribution for this $h$ is one-half of the value, approximately $0.004131$. Reporting the curvature and the quadratic contribution separately avoids a common factor-of-two error.

What the Representation Says

The score contribution $X_i\psi_\tau(U_i)$ is positive or negative according to whether the observation lies above or below the target quantile, with magnitude scaled by $\tau$ or $1-\tau$. Conditional on $X_i$, its expected value is zero because $P(U_i<0\mid X_i)=\tau$. The matrix $J_\tau$ converts the average score imbalance into coefficient units and accounts for the amount of conditional density near the target.

Under the stated moment conditions, the central limit theorem applied to the average score yields an asymptotically normal coefficient estimator. The linear representation is the key step: it identifies the random contributions whose large-sample variation drives the estimator. The exact covariance calculation and its estimation are separate tasks; here the emphasis is on deriving the representation and understanding the conditions that make it valid.

Common Mistakes and Limitations

Using a smooth least-squares derivative. Check loss is not everywhere differentiable. The correct local expansion uses the quantile score and the conditional density at the target, rather than a least-squares residual derivative.
Forgetting the density in $J_\tau$. The curvature is not generally just $E(XX^\top)$. Conditional density weights observations according to how much probability lies near the target quantile.
Reversing the score sign. With residual $U_i=Y_i-X_i^\top\beta_\tau$, the score used here is $X_i\{\tau-\mathbf{1}(U_i<0)\}$. The representation has the positive factor $J_\tau^{-1}$ multiplying its sum.
Assuming ties are harmless. An atom at the target can make the conditional probability jump instead of varying smoothly. In that case, the density-based quadratic expansion may fail, and the usual representation should not be applied without different conditions and analysis.
Treating the approximation as an exact finite-sample identity. The $o_p(1)$ remainder vanishes at the relevant asymptotic scale; it need not be small in a particular small sample.

Key takeaway: With a correctly specified conditional quantile, a smooth positive conditional density at that quantile, and nonsingular weighted design, the estimator's scaled error is asymptotically the inverse local curvature times the summed quantile scores.

Check Your Understanding

Use the definitions and expansion above to answer the following questions.

Why does the conditional density at the target quantile enter $J_\tau$, rather than only the regressor second moments?
In the intercept-only case, if $n=20$, $\tau=0.25$, and 3 residuals are below zero with no ties, what is the sum of the quantile scores?
What does nonsingularity of $J_\tau$ rule out for the local coefficient approximation?
Why can an atom at the target quantile invalidate the density-based quadratic expansion?
In the local objective expansion, why is the quadratic contribution one-half of $h^\top J_\tau h$?