? Murphy & Topell procedure with an example 

? MROZ.dat, Wooldridge, p. 544


Read  ; Nobs = 753 
	; Nvar = 22 
	; file = mroz.dat
      ; names = 3 $

Namelist ; X1 = nwifeinc, educ, exper, expersq, age$
namelist ; xx1= nwifeinc, educ, expersq, age $

? TOBIT WITH SAMPLE SELECTION
 Probit ; lhs = inlf
	  ; rhs = one, x1, hushrs, husage
	  ; hold $
 Select ; Tobit ; MLE 
	  ; lhs = hours
	  ; rhs = one, x1 $

?=====================================================================
? Two step estimation of a probit model when one of the RHS variables
? is the prediction from a first step linear regression model.
? Computations are based on Theorem 4.19 in Greene (1997, p. 142) and
? Murphy and Topel (1985).
?
? User must define the following:
? NAMELIST ; X = RHS in the probit model, not including fitted w
?          ; Z = RHS in the regression $
? CREATE   ; y = LHS in the probit model

?          ; w = LHS in the regression $
?=====================================================================
?
? Fit the regression and keep results.  Use ML rather than degrees of
? freedom corrected least squares results.  Need residuals, preciction,
? asymptotic covariance matrix and disturbance variance estimator.
? V1 is (e'e/N) * <X'X>

NAMELIST ; X = one, x1, hushrs, husage
         ; Z = one, x1 $
CREATE   ; y = inlf
         ; w = hours $

  REGRESS  ; Lhs = w ; Rhs = Z ; Keep = wfit ; Res = ew $
  MATRIX   ; V1  = {(N-Kreg)/N} * Varb $ (Use ML)
  CALC     ; sigmasq = sumsqdev/N $

?
? Now, fit probit model and correct asymptotic covariance matrix.
? Lambda is used to compute derivatives of log-likelihood.  Note
? dLog-L/dbeta = X'Lambda.
?
  NAMELIST ; XW = wfit,xx1,hushrs,husage,one $
  PROBIT   ; Lhs = y ; Rhs = XW ; Hold(IMR = lambda) $
?
? Two matrices C and R are of the form XW'[variable]Z where variable
? is computed from derivatives of log-likelihoods.  The V2 matrix in
? the computation is the asymptotic covariance matrix from the probit
? model, before the correction.

?
  CREATE   ; vc = B(1) * lambda^2 
           ; vr = lambda * ew/sigmasq $
  MATRIX   ; C = XW'[vc]Z
           ; R = XW'[vr]Z 
           ; V = C * V1 * C' - R * V1 * C' - C * V1 * R'
           ; V = Varb + Varb * V * Varb
           ; stat(B,V) $

  matrix; list; v $

? We get the same results given by Stata.