epoh前辈，谢谢您
可以看出，您很认真负责

epoh前辈，我运行这一段时，会提示错误
t=1.96
for (i in 1:ntaus){
#coefficient1
coef1mat[i,1]=taus
coef1mat[i,2]=coef[(2*(i-1)+1),1]
coef1mat[i,3]=coef[(2*(i-1)+1),2]
coef1mat[i,4]=coef[(2*(i-1)+1),1]+t*coef[(2*(i-1)+1),2]
coef1mat[i,5]=coef[(2*(i-1)+1),1]-t*coef[(2*(i-1)+1),2]
#coefficient2
coef2mat[i,1]=taus
coef2mat[i,2]=coef[2*i,1]
coef2mat[i,3]=coef[2*i,2]
coef2mat[i,4]=coef[2*i,1]+t*coef[2*i,2]
coef2mat[i,5]=coef[2*i,1]-t*coef[2*i,2]
}
coef1mat
coef2mat

错误于coef1mat[i, 1] = taus : 被替换的项目不是替换值长度的倍数

我哪里出错？

coef1mat[i,1]=taus[i]
coef2mat[i,1]=taus[i]

额，epoh老师，已经可以了
谢谢您的指导

epoh老师，又碰到疑问了
我看到一些文章运行出各分位点的系数后，如在0.1,0.2,0.3,0.4,0.5的系数值分别为coeff1，coeff2，coeff3，coeff4，coeff5，文章提到了通过统计量检验这些系数间是否存在显著差异。有些文章提到了系数的联合显著性检验（ie，有些文章提到：对FDI在不同分位点的各系数差异性的联合F检验），这些在R有相关命令吗？
附件中有一篇文章，其中我用红色标志了2个地方，您看在R中有相关命令可以实现吗？

function anova.rq()

There are two (as yet) distinct forms of the test.
In the first the fitted objects all have the same specified quantile (tau) and
the intent is to test the hypothesis that smaller models are adequate
relative to the largest specified model.

In the second form of the test the linear predictor of the fits
are all the same, but the specified quantiles (taus) are different

epoh前辈，anova.rq函数是和rq函数一起使用的
ie：
data(barro)
fit0 <- rq(y.net ~  lgdp2 + fse2 + gedy2 , data = barro)
fit1 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro)
fit2 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.75)
fit3 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.25)
anova(fit1,fit0)
anova(fit1,fit2,fit3)
anova(fit1,fit2,fit3,joint=FALSE)

我要的是anova(fit1,fit2,fit3,joint=FALSE)这个实现结果。我看了anova.rq内部函数，呵呵，好复杂啊
内部程序：
> anova.rq
function (object, ...)
{
    if (length(list(object, ...)) > 1) {
        return(anova.rqlist(object, ...))
    }
    stop("Anova is only defined (yet) for lists of rq objects")
}
> anova.rqlist
function (object, ..., test = "Wald", joint = TRUE, score = "wilcoxon",
    R = 200, trim = NULL)
{
    objects <- list(object, ...)
    responses <- as.character(lapply(objects, function(x) formula(x)[[2]]))
    sameresp <- responses == responses[1]
    if (!all(sameresp))
        stop("Models don't all have the same response variable")
    n <- length(objects[[1]]$y)
    models <- as.character(lapply(objects, function(x) formula(x)))
    nobjects <- length(objects)
    dimp <- lapply(objects, function(x) length(coef(x)))
    objects <- objects[order(-unlist(dimp))]
    mf <- model.frame(objects[[1]])
    models <- as.character(lapply(objects, function(x) formula(x)))
    taus <- unlist(lapply(objects, function(x) x$tau))
    if (is.matrix(coef(objects[[1]])))
        names <- lapply(objects, function(x) dimnames(coef(x))[[1]])
    else names <- lapply(objects, function(x) names(coef(x)))
    if (test == "Wald")
        objects <- lapply(objects, function(x) summary(x, se = "nid",
            covariance = TRUE))
    sametaus <- taus == taus[[1]]
    if (all(sametaus)) {
        Tn <- rep(0, nobjects - 1)
        ndf <- Tn
        ddf <- Tn
        pvalue <- Tn
        topnote <- paste("Model ", format(1:nobjects), ": ",
            models, sep = "", collapse = "\n")
        if (test == "anowar") {
            X1 <- model.matrix(objects[[1]], mf, contrasts = objects[[1]]$contrasts)
            y <- model.response(mf)
            weights <- model.weights(mf)
            tau <- taus[[1]]
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]]))
                  stop("Models aren't nested")
                mf <- model.frame(objects[[i]])
                X0 <- model.matrix(objects[[i]], mf, contrasts = objects[[i]]$contrasts)
                Htest <- rq.test.anowar(X0, X1, y, tau = tau,
                  R = R)
                ndf[i - 1] <- Htest$ndf
                Tn[i - 1] <- Htest$Tn
                ddf[i - 1] <- Htest$ddf
                pvalue[i - 1] <- Htest$pvalue
            }
            table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else if (test == "rank") {
            x1 <- model.matrix(objects[[1]], mf, contrasts = objects[[1]]$contrasts)
            y <- model.response(mf)
            weights <- model.weights(mf)
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]]))
                  stop("Models aren't nested")
                nullH <- is.na(match(names[[1]], names[[i]]))
                X1 <- as.matrix(x1[, nullH])
                mf <- model.frame(objects[[i]])
                X0 <- model.matrix(objects[[i]], mf, contrasts = objects[[i]]$contrasts)
                if (score == "tau")
                  tau <- taus[[1]]
                Htest <- rq.test.rank(X0, X1, y, score = score,
                  weights = weights, tau = tau, trim = trim)
                ndf[i - 1] <- Htest$ndf
                Tn[i - 1] <- Htest$Tn
                ddf[i - 1] <- Htest$ddf
                pvalue[i - 1] <- Htest$pvalue
            }
            table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else if (test == "Wald") {
            V <- lapply(objects, function(x) x$cov)
            coef <- lapply(objects, function(x) coef(x)[, 1])
            for (i in 2:nobjects) {
                if (!all(names[[i]] %in% names[[1]]))
                  stop("Models aren't nested")
                nullH <- is.na(match(names[[1]], names[[i]]))
                ndf[i - 1] <- sum(nullH)
                Tn[i - 1] <- t((coef[[1]])[nullH]) %*% solve((V[[1]])[nullH,
                  nullH], (coef[[1]])[nullH])/ndf[i - 1]
                ddf[i - 1] <- n - length(names[[1]])
                pvalue[i - 1] <- 1 - pf(Tn[i - 1], ndf[i - 1],
                  ddf[i - 1])
            }
            table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else stop("test only defined for anowar, Wald and rank")
    }
    else {
        m <- length(taus)
        for (i in 2:m) {
            if (!setequal(names[[i]], names[[1]]))
                stop("Models with common tau don't have same X")
        }
        if (names[[1]][1] != "(Intercept)")
            stop("Intercept required in common tau testing")
        Omega <- outer(taus, taus, pmin) - outer(taus, taus)
        J <- objects[[1]]$J
        p <- dim(J)[1]
        H <- array(unlist(lapply(objects, function(x) x$Hinv)),
            c(p, p, m))
        H <- matrix(aperm(H, c(1, 3, 2)), p * m, p) %*% t(chol(J))
        W <- (H %*% t(H)) * (kronecker(Omega, outer(rep(1, p),
            rep(1, p))))
        coef <- unlist(lapply(objects, function(x) coef(x)[,
            1]))
        if (joint) {
            D <- kronecker(diff(diag(m)), cbind(0, diag(p - 1)))
            ndf <- (p - 1) * (m - 1)
            Tn <- t(D %*% coef) %*% solve(D %*% W %*% t(D), D %*%
                coef)/ndf
            ddf <- n * m - (p - 1) * (m - 1)
            pvalue <- 1 - pf(Tn, ndf, ddf)
            nobjects <- 1
            tnote1 <- paste("Model: ", models[[1]], "\n", sep = "")
            tnote2 <- paste("Joint Test of Equality of Slopes: tau in { ",
                paste(taus, collapse = " "), " }\n")
            topnote <- paste(tnote1, tnote2, sep = "")
            table <- data.frame(ndf, ddf, Tn, pvalue)
        }
        else {
            Tn <- pvalue <- rep(0, p - 1)
            ndf <- m - 1
            ddf <- n * m - (m - 1)
            for (i in 2:p) {
                E <- matrix(0, 1, p)
                E[1, i] <- 1
                D <- kronecker(diff(diag(m)), E)
                Tn[i - 1] <- t(D %*% coef) %*% solve(D %*% W %*%
                  t(D), D %*% coef)/ndf
                pvalue[i - 1] <- 1 - pf(Tn[i - 1], ndf, ddf)
            }
            tnote1 <- paste("Model: ", models[[1]], "\n", sep = "")
            tnote2 <- paste("Tests of Equality of Distinct Slopes: tau in { ",
                paste(taus, collapse = " "), " }\n")
            topnote <- paste(tnote1, tnote2, sep = "")
            table <- data.frame(ndf, ddf, Tn, pvalue)
            dimnames(table)[[1]] <- names[[1]][2:p]
        }
    }
    x <- list(table = table, topnote = topnote)
    class(x) <- "anova.rq"
    return(x)
}

呵呵，看来rq.sfn.fit这个函数还是存在很大缺陷，主要还是用rq（）函数

哈哈!老兄,
应该说package 更新的速度,
不及你研究的"深度".
我稍瞄一下:
The standard errors (in parenthesis) are obtained after 1000 panel-bootstrap replications
我抽空,看看能否有解.
你有没这组DATA?

epoh前辈，见笑了，呵呵。我只是对方法很感兴趣。
其实这些估计分位的函数，包括quantreg包都是roger教授编写的。他专门写了一个rq.fit.panel（不属于R内部函数），其实就是调用rq.fit.sfn函数来估计面板数据分位，而据他主页上介绍rqss（）函数和rq（）函数也是可以实现这一功能了。不过前段时间您已经帮忙看了，他只能一次估计一个分位点，而不能同时估计多个分位点。

这篇文章是roger教授的学生写的一篇，这位老师一般都不共享数据的，所以我也没有数据。
额，标准误是通过panel-bootstrap产生的，这个您已经实现了。

哈哈，可惜epoh老师不是专门学习quantile regression的，要不我就可以向您请教疑问了。
我现在还堆积着很多问题没有得到解答，看的深度不够，基础不扎实啊。

前辈，您好，在论坛里看到您是R语言的高手，我问一下，为什么欲行面板分位数回归时总是出问题：
> library(SparseM)
> library(quantreg)
> library(Ecdat)
> data(Grunfeld)
> #panel data Grunfeld
> #n=10, T=20, N=200
> x1=Grunfeld$value
> x2=Grunfeld$capital
> X=cbind(x1,x2)
> y=Grunfeld$inv
> s= rep(1:10,rep(20,10))
> lambda=1
> w=c(0.1,0.2,0.3,0.4,0.5,0.4,0.3,0.2,0.1)
> taus=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
> ntaus=length(taus)
> K <- length(w)
> X <- as.matrix(X)
> p <- ncol(X)
> n <- length(levels(as.factor(s)))
> N <- length(y)
> Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
> Fidelity <- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)  
错误于cbind(deparse.level, ...) : args have differing numbers of rows
此外: 警告信息：
In as.matrix.csr(x, nrow = nrow, ncol = ncol, eps = eps) :
  ncol*nrow indivisable by length(x)
> Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))  
> D <- rbind(Fidelity,Penalty)   
错误于rbind(Fidelity, Penalty) : 找不到对象'Fidelity'
> yy <- 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,n))  
> rq.fit.sfn(D,yy,tau=taus,rhs=a)
请您指教
我计算过，一个是200*18的矩阵，一个是200*90的矩阵，为什么会提示说行不同呢？？请您指点一下

嘿嘿,一年了,有点忘了
不知ywh兄,还记不记得?

呵呵，刚才仔细看了一下， 是W的定义有问题

呵呵，epoh老师，楼主自己解决了。

DB file 是MicroTSP 的数据文件，在TSP中用fetch得到。
epoh老师，DB结尾的数据文件在TSP中如何打开呢？
好像是一个网页数据。

MicroTSP Databank Format
  http://econpy.googlecode.com/svn/trunk/pytrix/opendatabank.txt

?A simple example
?x.db  (x series:1,3,5,7,9)
?Undated Series
?Line 1 : starting index
?Line 2 : ending index
1
5
1
3
5
7
9
? save as x.db
??????
??????
?open x.db
fetch x;
print x;
? save as x.tsp

呵呵，谢谢epoh老师，可以显示
我开始用了dbprint命令，一直显示错误。

epoh 老师，我运行上边程序出现错误：

> xx1<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
错误于cbind(deparse.level, ...) : args have differing numbers of rows
此外: 警告信息：
In as.matrix.csr(x, nrow = nrow, ncol = ncol, eps = eps) :
  ncol*nrow indivisable by length(x)

哈哈，epoh老师太忙了
你可以直接利用rqpd 包了。

这个程序之前你也用过，那请问怎么回事？是哪里出问题了，你能修改过来吗？

你好，我发现这一段程序总是运行不过去，你之前运行过没有错误吗？为什么我的总出现以下错误？
急于求解：
> xx1<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
错误于cbind(deparse.level, ...) : args have differing numbers of rows
此外: 警告信息：
In as.matrix.csr(x, nrow = nrow, ncol = ncol, eps = eps) :
  ncol*nrow indivisable by length(x)

你的最后目的是什么？运行面板分位回归模型吗？
如果答案是，那你直接利用rqpd包就可以了，作者
已经把这些程序打包为package了，可以直接用

你是为了运行面板分位回归吗？
如果是，那用rqpd包就可以了

我不是很清楚这个包rqpd怎么用。。主要是我还实现惩罚项和工具变量

还有我不清楚为什么我运行那一段程序时总汇报那个错误，你当时没遇到吗？怎么解决的？

library(SparseM)
library(quantreg)
library(Ecdat)
data("Produc")  
x1=log(Produc$pcap)
x2=log(Produc$pc)
x3=log(Produc$emp)
x4=Produc$unemp
X=cbind(x1,x2,x3,x4)
y=log(Produc$gsp)
N<-length(y)
s= rep(1:48,rep(17,48))
X <- as.matrix(X)
p <- ncol(X)
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
lambda=1
w=c(0.1,0.25,0.5,0.25,0.1)
taus=c(0.1,0.25,0.5,0.75,0.9)
K <- length(w)
xx1<- cbind(as(w,"matrix.diag.csr") %x% X,as.matrix(w) %x% Z)
Penalty <- cbind(as.matrix.csr(0,n,K*p),lambda*as(n,"matrix.diag.csr"))
xx<-rbind(xx1,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(xx,y,a)