capture log close
log using sn.log , replace 

// do file to replicate SN results
// version 3.0.1 
// last edited David M. Drukker  30Sep2003


set matsize 800
version 8.1 

drop _all

use sndata

rename A a
replace VOTPOP = 1000*VOTPOP
sum 

capture program drop sn_ll

global doov DOOV
global tiebreak TIEBREAK
global a a

program define sn_ll
	version 8.1 
	args lnf dbar1 betar0 betad0 sep rho sed winst votpop tm xb1

	tempvar gamma alpha delta dbar Ef xb2 prob 

	gen double `gamma' = exp(`betar0' - `betad0') 

	gen double `alpha' = `gamma'/(`gamma' + 1)
	
	gen double `delta' = `gamma'/( (1/$doov) - (1-`gamma') )

	gen double `dbar'  = `dbar1'*(1-$a) + `delta'*$a

	gen double `Ef'    = ($tiebreak + `winst')*`rho'*`alpha'*	///
		(1-`alpha')*normden((`alpha'-`dbar')/`sed')/`sed' -	///
		`votpop'

	gen double `xb2'   = `xb1' + `betad0' + `tm'	
	gen double `prob'  =  (normden( (ln($ML_y1/`delta') - 		///
		`rho'*`Ef' - `xb2')/`sep')/`sep')/$ML_y1*(		///
		(normden( ( `delta'-`dbar')/`sed')/`sed')*		///
		`gamma'/((1-(1-`gamma')*$doov)^2)*(1-$a) +$a)

	qui replace `lnf' = ln(`prob')
end

ml model lf sn_ll (DOOP = TRIHEAT JULYEC2 PRESINC VP4 HOME4X ADAACA	///
	DEV3 DEV2 LEGIS ZECON41 DSOUTH64 SOUTH4 SOUTH64B WEST76 	///
	NCENT72 NEWENG64 ) /betar0 /betad0 /sep 	 		///
	/rho /sed /winst (votpop: VOTPOP, nocons) 			///
	(tm: DUMMY88 DUMMY84 DUMMY80 DUMMY76 DUMMY72 DUMMY68 DUMMY64	///
	DUMMY60 DUMMY56 DUMMY52,nocons)					///
	(xb1: GOVRACE RAIN PCAPREAL EDUCAT JIMCROW PCBLCK PCIN 		///
	OUTSTATE, nocons)

local sed 0.037847927  
local d0 0.29339712 
local cm1 0.0051469912 
local cm2 0.020902094 
local cm3 0.013200831 
local cm4  0.019494456  
local cm5 0.060477454   
local cm6 0.00035823595 
local cm7 0.0033246381  
local cm8   0.0025325467 
local cm9 0.00036951962  
local cm10 0.0063127442
local cm11 -0.14457178  
local cm12 0.07764634  
local cm13 -0.092488749 
local cm14 -0.071244138
local cm15 0.05687312   
local cm16 0.067957626 
local sep 0.11561017  
local betad0 -1.1301832
local betar0 -0.8589113  
local bgov 0.045194958 
local bpre -0.079465481  
local bjim -0.35111151
local bvinc 0.21899693  
local bvblk -1.0003932 
local bin -0.032955224 
local bedu 0.3623056
local tm1 0.04665423 
local tm2 0.12062019  
local tm3 0.032472619 
local tm4 0.1464833 
local tm5 0.16505689 
local tm6 0.017293131  
local tm7 0.32905255   
local tm8 0.26481688
local tm9 0.20499258 
local tm10 0.26819455 
local outst -0.17714668 
local rho 0.58661552
local zoe0 0.000014910489 
local winst 0.078760921
	
/* 
cm1
cm2 
cm3
cm4 
cm5 
cm6
cm7
cm8
cm9
cm10
cm11
cm12
cm13
cm14
cm15
cm16
d0
betar0
betad0
sep
rho
sed
winst
zoe0
tm1 
tm2 
tm3
tm4
tm5
tm6
tm7
tm8
tm9
tm10
bgov
bpre
bvinc
bedu
bjimmat 
bvblk
bin
outst
*/

mat b0 = [ `cm1', `cm2 ', `cm3', `cm4 ', `cm5 ', `cm6', `cm7', `cm8', 	///
	`cm9', `cm10', `cm11', `cm12', `cm13', `cm14', `cm15', `cm16', 	///
	`d0', `betar0', `betad0', `sep', `rho', `sed', `winst', `zoe0', ///
	`tm1', `tm2', `tm3', `tm4', `tm5', `tm6', `tm7', `tm8', `tm9', 	///
	`tm10', `bgov', `bpre', `bvinc', `bedu', `bjim', `bvblk', 	///
	`bin', `outst'] 
 
ml init b0, copy
ml max, difficult grad  nrtol(1e-8) gtol(1e-8) 
eret list
mat vhat = e(V)
mat symeigen X v = vhat
di "the condition number (in 2-norm) is " v[1,1]/v[1,42]
di "The eigenvalues of e(V) are "
mat list v


/* now do with with rescaled variables */

 
replace VOTPOP = VOTPOP/1000
replace JULYEC2 = JULYEC2/100
replace TRIHEAT = TRIHEAT/100
replace PRESINC = PRESINC/10
replace VP4 = VP4/10
replace ADAACA = ADAACA/100
replace LEGIS = LEGIS/100
replace DEV2 = DEV2/100
replace DEV3 = DEV3/100

mat b0[1,24] = b0[1,24]*1000
mat b0[1,1] = b0[1,1]*100
mat b0[1,2] = b0[1,2]*100
mat b0[1,3] = b0[1,3]*10
mat b0[1,4] = b0[1,4]*10
mat b0[1,6] = b0[1,6]*100
mat b0[1,9] = b0[1,9]*100

mat b0[1,7] = b0[1,7]*100
mat b0[1,8] = b0[1,8]*100


ml model lf sn_ll (DOOP = TRIHEAT JULYEC2 PRESINC VP4 HOME4X ADAACA	///
	DEV3 DEV2 LEGIS ZECON41 DSOUTH64 SOUTH4 SOUTH64B WEST76 	///
	NCENT72 NEWENG64 ) /betar0 /betad0 /sep 	 		///
	/rho /sed /winst (votpop: VOTPOP, nocons) 			///
	(tm: DUMMY88 DUMMY84 DUMMY80 DUMMY76 DUMMY72 DUMMY68 DUMMY64	///
	DUMMY60 DUMMY56 DUMMY52,nocons)					///
	(xb1: GOVRACE RAIN PCAPREAL EDUCAT JIMCROW PCBLCK PCIN 		///
	OUTSTATE, nocons)

 
ml init b0, copy
ml max, grad  nrtol(1e-8) gtol(1e-8)
eret list
mat vhat = e(V)
mat symeigen X v = vhat
di "the condition number (in 2-norm) is " v[1,1]/v[1,42]
di "The eigenvalues of e(V) are "
mat list v

mat b1 = e(b)

di b1[1,24]/1000
di b1[1,1]/100
di b1[1,2]/100
di b1[1,3]/10
di b1[1,4]/10
di b1[1,6]/100
di b1[1,9]/100
di b1[1,7]/100
di b1[1,8]/100


capture log close