Write your own estimator and bootstrap the standard errors

* Let's imagine that you have some new command that you would like to use and you do not know how to solve the math to calculate the standard errors.

* For instance, two stage quantile regression would be extremely challenging.
clear
set obs 1000

* There is an exogenous component to w which is z
gen z = rnormal()

* There is an endogenous component to w which is v
* It is drawn from a Cauchy distribution.
* See (previous post) for more information on drawing from non-standard distributions.
gen v = tan(_pi*(runiform()-.5))

gen w = z + v

* There is also a standard addative error in the model
* u1 is a portion of the error that is independent of w
* It is drawn from a Cauchy distribution.
gen u1 = tan(_pi*(runiform()-.5))

gen u = (v+u1*3)*5

* Thus y is generated this way
gen y = -14 + 5*w + u

* In a Cauchy distribution the median is defined but the mean and variance are not.
* Thus it is inappropriate and often infeasible to use OLS or many common estimators.

* Quantile regression however is theoretically justified.

* That said:
reg y w
* Looks like it is working fine.

* But just for the sake of doing the consistent thing:
qreg y w

* We can see that both estimators are heavily biased because of the endogeneity of w.

* We want to estimate a two stage quantile regression
qreg w z

predict w_hat

qreg y w_hat
* It looks pretty good but we need to estimate the standard errors because that generated from just qreg on the predicted values is ignoring the variance in the original prediction.

* We will write a short program:

cap program drop twostageqreg
program define twostageqreg

  qreg w z
    * First let's drop w_hat if it is already defined
    cap drop w_hat

    predict w_hat
  qreg y w_hat

end


* Now let's bootstrap
bs: twostageqreg

* Our two stage bootstrap looks like it is working pretty well.
* The standard errors are smaller than in the second stage estimate.
* This is somewhat surprising.

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!

Quantile Regression Fail

* Stata simulation to estimate performance of quantile regression at
* different quantiles

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 = 4-((ytile+30)/20)^.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"
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
 
  di _continue " `i'"
}

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

two (line beta ytile, sort) (line beta_hat ytile, sort), title(Results of Qreg)

* The quantile regression seems to be working mediocre
* at picking up some of the shape of the beta distribution.
* However, near the edges which people are often interested in
* it is failing.

* One might say that this is due to the tails having too few observations.
* However, increasing the sample size to 100,000 gives almost identical results!

* Ah, but hope is not lost.  Stay tuned for future releases of semi-parametric
* estimators!