``````
bandwidth.rq <- function(p, n, hs = TRUE, alpha = 0.05)
{
# Bandwidth selection for sparsity estimation two flavors:
#	Hall and Sheather(1988, JRSS(B)) rate = O(n^{-1/3})
#	Bofinger (1975, Aus. J. Stat)  -- rate = O(n^{-1/5})
# Generally speaking, default method, hs=TRUE is preferred.

x0 <- qnorm(p)
f0 <- dnorm(x0)
if(hs)
n^(-1/3) * qnorm(1 - alpha/2)^(2/3) *
((1.5 * f0^2)/(2 * x0^2 + 1))^(1/3)
else n^-0.2 * ((4.5 * f0^4)/(2 * x0^2 + 1)^2)^ 0.2
}

plot.rq.process <- function(x, nrow = 3, ncol = 2, ...)
{
## Function to plot estimated quantile regression  process
tdim <- dim(x\$sol)
p <- tdim[1] - 3
m <- tdim[2]
par(mfrow = c(nrow, ncol))
ylab <- dimnames(x\$sol)[[1]]
for(i in 1:p) {
plot(x\$sol[1,], x\$sol[3 + i,  ], xlab = "tau",
ylab = ylab[3 + i], type = "l")
}
par(oldpar)
}

print.rqs <- function (x, ...)
{
if (!is.null(cl <- x\$call)) {
cat("Call:\n")
dput(cl)
}
coef <- coef(x)
cat("\nCoefficients:\n")
print(coef, ...)
rank <- x\$rank
nobs <- nrow(residuals(x))
p <- nrow(coef)
rdf <- nobs - p
cat("\nDegrees of freedom:", nobs, "total;", rdf, "residual\n")
if (!is.null(attr(x, "na.message")))
cat(attr(x, "na.message"), "\n")
invisible(x)
}

"print.rq" <-
function(x, ...)
{
if(!is.null(cl <- x\$call)) {
cat("Call:\n")
dput(cl)
}
coef <- coef(x)
cat("\nCoefficients:\n")
print(coef, ...)
rank <- x\$rank
nobs <- length(residuals(x))
if(is.matrix(coef))
p <- dim(coef)[1]
else p <- length(coef)
rdf <- nobs - p
cat("\nDegrees of freedom:", nobs, "total;", rdf, "residual\n")
if(!is.null(attr(x, "na.message")))
cat(attr(x, "na.message"), "\n")
invisible(x)
}

print.summary.rqs <- function(x, ...) {
lapply(x, print.summary.rq)
invisible(x)
}

"print.summary.rq" <-
function(x, digits = max(5, .Options\$digits - 2), ...)
{
cat("\nCall: ")
dput(x\$call)
coef <- x\$coef
## df <- x\$df
## rdf <- x\$rdf
tau <- x\$tau
cat("\ntau: ")
print(format(round(tau,digits = digits)), quote = FALSE, ...)
cat("\nCoefficients:\n")
print(format(round(coef, digits = digits)), quote = FALSE, ...)
invisible(x)
}

"rq" <-
function (formula, tau = 0.5, data, subset, weights, na.action, method = "br",
model = TRUE, contrasts = NULL, ...)
{
call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action"), names(mf), 0)
mf <- mf[c(1,m)]
mf\$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval.parent(mf)
if(method == "model.frame")return(mf)
mt <- attr(mf, "terms")
weights <- as.vector(model.weights(mf))
Y <- model.response(mf)
X <- model.matrix(mt, mf, contrasts)
eps <- .Machine\$double.eps^(2/3)
Rho <- function(u,tau) u * (tau - (u < 0))
if(length(tau)>1){
if(any(tau < 0) || any(tau > 1))
stop("invalid tau:  taus should be >= 0 and <= 1")
if(any(tau == 0)) tau[tau == 0] <- eps
if(any(tau == 1)) tau[tau == 1] <- 1 -eps
coef <- matrix(0,ncol(X),length(tau))
rho <- rep(0,length(tau))
fitted <- resid <- matrix(0,nrow(X),length(tau))
for(i in 1:length(tau)){
z <- {if (length(weights))
rq.wfit(X, Y, tau = tau[i], weights, method, ...)
else rq.fit(X, Y, tau = tau[i], method, ...)
}
coef[,i] <- z\$coefficients
resid[,i] <- z\$residuals
rho[i] <- sum(Rho(z\$residuals,tau[i]))
fitted[,i] <- Y - z\$residuals
}
taulabs <- paste("tau=",format(round(tau,3)))
dimnames(coef) <- list(dimnames(X)[[2]],taulabs)
dimnames(resid) <- list(dimnames(X)[[1]],taulabs)
fit <- z
fit\$coefficients <-  coef
fit\$residuals <- resid
fit\$fitted.values <- fitted
if(method == "lasso") class(fit) <- c("lassorqs","rqs")
else class(fit) <- "rqs"
}
else{
process <- (tau < 0 || tau > 1)
if(tau == 0) tau <- eps
if(tau == 1) tau <- 1 -eps
fit <- {
if (length(weights))
rq.wfit(X, Y, tau = tau, weights, method, ...)
else rq.fit(X, Y, tau = tau, method, ...)
}
if(process)
rho <- list(x = fit\$sol[1,],y = fit\$sol[3,])
else {
dimnames(fit\$residuals) <- list(dimnames(X)[[1]],NULL)
rho <-  sum(Rho(fit\$residuals,tau))
}
if(method == "lasso") class(fit) <- c("lassorq","rq")
else class(fit) <- ifelse(process, "rq.process", "rq")
}
fit\$na.action <- attr(mf, "na.action")
fit\$formula <- formula
fit\$terms <- mt
fit\$xlevels <- .getXlevels(mt,mf)
fit\$call <- call
fit\$tau <- tau
fit\$weights <- weights
fit\$residuals <- drop(fit\$residuals)
fit\$rho <- rho
fit\$method <- method
fit\$fitted.values <- drop(fit\$fitted.values)

attr(fit, "na.message") <- attr(m, "na.message")
if(model) fit\$model <- mf
fit
}
"rq.fit" <-
function(x, y, tau = 0.5, method = "br", ...)
{
if(length(tau) > 1) {
tau <- tau[1]
warning("Multiple taus not allowed in rq.fit: solution restricted to first element")
}

fit <- switch(method,
fn = rq.fit.fnb(x, y, tau = tau, ...),
fnb = rq.fit.fnb(x, y, tau = tau, ...),
fnc = rq.fit.fnc(x, y, tau = tau, ...),
pfn = rq.fit.pfn(x, y, tau = tau, ...),
br = rq.fit.br(x, y, tau = tau, ...),
lasso = rq.fit.lasso(x, y, tau = tau, ...),
{
what <- paste("rq.fit.", method, sep = "")
if(exists(what, mode = "function"))
(get(what, mode = "function"))(x, y, ...)
else stop(paste("unimplemented method:", method))
}
)
fit\$fitted.values <- y - fit\$residuals
fit\$contrasts <- attr(x, "contrasts")
fit
}
"rqs.fit"<-
function(x, y, tau = 0.5, tol = 0.0001)
{
# function to compute rq fits for multiple y's
x <- as.matrix(x)
p <- ncol(x)
n <- nrow(x)
m <- ncol(y)
z <- .Fortran("rqs",
as.integer(n),
as.integer(p),
as.integer(m),
as.integer(n + 5),
as.integer(p + 2),
as.double(x),
as.double(y),
as.double(tau),
as.double(tol),
flag = integer(m),
coef = double(p * m),
resid = double(n),
integer(n),
double((n + 5) * (p + 2)),
double(n),
PACKAGE = "quantreg")
if(sum(z\$flag)>0){
if(any(z\$flag)==2)
warning(paste(sum(z\$flag==2),"out of",m,
"BS replications have near singular design"))
if(any(z\$flag)==1)
warning(paste(sum(z\$flag==1),"out of",m,"may be nonunique"))
}
return(t(matrix(z\$coef, p, m)))
}
"formula.rq" <-
function (x, ...)
{
form <- x\$formula
if (!is.null(form)) {
form <- formula(x\$terms)
environment(form) <- environment(x\$formula)
form
}
else formula(x\$terms)
}

"predict.rq" <-
function (object, newdata, type = "none", interval = c("none",
"confidence"), level = 0.95, na.action = na.pass, ...)
{
if (missing(newdata))
return(object\$fitted)
else {
tt <- terms(object)
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action,
xlev = object\$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object\$contrasts)
}
pred <- drop(X %*% object\$coefficients)
dots <- list(...)
if (length(dots\$se))
boot <- (dots\$se == "boot")
else boot <- FALSE
if (length(dots\$mofn))
mofn <- dots\$mofn
interval <- match.arg(interval)
if (!interval == "none") {
if (interval == "confidence") {
if (type == "percentile") {
if (boot) {
if(exists("mofn")) {# Rescale and recenter!!
n <- length(object\$fitted)
factor <- ifelse(mofn < n, sqrt(mofn/n), 1)
XB <- X %*% t(summary(object, cov = TRUE, ...)\$B)/factor
pl <- apply(XB, 1, function(x) quantile(x, (1 - level)/2))
pu <- apply(XB, 1, function(x) quantile(x, 1 - (1 - level)/2))
pl <- pred + factor * (pl - pred)
pu <- pred + factor * (pu - pred)
}
else {
XB <- X %*% t(summary(object, cov = TRUE, ...)\$B)
pl <- apply(XB, 1, function(x) quantile(x, (1 - level)/2))
pu <- apply(XB, 1, function(x) quantile(x, 1 - (1 - level)/2))
}
pred <- cbind(pred, pl, pu)
colnames(pred) <- c("fit", "lower", "higher")
}
else stop("Percentile method requires se = \"boot\".")
}
else if (type == "direct") {
if (boot)
stop("Direct method incompatible with bootstrap covariance matrix estimation")
Z <- rq.fit(object\$x, object\$y, tau = -1)\$sol
V <- summary(object, cov = TRUE, ...)
df <- V\$rdf
tfrac <- qt(1 - (1 - level)/2, df)
Vun <- V\$cov * V\$scale^2
tau <- object\$tau
bn <- tfrac * sqrt(diag(X %*% Vun %*% t(X)))
tauU <- pmin(tau + bn, 1 - 1/df)
tauL <- pmax(tau - bn, 1/df)
tauhat <- Z[1, ]
yhat <- X %*% Z[-(1:3), ]
n <- nrow(X)
pl <- yhat[cbind(1:n, cut(tauL, tauhat, label = FALSE))]
pu <- yhat[cbind(1:n, cut(tauU, tauhat, label = FALSE))]
pred <- cbind(pred, pl, pu)
colnames(pred) <- c("fit", "lower", "higher")
}
else {
V <- summary(object, cov = TRUE, ...)
df <- V\$rdf
tfrac <- qt((1 - level)/2, df)
sdpred <- sqrt(diag(X %*% V\$cov %*% t(X)))
pred <- cbind(pred, pred + tfrac * sdpred %o%
c(1, -1))
colnames(pred) <- c("fit", "lower", "higher")
}
}
else stop(paste("No interval method for", interval))
}
pred
}

"predict.rqs" <-
function (object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)
{
## with all defaults
if(missing(newdata) & !stepfun) return(object\$fitted)

## otherwise
tt <- delete.response(terms(object))
m <- if(missing(newdata)) model.frame(object) else model.frame(tt, newdata,
na.action = na.action, xlev = object\$xlevels)
if(!is.null(cl <- attr(tt, "dataClasses"))) .checkMFClasses(cl, m)
X <- model.matrix(tt, m, contrasts = object\$contrasts)
pred <- X %*% object\$coefficients

## return stepfun or matrix
if(stepfun) {
pred <- t(pred)
taus <- object\$tau
if(type == "Qhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) f[[i]] <- stepfun(taus, c(pred[1,i], pred[,i]))
} else
f <- stepfun(taus, c(pred[1,1], pred[,1]))
}
else if(type == "Fhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) {
o <- order(pred[,i])
f[[i]] <- stepfun(pred[o,i],c(taus[1],taus[o]))
}
} else
f <- stepfun(pred[,1],c(taus[1],taus))
}
else stop("type must be either 'Qhat' or 'Fhat'")
return(f)
} else {
return(pred)
}
}

"predict.rq.process" <-
function (object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)
{
if (missing(newdata))
return(object\$fitted)
else {
tt <- terms(object)
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action,
xlev = object\$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts = object\$contrasts)
}
pred <- t(X %*% object\$sol[-(1:3),])
if(stepfun){
taus <- object\$sol[1,]
if(type == "Qhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) f[[i]] <- stepfun(taus, c(pred[1,i], pred[,i]))
} else
f <- stepfun(taus, c(pred[1,1], pred[,1]))
}
else if(type == "Fhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) {
o <- order(pred[,i])
f[[i]] <- stepfun(pred[,i],c(taus[1],taus))
}
} else
f <- stepfun(pred[,1],c(taus[1],taus))
}
else stop("type must be either 'Qhat' or 'Fhat'")
return(f)
}
return(pred)
}
"rearrange"  <- function (f, xmin, xmax)
# Revised Version September 11 2007.
{
if (is.list(f))
lapply(f, rearrange)
else {
if (!is.stepfun(f))
stop("Only stepfuns can be rearranged.\n")
call	<- attributes(f)\$call;
right <- call[match("right",names(call))]=="TRUE()"
x 	<- knots(f)
n	<- length(x)
if(missing(xmin)) xmin = x[1]
if(missing(xmax)) xmax = x[n]
x	<- x[(x >= xmin) & (x <= xmax)]
x	<- c(xmin, x, xmax)
n	<- length(x)
y	<- f(x)
o 	<- ifelse(rep(right,n-1), order(y[-1])+1, order(y[-n]))
x 	<- cumsum(c(x[1], diff(x)[o - right]))
y 	<- y[o]
y 	<- c(y[1], y, max(y))
stepfun(x, y, right = right)
}

}

# Function to compute regression quantiles using original simplex approach
# of Barrodale-Roberts/Koenker-d'Orey.  There are several options.
# The options are somewhat different than those available for the Frisch-
# Newton version of the algorithm, reflecting the different natures of the
# problems typically solved.  Succintly BR for "small" problems, FN for
# "large" ones.  Obviously, these terms are conditioned by available hardware.
#
# Basically there are two modes of use:
# 1.  For Single Quantiles:
#
#       if tau is between 0 and 1 then only one quantile solution is computed.
#
#       if ci = FALSE  then just the point estimate and residuals are returned
#		If the column dimension of x is 1 then ci is set to FALSE since
#		since the rank inversion method has no proper null model.
#       if ci = TRUE  then there are two options for confidence intervals:
#
#               1.  if iid = TRUE we get the original version of the rank
#                       inversion intervals as in Koenker (1994)
#               2.  if iid = FALSE we get the new version of the rank inversion
#                       intervals which accounts for heterogeneity across
#                       observations in the conditional density of the response.
#                       The theory of this is described in Koenker-Machado(1999)
#               Both approaches involve solving a parametric linear programming
#               problem, the difference is only in the factor qn which
#               determines how far the PP goes.  In either case one can
#               specify two other options:
#                       1. interp = FALSE returns two intervals an upper and a
#                               lower corresponding to a level slightly
#                               above and slightly below the one specified
#                               by the parameter alpha and dictated by the
#                               essential discreteness in the test statistic.
#				interp = TRUE  returns a single interval based on
#                               linear interpolation of the two intervals
#                               returned:  c.values and p.values which give
#                               the critical values and p.values of the
#                               upper and lower intervals. Default: interp = TRUE.
#                       2.  tcrit = TRUE uses Student t critical values while
#                               tcrit = FALSE uses normal theory ones.
# 2. For Multiple Quantiles:
#
#       if tau < 0 or tau >1 then it is presumed that the user wants to find
#       all of the rq solutions in tau, and the program computes the whole
#	quantile regression solution as a process in tau, the resulting arrays
#	containing the primal and dual solutions, betahat(tau), ahat(tau)
#       are called sol and dsol.  These arrays aren't printed by the default
#       print function but they are available as attributes.
#       It should be emphasized that this form of the solution can be
#	both memory and cpu quite intensive.  On typical machines it is
#	not recommended for problems with n > 10,000.
#	In large problems a grid of solutions is probably sufficient.
#
rq.fit.br <-
function (x, y, tau = 0.5, alpha = 0.1, ci = FALSE, iid = TRUE,
interp = TRUE, tcrit = TRUE)
{
tol <- .Machine\$double.eps^(2/3)
eps <- tol
big <- .Machine\$double.xmax
x <- as.matrix(x)
p <- ncol(x)
n <- nrow(x)
ny <- NCOL(y)
nsol <- 2
ndsol <- 2
if (qr(x)\$rank < p)
stop("Singular design matrix")
if (tau < 0 || tau > 1) {
nsol <- 3 * n
ndsol <- 3 * n
lci1 <- FALSE
qn <- rep(0, p)
cutoff <- 0
tau <- -1
}
else {
if (p == 1)
ci <- FALSE
if (ci) {
lci1 <- TRUE
if (tcrit)
cutoff <- qt(1 - alpha/2, n - p)
else cutoff <- qnorm(1 - alpha/2)
if (!iid) {
h <- bandwidth.rq(tau, n, hs = TRUE)
bhi <- rq.fit.br(x, y, tau + h, ci = FALSE)
bhi <- coefficients(bhi)
blo <- rq.fit.br(x, y, tau - h, ci = FALSE)
blo <- coefficients(blo)
dyhat <- x %*% (bhi - blo)
if (any(dyhat <= 0)) {
pfis <- (100 * sum(dyhat <= 0))/n
warning(paste(pfis, "percent fis <=0"))
}
f <- pmax(eps, (2 * h)/(dyhat - eps))
qn <- rep(0, p)
for (j in 1:p) {
qnj <- lm(x[, j] ~ x[, -j] - 1, weights = f)\$resid
qn[j] <- sum(qnj * qnj)
}
}
else qn <- 1/diag(solve(crossprod(x)))
}
else {
lci1 <- FALSE
qn <- rep(0, p)
cutoff <- 0
}
}
z <- .Fortran("rqbr", as.integer(n), as.integer(p), as.integer(n +
5), as.integer(p + 3), as.integer(p + 4), as.double(x),
as.double(y), as.double(tau), as.double(tol), flag = as.integer(1),
coef = double(p), resid = double(n), integer(n), double((n +
5) * (p + 4)), double(n), as.integer(nsol), as.integer(ndsol),
sol = double((p + 3) * nsol), dsol = double(n * ndsol),
lsol = as.integer(0), h = integer(p * nsol), qn = as.double(qn),
cutoff = as.double(cutoff), ci = double(4 * p), tnmat = double(4 *
p), as.double(big), as.logical(lci1), PACKAGE = "quantreg")
if (z\$flag != 0)
warning(switch(z\$flag, "Solution may be nonunique", "Premature end - possible conditioning problem in x"))
if (tau < 0 || tau > 1) {
sol <- matrix(z\$sol[1:((p + 3) * z\$lsol)], p + 3)
dsol <- matrix(z\$dsol[1:(n * z\$lsol)], n)
vnames <- dimnames(x)[[2]]
dimnames(sol) <- list(c("tau", "Qbar", "Obj.Fun", vnames),
NULL)
return(list(sol = sol, dsol = dsol))
}
if (!ci) {
coef <- z\$coef
dual <- z\$dsol[1:n]
names(coef) <- dimnames(x)[[2]]
return(list(coefficients = coef, x = x, y = y, residuals = y - x %*% z\$coef,
dual = dual))
}
if (interp) {
Tn <- matrix(z\$tnmat, nrow = 4)
Tci <- matrix(z\$ci, nrow = 4)
Tci[3, ] <- Tci[3, ] + (abs(Tci[4, ] - Tci[3, ]) * (cutoff -
abs(Tn[3, ])))/abs(Tn[4, ] - Tn[3, ])
Tci[2, ] <- Tci[2, ] - (abs(Tci[1, ] - Tci[2, ]) * (cutoff -
abs(Tn[2, ])))/abs(Tn[1, ] - Tn[2, ])
Tci[2, ][is.na(Tci[2, ])] <- -big
Tci[3, ][is.na(Tci[3, ])] <- big
coefficients <- cbind(z\$coef, t(Tci[2:3, ]))
vnames <- dimnames(x)[[2]]
cnames <- c("coefficients", "lower bd", "upper bd")
dimnames(coefficients) <- list(vnames, cnames)
residuals <- y - drop(x %*% z\$coef)
return(list(coefficients = coefficients, residuals = residuals))
}
else {
Tci <- matrix(z\$ci, nrow = 4)
coefficients <- cbind(z\$coef, t(Tci))
residuals <- y - drop(x %*% z\$coef)
vnames <- dimnames(x)[[2]]
cnames <- c("coefficients", "lower bound", "Lower Bound",
"upper bd", "Upper Bound")
dimnames(coefficients) <- list(vnames, cnames)
c.values <- t(matrix(z\$tnmat, nrow = 4))
c.values <- c.values[, 4:1]
dimnames(c.values) <- list(vnames, cnames[-1])
p.values <- if (tcrit)
matrix(pt(c.values, n - p), ncol = 4)
else matrix(pnorm(c.values), ncol = 4)
dimnames(p.values) <- list(vnames, cnames[-1])
list(coefficients = coefficients, residuals = residuals,
c.values = c.values, p.values = p.values)
}
}

"rq.fit.fnb" <-
function (x, y, tau = 0.5, beta = 0.99995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method.  Set tau in (0,1)")
rhs <- (1 - tau) * apply(x, 2, sum)
d   <- rep(1,n)
u   <- rep(1,n)
wn <- rep(0,10*n)
wn[1:n] <- (1-tau) #initial value of dual solution
z <- .Fortran("rqfnb", as.integer(n), as.integer(p), a = as.double(t(as.matrix(x))),
c = as.double(-y), rhs = as.double(rhs), d = as.double(d),as.double(u),
beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1),PACKAGE= "quantreg")
if (z\$info != 0)
stop(paste("Error info = ", z\$info, "in stepy: singular design"))
coefficients <- -z\$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
list(coefficients=coefficients, tau=tau, residuals=residuals)
}

"rq.fit.fnc" <-
function (x, y, R, r, tau = 0.5, beta = 0.9995, eps = 1e-06)
{
n1 <- length(y)
n2 <- length(r)
p <- ncol(x)
if (n1 != nrow(x))
stop("x and y don't match n1")
if (n2 != nrow(R))
stop("R and r don't match n2")
if (p != ncol(R))
stop("R and x don't match p")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method.  Set tau in (0,1)")
rhs <- (1 - tau) * apply(x, 2, sum)
u <- rep(1, max(n1,n2)) #upper bound vector and scratch vector
wn1 <- rep(0, 9 * n1)
wn1[1:n1] <- (1 - tau) #store the values of x1
wn2 <- rep(0, 6 * n2)
wn2[1:n2] <- 1 #store the values of x2
z <- .Fortran("rqfnc", as.integer(n1), as.integer(n2), as.integer(p),
a1 = as.double(t(as.matrix(x))), c1 = as.double(-y),
a2 = as.double(t(as.matrix(R))), c2 = as.double(-r),
rhs = as.double(rhs), d1 = double(n1), d2 = double(n2),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn1 = as.double(wn1), wn2 = as.double(wn2), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
if (z\$info != 0)
stop(paste("Error info = ", z\$info, "in stepy2: singular design"))
coefficients <- -z\$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z\$it.count
list(coefficients=coefficients, tau=tau, residuals=residuals, it = it.count)
}
function (x, y, tau = 0.5, alpha = 3.2, lambda = 1, start = "rq", beta = 0.9995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method.  Set tau in (0,1)")
if(length(lambda) == 1)
lambda <- c(0,rep(lambda,p-1))
if(length(lambda) != p)
stop(paste("lambda must be either of length ",p," or length one"))
if(any(lambda < 0))
stop("negative lambdas disallowed")
R <- diag(lambda,nrow = length(lambda))
R <- R[which(lambda != 0),, drop = FALSE]
r <- rep(0,nrow(R))
X <- rbind(x,  R)
Y <- c(y, r)
N <- length(Y)
rhs <- (1 - tau) * apply(x, 2, sum) + apply(R,2,sum)
dscad <- function(x, a = 3.7, lambda = 2){
lambda *  sign(x) *  (abs(x) <= lambda) +
sign(x) * (a * lambda -  abs(x)) / (a - 1) *
(abs(x) <= a * lambda) * (abs(x) > lambda)
}
binit <- switch(start,
rq =   rq.fit.fnb(x, y, tau = tau)\$coef[-1],
lasso = rq.fit.lasso(x, y, tau = tau, lambda = lambda)\$coef[-1]
)
coef <- rep(.Machine\$double.xmax,p)
it <- 0
while(sum(abs(binit - coef[-1])) > eps){
it <- it + 1
d <- rep(1, N)
u <- rep(1, N)
wn <- rep(0, 10 * N)
wn[1:N] <- c(rep((1 - tau),n),rep(.5,nrow(R)))
binit <- coef[-1]
z <- .Fortran("rqfnb", as.integer(N), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-Y), vrhs = as.double(vrhs), d = as.double(d),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
coef <- -z\$wp[1:p]
}
if (z\$info != 0)
stop(paste("Error info = ", z\$info, "in stepy2: singular design"))
coefficients <- -z\$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z\$it.count
list(coefficients=coefficients, residuals=residuals, tau = tau,
lambda = lambda, it = it.count)
}

"rq.fit.lasso" <-
function (x, y, tau = 0.5, lambda = 1, beta = 0.9995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if(length(lambda) == 1)
lambda <- c(0,rep(lambda,p-1))
if(length(lambda) != p)
stop(paste("lambda must be either of length ",p," or length one"))
if(any(lambda < 0))
stop("negative lambdas disallowed")
R <- diag(lambda,nrow = length(lambda))
R <- R[which(lambda != 0),, drop = FALSE]
r <- rep(0,nrow(R))
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method.  Set tau in (0,1)")
X <- rbind(x, R)
Y <- c(y, r)
N <- length(Y)
rhs <- (1 - tau) * apply(x, 2, sum) + 0.5 * apply(R,2,sum)
d <- rep(1, N)
u <- rep(1, N)
wn <- rep(0, 10 * N)
wn[1:N] <- rep(0.5, N)
z <- .Fortran("rqfnb", as.integer(N), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-Y), rhs = as.double(rhs), d = as.double(d),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p), aa = double(p *
p), it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
if (z\$info != 0)
stop(paste("Error info = ", z\$info, "in stepy2: singular design"))
coefficients <- -z\$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z\$it.count
list(coefficients=coefficients, residuals=residuals, tau = tau,
lambda = lambda, it = it.count)
}

"rq.fit.pfn" <-
# This is an implementation (purely in R) of the preprocessing phase
# of the rq algorithm described in Portnoy and Koenker, Statistical
# Science, (1997) 279-300.  In this implementation it can be used
# as an alternative method for rq() by specifying method="pfn"
# It should probably be used only on very large problems and then
# only with some caution.  Very large in this context means roughly
# n > 100,000.  The options are described in the paper, and more
# explicitly in the code.  Again, it would be nice perhaps to have
# this recoded in a lower level language, but in fact this doesn't
# seem to make a huge difference in this case since most of the work
# is already done in the rq.fit.fnb calls.
#
function(x, y, tau = 0.5,  Mm.factor = 0.8,
max.bad.fixup = 3, eps = 1e-6)
{
#rq function for n large --
n <- length(y)
if(nrow(x) != n)
stop("x and y don't match n")
if(tau < 0 | tau > 1)
stop("tau outside (0,1)")
p <- ncol(x)
m <- round(((p + 1) * n)^(2/3))
not.optimal <- TRUE
while(not.optimal) {
if(m < n)
s <- sample(n, m)
else {
b <- rq.fit.fnb(x, y, tau = tau,  eps = eps)\$coef
break
}
xx <- x[s,  ]
yy <- y[s]
z <- rq.fit.fnb(xx, yy, tau = tau,  eps = eps)
xxinv <- solve(chol(crossprod(xx)))
band <- sqrt(((x %*% xxinv)^2) %*% rep(1, p))
#sqrt(h<-ii)
r <- y - x %*% z\$coef
M <- Mm.factor * m
lo.q <- max(1/n, tau - M/(2 * n))
hi.q <- min(tau + M/(2 * n), (n - 1)/n)
kappa <- quantile(r/pmax(eps, band), c(lo.q, hi.q))
sl <- r < band * kappa[1]
su <- r > band * kappa[2]
xx <- x[!su & !sl,  ]
yy <- y[!su & !sl]
if(any(sl)) {
glob.x <- c(t(x[sl,  , drop = FALSE]) %*% rep(
1, sum(sl)))
glob.y <- sum(y[sl])
xx <- rbind(xx, glob.x)
yy <- c(yy, glob.y)
}
if(any(su)) {
ghib.x <- c(t(x[su,  , drop = FALSE]) %*% rep(
1, sum(su)))
ghib.y <- sum(y[su])
xx <- rbind(xx, ghib.x)
yy <- c(yy, ghib.y)
}
z <- rq.fit.fnb(xx, yy, tau = tau,  eps = eps)
b <- z\$coef
r <- y - x %*% b
su.bad <- (r < 0) & su
sl.bad <- (r > 0) & sl
M) {
warning("Too many fixups:  doubling m")
m <- 2 * m
break
}
}
else not.optimal <- FALSE
}
}
coefficients <- b
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% b
return(list(coefficients=coefficients, tau=tau,
residuals=residuals))
}

"rq.wfit" <-
function(x, y, tau = 0.5, weights, method = "br",  ...)
{
if(any(weights < 0))
stop("negative weights not allowed")
if(length(tau) > 1) {
tau <- tau[1]
warning("Multiple taus not allowed in rq.wfit: solution restricted to first element")
}
contr <- attr(x, "contrasts")
wx <- x * weights
wy <- y * weights
fit <- switch(method,
fn = rq.fit.fnb(wx, wy, tau = tau, ...),
fnb = rq.fit.fnb(wx, wy, tau = tau, ...),
br = rq.fit.br(wx, wy, tau = tau, ...),
fnc = rq.fit.fnc(wx, wy, tau = tau, ...),
pfn = rq.fit.pfn(wx, wy, tau = tau, ...), {
what <- paste("rq.fit.", method, sep = "")
if(exists(what, mode = "function"))
(get(what, mode = "function"))(x, y, ...)
else stop(paste("unimplemented method:", method))
}
)
if(length(fit\$sol))
fit\$fitted.values <- x %*% fit\$sol[-(1:3),]
else
fit\$fitted.values <- x %*% fit\$coef
fit\$residuals <- y - fit\$fitted.values
fit\$contrasts <- attr(x, "contrasts")
fit\$weights <- weights
fit
}

"summary.rqs" <-
function (object, ...) {
taus <- object\$tau
xsum <- as.list(taus)
for(i in 1:length(taus)){
xi <- object
xi\$coefficients <- xi\$coefficients[,i]
xi\$residuals <- xi\$residuals[,i]
xi\$tau <- xi\$tau[i]
class(xi) <- "rq"
xsum[[i]] <- summary(xi, ...)
if(class(object)[1] == "dynrqs"){
class(xsum[[1]]) <- c("summary.dynrq", "summary.rq")
if(i == 1) xsum[[1]]\$model <- object\$model
}
}
class(xsum) <- "summary.rqs"
if(class(object)[1] == "dynrqs")
class(xsum) <- c("summary.dynrqs", "summary.rqs")
xsum
}
"logLik.rq" <- function(object,  ...){
n <- length(object\$residuals)
p <- length(object\$coefficients)
pen <- (length(object\$lambda) > 0)
tau <- object\$tau
fid <- object\$rho
val <- n * (log(tau * (1-tau)) - 1 - log(fid/n))
attr(val,"n") <- n
if(pen){
if(!hasArg(edfThresh)) edfThresh <- 0.0001
attr(val,"df") <- sum(abs(object\$coefficients) > edfThresh)
}
else  attr(val,"df") <- p
class(val) <- "logLik"
val
}
"logLik.rqs" <- function(object, ...){
n <- nrow(object\$residuals)
p <- nrow(object\$coefficients)
pen <- (length(object\$lambda) > 0)
tau <- object\$tau
fid <- object\$rho
val <- n * (log(tau * (1-tau)) - 1 - log(fid/n))
attr(val,"n") <- n
if(pen){
if(!hasArg(edfThresh)) edfThresh <- 0.0001
attr(val,"df") <- sum(abs(object\$coefficients) > edfThresh)
}
else  attr(val,"df") <- p
class(val) <- "logLik"
val
}
"AIC.rq" <- function(object, ... , k = 2){
v <- logLik(object)
if(k <= 0)
k <- log(attr(v,"n"))
val <- AIC(v, k = k)
attr(val,"edf") <- attr(v,"df")
val
}
"extractAIC.rq"  <- function(fit, scale, k=2, ...){
aic <- AIC(fit,k)
edf <- attr(aic, "edf")
c(edf, aic)
}

"AIC.rqs" <- function(object, ... , k = 2){
v <- logLik(object)
if(k <= 0)
k <- log(attr(v,"n"))
val <- AIC(v, k = k)
attr(val,"edf") <- attr(v,"df")
val
}

"summary.rq" <-
# This function provides  methods for summarizing the output of the
# rq command. In this instance, "summarizing" means essentially provision
# of either standard errors, or confidence intervals for the rq coefficents.
# Since the preferred method for confidence intervals is currently the
# rank inversion method available directly from rq() by setting ci=TRUE, with br=TRUE.
# these summary methods are intended primarily for comparison purposes
# and for use on large problems where the parametric programming methods
# of rank inversion are prohibitively memory/time consuming.  Eventually
# iterative versions of rank inversion should be developed that would
# employ the Frisch-Newton approach.
#
# Object is the result of a call to rq(), and the function returns a
# table of coefficients, standard errors, "t-statistics", and p-values, and, if
# covariance=TRUE a structure describing the covariance matrix of the coefficients,
# i.e. the components of the Huber sandwich.
#
# There are five options for "se":
#
#	1.  "rank" strictly speaking this doesn't produce a "standard error"
#		at all instead it produces a coefficient table with confidence
#		intervals for the coefficients based on inversion of the
#		rank test described in GJKP and Koenker (1994).
#	2.  "iid" which presumes that the errors are iid and computes
#		an estimate of the asymptotic covariance matrix as in KB(1978).
#	3.  "nid" which presumes local (in tau) linearity (in x) of the
#		the conditional quantile functions and computes a Huber
#		sandwich estimate using a local estimate of the sparsity.
#	4.  "ker" which uses a kernel estimate of the sandwich as proposed
#		by Powell.
#	5.  "boot" which uses a bootstrap method:
#		"xy"	uses xy-pair method
#		"wxy"	uses weighted (generalized) method
#		"pwy"	uses the parzen-wei-ying method
#		"mcmb"	uses the Markov chain marginal bootstrap method
#
#
function (object, se = NULL, covariance = FALSE, hs = TRUE, ...)
{
if(object\$method == "lasso")
stop("no inference for lasso'd rq fitting: try rqss (if brave)")
mt <- terms(object)
m <- model.frame(object)
y <- model.response(m)
x <- model.matrix(mt, m, contrasts = object\$contrasts)
wt <- as.vector(model.weights(object\$model))
tau <- object\$tau
eps <- .Machine\$double.eps^(1/2)
coef <- coefficients(object)
if (is.matrix(coef))
coef <- coef[, 1]
vnames <- dimnames(x)[[2]]
resid <- object\$residuals
n <- length(resid)
p <- length(coef)
rdf <- n - p
if (!is.null(wt)) {
resid <- resid * wt
x <- x * wt
y <- y * wt
}
if (is.null(se)) {
if (n < 1001 & covariance == FALSE)
se <- "rank"
else se <- "nid"
}
if (se == "rank") {
f <- rq.fit.br(x, y, tau = tau, ci = TRUE, ...)
}
if (se == "iid") {
xxinv <- diag(p)
xxinv <- backsolve(qr(x)\$qr[1:p, 1:p,drop=FALSE], xxinv)
xxinv <- xxinv %*% t(xxinv)
pz <- sum(abs(resid) < eps)
h <- max(p + 1, ceiling(n * bandwidth.rq(tau, n, hs = hs)))
ir <- (pz + 1):(h + pz + 1)
ord.resid <- sort(resid[order(abs(resid))][ir])
xt <- ir/(n - p)
sparsity <- rq(ord.resid ~ xt)\$coef[2]
cov <- sparsity^2 * xxinv * tau * (1 - tau)
scale <- 1/sparsity
serr <- sqrt(diag(cov))
}
else if (se == "nid") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1:  error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0:  error in summary.rq")
bhi <- rq.fit.fnb(x, y, tau = tau + h)\$coef
blo <- rq.fit.fnb(x, y, tau = tau - h)\$coef
dyhat <- x %*% (bhi - blo)
if (any(dyhat <= 0))
warning(paste(sum(dyhat <= 0), "non-positive fis"))
f <- pmax(0, (2 * h)/(dyhat - eps))
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)\$qr[1:p, 1:p,drop=FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "ker") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1:  error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0:  error in summary.rq")
uhat <- c(y - x %*% coef)
h <- (qnorm(tau + h) - qnorm(tau - h))*
min(sqrt(var(uhat)), ( quantile(uhat,.75)- quantile(uhat, .25))/1.34 )
f <- dnorm(uhat/h)/h
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)\$qr[1:p, 1:p,drop=FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "boot") {
B <- boot.rq(x, y, tau, ...)
cov <- cov(B)
serr <- sqrt(diag(cov))
}
if( se == "rank"){
coef <- f\$coef
}
else {
coef <- array(coef, c(p, 4))
dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value",
"Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if (rdf > 0)
2 * (1 - pt(abs(coef[, 3]), rdf))
else NA
}
object <- object[c("call", "terms")]
if (covariance == TRUE) {
object\$cov <- cov
if(se == "iid") object\$scale <- scale
if(se %in% c("nid", "ker")) {
object\$Hinv <- fxxinv
object\$J <- crossprod(x)
object\$scale <- scale
}
else if (se == "boot") {
object\$B <- B
}
}
object\$coefficients <- coef
object\$residuals <- resid
object\$rdf <- rdf
object\$tau <- tau
class(object) <- "summary.rq"
object
}

akj <- function(x,
z = seq(min(x), max(x), length = 2 * length(x)),
p = rep(1/length(x), length(x)),
h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0)
{
nx <- length(x)
stopifnot(is.numeric(x),
length(p) == nx,
any((iker1 <- as.integer(iker1)) == 0:1))
nz <- length(z)
if(is.unsorted(x))
x <- sort(x)
.Fortran("akj",
as.double(x),
as.double(z),
as.double(p),
iker1,
dens  = double(nz),
psi   = double(nz),
score = double(nz),
as.integer(nx),
as.integer(nz),
h = as.double(h),
as.double(alpha),
as.double(kappa),
double(nx),
PACKAGE = "quantreg")[c("dens", "psi", "score", "h")]
}

"lm.fit.recursive" <-
function(X, y, int = TRUE)
{
if(int)
X <- cbind(1, X)
p <- ncol(X)
n <- nrow(X)
D <- qr(X[1:p,  ])
A <- qr.coef(D, diag(p))
A[is.na(A)] <- 0
A <- crossprod(t(A))
Ax <- rep(0, p)
b <- matrix(0, p, n)
b[, p] <- qr.coef(D, y[1:p])
b[is.na(b)] <- 0
z <- .Fortran( "rls",
as.integer(n),
as.integer(p),
as.double(t(X)),
as.double(y),
b = as.double(b),
as.double(A),
as.double(Ax),
PACKAGE = "quantreg")
bhat <- matrix(z\$b, p, n)
return(bhat)
}

"rq.fit.hogg" <-
function (x, y, taus = c(.1,.3,.5), weights = c(.7,.2,.1),
R= NULL, r = NULL, beta = 0.99995, eps = 1e-06)
{
n <- length(y)
n2 <- nrow(R)
m <- length(taus)
p <- ncol(x)+m
if (n != nrow(x))
stop("x and y don't match n")
if (m != length(weights))
stop("taus and weights differ in length")
if (any(taus < eps) || any(taus > 1 - eps))
stop("taus outside (0,1)")
W <- diag(weights)
if(m == 1) W <- weights
x <- as.matrix(x)
X <- cbind(kronecker(W,rep(1,n)),kronecker(weights,x))
y <- kronecker(weights,y)
rhs <- c(weights*(1 - taus)*n, sum(weights*(1-taus)) * apply(x, 2, sum))
if(n2!=length(r))
stop("R and r of incompatible dimension")
if(ncol(R)!=p)
stop("R and X of incompatible dimension")
d <- rep(1, m*n)
u <- rep(1, m*n)
if(length(r)){
wn1 <- rep(0, 10 * m*n)
wn1[1:(m*n)] <- .5
wn2 <- rep(0,6*n2)
wn2[1:n2] <- 1
z <- .Fortran("rqfnc", as.integer(m*n), as.integer(n2), as.integer(p),
a1 = as.double(t(as.matrix(X))), c1 = as.double(-y),
a2 = as.double(t(as.matrix(R))), c2 = as.double(-r),
rhs = as.double(rhs), d1 = double(m*n), d2 = double(n2),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn1 = as.double(wn1), wn2 = as.double(wn2), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
}
else{
wn <- rep(0, 10 * m*n)
wn[1:(m*n)] <- .5
z <- .Fortran("rqfnb", as.integer(m*n), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-y), rhs = as.double(rhs), d = as.double(d), as.double(u),
beta = as.double(beta), eps = as.double(eps), wn = as.double(wn),
wp = double((p + 3) * p), it.count = integer(2), info = integer(1),
PACKAGE = "quantreg")
}
if (z\$info != 0)
warning(paste("Info = ", z\$info, "in stepy: singular design: iterations ", z\$it.count[1]))
coefficients <- -z\$wp[1:p]
if(any(is.na(coefficients)))stop("NA coefs:  infeasible problem?")
list(coefficients = coefficients, nit = z\$it.count, flag = z\$info)
}
``````