Heteroskedasticity-Robust Standard Errors in Median Regressions

Question

Does anyone know how to compute heteroskedasticity-robust standard errors in median regressions in R?

Assume the following example:

require("quantreg")

df <- iris
rq(Sepal.Length ~ Sepal.Width + Petal.Length, data = df)

I know that summary.rq offers multiple standard error options. But I do not know which option is the right one in my case. I would like to use the standard errors that Stata employs:

sysuse auto.dta
qreg price mpg trunk, vce(robust)

When running linear regressions in R, I correct standard errors with coeftest:

# Load packages
packs <- list("quantreg", "lmtest", "sandwich")
lapply(packs, require, character.only = T)
    
df <- iris
rq_outp <- lm(Sepal.Length ~ Sepal.Width + Petal.Length, data = df)
coeftest(rq_outp, vcov = vcovHC(rq_outp, type = "HC1"))

However, coeftest does apparently not work on rq objects.

Answer 1

There is no assumption of homoscedasticity in the quantile regression model. You need to assume independence, and you need assume linearity (or correct functional specification) of the quantile (median in your case) function, but you do not assume any particular distribution, and you do not assume homoscedasticity.

I can't say what stata might be doing here, but it does not seem right to me for this application. Whatever standard error methods are used (e.g., bootstrapping rows) should account for heteroscedasticity automatically.