/*
	To obtain confidence intervals for probit

		z = alpha + beta1*x1 + beta2*x2

	we use 	

		_b[x] +/- invt(13,.95)*_se[x]		(assuming _N=15)

	when we should use

		_b[x] +/- invnorm(.975)*_se[x]

	which is better?
*/

mac def ALPHA=.2
mac def BETA1=.4
mac def BETA2=1
mac def MODEL "probit y x1 x2"

capture program drop mkmodel /* n */
program define mkmodel
	local n `1'
	drop _all
	set obs `n'
	gen x1 = _n/`n' 
	gen x2 = (uniform() -.5)*2
	gen z = $ALPHA + $BETA1*x1 + $BETA2*x2
	gen p = normprob(z)
	gen y = p>=uniform()
end

capture program drop gens
program define gens
	gen alpha=.
	gen salpha=.
	gen beta1=.
	gen sbeta1=.
	gen beta2=.
	gen sbeta2=.
end

capture program drop sto
program define sto
	replace alpha = _b[_cons] in `1'
	replace salpha = _se[_cons] in `1'
	replace beta1 = _b[x1] in `1'
	replace sbeta1 = _se[x1] in `1'
	replace beta2 = _b[x2] in `1'
	replace sbeta2 = _se[x2] in `1'
end

capture program drop reports
program define reports
	report 95 `1'-3 $ALPHA alpha salpha 
	report 95 `1'-3 $BETA1 beta1 sbeta1
	report 95 `1'-3 $BETA2 beta2 sbeta2
end