Jonah - I don't have a link to these - but if you wanted to put them on
your site and "pretend" they were elsewhere, i.e. direct questions to me,
that'd be great (or just load this email as is if you want).

There are three files: probiv.do, probiv2.do (which has a "transformed"
rho),
and probivdata.do, that generates test data.  I have not implemented these
as stata commands: in each file, you need to enter the name of the dataset,
and then the names of the left hand side and right hand side variables for
each equation as local macros.  After that things should be "automatic".

later,
Deon



---------------------------- probiv.do starts here --------------------
capture log close
clear
set more 1
set matsize 300

log using probiv,replace

*--------------------------------------------------------*
* Writted by Deon Filmer                                 *
* Version 1 - Nov 8 1999                                 *
*                                                        *
* Following estimates the structure:                     *
* y1 = a*y2 + b *X        + eps1   (probit equation)     *
* y2 =        c1*X + c2*Z + eps2   (continuous)          *
* where (eps1,eps2) ~ jointly normal with correlation    *
* rho, eps1 normalized to have variance 1, eps2 has      *
* standard deviation sigma2                              *
*--------------------------------------------------------*

**********************************************************
*** Input dataset to be used *****************************
**********************************************************

use probivdata

**********************************************************
*** Input left hand side and right hand side variables ***
**********************************************************

* probit equation - replace contents of "" with
* lhs variable for pbdep1 and  rhs variables for pbind1
local pbdep1 "y1"
local pbind1 "y2 x1 x2"

* continuous equation 1- replace contents of "" with
* lhs variable for ctdep1 and  rhs variables for ctind1
local ctdep1 "y2"
local ctind1 "x1 x2 z1 z2"

**********************************************************
*****  You shouldn't need to change anything below  ******
**********************************************************

probit  `pbdep1'  `pbind1'
mat Mpb1 = e(b)

regress `ctdep1'  `ctind1'
mat Mct1    = e(b)
predict _alph2 , resid
local rmse2 = _result(9)


***********************************************************
* Rivers and Vuong test for exogeneity

probit  `pbdep1'  `pbind1' _alph2
test _alph2
if chiprob(1,_result(6)) < 0.05 {di "Reject exogeneity"}
if chiprob(1,_result(6)) >=0.05 {di "Cannot reject exogeneity"}

***********************************************************


capture program drop pbcont

program define pbcont
version 6
args lnf I1 I2 rho12 sigma2

tempvar eps2 index cindex lnf1 lnf2

qui gen double `eps2'  = $ML_y2 - `I2'
qui gen double `index' =  ( `I1' + `eps2'*`rho12'/`sigma2')
qui replace    `index' = `index' / sqrt(1-`rho12'^2)
qui gen double `cindex'= normprob(`index')
qui gen double `lnf1'  = .
qui replace    `lnf1'  = ln(`cindex')   if $ML_y1==1
qui replace    `lnf1'  = ln(1-`cindex') if $ML_y1==0
qui gen double `lnf2'  = -.5 * ln(2.0*_pi) -.5*ln(`sigma2'^2)
qui replace    `lnf2'  = `lnf2' -  .5*(`eps2'/`sigma2')^2.0
qui replace `lnf' = `lnf1' + `lnf2'

end

#delimit;
ml model lf pbcont
        (`pbdep1' = `pbind1')
        (`ctdep1' = `ctind1')
        /rho12
        /sigma2
;
#delimit cr

ml init Mpb1 Mct1 0 `rmse2' , copy
ml maximize

log close

--------------------- probiv.do ends here -------------------------------------


--------------------- probiv2.do starts here ----------------------------------
--------------------- (same as probiv.do but uses a transformation of rho) ----
capture log close
clear
set more 1
set matsize 300

log using probiv2,replace

*--------------------------------------------------------*
* Writted by Deon Filmer                                 *
* Version 1 - Nov 8 1999                                 *
*                                                        *
* Following estimates the structure:                     *
* y1 = a*y2 + b *X        + eps1   (probit equation)     *
* y2 =        c1*X + c2*Z + eps2   (continuous)          *
* where (eps1,eps2) ~ jointly normal with correlation    *
* rho, eps1 normalized to have variance 1, eps2 has      *
* standard deviation sigma2                              *
*--------------------------------------------------------*

**********************************************************
*** Input dataset to be used *****************************
**********************************************************

use probivdata

**********************************************************
*** Input left hand side and right hand side variables ***
**********************************************************

* probit equation - replace contents of "" with
* lhs variable for pbdep1 and  rhs variables for pbind1
local pbdep1 "y1"
local pbind1 "y2 x1 x2"

* continuous equation 1- replace contents of "" with
* lhs variable for ctdep1 and  rhs variables for ctind1
local ctdep1 "y2"
local ctind1 "x1 x2 z1 z2"

**********************************************************
*****  You shouldn't need to change anything below  ******
**********************************************************

probit  `pbdep1'  `pbind1'
mat Mpb1 = e(b)

regress `ctdep1'  `ctind1'
mat Mct1    = e(b)
predict _alph2 , resid
local rmse2 = _result(9)


***********************************************************
* Rivers and Vuong test for exogeneity

probit  `pbdep1'  `pbind1' _alph2
test _alph2
if chiprob(1,_result(6)) < 0.05 {di "Reject exogeneity"}
if chiprob(1,_result(6)) >=0.05 {di "Cannot reject exogeneity"}

***********************************************************


capture program drop pbcont

program define pbcont
version 6
args lnf I1 I2 trho12 sigma2

tempvar eps2 index cindex lnf1 lnf2 rho12

qui gen double `rho12' = (exp(2*`trho12')-1)/(exp(2*`trho12')+1)
qui gen double `eps2'  = $ML_y2 - `I2'
qui gen double `index' =  ( `I1' + `eps2'*`rho12'/`sigma2')
qui replace    `index' = `index' / sqrt(1-`rho12'^2)
qui gen double `cindex'= normprob(`index')
qui gen double `lnf1'  = .
qui replace    `lnf1'  = ln(`cindex')   if $ML_y1==1
qui replace    `lnf1'  = ln(1-`cindex') if $ML_y1==0
qui gen double `lnf2'  = -.5 * ln(2.0*_pi) -.5*ln(`sigma2'^2)
qui replace    `lnf2'  = `lnf2' -  .5*(`eps2'/`sigma2')^2.0
qui replace `lnf' = `lnf1' + `lnf2'

end

#delimit;
ml model lf pbcont
        (`pbdep1' = `pbind1')
        (`ctdep1' = `ctind1')
        /trho12
        /sigma2
;
#delimit cr

ml init Mpb1 Mct1 0 `rmse2' , copy
ml maximize

di (exp(2* [trho12]_cons)-1)/(exp(2*[trho12]_cons)+1)

log close
---------------------------- probiv2.do ends here -------------------------------


---------------------------- probivdata.do starts here --------------------------
---------------------------- (creates data to test probiv.do and probiv2.do) ----
capture log close
clear set more 1

log using probivdata , replace


clear

* following generates data for the the structure:
* y1 = a*y2 + b*x + eps1   (probit equation)
* y2 =        c*z + eps2   (continuous)


set seed 102167
set obs 1000

gen eps1 = invnorm(uniform())
gen eps2 = eps1 + 2*invnorm(uniform())

gen x1=uniform()
gen x2=uniform()

gen z1=uniform()
gen z2=uniform()

gen y2  =      x1 + 2 * x2  + z1 - z2 + eps2

gen y1s = y2 - x1 +     x2            + eps1

gen y1  = y1s>=0

save probivdata,replace

log close

---------------------------------- probivdata.do ends here ----------------