[R] Quantile regression: Discrepencies Between optimizer and rq()

[Roger Koenker]

Optim()  by default is using Nelder-Mead  which is an extremely poor way to
do linear programming, despite the fact that ?optim says that:  "It will work reasonably well for
non-differentiable functions."    I didn't check your coding of the objective function fully, but at the
very least you should explicitly pass the arguments y, x, and quant.  and you need to replace
what you call sample.cond.quantile by 0 in the definition of w.less.  


url:    www.econ.uiuc.edu/~roger            Roger Koenker
email    rkoenker at uiuc.edu            Department of Economics
vox:     217-333-4558                University of Illinois
fax:       217-244-6678                Urbana, IL 61801

On Jun 7, 2012, at 1:49 PM, Kevin Chang wrote:

> Hello Everyone,
> 
> 
> 
> I'm currently learning about quantile regressions. I've been using an
> optimizer to compare with the rq() command for quantile regression. 
> 
> When I run the code, the results show that my coefficients are consistent
> with rq(), but the intercept term can vary by a lot.
> 
> I don't think my optimizer code is wrong and suspects it has something to do
> with the starting values.
> 
> The results seems very sensitive to different starting values and I don't
> know how to make sense of it.
> 
> 
> 
> Advice from the community would be greatly appreciated.
> 
> 
> 
> Sincerely,
> 
> 
> Kevin Chang
> 
> 
> 
> ###################### CODE Below ###########################
> 
> 
> 
> library(quantreg)
> 
> data(engel)
> 
> y<-cbind(engel[,2])
> 
> x<-cbind(rep(1,length(engel[,1])),engel[,1])
> 
> x1<-cbind(engel[,1])
> 
> nn<-nrow(engel)
> 
> nn
> 
> 
> 
> bhat.ls<-solve(t(x)%*%x)%*%t(x)%*%y
> 
> #bhat.ls
> 
> 
> 
> # QUANTILES
> 
> quant=.25
> 
> 
> 
> fr.1=function(bhat.opt)
> 
> {
> 
>  uu=y-x%*%bhat.opt
> 
>  sample.cond.quantile=quantile(uu,quant)
> 
>  w.less=rep(0,nn)
> 
>  for(ii in 1:nn){if(uu[ii]<sample.cond.quantile) w.less[ii]=1}
> 
> 
> 
>  sum((quant-1)*sum((y-x%*%bhat.opt)*w.less) #negative residuals
> 
>      +quant*sum((y-x%*%bhat.opt)*(1-w.less))) #positive residuals
> 
> }
> 
> start<-c(0,0)
> 
> result=optim(start,fr.1)
> 
> bhat.cond=result$par
> 
> 
> 
> #Quantile Command Results
> 
> fit.temp=rq(y~x1,tau=quant)
> 
> fit.temp
> 
> 
> 
> #OPTIMIZER Results
> 
> bhat.cond
> 
> 
> 
> #OLS Command Results
> 
> mean=lm(y~x1)
> 
> mean
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.