/*   This program caculates the Robust standard errors as described by Wooldridge (JOE 1999) in the case of conditional fixed effects poisson models */
/*   Reporting ML coefficient estimates with the robust SE's is equivalent to reporting quasi-ML as described in Wooldridge (1999) */
/*   Robust_revise_new.do written in 03/2007 by Wes Yin; instructions as follows:
      1. Run the "xtpoisson depvar varlist fe if sample == X " regression 
	2. Immediately after, run robust.do program as "do robust groupi timet depvar varlist"
	3. Note that program can take a maximum of 17 RHS variables; this can be changed in the opening args & MACRO list
	4. groupi is the name of the xsection variable in the panel; timet is the name of the panel variable
	5. Note that the robust program keeps only the observations used in the regression (ie. e(sample) == 1)
      6. Note that this program applies only to Robust SE calculation for Conditional FE Poisson models w/ balanced panels
	7. Unfortunately, in order to rerun the program for a new regression, the program clears stata of data/matrices*/


args groupi timet depvar var1 var2 var3 var4 var5 var6 var7 var8 var9 var10 var11 var12 var13 var14 var15 var16 var17	/* Program for variables with three regressors (help to imbed?) */

cd U:\lotka\orphan\revise\analysis

keep if e(sample)									/* keeps only those */
save robust_temp, replace							/* to avoid saving over panel_skel.dta */

global depvar `depvar'								/* allows for no_clins, no_clins_subd */
global groupi `groupi'
global timet `timet'
global var1 `var1'							
global var2 `var2'							
global var3 `var3'							
global var4 `var4'							
global var5 `var5'							
global var6 `var6'							
global var7 `var7'							
global var8 `var8'							
global var9 `var9'							
global var10 `var10'							
global var11 `var11'							
global var12 `var12'							
global var13 `var13'							
global var14 `var14'							
global var15 `var15'							
global var16 `var16'							
global var17 `var17'							
global var18 `var18'							
global var19 `var19'							
global var20 `var20'							

/******* Creating Variables to be Inputted into Matrices ********/
sort $groupi $timet

matrix VC=e(V)									/* defines initial V-C matrix (includes poisson alpha coefficient) */
matrix b = e(b)	
bysort $groupi: egen n_i = sum($depvar) if e(sample)		/* total number of clinical trials for disease i over t */
predict xb_hat_it if e(sample), xb
gen mu_hat_it = exp( xb_hat_it) if e(sample)				/* predicted mu for each disease, every t; also numerator for p_it */ 
bysort $groupi: egen mu_hat_i = sum( mu_hat_it)  if e(sample)   /* i-specific denominator for p_it */
gen p_it =  mu_hat_it/ mu_hat_i
gen u_it = $depvar - p_it*n_i							/* this is u_hat_it, which forms Tx1 vector u_hat_i; needs to be made a matrix */

/* Variables for Del_p_it */							/* in Stata's editor, this will form the TxP matrix; will need ot be transposed */

global i = 1
global colsb = colsof(b)

while $i <= $colsb {
	gen mu_hat_it_var$i = ${var$i}*mu_hat_it  if e(sample)				
	bysort $groupi:  egen mu_hat_i_var$i = sum(mu_hat_it_var$i) if e(sample)
	gen del_p_it_var$i = (${var$i}*mu_hat_it*mu_hat_i - mu_hat_it*mu_hat_i_var$i) / ((mu_hat_i)^2)
	global i = $i + 1
}

sort $groupi $timet
bysort $groupi: gen id = 1 if _n==1
sort id $groupi
egen id2 = fill(1 2)
replace id2 = . if id == .
sort $groupi $timet
replace id = id2
drop id2
bysort $groupi: egen id2 = sum(id)
replace id = id2
drop id2											/* id is new disease specific count variable */

*ereturn list
global n = e(N_g)										/* $n counts the total number of diseases included in the e(sample)*/
sort $groupi $timet


/***** Inputting into Matrices ******/

global i =1
global p = 1
global t = 1
global colsb = colsof(b)
global T = e(g_max)
matrix D_0 = J($T, $colsb, 0)	
matrix D_1 = J($T, $colsb, 0)	


gen id2 = id
sort id $timet

while $i <= $n {
	while $p <= $colsb {
			while $t <= $T {							/* fill in for all t's first */
                  matrix D_$i[$t,$p] = del_p_it_var$p[$t]			/* value for disease i, year t (1981), variable p */
			qui global t = $t + 1
		}
		qui global p = $p + 1
		qui global t = 1
		sort id2 $timet
	}
	qui replace id2 = id2 + e(N_g) if id == $i				/* this sequence moves the ith disease that was just done to the bottom */
	global i = $i + 1
	global p = 1
	matrix D_$i = J($T, $colsb, 0)		 				/* create the P x T (time periods x col's of b) in Del P matrix */
	sort id2 $timet
	
}

/***** P, u and W matrix *****/
sort id $timet
gen id3 = id
global i = 1
global t = 1
matrix P_0 = J($T, 1, 0)	
matrix P_1 = J($T, 1, 0)	
matrix U_0 = J($T, 1, 0)	
matrix U_1 = J($T, 1, 0)	
matrix B_0 = J($colsb, $colsb, 0)
matrix B = J($colsb, $colsb, 0)
matrix A_0 = J($colsb, $colsb, 0)
matrix A = J($colsb, $colsb, 0)

while $i <= $n {
	while $t <= $T {								/* fill in for all t's first */
                	matrix P_$i[$t,1] = p_it[$t]				/* value for disease i, year t (1981) of P matrix */
                	matrix U_$i[$t,1] = u_it[$t]				/* value for disease i, year t (1981) of U matrix */
			qui global t = $t + 1
	}
	qui global t = 1
	qui replace id3 = id3 + e(N_g) if id == $i			/* this sequence moves the ith disease that was just done to the bottom */
	qui global i = $i + 1
	matrix P_$i = J($T, 1, 0)				 		/* create the T x 1 (time periods x 1) in P matrix */
	matrix U_$i = J($T, 1, 0)				 		/* create the T x 1 (time periods x 1) in U matrix */
	sort id3 $timet	
}

/***** Get B matrix *****/
global i = 1

while $i <= $n {
	matrix S_$i = D_$i'*inv(diag(P_$i))*U_$i
	matrix B_$i = S_$i*S_$i'
	matrix B = B + B_$i
	global i = $i + 1
}

/***** Get A matrix; should be MLE VC matrix *****/
global i = 1
sort id $timet
gen id4 = id

while $i <= $n {
	global ni = n_i	
	matrix A_$i = $ni*D_$i'*inv(diag(P_$i))*D_$i
	matrix A = A + A_$i
	qui replace id4 = id4 + e(N_g) if id == $i			/* this sequence moves the ith disease that was just done to the bottom */
	global i = $i + 1
	sort id4 $timet
}


/***** Compare SE's ******/
matrix Vmle = vecdiag(VC)
global p = 1
while $p <= $colsb {
	global v = Vmle[1,$p]
	matrix SEmle = [nullmat(SEmle), sqrt($v)]
	global p = $p + 1
}

matrix Ai = inv(A)
matrix R = Ai*B*Ai
matrix Rmle = vecdiag(R)
global p = 1
while $p <= $colsb {
	global v = Rmle[1,$p]
	matrix RSEmle = [nullmat(RSEmle), sqrt($v)]
	global p = $p + 1
}


/************IGNORE****** IGNORE ********IGNORE**********IGNORE**********IGNORE**************IGNORE*********IGNORE************/
/*** For random effects model ***/
/*local m = rowsof(VC)								/* `m' is # of rows (and columns) of VC */
local m1 = `m'-1
matrix V = VC[1..`m1', 1..`m1']					/* creates VC matrix "V" with the avar for only the regression parameters */
matrix list V*/

matrix list SEmle									/* MLE SE's */
matrix list RSEmle								/* Robust SE's */
clear