Note: This model could also be fit with sem, using maximum likelihood instead of a two-step method. 
You can find examples for recursive models fit with sem in the “Structural models 2: Dependencies between endogenous variables” section of [SEM] intro 5 — Tour of models.

Someone posed the following question:

        Y = a + bX + cZ + dW 

        W = e + fX + gZ + hQ 

ivregress will not let you do this and, moreover, if you believe W to be endogenous because it is part of a system, then you must include X and Z as instruments, or you will get biased estimates for b, c, and d.

Consider the system

        Y1 = a0 + a1*Y2 + a2*X1 + a3*X2 + e1               (1)

        Y2 = b0 + b1*Y1 + b2*X3 + b3*X4 + e2               (2)

Warning: Assume we are estimating structural equation (1); if X1 and X2 are exogenous, then they must be kept as instruments or your estimates will be biased. In a general system, such exogenous variables must be used as instruments for any endogenous variables when the instrumented value for the endogenous variables appears in an equation in which the exogenous variable also appears.

Consider the reduced forms of your two equations:

        Y1 = e0 + e1*X1 + e2*X2 + e3*X3 + e4*x4 + u1        (1r)

        Y2 = f0 + f1*X1 + f2*X2 + f3*X3 + f4*x4 + u2        (2r)

where e# and f# are combinations of the a# and b# coefficients from (1) and (2) and u1 and u2 are linear combinations of e1 and e2.

All exogenous variables appear in each equation for an endogenous variable. This is the nature of simultaneous systems, so efficiency argues that all exogenous variables be included as instruments for each endogenous variable.

Here is the real problem. Take (1): the reduced-form equation for Y2, (2r), clearly shows that Y2 is correlated with X2 (by the coefficient f2). If we do not include X2 among the instruments for Y2, then we will have failed to account for the correlation of Y2 with X2 in its instrumented values. Since we did not account for this correlation, when we estimate (1) with the instrumented values for Y2, the coefficient a3 will be forced to account for this correlation. This approach will lead to biased estimates of both a1 and a3.

For a brief reference, see Baltagi (2011). See the whole discussion of 2SLS, particularly the paragraph after equation 11.40, on page 265. (I have no idea why this issue is not emphasized in more books.)

Failing to include X4 affects only efficiency and not bias.

However, there is one case where it is not necessary to include X1 and X2 as instruments for Y2. That is when the system is triangular such that Y2 does not depend on Y1, but you believe it is weakly endogenous because the disturbances are correlated between the equations. You are still consistent here to do what ivregress does and retain X1 and X2 as instruments. They are, however, no longer required. Then you could do what you suggested and just regress on the predicted instruments from the first stage.

If you do use this method of indirect least squares, you will have to perform the adjustment to the covariance matrix yourself. Consider the structural equation

        y1 = y2 + x1 + e 

where you have an instrument z1 and you do not think that y2 is a function of y1.

The following example uses only z1 as an instrument for y2. Let’s begin by creating a dataset (containing made-up data) on y1, y2, x1, and z1:

Now we perform the first-stage regression and get predictions for the instrumented variable, which we must do for each endogenous right-hand-side variable.

Now we correct the variance–covariance by applying the correct mean squared error:

http://www.stata.com/support/faqs/statistics/instrumental-variables-regression/

本来我们完全可以使用 IVREG 这个简单的命令，但是如果工具变量只有一个，你们这里就给出了答案。

use  yourdata,clear

 regress x2  z1

predict double yhat 

regress y  yhat x1 

rename yhat yhold

rename x2 yhat

predict double res, residual

rename yhat x2                       

rename yhold yhat  

replace res = res^2  

summarize res

scalar realmse = r(mean)*r(N)/e(df_r) 

matrix bmatrix = e(b)

matrix Vmatrix = e(V)

matrix Vmatrix = e(V)*realmse/e(rmse)^2

scalar list realmse

matrix list Vmatrix

ereturn post bmatrix Vmatrix, noclear

ereturn display