The QuantileRegression Package

by Johannes Ludsteck (Johannes.Ludsteck@iab.de)

Quantile regression is a straightforward generalization of the median regression. Median regression is an alternative to least squares fitting. In a median regression absolute deviations of the data points from the regression hyperplane are minimized (instead of squared deviations as in least squares fitting). Since median regression puts smaller weights on large deviations than least squares regression, it is more robust against outliers (but less precise than least squares if the tails of the distribution of residuals are 'small enough').
This package implements median (which is the 50%-quantile) and quantile regression. The implementation is 'quick and dirty', i.e. no sophisticated input checking is done and the user is responsible for regression diagnostics. (An advantage of quantile regression is, however, that diagnostic checks are less important than with least squares fitting because of its higher robustness.)

Experienced users are encouraged to improve the package, bug reports and comments are welcome.

This loads the QuantileRegressionPackage and a standard package (employed to generate simulated data).

<< QuantileRegression`
<< Statistics`ContinuousDistributions`;

We generate simulated data for a small linear regression model with 5 regressors and 300 observations.
X is the Matrix of regressors, y the dependent variable and beta a vector of coefficients.

Our model does not contain a constant term. Note, however, that a constant term (i.e. a vector of ones) has to included in most regression problems and that the user is responsible for this by construction of the appropriate regressor matrix.

SeedRandom[123];
n = 300; p = 5;
beta = Table[Random[Real, {-2, 2}], {p}]
nd = NormalDistribution[0. , 1. ];
X = 
        Map[Prepend[#, 1.0] &, 
RandomArray[
            nd, {n, p - 1}]];
y = X . beta + RandomArray[nd, n];

Timing[{b, cov} = QuantileRegression[X, y]]

Here we extract the standard errors from the covariance matrix and compute asymptotic t-values from them

stddev = Sqrt[Tr[cov, List]]
tAs = b/stddev

Bootstrapping of the covariance matrix is extremely time-consuming. Note that the option below requires 200 re-estimates of the model.

However, since median and quantile regression techniques are often used to improve inference in presence of outliers and other problems stemming from 'perverse' residuals, it is not sensible to report asymptotic standard errors because they are derived under assumptions which are likely not met in such cases. In these cases we suggest to use bootstrap standard errors instead.

\(Timing[{b, covBS} = 
        QuantileRegression[X, 
          y, 
CovarianceMatrixEstimate -> {"Bootstrap", 200}]]\)

stddevBS = Sqrt[Tr[covBS, List]]
tBS = b/stddevBS

QuantileRegression is quite fast (by making use of interior point methods provided by the Mathematica functions LinearProgramming and DualLinearProgramming) and well suited for large problems. Here we compute an example with 10000 observations and 50 regressors.

SeedRandom[123];
n = 10000; p = 50;
beta = Table[Random[Real, {-2, 2}], {p}]
nd = NormalDistribution[0. , 1. ];
X = 
        Map[Prepend[#, 1.0] &, 
RandomArray[
            nd, {n, p - 1}]];
y = X . beta + RandomArray[nd, n];

Timing[
      b = QuantileRegression[X, y, 
          CovarianceMatrixEstimate -> None]]

The following example shows an important use of quantile regression. It delivers extremely interesting results in models with heteroscedastic errors and allows direct modelling of heteroscedasticity. Int the example below we generate heteroscedastic residuals and show that the slopes of the regression lines for different quantiles differ. The computation of different regression quantiles is now standard in many economic  applications. If, e.g. the dependent variable is income, we can analyse and test whether low income and high income owners respond differently to economic incentives.

n = 500; p = 2;
nd = NormalDistribution[0. , 1. ];
x = 7 RandomArray[nd, n];
X = 
        Map[Prepend[#, 1.0] &, 
RandomArray[
            nd, {n, p - 1}]];
beta = {0.5, 1.5};
resid = 
        RandomArray[nd, 
            n] (X[[All, 2]] + 5) 0.4;
y = X . beta + resid;

Timing[
      b05 = QuantileRegression[X, 
          y, 
RegressionQuantile -> 0.5, 
          CovarianceMatrixEstimate -> None]]
Timing[
      b09 = QuantileRegression[X, 
          y, 
RegressionQuantile -> 0.9, 
          CovarianceMatrixEstimate -> None]]
Timing[
      b03 = QuantileRegression[X, 
          y, 
RegressionQuantile -> 0.3, 
          CovarianceMatrixEstimate -> None]]

The plot shows that the slope of the regression line differs between quantiles

\(Show[
      ListPlot[Transpose[{X[\([All, 2]\)], 
            y}], 
DisplayFunction -> Identity], 

      p5 = Plot[
          b05[\([1]\)] + b05[\([2]\)]\ z, {z, \(-3\), 
            3}, 
DisplayFunction -> Identity], 

      p9 = Plot[
          b09[\([1]\)] + b09[\([2]\)]\ z, {z, \(-3\), 
            3}, 
DisplayFunction -> Identity], 

      p3 = Plot[
          b03[\([1]\)] + b03[\([2]\)]\ z, {z, \(-3\), 
            3}, 
DisplayFunction -> 
            Identity], 
DisplayFunction -> $DisplayFunction]\)