########################
### Replication code ###
########################

# Article: Between common responsibility and national interest:When do Europeans support a common European migration policy?

# Author: Philipp Lutz (University of Geneva & Vrije Universiteit Amsterdam)

### load libraries
library(foreign)
library(dplyr)
library(plyr)
library(ggplot2)
library(ggeffects)
library(ggpubr)
library(countrycode)
library(marginaleffects)
library(sf)
library(jtools)
library(grid)
library(scales)
library(RColorBrewer)
library(ggthemes)
library(eurostat) 
library(lmtest)
library(sandwich)
library(Hmisc)
library(sjPlot)
library(stargazer)
library(car)
library(viridis)

### load data
data_individual <- read.dta("data_individual.dta")
data_aggregate <- read.dta("data_aggregate.dta")
data_aggregate1 <- subset(data_aggregate, data_aggregate$year<2012)
data_aggregate2 <- subset(data_aggregate, data_aggregate$year>2012)
thresholds <- read.dta("thresholds.dta")

# Figure 1: Evolution of public support for an EU migration policy
data_individual <- subset(data_individual, !is.na(data_individual$cemp))
wtmeans <- data_individual %>% group_by(year) %>% mutate(cemp_mean = weighted.mean(cemp, weight))
wtsd <- data_individual %>% group_by(year) %>% mutate(cemp_sd = wtd.var(cemp, weight))
cemp_means <- aggregate(wtmeans$cemp_mean,by=list(wtmeans$year),FUN=mean, na.rm=TRUE, na.action="na.pass")
names(cemp_means)[1] <- "year"
names(cemp_means)[2] <- "cemp_mean"
cemp_sd <- aggregate(wtsd$cemp_sd,by=list(wtsd$year),FUN=mean)
names(cemp_sd)[1] <- "year"
names(cemp_sd)[2] <- "cemp_sd"
cemp_average <- merge(cemp_means,cemp_sd, by=c("year"),all=T)
cemp_average1 <- subset(cemp_average, cemp_average$year<2012)
cemp_average2 <- subset(cemp_average, cemp_average$year>2012)
jpeg("fig1.jpeg", width = 8, height = 5, units = "in", res=300)
qplot(cemp_average1$year,cemp_average1$cemp_mean) + geom_line() + xlab("Year") +
  geom_rect(aes(xmin = 1991.5, xmax = 1994.5, ymin = -Inf, ymax = Inf),fill = "gray", alpha = 0.01) +
  geom_rect(aes(xmin = 2003.5, xmax = 2006.5, ymin = -Inf, ymax = Inf),fill = "gray", alpha = 0.01) + 
  geom_rect(aes(xmin = 2012.5, xmax = 2019.5, ymin = -Inf, ymax = Inf),fill = "gray", alpha = 0.01) +
  annotate("text",x=1990, y=0.22, label=expression(underline("N countries")), size=3) +
  annotate("text",x=1993, y=0.22, label="12", size=3) +
  annotate("text",x=1999, y=0.22, label="15", size=3) +
  annotate("text",x=2005, y=0.22, label="25", size=3) +
  annotate("text",x=2009, y=0.22, label="27", size=3) +
  annotate("text",x=2016, y=0.22, label="28", size=3) +
  annotate("text",x=2020, y=0.22, label="27", size=3) +
  geom_vline(xintercept=c(1997.5,2000.5,2007.5), linetype="dashed", color = "#440154FF", size=0.25) + 
  annotate("text",x=1990, y=0.955, label=expression(underline("EU Treaty")), size=3) +
  annotate("text",x=1997.5, y=0.955, label="Amsterdam", size=3) +
  annotate("text",x=2000.5, y=0.955, label="Nice", size=3) +
  annotate("text",x=2007.5, y=0.955, label="Lisboa", size=3) + geom_hline(yintercept=c(0.25,0.5,0.75), linetype="dashed", color = "#440154FF", size=0.15) +
  scale_x_continuous("Year", c(1995,2000,2005,2010,2015,2020), labels=c("1995","2000","2005","2010","2015","2020")) +
  scale_y_continuous("Support EU migration policy (in %)", limits=c(0.22,0.955), c(0.25,0.5,0.75,1), labels=c("25%","50%","75%","100%")) +
  geom_errorbar(aes(x=cemp_average1$year, ymin=cemp_average1$cemp_mean-cemp_average1$cemp_sd, ymax=cemp_average1$cemp_mean+cemp_average1$cemp_sd), width=0.25, col="lightgray") + geom_errorbar(aes(x=cemp_average2$year, ymin=cemp_average2$cemp_mean-cemp_average2$cemp_sd, ymax=cemp_average2$cemp_mean+cemp_average2$cemp_sd), inherit.aes = F, width=0.25, col="lightgray") + 
  geom_point(data=cemp_average2, aes(x=year, y=cemp_mean)) +  geom_point(data=cemp_average1, aes(x=year, y=cemp_mean)) + geom_line(data=cemp_average2, aes(x=year, y=cemp_mean))  + theme_classic()
dev.off()
rm(wt_means,wtsd,cemp_means,cemp_means1,cemp_sd,cemp_sd1,cemp_average, cemp_average1, cemp_average2)


# Figure 2: Geographical distribution of support for an EU migration policy
means_year_country <- aggregate(data_individual$cemp, by=list(data_individual$year, data_individual$country),FUN=mean,na.rm=T)
names(means_year_country)[1] <- "year"
names(means_year_country)[2] <- "country"
names(means_year_country)[3] <- "value"
means_year_country$geo <- countrycode(means_year_country$country, origin = 'country.name', destination = 'iso2c')
means_year_country$geo[means_year_country$geo=="GB"] <- "UK"
means_year_country$geo[means_year_country$geo=="GR"] <- "EL"
means_year_country$value <- as.numeric(means_year_country$value)

brbg3 <- brewer.pal(11,"BrBG")[c(8,2,11)]
get_eurostat_geospatial(resolution = 10, nuts_level = 0, year = 2013)
SHP_0 <- get_eurostat_geospatial(resolution = 10, nuts_level = 0, year = 2013)
EU28 <- eu_countries %>% 
  select(geo = code, name)
SHP_28 <- SHP_0 %>% 
  select(geo = NUTS_ID, geometry) %>% 
  inner_join(EU28, by = "geo") %>% 
  arrange(geo) %>% 
  st_as_sf()
EU27 <- eu_countries %>% 
  filter(code != 'UK') %>% 
  select(geo = code, name)
SHP_27 <- SHP_0 %>% 
  select(geo = NUTS_ID, geometry) %>% 
  inner_join(EU27, by = "geo") %>% 
  arrange(geo) %>% 
  st_as_sf()

eb_shp <- means_year_country %>% 
  filter(year == 1995) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map1 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 1995",
    caption = ""
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))
eb_shp <- means_year_country %>% 
  filter(year == 2000) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map2 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 2000",
    caption = ""
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))
eb_shp <- means_year_country %>% 
  filter(year == 2005) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map3 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 2005",
    caption = ""
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))
eb_shp <- means_year_country %>% 
  filter(year == 2010) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map4 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 2010",
    caption = ""
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))
eb_shp <- means_year_country %>% 
  filter(year == 2015) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map5 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 2015",
    caption = ""
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))
eb_shp <- means_year_country %>% 
  filter(year == 2020) %>% 
  select(geo, value) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()
map6 <- eb_shp %>% 
  ggplot(aes(fill = value)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(midpoint=0.5, low="white", mid="gray80",
                       high="black", space ="Lab",
                       name = "Support in %",
                       breaks = pretty_breaks(10),
                       guide = guide_colorbar(
                         direction = "vertical", 
                         title.position = "top", 
                         label.position = "right",  
                         barwidth = unit(0.4, "cm"), 
                         barheight = unit(7, "cm"),  
                         ticks = TRUE,)) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(
    title = "",
    subtitle = "Year 2020",
    caption = "Data: Eurobarometer"
  ) +
  theme_void() +
  theme(legend.position = c(0.97, 0.50))

jpeg("fig2.jpeg", width = 11, height = 8, units = "in", res=300)
ggarrange(map1,map2,map3,map4,map5,map6, common.legend = TRUE, legend="right", hjust = -3, vjust = -3)
dev.off()

# Table 1: Panel regression estimates
model_base1 <- lm(cemp ~ exposure_absolute + exposure_relative + country, data=data_aggregate1) ##
coeftest(model_base1, vcov = vcovHC(model_base1, type="HC1"), cluster=~country) 
se_base1 <- sqrt(diag(vcovHC(model_base1, type = "HC1", cluster=~country)))
model_base2 <- lm(cemp ~ exposure_absolute + exposure_relative + country, data=data_aggregate2) ##
coeftest(model_base2, vcov = vcovHC(model_base2, type="HC1"), cluster=~country) 
se_base2 <- sqrt(diag(vcovHC(model_base2, type = "HC1", cluster=~country)))
model_base3 <- lm(cemp ~ exposure_absolute + exposure_relative + unemp + share_foreign + country, data=data_aggregate1) ## 
coeftest(model_base3, vcov = vcovHC(model_base3, type="HC1"), cluster=~country) 
se_base3 <- sqrt(diag(vcovHC(model_base3, type = "HC1", cluster=~country)))
model_base4 <- lm(cemp ~ exposure_absolute + exposure_relative + attach_eu + unemp + share_foreign + country, data=data_aggregate2) ## 
coeftest(model_base4, vcov = vcovHC(model_base4, type="HC1"), cluster=~country) 
se_base4 <- sqrt(diag(vcovHC(model_base4, type = "HC1", cluster=~country)))
stargazer(model_base1, model_base2,model_base3, model_base4,
          title="Panel regression estimates", 
          align=F, 
          dep.var.labels=c("Support EU migration policy"), 
          covariate.labels=c("Refugee arrivals (absolute)", "Refugee arrivals (relative)", "EU-attachment", "Unemployment rate", "Share of foreigners"),  
          no.space=TRUE, 
          dep.var.caption="", 
          se = list(se_base1,se_base2,se_base3,se_base4),
          model.numbers=T,
          star.cutoffs = c(0.05, 0.01, 0.001))


# Table 2: Difference-in-differences estimates
data_did <- subset(data_aggregate2, data_aggregate2$country!="Slovenia" & data_aggregate2$country!="Greece" & data_aggregate2$country!="Denmark" & data_aggregate2$country!="Croatia" & data_aggregate2$country!="Bulgaria") # 
data_did <- subset(data_did, data_did$year<2020)
data_did$`Crisis status` <- data_did$treat_lab
did_base <- lm(cemp ~ postcrisis*treatment + country*wave, data=data_did) #
coeftest(did_base, vcov = vcovHC(did_base, type="HC1"), cluster=~country) 
se_did <- sqrt(diag(vcovHC(did_base, type = "HC1", cluster=~country)))
stargazer(did_base,
          title="DiD estimates", 
          align=F, 
          dep.var.labels=c("Support EU migration policy"), 
          no.space=TRUE, 
          dep.var.caption="", 
          se = list(se_did),
          model.numbers=T,
          star.cutoffs = c(0.05, 0.01, 0.001))

# Figure 3: Heterogeneous effects of the refugee crisis by political orientation
data_did1 <- subset(data_individual, data_individual$country!="Slovenia" & data_individual$country!="Greece" & data_individual$country!="Denmark" & data_individual$country!="Croatia" & data_individual$country!="Bulgaria") #
data_did1 <- subset(data_did1,data_did1$year>2012 & data_did1$year<2020)
did_lrs <- lm(cemp ~ postcrisis*treat_lab*lrs, data=data_did1, weights=weight)
coeftest(did_lrs, vcov = vcovHC(did_lrs, type="HC1"), cluster=~country) 
se_did_lrs <- sqrt(diag(vcovHC(did_lrs, type = "HC1", cluster=~country)))
coefs <- coef(did_lrs)
pre_treat <- coefs[1] + coefs[3]
pre_untreat <- coefs[1] 
post_treat <- coefs[1] + coefs[3] + coefs[2]
post_untreat <- coefs[1] + coefs[2]
set_theme(base = theme_light())
p1 <- plot_model(did_lrs, type="pred", terms=c("postcrisis","treat_lab", "lrs [3.1,7.5]"),title="", grid.breaks=2, legend.title = "Crisis status", dot.size=2, line.size=1.1, color=c("#FDE725FF", "#440154FF"), fill=c("#FDE725FF", "#440154FF")) 
jpeg("fig3.jpeg", width = 10, height = 6, units = "in", res=300)
p1$data$facet <- factor(p1$data$facet, labels = c("Left-wing individuals (lrs = 3.1)", "Right-wing individuals (lrs = 7.5)"))
p1 + scale_y_continuous(name="Support EU migration policy", limits=c(0.64,0.88), breaks=c(0.65,0.7,0.75,0.8,0.85),labels=c("65 %", "70 %", "75 %", "80 %", "85 %")) +  scale_x_continuous(name="", limits=c(-0.1,1.1), breaks=c(0.1,0.9), labels=c("before\n crisis", "after\n crisis")) + theme_get() + 
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_line(linetype = 2), panel.grid.minor.y = element_blank() ) 
dev.off()




### Online Appendix

# Table A1: Descriptive statistics of the variables (1992-2012)
table_summary1 <- data.frame(data_aggregate1$cemp, data_aggregate1$exposure_absolute, data_aggregate1$exposure_relative, data_aggregate1$attach_eu, data_aggregate1$unemp, data_aggregate1$share_foreign)
stargazer(table_summary1, covariate.labels=c("Support EU migration policy", "Refugee arrivals (absolute)", "Refugee arrivals (relative)", "EU attachment", "Unemployment rate", "Share of foreigners"))
# Table A2: Descriptive statistics of variables (2014-2021)
table_summary2 <- data.frame(data_aggregate2$cemp, data_aggregate2$exposure_absolute, data_aggregate2$exposure_relative, data_aggregate2$attach_eu, data_aggregate2$unemp, data_aggregate2$share_foreign)
stargazer(table_summary2, covariate.labels=c("Support EU migration policy", "Refugee arrivals (absolute)", "Refugee arrivals (relative)", "EU attachment", "Unemployment rate", "Share of foreigners"))

# Figure A1: Histograms of the refugee exposure variables
jpeg("fig_A1.jpeg",   width = 12, height = 6, units = "in", res=300)
par (mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")   # 
par(mfrow = c(2,2))
hist(data_aggregate1$exposure_absolute, main="Refugee arrivals (absolute)" , xlab="Asylum requests per 100'000 inh. (log)", xlim=c(0,1), ylab="", border=F)
hist(data_aggregate2$exposure_absolute, main="Refugee arrivals (absolute)" , xlab="Asylum requests per 100'000 inh. (log)", xlim=c(0,1), ylab="", border=F)
hist(data_aggregate1$exposure_relative, main="Refugee arrivals (relative)", xlab="Share of asylum requests (relative to EU average), log", xlim=c(0,1), ylab="", border=F)
hist(data_aggregate2$exposure_relative, main="Refugee arrivals (relative)", xlab="Share of asylum requests (relative to EU average), log", xlim=c(0,1), ylab="", border=F)
dev.off()

# Figure A2: Scatter plot of country clusters
thresholds$cluster <- "non-exposed"
thresholds$cluster[thresholds$country=="Germany" | thresholds$country=="Finland" | thresholds$country=="Hungary" | thresholds$country=="Sweden" | thresholds$country=="Austria"] <- "exposed"
jpeg("fig_A2.jpeg",   width = 8, height = 6, units = "in", res=300)
ggplot(thresholds, aes(x = threshold1, y = threshold, color = cluster)) +
  geom_point(size = 2) + ylab("Arrivals as % of population (2015/2016)") + xlab("Arrivals as % of population (2015)") +
  theme_minimal()
dev.off()

# Figure A3: Component-plus-residual plots of panel regression models
jpeg("fig_A3.jpeg",   width = 20, height = 12, units = "in", res=300)
par (mar=c(4,4,0,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")   # 
par(mfrow=c(2,2))
crPlots(model_base1,  terms=~exposure_absolute, col.lines=c("#440154FF", "#FDE725FF"), grid.lines=2, main="2000-2012", ylab="Component+Residual (Policy support)", xlab="Refugee arrivals (absolute)")
text(0.5, 75, "2000-2012", cex=1.5, font=2)
crPlots(model_base1,  terms=~exposure_relative, col.lines=c("#440154FF", "#FDE725FF"), main="2000-2012", ylab="Component+Residual (Policy support)", xlab="Refugee arrivals (relative)")
text(0.5, 44, "2000-2012", cex=1.5, font=2)
crPlots(model_base2,  terms=~exposure_absolute, col.lines=c("#440154FF", "#FDE725FF"), main="2014-2021", ylab="Component+Residual (Policy support)", xlab="Refugee arrivals (absolute)")
text(0.6, 25, "2014-2021", cex=1.5, font=2)
crPlots(model_base2,  terms=~exposure_relative, col.lines=c("#440154FF", "#FDE725FF"), main="2014-2021", ylab="Component+Residual (Policy support)", xlab="Refugee arrivals (relative)")
text(0.5, 20, "2014-2021", cex=1.5, font=2)
dev.off()

# Figure A4: Predicted probabilities of policy support by number of refugee arrivals
pred_base1 <- plot_predictions(model_base1, condition=list("exposure_absolute", "exposure_relative"="threenum"), rug=T) + theme_minimal() + ggtitle("Period 2000-2012") + xlab("Countries' exposure to refugees (absolute)") + ylab("Support EU migration policy") + scale_color_viridis(discrete = TRUE)
pred_base2 <- plot_predictions(model_base2, condition=list("exposure_absolute", "exposure_relative"="threenum"),rug=T) + theme_minimal() + ggtitle("Period 2014-2021") + xlab("Countries' exposure to refugees (absolute)") + ylab("Support EU migration policy") + scale_color_viridis(discrete = TRUE)
jpeg("fig_A4.jpeg",   width = 15, height = 7, units = "in", res=300)
ggarrange(pred_base1, pred_base2, ncol=2, nrow=1)
dev.off()

# Table A4: Alternative panel model specifications
model_alt1 <- lm(cemp ~ exposure_absolute + country, data=data_aggregate1)
coeftest(model_alt1, vcov = vcovHC(model_alt1, type="HC1", cluster=~country)) 
se_alt1 <- sqrt(diag(vcovHC(model_alt1, type = "HC1", cluster=~country)))
model_alt2 <- lm(cemp ~ exposure_absolute + country, data=data_aggregate2) 
coeftest(model_alt2, vcov = vcovHC(model_alt2, type="HC1"), cluster=~country) 
se_alt2 <- sqrt(diag(vcovHC(model_alt2, type = "HC1", cluster=~country)))
data_aggregate1_lowcor <- subset(data_aggregate1,  data_aggregate1$country!="Slovenia" & data_aggregate1$country!="Slovakia" & data_aggregate1$country!="Malta" &  data_aggregate1$country!="Luxembourg" &  data_aggregate1$country!="Latvia" & data_aggregate1$country!="Ireland" & data_aggregate1$country!="Greece" & data_aggregate1$country!="Germany" & data_aggregate1$country!="Czech Republic" &  data_aggregate1$country!="Cyprus" & data_aggregate1$country!="Belgium")
data_aggregate2_lowcor <- subset(data_aggregate2,  data_aggregate2$country!="Slovenia" & data_aggregate2$country!="Spain" & data_aggregate2$country!="Sweden" &  data_aggregate2$country!="Luxembourg" &  data_aggregate2$country!="Italy" & data_aggregate2$country!="Ireland" & data_aggregate2$country!="Hungary" & data_aggregate2$country!="Denmark" & data_aggregate2$country!="Cyprus" &  data_aggregate2$country!="Croatia" & data_aggregate2$country!="Austria")
model_alt3 <- lm(cemp ~ exposure_absolute + exposure_relative + country, data=data_aggregate1_lowcor) #
coeftest(model_alt3, vcov = vcovHC(model_alt3, type="HC1"), cluster=~country) 
se_alt3 <- sqrt(diag(vcovHC(model_alt3, type = "HC1", cluster=~country)))
model_alt4 <- lm(cemp ~ exposure_absolute + exposure_relative + country, data=data_aggregate2_lowcor) #
coeftest(model_alt4, vcov = vcovHC(model_alt4, type="HC1"), cluster=~country) 
se_alt4 <- sqrt(diag(vcovHC(model_alt4, type = "HC1", cluster=~country)))
stargazer(model_alt1, model_alt2, model_alt3, model_alt4,
          title="Alternative panel model specifications", 
          align=F, 
          dep.var.labels=c("Support EU migration policy"), 
          no.space=TRUE, 
          dep.var.caption="", 
          se = list(se_alt1,se_alt2,se_alt3,se_alt4),
          model.numbers=T,
          star.cutoffs = c(0.05, 0.01, 0.001))

# Figure A5: Support for an EU migration policy by country group
data_individual <- subset(data_individual, !is.na(data_individual$cemp))
wtmeans1 <- data_individual %>% group_by(year,cgroup) %>% mutate(cemp_mean = weighted.mean(cemp, weight), na.rm = TRUE)
wtsd1 <- data_individual %>% group_by(year,cgroup) %>% mutate(cemp_sd = wtd.var(cemp, weight), na.rm = TRUE)
cemp_means1 <- aggregate(wtmeans1$cemp_mean,by=list(wtmeans1$year,wtmeans1$cgroup),FUN=mean, na.rm=T)
names(cemp_means1)[1] <- "year"
names(cemp_means1)[2] <- "group"
names(cemp_means1)[3] <- "cemp_mean"
cemp_sd1 <- aggregate(wtsd1$cemp_sd,by=list(wtsd1$year,wtsd1$cgroup),FUN=mean)
names(cemp_sd1)[1] <- "year"
names(cemp_sd1)[2] <- "group"
names(cemp_sd1)[3] <- "cemp_sd"
cemp_average_group <- merge(cemp_means1,cemp_sd1, by=c("year", "group"),all=T)
cemp_average_north1 <- subset(cemp_average_group, cemp_average_group$group=="North" & cemp_average_group$year<=2012)
cemp_average_south1 <- subset(cemp_average_group, cemp_average_group$group=="South" & cemp_average_group$year<=2012)
cemp_average_east1 <- subset(cemp_average_group, cemp_average_group$group=="East" & cemp_average_group$year<=2012)
cemp_average_east2 <- subset(cemp_average_group, cemp_average_group$group=="East" & cemp_average_group$year>2012)
cemp_average_north2 <- subset(cemp_average_group, cemp_average_group$group=="North" & cemp_average_group$year>2012)
cemp_average_south2 <- subset(cemp_average_group, cemp_average_group$group=="South" & cemp_average_group$year>2012)
jpeg("fig_A5.jpeg", width = 8, height = 5, units = "in", res=300)
qplot(cemp_average_north1$year,cemp_average_north1$cemp_mean) + geom_point(shape=17, color="#440154FF") + geom_line(color="#440154FF") + xlab("Year")  + 
  geom_point(data=cemp_average_south1, shape=15, aes(x=year, y=cemp_mean),color="#FDE725FF") +  geom_line(data=cemp_average_south1, aes(x=year, y=cemp_mean),color="#FDE725FF") +
  geom_point(data=cemp_average_east1, shape=16, aes(x=year, y=cemp_mean),color="gray30") +  geom_line(data=cemp_average_east1, aes(x=year, y=cemp_mean),color="gray30") +
  geom_point(data=cemp_average_south2, shape=15, aes(x=year, y=cemp_mean),color="#FDE725FF") +  geom_line(data=cemp_average_south2, aes(x=year, y=cemp_mean),color="#FDE725FF") +
  geom_point(data=cemp_average_east2, shape=16, aes(x=year, y=cemp_mean),color="gray30") +  geom_line(data=cemp_average_east2, aes(x=year, y=cemp_mean),color="gray30") +
  geom_point(data=cemp_average_north2, shape=17, aes(x=year, y=cemp_mean),color="#440154FF") +  geom_line(data=cemp_average_north2, aes(x=year, y=cemp_mean),color="#440154FF") +
  geom_vline(xintercept=c(1997.5,2000.5,2007.5), linetype="dashed", color = "gray", size=0.25) + 
  annotate("text",x=1992, y=0.95, label=expression(underline("EU Treaty")), size=3) +
  annotate("text",x=1997.5, y=0.95, label="Amsterdam", size=3) +
  annotate("text",x=2000.5, y=0.95, label="Nice", size=3) + 
  annotate("text",x=2007.5, y=0.95, label="Lisboa", size=3) +
  annotate("text",x=2021, y=0.87, label="South", size=3, color="#FDE725FF") +
  annotate("text",x=2021, y=0.805, label="North", size=3, color="#440154FF") +
  annotate("text",x=2021, y=0.64, label="East", size=3, color="gray50") +
  scale_x_continuous("Year", c(1995,2000,2005,2010,2015,2020), labels=c("1995","2000","2005","2010","2015","2020")) +
  scale_y_continuous("Support EU migration policy (in %)", c(0.25,0.5,0.75,1), labels=c("25%","50%","75%","100%")) + theme_classic()
dev.off()
rm(wtmeans1,wtsd1,cemp_average_group, cemp_average_east1, cemp_average_east2, cemp_average_north1, cemp_average_north2, cemp_average_south1,cemp_average_south2)

# Figure A6: Pre-treatment trend of support for an EU migration policy
means_cemp <- aggregate(data_aggregate1$cemp, by=list(data_aggregate1$wave, data_aggregate1$treatment), FUN=mean, na.rm=T) # calculate wave-means
jpeg("fig_A6.jpeg",   width = 8, height = 5, units = "in", res=300)
par (mar=c(3,5,1,0.5), mgp=c(0,.7,0), tck=-.01, las=1)
par(mfrow=c(1,1))
plot(means_cemp$Group.1[means_cemp$Group.2==1], means_cemp$x[means_cemp$Group.2==1], main="", type='l', lwd=2, col="#FDE725FF", xlab="", ylim=c(0,100),xaxt='n',yaxt='n', ylab="")
lines(means_cemp$Group.1[means_cemp$Group.2==0], means_cemp$x[means_cemp$Group.2==0], col="#440154FF", lwd=2)
points(means_cemp$Group.1[means_cemp$Group.2==0], means_cemp$x[means_cemp$Group.2==0], col="#440154FF", pch=16)
points(means_cemp$Group.1[means_cemp$Group.2==1], means_cemp$x[means_cemp$Group.2==1], col="#FDE725FF", pch=16)
axis(side=1, at=c(1,6,16,26,33), labels=c("1992", "1995", "2000", "2005", "2014"), cex=0.7)
axis(side=2, at=c(0, 20,40,60,80,100), labels=c("0 %", "20 %", "40 %", "60 %", "80 %", "100 %"))
legend(25,20, legend=c("exposed"), cex=0.85, box.lwd=0, bty='n', col="#FDE725FF",lwd=2)
legend(25,25, legend=c("not exposed"), cex=0.85, box.lwd=0, bty='n', col="#440154FF",lwd=2)
mtext("Year", side = 1, line=2,adj=0.5, cex = 1, las=1)
mtext("Public support EU migration policy", side = 2, line=3.5,adj=0.5, cex = 1, las=3)
dev.off()

# Figure A7: Country-specific trends in support for an EU migration policy
jpeg("fig_A7.jpeg",   width = 12, height = 10, units = "in", res=300)
grob <- grobTree(textGrob("Crisis peak", x=0.14,  y=0.95, hjust=0,
                          gp=gpar(col="red", fontsize=13, fontface="italic")))
data_did %>%
  ggplot(aes(x=wave, y=cemp, group=country, color=`Crisis status`)) +  geom_line(lwd=1) + xlab("Survey wave") + ylab("Support EU migration policy") +
  scale_y_continuous("Support EU migration policy (in %)", limits=c(0,100), c(20,40,60,80,100), labels=c("20%", "40%","60%","80%","100%")) + 
  scale_x_continuous(name="Year", breaks=c(34,36,38,40,42,44), labels=c("2014","2015", "2016","2017", "2018", "2019")) +
  scale_color_manual(values=c("#FDE725FF", "#440154FF")) + geom_vline(xintercept = c(35,36), linetype="dashed",  color = "gray", size=1) + annotation_custom(grob) + 
  theme_bw() + geom_text(data = data_did[data_did$wave == 44,], aes(label = country),size = 3, hjust = 0.5, vjust = 1, color="black")                                                                                                                                             
dev.off()

# Figure A8: Predicted probabilities of support for EU migration policy
dataset_did <- subset(data_individual, data_individual$year>2012 & data_individual$year<=2019)
did2a <- lm(cemp ~ wave2*treat_lab, data=dataset_did, weights = weight) ## 
jpeg("fig_A8.jpeg", width = 10, height = 6, units = "in", res=300)
plot_model(did2a, type="pred", terms=c("wave2","treat_lab"),title="", grid.breaks=2, legend.title = "Crisis status", dot.size=2, line.size=1.1, color=c("#FDE725FF", "#440154FF")) + geom_vline(xintercept=35.5, linetype="dashed") + scale_y_continuous(name="Support EU migration policy", limits=c(0.68,0.85), breaks=c(0.7,0.75,0.8),labels=c("70 %", "75 %", "80 %")) + scale_x_continuous(name="Year", breaks=c(34,36,38,40,42,44), labels=c("2014","2015", "2016","2017", "2018", "2019")) +  theme_blank()
dev.off()

# Table A5: DiD estimates with alternative model specifications
west <- subset(data_did, data_did$country=="Cyprus" | data_did$country=="Greece" | data_did$country=="Italy" | data_did$country=="Malta" | data_did$country=="Portugal" | data_did$country=="Spain" | data_did$country=="Austria" | data_did$country=="Belgium" | data_did$country=="Denmark" | data_did$country=="Finland" | data_did$country=="France" | data_did$country=="Germany" | data_did$country=="Ireland" | data_did$country=="Luxembourg" | data_did$country=="Netherlands" | data_did$country=="Sweden" | data_did$country=="United Kingdom")
did_west <- lm(cemp ~ postcrisis*treatment + country*wave, data=west) ##
coeftest(did_west, vcov = vcovHC(did_west, type="HC1"), cluster=~country) 
se_did_west <- sqrt(diag(vcovHC(did_west, type = "HC1", cluster=~country)))
# use continuous treatment variable of gradual exposure: variable 'threshold'
data_did <- merge(data_did, thresholds, by=c("country"), all=T)
did_gradual1 <- lm(cemp ~ postcrisis*threshold + country*wave, data=data_did) # + as.factor(wave)
coeftest(did_gradual1, vcov = vcovHC(did_gradual1, type="HC1"), cluster=~country) 
cov <- vcovHC(did_gradual1, type = "HC1", cluster=~country)
se_did_gradual1 <- sqrt(diag(cov))
did_gradual2 <- lm(cemp ~ postcrisis*threshold1 + country*wave, data=data_did) # + as.factor(wave)
coeftest(did_gradual2, vcov = vcovHC(did_gradual2, type="HC1"), cluster=~country) 
cov <- vcovHC(did_gradual2, type = "HC1", cluster=~country)
se_did_gradual2 <- sqrt(diag(cov))
# estimate model with transit countries
did_transit <- lm(cemp ~ postcrisis*treatment + country*wave, data=data_aggregate2)
coeftest(did_transit, vcov = vcovHC(did_transit, type="HC1"), cluster=~country) 
cov <- vcovHC(did_transit, type = "HC1", cluster=~country)
se_did_transit <- sqrt(diag(cov))
# add control variables to the model
did_control <- lm(cemp ~ postcrisis*treatment + attach_eu + unemp + share_foreign + country*wave, data=data_did)
coeftest(did_control, vcov = vcovHC(did_control, type="HC1"), cluster=~country) 
cov <- vcovHC(did_control, type = "HC1", cluster=~country)
se_did_control <- sqrt(diag(cov))
# export table
stargazer(did_west,did_transit,did_gradual1,did_gradual2,did_control,
          title="DiD estimates with alternative model specifications", 
          align=F, 
          dep.var.labels=c("Support EU migration policy"), 
          no.space=TRUE, 
          dep.var.caption="", 
          se = list(se_did_west, se_did_transit,se_did_gradual1, se_did_gradual2, se_did_control),
          model.numbers=T,
          star.cutoffs = c(0.05, 0.01, 0.001))