Showing posts with label Instrumental Variables. Show all posts
Showing posts with label Instrumental Variables. Show all posts

Wednesday, October 9, 2013

Weak Instruments with a Small Panel

This post is in response to a comment by Luis Enrique on my previous post:

"As somebody who regularly consumes cross-country empirical research based on IV
regressions with samples of 50-100, I found this quite alarming. But then most
of the papers I read will be panel, with T of let's say 50.

This question may reveal shocking ignorance, but if the number of observations
in a panel (N*T) is say 100 * 50, does that translate into a (very) safe
sample size?" - Luis

My response:

Thanks for asking. I am no expert on time series so my opinion should really
be regarded with a grain of salt. First off, the key to remember is that the
problem is the result of a weak instrument. As the first stage R-squared on
the 2SLS approaches 1 the IV estimator becomes the OLS estimator which has
the same properties of OLS.

However, as the first stage R-squared approaches zero we start suffering
serious problems with IV in terms of both bias and efficiency.

Based on the cross sectional simulation, I would say that given a similarly
 weak instrument, if each year's observation within each country is independent
then a panel of 100*50 should not be a problem. However, we expect explanatory
variables to be serially correlated and instruments to be often a 1 time only policy
changes which are extremely serially correlated which reduces the effective sample size
since each observation of the instrument can no longer be seen as independent.

Let's see how data generated in this manner may behave.

* First we define the original weakreg simulation in which there is 
* only cross sectional data.

cap program drop weakreg
program weakreg, rclass
  clear
  set obs `1' 
  * The first argument of the weakreg command is the number of 
  *  observations to draw.
  gen z = rnormal()
  gen w = rnormal()
  gen x = z*.2 + rnormal() + w
  gen u = rnormal()
  gen y = x + w + u*5
  reg y x
    return scalar reg_x = _b[x]
 return scalar reg_se_x = _se[x]
  ivreg y (x=z)
    return scalar ivreg_x = _b[x]
 return scalar iv_se_x = _se[x]
end 

* Now we define the original weakreg2 simulation in which there is 
* panel data.

cap program drop weakreg2
program weakreg2, rclass
  clear
  set obs `1' 
  * The first argument of the weakreg command is the number of 
  *  clusters to draw.
  gen id = _n

  * There is a one time policy change
  gen z = rnormal()

  * The second argument is the number of observations in each cluster
  expand `2' 
  
  bysort id: gen t = _n
  
  gen w = rnormal()
  
  * There is no policy effect before half way through the time period.
  replace z = 0 if (t < `2'/2)

  gen x = z*.2 + rnormal() + w
  
  gen u = rnormal()
  
  gen y = x + w + u*5 if t==1
  bysort id: replace y = x + w + u*12.5^.5 + u[_n-1]*12.5^.5 if t>1
  
  reg y x
    return scalar reg_x = _b[x]
 return scalar reg_se_x = _se[x]
  ivreg y (x=z), cluster(id)
    return scalar ivreg_x = _b[x]
 return scalar iv_se_x = _se[x]
end 

* Looking at first the cross sectional data:
simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(1000): weakreg 5000
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |      1000    1.490732    .0499934    1.29779   1.661524
    reg_se_x |      1000    .0500304    .0007235   .0481176   .0520854
     ivreg_x |      1000    .9938397    .3684834  -.2149411   1.970256
     iv_se_x |      1000    .3638622     .037906   .2568952   .5094048

*/

* We see everything seems to be working very well.

* Looking now at the case where there is 50 clusters and 100 observations
* in each of them.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(1000): weakreg2 50 100
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |      1000    1.493905    .0508205   1.269596   1.670965
    reg_se_x |      1000    .0502538    .0007578   .0478791   .0531899
     ivreg_x |      1000    1.021512    .7299972  -1.083777    3.64094
     iv_se_x |      1000    .7110581    .1806363   .3299663   1.566804
*/

* We can see that there is no huge bias in the ivregression though the
* standard errors are about twice that of the cross sectional data.

* Looking at first the cross sectional data:
simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(1000): weakreg 1000
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |      1000    1.483811    .1117863   1.130268     1.8281
    reg_se_x |      1000    .1121205    .0034219   .1026788      .1249
     ivreg_x |      1000    .9028068    .9626959  -8.379108   4.424785
     iv_se_x |      1000    .9215948    .7681726   .4671113   23.61619
*/

* We can see the IV estimator has a large variance and a significant bias
* though it seems to be doing well in general.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(1000): weakreg2 50 20
sum

/* 
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |      1000     1.49184    .1104733   1.148331   1.927401
    reg_se_x |      1000    .1125144    .0040137   .1009569   .1271414
     ivreg_x |      1000    .8209362    5.213284  -151.3602   17.03425
     iv_se_x |      1000    10.00779    257.6449   .5320991   8149.177
*/

* We can see the IV estimator though dealing with the same number of 
* observations has a huge variance and an even larger bias!
* The take away is that even small panels seem like they work so long
* as there is sufficient observations over time.
Formatted By Econometrics by Simulation

Tuesday, October 8, 2013

Finite Sample Properties of IV - Weak Instrument Bias

* There is no proof that an instrumental variables (IV) estimator is unbiased.

* In fact we know that in small enough samples the bias can be large.

* Let's see a simple setup with the endogeneity a result of omitted variable bias.

* Our instrument is valid, though biased because we are using a "small" sample and the instrument is weak.

clear
set obs 1000

gen z = rnormal()

gen w = rnormal()

gen x = z*.3 + rnormal() + w

gen u = rnormal()

gen y = x + w + u*5

reg y x

ivreg y (x=z)
* IVreg includes the true estimate in the confidence interval though the interval is quite wide.

* This is largely the result of z being a weak instrument for x
reg x z

* There is a conjecture that the IV estimator is biased in finite samples.

* In order to examine this bias we will run a monte carlo 
*  simulation to see how biased our estimates are at each level.

cap program drop weakreg
program weakreg, rclass
  clear
  set obs `1' 
  * The first argument of the weakreg command is the number of 
  *  observations to draw.
  gen z = rnormal()
  gen w = rnormal()
  gen x = z*.2 + rnormal() + w
  gen u = rnormal()
  gen y = x + w + u*5
  reg y x
    return scalar reg_x = _b[x]
 return scalar reg_se_x = _se[x]
  ivreg y (x=z)
    return scalar ivreg_x = _b[x]
 return scalar iv_se_x = _se[x]
end 

* With only 100 observations
simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(10000): weakreg 100
sum
/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |     10000    1.484123    .3619725   .1481085   2.739473
    reg_se_x |     10000     .357267    .0359624   .2261475   .5157685
     ivreg_x |     10000    14.45544    846.8431  -1318.675   78604.67
     iv_se_x |     10000    208878.2    1.93e+07   .6884335   1.92e+09
*/

* We can see the mean standard error estimate is much
* larger than the standard deviation of the estimates.

* In addition, the apparent bias of the IV is huge!
* Thus OLS is the better estimator in this case.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(10000): weakreg 300
sum

/*

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |     10000     1.48824    .2065498   .6367978   2.236391
    reg_se_x |     10000     .204807    .0117558   .1625977   .2537656
     ivreg_x |     10000     .883504    16.43395  -489.2355   1346.839
     iv_se_x |     10000     103.928    7065.172   .5418385   696229.8

*/
* Increasing the sample size to 300 vastly improves the IV estimator.
* Though it is now downward biased.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(10000): weakreg 500
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |     10000    1.490552     .159073   .9184232   2.088732
    reg_se_x |     10000    .1584161    .0071414   .1326629   .1896385
     ivreg_x |     10000     .841985    4.533252  -337.8417   41.04016
     iv_se_x |     10000    10.02729    672.1545   .5350561   66082.69
*/
* Increasing the sample size to 500 does not seem to improve the bias 
* of the IV estimator. Though the standard errors on average seem to be
* getting closer to the standard deviations of the estimators.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(10000): weakreg 750
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |     10000    1.491437    .1291091   1.012212   1.998231
    reg_se_x |     10000    .1292925    .0047437   .1129745   .1498492
     ivreg_x |     10000    .9714284    1.087322  -12.33198   7.718197
     iv_se_x |     10000    1.080065    .8625042   .4623444   39.95131
*/

* Increasing the sample size to 750 dramatically improves the IV estimator.
* It is still slightly biased but that is not a huge problem.
* Now the standard errors are working very well as well.
* The only problem would be the IV estimator still has such large variation
* that both the OLS estimator and the 0 coefficient would be included in
* most confidence intervals.

simulate reg_x=r(reg_x)     reg_se_x=r(reg_se_x)  /// 
         ivreg_x=r(ivreg_x) iv_se_x=r(iv_se_x)    /// 
   , rep(10000): weakreg 1000
sum

/*
    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       reg_x |     10000    1.488341    .1107187   1.064637   1.938462
    reg_se_x |     10000     .111871    .0035312   .1001198   .1255203
     ivreg_x |     10000    .9691499    .8924782  -5.659174    5.94314
     iv_se_x |     10000    .8812981    .2897863     .44236   6.674775
*/

* We can see that our primary gains from more observations is a smaller
* standard error.

Formatted By Econometrics by Simulation

Wednesday, February 20, 2013

Interpreting the Control Function Coefficient


* Is the control function coefficient a measure of the direction and size of the bias caused by endogeneity?

* Imagine the endogenous variable w being composed of three components: 1 endogenous portion, 2 exogenous portion correlated with z, 3 exogenous portion uncorrelated with z.

clear
set obs 10000

gen z = rnormal()
gen v = rnormal()

gen endogenous = rnormal()
gen exog_with_z = z
gen exog_without_z = rnormal()

gen w = endogenous + exog_with_z + exog_without_z

* Likewise we can think of the error u as composed of both an exogenous portion and an endogenous portion (correlated with part of w)

gen u = endogenous + rnormal()*3

gen y = 1*w + 3*u

reg y w
* We can see that OLS is clearly upward biased

ivreg y (w=z)
* Instrumental variables seems to be working well

* Now for the control function
reg w z

predict v_hat, resid

reg y w v_hat
* I was thinking that the control function coefficient could generally not be interpretted directly the sign of the bias but looking at this simulation it appears I was wrong.

* I will have to do some more thinking on this.

Thursday, February 7, 2013

2SLS with multiple endogenous variables


* I am wondering if when using 2SLS you must use a multivariate OLS in the reduced form or if you can just do each individual endogenous variable.

* Let's see!

clear
set obs 10000

* First generate the instruments
gen z1 = rnormal()
gen z2 = rnormal()

* Now the error.  u1 and u2 are the parts of w1 and w2 correlated with the error u.
gen u1 = rnormal()
gen u2 = rnormal()
gen u3 = rnormal()*3

gen w1 = 1*z1 + .5*z2 + rnormal() + 2*u1
gen w2 = -.5*z1 + 1.5*z2 + rnormal() - 2*u2 + .2*u1

gen u = u1 + u2 + u3
gen y = 4 + w1*1 + w2*1 + u

* We can see because u is correlated with w1 and w2, OLS will be biased and inconsistent.

reg y w1 w2

* Instrumental variable regression could solve this problem

ivreg y (w1 w2 = z1 z2)

* Let's see about our 2SLS (which should be the same as IV)

reg w1 z1 z2
predict w1_hat

reg w2 z1 z2
predict w2_hat

reg y w1_hat w2_hat
* It seems that 2SLS using separate regressors is producing the same results.

* This is probably because in SOLS if you use the same regressors then I think the coefficients are the same but the standard errors may be adjusted for cross equation correlations (I think I recall).

Monday, February 4, 2013

Potential IV Challenges even with RCTs


* You design an RCT with a single experimental binary treatment.

* But you think your instrument has the power to influence multiple endogenous response variables.

* We know we cannot use a single instrument to instrument for more than one variable jointly.

* Does it still work as a valid instrument if it affects other variables that can in turn affect our outcome?

clear
set obs 1000

gen z = rbinomial(1,.5)

gen end1 = rnormal()
gen end2 = rnormal()
* end1 and end2 are the endogenous component of x1 and x2 which will bias our results if we are not careful.

gen x1 = rnormal()+z+end1
gen x2 = rnormal()+z+end2

gen u = rnormal()*3

gen y = x1 + x2 + end1 + end2 + u

reg y x1 x2
* Clearly the OLS estimator is biased.

ivreg y (x1=z)
ivreg y (x2=z)
* Trying to use the IV individually only makes things worse.

ivreg y x1 (x2=z)
* Trying to control for x1 does not help.

* So what do we do?

* If we have control over how the instrument is being constructed then we think things through carefully making sure that if we are interested on the effect of x1 on y then our instrument only results in a direct effect and no additional affect on x2.

* Unfortunately, this task can not always be feasible.

* Another consideration we should take into account is potential variable responses to the istrument.

* For instance:

clear
set obs 1000

gen z = rbinomial(1,.5)

gen end1 = rnormal()
* end1 and end2 are the endogenous component of x1 and x2 which will bias our results if we are not careful.

gen g1 = rnormal()
gen grc = g1+1
  label var grc "Random Coefficient on the Instrument"

gen x1 = rnormal()+z*grc+(end1)

gen u = rnormal()*3

gen brc = g1+1
  label var brc "Random Coefficient on the Endogenous Variable"

gen y = x1*brc + end1 + u

reg y x1
* OLS is still biased.

ivreg y (x1=z)
* Unfortunatley, now too is IV.

* This is because, though our instrument is uncorrelated our errors, the response to the instrument is variable and that response may also be correlated with a variable response on the variable of interest.

* Thus, though corr(u,z)=0 and cor(x,z)!=0, we can still have problems.

* This particular problem is actually something that I am working on with Jeffrey Wooldridge and Yali Wang.

Sunday, August 12, 2012

Mata: program your own IV Estimator


* Mata: program your own IV Estimator

* First let see how we could program an IV estimator without matrices.
version 11

clear
set obs 1000

gen z = rnormal()

gen u= rnormal()

gen x = u + z + rnormal()

gen y = 2 + 2*x + 15*u

reg y x

foreach v in x y z {
  foreach vv in x y z {

  gen `v'`vv' = `v'*`vv'
  qui sum
  local `v'`vv' = r(mean)*r(N)
  }
}

di "IVhat = " (`xz'*(`zz')^-1*`zx')^-1 * (`xz'*(`zz')^-1*`zy')

ivreg y (x=z)
* Looks pretty close.  However, there is some divergence due to the above estimator not allowing for a constant.

ivreg y (x=z), nocon

* Now let's see how to do it in Mata.

mata

 x = st_data(.,("x"))
 z = st_data(.,("z"))
 y = st_data(.,("y"))

 IVhat1 = invsym((x'*z)*invsym(z'*z)*(z'*x))*(x'*z)*invsym(z'*z)*(z'*y)

"IV estimate without constant"
 IVhat1

"Create a constant vector of length x"
 cons = J(rows(x), 1, 1)

 X = (x, cons)
 Z = (z, cons)

 IVhat2 = invsym((X'*Z)*invsym(Z'*Z)*(Z'*X))*(X'*Z)*invsym(Z'*Z)*(Z'*y)

"IV estimate with constant"
 IVhat2

end
* End is necessary with mata do files to indicate when the mata subsystem is completed.

Wednesday, May 16, 2012

Testing the Rank Condition of IV Estimators


* One of the key assumptions of the IV estimators
* is that the rank E(Z'X)=K=# of elements of B.

* It is difficult to test the rank condition in
* general.  However following Wooldridge 2010 it
* is easy to test the rank condition of instrumental
* variables when there is only on endogenous variable
* xk

* xk = d1 + d2x2 + ... + dk-1 xk-1 + theta1*z1 + r

* In order for the rank condition to hold E(theta)!=0

*******************************************************
* Let's see this test in action:

clear
set obs 1000
set seed 100

gen x1 = rnormal()
gen x2 = rnormal()
gen x3 = rnormal()

* generate the unobserved factor V which contributes
* to create a correlation between xk and u.
gen V = rnormal()

* generate the instrumental variable
gen z = rnormal()

* generate xk the suspected endogenous variable.
gen xk = rnormal() + V + z

gen u = rnormal()*2 - 2*V

gen y = 1*x1 + 2*x2 + 4*x3 + xk + u

* Simulation End

reg y x1 x2 x3 xk
* We can see that xk is biased

* In order to test the rank condition of the IV estimator
* we estimate

reg xk x1 x2 x3 z
* z appears to be strongly significant.  Therefore we
* conclude that the rank condition is sustained.
* Therefore we can safely use the IV estimator:

ivreg y x1 x2 x3 (xk=z)
* So the IV estimator seems to be working pretty well.

* When might we expect the rank condition to fail?

*******************************************************
* When the instrumental variable is weak (has small
* explanatory power relative to xk):

clear
* Set the number of observations to generate to 1000.
set obs 1000
set seed 101

* generate the instrumental variable
gen z = rnormal()

gen x1 = rnormal()
gen x2 = rnormal()
gen x3 = rnormal()

* generate the unobserved factor V which contributes
* to create a correlation between xk and u.
gen V = rnormal()

* generate xk the suspected endogenous variable.
gen xk = 2*rnormal() + V + z*.1

reg xk x1 x2 x3 z
* We cannot see a statistical difference between
* theta and 0 (though we know one exists).

* note, just because our Rank condition test fails
* does not mean that the underlying rank condition
* fails.  In this case E(Z'X)!=0 therefore the
* rank condition does not fail in the population
* In order to see this trying increasing the number
* of observations.

* Generate the endogeneity between xk and yb
gen u = rnormal()*2 - 2*V

gen y = 1*x1 + 2*x2 + 4*x3 + xk + u

ivreg y x1 x2 x3 (xk=z)
* The IV estimator is simply too weak to get
* good estimates.

reg y x*
* Once again the OLS regression works pretty well even
* when the IV estimator is not working well.

Tuesday, May 15, 2012

The Weak Instrument Problem

* Instrumental variables are extremely useful in
* economics in general because they provide a
* mechanism to get around the issue that some
* important predictive variables of interest
* in economics might also be correlated with
* other unobservable factors that might affect
* the outcome.

* A typical example is years of education and wages:
* wage = alpha + beta*years_of_education + u

* The problem is that the number of years of
* education might be related to motivation which
* might also be related to u.  So
* corr(years_of_ed, u) != 0 and the
* normal OLS estimators will be biased
* because:

* Y = XB + u
* X'Y = X'XB + X'u
* (X'X)^-1 X'Y = B + (X'X)^-1 X'u

* If we assume that X'u are uncorrelated and
* the Xs are fixed then an unbiased estimate
* of B is B_hat:

* E((X'X)^-1 X'Y) = B + E((X'X)^-1 X'u)
* E((X'X)^-1 X'Y) = B
* B_hat = (X'X)^-1 X'Y

* But if corr(X,u)!=0 then:
* B_hat = (X'X)^-1 X'Y
* B_hat = (X'X)^-1 X'(XB + u)
* B_hat = (X'X)^-1 X'XB + (X'X)^-1 X'u
* B_hat = B + (X'X)^-1 X'u
* E(B_hat) = B + E((X'X)^-1 X'u)

* to show unbiasedness: Assume fixed Xs
* E(B_hat) = B + (X'X)^-1 E(X'u)
* E(X'u)!=0 if corr(X,u)!=0
* E(B_hat) = B

* to show consistency:
* plim(B_hat) = B + (X'X)^-1 plim(X'u)
* plim(X'u)!=0 if corr(X,u)!=0
* plim(B_hat) = B

* So in order to get around this economists
* often use the instrumental variables estimator.

* Imagine that there is some instrument Z
* which is uncorrelated with u and
* is correlated with X. Z must also be
* uncorrelated with B in the standard case.

* Y = XB + u
* Z'Y = Z'XB + Z'u
* (Z'X)^-1 Z'Y = (Z'X)^-1 Z'XB + (Z'X)^-1 Z'u
* (Z'X)^-1 Z'Y = B + (Z'X)^-1 Z'u

* Since we have assumed corr(Z'u) = 0
* Z'u=0 at the population level.

* Thus we get the IV estimator:
* (Z'X)^-1 Z'Y = B_hat
* B_hat = (Z'X)^-1 Z'Y

* To show consistency:
* B_hat = (Z'X)^-1 Z'Y
* B_hat = (Z'X)^-1 Z'(XB + u)
* B_hat = (Z'X)^-1 Z'XB + (Z'X)^-1 Z'u
* B_hat =  B + (Z'X)^-1 Z'u

* plim(B_hat) =  B + plim((Z'X)^-1 Z'u)
* plim(B_hat) =  B + (Z'X)^-1 plim(Z'u)
* plim(B_hat) =  B

* Note: the IV estimator is not unbiased but
* it is consistent!

* However, if corr(Z,u) != 0 even if it is small
* and the instrument is "weak" that is
* |Z'X| is small then we have a weak instrument
* problem.  This is because (Z'X)^-1 is 1 over
* a small number which will means that
* even a slight correlation in corr(Z,u) can cause
* a large bias in the IV estimator.

* Let's see how that works!

clear
set obs 1000

* First the standard IV estimator
gen Z1 = rnormal()

* However, Z is not correlated with u yet.
gen u1 = rnormal()

* So Z is weak and remember u1 has to be
* correlated with X in order to justify
* use of the IV estimator:
gen X1 = rnormal()*10 + Z1 + u1*2

reg X1 Z1

gen Y1 = 3*X1 + 5*u1

ivreg Y1 (X1=Z1)
* We can see that Z is a pretty effective
* instrument at estimating B even if
* the correlation between Z and X is small.

******************************************
* Now let's see what happens when there
* is a small correlation between Z and u

* Imagine there is some additional
* explanatory variable V which is unobserved
* but partially explains the instrument
* as well.

gen V = rnormal()

* First the standard IV estimator
gen Z2 = rnormal() - V

* However, Z is somewhat correlated with u.
gen u2 = 10*rnormal() + V/4

* So Z is weak and remember u2 has to be
* correlated with X in order to justify
* use of the IV estimator:
gen X2 = rnormal()*10 + Z2/4 + .1*u2

reg X2 Z2

gen Y2 = 3*X2 + 5*u2

ivreg Y2 (X2=Z2)
* In this case we get a better estimate of
* B by using the known biased OLS estimator
* than by using the instrumental variables
* estimator.

reg Y2 X2

* Even though the correlation between the
* instrument and the error is small
pwcorr Z2 u2, sig

******************************************
* This is primarily a function of the weakness
* of Z at explaining X.  See what happens
* if Z has more explanatory power

gen X3 = rnormal()*10 + 10*Z2 + .1*u2

reg X3 Z2

gen Y3 = 3*X3 + 5*u2

ivreg Y3 (X3=Z2)
* Now even though Z is slightly correlated
* with the error u.

reg Y2 X2

* Even though the correlation between
* Z and u is the same as previously,
* the strength of the instruments
* in explaining X wins out and gives us a
* better estimator than OLS.

Sunday, May 13, 2012

Average Treatement Effects and Correlated Random Coefficients


* Average Treatement Effects and Correlated Random Coefficients

* Wooldridge, J. (2003). Further results on instrumental
* variables estimation of average treatment effects in
* the correlated random coefficient model. Economics
* Letters, 79(2), 185-191. doi:10.1016/S0165-1765
* (02)00318-X

* This simulation follows the paper by Wooldridge which
* shows that it is possible to consistently estimate
* the population average of the correlated random
* coefficient (CRC) model with multiple treatment
* variables under the standard assumptions under
* which instrumental variables (IV) estimators are
* consistent.

* In order to understand CRC think of the following
* model. Equation (1):

* E(y|a,b,w)=a+wb=a + b1*w1 + b2*w2 ... + bg*wg     (1)

* b can depend on unobserved heterogeniety as well
* as w.  Writing the equation in error form (2):

* y = a + wb + e, E(e|a,b,w)=0

* b can vary for every individual observation i.

* Therefore we are not trying to estimate each individual
* b which is impossible, but instead we are attempting
* to estimate the E(b) or the average treatment effect (ATE).

* This in general is difficult when bj is correlated with
* unobserved heterogeniety of individual j.

* In order to tackle this problem Wooldridge introduces
* a instrumental variable which is redundant or ignorable
* in covariates x and instruments z ie.

* E(y|a,b,w,x,z)=E(y|a,b,w)                          (6)

* The next assumption is what seperates xs from zs.

* E(a|x,z)=E(a|x)=gamma0 + x*gamma                   (7)

* This says that the mean a can depend on the explanatory
* variables x but not on the instrumental variable z.

* E(bj|x,z)=E(b|x)=beta0 + (x-E(x))*deltaj

* This assumption says that the average coefficient can
* depend on the explanatory variables x but not on the
* instrumental variable z.

* We can write a = gamma0 + x*gamma + c , E(c|x,z)=0 (8)
* and b=beta0 + (x-E(x))*deltaj + vj, E(vj|x,z)=0    (9)

* Ultimately by substituting this back into (1) we get:

* y=gamma0 + xgamma + wbeta + w1(x-E(x))*delta1 +....
*  wG(x-E(x))*deltaG + c + wv + e                    (10)

* The composite error term is c + wv + e.  Under the
* assumptions thus far E(c|x,z)=E(e|x,z)=0.  However
* E(wv|x,z)!=0 because b is generally not a
* deterministic linear function of x.

* Let us simulate up to this point imagining that we do
* have our bs as deterministic linear functions of xs.

clear
set obs 10000

gen x1 = rnormal()
gen x2 = rnormal()

* We will force w to be uncorrelated with x.
gen w1 = rnormal()
gen w2 = rnormal()

* Each idividual has his/her own intercept or starting
* point c
gen c=rnormal()
gen a= -5 + 3*x1 - 2*x2 + c

* Let us first imagine b being a deterministic function
* of observables x.  Each idividual has a unique response
* w.
gen b1 =  1 + .5*x1 + -2*x2
gen b2 = -2 + 1.75*x1 + 3*x2

* Imagine w1 as being years of education, w2 as being
* offered a job right out of college and x1 as intelligence
* and x2 as GPA.  The interesting thing is having a high
* GPA might be correlated with years of experience
* but it also might help explain the effect years of education
* has on y (future income).

* Let us generate our reduced form error
gen e=rnormal()*10

* Now let's generate our y variables
gen y = a + b1*w1 + b2*w2+e

* First let us generate our w1(x-E(x)) variables:
* The average of observed x is of course not the E(x) but it is
* a consistent estimator of E(x).

sum x1
gen w1_x1 = w1*(x1-r(mean))
gen w2_x1 = w2*(x1-r(mean))

sum x2
gen w1_x2 = w1*(x2-r(mean))
gen w2_x2 = w2*(x2-r(mean))

* Now we can estimate all of our coefficients directly using
* equation 10.

* y=gamma0 + x gamma + w beta + w1(x-E(x))*delta1 +....
*  wG(x-E(x))*deltaG + c + wv + e                    (10)
reg y x1 x2 w1 w2 w1_x1 w1_x2 w2_x1 w2_x2

* One can see that in the case where b is completely
* linearly dependant on x the above CRC estimator works fine.

* Stay tuned for what happens when there is some error in b!

Tuesday, May 1, 2012

The Many Forms of Instrumental Variables


* The many forms of the Instrumental Variables estimator
* Stata Simulation and Estimation

* Imagine that you have the endogenous variable, years of 
* education, and you want to estimate the returns of education 
* in terms of earnings.  The problem is that diligence is 
* going to be correlated with education and it is also going to 
* be correlated with earnings. But diligence is unobservable 
* in your data set and perhaps from a feasibility perspective, 
* in life.

* Ideally you would like to have an experiment where you "give" 
* some people more education than others.  We have the next best 
* thing.  Let's imagine that we have education scholarship 
* lottery data which gives students one year of free education 
* upon completing that year and it is awarded completely RANDOMLY 
* among all potential students.

* The randomness is year in this application.  Formally:
* Y=XB+U
* The problem corr(X,U)!=0
* So, rather than solving the standard way:
* X'Y=X'XB+X'U
* B=(X'X)^-1 X'Y-(X'X)^-1 X'U  
* E(B)=B=E((X'X)^-1 X'Y) + 0  ---- because E(X'U)=0, OLS assumption

* We, use:
* Z'Y=Z'XB+Z'U
* B=(Z'X)^-1 Z'Y-(Z'X)^-1 Z'U  
* E(B)=B=E((Z'X)^-1 Z'Y) + 0  ---- because E(Z'U)=0, IV assumption
* It is not obvious from this formation, but IV is not an unbiased estimator,
* just a consistent one.

* For more details see Wooldridge, Book 1 Chapter 15

* Standard deviation of u
gl sdu=5

* Average effect of z on w
gl gamma11=4
gl gamma12=2
gl gamma21=4.4
gl gamma22=0

* Average effect of w on z
gl beta1=1
gl beta2=1

* Standard deviation of z
gl zsd1=1
gl zsd2=1

* Specify the correlation between the explanatory variables x1 and x2 and the error.
gl rho12=.5
gl rho13=.75

gl sdv1 = 1gl sdv2 = 1

drop _all
clear

set obs 10000

gen rv1=rnormal()
gen rv2=rnormal()
gen rv3=rnormal()

gen u =(rv1+rv2+rv3)*$sdu
gen v1 =$sdv1*(rv1*$rho12 + rv2*(1-$rho12)^.5)
gen v2 =$sdv2*(rv1*$rho13 + rv3*(1-$rho13)^.5)

gen z1 = rnormal()*$zsd1
gen z2 = rnormal()*$zsd2
  
gen x1=z1*$gamma11 + z2*$gamma12 + v1
gen x2=z1*$gamma21 + z2*$gamma22 + v2
    
gen y = x1*$beta1 + x2*$beta2 + u

* The most straight forward IV estimator is IV reg
ivreg y (x*=z*)

* An equivalent estimator is 2SLS
reg x1 z*
predict x1hat

reg x2 z*
predict x2hat

reg y x1hat x2hat
* The second stage errors need be adjusted for the first stage 
* being estimated.

* Also, an equivalent estimator is another 2 stage estimator 
* called the control function.
reg x1 z*
predict v1hat, residual

reg x2 z*
predict v2hat, residual

reg y x1 x2 v1hat v2hat