program define divprob, eclass
	version 6.0
	syntax varlist  [if] [in], ENDOG(varlist) IV(varlist) EXOG(varlist) [CLAssic]
	marksample touse
	markout `touse' `cluster', strok
	tokenize `varlist'
	local lhs `"`1'"'
				/* 
					touse:  sample
					tosumm:  lhs mapped to 0/1
				*/
	tempvar tosumm
	tempname b x
	
	qui ivprob `varlist' if `touse',endog(`endog') iv(`iv') exog(`exog')
	if e(N)==0 | e(N)==. { exit 2000 }

* e(ll), e(N), e(df_m), e(chi2), e(r2_p), etc. are obtained from -probit- call

	mat `b' = e(b)
	local varlist : colnames(`b')
	local i : word count `varlist'
	tokenize `varlist'
	if `"``i''"' != "_cons" { 
		di in red `"may not drop constant"'
		exit 399
	}
	local `i'
	local varlist `"`*'"'
	
				/* check if drops from sample	*/
	quietly replace `touse'=0 if !e(sample)

			/* obtain Xb evaluation at means, f(Xb) */
	quietly gen byte `tosumm'=`lhs'!=0 if `touse'
	quietly summ `tosumm' [`sweight'`exp'] if `touse'
	est scalar pbar = r(mean)
	matrix `x' = e(mns) * `b''
	est scalar xbar = `x'[1,1]
	if "`e(offset)'" != "" { 
		qui summ `e(offset)' if `touse'
		est scalar offbar = r(mean)
		est scalar xbar = e(xbar) + r(mean)
	}
	else	est scalar offbar = 0
				/* mark the dummies 			*/
	local i 1 
	while `"``i''"'!="" {
		capture assert ``i''==0 | ``i''==1 if `touse'
		if _rc==0 { est local dummy `"`e(dummy)' 1"' }
		else 	  { est local dummy `"`e(dummy)' 0"' }
		local i = `i' + 1
	}
	est local dummy `"`e(dummy)' 0"'
			
			/* save other info			*/
	est local wtype `"`weight'"'
	est local wexp `"`exp'"'

	/* double save in S_E_<stuff> and e()  */
	scalar S_E_pbar = e(pbar)
	scalar S_E_xbar = e(xbar)
	global S_E_dum `e(dummy)'
	global S_E_vl `"`lhs' `varlist'"'
	global S_E_if `"`if'"'
	global S_E_in `"`in'"'
	global S_E_wgt `e(wtype)'
	global S_E_exp `e(wexp)'
	global S_E_ll `e(ll)'
	global S_E_nobs `e(N)'
	global S_E_mdf `e(df_m)'
	est local cmd "dprobit"
	global S_E_cmd "dprobit"
	/* see farther down for saving of e(at) and S_E_at  */
	/* also see e(dfdx) and e(se_dfdx) */

	tempname mns b V Vhat v xb fxb c s ll ul z Z 
	local level 95

	scalar `Z' = invnorm(1-(1-`level'/100)/2)

	mat `b' = e(b)
	mat `V' = e(V)

	mat `mns' = e(mns)
	local ttl `"x-bar"'


*	probit, nocoef
	
	di _n in gr `"Probit estimates"' _col(57) `"Number of obs ="' /*
		*/ in ye %7.0f e(N)
	di in gr _col(57) `"`e(chi2type)' chi2("' in ye e(df_m) /*
		*/ in gr `")"' _col(71) /*
		*/ `"="' in ye %7.2f e(chi2)
	di in gr _col(57) `"Prob > chi2   ="' in ye /*
		*/ %7.4f chiprob(e(df_m),e(chi2))
	di in gr `"Log likelihood = "' in ye %10.0g e(ll) /*
		*/ _col(57) in gr `"Pseudo R2"' _col(71) `"="' /*
		*/ in ye %7.4f e(r2_p) _n


	di in gr _dup(78) `"-"'

        if `"`e(vcetype)'"'!="" {
		di in gr _col(10) `"|"' _col(26) `"`e(vcetype)'"'
	}

	if `"`classic'"'!="" {
		local cons `"_cons"'
	}
	_evlist
	tokenize `e(depvar)' `s(varlist)' `cons'
	sret clear
	local skip = 8-length(`"`1'"')

	di in gr /*
		*/ _skip(`skip') `"`1' |"' _col(17) `"dF/dx"' _col(25) /*
		*/ `"Std. Err."' _col(40) `"z"' _col(45) `"P>|z|"' /*
		*/ _col(55) `"`ttl'"' /*
		*/ _col(62) `"[    `level'% C.I.   ]"' _n /*
		*/ _dup(9) `"-"' `"+"' _dup(68) `"-"'

	Vhatdiag `mns' `b' `V' -> `Vhat'

	mat `xb' = `mns'*`b''
	scalar `xb' = `xb'[1,1] + e(offbar)
	scalar `fxb' = normd(`xb')

	local lhs `"`1'"'
	mac shift

	tempname dfdx se_dfdx
	local i 1
	while `"``i''"'!="" {
		local isdum : word `i' of `e(dummy)'
		if `isdum' & `"`classic'"'=="" {
			local anydum "true"
			local star `"*"'
			scalar `c' = /*
			*/ normprob(`xb'+(1-`mns'[1,`i'])*_b[``i'']) - /*
			*/ normprob(`xb'-(`mns'[1,`i'])*_b[``i'']) 
			vdummy `i' `mns' `b' `V' -> `v'
		}
		else {
			scalar `c' = _b[``i'']*`fxb'
			scalar `v' = `Vhat'[1,`i']
			local star `" "'
		}
		scalar `s' = sqrt(`v')
		scalar `ll' = `c'-`Z'*`s'
		scalar `ul' = `c'+`Z'*`s'
		scalar `z' = _b[``i'']/_se[``i'']

		local skip=8-length(`"``i''"')
		di in gr _skip(`skip') `"``i''`star'|  "' in ye /*
			*/ %9.0g `c' `"  "' /*
			*/ %9.0g `s' `" "' /*
			*/ %8.2f `z'  `"  "' /*
			*/ %6.3f 2*normprob(-abs(`z')) `"  "' /* 
			*/ %8.0g `mns'[1,`i'] `"  "' /*
			*/ %8.0g `ll' `" "' /*
			*/ %8.0g `ul'
		mat `dfdx' = nullmat(`dfdx') , `c'
		mat `se_dfdx' = nullmat(`se_dfdx') , `s'

*		`debug' `c' `s' `z' `mns'[1,`i'] `ll' `ul'
		local i=`i'+1
	}
	if "`e(offset)'"!="" { 
		di in gr %8s "`e(offset)'" " |   (offset)"
	}

	di in gr _dup(9) `"-"' `"+"' _dup(68) `"-"' _n /*
		*/ `"  obs. P |  "' in ye %9.0g e(pbar)
	di in gr `" pred. P |  "' in ye %9.0g normprob(e(xbar)) /*
		*/ in gr `"  (at x-bar)"'
	if `"`at'"'!="" {
		est scalar at = normprob(`xb')
		scalar S_E_at = e(at)
		di in gr `" pred. P |  "' in ye %9.0g /* 
		*/ e(at) in gr `"  (at x)"'
	}
	else	capture scalar drop S_E_at

	di in gr _dup(78) `"-"'
	if `"`anydum'"'==`"true"' {
		di in gr /*
		*/ `"(*) dF/dx"' /*
		*/ `" is for discrete change of dummy variable from 0 to 1"'
	}
	di in gr /*
     */ `"    z and P>|z| are the test of the underlying coefficient being 0"'

     tempname B
     matrix `B' = e(b)
     local names : colfullnames `B'
     local names : subinstr local names "_cons" ""
     mat colnames `dfdx' = `names'
     mat colnames `se_dfdx' = `names'

     est mat dfdx `dfdx'
     est mat se_dfdx `se_dfdx'
end


program define Vhatdiag /* mns b V -> res */
	args xbar b V ARROW res

	tempname tmp dgdb

	local k = colsof(`V')

	mat `tmp' = `xbar'*`b'' + e(offbar)
	mat `dgdb' = normd(`tmp'[1,1]) * ( I(`k') - `tmp'[1,1]*`b''*`xbar' )
	mat `res' = vecdiag(`dgdb'*`V'*`dgdb'')

end


program define vdummy /* i mns b V -> res */
	args i mns b V ARROW res

	tempname x xb fxb Delp Vhat
	mat `x' = `mns'
	mat `x'[1,`i']=1
	mat `xb' = `x'*`b'' + e(offbar)
	scalar `fxb' = normd(`xb'[1,1])
	mat `Delp' = `fxb' * `x'

	mat `x'[1,`i']=0
	mat `xb' = `x'*`b'' + e(offbar)
	scalar `fxb' = normd(`xb'[1,1])
	mat `xb' = `fxb'*`x' 
	mat `Delp' = `Delp' - `xb'

	mat `Vhat' = `Delp'*`V'*`Delp''
	scalar `res' = `Vhat'[1,1]
end

		     /*  1    2    3   4      5          6    7 */
program define Debug /* name  `c' `s' `z' `mns'[1,`i'] `ll' `ul' */
	local name `"`1'"'
	local c = `2'
	local s = `3'
	local z = `4'
	local m = `5' 
	local l = `6'
	local u = `7'
	tempname row
	mat `row' = (`c',`s',`z',`m',`l',`u')
	mat `name' = nullmat(`name') \ `row'
end
exit

Dummies

        Dp = F(x2*b) - F(x1*b)

        d(Dp)/db = f(x2*b)*x2 - f(x1*b)*x1

        Vhat(Dp) = [d(Dp)/db)]' V [d(Dp)/db]


Continuous 

        p = F(x*b)

        dp/dx = f(x*b)*x = g(b)

        dg/db = f'(x*b)b'x + f(x*b)I

        Vhat  = [dg/db] V [dg/db]'