# Star Wars at Central Banks Script - Create p-curve outputs
# Last Updated: 27 January 2021

library(ggplot2)
library(ragg)

set.seed(1) # this is to ensure reproducibility of the pseudo-random sampling of test statistics
source("Programs/D_p-curve_analysis/D1_construct_inputs.R")
source("Programs/D_p-curve_analysis/D2_udf.R")
source("Programs/D_p-curve_analysis/D3_apply_method.R")
source("Programs/E_paper-formatted_graphs/starwars_graph_theme.R") # ggplot theme for graphs


# Create main graph for the paper (Figure 4)
groupings <- c("Brodeur sample", "Full sample", "Reverse-causal research", "Data-driven research", "Restricted sample")
full_curve <- map2(full_curve, groupings, function (x, y) mutate(.data = x, grouping = y)) # identification of each sub-group in each dataframe is necessary to draw this graph

main_pcurve_data <- rbind(full_curve[["full_cb"]], full_curve[["full_brodeur"]]) %>%
  mutate(bins = as.numeric(as.character(bins)),
         grouping = factor(grouping, levels = c("Full sample", "Brodeur sample")),
         `Observed p-curve` = `Observed p-curve`/100)

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

Figure4 <- ggplot(main_pcurve_data, aes(x = bins)) +
  geom_col(aes(y = `Observed p-curve`), fill = "#F8766D", na.rm = TRUE, size = 1) +
  theme_starwars() +
  facet_grid(cols = vars(grouping)) +
  geom_text(aes(label = groupings, x = 0.025, y = c(0.9, 0.94)),
            data = data.frame(grouping = unique(main_pcurve_data[["grouping"]]), groupings = c("Central bank discussion\npapers (2000-2019)", "Top journals (2005-2011)")),
            size = 5.5) +
  ylab("Share of significant p-values") +
  xlab("p-value") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1), sec.axis = sec_axis(~., name = "Share of significant p-values")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 0.05), breaks = figure3_breaks) +
  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")))

# Appendix graph 1 - subgroup chart - alpha = 0.05 (exploratory, data-driven, and total minus exploratory and data driven)
groupings <- c("Brodeur sample", "Full sample", "Reverse-causal research", "Data-driven models", "Preferred sample")
full_curve <- map2(full_curve, groupings, function (x, y) mutate(.data = x, grouping = y)) # identification of each sub-group in each dataframe is necessary to draw this graph

ag1_data <- rbind(full_curve[["exploratory_cb"]], full_curve[["datadriven_cb"]], full_curve[["total_less_exp_and_dd"]]) %>%
  mutate(bins = as.numeric(as.character(bins)),
         grouping = factor(grouping, 
                           levels = c("Data-driven models", 
                                      "Reverse-causal research", 
                                      "Preferred sample")),
         `Observed p-curve` = `Observed p-curve`/100)

ag1 <- ggplot(ag1_data, aes(x = bins)) +
  geom_line(aes(y = `Observed p-curve`, color = grouping), na.rm = TRUE, size = 1) +
  geom_point(aes(y = `Observed p-curve`, colour = grouping), na.rm = TRUE, size = 3) +
  theme_starwars() +
  ylab("Share of significant p-values") +
  xlab("p-value") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1), sec.axis = sec_axis(~., name = "Share of significant p-values")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 0.05)) +
  theme(legend.position = c(0.97, 1.03)) +
  scale_color_manual(values = c("blue4", "gold1", "hotpink2")) +
  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")))

# Appendix graph 2 - subgroup chart - alpha = 0.01 (exploratory, data-driven, and total minus exploratory and data driven)
triple_star_curve <- map2(triple_star_curve, groupings, function (x, y) mutate(.data = x, grouping = y)) # identification of each sub-group in each dataframe is necessary to draw this graph

ag2_data <- rbind(triple_star_curve[["exploratory_cb"]], triple_star_curve[["datadriven_cb"]], triple_star_curve[["total_less_exp_and_dd"]]) %>%
  mutate(bins = as.numeric(as.character(bins)),
         grouping = factor(grouping, 
                           levels = c("Data-driven models", 
                                      "Reverse-causal research", 
                                      "Preferred sample")),
         `Observed p-curve` = `Observed p-curve`/100)

ag2 <- ggplot(ag2_data, aes(x = bins)) +
  geom_line(aes(y = `Observed p-curve`, color = grouping), na.rm = TRUE, size = 1) +
  geom_point(aes(y = `Observed p-curve`, colour = grouping), na.rm = TRUE, size = 3) +
  theme_starwars() +
  ylab("Share of significant p-values") +
  xlab("p-value") +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1), sec.axis = sec_axis(~., name = "Share of significant p-values")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 0.01), breaks = pretty(ag2_data[["bins"]], n = 5)) +
  theme(legend.position = c(0.97, 1.03)) +
  scale_color_manual(values = c("blue4", "gold1", "orchid3")) +
  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
ggsave("Figures/RDP_Figure4.png", Figure4, device = agg_png, scaling = 1.3,
       units = "mm", width = 260, height = 161, dpi = 500)

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

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