Labour Event Study with Control Counties

Functions

## Functions
with.overlay <- function(est, s) { attr(est,'overlay') = s; est }

estimators <- function(sdid,sc,s) {
  estimator.list = list(with.overlay(sdid, s), sc, sdid)
  names(estimator.list)=c('SDID+fixed effect', 'synth. control', 'synth. diff-in-diff')
  estimator.list
}

plot.estimators <- function(ests, alpha.multiplier) {
  p = synthdid_plot(ests, se.method='none',
                    alpha.multiplier=alpha.multiplier, facet=rep(1,length(ests)),
                    trajectory.linetype = 1, effect.curvature=-.4,
                    trajectory.alpha=.5, effect.alpha=.5, diagram.alpha=1)
  suppressMessages(p + scale_alpha(range=c(0,1), guide='none'))
}

Civil Labor

df2 <- read.csv("Civil Labor Ballance.csv") %>%
    filter(COUNTY != "jefferson parish", COUNTY != "orleans parish",
        COUNTY != "plaquemines parish", COUNTY != "st. bernard parish",
        COUNTY != "st. charles parish", COUNTY != "st. helena parish",
        COUNTY != "st. james parish") %>%
    filter(STATE != "FL" & STATE != "AL" & STATE != "MS") %>%
    mutate(DATE = as.Date(DATE),
        treat = ifelse(STATE == "LA" & DATE >= as.Date("2005-08-01"), 1, 0),
        year = as.numeric(format(DATE, "%Y"))) %>%
    group_by(COUNTY, STATE) %>%
    filter(!any(is.na(Civil.Labor)))

Balancing Data and Calculating Unit Weights

dt <- seq.Date(from = as.Date("2000-01-01"),
    to = as.Date("2015-12-01"), by = "month")

df.fix <- data.frame(FIPS = rep(unique(df2$FIPS), length(dt))) %>%
    arrange(FIPS) %>%
    mutate(DATE = rep(dt, length(unique(df2$FIPS))))

bal_df <- left_join(df.fix, df2, by = c("FIPS", "DATE")) %>%
    select(COUNTY, FIPS, STATE, DATE, treat, Civil.Labor) %>%
    group_by(FIPS) %>%
    filter(!any(is.na(Civil.Labor))) %>%
    ungroup()

bal_df$FIPS <- factor(bal_df$FIPS)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- log(bal_df$Civil.Labor)

setup <- panel.matrices(bal_df,
    unit = "FIPS", time = "DATE", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)

unit.weights <- synthdid_controls(sdid.df, weight.type = 'omega', mass = 1)
unit.df <- data.frame(sdid = round(unit.weights[rev(1:nrow(unit.weights)), ], 3)) %>%
    rownames_to_column(var = "FIPS")

head(unit.df)

##    FIPS sdid
## 1  8091    0
## 2  6065    0
## 3 40033    0
## 4 40099    0
## 5 48259    0
## 6 46071    0

Subsetting data by Weight Counties (SDID > 0)

map_df <- unit.df %>%
    mutate(sdid_mp = ifelse(sdid > 0, 1, 0))

cnty_fips <- df2 %>%
    distinct(COUNTY, FIPS)

map_sdid <- map_df %>%
  filter(sdid_mp == 1) %>%
  merge(cnty_fips, by = "FIPS")

regdf <- df2 %>%
    filter(STATE == "LA" | FIPS %in% map_sdid$FIPS)

df3 <- regdf %>%
    filter(DATE < as.Date("2015-08-01")) %>%
    mutate(treat = as.numeric(STATE == "LA"))

Civil Labor Regression Model

## Makes Date Dummy
dloop <- data.frame(DATE = unique(df3$DATE))

dloop <- dloop %>%
  mutate(year = format(DATE,"%Y"),
         month=str_pad(string = format(DATE,"%m"),width = 2,side = "left",pad = 0))
head(dloop)

##         DATE year month
## 1 2000-01-01 2000    01
## 2 2000-02-01 2000    02
## 3 2000-03-01 2000    03
## 4 2000-04-01 2000    04
## 5 2000-05-01 2000    05
## 6 2000-06-01 2000    06

## YEAR DUMMY
for (i in 1:nrow(dloop)) {
  df3[paste0("d.", dloop$year[i],".",dloop$month[i])] <- as.numeric(df3$DATE == unique(df3$DATE)[i])
}

## RUN MODEL
did.reg <- felm(log(Civil.Labor) ~
           treat:(d.2000.01+d.2000.02+d.2000.03+d.2000.04+d.2000.05+d.2000.06+d.2000.07+d.2000.08+d.2000.09+d.2000.10+d.2000.11+d.2000.12+
           d.2001.01+d.2001.02+d.2001.03+d.2001.04+d.2001.05+d.2001.06+d.2001.07+d.2001.08+d.2001.09+d.2001.10+d.2001.11+d.2001.12+
           d.2002.01+d.2002.02+d.2002.03+d.2002.04+d.2002.05+d.2002.06+d.2002.07+d.2002.08+d.2002.09+d.2002.10+d.2002.11+d.2002.12+
           d.2003.01+d.2003.02+d.2003.03+d.2003.04+d.2003.05+d.2003.06+d.2003.07+d.2003.08+d.2003.09+d.2003.10+d.2003.11+d.2003.12+
           d.2004.01+d.2004.02+d.2004.03+d.2004.04+d.2004.05+d.2004.06+d.2004.07+d.2004.08+d.2004.09+d.2004.10+d.2004.11+d.2004.12+
           d.2005.01+d.2005.02+d.2005.03+d.2005.04+d.2005.05+d.2005.06+d.2005.07+d.2005.09+d.2005.10+d.2005.11+d.2005.12+
           d.2006.01+d.2006.02+d.2006.03+d.2006.04+d.2006.05+d.2006.06+d.2006.07+d.2006.08+d.2006.09+d.2006.10+d.2006.11+d.2006.12+
           d.2007.01+d.2007.02+d.2007.03+d.2007.04+d.2007.05+d.2007.06+d.2007.07+d.2007.08+d.2007.09+d.2007.10+d.2007.11+d.2007.12+
           d.2008.01+d.2008.02+d.2008.03+d.2008.04+d.2008.05+d.2008.06+d.2008.07+d.2008.08+d.2008.09+d.2008.10+d.2008.11+d.2008.12+
           d.2009.01+d.2009.02+d.2009.03+d.2009.04+d.2009.05+d.2009.06+d.2009.07+d.2009.08+d.2009.09+d.2009.10+d.2009.11+d.2009.12+
           d.2010.01+d.2010.02+d.2010.03+d.2010.04+d.2010.05+d.2010.06+d.2010.07+d.2010.08+d.2010.09+d.2010.10+d.2010.11+d.2010.12+
           d.2011.01+d.2011.02+d.2011.03+d.2011.04+d.2011.05+d.2011.06+d.2011.07+d.2011.08+d.2011.09+d.2011.10+d.2011.11+d.2011.12+
           d.2012.01+d.2012.02+d.2012.03+d.2012.04+d.2012.05+d.2012.06+d.2012.07+d.2012.08+d.2012.09+d.2012.10+d.2012.11+d.2012.12+
           d.2013.01+d.2013.02+d.2013.03+d.2013.04+d.2013.05+d.2013.06+d.2013.07+d.2013.08+d.2013.09+d.2013.10+d.2013.11+d.2013.12+
           d.2014.01+d.2014.02+d.2014.03+d.2014.04+d.2014.05+d.2014.06+d.2014.07+d.2014.08+d.2014.09+d.2014.10+d.2014.11+d.2014.12+
           d.2015.01+d.2015.02+d.2015.03+d.2015.04+d.2015.05+d.2015.06+d.2015.07) |
           COUNTY + DATE | 0 | COUNTY,
              data = df3)

## Export Coeffiecent and test for Heteroscedasticity SE
coef_df <- coeftest(did.reg)

coef_df <- coef_df[,] %>%
    as_tibble() %>%
    mutate(variables = rownames(coef_df))

coef_df$date <- as.Date(paste0(coef_df$variables, ".01"),
    format = "treat:d.%Y.%m.%d")
 
event_df <- data.frame(date=seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by = "month"))
event_df <- left_join(event_df, coef_df[, c(1, 2, 6)])

## Joining with `by = join_by(date)`

event_df[is.na(event_df)] <- 0

colnames(event_df) <- c("time","coef","se")

## Makes Confident Interval
event_df$ci_upper <- event_df$coef + 1.96 * event_df$se
event_df$ci_lower <- event_df$coef - 1.96 * event_df$se
head(event_df,3)

##         time        coef          se     ci_upper    ci_lower
## 1 2000-01-01 -0.02199253 0.009273899 -0.003815685 -0.04016937
## 2 2000-02-01 -0.03010525 0.008785003 -0.012886642 -0.04732386
## 3 2000-03-01 -0.02473713 0.008821614 -0.007446771 -0.04202750

dates <- seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by='month')
# Define start and end dates for plot
start_date <- as.Date("2004-08-01")
end_date <- as.Date("2015-07-01")

# Create sequence of dates for x-axis
dates_seq <- seq.Date(from = start_date, to = end_date, by = "month")

# Subset event_df to only include dates after or equal to start_date
event_df_sub <- event_df[dates >= start_date, ]

pl <- ggplot(event_df, aes(x = dates, y = coef)) +
  geom_point() +
  geom_line()+
  #geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = as.Date("2005-08-01"), linetype = "dashed", color = "blue") +
  scale_x_date(date_labels = "%b %Y", limits = c(start_date, end_date), breaks = "1 year") +
  labs(title = "Civil Labor") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14, face = "bold"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text = element_text(size = 12),
    legend.text = element_text(size = 12),
    plot.caption = element_text(size = 12, face = c("italic", "bold")),
    axis.line = element_line(color = "black", size = 0.5),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_line(color = "black", size = 0.1))
pl

Employee Person

df4 <- read.csv("Merge Employee Person Full-Location Ballance.csv") %>%
    filter(COUNTY != "Jefferson Parish", COUNTY != "Orleans parish",
        COUNTY != "Plaquemines Parish", COUNTY != "St. Bernard parish",
        COUNTY != "St. Charles Parish", COUNTY != "St. Helena Parish",
        COUNTY != "St. James Parish") %>%
    filter(STATE != "FL" & STATE != "AL" & STATE != "MS") %>%
    mutate(DATE = as.Date(DATE),
        treat = ifelse(STATE == "LA" & DATE >= as.Date("2005-08-01"), 1, 0),
        year = as.numeric(format(DATE, "%Y"))) %>%
    group_by(COUNTY, STATE) %>%
    filter(!any(is.na(Employee)))

Balancing Data and Calculating Unit Weights

dt <- seq.Date(from = as.Date("2000-01-01"),
    to = as.Date("2015-12-01"), by = "month")

df.fix <- data.frame(FIPS = rep(unique(df4$FIPS), length(dt))) %>%
    arrange(FIPS) %>%
    mutate(DATE = rep(dt, length(unique(df4$FIPS))))

bal_df <- left_join(df.fix, df4, by = c("FIPS", "DATE")) %>%
    select(COUNTY, FIPS, STATE, DATE, treat, Employee) %>%
    group_by(FIPS) %>%
    filter(!any(is.na(Employee))) %>%
    ungroup()

bal_df$FIPS <- factor(bal_df$FIPS)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- log(bal_df$Employee)

setup <- panel.matrices(bal_df,
    unit = "FIPS", time = "DATE", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)

unit.weights <- synthdid_controls(sdid.df, weight.type = 'omega', mass = 1)
unit.df <- data.frame(sdid = round(unit.weights[rev(1:nrow(unit.weights)), ], 3)) %>%
    rownames_to_column(var = "FIPS")

head(unit.df)

##    FIPS sdid
## 1 10005    0
## 2 13043    0
## 3 13101    0
## 4 13117    0
## 5 13151    0
## 6 13177    0

Subsetting data by Weight Counties (SDID > 0)

map_df <- unit.df %>%
    mutate(sdid_mp = ifelse(sdid > 0, 1, 0))

cnty_fips <- df4 %>%
    distinct(COUNTY, FIPS)

map_sdid <- map_df %>%
  filter(sdid_mp == 1) %>%
  merge(cnty_fips, by = "FIPS")

regdf <- df4 %>%
    filter(STATE == "LA" | FIPS %in% map_sdid$FIPS)

df5 <- regdf %>%
    filter(DATE < as.Date("2015-08-01")) %>%
    mutate(treat = as.numeric(STATE == "LA"))

Employee Person Regression Model

## Makes Date Dummy
dloop <- data.frame(DATE = unique(df5$DATE))

dloop <- dloop %>%
  mutate(year = format(DATE,"%Y"),
         month=str_pad(string = format(DATE,"%m"),width = 2,side = "left",pad = 0))
head(dloop)

##         DATE year month
## 1 2000-01-01 2000    01
## 2 2000-02-01 2000    02
## 3 2000-03-01 2000    03
## 4 2000-04-01 2000    04
## 5 2000-05-01 2000    05
## 6 2000-06-01 2000    06

## YEAR DUMMY
for (i in 1:nrow(dloop)) {
  df5[paste0("d.", dloop$year[i],".",dloop$month[i])] <- as.numeric(df5$DATE == unique(df5$DATE)[i])
}

## RUN MODEL
did.reg <- felm(log(Employee) ~
           treat:(d.2000.01+d.2000.02+d.2000.03+d.2000.04+d.2000.05+d.2000.06+d.2000.07+d.2000.08+d.2000.09+d.2000.10+d.2000.11+d.2000.12+
           d.2001.01+d.2001.02+d.2001.03+d.2001.04+d.2001.05+d.2001.06+d.2001.07+d.2001.08+d.2001.09+d.2001.10+d.2001.11+d.2001.12+
           d.2002.01+d.2002.02+d.2002.03+d.2002.04+d.2002.05+d.2002.06+d.2002.07+d.2002.08+d.2002.09+d.2002.10+d.2002.11+d.2002.12+
           d.2003.01+d.2003.02+d.2003.03+d.2003.04+d.2003.05+d.2003.06+d.2003.07+d.2003.08+d.2003.09+d.2003.10+d.2003.11+d.2003.12+
           d.2004.01+d.2004.02+d.2004.03+d.2004.04+d.2004.05+d.2004.06+d.2004.07+d.2004.08+d.2004.09+d.2004.10+d.2004.11+d.2004.12+
           d.2005.01+d.2005.02+d.2005.03+d.2005.04+d.2005.05+d.2005.06+d.2005.07+d.2005.09+d.2005.10+d.2005.11+d.2005.12+
           d.2006.01+d.2006.02+d.2006.03+d.2006.04+d.2006.05+d.2006.06+d.2006.07+d.2006.08+d.2006.09+d.2006.10+d.2006.11+d.2006.12+
           d.2007.01+d.2007.02+d.2007.03+d.2007.04+d.2007.05+d.2007.06+d.2007.07+d.2007.08+d.2007.09+d.2007.10+d.2007.11+d.2007.12+
           d.2008.01+d.2008.02+d.2008.03+d.2008.04+d.2008.05+d.2008.06+d.2008.07+d.2008.08+d.2008.09+d.2008.10+d.2008.11+d.2008.12+
           d.2009.01+d.2009.02+d.2009.03+d.2009.04+d.2009.05+d.2009.06+d.2009.07+d.2009.08+d.2009.09+d.2009.10+d.2009.11+d.2009.12+
           d.2010.01+d.2010.02+d.2010.03+d.2010.04+d.2010.05+d.2010.06+d.2010.07+d.2010.08+d.2010.09+d.2010.10+d.2010.11+d.2010.12+
           d.2011.01+d.2011.02+d.2011.03+d.2011.04+d.2011.05+d.2011.06+d.2011.07+d.2011.08+d.2011.09+d.2011.10+d.2011.11+d.2011.12+
           d.2012.01+d.2012.02+d.2012.03+d.2012.04+d.2012.05+d.2012.06+d.2012.07+d.2012.08+d.2012.09+d.2012.10+d.2012.11+d.2012.12+
           d.2013.01+d.2013.02+d.2013.03+d.2013.04+d.2013.05+d.2013.06+d.2013.07+d.2013.08+d.2013.09+d.2013.10+d.2013.11+d.2013.12+
           d.2014.01+d.2014.02+d.2014.03+d.2014.04+d.2014.05+d.2014.06+d.2014.07+d.2014.08+d.2014.09+d.2014.10+d.2014.11+d.2014.12+
           d.2015.01+d.2015.02+d.2015.03+d.2015.04+d.2015.05+d.2015.06+d.2015.07) |
           COUNTY + DATE | 0 | COUNTY,
              data = df5)

## Export Coeffiecent and test for Heteroscedasticity SE
coef_df <- coeftest(did.reg)

coef_df <- coef_df[,] %>%
    as_tibble() %>%
    mutate(variables = rownames(coef_df))

coef_df$date <- as.Date(paste0(coef_df$variables, ".01"),
    format = "treat:d.%Y.%m.%d")
 
event_df <- data.frame(date=seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by = "month"))
event_df <- left_join(event_df, coef_df[, c(1, 2, 6)])

## Joining with `by = join_by(date)`

event_df[is.na(event_df)] <- 0

colnames(event_df) <- c("time","coef","se")

## Makes Confident Interval
event_df$ci_upper <- event_df$coef + 1.96 * event_df$se
event_df$ci_lower <- event_df$coef - 1.96 * event_df$se
head(event_df,3)

##         time        coef          se    ci_upper    ci_lower
## 1 2000-01-01 -0.02606679 0.009502501 -0.00744189 -0.04469169
## 2 2000-02-01 -0.03020137 0.009043699 -0.01247572 -0.04792702
## 3 2000-03-01 -0.02978123 0.009207189 -0.01173514 -0.04782732

dates <- seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by='month')
# Define start and end dates for plot
start_date <- as.Date("2004-08-01")
end_date <- as.Date("2015-07-01")

# Create sequence of dates for x-axis
dates_seq <- seq.Date(from = start_date, to = end_date, by = "month")

# Subset event_df to only include dates after or equal to start_date
event_df_sub <- event_df[dates >= start_date, ]

pl1 <- ggplot(event_df, aes(x = dates, y = coef)) +
  geom_point() +
  geom_line()+
  #geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = as.Date("2005-08-01"), linetype = "dashed", color = "blue") +
  scale_x_date(date_labels = "%b %Y", limits = c(start_date, end_date), breaks = "1 year") +
  labs(title = "Employee Person") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14, face = "bold"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text = element_text(size = 12),
    legend.text = element_text(size = 12),
    plot.caption = element_text(size = 12, face = c("italic", "bold")),
    axis.line = element_line(color = "black", size = 0.5),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_line(color = "black", size = 0.1))
pl1

Unemployment Rate

df6 <- read.csv("Merge UnemRate ballance V2 FULL-Location.csv") %>%
    filter(COUNTY != "Jefferson Parish", COUNTY != "Orleans parish",
        COUNTY != "Plaquemines Parish", COUNTY != "St. Bernard parish",
        COUNTY != "St. Charles Parish", COUNTY != "St. Helena Parish",
        COUNTY != "St. James Parish") %>%
    filter(STATE != "FL" & STATE != "AL" & STATE != "MS") %>%
    mutate(DATE = as.Date(DATE),
        treat = ifelse(STATE == "LA" & DATE >= as.Date("2005-08-01"), 1, 0),
        year = as.numeric(format(DATE, "%Y"))) %>%
    group_by(COUNTY, STATE) %>%
    filter(!any(is.na(Unemployee.Rate)))

Balancing Data and Calculating Unit Weights

dt <- seq.Date(from = as.Date("2000-01-01"),
    to = as.Date("2015-12-01"), by = "month")

df.fix <- data.frame(FIPS = rep(unique(df6$FIPS), length(dt))) %>%
    arrange(FIPS) %>%
    mutate(DATE = rep(dt, length(unique(df6$FIPS))))

bal_df <- left_join(df.fix, df6, by = c("FIPS", "DATE")) %>%
    select(COUNTY, FIPS, STATE, DATE, treat, Unemployee.Rate) %>%
    group_by(FIPS) %>%
    filter(!any(is.na(Unemployee.Rate))) %>%
    ungroup()

bal_df$FIPS <- factor(bal_df$FIPS)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- log(bal_df$Unemployee.Rate)

setup <- panel.matrices(bal_df,
    unit = "FIPS", time = "DATE", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)

unit.weights <- synthdid_controls(sdid.df, weight.type = 'omega', mass = 1)
unit.df <- data.frame(sdid = round(unit.weights[rev(1:nrow(unit.weights)), ], 3)) %>%
    rownames_to_column(var = "FIPS")

head(unit.df)

##    FIPS sdid
## 1 10005    0
## 2 13013    0
## 3 13015    0
## 4 13023    0
## 5 13045    0
## 6 13047    0

Subsetting data by Weight Counties (SDID > 0)

map_df <- unit.df %>%
    mutate(sdid_mp = ifelse(sdid > 0, 1, 0))

cnty_fips <- df6 %>%
    distinct(COUNTY, FIPS)

map_sdid <- map_df %>%
  filter(sdid_mp == 1) %>%
  merge(cnty_fips, by = "FIPS")

regdf <- df6 %>%
    filter(STATE == "LA" | FIPS %in% map_sdid$FIPS)

df7 <- regdf %>%
    filter(DATE < as.Date("2015-08-01")) %>%
    mutate(treat = as.numeric(STATE == "LA"))

Unemployment Rate Regression Model

## Makes Date Dummy
dloop <- data.frame(DATE = unique(df7$DATE))

dloop <- dloop %>%
  mutate(year = format(DATE, "%Y"),
         month=str_pad(string = format(DATE,"%m"),width = 2,side = "left",pad = 0))

## YEAR DUMMY
for (i in 1:nrow(dloop)) {
  df7[paste0("d.", dloop$year[i],".",dloop$month[i])] <- as.numeric(df7$DATE == unique(df7$DATE)[i])
}

## RUN MODEL
did.reg <- felm(log(Unemployee.Rate) ~
           treat:(d.2000.01+d.2000.02+d.2000.03+d.2000.04+d.2000.05+d.2000.06+d.2000.07+d.2000.08+d.2000.09+d.2000.10+d.2000.11+d.2000.12+
           d.2001.01+d.2001.02+d.2001.03+d.2001.04+d.2001.05+d.2001.06+d.2001.07+d.2001.08+d.2001.09+d.2001.10+d.2001.11+d.2001.12+
           d.2002.01+d.2002.02+d.2002.03+d.2002.04+d.2002.05+d.2002.06+d.2002.07+d.2002.08+d.2002.09+d.2002.10+d.2002.11+d.2002.12+
           d.2003.01+d.2003.02+d.2003.03+d.2003.04+d.2003.05+d.2003.06+d.2003.07+d.2003.08+d.2003.09+d.2003.10+d.2003.11+d.2003.12+
           d.2004.01+d.2004.02+d.2004.03+d.2004.04+d.2004.05+d.2004.06+d.2004.07+d.2004.08+d.2004.09+d.2004.10+d.2004.11+d.2004.12+
           d.2005.01+d.2005.02+d.2005.03+d.2005.04+d.2005.05+d.2005.06+d.2005.07+d.2005.09+d.2005.10+d.2005.11+d.2005.12+
           d.2006.01+d.2006.02+d.2006.03+d.2006.04+d.2006.05+d.2006.06+d.2006.07+d.2006.08+d.2006.09+d.2006.10+d.2006.11+d.2006.12+
           d.2007.01+d.2007.02+d.2007.03+d.2007.04+d.2007.05+d.2007.06+d.2007.07+d.2007.08+d.2007.09+d.2007.10+d.2007.11+d.2007.12+
           d.2008.01+d.2008.02+d.2008.03+d.2008.04+d.2008.05+d.2008.06+d.2008.07+d.2008.08+d.2008.09+d.2008.10+d.2008.11+d.2008.12+
           d.2009.01+d.2009.02+d.2009.03+d.2009.04+d.2009.05+d.2009.06+d.2009.07+d.2009.08+d.2009.09+d.2009.10+d.2009.11+d.2009.12+
           d.2010.01+d.2010.02+d.2010.03+d.2010.04+d.2010.05+d.2010.06+d.2010.07+d.2010.08+d.2010.09+d.2010.10+d.2010.11+d.2010.12+
           d.2011.01+d.2011.02+d.2011.03+d.2011.04+d.2011.05+d.2011.06+d.2011.07+d.2011.08+d.2011.09+d.2011.10+d.2011.11+d.2011.12+
           d.2012.01+d.2012.02+d.2012.03+d.2012.04+d.2012.05+d.2012.06+d.2012.07+d.2012.08+d.2012.09+d.2012.10+d.2012.11+d.2012.12+
           d.2013.01+d.2013.02+d.2013.03+d.2013.04+d.2013.05+d.2013.06+d.2013.07+d.2013.08+d.2013.09+d.2013.10+d.2013.11+d.2013.12+
           d.2014.01+d.2014.02+d.2014.03+d.2014.04+d.2014.05+d.2014.06+d.2014.07+d.2014.08+d.2014.09+d.2014.10+d.2014.11+d.2014.12+
           d.2015.01+d.2015.02+d.2015.03+d.2015.04+d.2015.05+d.2015.06+d.2015.07) |
           COUNTY + DATE | 0 | COUNTY,
              data = df7)

## Export Coeffiecent and test for Heteroscedasticity SE
coef_df <- coeftest(did.reg)

coef_df <- coef_df[,] %>%
    as_tibble() %>%
    mutate(variables = rownames(coef_df))

coef_df$date <- as.Date(paste0(coef_df$variables, ".01"),
    format = "treat:d.%Y.%m.%d")
 
event_df <- data.frame(date=seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by = "month"))
event_df <- left_join(event_df, coef_df[, c(1, 2, 6)])

## Joining with `by = join_by(date)`

event_df[is.na(event_df)] <- 0

colnames(event_df) <- c("time","coef","se")

## Makes Confident Interval
event_df$ci_upper <- event_df$coef + 1.96 * event_df$se
event_df$ci_lower <- event_df$coef - 1.96 * event_df$se
head(event_df,3)

##         time        coef         se     ci_upper    ci_lower
## 1 2000-01-01 -0.01626689 0.02093174  0.024759325 -0.05729311
## 2 2000-02-01 -0.12875379 0.02240263 -0.084844632 -0.17266295
## 3 2000-03-01 -0.04782277 0.02159325 -0.005500005 -0.09014553

dates <- seq.Date(as.Date("2000-01-01"),as.Date("2015-07-01"),by='month')
# Define start and end dates for plot
start_date <- as.Date("2004-08-01")
end_date <- as.Date("2015-07-01")

# Create sequence of dates for x-axis
dates_seq <- seq.Date(from = start_date, to = end_date, by = "month")

# Subset event_df to only include dates after or equal to start_date
event_df_sub <- event_df[dates >= start_date, ]

pl2 <- ggplot(event_df, aes(x = dates, y = coef)) +
  geom_point() +
  geom_line()+
  #geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = as.Date("2005-08-01"), linetype = "dashed", color = "blue") +
  scale_x_date(date_labels = "%b %Y", limits = c(start_date, end_date), breaks = "1 year") +
  labs(title = "Unemployment Rate") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(size = 14, face = "bold"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text = element_text(size = 12),
    legend.text = element_text(size = 12),
    plot.caption = element_text(size = 12, face = c("italic", "bold")),
    axis.line = element_line(color = "black", size = 0.5),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_line(color = "black", size = 0.1))
pl2

All plots

layout_matrix <- rbind(c(1, 1, 2, 2), c(4, 3, 3, 4))
grid.arrange(pl, pl1, pl2, layout_matrix = layout_matrix)