******************************************
* This is triprob_ml.ado
* February, 22 2002, beta version
* Antoine Terracol
* terracol@univ-paris1.fr
*
* -ml- code for -triprobit-
******************************************

cap prog drop triprob_ml
program define triprob_ml
version 7.0
args lnf theta1 theta2 theta3 theta4 theta5 theta6
******************
* theta1=X1Beta1
* theta2=X2Beta2
* theta3=X3Beta3
* theta4=rho12
* theta5=rho13
* theta6=rho23
*****************

tempvar rho12 rho13 rho23 triv1 triv2 

qui gen double `rho12'=[exp(2*`theta4')-1]/[exp(2*`theta4')+1]
qui gen double `rho13'=[exp(2*`theta5')-1]/[exp(2*`theta5')+1]
qui gen double `rho23'=[exp(2*`theta6')-1]/[exp(2*`theta6')+1]



********************
*1st trivariate CDF
********************
tempvar Prod B1 eps1 B2 eps2 B3

local R12=`rho12'
local R13=`rho13'
local R23=`rho23'
if (`R12'==. | `R13'==. | `R23'==.) {
						qui replace `lnf'=.
						exit
						}

matrix Sig = (1    , `R12' , `R13' \ /*
		*/`R12' ,   1   , `R23' \ /*
            */`R13' , `R23' ,  1    )


/*Check that the var-covar matrix is d.p.*/
matrix symeigen X V = Sig

if V[1,3] <= 0 {
		 qui gen double `triv1' =.
		 exit
		 }

matrix L=cholesky(Sig)

local l11=L[1,1]
local l21=L[2,1]
local l22=L[2,2]
local l33=L[3,3]
local l32=L[3,2]
local l31=L[3,1]



qui gen double `Prod'=0 
qui gen double  `B1'=(-`theta1'/`l11') 
qui gen double  `eps1' = 0 
qui gen double  `B2' = 0 
qui gen double  `eps2' = 0 
qui gen double  `B3' = 0 
set seed $S_seed
local repl=$D_Draws
local r=1
while `r'<=`repl' {
			qui replace `eps1'=(invnorm(normprob(`B1')*uniform())) 
	      	qui replace `B2'=[-`theta2'-(`l21'*`eps1') ]/`l22' 
			qui replace `eps2'=(invnorm(normprob(`B2')*uniform())) 
	      	qui replace `B3'=[-`theta3'-(`l31'*`eps1' )-(`l32'*`eps2')]/`l33' 
			qui replace `Prod'=`Prod'+[normprob(`B1')*normprob(`B2')*normprob(`B3')]
			local r=`r'+1
		}


qui gen double `triv1' = `Prod'/`repl' 

******************
*2nd trivariate CDF
*****************;

local R12=`rho12'
local R13=`rho13'
local R23=`rho23'
if (`R12'==. | `R13'==. | `R23'==.) {
						qui replace `lnf'=.
						exit
						}

matrix Sig = (1    , `R12' , `R13' \ /*
		*/`R12' ,   1   , `R23' \ /*
            */`R13' , `R23' ,  1    )

/*Check that the var-covar matrix is d.p.*/
matrix symeigen X V = Sig
if V[1,3] <= 0 {
		 	qui gen double `triv2' =.
			exit
			}
matrix L=cholesky(Sig)

local l11=L[1,1]
local l21=L[2,1]
local l22=L[2,2]
local l33=L[3,3]
local l32=L[3,2]
local l31=L[3,1]



qui replace `Prod'=0 
qui replace  `B1'=(`theta1'/`l11') 
qui replace  `eps1' = 0 
qui replace  `B2' = 0 
qui replace  `eps2' = 0 
qui replace  `B3' = 0 

local r=1
while `r'<=`repl' {			
			qui replace `eps1'=(invnorm(normprob(`B1')*uniform())) 
	      	qui replace `B2'=[`theta2'-(`l21'*`eps1') ]/`l22' 
			qui replace `eps2'=(invnorm(normprob(`B2')*uniform())) 
	      	qui replace `B3'=[`theta3'-(`l31'*`eps1' )-(`l32'*`eps2')]/`l33' 
			qui replace `Prod'=`Prod'+[normprob(`B1')*normprob(`B2')*normprob(`B3')]
			local r=`r'+1
		}



qui gen double `triv2' = `Prod'/`repl' 

******************************************************************************
qui replace `lnf'=ln(`triv1') if $ML_y1==0 & $ML_y2==0 & $ML_y3==0
qui replace `lnf'=ln(binorm(-`theta1',-`theta2',`rho12')-`triv1') if $ML_y1==0 & $ML_y2==0 & $ML_y3==1
qui replace `lnf'=ln(binorm(-`theta1',-`theta3',`rho13')-`triv1') if $ML_y1==0 & $ML_y2==1 & $ML_y3==0
qui replace `lnf'=ln(binorm(`theta2',`theta3',`rho23')-`triv2') if $ML_y1==0 & $ML_y2==1 & $ML_y3==1
qui replace `lnf'=ln(binorm(-`theta2',-`theta3',`rho23')-`triv1') if $ML_y1==1 & $ML_y2==0 & $ML_y3==0
qui replace `lnf'=ln(binorm(`theta1',`theta3',`rho13')-`triv2') if $ML_y1==1 & $ML_y2==0 & $ML_y3==1
qui replace `lnf'=ln(binorm(`theta1',`theta2',`rho12')-`triv2') if $ML_y1==1 & $ML_y2==1 & $ML_y3==0
qui replace `lnf'=ln(`triv2') if $ML_y1==1 & $ML_y2==1 & $ML_y3==1


end