```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r,warning = F,message = F}
library(dplyr)
library(ggplot2)
library(plm)
library(did)
```

#######################################################################################################

In this part of the exercise we show examples of how the basic difference-in-differences estimator that was discussed in the lecture can be implemented in R. This is done using a single simulated dataset, for demonstration purposes. We show both the standard approach (using base R) and the more contemporary approach using the did-package. In a first step, let's load the dataset.

## 1. Load the dataset

```{r}
# LOAD DATASET
df <- readRDS("simdata.rds")
```


## 2. Case 1: The 2x2-DiD

Let's assume that the intervention happened after period 7, and for simplicity, let's only consider the 2 x 2 case, which means that we only have one pre- (7) and one post-treatment period (8). As always, we have two ways of estimating the treatment effect: (i) simple difference in means, and (ii) regression.

```{r}

# SUBSET DATA TO TWO PERIODS (PRE & POST)
temp <- df %>% subset(., t %in% c(7, 8)) %>% arrange(id, t)

# 1. SIMPLE DIFFERENCE IN MEANS
pre_treated <- temp %>% subset(., t == 7) %>% subset(., g == 1) %>% pull(y) %>% mean(., na.rm = T)
pre_untreated <- temp %>% subset(., t == 7)  %>% subset(., g == 0) %>% pull(y) %>% mean(., na.rm = T)
post_treated <- temp %>% subset(., t == 8)  %>% subset(., g == 1) %>% pull(y) %>% mean(., na.rm = T)
post_untreated <- temp %>% subset(., t == 8)  %>% subset(., g == 0) %>% pull(y) %>% mean(., na.rm = T)

(pre_untreated-post_untreated)-(pre_treated-post_treated)

###

# 2. REGRESSION
temp <- temp %>% mutate(t8 = as.numeric(t == 8))
reg <- lm(y ~ g*t8, data = temp)
summary(reg)

reg <- lm(y ~ g + t8 + g:t8, data = temp)
summary(reg)

reg <- lm(y ~ g + t8 + d1, data = temp)
summary(reg)
```

Consider then the same estimation using a fixed-effects model with fixed effects for the individual units (and times). This does require panel data for (at least some) units.

```{r}
## ESTIMATION WITH TWO-WAY FIXED EFFECTS

# A) WITH BASE R
reg <- lm(y ~ d1 + factor(id) + factor(t), data = temp)
summary(reg)[['coefficients']][1:10,]

# B) WITH PLM PACKAGE
reg <- plm(y ~ d1 + factor(t), data = temp, index = c("id", "t"), model = "within", effect = "individual")
summary(reg)
```


## 3. Case 2: DiD with multiple pre- and post-periods

In addition to these 2x2-DiD setup, having more pre- and post-treatment periods helps us to assess whether we find evidence that the parallel trend assumption is likely to be valid. First, we simply plot the outcomes for both treated and untreated groups over time. For this, we use ggplot.

```{r}
plot_df <- df %>% 
  group_by(t, g, ) %>% summarize(y = mean(y, na.rm = T)) %>% ungroup() %>%
  mutate(group = if_else(g == 0, "Control", "Treatment"))

ggplot(plot_df, aes(x = t, y = y, group = group)) +
  geom_vline(xintercept = 7.5, color = 'black') +
  geom_segment(x = 6.4, y = 20, xend = 7.4, yend = 20,
               arrow = arrow(length = unit(.1, "inches"))) +
  annotate("text", x = .62, y = 20, label = "Intervention (Treatment)", hjust = -.6) +
  geom_line(aes(color = group)) + geom_point(aes(color = group)) +
  scale_x_continuous(breaks = 1:10) +
  labs(x = "Time Period", color = "Experimental Group")
```

As we can see, the lines mirror more or less each other in the pre-treatment period, which seems to support that the parallel trends assumption is likely to be valid. Beyond this visual test, we also want to conduct an empirical test. For this, the did-package provides all relevant functions we need. The main functions is a pre-processing function and a function that estimates the ATT. The pre-processing function, called pre_process_did, requires as inputs the outcome (y), the variable name of the time period (t), the variable name of the unit (id), and the time of treatment which is equal to 8 if the unit is part of the treatment group, and 0 otherwise.

```{r}
# A) PRE-PROCESS
temp <- df %>%
  mutate(treat_t = case_when(g == 1 ~ 8, g == 0 ~ 0)) %>%
  pre_process_did(yname = "y", tname = "t", idname = "id", gname = "treat_t", allow_unbalanced_panel = T, data = .)

# B) ESTIMATION OF ATT
reg <- temp[['data']] %>%
  att_gt(yname = "y", tname = "t", idname = "id", gname = "treat_t", 
         control_group = c("nevertreated"), est_method = 'dr', 
         allow_unbalanced_panel = T, data = .)
reg
```

As the results show, we do not find any significant coefficients prior treament which provides (some) evidence that the parallel trends assumption is likely to be true. We can also plot these results with a lead-lag plot.

```{r}
ggdid(reg)
```

As expected the (visual and empirical) placebo tests all produce null results, indicating no differences in the trends in the pre-intervention periods, a result that would be helpful for convincing ourselves that the parallel trends assumption is reasonable.

The did-package has many different input parameters that can be used to specify the estimation approach in a specific way. This is particularly useful when we deal with staggered-DiD when we want to incorporate the advances of the new DiD literature on heterogenous treatment effects. We encourage you play around with the different parameters from the did-package and try to understand what they mean and what they do in the background.


## 4. Conclusion

This exercise was for illustrative purposes to show that mathematics and minimum data requirements necessary to obtain a DiD estimate is relatively simple. However, whether DiD provides a good ATT estimate depends solely on the validity of the parallel trends assumption. Unfortunately, we will never know whether PTA is valid because it is a counterfactual assumption. However, conducting pre-intervention placebo tests is the best (indirect) way to convince ourself whether PTA holds. 
