The Sandwich Variance for Quantile Regression

Why Quantile Regression Needs a Sandwich Variance

A fitted conditional quantile describes a location in the outcome distribution, but its uncertainty depends on more than the spread of the observations. In linear quantile regression, the slope and intercept are informed by observations near the target conditional quantile. How densely observations are packed around that quantile affects how quickly the fitted coefficients respond to sampling variation. The regressors and the quantile score also determine how much variation enters the estimating equation.

The resulting covariance has a sandwich form: a matrix describing local sensitivity surrounds a matrix describing score variation. This structure matters especially when the conditional density at the target quantile changes with the regressors. In that case, a single common density cannot generally represent the uncertainty of every coefficient.

Model, Score, and Regularity Conditions

Definition: For a quantile level $0<\tau<1$, suppose the conditional quantile of $Y_i$ given the regressor vector $X_i$ is linear: $$ Q_{Y_i}(\tau\mid X_i)=X_i'\beta_\tau. $$ The residual at this quantile is $U_i=Y_i-X_i'\beta_\tau$, and the quantile score is $\psi_\tau(U_i)=\tau-\mathbf{1}\{U_i<0\}$. The conditional-density bread is $$ J_\tau=E\!\left[f_{Y\mid X}(X_i'\beta_\tau\mid X_i)X_iX_i'\right], $$ and the quantile-score meat is $$ \Omega_\tau=E\!\left[\psi_\tau(U_i)^2X_iX_i'\right]. $$

Here $f_{Y\mid X}(y\mid X_i)$ is the conditional density of the outcome, evaluated at the fitted conditional quantile. The bread describes the local change in the expected score when the coefficient vector changes. The meat is the variance contribution from the score. Both are matrices with the same dimensions as $X_iX_i'$.

The standard result below assumes independent, identically distributed observations; a correctly specified linear conditional quantile; no atom at the conditional quantile; a conditional density that is positive, finite, and sufficiently smooth near that quantile; and enough finite moments of the regressors for the relevant expectations to exist. The matrix $J_\tau$ must be nonsingular, which requires that the regressors identify all coefficient directions. These conditions allow the score to be locally approximated by a linear function of the coefficient error and its sample variation to satisfy a central limit theorem.

Correct specification gives $P(U_i<0\mid X_i)=\tau$. Consequently, conditional on $X_i$, the score is centered and its squared value has expectation $\tau(1-\tau)$. The meat therefore simplifies to

\Omega_\tau=\tau(1-\tau)E[X_iX_i'].

This simplification applies under the stated conditional-quantile model. The bread does not usually simplify in the same way: if the conditional density differs across regressor values, its value remains inside the expectation defining $J_\tau$.

The Asymptotic Covariance

Theorem: Under the regularity conditions above, the quantile-regression estimator $\hat{\beta}_\tau$ satisfies $$ \sqrt{n}(\hat{\beta}_\tau-\beta_\tau)\ \xrightarrow{d}\ N\!\left(0,\ J_\tau^{-1}\Omega_\tau J_\tau^{-1}\right). $$ Thus the asymptotic covariance matrix of $\hat{\beta}_\tau$ itself is $n^{-1}J_\tau^{-1}\Omega_\tau J_\tau^{-1}$.

The covariance of the scaled estimator, $\sqrt{n}(\hat{\beta}_\tau-\beta_\tau)$, is the matrix inside the normal distribution. For the unscaled coefficient estimator, divide that matrix by $n$. The two inverse bread matrices appear on opposite sides of the meat because score variation is translated into coefficient variation through the sensitivity of the estimating equation.

Why the Sandwich Appears

The quantile-regression score equation is the sample condition $\sum_{i=1}^n X_i\psi_\tau(Y_i-X_i'\hat{\beta}_\tau)=0$, with the usual qualification that an optimizer may satisfy the nonsmooth first-order condition through a subgradient. The asymptotic linearization of this equation yields

\sqrt{n}(\hat{\beta}_\tau-\beta_\tau) = J_\tau^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^n X_i\psi_\tau(U_i)+o_p(1).

The key steps in the covariance calculation are as follows.

Identify score variation.
Each observation contributes $X_i\psi_\tau(U_i)$. Under correct specification its expectation is zero, and its covariance is $\Omega_\tau$.

Apply the central limit theorem.
For independent observations with finite second moments, $n^{-1/2}\sum_i X_i\psi_\tau(U_i)$ converges to a mean-zero normal vector with covariance $\Omega_\tau$.

Translate score variation into coefficient variation.
The linearization premultiplies that score sum by $J_\tau^{-1}$. Therefore its covariance is $J_\tau^{-1}\Omega_\tau(J_\tau^{-1})'$. Since $J_\tau$ is symmetric, this is $J_\tau^{-1}\Omega_\tau J_\tau^{-1}$.

This argument also explains why the density belongs in the bread. Near the target quantile, a small change in a fitted value changes the probability of an observation falling below that value at a rate determined by the conditional density. The derivative of that probability supplies the factor $f_{Y\mid X}(X_i'\beta_\tau\mid X_i)$ in $J_\tau$.

Worked Calculations

Worked Example: A Binary Regressor With Different Conditional Densities

Suppose $X_i=(1,D_i)'$, where $D_i$ is zero or one with probability one-half each. At the conditional median, assume the conditional densities are $0.25$ outcome-units$^{-1}$ when $D_i=0$ and $0.10$ outcome-units$^{-1}$ when $D_i=1$. Assume the conditional median is linear in $D_i$, observations are independent, and there is no atom at either median. Calculate the asymptotic covariance of the coefficient estimator.

At $\tau=0.5$, the conditional-quantile model is correctly specified and the score meat is $0.25E[X_iX_i']$. The regressor moment is

E[X_iX_i'] =\frac{1}{2} \begin{pmatrix}1&0\\0&0\end{pmatrix} +\frac{1}{2} \begin{pmatrix}1&1\\1&1\end{pmatrix} = \begin{pmatrix}1&0.5\\0.5&0.5\end{pmatrix}.

The assumptions on the conditional densities make the bread

J_{0.5} =\frac{1}{2}(0.25) \begin{pmatrix}1&0\\0&0\end{pmatrix} +\frac{1}{2}(0.10) \begin{pmatrix}1&1\\1&1\end{pmatrix} = \begin{pmatrix}0.175&0.050\\0.050&0.050\end{pmatrix}.

The meat is

\Omega_{0.5} =0.25 \begin{pmatrix}1&0.5\\0.5&0.5\end{pmatrix} = \begin{pmatrix}0.250&0.125\\0.125&0.125\end{pmatrix}.

The determinant of $J_{0.5}$ is $0.175(0.050)-0.050^2=0.00625$, so

J_{0.5}^{-1} =\frac{1}{0.00625} \begin{pmatrix}0.050&-0.050\\-0.050&0.175\end{pmatrix} = \begin{pmatrix}8&-8\\-8&28\end{pmatrix}.

Multiplying out the sandwich gives

J_{0.5}^{-1}\Omega_{0.5}J_{0.5}^{-1} = \begin{pmatrix}8&-8\\-8&58\end{pmatrix}.

This is the asymptotic covariance of the coefficients after scaling by $\sqrt{n}$; the asymptotic covariance of the unscaled estimator is this matrix divided by $n$. The intercept and slope covariance is negative in this parameterization, while the slope variance is especially large. The result reflects the relatively low density at the conditional median in the $D_i=1$ group: a lower density means that small changes in the fitted quantile correspond to larger changes in the score.

Worked Example: An Intercept-Only Quantile

Consider estimating the $\tau=0.25$ quantile from independent observations with no regressors other than an intercept. Assume the population quantile has density $0.20$ outcome-units$^{-1}$, and the density is positive and smooth near that quantile. Find the asymptotic variance for a sample of $n=100$ and its standard error.

With an intercept, $X_i=1$. The bread is $J_{0.25}=0.20$, and the correctly specified score meat is $\Omega_{0.25}=0.25(1-0.25)=0.1875$. Thus the asymptotic variance of the scaled estimator is

J_{0.25}^{-1}\Omega_{0.25}J_{0.25}^{-1} =\frac{0.1875}{0.20^2} =4.6875.

For $n=100$, the asymptotic variance of the sample quantile estimator is $4.6875/100=0.046875$ outcome-units squared. Its standard error is therefore

SE=\sqrt{0.046875}=0.2165\text{ outcome units}.

The score variation alone is $0.1875$, but the density rescales it by the inverse density on both sides. Omitting that rescaling would give the wrong variance and the wrong units.

Worked Example: A Confidence Interval for the Binary-Regressor Slope

Use the covariance calculation from the first example for a hypothetical sample of $n=200$, with the same design proportions and conditional densities. Suppose the estimated conditional-median slope is $0.80$ outcome units per unit of $D_i$. Calculate an approximate 95% confidence interval, treating the bread and meat values used above as known for this illustration.

The estimated asymptotic variance of the slope is the slope diagonal element, $58$, divided by $n$. Hence its standard error is

SE(\hat{\beta}_{0.5,1})=\sqrt{\frac{58}{200}} =\sqrt{0.29}=0.5385\text{ outcome units}.

Using the standard normal critical value $1.96$, the interval is

0.80\pm1.96(0.5385) =0.80\pm1.0555 =[-0.2555,\ 1.8555]\text{ outcome units}.

Under the model and large-sample approximation, this procedure has 95% coverage across repeated samples. The interval includes zero, so it does not rule out a zero conditional-median difference at the 5% two-sided level. This is not a statement that there is a 95% probability that the fixed population slope lies in this particular interval.

Estimating and Interpreting the Sandwich

In practice, the matrices are replaced by sample-based estimates. A common plug-in form is

\widehat{\operatorname{Var}}(\hat{\beta}_\tau) = \frac{1}{n}\hat{J}_\tau^{-1}\hat{\Omega}_\tau\hat{J}_\tau^{-1}.

Here $\hat{J}_\tau$ estimates the expectation of the conditional-density-weighted regressor products, and $\hat{\Omega}_\tau$ estimates the variation in the quantile scores. Under the correctly specified model, a population expression for the meat is $\tau(1-\tau)E[X_iX_i']$; its sample counterpart uses the observed regressor products. A residual-score version instead forms the sample average of $X_iX_i'\psi_\tau(\hat{U}_i)^2$. These are large-sample approaches to the score variation, not a reason to ignore the bread.

The bread and meat play distinct roles. The meat records score variation. The bread describes how responsive the score is to coefficient changes. When conditional densities vary with $X_i$, the bread weights different regressor patterns differently. In general, replacing it with one common density can produce incorrect standard errors even if the score meat is calculated correctly.

A further check is dimensional: a density has inverse-outcome units, so $J_\tau^{-1}$ contributes outcome units per regressor scale. The complete covariance must then have the corresponding squared coefficient units. This is a useful way to detect a formula that has accidentally omitted one of the bread factors.

Common Mistake

Common mistake: Using only $\tau(1-\tau)E[X_iX_i']$ as the coefficient covariance, or inserting a heteroskedasticity-robust score meat while leaving out the density bread. The expression $\tau(1-\tau)E[X_iX_i']$ is the meat under correct specification; it is not the full coefficient covariance. The complete asymptotic covariance is $J_\tau^{-1}\Omega_\tau J_\tau^{-1}$. Both inverse bread factors are required because score fluctuations must be translated into coefficient fluctuations.

The sandwich is useful because it separates two sources of uncertainty without treating them as interchangeable. For quantile regression, a robust variance calculation must account for the score and for the conditional density at the target quantile. A plausible-looking score variance alone does not supply valid coefficient uncertainty.

Check Your Understanding

Answer each question using the definitions and result above.

In words, what does the conditional-density bread measure, and what does the quantile-score meat measure?
Under a correctly specified conditional quantile model with no atom at the quantile, simplify $E[\psi_\tau(U_i)^2X_iX_i']$.
Why is the asymptotic covariance of $\hat{\beta}_\tau$ divided by $n$, while the covariance in the limit distribution of $\sqrt{n}(\hat{\beta}_\tau-\beta_\tau)$ is not?
If the conditional density at the target quantile is lower for one regressor group, what does that generally imply about the coefficient uncertainty associated with that group?
What is wrong with treating $\tau(1-\tau)E[X_iX_i']$ as the complete asymptotic covariance matrix?