2SLS in Mathematica
2SLS data setup. Note that there is RV u that appears in both x and in the error term. There is also an IV, z, that affects x but not u. n = 10000; z = Table[Random[NormalDistribution[0, 1]], {n}]; B0 = 1; B1 = 2; gamma0 = -2; gamma1 = 4; u = Table[Random[NormalDistribution[0, 1]], {n}]; x = u + gamma0 + gamma1 * z + Table[Random[NormalDistribution[0, 1]], {n}]; e = 5*u + Table[Random[NormalDistribution[0, 1]], {n}]; y = B0 + B1*x + e;
Note that the real coefficient on x is 2, but the estimated coefficient biased upwards. Now we can do the first stage:
First stage
iota = Table[1, {n}]; Z = Transpose[{iota, z}]; Gammahat = Inverse[Transpose[Z].Z].Transpose[Z].x; xhat = Z.Gammahat; Xhat = Transpose[{iota, xhat}]; Bhat = Inverse[Transpose[Xhat].Xhat].Transpose[Xhat].y
and now the coefficient estimates are close to the true values: