# Star Wars at Central Banks Script - Z-curve charts
# Last Updated: 27 January 2021

library(tidyverse)
library(haven)
library(ggplot2)
library(vctrs)
library(ragg)
source("Programs/E_paper-formatted_graphs/starwars_graph_theme.R") # ggplot theme for graphs
source("Programs/E_paper-formatted_graphs/adjust_kernal_boundaries.R") # extra code for fixing the boundaries of the kernel estimates

starwars_data <- read_dta("Data/Final/final_stars_supp.dta")

# Figure 1
Figure1_data <- starwars_data %>%
  mutate(dataset = case_when(D_cb == 1 ~ "central_bank",
                             D_topJ == 1 ~ "top_journals",
                             TRUE ~ NA_character_)) %>%
  select(dataset, journal_id, eye_catcher, exploratory, data_driven, peer, 
         data_availability, codes_availability, t_stat_sm) %>%
  filter(t_stat_sm < 10.0 & t_stat_sm >= 0 & !is.na(dataset))

facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0

Figure1 <- ggplot(data = Figure1_data, aes(x = t_stat_sm)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.04999471, center = 0.024997355, fill = "#F8766D") +
  stat_density(size = 0.8, kernel = 'epanechnikov', bw = 0.2, geom = "line", position = "identity", bounds = c(0, 11)) +
  facet_grid(cols = vars(dataset)) +
  geom_text(aes(label = groupings, x = 5, y = c(0.3625, 0.3775)),
            data = data.frame(dataset = unique(Figure1_data[["dataset"]]), 
                              groupings = c("Central bank discussion papers\n(2000-2019)", 
                                            "Top journals (2005-2011)")),
            size = 5.5) +
  theme_starwars() +
  ylab("Density") +
  xlab("|z-score|") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.4), sec.axis = sec_axis(~., name = "Density")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 10),  breaks = z10_breaks, labels = z10_labels) +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust = 1, margin=unit(c(0,5,0,0), "mm")),
        axis.title.y.right = element_text(size = rel(1), angle = 90, vjust = 1, margin=unit(c(0,0,0,5), "mm")))

# Figure 5
facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0

data_driven_data <- starwars_data %>%
  filter(D_dataDriven_cb == 1) %>%
  mutate(analysis_type = "data_driven")

reverse_causal <- starwars_data %>%
  filter(D_explore_cb == 1) %>%
  mutate(analysis_type = "reverse_causal")

neither <- starwars_data %>%
  filter(D_noExpData_cb == 1) %>%
  mutate(analysis_type = "neither")
  
Figure5_data <- rbind(data_driven_data, reverse_causal, neither) %>%
  select(analysis_type, journal_id, eye_catcher, exploratory, data_driven, peer, 
         data_availability, codes_availability, t_stat_sm) %>%
  filter(t_stat_sm < 10.0 & t_stat_sm >= 0) %>%
  mutate(analysis_type = factor(analysis_type, levels = c("data_driven",
                                                             "reverse_causal",
                                                             "neither")))

Figure5 <- ggplot(data = Figure5_data, aes(x = t_stat_sm)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, center = 0.025, fill = "#F8766D") +
  stat_density(size = 0.8, kernel = 'epanechnikov', bw = 0.2, geom = "line", position = "identity", bounds = c(0, 11)) +
  facet_grid(cols = vars(analysis_type)) +
  geom_text(aes(label = groupings, x = 5, y = c(0.565, 0.565, 0.52)),
            data = data.frame(analysis_type = unique(Figure5_data[["analysis_type"]]), 
                              groupings = c("Data-driven model selection", 
                                            "Reverse causal research",
                                            "Neither data-driven model\nselection or\nreverse causal research")),
            size = 5.5) +
  theme_starwars() +
  ylab("Density") +
  xlab("|z-score|") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.6), breaks = c(0, 0.2, 0.4, 0.6), 
                     sec.axis = sec_axis(~., name = "Density", breaks = c(0, 0.2, 0.4, 0.6))) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 10),  breaks = z10_breaks, labels = z10_labels) +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust = 1, margin=unit(c(0,5,0,0), "mm")),
        axis.title.y.right = element_text(size = rel(1), angle = 90, vjust = 1, margin=unit(c(0,0,0,5), "mm")))

#Figure 6
Figure6_data <- starwars_data %>%
  filter(D_control == 1 & t_stat_sm < 10.0 & t_stat_sm >= 0) %>%
  select(journal_id, eye_catcher, exploratory, data_driven, peer, 
         data_availability, codes_availability, t_stat_sm)

Figure6 <- ggplot(data = Figure6_data, aes(x = t_stat_sm)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, center = 0.025, fill = "#F8766D") +
  stat_density(size = 0.8, kernel = 'epanechnikov', bw = 0.2, geom = "line", position = "identity", bounds = c(0, 11)) +
  theme_starwars() +
  ylab("Density") +
  xlab("|z-score|") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.4), breaks = c(0, 0.1, 0.2, 0.3, 0.4), 
                     sec.axis = sec_axis(~., name = "Density", breaks = c(0, 0.1, 0.2, 0.3, 0.4))) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 10),  breaks = c(0, 1.96, 10), labels = c("0", "1.96", "10")) +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust = 1, margin=unit(c(0,5,0,0), "mm")),
        axis.title.y.right = element_text(size = rel(1), angle = 90, vjust = 1, margin=unit(c(0,0,0,5), "mm")))

#Figure 8
facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0

Figure8_data <- starwars_data %>%
  mutate(central_bank = case_when(D_rba == 1 ~ "RBA",
                             D_rbnz == 1 ~ "RBNZ",
                             D_minn == 1 ~ "minn_fed",
                             TRUE ~ NA_character_)) %>%
  select(central_bank, D_noExpData_cb, journal_id, eye_catcher, exploratory, data_driven, peer, 
         data_availability, codes_availability, t_stat_sm) %>%
  filter(t_stat_sm < 10.0 & t_stat_sm >= 0 & D_noExpData_cb == 1 & !is.na(central_bank))

Figure8 <- ggplot(data = Figure8_data, aes(x = t_stat_sm)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, center = 0.025, fill = "#F8766D") +
  stat_density(size = 0.8, kernel = 'epanechnikov', bw = 0.2, geom = "line", position = "identity", bounds = c(0, 11)) +
  facet_grid(cols = vars(central_bank)) +
  geom_text(aes(label = groupings, x = 5, y = c(0.57, 0.57, 0.57)),
            data = data.frame(central_bank = unique(Figure8_data[["central_bank"]]), 
                              groupings = c("Minneapolis Fed", 
                                            "Reserve Bank of Australia",
                                            "Reserve Bank of New Zealand")),
            size = 5.5) +
  theme_starwars() +
  ylab("Density") +
  xlab("|z-score|") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.6), breaks = c(0, 0.2, 0.4, 0.6), 
                     sec.axis = sec_axis(~., name = "Density", breaks = c(0, 0.2, 0.4, 0.6))) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 10),  breaks = z10_breaks,labels = z10_labels) +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust = 1, margin=unit(c(0,5,0,0), "mm")),
        axis.title.y.right = element_text(size = rel(1), angle = 90, vjust = 1, margin=unit(c(0,0,0,5), "mm")))

#Figure 9
Figure9_data <- starwars_data %>%
  filter(D_noExpData_cb == 1 & t_stat_sm < 10.0 & t_stat_sm >= 0 & D_pub_cb == 1) %>%
  select(journal_id, eye_catcher, exploratory, data_driven, peer, 
         data_availability, codes_availability, t_stat_sm)

Figure9 <- ggplot(data = Figure9_data, aes(x = t_stat_sm)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, center = 0.025, fill = "#F8766D") +
  stat_density(size = 0.8, kernel = 'epanechnikov', bw = 0.2, geom = "line", position = "identity", bounds = c(0, 11)) +
  theme_starwars() +
  ylab("Density") +
  xlab("|z-score|") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.6), breaks = c(0, 0.2, 0.4, 0.6), 
                     sec.axis = sec_axis(~., name = "Density", breaks = c(0, 0.2, 0.4, 0.6))) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 10),  breaks = c(0, 1.96, 10), labels = c("0", "1.96", "10")) +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust = 1, margin=unit(c(0,5,0,0), "mm")),
        axis.title.y.right = element_text(size = rel(1), angle = 90, vjust = 1, margin=unit(c(0,0,0,5), "mm")))

# Save RDP charts
facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0
ggsave("Figures/RDP_Figure1.png", Figure1, device = agg_png, scaling = 1.3,
       units = "mm", width = 360, height = 161, dpi = 500)

facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0

ggsave("Figures/RDP_Figure5.png", Figure5, device = agg_png, scaling = 1.3,
       units = "mm", width = 360, height = 161, dpi = 500)

ggsave("Figures/RDP_Figure6.png", Figure6, device = agg_png, scaling = 1.3,
       units = "mm", width = 220, height = 161, dpi = 500)

facet_number <- 0 #this is necessary because of the hack way I have implemented custom breaks
label_number <- 0

ggsave("Figures/RDP_Figure8.png", Figure8, device = agg_png, scaling = 1.3,
       units = "mm", width = 360, height = 161, dpi = 500)

ggsave("Figures/RDP_Figure9.png", Figure9, device = agg_png, scaling = 1.3,
       units = "mm", width = 220, height = 161, dpi = 500)



