## Wednesday, May 9, 2012

### Oh wait! Quantile regression wins!

* Stata simulation to estimate performance of quantile regression at
* different quantiles compared with other alternatives.

clear
set obs 10000
set seed 101

gen x1 = runiform()
* Single exogenous x

gen u = rnormal()
* Single uncorrelated error u

gen beta=1
label var beta "True Beta"
* At first beta is = 1

gen y0 = u

gen y = beta*x1 + y0
* Generate a starting point for y

forv i=1/20 {
cap drop ytile
* Removes old ytile from the data
xtile ytile=y, nq(100)
* Generates a variable that orders all of the
* y variables into 100 quantiles
qui replace beta = 8-((200+ytile)*.21)^.5
* beta goes from  4-31/20)^.5 to 4-(133/20)^.5 which is
* a range of 1.25 to 2.6 ish.
replace y = beta*x1 + u + 190

di "`i'"
}
* This repetitive routine makes it so that beta is smoothly related to y

replace y = y-190
* Make y closer a small range so that the beta hat seems more effective

line beta y, sort
* We can see the beta function is changing effect accross the different
* quantiles of y.  We can see that if x1 is a policy variable then because
* x1 is greater for lower ys and smaller for higher ys.  Then the effect of
* sign of the betas is to keep the variance in y constant or even decrease
* it.  It is possible for the variance in y to still be larger if the
* variances of the xs are larger since in effect we are adding two random
* variables together.

sum y y0
gen beta_hat=.
label var beta_hat "Estimated Beta from qreg"
gen beta_hat1=.
label var beta_hat1 "Estimated Beta from reg at quantile"
gen beta_hat2=.
label var beta_hat2 "Estimated Beta from reg at quantile with overlap"
gen beta_hat3=.
label var beta_hat3 "Estimated Beta from qreg at quantile with overlap"

forv i =1/100 {
cap qreg y x1, quantile(`=`i'/100')
* This tells stata to do a quantile regression at every point
* the 'cap' tells it not to stop if there is an error.

if _rc==0 qui replace beta_hat=_b[x1] if ytile==`i'
* _rc==0 then it means there was no error in the previous quantile regression

* This does a regression for each of the 100 quantiles
cap reg y x1 if ytile==`i'
* This tells stata to do a regression only in the specific quantile
* the 'cap' tells it not to stop if there is an error.

if _rc==0 qui replace beta_hat1=_b[x1] if ytile==`i'
* _rc==0 then it means there was no error in the previous regression

* This does a regression for each of the 100 quantiles but allows for
* overlap.  We expect that the estimated betas will not be independent
* because there is overlap in the sample.
cap reg y x1 if ytile<= `i'+15 & ytile>= `i'-15
* This tells stata to do a regression only in the specific quantile
* the 'cap' tells it not to stop if there is an error.

if _rc==0 qui replace beta_hat2=_b[x1] if ytile==`i'
* _rc==0 then it means there was no error in the previous regression

cap qreg y x1 if ytile<= `i'+15 & ytile>= `i'-15
* This tells stata to do a regression only in the specific quantile
* the 'cap' tells it not to stop if there is an error.

if _rc==0 qui replace beta_hat3=_b[x1] if ytile==`i'
* _rc==0 then it means there was no error in the previous regression

di _continue " `i'"
}

two (line beta y, sort) (line beta_hat y, sort) (line beta_hat1 y, sort) (line beta_hat2 y, sort), title(Results of Qreg)

two (connected beta ytile, sort) (line beta_hat ytile, sort) (line beta_hat1 ytile, sort) ///
(line beta_hat2 ytile, sort) (line beta_hat3 ytile, sort), title(Qreg against semi-parametric alternatives) legend(on size(vsmall))

* Compared to other estimators, quantile regreesion is working awesome!