---
title: "tutorial9_lars"
format: html
editor: visual
---

## Setup

```{r}
library(pacman)
p_load(ggplot2, data.table, fixest, haven)
```

## Code für die Simulation

```{r}
# Set seed for reproducibility
set.seed(123)

# Parameters
n_groups <- 100   # Number of groups
n_periods <- 10   # Number of time periods
n <- n_groups * n_periods

# Simulate groups and time periods
data <- data.table(
  group = rep(1:n_groups, each = n_periods),
  time = rep(1:n_periods, times = n_groups)
)

# Simulate fixed effects
data[, group_fe := rnorm(1 , sd = 100), by = group]
data[, time_fe := rnorm(1, sd = 100), by = time]

# Simulate treatment variable and dependent variable
data[, treatment := sample(0:1, n, replace = TRUE)]
data[, epsilon := rnorm(n)]
data[, y := 5 + 2 * treatment + group_fe + time_fe + epsilon]
```

Wir simulieren unsere Daten mit dem bereitgestellten Code.

## Schätzen eines klassischen linearen Regressionsmodells

```{r}
lm_model <- lm(y ~ treatment, data = data)
summary(lm_model)
```

## Und jetzt, zum Vergleich, mit Fixed Effects

```{r}
# Kleinste Quadrate Schätzer mit Fixed Effects
fe_model <- feols(y ~ treatment | group + time, data = data)
summary(fe_model)
```

Wenn wir für die Effekte von `group`und `time` kontrollieren, sehen wir, dass der Effekt von `treatment` höher eingeschätzt wird, und der Standardfehler deutlich kleiner ist, wir uns also deutlich sicherer in unserer Schätzung sind. Zudem wird dies als hochsignifikant eingestuft.

## Zum Vergleich mit einer anderen Standardabweichung

```{r}
# Set seed for reproducibility
set.seed(123)

# Simulate groups and time periods
new_data <- data.table(
  group = rep(1:n_groups, each = n_periods),
  time = rep(1:n_periods, times = n_groups)
)

# Simulate fixed effects
new_data[, group_fe := rnorm(1, sd = 5), by = group]
new_data[, time_fe := rnorm(1, sd = 200), by = time]

# Simulate treatment variable and dependent variable
new_data[, treatment := sample(0:1, n, replace = TRUE)]
new_data[, epsilon := rnorm(n)]
new_data[, y := 5 + 2 * treatment + group_fe + time_fe + epsilon]

# Fixed Effects Model
new_fe_model <- feols(y ~ treatment | group + time, data = new_data)
summary(new_fe_model)
```

Der Code simuliert uns einen neuen Datensatz, bei dem die Standardabweichung des `time`-Effekts doppelt so hoch ist, während die Standardabweichung des `group`-Effekts deutlich geringer ist. Diese Veränderung ist für unsere Schätzung mit Hilfe von Fixed-Effects Modellen nicht relevant, da wir die Effekte von `group` und `time` kontrollieren. Der Effekt von `treatment` bleibt also unverändert. Dies kann gerne auch mit anderen Werten ausprobiert werden, das Ergebnis bleibt gleich. Wichtig ist, dass für die Simulation der Daten, sofern sie neu generiert werden, der gleiche Seed verwendet wird, andernfalls erhalten wir andere Werte für unseren `Treatment`-Effekt.

# Zum Simpson's Paradox

Das hier angesprochene Paradoxon lässt sich sehr gut visuell beobachten, insbesondere auf der Wikipedia-Seite zum \`Simpson's paradox' gibt es diese schöne Animation: https://en.wikipedia.org/wiki/Simpson%27s_paradox#/media/File:Simpsons_paradox\_-\_animation.gif Sie verdeutlicht, dass Fixed Effects den wahren Effekt, teils erheblich, verschleiern können, und die wahre kausale Richtung zum Teil entgegen der naiven Schätzung wirken kann.

### Difference-In-Differences

```{r}
# Data.Table anlegen
did_data <- data.table(
    Gruppe = c("Treatment", "Treatment", "Control", "Control"),
    Zeit = c("Vorher", "Nacher", "Vorher", "Nacher"),
    Wert = c(70, 85, 75, 80)
    )

# Kopie der Daten erstellen, als Sicherungskopie (und Replizierbarkeit später)
did_data_backup <- copy(did_data)

# Daten visualisieren
did_plot <- ggplot(did_data, aes(x = Zeit, y = Wert, fill = Gruppe)) +
    geom_bar(stat = "identity", position = "dodge") +
    labs(title = "Differences-In-Differences Visualisierung", x = "Zeitpunkt", y = "Wert") +
    theme_minimal()
did_plot
```

Wir sehen, dass die Werte für beide Gruppen nach dem Treatment ansteigen, aber der Anstieg bei der Treatment-Gruppe deutlich stärker ist. Dies ist der Effekt, den wir mit einem Differences-In-Differences Schätzer messen wollen.

## Difference-In-Differences Schätzer

```{r}
# Wir nutzen unsere Sicherungskopie der Daten, sodass wir diesen Code-Chunk immer wieder ausführen können, ohne dass wir unsere Daten dauerhaft verändern und Replizierbarkeit verloren geht.
did_data <- copy(did_data_backup)

# Zeitwerte in numerische Werte übersetzen
did_data[, Zeit := ifelse(Zeit == "Vorher", 0, 1)]

# Gruppenwerte in numerische Werte übersetzen
did_data[, Gruppe := ifelse(Gruppe == "Treatment", 1, 0)]

# Berechnung des Effekts mit unserer "normalen" Formel über die Differenzen
diff1 <- did_data[Gruppe == 1 & Zeit == 1]$Wert - did_data[Gruppe == 1 & Zeit == 0]$Wert
diff2 <- did_data[Gruppe == 0 & Zeit == 1]$Wert - did_data[Gruppe == 0 & Zeit == 0]$Wert

diff_effect <- diff1 - diff2
diff_effect
```

Diff-In-Diff Effekt beträgt 10.


## Jetzt als Lineare Regression

```{r}
# Erstellen unseres Interaktionsterms
did_data[, Interaktion := Gruppe * Zeit]
# Wird nur zu 1, wenn wir in Treatment-Gruppe und nach Verabreichung des Treatments

# Lineares Regressionsmodell schätzen
fit_model <- lm(Wert ~ Gruppe + Zeit + Interaktion, data = did_data)
summary(fit_model)
```

Wir sehen, dass unsere Diff-in-Diff Schätzung mit Hilfe eines linearen Regressionsmodells den gleichen Effekt schätzt, wie unsere manuelle Berechnung. Dieser ist in unserem Interaktionsterm wiederzufinden, da dieser genau den Fall des isolierten Treatment-Effekts darstellt. Zudem schätzen wir auch zeitliche Effekte in unserer `post`-Variable und Gruppeneffekte in unserer `group`Variable, welch hier aber keine Bedeutung besitzen.


### Agrawal und Goldfarb - Ein Beispiel aus der Forschung

Kernfoschungsfrage: Wie beeinflusst die Einführung von `bothbitnet` die Kollaboration in der Zukunft (hier gemessen an der Anzahl an Co-Autoren)?

Was bedeutet die Variablen:
i,j = Universitäten
t = Zeitpunkt
$\mu_t$ = Fixed Effects der Zeitpunkte
$\phi_ij$ = Fixed Effects der Universitäten

```{r}
# Laden der Daten
goldfarb <- haven::read_dta("../input/bitnet_paired_final.dta")
goldfarb <- as.data.table(goldfarb)

# Identifizieren der relevanten Variablen
colnames(goldfarb)

# DiD-Effekt berechnen
# lm(dcoauths1_2 ~ bothbitnet + as.factor(year) + as.factor(gpinst), data = goldfarb)
# führt (bei mir) zur Fehlermeldung, da zu viele Paare existieren und zu viel Speicher benötigt wird

# Jetzt mit Fixed Effects
feols(dcoauths1_2 ~ bothbitnet | gpinst + year , data = goldfarb)
```

In diesem Beispiel wird ein Untersuchung von Agrawal und Goldfarb (2008) nachvollzogen. Hierbei wird der Effekt von `bothbitnet` auf `dcoauths1_2` untersucht. Es wird ein einfaches lineares Regressionsmodell geschätzt, sowie ein Fixed Effects Modell. Die Berechnung mit klassischem `lm`ist allerdings so groß, dass der Code nur als Beispiel dient, und hier nicht ausgeführt wird. Mit Fixed Effects funktioniert die Schätzung allerdings problemlos und ist äußerst schnell. Der errechnete Wert, 0.000852, ist hochsignifikant und zeigt, dass `bothbitnet` einen positiven Effekt auf `dcoauths1_2` hat. Im Paper wird der Effekt exakt mit dem gleichen Wert angegeben.
