******************************************
* This is the Stata do-file for the paper:
* Andrew Leigh, "Estimating Teacher Effectiveness From Two-Year Changes in Students’ Test Scores", Economics of Education Review, 2010, vol 29, pp.480-488
* Feel free to use or modify the do-file, so long as you cite that paper.
* The data in the paper are confidential. Requests to access them should be directed to the Queensland Department of Education
******************************************

version 10
clear
set mem 60m
set matsize 3000
program drop _all
set more off
cd <DATA DIRECTORY>

/*

cd <DATA DIRECTORY>

* Here's a run-down of the key variables:
* A file: record_id, education_dept_id, centre_code, date_of_birth, test score, year
* B file: education_dept_id, roll_class, centre_code, date_of_birth, student_id
* C file: centre_code, roll_class, payroll_number
* D file: Payroll_Number, 4_years, Masters, Date_of_Birth (of teacher), Gender, Start_Date, Suitability_Rating, Status

* Matching
* 2001-3 - 2003-5
* 2002-3 - 2004-5
* 2001-5 - 2003-7
* 2002-5 - 2004-7

* The B files (2003) is missing variable names
*use "b qry individual student 2003", clear
*for num 1/10 \ any year collection_date centre_code student_id date_of_birth year_level education_dept_id date_enrolled roll_class time_series: ren colX Y
*save, replace

* Get staff characteristics ready to merge on
use "d qry staff characteristics.dta", clear
* We recode one miscoded starting date
replace date_of_birth=d(04Nov1924) if date_of_birth==d(04Nov2024)
gen masters=0 if masters_=="No"
replace masters=1 if masters_=="Yes"
gen fouryr=0 if _4_year=="No"
replace fouryr=1 if _4_year=="Yes"
gen birthyear=year(date_of_birth)
gen startyear=year(start_date)
* Rating=T4 means not assessed, so I code it as missing.
gen rating=.
for num 1/4: replace rating=X if suitability=="SX"
drop masters_ _4_year 
gen female=1 if gender=="F"
replace female=0 if gender=="M"
ren payroll teacher
sort teacher
save temp_d, replace

cd <DATA DIRECTORY>
* Match teacher payroll numbers to a classroom identifier (roll_class & centre_code)
*for num 3 3 5 5 5 5 7 7 \ num 2001/2002 2001/2004 2003/2004: use "a qry year X results Y.dta", clear \ gen grade=X \ sort education_dept_id \ save, replace
*for num 3 3 5 5 5 5 7 7 \ num 2001/2002 2001/2004 2003/2004: use "a qry year X results Y.dta", clear \ destring centre_code, force replace \ ren date_of_birth dob \ sort dob centre \ save temp_a_X_Y, replace
*for num 2001/2004: use "c qry student roll class X", clear \ destring centre_code, force replace \ sort centre_code roll_class \ save temp_c_X, replace

#delimit ;
for num 3 3 5 5 5 5 7 7 \ num 2001/2002 2001/2004 2003/2004: use "a qry year X results Y.dta", clear \ 
destring centre_code, force replace \ ren date_of_birth dob \ destring education_dept_id, force replace \ 
\ gen temp1=1 \ bysort education_dept_id dob: egen temp2=count(temp1) \ sum literacy_o numeracy_o if temp2>1 \
drop if literacy_o==0 & numeracy_o==0 & temp2>1 \ drop temp1 temp2 \
\ gen temp1=1 \ bysort education_dept_id dob: egen temp2=count(temp1) \
drop if literacy_o==0 & temp2>1 \ drop temp1 temp2 \
\ gen temp1=1 \ bysort education_dept_id dob: egen temp2=count(temp1) \
drop if numeracy_o==0 & temp2>1 \ drop temp1 temp2 \
sort education_dept_id dob centre_code \ save temp_a_X_Y, replace;

for num 2001/2004: use "b qry individual student X", clear \ bysort centre_code roll_class: egen class_sizeX=count(year) \ 
sort centre_code roll_class \ merge centre_code roll_class using temp_c_X, keep(payroll_number) nokeep \ 
tab _merge \ drop _merge collection_date time_series \ gen dob=date(date_of_birth,"DMY") \ 
replace educ=substr(educ,1,8) \ destring education_dept_id, force replace \ 
sort education_dept_id dob centre_code \ save temp_b_X, replace;

#delimit cr

* Match students' two test scores
program drop _all
program define matchscores
use temp_a_$gy1.dta, clear
for any literacy numeracy: recode X_overall 0=. \ egen X_z$y1=std(X_overall) \ drop X_overall
keep grade education_dept_id dob centre *_z$y1
ren centre_code cc_check$y1
sort education_dept_id dob cc_check
merge education_dept_id dob using temp_a_$gy2, nokeep
keep if _merge==3
ren centre_code cc_check$y3
gen cc_check$y2=cc_check$y1
for any literacy numeracy: recode X_overall 0=. \ egen X_z$y3=std(X_overall) \ drop X_overall
keep grade education_dept_id dob cc* *_z$y3 *_z$y1
*for any literacy numeracy: gen dX=X_z$y3-X_z$y1
for num $y1 $y2 $y3: sort education_dept_id dob cc_checkX \ merge education_dept_id dob using temp_b_X, keep(roll_class payroll_number centre_code class_size*) nokeep \ tab _merge \ drop _merge \ ren payroll_number teacherX \ ren centre_code centre_codeX \ replace teacherX=. if centre_codeX~=cc_check$y1 & centre_codeX~=cc_check$y3
*list education_dept_id centre_code* in 1/100
drop cc_check*
gen studentid=_n*10+$iter
expand 3
bysort studentid: egen weight=seq()
for var class_size*: recode X 40/max=40 .=.
gen class_size_mean=(class_size$y1+class_size$y2+class_size$y3)/3
gen class_size=.
for num 1/3 \ num $y1 $y2 $y3: replace class_size=class_sizeY if weight==X \ drop class_sizeY 
gen centre_code=.
for num 1/3 \ num $y1 $y2 $y3: replace centre_code=centre_codeY if weight==X \ drop centre_codeY 
gen teacher=.
for num 1/3 \ num $y1 $y2 $y3: replace teacher=teacherY if weight==X \ drop teacherY 
replace grade=weight+$g1-1
for num 1 2 \ num $y1 $y3: ren literacy_zY literacyX \ ren numeracy_zY numeracyX
gen year=weight
recode year 1=$y1 2=$y2 3=$y3
recode weight 1=.5 2=1 3=.5
sort studentid
save "temp_$iter.dta", replace
end

for num 1/3 \ num 2001/2003: global yX=Y
for num 1 2 \ num 3 5: global gX=Y
for X in num 1 2 \ Y in num 3 5 \ Z in num 2001 2003: global gyX="Y_Z"
global iter 1
matchscores

for num 1/3 \ num 2002/2004: global yX=Y
for num 1 2 \ num 3 5: global gX=Y
for X in num 1 2 \ Y in num 3 5 \ Z in num 2002 2004: global gyX="Y_Z"
global iter 2
matchscores

for num 1/3 \ num 2001/2003: global yX=Y
for num 1 2 \ num 5 7: global gX=Y
for X in num 1 2 \ Y in num 5 7 \ Z in num 2001 2003: global gyX="Y_Z"
global iter 3
matchscores

* Note: education_dept_id is missing for all records in the 2004 grade 7 file
for num 1/3 \ num 2002/2004: global yX=Y
for num 1 2 \ num 5 7: global gX=Y
for X in num 1 2 \ Y in num 5 7 \ Z in num 2002 2004: global gyX="Y_Z"
global iter 4
matchscores

*/

/*
***************************
* Teacher quality
***************************
clear
set mem 40m
set more off
cd <DATA DIRECTORY>
set matsize 1000
use temp_1, clear
for num 2/3: append using temp_X

* Generating grade*year fixed effects
egen gradeyr=group(grade year)

egen tag=tag(studentid year)
tab grade year if tag
drop tag


bysort studentid: egen seq=rank(year)

* Assigning relevant score to each obs
gen literacy=literacy2 if seq==3
gen numeracy=numeracy2 if seq==3
replace literacy=literacy1 if seq==1
replace numeracy=numeracy1 if seq==1
replace literacy=0.5*(literacy1+literacy2) if seq==2
replace numeracy=0.5*(numeracy1+numeracy2) if seq==2
drop literacy1 literacy2 numeracy1 numeracy2

* Dropping extraneous variables
keep numeracy literacy roll_class class_siz* grade year gradeyr centre_code teacher weight studentid seq
compress

save temp_teacher, replace
*/

/*
clear
*set maxvar 2048
set matsize 50
set mem 1300m
set more off
cd <DATA DIRECTORY>
use temp_teacher, clear

* Using a two-way fixed effects model
* Source: T. Cornelissen. "The Stata command felsdvreg to fit a linear model with two high-dimensional fixed effects". Stata Journal 8: 170–189.
*xi: felsdvreg literacy class_size i.gradeyr, ftest cons ivar(studentid) jvar(teacher) feff(fe_lit) normalize feffse(se_lit) peff(temp1) xb(temp2) res(temp3) mover(temp4) group(temp5) mnum(temp6) pobs(temp7)
*drop temp*

* Dropping out cases where a student only has one literacy or one numeracy observation
for any literacy numeracy: bysort studentid: egen count=count(X) \ tab count \ replace X=. if count<3 \ drop count

* De-meaning the data to take out student effects (follows http://www.stata.com/statalist/archive/2003-01/msg00688.html)
for any literacy numeracy: gen X_raw=X
for any literacy numeracy: bysort studentid:egen X_m=mean(X) \ egen X_m2=mean(X) \ gen temp=X-X_m+X_m2 \ replace X=temp \ drop temp X_m X_m2

* For Table 1
for any literacy numeracy: xi: xtreg X class_size i.gradeyr, i(teacher) fe \ tab year grade if e(sample)

* Table 2, Panel A (excluding non-test years)
* Without student FE
xi: xtreg literacy_raw class_size i.gradeyr if seq~=2, i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster replace bdec(4) rdec(4) se ct("Panel A - Lit - No SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
xi: xtreg numeracy_raw class_size i.gradeyr if seq~=2, i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel A - Num - No SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
* With student FE
xi: xtreg literacy class_size i.gradeyr if seq~=2, i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel A - Lit - SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
codebook centre_code studentid if e(sample)
sum class_size if e(sample)
xi: xtreg numeracy class_size i.gradeyr if seq~=2, i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel A - Num - SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
codebook centre_code studentid if e(sample)
sum class_size if e(sample)

* Table 2, Panel B (interpolating non-test years)
* Without student FE
xi: xtreg literacy_raw class_size i.gradeyr , i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel B - Lit - No SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
xi: xtreg numeracy_raw class_size i.gradeyr , i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel B - Num - No SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
* With student FE
xi: xtreg literacy class_size i.gradeyr , i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel B - Lit - SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
predict fe_lit, u
codebook centre_code studentid if e(sample)
sum class_size if e(sample)
xi: xtreg numeracy class_size i.gradeyr , i(teacher) fe
outreg using table_student.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se ct("Panel B - Num - SFE") addstat("F-test for joint significance of teacher fixed effect terms",e(F_f))
predict fe_num, u
codebook centre_code studentid if e(sample)
sum class_size if e(sample)

* Check - is it similar without class size or grade*year FE?
log using robustness,replace
xi: xtreg literacy i.gradeyr, i(teacher) fe
xi: xtreg literacy class_size , i(teacher) fe
xi: xtreg literacy , i(teacher) fe
xi: xtreg numeracy i.gradeyr, i(teacher) fe
xi: xtreg numeracy class_size , i(teacher) fe
xi: xtreg numeracy , i(teacher) fe
log close

* Estimating fixed effects with data aggregated to the grade-year-school level
egen gradeyrschool=group(gradeyr centre_code)
xi: xtreg literacy class_size , i(gradeyrschool) fe
predict fe_lit_gy, u
xi: xtreg numeracy class_size , i(gradeyrschool) fe
predict fe_num_gy, u
save teacher_fe_gy,replace
drop gradeyrschool

/*
* Saving the fixed effects and their standard errors (asterisked out for faster running)
* Uses Nichols, Austin. 2008.  fese: Stata module calculating standard errors for fixed effects. http://ideas.repec.org/c/boc/bocode/s456914.html
drop fe_*
xi: fese literacy class_size i.gradeyr, a(teacher) oonly
ren b fe_lit
ren se se_lit
xi: fese numeracy class_size i.gradeyr, a(teacher) oonly
ren b fe_num
ren se se_num
for var fe_lit se_lit fe_num se_num: la var X "Estimated with linear interpolation for non-tested year"
xi: fese literacy class_size i.gradeyr if seq~=2, a(teacher) oonly
ren b fe_lit2
ren se se_lit2
xi: fese numeracy class_size i.gradeyr if seq~=2, a(teacher) oonly
ren b fe_num2
ren se se_num2
for var fe_lit2 se_lit2 fe_num2 se_num2: la var X "Estimated by omitting non-tested year"
*/

*slash-asterisk
* Test for streaming within schools, by looking at whether a student's relative rank in grade in year t predicts relative rank of her class (relative to grade) in year t+1
tsset studentid year
for any literacy numeracy: bysort centre_code gradeyr: egen mean1X=mean(X_raw) \ bysort centre_code gradeyr: egen count1X=count(X_raw) \ replace mean1X=((mean1X*count1X)-X_raw)/(count1X-1) \ gen rank1X=X_raw-mean1X \ drop mean1X count1X \ la var rank1X "Student's within-grade rank in X in year t (z-score diff)"
for any literacy numeracy: bysort centre_code gradeyr: egen mean2X=mean(X_raw) \ bysort centre_code gradeyr: egen count2X=count(X_raw) \ replace mean2X=((mean2X*count2X)-X_raw)/(count2X-1) 
for any literacy numeracy: bysort centre_code gradeyr roll_class: egen mean3X=mean(X_raw) \ bysort centre_code gradeyr roll_class: egen count3X=count(X_raw) \ replace mean3X=((mean3X*count3X)-X_raw)/(count3X-1)  \ tsset studentid year \ gen rank2X=f2.mean3X-f2.mean2X \ la var rank2X "Student's avg class score in year t+2 relative to avg grade score in year t+2 (z-score diff)"
sum rank*,d
for any literacy numeracy: reg rank2X rank1X
* Interpretation: if a student is 0.1SD above her class mean in year t (ie. she is at the 75th percentile), 
* then she will be placed in a class in year t+2 that is 0.003 to 0.006 SD below the grade average.

* Test for non-random assignment of students across classes
* Does a student's relative rank in grade in year t predict her teacher's characteristics in year t+2?
sort teacher
merge teacher using temp_d, nokeep
gen test_year=year
gen age=test_year-birthyear
gen experience=test_year-startyear
replace masters=1 if fouryr==1
for num 1/2: gen ratingX=1 if rating==X \ recode ratingX .=0 if rating~=.
gen rating34=1 if rating==3 | rating==4 
recode rating34 .=0 if rating~=. 
tsset studentid year
for any literacy numeracy: gen rank_X_lag=l2.rank1X \ gen X_lag=l2.X_raw
egen gradeyrschool=group(gradeyr centre_code)
egen numteachers=nvals(teacher), by(gradeyrschool) 
sum numteachers,d
for any literacy numeracy: areg rank_X_lag age experience masters rating1 rating2 rating34 if numteachers>=2 & numteachers~=.,a(gradeyrschool) \ areg X_lag age experience masters rating1 rating2 rating34 if numteachers>=2 & numteachers~=.,a(gradeyrschool)
drop numteachers test_year gradeyrschool rank_literacy_lag rank_numeracy_lag literacy_lag numeracy_lag mean* count* rank*  date_of_birth gender start_date suitability_rating status masters fouryr birthyear startyear rating* female age experience temp* _merge
*asterisk-slash

bysort teacher: egen count_num=count(numeracy)
bysort teacher: egen count_lit=count(literacy)
bysort teacher: egen test_year=mean(year)
egen tag=tag(teacher)
keep if tag
keep centre_code teacher fe_* se_* count_* test_year
sort teacher
save teacher_fe, replace



/*
*************************************************************************************************************************************
* Explicitly estimating 'lagged teacher' effects (note: there are 5906 lagged teachers, so felsdvreg doesn't run)
*************************************************************************************************************************************
use teacher_fe, clear
tsset studentid year
gen lag_teacher=l.teacher if seq==3
drop if lag_teacher==.
tsset, clear
* This involves separately estimating two sets of fixed effects (current and lagged teachers). We therefore need to use felsdvreg 
*xi: xtreg literacy class_size i.gradeyr i.lag_teacher, i(teacher) fe
*xi: felsdvreg literacy class_size i.gradeyr, ftest cons ivar(teacher) jvar(lag_teacher) feff(fe_lit_lag) normalize feffse(se_lit_lag) peff(temp1) xb(temp2) res(temp3) mover(temp4) group(temp5) mnum(temp6) pobs(temp7)
*drop cons temp*
*xi: felsdvreg numeracy class_size i.gradeyr, ftest cons ivar(teacher) jvar(lag_teacher) feff(fe_num_lag) normalize feffse(se_num_lag) peff(temp1) xb(temp2) res(temp3) mover(temp4) group(temp5) mnum(temp6) pobs(temp7)
*keep studentid year fe_lit_lag se_lit_lag fe_num_lag se_num_lag 

* De-meaning the data to take out *teacher* effects 
for any literacy numeracy: bysort teacher:egen X_m=mean(X) \ egen X_m2=mean(X) \ gen temp=X-X_m+X_m2 \ replace X=temp \ drop temp X_m X_m2
xi: fese literacy class_size i.gradeyr, a(lag_teacher) oonly
ren b fe_lit_lag
ren se se_lit_lag
xi: fese numeracy class_size i.gradeyr, a(lag_teacher) oonly
ren b fe_num_lag
ren se se_num_lag
egen tag=tag(lag_teacher)
keep if tag
keep fe_lit_lag se_lit_lag fe_num_lag se_num_lag lag_teacher
ren lag_teacher teacher
sort teacher
save temp_lagged_teacher_effects,replace
use teacher_fe, clear
* Aggregating lagged teacher effects into regular teacher effects
sort teacher
merge teacher using temp_lagged_teacher_effects
drop _merge
for any lit num: gen fe_X3=(fe_X_lag+fe_X2)/2 \ gen se_X3=(se_X_lag+se_X2)/2 \ replace fe_X3=fe_X2 if fe_X3==. \ replace se_X3=se_X2 if se_X3==. \ replace fe_X3=fe_X_lag if fe_X3==. \ replace se_X3=se_X_lag if se_X3==.
for var fe_???3 se_???3: la var X "Estimated by explicitly modelling non-tested years"
sort teacher
save teacher_fe_lagged, replace
sum fe_*  
for any lit num: metareg fe_X3, wsse(se_X3) bsest(eb) 
*/

*/

*****************************************************************************************
* Looking at which teacher demographics correlate with higher value-added
*****************************************************************************************
set more off
cd <DATA DIRECTORY>
use teacher_fe, clear

egen tag=tag(teacher)
keep if tag
drop tag

merge teacher using temp_d, nokeep
tab _merge
drop _merge

* Getting a precise estimate of age and experience
* (This exploits the fact that we have precise birth dates and start dates)
gen temp1=doy(date_of_birth)
gen age=test_year-birthyear+(temp1/365)
gen temp2=doy(start_date)
gen experience=test_year-startyear+(temp2/365)
drop temp*
sum experience age, d

* Generating weights
*for any lit num: gen w_X=count_X
*gen w_both=w_lit+w_num
for any lit num: gen w_X=1/se_X
gen w_both=(w_lit+w_num)/2

* De-meaning literacy and numeracy scores
for var fe_*: sum X \ replace X=X-r(mean)

* Note that fouryr=1 means more than 4-year trained. So we code them as masters=1
replace masters=1 if fouryr==1

* Keeping only teachers with nonmissing obs
for any num lit: reg fe_X masters female age experience [aw=w_X], r cl(teacher) \ drop if e(sample)~=1

* EQ rating
for num 1/2: gen ratingX=1 if rating==X \ recode ratingX .=0 if rating~=.
gen rating34=1 if rating==3 | rating==4 
recode rating34 .=0 if rating~=. 

* Summary statistics (unweighted)
for var masters female age experience rating1 rating2 rating34: format X %9.3f
sum masters female age experience rating1 rating2 rating34, sep(0) format 

log using metareg, replace
sum fe_*  
for any lit num: metareg fe_X, wsse(se_X) bsest(eb) 
for any lit num: metareg fe_X2, wsse(se_X2) bsest(eb) 
* True SD is the square root of tau-squared, ie: <di tau2^.5>
log close

* What do these mean in SD terms? ("X" is the SD on teacher FE term)
* If we move from a teacher at the 25th to a teacher at the 75th percentile, 
* what fraction of an SD do student scores go up? (eg. 7 means 1/7th of an SD)
for num 0.1: di 1/(X*(invnorm(.75)-invnorm(.25)))
* If we move from a teacher at the 10th to a teacher at the 90th percentile, 
* what fraction of an SD do student scores go up? 
for num 0.1: di 1/(X*(invnorm(.9)-invnorm(.1)))

* Kernel density plot of literacy and numeracy gains
*kdensity fe_lit2 if fe_num2<2 & se_num2<2,addplot(kdensity fe_num2 if fe_num2<2 & se_num2<2) title("Figure 1: Dispersion of Teacher Fixed Effects (unadjusted)") legend(label(1 "Literacy") label(2 "Numeracy"))

* How correlated are literacy and numeracy gains?
pwcorr fe_num fe_lit [aw=w_both]
reg fe_num fe_lit [aw=w_both]
la var fe_num "Average change in numeracy (SD)"
la var fe_lit "Average change in literacy (SD)"
set scheme s1mono
*twoway scatter fe_num fe_lit if fe_num<2 & se_num<2 || lfit fe_num fe_lit if fe_num<2 & se_num<2 [aw=w_both], title("Figure 1: Are Good Numeracy Teachers" "Also Good Literacy Teachers?") subtitle("Each dot represents one teacher") legend(off) caption("Numeracy fixed effect = 0.004 + 0.388 * Literacy fixed effect (t=40.5)", ring(0) pos(5) size(medsmall)) xlabel(-2 0 2) ylabel(-2 0 2) yti("Average change in numeracy (SD)")
twoway scatter fe_num fe_lit if fe_num<2 & se_num<2 || lfit fe_num fe_lit if fe_num<2 & se_num<2 [aw=w_both], legend(off) title("") caption("Numeracy fixed effect = 0.004 + 0.388 * Literacy fixed effect (t=40.5)", ring(0) pos(5) size(medsmall)) xlabel(-2 0 2) ylabel(-2 0 2) yti("Average change in numeracy (SD)")
graph export fig1.eps,replace
graph export fig1.wmf,replace
graph save fig1.gph,replace

* Regressions: Including variables singly, then all together
for any masters: reg fe_num X [aw=w_num], r cl(teacher) \ outreg using table_num.doc, coefastr nocons bracket 3aster replace bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any female age : reg fe_num X [aw=w_num], r cl(teacher) \ outreg using table_num.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any experience : reg fe_num X [aw=w_num], r cl(teacher) \ outreg using table_num.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any num: reg fe_X masters female age experience [aw=w_X], r cl(teacher) \ outreg using table_X.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any masters: reg fe_lit X [aw=w_lit], r cl(teacher) \ outreg using table_lit.doc, coefastr nocons bracket 3aster replace bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any female age : reg fe_lit X [aw=w_lit], r cl(teacher) \ outreg using table_lit.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any experience : reg fe_lit X [aw=w_lit], r cl(teacher) \ outreg using table_lit.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any lit: reg fe_X masters female age experience [aw=w_X], r cl(teacher) \ outreg using table_X.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
* Including EQ rating
for any rating: reg fe_num X1 X2 [aw=w_num], r cl(teacher) \ outreg using table_num.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any rating: reg fe_lit X1 X2 [aw=w_lit], r cl(teacher) \ outreg using table_lit.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))
for any num lit: reg fe_X masters female age experience rating1 rating2 [aw=w_X], r cl(teacher) \ outreg using table_X.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se addstat("Teachers",e(N_clust))

* Nonlinear experience - locally weighted regressions
set scheme s1mono
format experience %9.0f
*for any lit: twoway lowess fe_X experience if exp<25,adjust title("Figure 2: Literacy and Teacher Experience") note("Note: Results are based on locally weighted regressions.") legend(off) ytitle("Average change in literacy") xtitle("Years of experience") 
*for any lit: lpoly fe_X experience if experience<=30 [aweight = w_X], noscatter ci legend(off) title("Figure 2: Literacy and Teacher Experience") note("") xtitle("Years of experience")
for any lit: lpoly fe_X experience if experience<=30 [aweight = w_X], noscatter ci legend(off) title("") note("") xtitle("Years of experience")
graph export fig2.eps,replace
graph export fig2.wmf,replace
graph save fig2, replace
*for any num: twoway lowess fe_X experience if exp<25,adjust title("Figure 3: Numeracy and Teacher Experience") note("Note: Results are based on locally weighted regressions.") legend(off) ytitle("Average change in numeracy") xtitle("Years of experience") 
*for any num: lpoly fe_X experience if experience<=30 [aweight = w_X], noscatter ci legend(off) title("Figure 3: Numeracy and Teacher Experience") note("") xtitle("Years of experience")
for any num: lpoly fe_X experience if experience<=30 [aweight = w_X], noscatter ci legend(off) title("") note("") xtitle("Years of experience")
graph save fig3, replace
graph export fig3.eps,replace
graph export fig3.wmf,replace


***********************************************
* Analysing teacher demographics - aggregating to grade-year level
***********************************************
use teacher_fe_gy, clear
gen test_year=year
drop weight
sort teacher
merge teacher using temp_d, nokeep
tab _merge
drop _merge

* Getting a precise estimate of age and experience
* (This exploits the fact that we have precise birth dates and start dates)
gen temp1=doy(date_of_birth)
gen age=test_year-birthyear+(temp1/365)
gen temp2=doy(start_date)
gen experience=test_year-startyear+(temp2/365)
drop temp*
sum experience age, d

* Generating weights
bysort gradeyrschool: egen weight=count(fe_num_gy)
replace weight=1

* Note that fouryr=1 means more than 4-year trained. So we code them as masters=1
replace masters=1 if fouryr==1

* EQ rating
for num 1/2: gen ratingX=1 if rating==X \ recode ratingX .=0 if rating~=.
gen rating34=1 if rating==3 | rating==4 
recode rating34 .=0 if rating~=. 

* Summary statistics (unweighted)
collapse fe_num_gy fe_lit_gy masters female age experience rating1 rating2 rating34 weight [aw=weight], by(gradeyrschool)

* Regressions: Including variables singly, then all together
for any masters: reg fe_num X [aw=weight ], r \ outreg using table_num_gy.doc, coefastr nocons bracket 3aster replace bdec(4) rdec(4) se 
for any female age : reg fe_num X [aw=weight ], r \ outreg using table_num_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any experience : reg fe_num X [aw=weight ], r \ outreg using table_num_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any num: reg fe_X masters female age experience [aw=weight ], r \ outreg using table_X_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any masters: reg fe_lit X [aw=weight ], r \ outreg using table_lit_gy.doc, coefastr nocons bracket 3aster replace bdec(4) rdec(4) se 
for any female age : reg fe_lit X [aw=weight ], r \ outreg using table_lit_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any experience : reg fe_lit X [aw=weight ], r \ outreg using table_lit_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any lit: reg fe_X masters female age experience [aw=weight ], r \ outreg using table_X_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
* Including EQ rating
for any rating: reg fe_num X1 X2 [aw=weight ], r \ outreg using table_num_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any rating: reg fe_lit X1 X2 [aw=weight ], r \ outreg using table_lit_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 
for any num lit: reg fe_X masters female age experience rating1 rating2 [aw=weight ], r \ outreg using table_X_gy.doc, coefastr nocons bracket 3aster append bdec(4) rdec(4) se 

***********************************************
* Looking at gains by SD
***********************************************
cd <DATA DIRECTORY>
use "a qry year 3 results 2001.dta", clear
append using "a qry year 3 results 2002.dta"
append using "a qry year 5 results 2001.dta"
append using "a qry year 5 results 2002.dta"
append using "a qry year 5 results 2003.dta"
append using "a qry year 5 results 2004.dta"
append using "a qry year 7 results 2003.dta"
append using "a qry year 7 results 2004.dta"
for any literacy numeracy: recode X_overall 0=. 
collapse (mean) literacy_mean=literacy_overall numeracy_mean=numeracy_overall (sd) literacy_sd=literacy_overall numeracy_sd=numeracy_overall, by(grade year)
* Mean change
for any literacy numeracy: gen X_change=.
for num 2001/2002: gen temp1=literacy_mean if grade==3 & year==X \ egen temp2=max(temp1) \ replace literacy_change=literacy_mean-temp2 if grade==5 & year==X+2 \ drop temp*
for num 2001/2002: gen temp1=numeracy_mean if grade==3 & year==X \ egen temp2=max(temp1) \ replace numeracy_change=numeracy_mean-temp2 if grade==5 & year==X+2 \ drop temp*
for num 2001/2002: gen temp1=literacy_mean if grade==5 & year==X \ egen temp2=max(temp1) \ replace literacy_change=literacy_mean-temp2 if grade==7 & year==X+2 \ drop temp*
for num 2001/2002: gen temp1=numeracy_mean if grade==5 & year==X \ egen temp2=max(temp1) \ replace numeracy_change=numeracy_mean-temp2 if grade==7 & year==X+2 \ drop temp*
* Mean SD
sum *sd
* Mean change
sum *change