Running 'panel.R' [2s]
  Running 'rq.R' [148s]
  Running 'rq.fit.panel.R' [165s]
  Running 'run-demos.R' [191s]
Running the tests in 'tests/rq.R' failed.
Complete output:
  > ## This is just from ?anova.rq (extended)
  > library(quantreg)
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > data(barro)
  > fit0 <- rq(y.net ~  lgdp2 + fse2 + gedy2 , data = barro)
  > fit1 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro)
  > 
  > a01 <- anova(fit1,fit0)
  > a01
  Quantile Regression Analysis of Deviance Table
  
  Model 1: y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2
  Model 2: y.net ~ lgdp2 + fse2 + gedy2
    Df Resid Df F value    Pr(>F)    
  1  2      155  18.879 4.596e-08 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  > 
  > fit2 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,
  + 	   tau = 0.75)
  > fit3 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,
  + 	   tau = 0.25)
  > 
  > a123  <- anova(fit1,fit2, fit3)
  Warning message:
  In summary.rq(x, se = se, R = R, covariance = TRUE) : 1 non-positive fis
  > a.123 <- anova(fit1,fit2, fit3, joint=FALSE)
  Warning message:
  In summary.rq(x, se = se, R = R, covariance = TRUE) : 1 non-positive fis
  > a.123
  Quantile Regression Analysis of Deviance Table
  
  Model: y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2
  Tests of Equality of Distinct Slopes: tau in {  0.5 0.75 0.25  }
  
         Df Resid Df F value  Pr(>F)  
  lgdp2   2      481  1.0656 0.34535  
  fse2    2      481  2.6398 0.07241 .
  gedy2   2      481  0.7862 0.45614  
  Iy2     2      481  0.0447 0.95632  
  gcony2  2      481  0.0653 0.93675  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  > 
  > AE <- function(x,y) all.equal(x, y, tol = 1e-5)
  > ##                                  ---------- {giving a bit more digits below}
  > stopifnot(
  + 	  AE(100 * unname(coef(fit0)),
  + 	     c(-0.74679759, 0.46539963, 0.15902838, -36.619915))
  + 	  ,
  + 	  AE(unlist(a01$table),
  + 	     c(ndf=2, ddf=155, Tn = 18.878717, pvalue= 4.6e-08))
  + 	  ,
  + 	  AE(100* unname(coef( fit2 )),
  + 	     c(13.103018, -1.4885239, -0.026452369,
  + 	       0.3999839, 14.526663, -13.504643))
  + 	  ,
  + 	  AE(100* unname(coef( fit3 )),
  + 	     c(6.0860719, -0.88350554, 0.24596781,
  + 	       -14.962498, 15.592489, -15.861804))
  + 	  ,
  + 	  AE(unlist(a123$table),
  + 	     c(ndf = 10, ddf = 473, Tn = 1.80385526, pvalue=0.0575117558))
  + 	  ,
  + 	  AE(a.123$table[,"Tn"],
  + 	     c(1.0655561, 2.6398508, 0.78623238, 0.04467014, 0.065344348))
  + 	  )
  > 
  > proc.time()
     user  system elapsed 
     1.45    0.14    1.57 
Running the tests in 'tests/rq.fit.panel.R' failed.
Complete output:
  > rq.fit.panel <- function(X,y,s,w=c(.25,.5,.25),taus=(1:3)/4,lambda = 1){
  + # prototype function for panel data fitting of QR models
  + # the matrix X is assumed to contain an intercept
  + # the vector s is a strata indicator assumed (so far) to be a one-way layout
  + # NB:  
  + # 0.  This is an altered version from that originally posted -- the definition
  + #     of the rhs vector now incorporates a factor of 1/2 for the penalty.
  + # 1.  The value of the shrinkage parameter lambda is an open research problem in
  + #   the simplest homogneous settings it should be the ratio of the scale parameters
  + #   of the fixed effects and the idiocyncratic errors
  + # 2.  On return the coefficient vector has m*p + n elements where m is the number
  + # quantiles being estimated, p is the number of colums of X, and n is the
  + # number of distinct values of s.  The first m*p coefficients are the 
  + # slope estimates, and the last n are the "fixed effects"
  + # 3.  Like all shrinkage (regularization) estimators, asymptotic inference is somewhat
  + # problematic... so the bootstrap is the natural first resort.
  +   require(SparseM)
  +   require(quantreg)
  +   K <- length(w)
  +   if(K != length(taus))
  +     stop("length of w and taus must match")
  +   X <- as.matrix(X)
  +   p <- ncol(X)
  +   n <- length(levels(as.factor(s)))
  +   N <- length(y)
  +   if(N != length(s) || N != nrow(X))
  +     stop("dimensions of y,X,s must match")
  +   Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
  +   Fidelity <- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
  +   Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
  +   D <- rbind(Fidelity,Penalty)
  +   y <- c(w %x% y,rep(0,n))
  +   a <- c((w*(1-taus)) %x% (t(X)%*%rep(1,N)),
  +     sum(w*(1-taus)) * (t(Z) %*% rep(1,N)) + lambda * rep(1/2,n))
  +   rq.fit.sfn(D,y,rhs=a)
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.25    0.01    0.21 
Running the tests in 'tests/run-demos.R' failed.
Complete output:
  > library(quantreg)
  Loading required package: SparseM
  
  Attaching package: 'SparseM'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > 
  > (dDIR <- system.file("demo", package = "quantreg"))
  [1] "D:/RCompile/CRANpkg/lib/4.4/quantreg/demo"
  > set.seed(1) # since some demos randomly generate
  > 
  > cat("Running demos from package 'quantreg' : \n\n")
  Running demos from package 'quantreg' : 
  
  > for(ff in list.files(dDIR, pattern="\\.R$", full.names = TRUE)) {
  +    f <- basename(ff)
  +    cat("\n", f," :\n", paste(rep.int("-", nchar(f)), collapse=''),
  +        "\n", sep='')
  + 
  +    source(ff, echo = TRUE)
  + }
  
  Frank.R :
  -------
  
  > vFrank <- function(x, df, delta, u) -log(1 - (1 - 
  +     exp(-delta))/(1 + exp(-delta * pt(x, df)) * ((1/u) - 1)))/delta
  
  > FrankModel <- function(x, delta, mu, sigma, df, tau) {
  +     z <- qt(vFrank(x, df, delta, u = tau), df)
  +     mu + sigma * z
  + }
  
  > n <- 200
  
  > df <- 8
  
  > delta <- 8
  
  > set.seed(1989)
  
  > x <- sort(rt(n, df))
  
  > v <- vFrank(x, df, delta, u = runif(n))
  
  > y <- qt(v, df)
  
  > plot(x, y, pch = "o", col = "blue", cex = 0.25)
  
  > Dat <- data.frame(x = x, y = y)
  
  > us <- c(0.25, 0.5, 0.75)
  
  > for (i in 1:length(us)) {
  +     v <- vFrank(x, df, delta, u = us[i])
  +     lines(x, qt(v, df))
  + }
  
  > cfMat <- matrix(0, 3, length(us))
  
  > trace <- TRUE
  
  > trace <- FALSE
  
  > for (i in 1:length(us)) {
  +     tau <- us[i]
  +     cat("tau = ", format(tau), ".. ")
  +     fit <- nlrq(y ~ FrankModel(x, delta, mu, sigma, df = 8, t .... [TRUNCATED] 
  tau =  0.25 .. 
  tau =  0.5 .. 
  tau =  0.75 .. 
  
  > colnames(cfMat) <- names(coef(fit))
  
  > cfMat
          delta          mu     sigma
  [1,] 14.86275 -0.20513700 0.9135934
  [2,] 16.65111  0.03327799 0.9578400
  [3,] 11.99822  0.11557829 0.9464470
  
  KMvCRQ.R :
  --------
  
  > if (requireNamespace("survival", quietly = TRUE)) {
  +     n <- 100
  +     y <- rchisq(n, 3)
  +     c <- rchisq(n, 5)
  +     Y <- pmin(y, c)
  +     d <-  .... [TRUNCATED] 
  
  MCV.R :
  -----
  
  > MCV <- function(lambdas, formula, data, tau = 0.5, 
  +     k = 10) {
  +     F <- Munge(formula, lambdas = lambdas)
  +     f <- rqss(F, data, tau = tau) .... [TRUNCATED] 
  
  > set.seed(1729)
  
  > n <- 200
  
  > x <- sort(runif(n, 0, 20))
  
  > g0 <- function(x, tau) log(x) + 0.2 * (log(x))^3 + 
  +     log(x) * qnorm(tau)/4
  
  > y <- g0(x, runif(n))
  
  > D <- data.frame(y = y, x = x)
  
  > lams <- mcvs <- seq(0.02, 5, by = 0.2)
  
  > for (i in 1:length(mcvs)) mcvs[i] <- MCV(lams[i], 
  +     y ~ qss(x, lambda = lambdas[1]), D)
  
  > par(mfrow = c(1, 2))
  
  > plot(lams, mcvs, cex = 0.5, lwd = 2, type = "l", xlab = expression(lambda), 
  +     ylab = expression(MCV(lambda)))
  
  > lambdastar <- lams[which.min(mcvs)]
  
  > plot(x, y, cex = 0.5, col = "grey")
  
  > f <- rqss(y ~ qss(x, lambda = lambdastar), data = D)
  
  > plot(f, add = TRUE, lwd = 2)
  
  > lines(x, g0(x, 0.5), col = "red", lwd = 2)
  
  > text(10, 1, bquote(lambda == ~.(lambdastar)))
  
  Mammals.R :
  ---------
  
  > require(quantreg)
  
  > data(Mammals)
  
  > attach(Mammals)
  
  > x <- log(weight)
  
  > xx <- unique(x[order(x)])
  
  > y <- log(speed)
  
  > plot(x, y, xlab = "Weight in log(Kg)", ylab = "Speed in log(Km/hour)", 
  +     type = "n")
  
  > points(x[hoppers], y[hoppers], pch = "h", col = "red")
  
  > points(x[specials], y[specials], pch = "s", col = "blue")
  
  > others <- (!hoppers & !specials)
  
  > points(x[others], y[others], col = "black", cex = 0.75)
  
  > taus <- c(0.5, 0.9)
  
  > for (i in 1:length(taus)) {
  +     fit <- rqss(y ~ qss(x, lambda = 1, constraint = "C"), tau = taus[i])
  +     plot(fit, title = "Running Speed of Mam ..." ... [TRUNCATED] 
  
  > legend(4, 2, c("Median", "0.9 Quantile"), lty = 1, 
  +     col = 1:2, lwd = 1.5)
  
  > plot(fit, title = "Running Speed of Mammals", band = "both", 
  +     col = i, lwd = 1.5)
  
  > xy <- fit$qss[[1]]$xyz
  
  > xx <- xy[, 1]
  
  > yhat <- fit$coef[1] + xy[, 2]
  
  > g <- diff(yhat)/diff(xx)
  
  > plot(xx[-1], g, main = "Fitted Slopes of Running Speed", 
  +     xlab = "Weight in log(Kg)", ylab = "dlog(Speed) /dlog(Weight)")
  
  Mel.R :
  -----
  
  > require(splines)
  Loading required package: splines
  
  > if (interactive()) {
  +     oldpar <- par(ask = TRUE)
  +     data(MelTemp)
  +     x <- MelTemp[-3650]
  +     y <- MelTemp[-1]
  +     s <- (x < 40)
  +      .... [TRUNCATED] 
  
  Mel2.R :
  ------
  
  > require(splines)
  
  > if (interactive()) {
  +     oldpar <- par(ask = TRUE)
  +     data(MelTemp)
  +     x <- MelTemp[-3650]
  +     y <- MelTemp[-1]
  +     s <- (x < 40)
  +      .... [TRUNCATED] 
  
  Panel.R :
  -------
  
  > library(quantreg)
  
  > rq.fit.panel <- function(X, y, s, w = c(1/3, 1/3, 
  +     1/3), taus = c(0.25, 0.5, 0.75), lambda = 0) {
  +     require(SparseM)
  +     require(quantre .... [TRUNCATED] 
  
  > n <- 3
  
  > T <- 50
  
  > nT <- n * T
  
  > u1 <- rnorm(T)
  
  > u2 <- rnorm(T)
  
  > u3 <- rnorm(T)
  
  > x1 <- rnorm(T, 1, 0.85)
  
  > x2 <- rnorm(T, 4, 1)
  
  > x3 <- rnorm(T, 7, 1)
  
  > beta1 <- 1
  
  > beta2 <- 0.5
  
  > y1 <- 0 + beta1 * x1 + (beta2 * x1) * u1
  
  > y2 <- 4 + beta1 * x2 + (beta2 * x2) * u2
  
  > y3 <- 8 + beta1 * x3 + (beta2 * x3) * u3
  
  > plot(c(0, 9), c(0, 25), type = "n", xlab = expression(x[it]), 
  +     ylab = expression(y[it]))
  
  > points(x1, y1, pch = 15)
  
  > points(x2, y2, pch = 15, col = "blue")
  
  > points(x3, y3, pch = 15, col = "red")
  
  > legend(1, 17, paste("i = ", 1:3, sep = ""), pch = 15, 
  +     col = c("black", "blue", "red"))
  
  > ya <- c(y1, y2, y3)
  
  > xa <- c(x1, x2, x3)
  
  > taus <- c(0.25, 0.5, 0.75)
  
  > xx <- seq(min(xa), max(xa), 0.25)
  
  > f <- coef(rq(ya ~ xa, tau = taus))
  
  > yy <- cbind(1, xx) %*% f
  
  > for (i in 1:3) lines(xx, yy[, i], col = "grey")
  
  > s <- rep(1:n, rep(T, n))
  
  > fp <- rq.fit.panel(xa, ya, s)$coef
  
  > bhat <- fp[1:3]
  
  > fehat <- fp[4:6]
  
  > xx1 <- seq(min(x1), max(x1), 0.25)
  
  > for (i in 1:3) {
  +     yy1 <- fehat[1] + bhat[i] * xx1
  +     lines(xx1, yy1, col = "black")
  + }
  
  > xx2 <- seq(min(x2), max(x2), 0.25)
  
  > for (i in 1:3) {
  +     yy2 <- fehat[2] + bhat[i] * xx2
  +     lines(xx2, yy2, col = "blue")
  + }
  
  > xx3 <- seq(min(x3), max(x3), 0.25)
  
  > for (i in 1:3) {
  +     yy3 <- fehat[3] + bhat[i] * xx3
  +     lines(xx3, yy3, col = "red")
  + }
  
  Polson.R :
  --------
  
  > dgp <- function(n) {
  +     x <- 0:n/n
  +     y <- rnorm(n + 1, 5 * sin(2 * pi * x), 0.5 + exp(1.5 * sin(4 * 
  +         pi * x)))
  +     data.frame(x = .... [TRUNCATED] 
  
  > D <- dgp(1000)
  
  > plot(D$x, D$y, cex = 0.5)
  
  > taus <- 1:9/10
  
  > for (i in 1:length(taus)) plot(rqss(y ~ qss(x, lambda = 1/10), 
  +     tau = taus[i], data = D), rug = FALSE, add = TRUE)
  
  RB-r.R :
  ------
  
  > if (require(MASS)) {
  +     require(quantreg)
  +     X <- cbind(1, matrix(rnorm(500), 100, 5))
  +     y <- rnorm(100)
  +     R <- matrix(rnorm(18), 3, 6 .... [TRUNCATED] 
  Loading required package: MASS
  
  arqss.R :
  -------
  
  > n <- 2000
  
  > x <- 1:n/n
  
  > noise <- rgamma(n, 3, 1)
  
  > g0 <- function(x) sin(10 * x)
  
  > y <- g0(x) + noise
  
  > arqss <- function(x, y, tau, g0 = NULL) {
  +     g <- function(lam, y, x, tau) AIC(rqss(y ~ qss(x, lambda = lam), 
  +         tau = tau), k = -1)
  +    .... [TRUNCATED] 
  
  > arqss(x, y, 0.5, g0)
  
  cobar.R :
  -------
  
  > require(quantreg)
  
  > if (requireNamespace("interp")) {
  +     do.rgl <- interactive() && require(rgl)
  +     data(CobarOre)
  +     fit <- rqss(z ~ qss(cbind(x, y), lambda = .... [TRUNCATED] 
  Loading required namespace: interp
  
  combos.R :
  --------
  
  > H <- combos(20, 3)
  
  > if (!require("rgl", quietly = TRUE)) {
  +     warning("The package rgl is needed for plotting")
  + } else {
  +     if (interactive()) {
  +         plot3 .... [TRUNCATED] 
  
  cpoint.R :
  --------
  
  > x <- runif(200, 0, 10)
  
  > u <- rnorm(200)/5
  
  > y <- 1 + x - 0.5 * (x - 3) * (x > 3) + u
  
  > plot(y ~ x, cex = 0.5, col = "grey")
  
  > z <- rqss(y ~ qss(x, lambda = 10), tau = 0.5)
  
  > plot(z, col = "dark blue")
  
  > eps <- 1e-05
  
  > Nz <- abs(z$resid[1:200]) < eps
  
  > Nj <- abs(z$resid[201:398]) > eps
  
  > xx <- z$qss[[1]]$xyz[, 1]
  
  > yy <- z$coef[1] + z$qss[[1]]$xyz[, 2]
  
  > xj <- xx[3:200]
  
  > yj <- yy[3:200]
  
  > points(xx[Nz], yy[Nz], col = "green")
  
  > points(xj[Nj], yj[Nj], col = "red")
  
  > print(paste("Number of zero residuals:  ", sum(Nz)))
  [1] "Number of zero residuals:   4"
  
  > print(paste("Number of jumps in slope:  ", sum(Nj)))
  [1] "Number of jumps in slope:   2"
  
  > legend(6, 3, c("Derivative Jumps", "Zero residuals"), 
  +     pch = "o", col = c("red", "green"))
  
  crquis.R :
  --------
  
  > if (requireNamespace("survival", quietly = TRUE)) {
  +     data(uis)
  +     Surv <- survival::Surv
  +     fit <- crq(Surv(log(TIME), CENSOR) ~ ND1 + ND .... [TRUNCATED] 
  
  engel1.R :
  --------
  
  > data(engel)
  
  > plot(foodexp ~ income, data = engel, cex = 0.5, col = "blue", 
  +     xlab = "Household Income", ylab = "Food Expenditure")
  
  > z <- rq(foodexp ~ income, tau = 0.5, data = engel)
  
  > abline(z, col = "dark blue")
  
  > abline(lm(foodexp ~ income, data = engel), lty = 2, 
  +     col = "red")
  
  > taus <- c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
  
  > nt <- length(taus)
  
  > for (i in 1:length(taus)) {
  +     abline(rq(foodexp ~ income, tau = taus[i], data = engel), 
  +         col = "gray")
  + }
  
  > legend("bottomright", c("L1 (tau = .50)", "OLS", paste("tau= ", 
  +     formatC(rev(taus)))), col = c("dark blue", "red", rep("gray", 
  +     nt)), lt .... [TRUNCATED] 
  
  engel2.R :
  --------
  
  > data(engel)
  
  > x.poor <- quantile(engel[, "income"], 0.1)
  
  > x.rich <- quantile(engel[, "income"], 0.9)
  
  > z <- rq(foodexp ~ income, tau = -1, data = engel)
  
  > ps <- z$sol["tau", ]
  
  > coefs <- z$sol[4:5, ]
  
  > qs.poor <- c(c(1, x.poor) %*% coefs)
  
  > qs.rich <- c(c(1, x.rich) %*% coefs)
  
  > par(mfrow = c(1, 2))
  
  > plot(c(ps, ps), c(qs.poor, qs.rich), type = "n", xlab = expression(tau), 
  +     ylab = "quantile")
  
  > plot(stepfun(ps, c(qs.poor[1], qs.poor)), do.points = FALSE, 
  +     add = TRUE)
  
  > plot(stepfun(ps, c(qs.poor[1], qs.rich)), do.points = FALSE, 
  +     add = TRUE, col.hor = "gray", col.vert = "gray")
  
  > ps.wts <- (c(0, diff(ps)) + c(diff(ps), 0))/2
  
  > ap <- akj(qs.poor, z = qs.poor, p = ps.wts)
  
  > ar <- akj(qs.rich, z = qs.rich, p = ps.wts)
  
  > plot(c(qs.poor, qs.rich), c(ap$dens, ar$dens), type = "n", 
  +     xlab = "Food Expenditure", ylab = "Density")
  
  > lines(qs.rich, ar$dens, col = "gray")
  
  > lines(qs.poor, ap$dens, col = "black")
  
  > legend("topright", c("poor", "rich"), lty = c(1, 1), 
  +     col = c("black", "gray"))
  
  hinged.R :
  --------
  
  > require(quantreg)
  
  > if (requireNamespace("interp")) {
  +     do.rgl <- interactive() && require(rgl)
  +     n <- 1000
  +     x <- runif(n)
  +     y <- runif(n)
  +     z <- - .... [TRUNCATED] 
  NULL
  
  panelfig.R :
  ----------
  
  > library(quantreg)
  
  > rq.fit.panel <- function(X, y, s, w = c(1/3, 1/3, 
  +     1/3), taus = c(0.25, 0.5, 0.75), lambda = 0) {
  +     require(SparseM)
  +     require(quantre .... [TRUNCATED] 
  
  > n <- 3
  
  > T <- 50
  
  > nT <- n * T
  
  > u1 <- rnorm(T)
  
  > u2 <- rnorm(T)
  
  > u3 <- rnorm(T)
  
  > x1 <- rnorm(T, 1, 0.85)
  
  > x2 <- rnorm(T, 4, 1)
  
  > x3 <- rnorm(T, 7, 1)
  
  > beta1 <- 1
  
  > beta2 <- 0.5
  
  > y1 <- 0 + beta1 * x1 + (beta2 * x1) * u1
  
  > y2 <- 4 + beta1 * x2 + (beta2 * x2) * u2
  
  > y3 <- 8 + beta1 * x3 + (beta2 * x3) * u3
  
  > plot(c(0, 9), c(0, 25), type = "n", xlab = expression(x[it]), 
  +     ylab = expression(y[it]))
  
  > points(x1, y1, pch = 15)
  
  > points(x2, y2, pch = 15, col = "blue")
  
  > points(x3, y3, pch = 15, col = "red")
  
  > legend(1, 17, paste("i = ", 1:3, sep = ""), pch = 15, 
  +     col = c("black", "blue", "red"))
  
  > ya <- c(y1, y2, y3)
  
  > xa <- c(x1, x2, x3)
  
  > taus <- c(0.25, 0.5, 0.75)
  
  > xx <- seq(min(xa), max(xa), 0.25)
  
  > f <- coef(rq(ya ~ xa, tau = taus))
  
  > yy <- cbind(1, xx) %*% f
  
  > for (i in 1:3) lines(xx, yy[, i], col = "grey")
  
  > s <- rep(1:n, rep(T, n))
  
  > fp <- rq.fit.panel(xa, ya, s)$coef
  
  > bhat <- fp[1:3]
  
  > fehat <- fp[4:6]
  
  > xx1 <- seq(min(x1), max(x1), 0.25)
  
  > for (i in 1:3) {
  +     yy1 <- fehat[1] + bhat[i] * xx1
  +     lines(xx1, yy1, col = "black")
  + }
  
  > xx2 <- seq(min(x2), max(x2), 0.25)
  
  > for (i in 1:3) {
  +     yy2 <- fehat[2] + bhat[i] * xx2
  +     lines(xx2, yy2, col = "blue")
  + }
  
  > xx3 <- seq(min(x3), max(x3), 0.25)
  
  > for (i in 1:3) {
  +     yy3 <- fehat[3] + bhat[i] * xx3
  +     lines(xx3, yy3, col = "red")
  + }
  
  predemo.R :
  ---------
  
  > par(ask = TRUE)
  
  > sub <- "Classical Gaussian Version"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > ss = seq(0, 27, len = 100)
  
  > flm <- lm(dist ~ speed + I(speed^2), data = cars)
  
  > cflm <- predict(flm, list(speed = ss), data = cars, 
  +     interval = "confidence")
  
  > pflm <- predict(flm, list(speed = ss), data = cars, 
  +     interval = "prediction")
  
  > matlines(ss, pflm, lty = c(1, 2, 2), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, cflm, lty = c(1, 2, 2), col = c("black", 
  +     "blue", "blue"))
  
  > fl = rq(dist ~ speed + I(speed^2), data = cars, tau = 0.25)
  
  > fu = rq(dist ~ speed + I(speed^2), data = cars, tau = 0.75)
  
  > sub <- "Covariance matrix estimation using nid"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     se = "nid")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     se = "nid")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > sub <- "Covariance matrix estimation using xy-bootstrap"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     se = "boot", bsmethod = "xy")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     se = "boot", bsmethod = "xy")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > sub <- "Covariance matrix estimation using weighted xy-bootstrap"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     se = "boot", bsmethod = "wxy")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     se = "boot", bsmethod = "wxy")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > sub <- "Percentile method using xy-bootstrap"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     type = "percentile", se = "boot", bsmethod = "xy")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     type = "percentile", se = "boot", bsmethod = "xy")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > sub <- "Direct method using iid covariance matrix"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     type = "direct", se = "iid")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     type = "direct", se = "iid")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > sub <- "Direct method using nid covariance matrix"
  
  > plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
  +     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
  +     ylab = "stopp ..." ... [TRUNCATED] 
  
  > pl <- predict(fl, list(speed = ss), interval = "confidence", 
  +     type = "direct", se = "nid")
  
  > pu <- predict(fu, list(speed = ss), interval = "confidence", 
  +     type = "direct", se = "nid")
  
  > matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
  +     "red", "red"))
  
  > matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
  +     "red", "red"))
  
  > par(ask = FALSE)
  
  rqsslasso.R :
  -----------
  
  > n <- 100
  
  > p <- 9
  
  > q <- 3
  
  > beta <- c(rep(1, q), rep(0, p - q))
  
  > w <- matrix(rnorm(n * p), n, p)
  
  > x <- runif(n, 0, 10)
  
  > z <- runif(n, 0, 10)
  
  > y <- w %*% beta + sin(x) + (z^2)/50 + rnorm(n)/5
  
  > d <- data.frame(w, x, y, z)
  
  > f <- rqss(y ~ w + qss(x, lambda = 3) + qss(z, lambda = 2), 
  +     method = "lasso", lambda = 3, data = d)
  
  > plot(f, bands = "both", bcol = c("lightsteelblue", 
  +     "lightsteelblue4"))
  
  stack.R :
  -------
  
  > require(quantreg)
  
  > data(stackloss)
  
  > logLik.rq.process <- function(fit) {
  +     y <- model.response(model.frame(fit))
  +     fhat <- predict(fit, type = "fhat")
  +     fy <- mapply(functi .... [TRUNCATED] 
  
  > f0 <- rq(stack.loss ~ 1, tau = -1)
  
  > f1 <- rq(stack.loss ~ stack.x, tau = -1)
  
  > l0 <- logLik(f0)
  
  > l1 <- logLik(f1)
  
  > f0 <- rq(stack.loss ~ 1, tau = 1:19/20)
  
  > f1 <- rq(stack.loss ~ stack.x, tau = 1:19/20)
  
  > l0 <- logLik(f0)
  
  > l1 <- logLik(f1)
  
  subset.R :
  --------
  
  > require(quantreg)
  
  > n <- 200
  
  > x <- sort(rchisq(n, 4))
  
  > z <- rnorm(n)
  
  > s <- sample(1:n, n/2)
  
  > y <- log(x) + rnorm(n)/5
  
  > D = data.frame(y = y, x = x, z = z, s = (1:n) %in% 
  +     s)
  
  > plot(x, y, cex = 0.5, col = "grey")
  
  > points(x[s], y[s], col = "pink", cex = 0.5)
  
  > lam = 0.2
  
  > f0 <- rqss(y ~ qss(x, lambda = lam) + z, subset = s)
  
  > f1 <- rqss(y ~ qss(x, lambda = lam) + z, subset = s, 
  +     data = D)
  
  > plot(f0, add = TRUE, col = 2, lwd = 3)
  
  > plot(f1, add = TRUE, col = 4, lwd = 3)
  Warning messages:
  1: In summary.rq(x, se = se, R = R, covariance = TRUE) :
    5 non-positive fis
  2: In survfit.coxph(x) :
    the model contains interactions; the default curve based on columm means of the X matrix is almost certainly not useful. Consider adding a newdata argument.
  3: In pt(abs(ptab[, 3]), n - edf) : NaNs produced
  4: In pf(ntab[i - 1, 4], ntab[i - 1, 1], n - edf) : NaNs produced
  5: In summary.rq(object, cov = TRUE, ...) : 2 non-positive fis
  6: In summary.rq(object, cov = TRUE, ...) : 2 non-positive fis
  > 
  > 
  > proc.time()
     user  system elapsed 
    31.48    3.57   35.07