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: