带置信区间的画图问题

仔细把底下三个程序看完
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")

        

epoh前辈，谢谢您了
您真是太厉害了，这应该和R内部的程序一样了
不过您还是使用rq函数，而不是在rq.fit.panel运行基础上的
谢谢

不太懂你的题意