*! version 1.2.0  07oct1996
program define biprob
	version 5.0
        local options "Level(integer $S_level)"
        if (substr("`1'",1,1)=="," | "`*'"=="") {
                if "$S_E_cmd"~="biprob" {
                        error 301
                }
                parse "`*'"
                if "`level'" != "" {
                        if `level' < 10 | `level' > 99 {
                                local level = 95
                        }
                }
        }
        else {
		local varlist   "req ex min(2)"
		local if	"opt"
		local in	"opt"
		local options	"Robust SCore(string) TRace noLOg"
		local options	"`options' Level(integer $S_level)"
		local options	"`options' noCHI CLuster(string) *"
		parse "`*'"

		if "`score'" != "" {
			if "`robust'" == "" {
				local robust "robust"
			}
			local n : word count `score'
			if `n' != 2 {
				noi di in red "need 2 arguments in score()"
				exit 198
			}
			local sc1nam : word 1 of `score'
			local sc2nam : word 2 of `score'
			confirm new var `sc1nam'
			confirm new var `sc2nam'
		}
		parse "`varlist'", parse(" ")
		local dep1 "`1'"
		local dep2 "`2'"
		mac shift 2
		local ind1 "`*'"
		local ind2 "`*'"
		eq `dep1' `ind1'
		eq `dep2' `ind2'
		local eq1 "`dep1'"
		local eq2 "`dep2'"
		
		quietly {
			tempvar touse
			mark `touse' `if' `in' 
			markout `touse' `dep1' `ind1' `dep2' 

			if "`cluster'" != "" {
				unabbrev `cluster'
				if $S_2 > 1 {
					noi di in red /*
				*/ "only one varname allowed in cluster()"
					exit 198
				}
				local cluster "$S_1"
				local clarg ", cluster(`cluster')"
				local robust "robust"
				confirm var `cluster'
				markout `touse' `cluster', strok
			}

			global SC_1
			global SC_2
			global SC_3
			global S_y1 "`dep1'"
			global S_y2 "`dep2'"
			global S_X1 "`ind1'"
			global S_X2 "`ind2'"

			if "`chi'" == "" {
				/* Fit constant only model */
				probit `dep1' if `touse'
				tempname b0 b1 b2 b3 eq01 eq02
				mat `b1' = get(_b)
				mat coleq `b1' = `eq01'
				probit `dep2' if `touse'
				mat `b2' = get(_b)
				mat coleq `b2' = `eq02'
				mat `b3' = (0)
				mat colnames `b3' = rho:_cons
				mat `b0' = `b1',`b2'
				mat `b0' = `b0',`b3'
				eq rho:
				eq `eq01' : `dep1'
				eq `eq02' : `dep2'

				ml begin
				tempvar mltouse 
				ml sample `mltouse' if `touse', noauto
				ml function bipr_l2
				ml method deriv2
				ml model b = `eq01' `eq02' rho, /*
					*/ depv(110) from(`b0')
				if "`log'" == "nolog" {
					local prel "qui"
				}
				else    local prel "noi"
				`prel' di 
				`prel' di in gr "Fitting constant only model"
				`prel' ml max f V, `trace' ltol(1e-4)
				local ll0 = scalar(f)
				eq drop `eq01'
				eq drop `eq02'
				local llarg "lf0(`ll0') pr2"
			}
			else {
				local llarg ""
			}

			/* Fit full model */
		
			probit `dep1' `ind1' if `touse'
			tempname b0 b1 b2 b3
			mat `b1' = get(_b)
			mat coleq `b1' = `dep1'
			probit `dep2' `ind2'  if `touse'
			mat `b2' = get(_b)
			mat coleq `b2' = `dep2'
			mat `b3' = (0)
			mat colnames `b3' = rho:_cons
			mat `b0' = `b1',`b2'
			mat `b0' = `b0',`b3'
			eq rho:

			if "`robust'" != "" {
				tempvar s1 s2 s3
				gen double `s1' = .
				gen double `s2' = .
				gen double `s3' = .
				global SC_1 "`s1'"
				global SC_2 "`s2'"
				global SC_3 "`s3'"
			}
			ml begin
			tempvar mltouse 
			ml sample `mltouse' if `touse', noauto
			ml function bipr_l2
			ml method deriv2
			ml model b = `eq1' `eq2' rho, depv(110) from(`b0')
			if "`log'" == "nolog" {
				local prel "qui"
			}
			else    local prel "noi"
			`prel' di 
			`prel' di in green "Fitting full model"
			`prel' ml max f V, `trace' `options'
			local title "Bivariate probit regression"
			ml post biprob, `llarg' title(`title')
			if "`robust'" != "" {
				_robust `s1' `s2' `s3' `clarg'
				if "`score'" != "" {
					rename `s1' `sc1nam'
					rename `s2' `sc2nam'
				}
			}
		}
		if "`cluster'" != "" {
			global S_E_cvn "`cluster'"
		}
		global S_E_cmd "biprob"
	}
	ml mlout biprob, level(`level') 
end