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: