*! version 2.6.1  03oct2000 
program define biprobit2, eclass byable(onecall)
	version 6.0
	if replay() {
		if "`e(cmd)'" != "biprobit" {
			noi di in red "results for biprobit not found"
			exit 301
		}
		if _by() { error 190 } 
		Display `0'
		exit `rc'
	}
	if _by() {
		by `_byvars'`_byrc0': Estimate `0'
	}
	else	Estimate `0'
end

program define Estimate, eclass byable(recall)

	/* two syntax, handle one at a time */

	gettoken first : 0, match(paren)
	
	if "`paren'" == "" {
		/* syntax 1, bivariate probit models */

		gettoken dep1 0:0, parse (" =,[")
		unab dep1 : `dep1'
		confirm variable `dep1'
		gettoken dep2 0:0, parse (" =,[")
		unab dep2 : `dep2'
		confirm variable `dep2'
		gettoken junk left :0, parse ("=")
		if "`junk'" == "=" {
			local 0 "`left'"
		}
		syntax [varlist(default=none)] [if] [in] [pw fw iw] /*
			*/[, Robust CLuster(varname) SCore(string)  /*
			*/ offset1(varname) offset2(varname) PARtial /*
			*/ noCONstant noSKIP Level(integer $S_level) /* 
			*/ noLog MLOpts(string) FROM(string) 	/*
			*/ CIter(integer 16000) * ]  

		if _by() {
			_byoptnotallowed score() `"`score'"'
		}

		local ind1 `varlist'
		local ind2 `varlist'
		local title "Bivariate probit regression"
                local dep1n "`dep1'"
                local dep2n "`dep2'"
                local e_neq = 2
		local nc1 `constant'
		local nc2 `constant'
		local option0 `options'
		marksample touse
		markout `touse' `dep1' `dep2' `offset1' `offset2', strok
	}
	else {
		/* syntax 2, seemingly unrelated bivariate probit model */
		

					/* get first equation */
		gettoken first 0:0, parse(" ,[") match(paren)
		local left "`0'"
		local junk: subinstr local first ":" ":", count(local number)
		if "`number'" == "1" {
			gettoken dep1n first: first, parse(":")
			gettoken junk first: first, parse(":")
		}
		local first : subinstr local first "=" " "
		gettoken dep1 0: first, parse(" ,[") 
		unab dep1: `dep1' 
		confirm variable `dep1'
		if "`dep1n'" == "" {
			local dep1n "`dep1'"
		}
		syntax [varlist(default=none)] [, OFFset(varname) noCONstant]
		local ind1 `varlist'
		local offset1 `offset' 
		local nc1 `constant'

					/* get second equation */
		local 0 "`left'"
		gettoken second 0:0, parse(" ,[") match(paren)
		if "`paren'" != "(" {
			dis in red "two equations required"
			exit 110
		}
		local left "`0'"
		local junk : subinstr local second ":" ":", count(local number)
		if "`number'" == "1" {
			gettoken dep2n second: second, parse(":")
			gettoken junk second: second, parse(":")
		}
		local second : subinstr local second "=" " "
		gettoken dep2 0: second, parse(" ,[") 
		unab dep2: `dep2' 
		confirm variable `dep2'
		if "`dep2n'" == "" {
			local dep2n "`dep2'"
		}
		syntax [varlist(default=none)] [, OFFset(varname) noCONstant]
		local ind2 `varlist'
		local offset2 `offset' 
		local nc2 `constant'

					/* remain options */
		local 0 "`left'"
		syntax [if] [in] [pw fw iw] [, Robust Cluster(varname)    /*
			*/ SCore(string)  /*
			*/ PARtial noSKIP Level(integer $S_level)  	   /*
			*/ noLOG MLOpts(string) FROM(string) * ]

		local title "Seemingly unrelated bivariate probit"
		local option0 `options'
		marksample touse
		markout `touse' `dep1' `dep2' `ind1' `ind2' `offset1'	  /*
          		*/ `offset2'
	}
	
					
	if "`partial'" != "" {
		local title "Partial observability bivariate probit"
	}        
	
	local wtype `weight'
        local wtexp `"`exp'"'
        if "`weight'" != "" { local wgt `"[`weight'`exp']"'  }

        if "`weight'" == "pweight" | "`cluster'" != "" {
                local robust "robust"
        }
	if "`cluster'" ! = "" { local clopt "cluster(`cluster')" }
	mlopts stdopts, `option0'

					/* test collinearity */
	_rmcoll `ind1' `wgt' if `touse', `nc1'
	local ind1 "`r(varlist)'"
	_rmcoll `ind2' `wgt' if `touse', `nc2'
	local ind2 "`r(varlist)'"

	if "`level'" != "" {
		local level "level(`level')"
	}

	local scvar `"`score'"'
	if "`scvar'" != "" {
                local scvar1 : word 1 of `scvar'
                local scvar2 : word 2 of `scvar'
                local scvar3 : word 3 of `scvar'
                local scvar4 : word 4 of `scvar'
                if "`scvar4'" != "" {
    		noi di in red "may only request up to 3 score variables"
             	   exit 198 
                }
                tempvar sc1 sc2 sc3
                local scopt "score(`sc1' `sc2' `sc3')"
        }

	if "`offset1'" != "" { local offo1 "offset(`offset1')" }
        if "`offset2'" != "" { local offo2 "offset(`offset2')" }


	local skip = cond("`skip'"!="","","skip")

	qui {
		if "`log'" == "" {
                        local log "noisily"
                }
                else    local log "quietly"

		count if `touse' 
		local N = r(N)
		count if `dep1'==0 & `touse'
		local d10 = r(N)
		count if `dep2'==0 & `touse'
		local d20 = r(N)
		if `d10' == 0 {
			di in red "`dep1' is never zero"
			exit 2000
		}
		else if `d10' == `N' {
			di in red "`dep1' is always zero"
			exit 2000
		}
		else if `d20' == 0 {
			di in red "`dep2' is never zero"
			exit 2000
		}
		else if `d20' == `N' {
			di in red "`dep2' is always zero"
			exit 2000
		}


		if "`partial'" != "" {
			tempvar p1
			gen `p1' = (`dep1'!=0)*(`dep2'!=0)
			local ddd1 "`p1'"
			local ddd2 "`p1'"
		}
		else {
			local ddd1 "`dep1'"
			local ddd2 "`dep2'"
		}


		if "`from'" == "" {
			local log0 `log'
			if "`partial'" != "" {
				local log qui
			}
			`log' di in gr _n /*
		        */	"Fitting comparison equation 1:" 
			capture `log' probit `ddd1' `ind1' `wgt' if `touse', /*
				*/`offo1'   `nc1' asis nocoef iter(`citer')
			if _rc == 0 {
				tempname b1
				mat `b1' = get(_b)
				local ll1 = e(ll)
				mat coleq `b1' = `dep1n'
			}
			`log' di in gr _n /*
			*/ "Fitting comparison equation 2:" 
			capture `log' probit `ddd2' `ind2' `wgt' if `touse', /*
				*/`offo2'   `nc2' asis nocoef iter(`citer')
			if _rc == 0 {
				tempname b2
				mat `b2' = get(_b)
				local ll2 = e(ll)
				mat coleq `b2' = `dep2n'
			}
			local llp = `ll1' + `ll2'
			local log `log0'	

			if "`b1'`b2'" != "" {
				tempname from
				if "`b1'" != "" & "`b2'" != "" { 
					mat `from' = `b1' , `b2'
					if "`partial'" == "" {
					`log' di in gr _n /*
					*/ "Comparison:    log likelihood = " /*
                                        */ in ye %10.0g `llp'
					}
				}
				else if "`b1'" != "" {
					mat `from' = `b1'
					local lrtest "nolrtest"
				}
				else {
					mat `from' = `b2'
					local lrtest "nolrtest"
				}
				local getvals 1
			}
			if "`partial'" != "" {
				`log' di in gr _n /*
				*/	"Fitting comparison model: "
				#delimit ;
				`log' ml model lf bip0_lf  	
				   (`dep1n' : `dep1' = `ind1', `nc1' `offo1')
				   (`dep2n' : `dep2' = `ind2', `nc2' `offo2')
				   if `touse' `wgt', collinear missing max
				   nooutput nopreserve init(`from');
				#delimit cr
				local llp = e(ll)
			}	

		}

		if "`nc1'`nc2'" != ""   { local skip "skip" }
		if "`ind1'`ind2'" == "" { local skip "skip" }
		if "`robust'" != ""     { local skip "skip" }

		tempname a0
		if "`skip'" == "" {
			tempname a from0
			capture probit `ddd1' if `touse' `wgt', `offo1'   asis
			if _rc == 0 {
				tempname cb1
				mat `cb1' = get(_b)
				mat coleq `cb1' = `dep1n'
			}
			capture probit `ddd2' if `touse' `wgt', `offo2'   asis
			if _rc == 0 {
				tempname cb2
				mat `cb2' = get(_b)
				mat coleq `cb2' = `dep2n'
			}
			mat `a' = (-.3)
			mat colnames `a' = athrho:_cons
			if "`cb1'" != "" & "`cb2'" != "" {
				mat `from0' = `cb1',`cb2'
			}
			else if "`cb1'" != "" { mat `from0' = `cb1' }
			else if "`cb2'" != "" { mat `from0' = `cb2' }
			mat `from0' = `from0',`a'

			`log' di in gr _n "Fitting constant-only model:"
			if "`partial'" != "" {
				#delimit ;
				`log' ml model d2 bipp_lf 
					(`dep1n': `dep1' = , `nc' `offo1'  )
					(`dep2n': `dep2' = , `nc' `offo2'  )
					/athrho	
					if `touse' `wgt', collinear missing 
					max nooutput nopreserve wald(0) 
					init(`from0') search(off) `lrtest'
					`mlopts' `stdopts' `level' ;

				#delimit cr
			}
			else {
				#delimit ;
				`log' ml model d2 bipr_lf 
					(`dep1n': `dep1' = , `nc' `offo1'    )
					(`dep2n': `dep2' = , `nc' `offo2'  )
					/athrho	
					if `touse' `wgt', collinear missing 
					max nooutput nopreserve wald(0) 
					init(`from0') search(off) `lrtest'
					`level' `mlopts' `stdopts' ;
				#delimit cr
			}
			local cont "continue"
			if "`getvals'" != "" {
				mat `a0' = get(_b)
				mat `a0' = `a0'[1,3]
				mat colnames `a0' = athrho:_cons
				mat `from' = `from',`a0'
			}
		}
		else {
			if "`getvals'" != "" {
				mat `a0' = (0)
				mat colnames `a0' = athrho:_cons
				mat `from' = `from',`a0'
			}
			local cont "wald(2)"
		}
	}

	`log' di in gr _n "Fitting full model:"
	if "`partial'" != "" {
		#delimit ;
		`log' ml model d2 bipp_lf 
			(`dep1n': `dep1' = `ind1', `nc1' `offo1'  )
			(`dep2n': `dep2' = `ind2', `nc2' `offo2'  )
			/athrho	
			if `touse' `wgt', collinear missing max nooutput 
			nopreserve `cont' title(`title') `scopt' `robust' 
			`clopt' init(`from') search(off) `lrtest' 
			`level' `mlopts' `ll0' `stdopts';
		#delimit cr
	}
	else {
		#delimit ;
		`log' ml model d2 bipr_lf 
			(`dep1n': `dep1' = `ind1', `nc1' `offo1'  )
			(`dep2n': `dep2' = `ind2', `nc2' `offo2'  )
			/athrho	
			if `touse' `wgt', collinear missing max nooutput 
			nopreserve `cont' title(`title') `scopt' `robust' 
			`clopt' init(`from') search(off) `lrtest' 
			`level' `mlopts' `ll0' `stdopts';
		#delimit cr
	}

	local r = [athrho]_b[_cons]
	est scalar rho = (exp(2*`r')-1)/(exp(2*`r')+1)
        if "`scvar'" != "" { 
		rename `sc1' `scvar1' 
		if "`scvar2'" != "" { rename `sc2' `scvar2' }
		if "`scvar3'" != "" { rename `sc3' `scvar3' }
	}

	if "`llp'" != "" & "`robust'`lrtest'" == "" {
		est scalar ll_c = `llp'
		est scalar chi2_c = abs(-2*(e(ll_c)-e(ll)))
		est local chi2_ct "LR"
	}
	else {
		test [athrho]_cons
		est scalar chi2_c = r(chi2)
		est local chi2_ct "Wald"
	}
	est local predict "bipr_p"
        est local cmd "biprobit"

        Display ,`dispopt'
	exit `e(rc)'
end

program define Display
	syntax [,Level(int $S_level)]

	ml di, level(`level') neq(2) plus
	DispAl `level'
	DispLr
end

program define DispAl
	local level = `1'
	_diparm athrho, level(`level')
        di in smcl in gr "{hline 13}{c +}{hline 64}"
	_diparm athrho, level(`level') tanh label("rho")
	di in smcl in gr "{hline 13}{c BT}{hline 64}"
end

program define DispLr
	if "`e(ll_c)'" == "" { 
		di in green "Wald test of rho=0:                 " _c
		di in green "chi2(" in ye "1" in gr ") = " /*
			*/ in ye %8.0g e(chi2_c) _c
		di in green "    Prob > chi2 = " in ye %6.4f /*
			*/ chiprob(1,e(chi2_c))
		exit 
	}
	di in green "Likelihood ratio test of rho=0:     " _c
        di in green "chi2(" in ye "1" in gr ") = " /*
                */ in ye %8.0g e(chi2_c) _c
        di in green "    Prob > chi2 = " in ye %6.4f /*
                */ chiprob(1,e(chi2_c))
end