R help - Jun 2012 - Quantile regression: Discrepencies Between optimizer and rq()

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]]

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.