Solución Tarea 2

Author

Emiliano Ramírez y Mauricio Romero

Published

September 24, 2023

Code

rm(list = ls())
library(pacman)
p_load(tidyverse, data.table, readr, tidylog, readxl,
       skimr, Hmisc, lubridate, foreign, ggtext,
       gridExtra, modelsummary, sf, stringi,
       did, haven, fixest
       )

# tidyverse for obious reasons
# data.table
# readr for reading files easily
# tidylog for captions about data wrangling commands
# readxl for reading excel files
# skimr for skim_without_charts function
# Hmisc for describe function 
# lubridate for date processing 
# foreign for dbf databases reading 
# ggtext for text format for Gini plot 
# modelsummary for pretty tables
# sf for shapefiles processing
# stringi for accent string processing 
# did for difference in differences functions (from Callaway and Sant'Anna)
# haven for reading dta files (dataframe format from Stata)
# fixest for feols function for fixed effects regressions

Problema 1

1. Preguntas de incisos.

a)

Un mérito importante de optar por este diseño es que el tamaño de muestra para detectar un efecto de tamaño razonable es menor pues el tratamiento es asignado directamente a las unidades. Un problema importante que surge de este tipo de asignación es que se corre un gran riesgo de que SUTVA no pueda cumplirse ya que el hecho de que maestros de una misma escuela tenga unos tratamiento y otros no puede generar cambios de comportamiento en los maestros y, por ende, el tratamiento tendría ingerencia no solo en las unidades asignadas a él si no en unidades que no fueron asignadas al tratamiento. Es decir, habría externalidades del tratamiento.

b)

Si clusterizamos el asignamiento del tratamiento a nivel escuela podríamos no preocuparnos tanto de que no se cumpla SUTVA. No obstante, para que podamos detectar un efecto de tamaño razonable será necesario que el número de clusters aumente y que el número de observaciones por cluster también lo haga, lo que implica una mayor complejidad en el diseño de la intervención.

c)

Generalmente en este tipo de intervenciones, suena más razonable optar por un diseño de asignación del tratamiento a clusters ya que es importante que, si queremos obtener estiamciones causales, no sea un problema suponer que se cumple SUTVA.

d)

El hecho de tener una muestra más grande de escuelas implica que se puede tener un mayor número de clusters y, por lo tanto, esto nos beneficiaría en terminos de poder estadístico aumentándolo, todo lo demás constante.

e)

El parámetro indicado para calibrar y cumplir con las condiciones del gobierno federal es el mínimo efecto detectable y este depende de parámetros como el número de observaciones, la varianza de la variable de respuesta, el nivel de signifacncia de la prueba, y el poder de la prueba que se desee y la fracción de la población que desea sea asignada al tratamiento.

f)

Una mayor varianza poblaciónal de la variable de respuesta implica un menor poder estadístico de la prueba en cuestión.

g)

2. Ejercicio práctico

El poder estadístico se maximiza cuando el tamaño de los grupos es uniforme.

Primero identificamos las bases de nivel estudiante, nivel escuela y de los puntjes de la prueba estandarizada que se realizó en 2014. Estas bases son las sigueintes:

School.dta (donde encontramos la información de la asignación del tratamiento)
Student_PSLE_2014.dta (donde están los puntajes de las escuelas estandarizadas)
Student_School_House_Teacher_Char.dta (donde se pueden encontrar algunos controles para la especificación deseada que se escoja).

Code

school <- read_dta("tarea2/p1/raw_data/QJE_Replication/Data/School.dta")

puntajes_endline <- read_dta("tarea2/p1/raw_data/QJE_Replication/Data/Student_PSLE_2014.dta")

puntajes_baseline <- read_dta("tarea2/p1/raw_data/QJE_Replication/Data/Student_PSLE_2013.dta")

Seleccionamos las variables que nos servirán

Code

school <- school %>% 
  select(SchoolID, treatment, treatarm, DistrictID, factyn9_T7, tbgdno_math1_T7,
         tbgdno_math2_T7, tbgdno_math3_T7, tbgdno_math4_T7, tbgdno_math5_T7,
         tbgdno_math6_T7, tbgdno_math7_T7, s1841_T7, s1851_T7, s1842_T7, s1852_T7,
         s1843_T7, s1853_T7, s1844_T7, s1854_T7, s1845_T7, s1855_T7, s1846_T7,
         s1856_T7, s1847_T7, s1857_T7, studentsTeacherRatio_T7, s108_T1) %>%
  rename(urban=s108_T1,
         library=factyn9_T7)

#Create school controls
school <- school %>%
  mutate(treatment_1=ifelse(treatment=='COD', 1,0),
         treatment_2=ifelse(treatment=="CG", 1,0),
         treatment_3=ifelse(treatment=="Both",1,0),
         control=ifelse(treatment=="Control",1,0),
         log_mathbooks=log(tbgdno_math1_T7+
         tbgdno_math2_T7+ tbgdno_math3_T7+ tbgdno_math4_T7+ tbgdno_math5_T7+
         tbgdno_math6_T7+ tbgdno_math7_T7),
         num_boys=s1841_T7+s1842_T7+s1843_T7+ s1844_T7+s1845_T7+ s1846_T7+s1847_T7,
         num_girls=s1851_T7+s1852_T7+s1853_T7+s1854_T7+s1855_T7+s1856_T7+s1857_T7,
         ratio_girls=num_girls/(num_boys+num_girls)
         ) %>%
  select(SchoolID, DistrictID, treatment_1, treatment_2, treatment_3, control, log_mathbooks,
         num_boys, num_girls, ratio_girls, urban, library, studentsTeacherRatio_T7)
  

puntajes_endline <- puntajes_endline %>% 
  select(CAND, SX, AVERAGE, SchoolID) %>%
  rename(SEX=SX, 
         AVERAGE_ENDLINE=AVERAGE) %>% 
  mutate(SEX=ifelse(SEX=='F', 1,0),
         AVERAGE_ENDLINE_STD=(AVERAGE_ENDLINE - mean(AVERAGE_ENDLINE))/sd(AVERAGE_ENDLINE)) 



puntajes_baseline <- puntajes_baseline %>% 
  select(CAND, AVERAGE) %>%
  rename(AVERAGE_BASELINE=AVERAGE) %>%
  mutate(AVERAGE_BASELINE_STD=(AVERAGE_BASELINE - mean(AVERAGE_BASELINE))/sd(AVERAGE_BASELINE))

Unimos base de datos para estimaciones

Code

puntajes <- inner_join(puntajes_endline, puntajes_baseline, by=c("CAND"))

puntajes <- inner_join(puntajes, school, by=c("SchoolID"))

Estimamos por MCO la siguiente especificación

\[ Z_{i s d t}= \delta_0+\delta_1 \text {Treatment 1}_s+\delta_2 \text {Treatment 2}_s+\delta_3 \text {Treatment 3}_s +\gamma_z Z_{i s d, t=0}+\gamma_d+X_i \delta_4+X_s \delta_5+\varepsilon_{i s d t} \]

donde \(Z_{isd,t}\) es el puntaje estandarizado del examen baseline que se llevo a acabo antes de la intervención, \(\gamma_d\) es el efecto fijo del distrito de la escuela, \(X_i\) es el género de la unidad y \(X_s\) son los controles a nivel escuela que escogimos.

Code

p1_MCO <- feols(AVERAGE_ENDLINE_STD ~ treatment_1 + treatment_2 + treatment_3 + 
                  AVERAGE_BASELINE_STD + SEX  + 
                  ratio_girls + urban + library +
                  studentsTeacherRatio_T7| DistrictID,
  cluster = "DistrictID",
  data = puntajes
)

f <- function(x) format(round(as.numeric(x),0), nsmall=0, big.mark=",")  # formato de coma

gm <- list(
  list("raw" = "nobs", "clean" = "N", "fmt" = f),
  list("raw" = "FE: DistrictID", "clean" = "FE: Distrito", "fmt"= f),
  list("raw" = "vcov.type", "clean" = "Clustered error", "fmt"=f))

#reportamos con modelsummary la regresión estimada
modelsummary(p1_MCO, stars=TRUE, statistic = "{std.error} ({p.value})",
             coef_omit = c(-1,-2,-3), gof_map = gm)

	(1)
treatment_1	0.022
	0.074 (0.773)
treatment_2	−0.002
	0.058 (0.977)
treatment_3	0.184*
	0.067 (0.025)
N	12,803
FE: Distrito	X
Clustered error	by: DistrictID
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Ahora estimemos por inferencia aleatorizada. Escogeremos como test statistic a evaluar la diferencia de medias simple.

Para ello veamos si existen outliers en la distribución de puntajes del examen estandarizado endline. La densidad se ve bien y no hay colas muy pesadas por lo que no se aprecian outliers y podemos usar como test statistic la diferencia de medias simple.

Code

ggplot(puntajes, aes(x=AVERAGE_ENDLINE_STD))+
  geom_density(color="darkblue", fill="lightblue") +
  geom_vline(aes(xintercept=mean(AVERAGE_ENDLINE_STD)),
            color="blue", linetype="dashed", size=1) +
  annotate(x=mean(puntajes$AVERAGE_ENDLINE_STD),y=+Inf,
           label=paste('media ',round(mean(puntajes$AVERAGE_ENDLINE_STD),2)),
           vjust=10, hjust=1.2, geom="text", size=3) +
  labs(x='Puntaje promedio estandarizado') +
  theme_minimal()

Para efecto de tratamiento 1

Code

Y <- as.numeric(puntajes$AVERAGE_ENDLINE_STD)

treatment_1 <- as.numeric(puntajes$treatment_1)
treatment_2 <- as.numeric(puntajes$treatment_2)
treatment_3 <- as.numeric(puntajes$treatment_3)


real_diff_1=mean(Y[treatment_1==1])-mean(Y[treatment_1==0]) 
placebo_diff_1=NULL

real_diff_2=mean(Y[treatment_2==1])-mean(Y[treatment_2==0]) 
placebo_diff_2=NULL

real_diff_3=mean(Y[treatment_3==1])-mean(Y[treatment_3==0]) 
placebo_diff_3=NULL

num_rep=30000

#for treatment effect 1 
for(r in 1:num_rep){
  placebo_treatment1=sample(treatment_1) 
  placebo_diff_1=c(placebo_diff_1, mean(Y[placebo_treatment1==1])-mean(Y[placebo_treatment1==0]))
}

pvalue1=mean(placebo_diff_1 > real_diff_1) 

hist(placebo_diff_1 ,freq=F,las=1,xlab="Difference",
     main=paste("P-value for treatment 1 effect",format(pvalue1 ,digits=2))) 
abline(v=real_diff_1 ,col=rgb(1,0,0,1),lwd=3,lty=2)

Para efecto de tratamiento 2

Code

#for treatment effect 2 
for(r in 1:num_rep){
  placebo_treatment2=sample(treatment_2) 
  placebo_diff_2=c(placebo_diff_2, mean(Y[placebo_treatment2==1])-mean(Y[placebo_treatment2==0]))
}

pvalue2=mean(placebo_diff_2 > real_diff_2) 

hist(placebo_diff_2 ,freq=F,las=1,xlab="Difference",
     main=paste("P-value for treatment 2 effect",format(pvalue2 ,digits=2))) 
abline(v=real_diff_2 ,col=rgb(1,0,0,1),lwd=3,lty=2)

Para efecto de tratamiento 3

Code

#for treatment effect 3 
for(r in 1:num_rep){
  placebo_treatment3=sample(treatment_3) 
  placebo_diff_3=c(placebo_diff_3, mean(Y[placebo_treatment3==1])-mean(Y[placebo_treatment3==0]))
}

pvalue3=mean(placebo_diff_3 > real_diff_3) 

hist(placebo_diff_3 ,freq=F,las=1,xlab="Difference",
     main=paste("P-value for treatment 3 effect",format(pvalue3 ,digits=2))) 
abline(v=real_diff_3 ,col=rgb(1,0,0,1),lwd=3,lty=2)

Problema 2

Primero leemos los datos que se nos proporcionaron

Code

df <- haven::read_dta("tarea2/p2/raw_data/medicaid_did.dta")

Non-staggered DiD

Procesamos los datos como se nos indica en el taller para realizar el ejercicio. Esto nos permitirá tener puras observaciones tratadas en el año 2014 o no tratadas en lo absoluto.

Code

# tiramos años estrictamente mayores a 2015 y tiramos los tres estados que expandieron la ayuda en 2015
df_p1did <- df %>%
  filter(year <= 2015) %>%
  filter(is.na(yexp2) | yexp2 < 2015)

1. Estimamos los coeficintes de la regresión two way fixed effects para esta submuestra

Expresemos la especificación de la regresión que nos dará estos coeficientes.

\[ Y_{it} = \alpha_i + \lambda_t + \sum_{s \neq 2013} \mathbb{I}_{\{s=t\}} \times D_i \times \beta_s + u_{it} \]

Donde \(Y_{it}\) es la variable de respuesta que utilizaremos de la unidad \(i\) en el tiempo \(t\) , en este caso dins, \(\alpha_i\) es el efecto fijo del estado, \(\lambda_t\) es el efecto fijo del año (ya tenemos los dos efectos fijos clásicos en nuestra especificación), \(\mathbb{I}_{\{s \neq 2013\}}\) es una indicadora que vale uno cuando el año es igual a \(t\) (note que la suma se salta el año 2013, esto es porque estamos tomando como año de referencia al año 2013, por lo que no se estima su coeficiente para evitar colinealidad), \(D_i\) es la indicadora que vale uno cuando la unidad \(i\) fue tratada en el año 2014, \(\beta_s\) es el coeficiente que queremos estimar para el año \(s\) y \(u_{it}\) es el componente no explicado de la unidad \(i\) en el tiempo \(t\).

Recordemos que los coeficientes estimados \(\hat{\beta}\) de TWFE son equivalentes a las estimaciones DiD entre unidades tratadas y no tratadas entre los períodos \(s\) y 2013.

Code

#Creemos la variable indicadora del tratamiento en 2014
df_p1did <- df_p1did %>%
  mutate(D = ifelse(is.na(yexp2), 0, 1))

Estimemos nuestro modelo que especificamos arriba con la función feols de fixest. Clusterizamos a nivel estado los errores estandar.

Code

reg_nonStag <- feols(dins ~ i(year, D, ref = 2013) | stfips + year,
  cluster = "stfips",
  data = df_p1did
)

#obtenemos nombres de coeficientes y estadisticos

f <- function(x) format(round(as.numeric(x),0), nsmall=0, big.mark=",")  # formato de coma

gm <- list(
  list("raw" = "nobs", "clean" = "N", "fmt" = f),
  list("raw" = "FE: stfips", "clean" = "FE: estado", "fmt"= f),
  list("raw" = "FE: year", "clean" = "FE: año", "fmt"= f),
  list("raw" = "vcov.type", "clean" = "Clustered error", "fmt"=f))

#reportamos con modelsummary la regresión estimada
modelsummary(reg_nonStag, stars=TRUE, gof_map = gm)

	(1)
year = 2008 × D	−0.010
	(0.008)
year = 2009 × D	−0.013+
	(0.007)
year = 2010 × D	−0.002
	(0.007)
year = 2011 × D	−0.006
	(0.007)
year = 2012 × D	−0.006
	(0.006)
year = 2014 × D	0.042***
	(0.008)
year = 2015 × D	0.069***
	(0.011)
N	304
FE: estado	X
FE: año	X
Clustered error	by: stfips
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

2. Leads & lags

Obtenemos plot de leads & lags

Code

iplot(reg_nonStag, xlab = "año", ylab="Estimación e IC al 95%", main='Event study: Expansión Medicaid en 2014')

Staggered DiD

1. Estimación de ATT(g,t)

La especificación del modelo sería la sigueinte

\[ Y=\alpha_1^{g, t}+\alpha_2^{g, t} \cdot G_g+\alpha_3^{g, t} \cdot 1\{T=t\}+\beta^{g, t} \cdot\left(G_g \times 1\{T=t\}\right)+\epsilon^{g, t} \]

Donde \(\alpha_{1}^{g,t}\) es el efecto fijo de la unidad que pertenece al grupo que fue tratado en el tiempo \(g\) observada en el tiempo \(t\), \(G_g\) es una indicadora de que vale 1 si la unidad fue tratada en tiempo \(g\) y \(\epsilon^{g,t}\) es el termino de error. el ATT(g,t) corresponde al coeficiente \(\beta^{g,t}\).

Estimación cuando unidades de control todavía no tratados.

Code

reg_stag1 <- att_gt(
  yname = "dins",
  tname = "year",
  idname = "stfips",
  gname = "yexp2",
  data = df %>% mutate(yexp2 = ifelse(is.na(yexp2), 3000, yexp2)),
  control_group = "notyettreated",
  allow_unbalanced_panel = TRUE,
  print_details = FALSE,
  alp = 0.05
)

modelsummary(reg_stag1)

	(1)
ATT(2014,2009)	−0.007
	(0.005)
ATT(2014,2010)	0.011
	(0.006)
ATT(2014,2011)	0.002
	(0.005)
ATT(2014,2012)	0.002
	(0.006)
ATT(2014,2013)	0.001
	(0.007)
ATT(2014,2014)	0.047
	(0.008)
ATT(2014,2015)	0.069
	(0.010)
ATT(2014,2016)	0.079
	(0.011)
ATT(2014,2017)	0.073
	(0.011)
ATT(2014,2018)	0.074
	(0.013)
ATT(2014,2019)	0.080
	(0.011)
ATT(2015,2009)	0.007
	(0.013)
ATT(2015,2010)	−0.024
	(0.011)
ATT(2015,2011)	−0.003
	(0.006)
ATT(2015,2012)	0.000
	(0.003)
ATT(2015,2013)	−0.010
	(0.011)
ATT(2015,2014)	−0.002
	(0.007)
ATT(2015,2015)	0.049
	(0.022)
ATT(2015,2016)	0.053
	(0.015)
ATT(2015,2017)	0.068
	(0.010)
ATT(2015,2018)	0.066
	(0.013)
ATT(2015,2019)	0.074
	(0.012)
ATT(2016,2009)	0.002
	(0.006)
ATT(2016,2010)	0.051
	(0.008)
ATT(2016,2011)	−0.026
	(0.039)
ATT(2016,2012)	−0.028
	(0.045)
ATT(2016,2013)	0.044
	(0.045)
ATT(2016,2014)	−0.032
	(0.049)
ATT(2016,2015)	0.041
	(0.019)
ATT(2016,2016)	0.032
	(0.010)
ATT(2016,2017)	0.037
	(0.009)
ATT(2016,2018)	0.065
	(0.013)
ATT(2016,2019)	0.083
	(0.009)
ATT(2017,2009)	0.010
	(0.003)
ATT(2017,2010)	−0.016
	(0.003)
ATT(2017,2011)	0.002
	(0.003)
ATT(2017,2012)	0.018
	(0.003)
ATT(2017,2013)	0.001
	(0.004)
ATT(2017,2014)	−0.006
	(0.006)
ATT(2017,2015)	0.004
	(0.006)
ATT(2017,2016)	0.040
	(0.006)
ATT(2017,2017)	0.047
	(0.005)
ATT(2017,2018)	0.068
	(0.005)
ATT(2017,2019)	0.065
	(0.004)
ATT(2019,2009)	0.019
	(0.005)
ATT(2019,2010)	−0.027
	(0.009)
ATT(2019,2011)	−0.024
	(0.003)
ATT(2019,2012)	0.002
	(0.025)
ATT(2019,2013)	0.013
	(0.007)
ATT(2019,2014)	0.000
	(0.007)
ATT(2019,2015)	−0.011
	(0.009)
ATT(2019,2016)	−0.018
	(0.028)
ATT(2019,2017)	0.009
	(0.019)
ATT(2019,2018)	0.006
	(0.006)
ATT(2019,2019)	0.037
	(0.005)
Num.Obs.	46
Std.Errors	by: stfips
ngroup	5
ntime	12
control.group	notyettreated
est.method	dr

Code

did::ggdid(reg_stag1) + ggtitle('Event Study por Grupo | not yet treated') + 
labs(subtitle = "Estados expanden Medicaid en diferentes años") + 
  theme(plot.title = element_text(color = "black", hjust=0.5),
        plot.subtitle=element_text(hjust=0.5)) + 
  scale_x_discrete(name ="Periodos", limits=c("2009","2010", "2011", "2012", "2013",
                                              "2014", "2015","2016", "2017", "2018", "2019"))

Estimación cuando unidades de control como nunca tratadas.

Code

reg_stag2 <- att_gt(
  yname = "dins",
  tname = "year",
  idname = "stfips",
  gname = "yexp2",
  data = df %>% mutate(yexp2 = ifelse(is.na(yexp2), 3000, yexp2)),
  control_group = "nevertreated",
  allow_unbalanced_panel = TRUE,
  print_details = FALSE,
  alp = 0.05
)

modelsummary(reg_stag2)

	(1)
ATT(2014,2009)	−0.004
	(0.006)
ATT(2014,2010)	0.011
	(0.007)
ATT(2014,2011)	−0.005
	(0.005)
ATT(2014,2012)	0.000
	(0.006)
ATT(2014,2013)	0.006
	(0.006)
ATT(2014,2014)	0.042
	(0.008)
ATT(2014,2015)	0.069
	(0.010)
ATT(2014,2016)	0.078
	(0.010)
ATT(2014,2017)	0.071
	(0.011)
ATT(2014,2018)	0.073
	(0.012)
ATT(2014,2019)	0.080
	(0.010)
ATT(2015,2009)	0.006
	(0.013)
ATT(2015,2010)	−0.017
	(0.012)
ATT(2015,2011)	−0.008
	(0.006)
ATT(2015,2012)	0.000
	(0.005)
ATT(2015,2013)	−0.004
	(0.011)
ATT(2015,2014)	−0.006
	(0.006)
ATT(2015,2015)	0.053
	(0.023)
ATT(2015,2016)	0.053
	(0.014)
ATT(2015,2017)	0.067
	(0.009)
ATT(2015,2018)	0.066
	(0.012)
ATT(2015,2019)	0.074
	(0.011)
ATT(2016,2009)	0.002
	(0.008)
ATT(2016,2010)	0.054
	(0.010)
ATT(2016,2011)	−0.030
	(0.064)
ATT(2016,2012)	−0.027
	(0.075)
ATT(2016,2013)	0.048
	(0.077)
ATT(2016,2014)	−0.034
	(0.080)
ATT(2016,2015)	0.041
	(0.019)
ATT(2016,2016)	0.032
	(0.010)
ATT(2016,2017)	0.036
	(0.009)
ATT(2016,2018)	0.064
	(0.014)
ATT(2016,2019)	0.083
	(0.010)
ATT(2017,2009)	0.010
	(0.005)
ATT(2017,2010)	−0.010
	(0.007)
ATT(2017,2011)	−0.004
	(0.004)
ATT(2017,2012)	0.017
	(0.005)
ATT(2017,2013)	0.007
	(0.005)
ATT(2017,2014)	−0.010
	(0.005)
ATT(2017,2015)	0.007
	(0.007)
ATT(2017,2016)	0.038
	(0.006)
ATT(2017,2017)	0.048
	(0.005)
ATT(2017,2018)	0.070
	(0.005)
ATT(2017,2019)	0.065
	(0.004)
ATT(2019,2009)	0.018
	(0.007)
ATT(2019,2010)	−0.020
	(0.010)
ATT(2019,2011)	−0.029
	(0.004)
ATT(2019,2012)	0.002
	(0.016)
ATT(2019,2013)	0.018
	(0.008)
ATT(2019,2014)	−0.004
	(0.006)
ATT(2019,2015)	−0.006
	(0.009)
ATT(2019,2016)	−0.016
	(0.021)
ATT(2019,2017)	0.009
	(0.017)
ATT(2019,2018)	0.006
	(0.006)
ATT(2019,2019)	0.037
	(0.006)
Num.Obs.	46
Std.Errors	by: stfips
ngroup	5
ntime	12
control.group	nevertreated
est.method	dr

Code

did::ggdid(reg_stag2) + ggtitle('Event Study por Grupo | never treated') + 
labs(subtitle = "Estados expanden Medicaid en diferentes años") + 
  theme(plot.title = element_text(color = "black", hjust=0.5),
        plot.subtitle=element_text(hjust=0.5)) + 
  scale_x_discrete(name ="Periodos", limits=c("2009","2010", "2011", "2012", "2013",
                                              "2014", "2015","2016", "2017", "2018", "2019"))

2. Construyamos a mano los coeficientes ATT(2015, 2013), ATT(2013, 2013) y ATT(2016, 2014).

Si se revisa con los coeficientes de summary(notyettreated) son los mismos.

ATT(2015, 2013)

Code

#ATT(2015,2013)

#Creamos un df auxiliar con las medias condicionales de la var. de respuesta de 2012 y 2013 
# Para ello nos quedamos con observaciones de los años 2012 y 2013 y 
# creamos dummy de si la unidad fue tratada en 2015.

ATT_15_13 <- df %>%
  filter(year %in% c(2012, 2013)) %>%
  mutate(treated = case_when(
    yexp2 == 2015 ~ 1,
    is.na(yexp2) | yexp2 != 2015 ~ 0
  )) %>%
  group_by(treated, year) %>%
  summarise(dins = mean(dins)) %>% ungroup() 

round(with(ATT_15_13,
  {(dins[year == 2013 & treated == 1] - dins[year == 2012 & treated == 1]) -
    (dins[year == 2013 & treated == 0] - dins[year == 2012 & treated == 0])}
  ),4)

[1] -0.01

ATT(2017, 2017)

Code

ATT_17_17 <- df %>%
  filter(year %in% c(2016, 2017)) %>%
  mutate(treated = case_when(
    yexp2 == 2017 ~ 1,
    is.na(yexp2) | yexp2 > 2017 ~ 0
  )) %>%
  group_by(treated, year) %>%
  summarise(dins = mean(dins)) %>% ungroup() %>%
  filter(!is.na(treated))

round(with(ATT_17_17,
  {(dins[year == 2017 & treated == 1] - dins[year == 2016 & treated == 1]) -
    (dins[year == 2017 & treated == 0] - dins[year == 2016 & treated == 0])}
  ),4)

[1] 0.0471

ATT(2016, 2014)

Code

ATT_16_14 <- df %>%
  filter(year %in% c(2013, 2014)) %>%
  mutate(treated = case_when(
    yexp2 == 2016 ~ 1,
    is.na(yexp2) | yexp2 > 2016 ~ 0
  )) %>%
  group_by(treated, year) %>%
  summarise(dins = mean(dins)) %>% ungroup() %>%
  filter(!is.na(treated))

round(with(ATT_16_14,
  {(dins[year == 2014 & treated == 1] - dins[year == 2013 & treated == 1]) -
    (dins[year == 2014 & treated == 0] - dins[year == 2013 & treated == 0])}
  ),4)

[1] -0.0329

3. Event study de agregación dinámica.

La agregación dinámica está basada en la siguiente especificación TWFE:

\[ Y_{i, t}=\alpha_t+\alpha_g+\sum_{e=-K}^{-2} \delta_e^{\text {anticip }} \cdot D_{i, t}^e+\sum_{e=0}^L \beta_e \cdot D_{i, t}^e+v_{i, t} \]

Donde \(D_{i,t}^e = \mathbb{I}_{\{ t- G_i = e \}}\) es una indicadora que vale 1 cuando la unidad \(i\) está desfasada \(e\) períodos del tratamiento inicial en el tiempo \(t\). Se suele enfocarse en los coeficientes \(\beta_e, \, e \geq0\) pues estos parámetros son típicamente interpretados como la medición del efecto de participar en el tratamiento después de \(e\) períodos de tiempo del tratamiento inicial.

Code

dynamic_gg <- aggte(reg_stag1,
  type = "dynamic",
  min_e = -5, max_e = 5
)

ggdid(dynamic_gg) + ggtitle('Event Study Agregación Dinámica') + 
  theme(plot.title = element_text(color = "black", hjust=0.5),
        plot.subtitle=element_text(hjust=0.5)) + 
  scale_x_discrete(name ="Períodos", limits=seq(-5,5,1) )

4. Estimación de regresión por MCO.

La siguiente especificación

\[ Y_{it}=\alpha_i + \lambda_t + D_{it}\beta + \epsilon_it \]

es la versión estática del TWFE. Con ella se puede hacer la agregación simple que toma un promedio ponderado de los ATT(g,t), ponderando proporcionalmente por tamaño de grupo.

Estimemos la regresión con MCO. Para ello tenemos crear una dummy que indique si el período es después del tratamiento.

Code

# Creamos variable post tratamiento
df <- df %>% mutate(post_treated = !is.na(yexp2) & year >= yexp2)

# Corremos TWFE estático con errores estándar clusterizados a nivel estado
reg_static <- feols(dins ~ post_treated | stfips + year, data = df, cluster = "stfips")

gm <- list(
  list("raw" = "nobs", "clean" = "N", "fmt" = f),
  list("raw" = "FE: stfips", "clean" = "FE: estado", "fmt"= f),
  list("raw" = "FE: year", "clean" = "FE: año", "fmt"= f),
  list("raw" = "vcov.type", "clean" = "Clustered error", "fmt"=f))

#reportamos con modelsummary la regresión estimada
modelsummary(reg_static, stars=TRUE, gof_map = gm, coef_rename = c("post_treatedTRUE" = "post treated beta"))

	(1)
post treated beta	0.070***
	(0.007)
N	552
FE: estado	X
FE: año	X
Clustered error	by: stfips
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Observemos cómo es el coeficiente de la agregación simple.

Code

aggte(reg_stag1, type = "simple")


Call:
aggte(MP = reg_stag1, type = "simple")

Reference: Callaway, Brantly and Pedro H.C. Sant'Anna.  "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015> 


   ATT    Std. Error     [ 95%  Conf. Int.]  
 0.068        0.0074     0.0534      0.0826 *


---
Signif. codes: `*' confidence band does not cover 0

Control Group:  Not Yet Treated,  Anticipation Periods:  0
Estimation Method:  Doubly Robust

De hecho, son muy parecidos. Para explicar esto se pude usar la descomposición de Bacon que nos arroja los pesos con los cuales son los pesos que el TWFE pone en las comparaciones con unidades que nunca son tratadas y con las que todavía no son tratadas.