"trq.fit"<-
function(x, y, a1 = 0.1, a2, z, int = T, method = "primal", tol = 0.0001)
{
#compute trimmed regression quantiles
#z is the rq strcture 
	if(missing(a2)) a2 <- a1
	if(a1 < 0 | a2 < 0)
		stop("trimming proportion negative.")
	if(a1 + a2 - 1 > tol)
		stop("trimming proportion greater than 1.")
	if(method!="primal" & method!="dual")
		stop("invalid method: should be 'primal' or 'dual'.")
	x <- as.matrix(x)
	if(missing(z))
		z <- rq.fit.br(x, y, tau = -1)
	p <- z$sol[1,  ]
	q <- matrix(z$sol[ - c(1:3),  ], nrow(z$sol) - 3, ncol(z$sol))
	n <- nrow(z$dsol)
	s <- NULL
	if(length(dimnames(x)[[2]]) == 0)
		dimnames(x) <- list(NULL, paste("X", 1:(nrow(q) - 1), sep = "")
			)
	if(int) {
		x <- cbind(1, x)
		dimnames(x)[[2]][1] <- "Intercept"
	}
	xbar <- apply(x, 2, "mean")
	xxinv <- solve(t(x) %*% x)
	if(abs(a1 + a2 - 1) <= tol) {
#single quantile case
		i <- sum(p < a1)
		s$coef <- q[, i]
		names(s$coef) <- dimnames(x)[[2]]
		s$resid <- y - x %*% s$coef
		PI <- 3.14159
		x0 <- qnorm(a1)
		d0 <- (1/sqrt(2 * PI)) * exp( - (x0^2/2))
		d0 <- ((4.5 * d0^4)/(2 * x0^2 + 1)^2)^0.2
		d <- d0 * (length(s$resid) - length(s$coef))^(-0.2)
		if(d > min(a1, 1 - a1))
			d <- min(a1, 1 - a1)
		s$d <- d
		i <- sum(p < a1 + d)
		j <- sum(p < a1 - d)
		shat <- as.numeric(xbar %*% t(q[, i] - q[, j]))/(2 * d)
		s$int <- int
		s$v <- a1 * (1 - a1) * shat^2
		s$cov <- s$v * xxinv
	}
	else {
#real trimming
		p1 <- p[-1]
		f <- 1/(1 - a1 - a2)
		d <- pmax((pmin(p1, 1 - a2) - c(a1, pmax(p1[ - length(p1)],
			a1))), 0)
		if(method == "primal") {
			s$coef <- q[, 1:length(p1)] %*% d * f
			s$resid <- y - x %*% s$coef
			s$int <- int
		}
		else {
#Jureckova-Gutenbrunner trimmed least squares
			i <- max(1, sum(p < a1))
			g <- (z$dsol[, i + 1] - z$dsol[, i])/(p[i + 1] - p[
				i])
			wa <- z$dsol[, i] + (a1 - p[i]) * g
			j <- sum(p < 1 - a2)
			g <- (z$dsol[, j + 1] - z$dsol[, j])/(p[j + 1] - p[
				j])
			wb <- z$dsol[, j] + (1 - a2 - p[j]) * g
			wt <- wa - wb
			if(min(wt) < 0)
				warning("some weights negative!")
			s <- lsfit(x, y, abs(wt), int = F)
		}
#now compute covariance matrix estimate
		mu <- xbar %*% s$coef
		v <- d %*% (z$sol[2, 1:length(d)] - mu)^2
		k <- qrq(z, c(a1, a2)) - mu
		v <- v + a1 * k[1]^2 + a2 * k[2]^2 + (a1 * k[1] + a2 * k[2])^2
		names(s$coef) <- dimnames(x)[[2]]
		s$v <- as.vector(f^2 * v)
		s$cov <- s$v * xxinv
	}
	class(s) <- "trq"
	s
}
"print.trq"<-
function(object, digits = 4)
{
	n <- length(object$resid)
	p <- length(object$coef)
	options(warn = -1)
	if(object$int) {
		df.num <- p - 1
		fstat <- c(t(object$coef[-1]) %*% solve(object$cov[-1, -1]) %*%
			(object$coef[-1]))/df.num
	}
	else {
		df.num <- p
		fstat <- t(object$coef) %*% solve(object$cov) %*% (object$coef
			)/df.num
	}
	pvalue <- 1 - pf(fstat, df.num, (n - p))
	regstat <- c(sqrt(object$v), n, fstat, df.num, (n - p), pvalue)
	names(regstat) <- c("rse", "n", "F.stat", "df.num", "df.den", "p.value"
		)
	err <- sqrt(diag(object$cov))
	tstat <- c(object$coef/err)
	tabcoef <- cbind(object$coef, err, tstat, 2 * (1 - pt(abs(tstat),
		n - p)))
	dimnames(tabcoef) <- list(names(object$coef), c("coef", "std.err",
		"t.stat", "p.value"))
	options(warn = 0)
	print(round(tabcoef, digits))
	cat(paste("Winsorized Standard Error of Regression= ", format(round(
		sqrt(object$v), digits)), "\n", "N = ", format(n), 
		",  F-statistic = ", format(round(fstat, digits)), " on ",
		format(df.num), " and ", format((n - p)), " df, ", "p-value = ",
		format(round(pvalue, digits)), "\n\n", sep = ""))
	invisible(list(summary = regstat, coef.table = tabcoef))
}
"qrq"<-
function(s, a)
{
#computes linearized quantiles from rq data structure
#v is the rq structure e.g. rq(x,y)
#a is a vector of quantiles required
        if(min(a) < 0 | max(a) > 1) stop("alphas out of range [0,1]")
        r <- s$sol[1,  ]
        q <- s$sol[2,  ]
        q <- c(q[1], q)
        J <- length(r)
        r <- c(0, (r[1:J - 1] + r[2:J])/2, 1)
        u <- rep(0, length(a))
        for(k in 1:length(a)) {
                i <- sum(r < a[k])
                w <- (a[k] - r[i])/(r[i + 1] - r[i])
                u[k] <- w * q[i + 1] + (1 - w) * q[i]
        }
        u
}
