---
title: "Problem Set 7"
author: "H. Busshoff & P. Heinemann"
date: "17 Juni 2018"
output:
  word_document: default
  pdf_document: default
  html_document:
    df_print: paged
---

#1 Background

#1

Police presence and crime are likely to be simultaneously determined. Cities with high crime rates are likely to employ more police. Therefore, the correlation between police presence and crime cannot easily be attributed to a one-directional causal link.

# 2

Terrorist attack on July 18th, 1994 on the main Jewish center in Argentina. A week later, the government assigned police protection to every Jewish and Muslim building in the country. The authors used data on motor vehicle theft per block in different neighbourhoods before and after the attack and information about the location of Jewish institutions in these neighbourhoods.

#3

The different groups are different blocks, some of which received treatment (blocks where Jewish institutions are located and thus police protection has been assigned after the attack) and some of which didn't (where no such institution is located).
Expectations are that being treated or being close to a treated block should reduce car theft in a block and that the effect reduces with increasing distance to a treated block.


#4

Common trend assumption: In the absence of more police, the development of thefts in the different blocks would have been parallel.
Given this assumption, the difference between how strongly theft changed in treated blocks and how strongly theft changed in non-treated blocks after the attack can be attributed to the causal effect of police presence on car theft.

# 5

Find evidence that police lowers crime. Relative to the control group, car thefts fall by 75 percent in the blocks in which the protected institutions are situated. Effect is local.  Police presence in a given block does not reduce car theft one or two blocks away.


#2 Empirics

Importing Data

```{r}
df <- read.csv("C:/Users/Hannah/Desktop/ps7_crimerates.csv")
```


Loading packages

```{r}
library(tidyr)
library(dplyr)
library(psych)
library(skimr)
library(ggplot2)
library(mice)
library(miceadds)
```

Let's have a brief look at data

```{r}
skim(df)
```


#1

Recoding variables to distinguish relevant dimensions.

```{r}
df$D = 0
df$D[df$distance == 0] =  1
df$T = 0
df$T[df$week >= 16] = 1
```

#2

Computing all conditional means

```{r}
e = data.frame(matrix(NA, ncol = 2, nrow = 2))

for (i in 0:1) {
  for (j in 0:1)
    e[i+1, j+1] = mean(df$NumCarThefts[df$D == i & df$T == j])
  }
```

Deriving the causal estimates

```{r}
treatment = e[2,2] - e[2,1] - (e[1,2] - e[1,1])
```

The treatment induces the relative number of auto thefts to decrease by approx. 0.018 pts.

#3

Alternative to compute causal effect via a linear regression

```{r}
lm1 = lm.cluster(data = df, NumCarThefts ~ D + T + D:T, cluster="blockID")
summary(lm1)
```
Gives same result as it should be.

#4

```{r message=TRUE, warning=FALSE, paged.print=FALSE}
df$block = as.factor(df$blockID)
lm2 = lm.cluster(data = df, NumCarThefts ~ D + T + D:T + block, cluster="block")
coef(lm2)[c("(Intercept)","D","T","D:T")]
```

Effect of interaction term does not change.

#5

They will be already subsumed in the fixed effects as they do not change over time as the following shows.


```{r}

x = df %>%
  group_by(blockID, week) %>%
  summarise (
    bank = mean(df$dumBank),
    inst = mean(df$dumPublicInstitution)
    )
y = x %>%
  group_by(blockID) %>%
  summarise(
  varbank = var(bank),
  varinst = var(inst)
  )

skim(y)
```

# 6

```{r}
df$weekD = as.factor(df$week)
lm3 = lm.cluster(data = df, NumCarThefts ~ D + T + D:T + block + weekD, cluster="block")
coef(lm3)[c("(Intercept)","D","T","D:T")]
```

No, coefficient does not change.

# 7

Exploring the common trends assumption

```{r message=FALSE, warning=FALSE}
y = df %>%
  group_by(D, week) %>%
summarise(
  y = mean(NumCarThefts))

y$D = as.factor(y$D)

ggplot(data = y, aes(x = week, y = y, color = D)) +
  geom_line(aes(group = D))


```


# 8

Testing the common trend assumption via a regression

```{r}
reg4 = lm(df$NumCarThefts ~ D + T + D:T + week + week:D, data = df)
summary(reg4)
```

Interaction week and dummy variable not significant. Therefore, conclude no violation of the common trend assumption.

Alternatively, use an interaction of week indicators and dummy

```{r}
reg5 = lm(df$NumCarThefts ~ T + D + weekD + weekD:D, data = df)

anova(reg5)
```

No evidence that common trend assumption is violated.

# Part 3: Synthetic controls

# 1

Importing data

```{r}
df1 <- read.csv("C:/Users/Hannah/Downloads/ps7_birthrates.csv")

library(ggplot2)
library(zoo)
df1$date = with(df1, interaction(df1$year, df1$quarter))
df1$date1 = as.Date(as.yearqtr(df1$date, format = "%Y.%q"))

df1$date2 = as.numeric(df1$date1)
ggplot(data = df1, aes(x = date1, y = birthrate, color = countyname)) +
  geom_line(aes(group = countyname)) +
  theme(axis.text.x = element_text(angle = 90)) + geom_point() +
  geom_vline(xintercept = as.numeric(as.Date("1995-04-01")), linetype=4) +
  geom_vline(xintercept = as.numeric(as.Date("1996-01-01")), linetype = 4)

```

No, common trend assumtion seems to be violated.


# 2

```{r}
df1$D = 0
df1$D[df1$countyname == 'oklahoma'] = 1
df1$D = as.factor(df1$D)

x = df1 %>%
  group_by(D, date1) %>%
  summarise(
    births = mean(birthrate))

ggplot(data = x, aes(x = date1, y = births, color = D)) +
  geom_line(aes(group = D)) +
  theme(axis.text.x = element_text(angle = 90)) + geom_point() +
  geom_vline(xintercept = as.numeric(as.Date("1995-04-01")), linetype=4) +
  geom_vline(xintercept = as.numeric(as.Date("1996-01-01")), linetype = 4)
```

Yes, looks much better now.

#3

Load package for synthetic controls

```{r}
library(Synth)

df1$countyname = as.factor(df1$countyname)
df1$SD = as.numeric(df1$countyname)
df1$countyname = as.character(df1$countyname)

a = unique(df1$date2[df1$year == 1992 & df1$quarter == 1])
b = unique(df1$date2[df1$year == 1995 & df1$quarter == 4])
c = unique(df1$date2[df1$year == 1990 & df1$quarter == 1 ])

d = unique(df1$date2[df1$date2 <= b])
e = unique(df1$date2[df1$date2 >= a & df1$date2 <= b])
f = unique(df1$date2[df1$year == 1990 & df1$quarter == 1])

dataprepout1 = dataprep(foo = df1, predictors = "schoolenrollment", predictors.op = "mean", dependent = "birthrate", unit.variable = "SD", unit.names.variable = "countyname", time.variable = "date2", treatment.identifier = 7, controls.identifier = c(1:6, 8:12), time.optimize.ssr = c(d), time.predictors.prior = f)


synth.out1 <- synth(data.prep.obj = dataprepout1,
                    method = "BFGS")


synth.tables <- synth.tab(dataprep.res = dataprepout1,
                           synth.res = synth.out1
                           )

synth.tables
```


```{r}
x  = synth.tables$tab.w

df1$weight = 0

for (ii in c(1:6, 8:12)) {
  df1$weight[df1$SD == ii] = x$w.weights[x$unit.numbers == ii]
}

df1$y = df1$birthrate * df1$weight
df1$y[df1$D == 1] = df1$birthrate[df1$D == 1]

y = df1 %>%
  group_by(date1, D) %>%
  summarise(
    births = sum(y))

ggplot(data = y, aes(x = date1, y = births, color = D)) +
  geom_line(aes(group = D)) +
  geom_vline(xintercept = as.numeric(as.Date("1996-01-01")), linetype=4)


```

# 5


Assessing significance of the effect via randomization inference - we take all counties (without weighting scheme) as units of observations:

```{r}
df1$T = 0
df1$T[as.numeric(as.Date(df1$date1)) >= as.numeric(as.Date("1996-01-01"))] = 1

fisher = df1 %>%
group_by(countyname) %>%
summarise(
pretreatment = mean(birthrate[T==0]),
posttreatment = mean(birthrate[T==1])
)

fisher$delta = fisher$posttreatment - fisher$pretreatment

#Under H0: Delta(0) = Delta(1) for all counties

a = sum(fisher$delta)

#There are 12 possible permutations in the setting (N over NT).

j = 12

X = data.frame(matrix(NA, ncol = 3, nrow = j))
names(X) = c("Treatment Effect H0", "Absolute Value", "Relative Size")

ref = fisher$delta[7] - 1/11*sum(fisher$delta[c(1:6, 8:12)])

for (i in c(1:6, 8:12)) {
  b = fisher$delta[i]
  X[i, 1] = b - (a - b)*1/11
  X[i, 2] = abs(X[i, 1] )
  if (X[i,2] >= ref ) {
    X[i,3] = 1
  } else {
    X[i,3] = 0
  }
  }

p = (1 + sum(X$V3[c(1:6, 8:12)]))/12
```

Effect (weakly) significant - computed p-value is 0.083 percent.