/**************************************************************************************
Stata script for "Regression Models for Binary Dependent Variables" and
				 http://www.indiana.edu/~statmath/stat/all/cdvm/index.html

Created on 09/17/2005
Last Modified on September 15, 2009

Written by Hun Myoung Park, Ph.D. 
	Statistical Software Analyst and Programmer
	UITS Center for Statistical and Mathematical Computing, Indiana University
	535 West Michigan Street IT417 Indianapolis, IN 46202
	317-274-0573 kucc625@indiana.edu http://www.masil.org

University Information Technology Services, Indiana University
Center for Statistical and Mathematical Computing
http://www.indiana.edu/~statmath

First, change the default directory 
**************************************************************************************/

cd c:\temp\cdvm

net from http://www.indiana.edu/~jslsoc/stata/
net install spost9_ado, replace
net get spost9_do, replace



use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear

sum trust belief educate income age male www, sep(10)
tab trust male, miss
tab trust www, miss
tab male www, miss
tab belief male, miss
tab belief www, miss





// Binary logit model 
use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear
logistic trust educate income age male www

logit

logit trust educate income age male www

predict resid, residual

test educate

mfx, dydx at(mean educate=16 male=0 www=1)

fitstat

di 2*(-326.755 - (-357.926))


// factor changes in the odds
listcoef, help 

listcoef, reverse

listcoef, percent help


// predicted probability 
prtab male www, x(educate=16 male=0 www=1) rest(mean)

prvalue, x(educate=16 male=0 www=1) rest(mean)


// marginal effects and discrete changes
prchange, x(educate=16 male=0 www=1) rest(mean)



// Draw a plot of predicted probability
// use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear

quietly logit trust educate income age male www

global points = 20
global group1 = "Logit_ed11"
global group2 = "Logit_ed10"
global group3 = "Logit_ed01"
global group4 = "Logit_ed00"

prgen educate, from(0) to(20) ncases($points) x(male=1 www=1) rest(mean) gen($group1)
prgen educate, from(0) to(20) ncases($points) x(male=1 www=0) rest(mean) gen($group2)
prgen educate, from(0) to(20) ncases($points) x(male=0 www=1) rest(mean) gen($group3)
prgen educate, from(0) to(20) ncases($points) x(male=0 www=0) rest(mean) gen($group4)

rename ${group1}x x

gen bygroup = 1 if x !=.
replace bygroup = 0 if (_n > $points) & (x== .) 
replace x = x[_n-$points] if _n > $points

forvalues v = 1/4 {
	gen cat`v' = ${group`v'}p1
}
drop if _n > $points*2

gen y1 = cat1
gen y2 = cat3
replace y1=cat2[_n-$points] if y1==.
replace y2=cat4[_n-$points] if y2==.

label define bygroup 1 "WWW Users" 0 "WWW Non-users"
label values bygroup bygroup

label var y1 "Men"
label var y2 "Women"
label var x "Education (Years)"

global min = 0
global max = 1
global step = .2
twoway ///
	(connected y1 x, lcolor(maroon) lpattern(solid) lwidth(medthick) connect(direct) msymbol(dh) mcolor(maroon)  ) ///
	(connected y2 x, lcolor(blue) lpattern(solid) lwidth(medium) connect(direct) msymbol(X) mcolor(blue) ) if _n<=40, ///
	by(bygroup) xscale(range(0 20)) xlabel(1(2)20, grid) yscale(range($min $max)) yline(.5) ylabel($min($step)$max, grid) ///    
    legend(region(lpattern(blank))) xtitle("Education (Years)") ytitle("Predicted Probabilities") 	








// Binary probit model 

use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear

probit trust educate income age male www

di _pi/sqrt(3)*.0907207

di 1.7*.0907207

probit trust educate income age male www, robust

fitstat


// standardized effects
listcoef, help


// marginal effects and discrete changes
mfx, at(mean educate=16 male=0 www=1)

prchange, x(educate=16 male=0 www=1) rest(mean)


// predicted probabilities
prtab male www, x(educate=16 male=0 www=1) rest(mean)

prvalue, x(educate=16 male=0 www=1) rest(mean)


// Draw a plot of predicted probabilities
//use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear

probit trust educate income age male www

global points = 20
global group1 = "Probit_ed11"
global group2 = "Probit_ed10"
global group3 = "Probit_ed01"
global group4 = "Probit_ed00"

prgen educate, from(0) to(20) ncases($points) x(male=1 www=1) rest(mean) gen($group1)
prgen educate, from(0) to(20) ncases($points) x(male=1 www=0) rest(mean) gen($group2)
prgen educate, from(0) to(20) ncases($points) x(male=0 www=1) rest(mean) gen($group3)
prgen educate, from(0) to(20) ncases($points) x(male=0 www=0) rest(mean) gen($group4)

rename ${group1}x x

gen bygroup = 1 if x !=.
replace bygroup = 0 if (_n > $points) & (x== .) 
replace x = x[_n-$points] if _n > $points

forvalues v = 1/4 {
	gen cat`v' = ${group`v'}p1
//	replace cat`v'=${group2}s`v'[_n-$points] if cat`v'==.
}
drop if _n > $points*2

gen y1 = cat1
gen y2 = cat3
replace y1=cat2[_n-$points] if y1==.
replace y2=cat4[_n-$points] if y2==.

label define bygroup 1 "WWW Users" 0 "WWW Non-users"
label values bygroup bygroup

label var y1 "Men"
label var y2 "Women"
label var x "Education (Years)"

global min = 0
global max = 1
global step = .2
twoway ///
	(connected y1 x, lcolor(maroon) lpattern(solid) lwidth(medthick) connect(direct) msymbol(dh) mcolor(maroon)  ) ///
	(connected y2 x, lcolor(blue) lpattern(solid) lwidth(medium) connect(direct) msymbol(X) mcolor(blue) ) if _n<=40, ///
	by(bygroup) xscale(range(0 20)) xlabel(1(2)20, grid) yscale(range($min $max)) yline(.5) ylabel($min($step)$max, grid) ///    
    legend(region(lpattern(blank))) xtitle("Education (Years)") ytitle("Predicted Probabilities") 	









// Bivariate probit model
use http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.dta, clear

quietly probit www educate income age male 

biprobit (trust = educate income age male) (www = educate income age male) 

// to get AIC and BIC
estat ic 


// conditional marginal effects
mfx, predict(pcond1) at(mean educate=16 male=0 www=1)

// marginal effects
mfx, predict(pmarg1) at(mean educate=16 male=0 www=1)



// recursive bivariate probit model
biprobit (trust = educate income age male www) (www = educate income age male) 

// to get AIC and BIC
estat ic 

// conditional marginal effects
mfx, predict(pcond1) at(mean educate=16 male=0 www=1)

// marginal effects
mfx, predict(pmarg1) at(mean educate=16 male=0 www=1)

mfx, predict(pmarg2) at(mean educate=16 male=0 www=1)





// compute direct and indirect effects manually
biprobit (trust = educate income age male www) (www = educate income age male) 

global rho=e(rho) 	// correlation coefficient of disturbances
global n1 = 6 		// the number of parameters in equation 1
global n2 = 5 		// the number of parameters in equation 2

// get the reference points 
//tabstat $rhs, stat(mean) col(variable) save
tabstat $rhs, stat(mean) col(stat) save

matrix ref1 = r(StatTotal),I(1) // reference points for equation 1
//matrix ref1[1,1]=16 			// education (college graduation)
//matrix ref1[1,4]=0 			// female
//matrix ref1[1,5]=1 			// WWW use
matrix ref0 = ref1[1,1..$n2]	// for display

matrix ref2 = ref1[1,1..$n2] // reference points for equation 2
matrix ref2[1,$n2]=1

// get parameter estimates 
matrix b0=e(b)

matrix b1=b0[1,1..$n1] // parameter estimates for equation 1
matrix b2=b0[1,$n1+1..$n1+$n2] // parameter estimates for equation 2

matrix xb1=b1*ref1' 	// compute xb1 of equation 1
matrix xb2=b2*ref2' 	// compute xb2 of equation 2

global xb1=xb1[1,1]		// put xb1 into a global macro for computation
global xb2=xb2[1,1]		// put xb1 into a global macro for computation

// compute the predicted probability at the reference points
di binormal($xb1, $xb2, $rho)/normal($xb2)

// compute direct effects 
global g1=normalden($xb1)*normal(($xb2-($rho)*$xb1)/sqrt(1-($rho)^2))
matrix directE=$g1/normal($xb2)*b1
matrix directE=directE[1,1..$n2]

// compute indirect effects
global g2=normalden($xb2)*normal(($xb1-($rho)*$xb2)/sqrt(1-($rho)^2))
matrix indirectE=($g2/normal($xb2)-(binormal($xb1,$xb2,$rho)*normalden($xb2))/(normal($xb2)^2))*b2
matrix indirectE[1,$n2]=0

// compute overall effects
matrix Overall=directE+indirectE

//mfx, predict(pcond1) at(mean www=1)
//mfx, predict(pcond1) at(mean www=0)

// to compute discrete changes
// list the location of binary variables
foreach bin of numlist 4 5 {
matrix ref10=ref1
matrix ref10[1,`bin']=0

matrix ref11=ref1
matrix ref11[1,`bin']=1

matrix ref20=ref2
if `bin' != 5 { 			// if the endogenous variable, do not change the value
matrix ref20[1,`bin']=0 
}

matrix ref21=ref2
if `bin' != 5 {	 			// if the endogenous variable, do not change the value
matrix ref21[1,`bin']=1 
}

matrix xb10=b1*ref10' 	// compute xb1 of equation 1 when d=0
matrix xb20=b2*ref20' 	// compute xb2 of equation 2 when d=0
matrix xb11=b1*ref11' 	// compute xb1 of equation 1 when d=1
matrix xb21=b2*ref21' 	// compute xb2 of equation 2 when d=1
global xb10=xb10[1,1]	// put xb10 to a global macro xb10
global xb20=xb20[1,1]
global xb11=xb11[1,1]
global xb21=xb21[1,1]

//di binormal($xb11, $xb21, $rho)/normal($xb21) // predicted probability when d=1
//di binormal($xb10, $xb20, $rho)/normal($xb20) // predicted probability when d=0

// discrete changes
matrix dchange = binormal($xb11, $xb21, $rho)/normal($xb21)- binormal($xb10, $xb20, $rho)/normal($xb20)

//di binormal($xb11, $xb21, $rho)/normal($xb21)- binormal($xb10, -($xb21), -($rho))/(1-normal($xb21)) 

matrix Overall[1,`bin']=dchange

} // foreach

// join matrices vertically
matrix Marginal=ref0\directE\indirectE\Overall 

// add row and column names
matrix coleq Marginal="" // remove equation names
matrix colnames Marginal = Education Income Age Male WWW
matrix rownames Marginal = Reference Direct Indirect Overall

// display marginal effects and discrete changes (direct and indirect effects)
matrix list Marginal 





// to compare with LIMDEP (Section 4.4)
mfx, predict(pcond1) at(mean male=.450596 www=.785349)

// end of Stata script