Quantile regression in Python: qreg, bootstrap and simultaneous quantiles

OLS describes how the average outcome changes with the regressors. Often the question is about other parts of the distribution: does schooling raise low wages as much as high ones, or does a policy reduce the risk of very low birth weight? Quantile regression estimates how a chosen quantile of the outcome, such as the median or the 90th percentile, shifts with the regressors. Fitting several quantiles shows whether a variable moves the whole distribution or changes its spread.

What quantile regression estimates

The model states that the τ-th conditional quantile of y is linear in x: Q_τ(y | x) = x′β(τ). Koenker and Bassett (1978) showed that β(τ) can be estimated by minimizing an asymmetrically weighted sum of absolute residuals:

β̂(τ) = argmin_b Σᵢ ρ_τ(yᵢ − xᵢ′b),   ρ_τ(u) = u · (τ − 1[u < 0])

At τ = 0.5 the loss is the sum of absolute deviations, and the estimator is median regression. Residuals above the fitted line get weight τ and those below get weight 1 − τ, so at τ = 0.9 the fitted line sits near the top of the conditional distribution.

The slope β(τ) says how the τ-th quantile of the conditional distribution moves when x changes. If all slopes are equal, x shifts the distribution without changing its shape. If the slope rises with τ, x also widens the spread. Koenker and Hallock (2001) give an accessible introduction with economic examples. Two points of interpretation matter:

• β(0.9) describes the 90th percentile of the outcome among people with the same x. It is not the effect for a particular person at the 90th percentile, unless people would keep their rank in the distribution.
• The median is less affected by extreme outcomes than the mean, but quantile regression is not protected against outliers in the regressors.

The estimator does not need an assumption on the error distribution, but its standard errors do. The classical (iid) covariance assumes that the error density at the quantile is the same for every x. The robust sandwich allows it to vary with x, using the fitted-quantile density estimate of Hendricks and Koenker (1992). The cluster-robust version of Parente and Santos Silva (2016) allows dependence within clusters. A bootstrap avoids estimating the density altogether.

In OpenEconometrics

oe.qreg fits one quantile with analytic standard errors. The call oe.qreg(data=df, y="lwage", x=["educ"], quantile=0.5, covariance="robust") corresponds to the Stata command qreg lwage educ, quantile(.5) vce(robust).

• covariance="nonrobust" (the default, Stata's vce(iid)), "robust", or "cluster" with cluster= for one cluster column.
• density= chooses how the error density is estimated: "fitted", "residual" (iid only) or "kernel", which is the default for the cluster covariance. bandwidth= is "hsheather" (the default), "bofinger" or "chamberlain", and kernel= selects the kernel.
• weights= with weight_type="aweight", "fweight" or "pweight". Sampling weights select the robust covariance.

The linear program is solved by an interior-point method and then certified as an exact vertex solution. When the minimizer is not unique, as often happens with discrete outcomes, one optimal solution is reported and a warning is recorded.

oe.bsqreg fits one quantile with pairs-bootstrap standard errors (reps=, seed=). oe.sqreg fits several quantiles on the same bootstrap samples, so the result holds the covariance between quantiles. Its terms are named by equation, such as q10:educ and q90:educ, and oe.test and oe.lincom accept those names for cross-quantile tests. The model in the example corresponds to sqreg lwage educ, quantiles(.1 .5 .9) reps(200) followed by test [q90]educ = [q10]educ. With cluster=, both bootstrap functions resample whole clusters. oe.iqreg estimates the difference between two quantile regressions directly. All tests use Student t with N − K degrees of freedom. The quantile family guide gives the density, bandwidth and bootstrap formulas.

Supported options

Options of oe.qreg

Options of oe.bsqreg

Options of oe.sqreg

Options of oe.iqreg

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 2,000 workers with 8 to 20 years of schooling. Log wages follow a location-scale model: lwage = 1.0 + 0.08 · educ + (0.3 + 0.02 · educ) · e, with e standard normal. Schooling raises the mean and also widens the spread, so the true slope of the τ-th quantile is 0.08 + 0.02 · z_τ, where z_τ is the standard normal quantile: 0.0544 at τ = 0.1, 0.08 at the median and 0.1056 at τ = 0.9.

from statistics import NormalDist

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=2000, seed=5):
    """Log wages whose spread widens with schooling (a location-scale model)."""
    rng = np.random.default_rng(seed)
    educ = rng.integers(8, 21, size=n)
    lwage = 1.0 + 0.08 * educ + (0.3 + 0.02 * educ) * rng.normal(size=n)
    return pd.DataFrame({"lwage": lwage, "educ": educ})


def true_slope(tau):
    """Slope of the tau-th conditional quantile: 0.08 + 0.02 * z_tau."""
    return 0.08 + 0.02 * NormalDist().inv_cdf(tau)


if __name__ == "__main__":
    df = simulate()
    median = oe.qreg(data=df, y="lwage", x=["educ"], quantile=0.5, covariance="robust")
    print(median.summary())
    ols = oe.ols(data=df, y="lwage", x=["educ"], covariance="HC1")
    print(f"OLS slope (conditional mean): {ols.coefficients[1].estimate:.4f}")
    print()
    sq = oe.sqreg(data=df, y="lwage", x=["educ"], quantiles=[0.1, 0.5, 0.9],
                  reps=200, seed=1)
    coef = {c.term: c for c in sq.coefficients}
    for tau, eq in [(0.1, "q10"), (0.5, "q50"), (0.9, "q90")]:
        c = coef[f"{eq}:educ"]
        print(f"tau = {tau}: slope {c.estimate:.4f} (bootstrap se {c.std_error:.4f}), "
              f"true {true_slope(tau):.4f}")
    gap = oe.lincom(sq, {"q90:educ": 1, "q10:educ": -1})
    print(f"q90 - q10 slope: {gap['estimate']:.4f}, 95% CI [{gap['ci_low']:.4f}, "
          f"{gap['ci_high']:.4f}], true {true_slope(0.9) - true_slope(0.1):.4f}")
    equal = oe.test(sq, {"q90:educ": 1, "q10:educ": -1})
    print(f"Test of equal slopes: F({equal['df_num']}, {equal['df_denom']:.0f}) = "
          f"{equal['statistic']:.2f}, p = {equal['p_value']:.3g}")

With openecon 0.3.20a1 from PyPI the script prints:

Median regression — lwage
Observations: 2000  |  Covariance: robust  |  Confidence: 95%

Term        Estimate  Std. error        t     P>|stat|   CI lower   CI upper
---------  ---------  ----------  -------  -----------  ---------  ---------
Intercept    1.13792   0.0505376  22.5162  3.0693e-100     1.0388    1.23703
educ       0.0716132  0.00377604  18.9652  6.94202e-74  0.0642079  0.0790186

quantile: 0.5  |  pseudo_r_squared: 0.116624  |  sum_adev: 458.49  |  sum_rdev: 519.02
raw_quantile: 2.1006  |  sparsity: 1.28467  |  density: 0.77841  |  bandwidth: 0.0771127
iterations: 10  |  df_model: 1  |  df_resid: 1998
OLS slope (conditional mean): 0.0803

tau = 0.1: slope 0.0633 (bootstrap se 0.0062), true 0.0544
tau = 0.5: slope 0.0716 (bootstrap se 0.0038), true 0.0800
tau = 0.9: slope 0.1038 (bootstrap se 0.0059), true 0.1056
q90 - q10 slope: 0.0406, 95% CI [0.0251, 0.0560], true 0.0513
Test of equal slopes: F(1, 1998) = 26.53, p = 2.86e-07

Reading the output

• The median slope misses the truth in this sample. Median regression estimates 0.0716 with a 95% interval from 0.0642 to 0.0790, which excludes the true 0.08: the estimate lies 2.2 standard errors below it. A 95% interval misses in about one sample out of twenty, and this is such a sample. The errors are normal and symmetric here, so the conditional mean and median coincide, and OLS, which is more precise with normal errors, lands at 0.0803.
• The slopes rise across the distribution. The bootstrap estimates are 0.0633 at τ = 0.1, 0.0716 at the median and 0.1038 at τ = 0.9, against true values of 0.0544, 0.08 and 0.1056. Schooling raises the 90th percentile of log wages much more than the 10th.
• The difference is estimated with its joint covariance. Because oe.sqreg fits all quantiles on the same 200 resamples, oe.lincom can use the covariance between the q90 and q10 slopes. The gap is 0.0406 with a 95% interval from 0.0251 to 0.0560, which contains the true 0.0513.
• Equal slopes are rejected. The Wald test gives F(1, 1998) = 26.53 with p = 2.86e-07, so schooling changes the spread of wages, not only their level. One OLS slope would hide this.
• The header reports the check-loss fit. The pseudo R² of 0.117 compares the minimized sum of check losses (458.49) with the sum about the raw median (519.02). The sparsity of 1.285 is the estimated inverse density of the errors at the median, which the iid covariance uses; the robust covariance shown here estimates a density for each observation instead.

Limitations

• Cluster-robust covariance for oe.qreg allows at most 1 clustering dimension.
• oe.bsqreg accepts no weights.
• oe.sqreg accepts no weights.
• oe.iqreg accepts no weights.
• The analytic covariances of oe.qreg need τ ± h inside (0, 1), where h is the bandwidth. Extreme quantiles in small samples raise bandwidth_out_of_range; use the bootstrap functions instead.
• density="residual" gives a single sparsity and is accepted only with the nonrobust covariance.
• The Hall–Sheather and Chamberlain bandwidths fix their normal quantile at the 5% level, so analytic standard errors do not change with alpha.
• The cluster covariance of oe.qreg applies no finite-sample correction. The default of 20 bootstrap replications in oe.bsqreg and oe.sqreg gives noisy standard errors; use a few hundred for reported results.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Hendricks, W. and Koenker, R. (1992). Hierarchical spline models for conditional quantiles and the demand for electricity. Journal of the American Statistical Association, 87(417), 58–68. doi:10.1080/01621459.1992.10475175
• Koenker, R. and Bassett, G. (1978). Regression quantiles. Econometrica, 46(1), 33–50. doi:10.2307/1913643
• Koenker, R. and Hallock, K. F. (2001). Quantile regression. Journal of Economic Perspectives, 15(4), 143–156. doi:10.1257/jep.15.4.143
• Parente, P. M. D. C. and Santos Silva, J. M. C. (2016). Quantile regression with clustered data. Journal of Econometric Methods, 5(1), 1–15. doi:10.1515/jem-2014-0011