function rq.fit.panel()
is now superseded by the use of the function rqss() in the R quantreg package,
or the use of method = "lasso" for the function rq() in that package.

Package‘quantreg’Version 4.57  Date 2011-2-26  

#############
library(quantreg)
data(engel)
attach(engel)
summary(rq(foodexp~income,tau = 1:49/50,data=engel))

Call: rq(formula = foodexp ~ income, tau = 1:49/50, data = engel)
tau: [1] 0.02
Coefficients:
            coefficients lower bd  upper bd
(Intercept) 108.60499    103.39624 141.26345
income        0.34664      0.28526   0.34664

Call: rq(formula = foodexp ~ income, tau = 1:49/50, data = engel)
tau: [1] 0.04
............

#Example of plotting of coefficients and their confidence bands
plot(summary(rq(foodexp~income,tau = 1:49/50,data=engel)))

##epoh ,哈哈，谢谢啊
我看过rqss（）函数帮助，但是找不到识别面板结构的程序，还不知道怎么使用。

panel data:
  package "Ecdat"
  Produc Us States Production
  a panel of 48 observations from 1970 to 1986
  number of observations : 816

##############
library(Ecdat)
data("Produc")  
#n=48, T=17, N=816
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)
s= rep(1:48,rep(17,48))
w=1

####case 1
####use function rq.fit.panel()
require(SparseM)
require(quantreg)

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))
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
rq.fit.sfn(xx,y,tau=0.5)
>$coef
[1] -0.001856897  0.227955745  0.806905899 -0.003254027  2.348518722
[6]  2.531511125  2.416402785  2.646310826  2.554288408  2.598468250
[11]  2.540988532  2.473361510  2.401692284  2.497764134  2.537963266
[16]  2.425204261  2.482897448  2.510841856  2.553142722  2.703008416
[21]  2.404225489  2.537962766  2.486090983  2.567535469  2.455404207
[26]  2.403301315  2.463135907  2.545757570  2.474922409  2.515871136
[31]  2.443668078  2.580372182  2.643673666  2.605560174  2.384878447
[36]  2.541337214  2.468143678  2.587964518  2.509402455  2.432777812
[41]  2.495204924  2.258353874  2.465798076  2.372996183  2.600392620
[46]  2.471721790  2.457974490  2.519673773  2.590197854  2.472063213
[51]  2.470425047  2.862672642


####case 2
####use function rqss()
tt=as.matrix(xx)
fit=rqss(y~tt-1)
summary(fit)
>
Formula:
y ~ tt - 1

Parametric coefficients:
      Estimate    Std. Error   t value    Pr(>|t|)   
tt1  -0.001857   0.040474  -0.046   0.9634   
tt2   0.227956   0.033949   6.715 3.68e-11 ***
tt3   0.806906   0.041542  19.424  < 2e-16 ***
tt4  -0.003254   0.001402  -2.321   0.0205 *  
tt5   2.348519   0.227838  10.308  < 2e-16 ***
tt6   2.531511   0.226443  11.179  < 2e-16 ***
tt7   2.416403   0.216806  11.145  < 2e-16 ***
tt8   2.646311   0.254619  10.393  < 2e-16 ***
tt9   2.554288   0.222018  11.505  < 2e-16 ***
tt10  2.598468   0.220106  11.806  < 2e-16 ***
tt11  2.540989   0.210428  12.075  < 2e-16 ***
tt12  2.473362   0.234957  10.527  < 2e-16 ***
tt13  2.401692   0.226512  10.603  < 2e-16 ***
tt14  2.497764   0.208589  11.975  < 2e-16 ***
tt15  2.537963   0.241893  10.492  < 2e-16 ***
tt16  2.425204   0.227976  10.638  < 2e-16 ***
tt17  2.482897   0.232136  10.696  < 2e-16 ***
tt18  2.510842   0.231957  10.825  < 2e-16 ***
tt19  2.553143   0.231057  11.050  < 2e-16 ***
tt20  2.703008   0.244407  11.059  < 2e-16 ***
tt21  2.404225   0.201550  11.929  < 2e-16 ***
tt22  2.537963   0.228998  11.083  < 2e-16 ***
tt23  2.486091   0.220021  11.299  < 2e-16 ***
tt24  2.567535   0.239121  10.737  < 2e-16 ***
tt25  2.455404   0.235207  10.439  < 2e-16 ***
tt26  2.403301   0.223782  10.739  < 2e-16 ***
tt27  2.463136   0.226028  10.897  < 2e-16 ***
tt28  2.545758   0.227966  11.167  < 2e-16 ***
tt29  2.474922   0.232918  10.626  < 2e-16 ***
tt30  2.515871   0.211395  11.901  < 2e-16 ***
tt31  2.443668   0.201243  12.143  < 2e-16 ***
tt32  2.580372   0.224698  11.484  < 2e-16 ***
tt33  2.643674   0.222489  11.882  < 2e-16 ***
tt34  2.605560   0.253096  10.295  < 2e-16 ***
tt35  2.384878   0.221027  10.790  < 2e-16 ***
tt36  2.541337   0.225034  11.293  < 2e-16 ***
tt37  2.468144   0.238798  10.336  < 2e-16 ***
tt38  2.587965   0.227848  11.358  < 2e-16 ***
tt39  2.509402   0.221705  11.319  < 2e-16 ***
tt40  2.432778   0.239579  10.154  < 2e-16 ***
tt41  2.495205   0.191899  13.003  < 2e-16 ***
tt42  2.258354   0.211827  10.661  < 2e-16 ***
tt43  2.465798   0.222536  11.080  < 2e-16 ***
tt44  2.372996   0.232450  10.209  < 2e-16 ***
tt45  2.600393   0.251702  10.331  < 2e-16 ***
tt46  2.471722   0.214506  11.523  < 2e-16 ***
tt47  2.457974   0.199837  12.300  < 2e-16 ***
tt48  2.519674   0.226890  11.105  < 2e-16 ***
tt49  2.590198   0.241976  10.704  < 2e-16 ***
tt50  2.472063   0.224931  10.990  < 2e-16 ***
tt51  2.470425   0.229340  10.772  < 2e-16 ***
tt52  2.862673   0.230734  12.407  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


  Quantile Fidelity at tau = 0.5  is      10.8697
  Effective Degrees of Freedom = 52        Sample Size = 816

#（1）先回归
   lmfit <- lm( dat$Sunday ~ dat$Daily, data=dat )
# （2）  计算置信区间（95% ）---------------------------
     Ypred <- predict( lmfit, interval = "confidence", newdata= dat, level = 0.95 )
#（3） 画图
        plot(dat$Daily, dat$Sunday,
              xlab="Daily", ylab="Sunday",
             main="Figure:  Yhat Predict confidence")
        lines ( Ypred[ ,1 ] ~ dat$Daily, col="red", lwd=1 )
        lines ( Ypred[ ,2 ] ~ dat$Daily, col="gray", lwd=1 )
        lines ( Ypred[ ,3 ] ~ dat$Daily, col="gray", lwd=1 )

epoh，谢谢您了
呵呵，感谢楼上兄弟

收藏了。。。。。。。。。。。。。。

谢谢各位了，谢谢epoh老师

epoh老师，您好，能否帮我看看这个该如何实现？
您编写的case1，我修改了一下可以实现多个分位同时估计
case1：
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,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))
##rhs = (1 - tau) * c(t(a) %*% rep(1, length(y)))
rq.fit.sfn(xx,y,a)
case2： 我修改了一点
library(SparseM)
library(quantreg)
library(Ecdat)
data("Produc")  
#n=48, T=17, N=816
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)
s= rep(1:48,rep(17,48))
X <- as.matrix(X)
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
w<-1
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)  tt<-as.matrix(xx)
taus<-c(0.1,0.25,0.75,0.5,0.9)
for(j in 1:length(taus)){
fit<-rqss(y~tt-1,tau=taus[j])
}
但是case不能实现多个分位同时运行
就假如我一次一次运行，如现在运行5个分位点，就5次，每一次都估计出x1-x4系数，而每一次的固定效应结果不同。
在rq.fit.sfn中，我同时估计5个分位点，得出的只有一组固定效应。那两者肯定不一样
在rqss中能实现5个分位同时估计吗？

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,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))
##############case1
rq.fit.sfn(xx,y)
################
$coef
[1]  2.836199e-01  4.602578e-01  4.218402e-01 -1.293751e-02  2.836199e-01
[6]  4.602578e-01  4.218402e-01 -1.293751e-02  2.836199e-01  4.602578e-01
[11]  4.218402e-01 -1.293751e-02  2.836199e-01  4.602578e-01  4.218402e-01
[16] -1.293751e-02  2.836199e-01  4.602578e-01  4.218402e-01 -1.293751e-02
[21] -1.412773e-01  6.149486e-02  2.992886e-02 -4.000146e-02  1.343490e-01
[26]  2.414043e-01  2.334016e-01 -4.028156e-02 -4.768962e-02  1.829905e-01
[31] -5.174577e-02 -6.450664e-02 -5.472290e-02 -4.615042e-02  3.985560e-02
[36]  2.038454e-10  1.901025e-01  7.812054e-02  1.347637e-01  2.923434e-11
[41] -8.490129e-02 -3.631191e-02  2.705989e-02  1.312208e-03 -1.009822e-01
[46]  1.433116e-01  2.792632e-01  1.822704e-01  1.587613e-01 -5.398773e-02
[51]  3.816106e-10  1.061240e-02 -9.002397e-02  7.444380e-02  9.299111e-02
[56] -1.201060e-01  4.259159e-01 -5.967792e-02 -4.599368e-03 -1.326447e-01
[61] -1.550668e-01  1.112986e-01  2.706360e-01  9.057896e-02 -1.589800e-02
[66] -2.147877e-02 -3.657946e-03  2.225002e-01


###############case2
tt<-as.matrix(xx)
fit<-rqss(y~tt-1)
summary(fit)
Parametric coefficients:
       Estimate Std. Error t value Pr(>|t|)   
tt1   2.836e-01  1.726e-02  16.433  < 2e-16 ***
tt2   4.603e-01  1.473e-02  31.244  < 2e-16 ***
tt3   4.218e-01  1.098e-02  38.425  < 2e-16 ***
tt4  -1.294e-02  1.002e-03 -12.908  < 2e-16 ***
tt5   2.836e-01  1.686e-02  16.827  < 2e-16 ***
tt6   4.603e-01  1.441e-02  31.933  < 2e-16 ***
tt7   4.218e-01  1.077e-02  39.172  < 2e-16 ***
tt8  -1.294e-02  9.440e-04 -13.705  < 2e-16 ***
tt9   2.836e-01  1.629e-02  17.411  < 2e-16 ***
tt10  4.603e-01  1.413e-02  32.565  < 2e-16 ***
tt11  4.218e-01  1.048e-02  40.260  < 2e-16 ***
tt12 -1.294e-02  1.007e-03 -12.850  < 2e-16 ***
tt13  2.836e-01  1.686e-02  16.827  < 2e-16 ***
tt14  4.603e-01  1.441e-02  31.933  < 2e-16 ***
tt15  4.218e-01  1.077e-02  39.172  < 2e-16 ***
tt16 -1.294e-02  9.440e-04 -13.705  < 2e-16 ***
tt17  2.836e-01  1.726e-02  16.433  < 2e-16 ***
tt18  4.603e-01  1.473e-02  31.244  < 2e-16 ***
tt19  4.218e-01  1.098e-02  38.425  < 2e-16 ***
tt20 -1.294e-02  1.002e-03 -12.908  < 2e-16 ***
tt21 -1.413e-01  8.672e-03 -16.291  < 2e-16 ***
tt22  6.149e-02  9.322e-03   6.596 4.78e-11 ***
tt23  2.993e-02  8.580e-03   3.488 0.000492 ***
tt24 -4.000e-02  1.679e-02  -2.382 0.017265 *  
tt25  1.343e-01  9.265e-03  14.501  < 2e-16 ***
tt26  2.414e-01  1.648e-02  14.649  < 2e-16 ***
tt27  2.334e-01  1.483e-02  15.738  < 2e-16 ***
tt28 -4.028e-02  1.087e-02  -3.707 0.000213 ***
tt29 -4.769e-02  9.713e-03  -4.910 9.49e-07 ***
tt30  1.830e-01  9.089e-03  20.133  < 2e-16 ***
tt31 -5.175e-02  1.210e-02  -4.278 1.94e-05 ***
tt32 -6.451e-02  8.223e-03  -7.845 5.55e-15 ***
tt33 -5.472e-02  1.170e-02  -4.679 2.98e-06 ***
tt34 -4.615e-02  1.009e-02  -4.576 4.90e-06 ***
tt35  3.986e-02  1.140e-02   3.497 0.000477 ***
tt36  2.038e-10  4.322e-03   0.000 1.000000   
tt37  1.901e-01  9.868e-03  19.264  < 2e-16 ***
tt38  7.812e-02  1.591e-02   4.910 9.50e-07 ***
tt39  1.348e-01  1.544e-02   8.729  < 2e-16 ***
tt40  2.923e-11  6.683e-03   0.000 1.000000   
tt41 -8.490e-02  1.161e-02  -7.316 3.09e-13 ***
tt42 -3.631e-02  1.556e-02  -2.333 0.019698 *  
tt43  2.706e-02  1.056e-02   2.563 0.010407 *  
tt44  1.312e-03  7.382e-03   0.178 0.858914   
tt45 -1.010e-01  1.695e-02  -5.956 2.82e-09 ***
tt46  1.433e-01  1.238e-02  11.579  < 2e-16 ***
tt47  2.793e-01  1.128e-02  24.764  < 2e-16 ***
tt48  1.823e-01  1.111e-02  16.407  < 2e-16 ***
tt49  1.588e-01  1.160e-02  13.682  < 2e-16 ***
tt50 -5.399e-02  1.946e-02  -2.774 0.005568 **
tt51  3.816e-10  5.112e-03   0.000 1.000000   
tt52  1.061e-02  1.576e-02   0.674 0.500640   
tt53 -9.002e-02  1.151e-02  -7.820 6.66e-15 ***
tt54  7.444e-02  9.354e-03   7.958 2.22e-15 ***
tt55  9.299e-02  1.008e-02   9.227  < 2e-16 ***
tt56 -1.201e-01  1.349e-02  -8.906  < 2e-16 ***
tt57  4.259e-01  1.282e-02  33.218  < 2e-16 ***
tt58 -5.968e-02  9.418e-03  -6.337 2.61e-10 ***
tt59 -4.599e-03  1.166e-02  -0.395 0.693160   
tt60 -1.326e-01  1.287e-02 -10.304  < 2e-16 ***
tt61 -1.551e-01  9.147e-03 -16.952  < 2e-16 ***
tt62  1.113e-01  1.079e-02  10.312  < 2e-16 ***
tt63  2.706e-01  1.195e-02  22.653  < 2e-16 ***
tt64  9.058e-02  1.167e-02   7.759 1.09e-14 ***
tt65 -1.590e-02  1.465e-02  -1.085 0.277769   
tt66 -2.148e-02  1.096e-02  -1.959 0.050195 .  
tt67 -3.658e-03  1.035e-02  -0.354 0.723707   
tt68  2.225e-01  1.932e-02  11.514  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


  Quantile Fidelity at tau = 0.5  is      17.5209
  Effective Degrees of Freedom = 248       Sample Size = 4128

##epoh，我对您程序还要疑惑，为什么在case1中，函数rq.fit.sfn(xx,y)不加入
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内部是rq.fit.sfn（D，x，a）
在您程序中，为什么没有指明是这5个分位点，而就直接可以？

我试加了，程序也是可以运行的，只是结果不一样
> 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,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)
$coef
[1]  0.3479225632  0.4579939079  0.3419116807 -0.0132019291  0.3479225632
[6]  0.4579939079  0.3419116807 -0.0132019291  0.3479225632  0.4579939079
[11]  0.3419116807 -0.0132019291  0.3479225632  0.4579939079  0.3419116807
[16] -0.0132019291  0.3479225632  0.4579939079  0.3419116807 -0.0132019291
[21] -0.1333094818 -0.0001001004  0.0005137676 -0.0406687377  0.0783307839
[26]  0.1897515687  0.1388376698 -0.0243480092 -0.0477144706  0.1091862206
[31] -0.0495600876 -0.0631116318 -0.0821633242 -0.0649543957 -0.0088250215
[36] -0.0295872475  0.1433398180  0.0385951674  0.1092909477 -0.0157361471
[41] -0.1005352958 -0.0707680805  0.0030014620 -0.0804312417 -0.1401517851
[46]  0.0825749965  0.2307316486  0.1470793361  0.0835593363 -0.0659401370
[51]  0.0033954589 -0.0757982221 -0.0809079723  0.0426541558  0.0304444640
[56] -0.1109061177  0.3313202668 -0.0641237261 -0.0821561387 -0.1284287320
[61] -0.1464990138  0.0452016648  0.1755758956  0.0560798608 -0.0777549839
[66] -0.0559941815 -0.0202300682  0.0801947593

您觉得我哪里理解错了？

糊涂了,封掉a没发觉.
这问题,我再仔细看一下function rqss().

谢谢epoh指导
rqss（）函数也有用到方法lasso，所以也是加了惩罚项，因此才可以实现rq.panel.fit功能（他就是通过引入固定效应惩罚项的）。
我看了很久，没有看懂

epoh老师，如果rqss（）可以实现多个分位同时估计，那也直接给出了标准误和系数，能不能直接画出第一帖附件中的图形？

谢谢您

仔细把底下三个程序看完
function quantreg()
  function rqss()
  function sfn()

  1.function rq()  method = "lasso"
     清楚的告诉你
     if(object$method == "lasso")
        stop("no inference for lasso'd rq fitting: try rqss (if brave)")

  2.看来panel data还是只能适用function rq.fit.sfn()
    a = structure of the design matrix A' stored in csr format
    毕竟是一个  "design matrix"--------s,taus,w

    不过只估计出coefficients & resid
    楼主需要的统计量,可能要用bootstrap
.

##epoh前辈，根据作者说明是rq和rqss都可以实现rq.fit.sfn功能的了
我都试过了
1、rq函数选择lasso时根本不可行，结果和您一样
2、rqss（）函数我觉得应该可以的，从您第一次用一个分位点试验的结果可以看出，他的结果和rq.fit.sfn结果是一样的，所以应该是可以的，我只是不知道怎么修改成同时估计多个分位点。
3、因为rq.fit.sfn函数要经过bootstrap抽样，感觉不是很稳定

前辈，若使用rq.fit.sfn函数，再bootstrap得出标准误和t值，可以计算其置信区间得到 帖子中的图像吗？

呵呵，那个bootstrap实现统计量的程序已经让您编写了，我只是觉得这个程序不稳定了

rqss（）函数只适用一个分位点
无法估计同时估计多个分位点
最主要是在程序中summary.rqss()
已经写明
if(tau + h > 1)
    stop("tau + h > 1:  error in summary.rq")

ex:tau=0.5
     h=0.1
      if(tau + h > 1)
      stop("tau + h > 1:  error in summary.rq")
      #OK
   
    tau=c(0.5,0.75)
      if(tau + h > 1)
      stop("tau + h > 1:  error in summary.rq")
      #PROBLEM
   
  function rqss()呼叫
  rqss.fit(X, Y, tau = tau, rhs = rhs, method = "sfnc")
  再呼叫
  fit <- sfn = rq.fit.sfn(x, y, tau = tau,  nsubmax = nsubmax, nnzlmax = nnzlmax, tmpmax = tmpmax, ...)
  再由fortran srqfnb.o <- .Fortran("srqfn",....)

  所产生的object fit 就是提供 coef &resid

  使用rq.fit.sfn函数，再bootstrap得出标准误和t值，
  有了标准误,就可以计算其置信区间,得到帖子中的图像.

epoh前辈，谢谢您了
画图程序我自己看看，再向您请教

rqss不能实现多个分位点，那就没有什么意义了，因为quantile regression要的就是对几个分位进行加权处理

Stata:
command qreg:
  qreg depvar indepvars , quantile(choice of quantile)

  The standard errors estimated by qreg are valid for the
  location model only. To obtain better standard errors
  estimates (given with the bootstrap), you should use
bsqreg instead.

  To obtain simultaneously several quantile regressions, use
  sqreg:
  sqreg depvar indepvars , quantiles(choice of quantiles)

  N.B: The standard errors estimates provided by sqreg
         Are computed with the bootstrap.

epoh前辈，谢谢您
呵呵，可能我刚才说的用bootstrap结果不稳定，是不正确的
那个stata命令我看过了，如果是用截面数据，那很完美，图形也非常漂亮

可惜stata目前还未有估计面板分位的命令


呵呵，前辈，您能否把求置信区间的命令也写出来？

Grunfeld   Grunfeld Investment Data
  a panel of 10 observations from 1935 to 1954
  number of observations : 200
##########
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))
w=1
n <- length(levels(as.factor(s)))
Z <- as.matrix.csr(model.matrix(~as.factor(s)-1))
xx<- cbind(as(w,"matrix.diag.csr") %x% X,w %x% Z)
####################
tt=as.matrix(xx)
taus=seq(0.1,0.9,0.01)
ntau=length(taus)  #81
coef1mat=matrix(NA,ntau,5)
coef2mat=matrix(NA,ntau,5)
rq.fit=rq(y~tt-1,taus)
sum=summary(rq.fit,"ker")
#If n is large(>1000), t=1.96
t=1.96
for (i in 1:ntau){
#coefficient1
coef1mat[i,1]=sum[]$tau
coef1mat[i,2]=sum[]$coefficients[1,1]
coef1mat[i,3]=sum[]$coefficients[1,2]
coef1mat[i,4]=sum[]$coefficients[1,1]+t*sum[]$coefficients[1,2]
coef1mat[i,5]=sum[]$coefficients[1,1]-t*sum[]$coefficients[1,2]
#coefficient2
coef2mat[i,1]=sum[]$tau
coef2mat[i,2]=sum[]$coefficients[2,1]
coef2mat[i,3]=sum[]$coefficients[2,2]
coef2mat[i,4]=sum[]$coefficients[2,1]+t*sum[]$coefficients[2,2]
coef2mat[i,5]=sum[]$coefficients[2,1]-t*sum[]$coefficients[2,2]
}
coef1mat
coef2mat
###plot coef1,coef2 95%CI
par(mfrow = c(1, 2))
####coefficient1
upperbd1=coef1mat[,4]
lowerbd1=coef1mat[,5]
coef1=coef1mat[,2]
plot(taus,upperbd1,type = 'n',  ylim = c(0.00, 0.20),
     ylab = 'Coef1 Value', xlab = 'Quantile')
lines(taus, lowerbd1, col = 'grey')
lines(taus, upperbd1, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd1, rev(lowerbd1)),
     col = "lightblue", border = 'black')
points(taus, coef1,col='blue',pch=21)
lines(taus, coef1,col='blue')
title("gross Investment vs value of the firm")
####coefficient2
upperbd2=coef2mat[,4]
lowerbd2=coef2mat[,5]
coef2=coef2mat[,2]
plot(taus,upperbd2,type = 'n',  ylim = c(0.00, 0.50),
     ylab = 'Coef2 Capital', xlab = 'Quantile')
lines(taus, lowerbd2, col = 'grey')
lines(taus, upperbd2, col = 'grey')
polygon(c(taus, rev(taus)), c(upperbd2, rev(lowerbd2)),
     col = "lightblue", border = 'black')
points(taus, coef2,col='blue',pch=21)
lines(taus, coef2,col='blue')
title("gross Investment vs stock of plant and equipment")