Write your own System IV estimator in R (and SOLS)

# Specify a variable to hold number of observations
obs = 10000

# Create independent variables
z1 = rnorm(obs)
z2 = rnorm(obs)

# Create error
u1 = rnorm(obs)
u2 = rnorm(obs)

# Create endogenous variable
x1 = .5*z1  -   z2 + u1 + rnorm(obs)
x2 = 1.5*z1 - 2*z2 - u2 + rnorm(obs)
x3 = rnorm(obs)

# Create dependent variable
y1 =  5 + 2*x1 - x2 +   x3 + u1*5 - u2*2
y2 = -5 - 2*x1 + x2 - 2*x3 + u2*5 - u1*1

# First let's attempt a nieve SOLS
# Y = XB
# X'Y = X'XB
# (X'X)^-1(X'Y) = Bhat

X = cbind(x1,x2,x3,1)
Y = cbind(y1,y2)

A = solve(t(X)%*%X)
# solve simply finds the inverse of its argument
B = t(X)%*%Y

SOLS = A%*%B
SOLS

# We can see the coefficient on x is upwards biased.

# Now let's construct the system IV

# Y = XB
# Z'Y = Z'XB
# (Z'X)'Z'Y = (Z'X)'Z'XB
# (X'Z)Z'Y  = (X'Z)Z'XB
# ((X'Z)Z'X)^-1 (X'Z)Z'Y = B
# ((X'Z Z'X)^-1 (X'Z Z'Y) = IVhat
# ((X'Z Z'X)^-1 (X'Z Z'Y) = IVhat
# C * D = IVhat

X = cbind(x1, x2, x3, 1)

# It is important to remember that any explanatory variables that are not instrumented for must be included in the instrument.  Thus x3 in Z.
Z = cbind(x3, z1, z2, 1)

C = solve(t(X)%*%Z%*%t(Z)%*%X)
D = t(X)%*%Z%*%t(Z)

SIV = C%*%D%*%Y
SIV

# System IV seems to be working pretty well
# Alternative non-system IV on the two equations would be

IV1 = C%*%D%*%y1
IV1

IV2 = C%*%D%*%y2
IV2

# We can see that there is no difference in coefficient estimates between seperate IV estimates and SIV.  This, I believe, is because like in the SOLS model if the regressors are the same in both equations and there is no cross equation restrictions then the system produces identical results to the non-system estimates.  This however is not true when estimating standard errors.

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.

A note on correlated variables

* If A is correlated with B and B with C, does A need to therefore be correlated with C?

* No.
clear
set obs 10000

gen A = rnormal()
gen C = rnormal()

gen B = A+C

corr A B C

Selection and Bias

* Selection bias is often regarded as a severe condition that if diagnosed, significantly limits the believability of any study.  In this simulation I will show that while selection bias can bias coefficients the type and severity of the selection can create widely different results.

* Generally speaking it is a neccessary assumption that selection is random in order to achieve consistent estimates.

* In this post I will look at several different potential selection mechanisms and examine how these different mechanisms may bias estiamtes.

* I: Selection is based on one or more explanatory variables.  This might be the case if you were wondering what the returns to years of education were but people with more years of education were less likely to respond to the survey.

  clear
  set obs 100000
  * Set the number of observations that will be generated.

  gen x1 = rnormal()
  gen x2 = runiform()
  gen u = runiform()

  gen s = rbinomial(1,x2)
  * Selection is based on x2.  The more years of education in terms of a 0 to 1 scale that a person has the more likely the person is to opt out of the survey.

  gen y = 1 + 1.5*x1 + 2*x2 + 10*(u-.5) if s == 1
  * y is some income index

  reg y x1 x2
  * When selection is based on years experience alone there is no detectable bias in estimates.

  * This type of selection is called "selection at random" in the statics literature.  This is due to selection not being a source of bias when it is based on observable characteristics.  However, this is a poor name for selection because even when based on observable characteristics, selection is not at random and does not have many of the convenient properties of selection truly random, such as the sample being an unbiased estimator of the mean x2.

sum x2 if s==1

* II: Selection is based on the unobserved shock term u.  Imagine that some random of the target population gets hired by private corporations that pay more (thus larger u) but require that employees do not respond to unapproved surveys.

  clear
  set obs 100000

  gen x1 = rnormal()
  gen x2 = runiform()
  gen u = runiform()

  gen s = rbinomial(1,1-u)

  gen y = 1 + 1.5*x1 + 2*x2 + 10*(u-.5) if s == 1

  reg y x1 x2
  * Once again there is no bias on the coefficients of x1 and x2 however there is a bias in the constant estimate.  This is because the sample that responds to the survey on average has lower expected wage than those who do not respond.

  * Thus E(u|s=1)>0

* III: The consistency of the results for mechanism I is a bit misleading.  For one, we are assuming that x2 (education) is uncorrelated with unobserved error u (such as ability, family background, or geographic heterogeneity).  All of these factors are likely to cause difficulties when estimating the returns to education.  However, these problems exist absent of selection.  What we are interested in is how does selection bias estimates.  If we relax the implicit assumption of constant returns to education then it is easier to see how selection biases these results.

  clear
  set obs 100000

  gen x1 = rnormal()
  gen x2 = runiform()
  gen u = rnormal()
  gen v = runiform()

  * Imagine now that there is some unobserved heterogeneity v (hard-workingness) which affects the returns to education.
  gen r = (2*v)

  sum r
  * The expected value of r is 1.

  gen y = 1 + 1.5*x1 + 2*x2*r + 10*u

  reg y x1 x2
  * Without selection, the random coefficient (r) does not bias the average partial effect of x2 on y.
  * One important reason is because cor(x2,r)=cov(x2,r)=0 and
  *  cov(x2,r) = E(x2*r) - E(x2)E(r) =  E(x2*r) - E(x2) = 0
  * => E(x2*r) = E(x2)

  gen s = rbinomial(1,v)
  * Let's imagine that people who are hard working, are more likely to make the time to answer surveys.

  reg y x1 x2 if s == 1
  * Section on v however does bias the coefficient on x2.

  * This is because E(r|s=1)>1

* IV: Finally we will exam is selection based on outcome variable y.  Imagine that people who have less y (income) are less likely to be willing to respond to the survey.  This will cause selection to be correlated with the explanatory variables x1 and x2 as well as the error u.

  clear
  set obs 100000

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

  gen y = 1 + 1.5*x1 + 2*x2 + 7.5*u

  sum y

  * Now we will create an index of y from 0 to 1 from which selection will occur
  gen yp = (y-r(min))/(r(max)-r(min))

  gen s = rbinomial(1,yp)

  reg y x1 x2 if s == 1
  * Section on v however does bias the estimates even if there exists an instrumental variable.

  * This is because Corr(u,s)>0 and Corr(x, s)>0.  This alone does not imply Corr(x,u|s)!=0 however in this case it happens to be the case.

  cor u x1 x2 if s == 1

  * In summary, selection bias can bias the results of OLS estimates if not taken into consideration (II,III,IV).  However, this bias may be small enough as to be trivial in many cases (II, IV).  In contrast when the standard assumption on the constant nature of the coefficients is relaxed then selection can become much more of a biasing factor (III).

on Vacation

Dear Readers, I am on vacation till the 9th of August.  I will continue to post as the spirit takes me but probably not daily posts again until after that.

3 ways to write a probit command

* Probit is a commonly used command to model binary outcome variables.

* The maximum likelihood sytax in Stata is specific
* I will first define three different equivalent ml programs.  Then we will solve them to see how they perform.

cap program drop myprobit1
program define myprobit1
  * Tell Stata that this code is written in version 11
  version 11
  * Define the arguments input in the maximum likelihood function
  * Stata automatically generates the first argument as a temporary variable name
  * that is used to store the natural log of the likelihood value of the likelihood function.
  args ln_likehood xb
  * The probit tries to directly maximize the probability of seeing the outcome (either y==1 or y==0)
  * If y==1 then we maximize the probability of seeing a positive outcome
  qui replace `ln_likehood' = ln(  normal(`xb')) if $ML_y==1
  * If y==0 then we want to maximize the probability of seeing a negative outcome
  qui replace `ln_likehood' = ln(1-normal(`xb')) if $ML_y==0
end

cap program drop myprobit2
program define myprobit2
  version 11
  args ln_likehood xb
  * Because of the symetry of the normal CDF the following set of expressions is equivalent.
  qui replace `ln_likehood' = ln(normal( `xb' )) if $ML_y==1
  qui replace `ln_likehood' = ln(normal(-`xb' )) if $ML_y==0
end

cap program drop myprobit3
program define myprobit3
  version 11
  args ln_likehood xb
  * This is a somewhat opeque way of writing the same thing as above.
  * I do not prefer it to the first two forms except that it is likely to run slightly faster
  * since it is a single command.
  qui replace `ln_likehood' = ln(normal((-1)^(1-$ML_y)*`xb'))
end


clear
set obs 1000

gen x1 = rnormal()*.5
gen x2 = rnormal()*.5
gen u = rnormal()*(.5^.5)

gen p=normal(x1+x2+u)

gen y=rbinomial(1,p)

probit y x1 x2
* This is the benchmark

ml model lf myprobit1 (y=x1 x2)
ml maximize
* Basically works the same except there is no pseudo r2 reported.

timer on 1
ml model lf myprobit2 (y=x1 x2)
ml maximize
timer off 1

timer on 2
ml model lf myprobit3 (y=x1 x2)
ml maximize
timer off 2

timer list

di "myprobit2 is " r(t1)/r(t2) " times slower than myprobit3"

* The following is a brief exploration attempting to run two probits at once.
* It shows how ml can optimize two equations jointly
cap program drop mytwoprobit1
program define mytwoprobit1
  version 11
  args ln_likehood xb1 xb2
  qui replace `ln_likehood' = ln(normal((-1)^(1-$ML_y1)*`xb1')) ///
                            + ln(normal((-1)^(1-$ML_y2)*`xb2'))
end

* This is equivalent to the previous except that it uses the binormal distribution rather than two seperate normals.  I tried specifying the correlation as a third equation but that did not work.  I will need to play around a bit more to figure out how to program a biprobit.
cap program drop mytwoprobit2
program define mytwoprobit2
  version 11
  args ln_likehood xb1 xb2
  qui replace `ln_likehood' = ln(binormal((-1)^(1-$ML_y1)*`xb1', ///
                                          (-1)^(1-$ML_y2)*`xb2', ///
 0))

end

clear
set obs 1000

gen x1 = rnormal()*.5
gen x2 = rnormal()*.5

gen p1=normal(x1-x2)
gen p2=normal(-x1+x2)

gen y1=rbinomial(1,p1)
gen y2=rbinomial(1,p2)

ml model lf mytwoprobit1 (y1:y1=x1 x2) (y2:y2=x1 x2)
ml maximize

ml model lf mytwoprobit2 (y1:y1=x1 x2) (y2:y2=x1 x2)
ml maximize

The Joy of Logic - Stata logical operators

* Both Stata and R handle logical operators in a similar fassion

* There are two types of operators & (and) and | (or)

if (1==2 & 2==2) di "both 1 = 2 and 2 = 2"
* This message is not displayed

if (1==2 | 2==2) di "either 1 = 2 or 2 = 2 (obviously 2==2)"

* Operators can easy switch values
if !(1==2 | 2==2) di "neither 1 = 2 nor 2 = 2"
  * Is not displayed

* You might be wondering what the value of these operators are

di "True: " 1==1 ", False: " 3==2

* Logical operators can also take order of operations.

* Thus

if (1==1|2==3)|(3==3&3==2) di "Will display because either 1==1 or 2==3 (obviously 1==1) or 3==3 and 3==2 (and 3==3)"
if (1==1|2==3|3==3)&(3==2) di "Will not display because either 1==1, 2==3, or 3==3 (we know 1==1, and 3==3) but 3 does not equal 2."

* Order of operations matters.

* Any logical operator combination can be specified two ways.
if (1==1|2==3) di "Straightforward method"
if !(1!=1&2!=3) di "Alternative display (more clunky)"