**********************************
**********************************
**********************************
* Part 1: PS Estimation
**********************************
**********************************
**********************************

** Install 
search pscore
search pbalchk
search dr
ssc install psmatch2


clear

** Change to your directory! **
cd "/Users/albertoguzman-alvarez/Desktop/Workshop"
use "shadishclark.dta", clear
set more off

** Store the covariates in a global	
global covariates "vocabpre numbmath likemath likelit preflit pagree pconsc pemot pintell mars beck cauc afram other age male married momdegr daddegr credit majormi actcomp hsgpaar collgpaa mathpre pextra"

** Global for continuous variables
global c_cov "vocabpre numbmath likemath likelit pagree pconsc pemot pintell mars beck age credit actcomp hsgpaar collgpaa mathpre pextra"


** Create treatment indicator
codebook vm
gen treatment = . 
	replace treatment = 1 if vm == 1
	replace treatment = 0 if vm == 2
	
	
**********************************
* Assessing Imbalance *
**********************************
help pbalchk


* Check covariate balance
*standardized mean difference, measure of difference between group1 and group2, standardized so that 
*the difference is measurement independent ( < .2 considered good)

pbalchk treatment $covariates, nostandardize
pbalchk treatment $covariates


// ** Covariate balance checks	
// foreach i of global covariates {
// 	reg `i' treatment 
// 	estimates store `i', title (`i')
// 	}

// foreach i of global covariates  {
// 	outreg2 [`i'] using covariate_balance_pre.xls, dec(2) rdec(3) append cttop(full) alpha(0.10, 0.05,    0.01, 0.001) ///
// 	symbol(~, *, **, ***) keep(treatment) ctitle(`var') label
// }
	

** Look at distribution of continuous covariates by treatment condition**
foreach var of varlist $c_cov {		
	graph twoway ///
		(kdensity `var' if treatment  == 1) ///
		(kdensity `var' if treatment  == 0), ///
		legend(label(1 "Treatment") label(2 "Control")) ///
		ytitle(Kernel density) ///
		title(`var') ///
		saving(`var'.gph, replace)
}
	
graph combine ///
vocabpre.gph numbmath.gph likemath.gph likelit.gph ///
pagree.gph   pconsc.gph   pemot.gph    pintell.gph ///
mars.gph     beck.gph     age.gph      credit.gph  ///
actcomp.gph  hsgpaar.gph  collgpaa.gph mathpre.gph ///
pextra.gph, ///
title("Kernel Densities by Treatment Condition") ///
saving(k_pre.gph, replace)

**********************************
* Calculate PS *
**********************************

*** calculate propensity score
help pscore

pscore treatment $covariates, pscore(pscore) 

*region of common support
psgraph, treated(treatment)pscore(pscore)

* Iterate AND/OR Trim
drop pscore 
gen numbmath_2 = numbmath^2
pscore treatment $covariates, pscore(pscore) 

* balance after 
psgraph, treated(treatment)pscore(pscore)

**********************************
* Stratify *
**********************************

* Stratify ps score or use pscore blocks
xtile pq = pscore, n(3)
tab pq treatment, ro

* Check balance
pbalchk treatment $covariates, strata(pq) 

* Naive estimate
reg mathall treatment 

* Regress with strata 
help xi

xi: reg mathall treatment i.pq



**********************************
* IPW *
**********************************
 
* Generate IPW
qui dr mathall treatment $covariates, genvars
br iptwt 

* Check balance 
pbalchk treatment $covariates, wt(iptwt)

* Regress with IPW
reg mathall treatment [pweight = iptwt]

**********************************
* Greedy Matching *
**********************************

help psmatch2

**********************************
*matching 1:1 Nearest neighbor
**********************************

psmatch2 treatment, ///
	outcome(mathall) pscore(pscore) neighbor(1) ate

*test balance
pstest $covariates, both graph
pstest $covariates, both hist

*treatment effect, psmatch2 doesnt account for estimation of PS, fine tune above
teffects psmatch (mathall) (treatment $covariates), nneighbor(1) ate

**********************************
*matching 1:1 Nearest neighbor with caliper
**********************************

psmatch2 treatment, ///
	outcome(mathall) pscore(pscore) caliper(.1) neighbor(1) ate

*test balance
pstest $covariates, both graph
pstest $covariates, both hist

*treatment effect
teffects psmatch (mathall) (treatment $covariates), nneighbor(1) caliper(.1) ate

**********************************
*matching 1:M with caliper & replacement
**********************************

psmatch2 treatment, ///
	outcome(mathall) pscore(pscore) caliper(.1) neighbor(2) ate

*test balance
pstest $covariates, both graph
pstest $covariates, both hist

*treatment effect
teffects psmatch (mathall) (treatment $covariates), nneighbor(2) caliper(.1) ate

*compare to naive
reg mathall treatment 
reg mathall treatment $covariates