******************************************
*           -triprobit-
* This is triprobit.ado
* February, 25 2002, beta version
* Antoine Terracol
* terracol@univ-paris1.fr
*
* (Seemingly unrelated) trivariate probit
* 	using the GHK simulator
******************************************


cap prog drop triprobit
program define triprobit
version 7.0

/*if no arguments: redisplay results*/
		if "`0'"=="" {
				if e(cmd)!="triprobit" {
								display as error "results from triprobit not found"
								exit
								}
				else {
					di ""
					di in gr "trivariate probit, GHK simulator, " in ye e(draws) in gr " draws"
					ml display
					di in smcl in gr "{hline 78}"
					di in gr "rho12= " in ye %10.0g e(rho12)  in gr "  Std. Err.= " in ye %10.0g e(rho12_se) in gr "  z= " in gr %10.0g in ye e(rho12)/e(rho12_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho12)/e(rho12_se))))*2
					di in gr "rho13= " in ye %10.0g e(rho13)  in gr "  Std. Err.= " in ye %10.0g e(rho13_se) in gr "  z= " in gr %10.0g in ye e(rho13)/e(rho13_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho13)/e(rho13_se))))*2
					di in gr "rho23= " in ye %10.0g e(rho23)  in gr "  Std. Err.= " in ye %10.0g e(rho23_se) in gr "  z= " in gr %10.0g in ye e(rho23)/e(rho23_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho23)/e(rho23_se))))*2
					di in smcl in gr "{hline 78}"
					di in green "LR test of rho12=rho13=rho23=0: " _c
					di in green "chi2(" in ye "3" in gr ") = " in ye %10.0g e(chi2_c) _c
					di in green "  Prob > chi2 = " in ye %10.0g chiprob(3,e(chi2_c)) _c
					exit
					}
			}

/*else estimate*/
		else triprob `0'

end




cap prog drop triprob
program define triprob, eclass
version 7.0

/*Parsing syntax (hacked from biprobit.ado)*/

/* get first equation*/
		gettoken first : 0, match(paren)
		gettoken first 0:0, parse(" ,[") match(paren)
		local left "`0'"
		local junk: subinstr local first ":" ":"
		local first : subinstr local first "=" " "
		gettoken dep1 0: first, parse(" ,[")
		unab dep1: `dep1'
		confirm variable `dep1'
		capture assert `dep1'==1 | `dep1'==0
		if _rc==9 {
				di in red "depvar1 should be binary (0 or 1)"
				exit
				}
		syntax [varlist(default=none)] [, noCONstant]
		local ind1 `varlist'
		local nc1 `constant'

/* get second equation*/
		local 0 "`left'"
		gettoken second 0:0, parse(" ,[") match(paren)
		if "`paren'" != "(" {
			dis in red "three equations required"
			exit 110
		}
		local left "`0'"
		local junk : subinstr local second ":" ":"
		local second : subinstr local second "=" " "
		gettoken dep2 0: second, parse(" ,[")
		unab dep2: `dep2'
		confirm variable `dep2'
		capture assert `dep2'==1 | `dep2'==0
		if _rc==9 {
				di in red "depvar2 should be binary (0 or 1)"
				exit
				}
		syntax [varlist(default=none)] [, noCONstant]
		local ind2 `varlist'
		local nc2 `constant'

/* get third equation*/
		local 0 "`left'"
		gettoken third 0:0, parse(" ,[") match(paren)
		if "`paren'" != "(" {
			dis in red "three equations required"
			exit 110
		}
		local left "`0'"
		local junk : subinstr local third ":" ":"
		local third : subinstr local third "=" " "
		gettoken dep3 0: third, parse(" ,[")
		unab dep3: `dep3'
		confirm variable `dep3'
		capture assert `dep3'==1 | `dep3'==0
		if _rc==9 {
				di in red "depvar3 should be binary (0 or 1)"
				exit
				}
		syntax [varlist(default=none)] [, noCONstant]
		local ind3 `varlist'
		local nc3 `constant'

/* remaining options */
		local 0 "`left'"
		syntax [if] [in] [aw pw fw iw] [, DRaws(integer 25) Robust CLuster(varname) /*
							*/ SEED(integer 123456789)  noINIT MLOpts(string) /*
							*/ FROM(string) *]
		global D_Draws=`draws'
		global S_seed=`seed'
		local option0 `options'
		marksample touse
		markout `touse' `dep1' `dep2' `dep3' `ind1' `ind2' `ind3'
		if "`from'"!="" {
					local init="noinit"
					local initfrom="init `from'"
					}

/*Run comparison and initialization probits*/
			qui probit `dep1' `ind1' if `touse' [`weight' `exp'] ,`nc1'
			local ll1=e(ll)
			matrix I1=e(b)
			matrix coleq I1 = `dep1'
			qui probit `dep2' `ind2' if `touse' [`weight' `exp'] ,`nc2'
			local ll2=e(ll)
			matrix I2=e(b)
			matrix coleq I2 = `dep2'
			qui probit `dep3' `ind3' if `touse' [`weight' `exp'] ,`nc3'
			local ll3=e(ll)
			matrix I3=e(b)
			matrix coleq I3 = `dep3'
			matrix Init=I1,I2,I3
		      local llp=`ll1'+`ll2'+`ll3'


/*cluster option*/
			if "`cluster'"!="" {
						local kluster="cluster(`cluster')"
						}

/*display some info*/
			di ""
			di in gr "trivariate probit, GHK simulator, " in ye $D_Draws in gr " draws"
			di ""
			di in gr "Comparison log likelihood = " in ye %10.0g `llp'


/*Estimate the bastard*/
			ml model lf triprob_ml (`dep1':`dep1'=`ind1' ,`nc1') (`dep2':`dep2'=`ind2' ,`nc2') (`dep3':`dep3'=`ind3' ,`nc3') /*
							*/ / athrho12 / athrho13 / athrho23 if `touse' [`weight' `exp'], wald(-6) `robust' `kluster'
			if "`init'"=="" {
						ml init Init , noskip
						}
			if "`initfrom'"!="" {
						ml `initfrom'
						}
			ml max ,`options'


/*delta method to get rhos' Std. Errs*/
			local ath12=[athrho12][_cons]
			local athvar12=([athrho12]_se[_cons])^2
			local rho12=(exp(2*`ath12')-1)/(exp(2*`ath12')+1)
			local var12=[((4*exp(2*`ath12'))/((exp(2*`ath12')+1)^2))^2]*`athvar12'
			local se12=sqrt(`var12')

			local ath13=[athrho13][_cons]
			local athvar13=([athrho13]_se[_cons])^2
			local rho13=(exp(2*`ath13')-1)/(exp(2*`ath13')+1)
			local var13=[((4*exp(2*`ath13'))/((exp(2*`ath13')+1)^2))^2]*`athvar13'
			local se13=sqrt(`var13')

			local ath23=[athrho23][_cons]
			local athvar23=([athrho23]_se[_cons])^2
			local rho23=(exp(2*`ath23')-1)/(exp(2*`ath23')+1)
			local var23=[((4*exp(2*`ath23'))/((exp(2*`ath23')+1)^2))^2]*`athvar23'
			local se23=sqrt(`var23')

/*post results*/
			estimate scalar rho12 = `rho12'
			estimate scalar rho13 = `rho13'
			estimate scalar rho23 = `rho23'
			estimate scalar rho12_se = `se12'
			estimate scalar rho13_se = `se13'
			estimate scalar rho23_se = `se23'
			estimate scalar ll_c = `llp'
			est scalar chi2_c = abs(-2*(e(ll_c)-e(ll)))
			estimate local draws = `draws'
			estimate local seed = `seed'
			estimates local cmd  triprobit
			estimate local depvar "`dep1' `dep2' `dep3'"


/*display results*/
			di in smcl in gr "{hline 78}"
			di in gr "rho12= " in ye %10.0g e(rho12)  in gr "  Std. Err.= " in ye %10.0g e(rho12_se) in gr "  z= " in gr %10.0g in ye e(rho12)/e(rho12_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho12)/e(rho12_se))))*2
			di in gr "rho13= " in ye %10.0g e(rho13)  in gr "  Std. Err.= " in ye %10.0g e(rho13_se) in gr "  z= " in gr %10.0g in ye e(rho13)/e(rho13_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho13)/e(rho13_se))))*2
			di in gr "rho23= " in ye %10.0g e(rho23)  in gr "  Std. Err.= " in ye %10.0g e(rho23_se) in gr "  z= " in gr %10.0g in ye e(rho23)/e(rho23_se) in gr "  Pr>|z|= " in ye %10.0g (1-norm(abs(e(rho23)/e(rho23_se))))*2
			di in smcl in gr "{hline 78}"
			di in green "LR test of rho12=rho13=rho23=0: " _c
			di in green "chi2(" in ye "3" in gr ") = " in ye %10.0g e(chi2_c) _c
			di in green "  Prob > chi2 = " in ye %10.0g chiprob(3,e(chi2_c)) _c


/*how to get rid of these?*/
			macro drop S_seed D_draws


end