# Star Wars at Central Banks Script - Create hypothetical distribution charts
# Last Updated: 29 September 2020

library(tidyverse)
library(ggplot2)
library(ragg)
source("Programs/E_paper-formatted_graphs/starwars_graph_theme.R") # ggplot theme for graphs

### Z-score chart - RDP Figure 2
# Generate hypothetical data
hypothetical_z <- tibble(z_score = rep(seq(0.5, 4.5, by = 1), times = 3),
                         distribution_type = c(rep("P[z|disseminated]", times = 5), 
                                               rep("P[disseminated|z]", times = 5),
                                               rep("P[z]", times = 5)),
                         probability = c(0.3,0.2,0.3,0.1,0.1,
                                         0.15,0.15,0.3,0.3,0.3,
                                         0.6/1.5,0.4/1.5,0.3/1.5,0.1/1.5,0.1/1.5))

# Set up the bias series to appear in correct facet order when plotting
hypothetical_z[["distribution_type"]] <- factor(hypothetical_z$distribution_type, 
                                                levels = c("P[z|disseminated]",
                                                           "P[disseminated|z]",
                                                           "P[z]"))

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

Figure2 <- ggplot(data = hypothetical_z) +
  geom_col(mapping = aes(x = z_score, y = probability)) +
  theme_starwars() +
  facet_grid(cols = vars(distribution_type)) +
  aes(fill = distribution_type) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5), sec.axis = sec_axis(~., name = "Probability")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 5), breaks = figure2_breaks) +
  xlab("|z-score|") +
  ylab("Probability") +
  geom_text(data = tibble(distribution_type = (unique(hypothetical_z[["distribution_type"]])), 
                          label = (c("P[z|disseminated]", "P[disseminated|z]", "P[z]"))),
            mapping = aes(label = label, x = 2.5, y = 0.47),
            size = 5.5,
            parse = FALSE) +
  theme(legend.position = "none",
        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")))

ggsave("Figures/RDP_Figure2.png", Figure2, device = agg_png, scaling = 1.3,
       units = "mm", width = 350, height = 161, dpi = 500)

### P-curve chart
# Generate hypothetical data
hypothetical_p <- tibble(p_value = rep(seq(0.005, 0.045, by = 0.01), times = 3),
                         distribution_type = c(rep("Ambiguous P[p]; form 1", times = 5), 
                                               rep("Ambiguous P[p]; form 2", times = 5),
                                               rep("Biased P[p]; form 3", times = 5)),
                         probability = c(0.2,0.2,0.2,0.2,0.2,
                                         0.3,0.25,0.2,0.15,0.1,
                                         0.1,0.15,0.2,0.25,0.3))

# Set up the bias series to appear in correct facet order when plotting
hypothetical_p[["distribution_type"]] <- factor(hypothetical_p$distribution_type, 
                                           levels = c("Ambiguous P[p]; form 1",
                                                      "Ambiguous P[p]; form 2",
                                                      "Biased P[p]; form 3"))

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

Figure3 <- ggplot(data = hypothetical_p) +
  geom_col(mapping = aes(x = p_value, y = probability), fill = "#619CFF") +
  theme_starwars() +
  facet_grid(cols = vars(distribution_type)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.4), sec.axis = sec_axis(~., name = "Probability")) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 0.05), breaks = figure3_breaks) +
  xlab("p-value") +
  ylab("Probability") +
  geom_text(data = tibble(distribution_type = (unique(hypothetical_p[["distribution_type"]])), 
                              label = (c("Ambiguous P[p];\n form 1", "Ambiguous P[p];\n form 2", "Biased P[p];\n form 3"))),
            mapping = aes(label = label, x = 0.025, y = 0.365),
            size = 5.5,
            parse = FALSE) +
  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")))
 
ggsave("Figures/RDP_Figure3.jpg", Figure3, device = agg_jpeg, scaling = 1.3,
       units = "mm", width = 350, height = 161, dpi = 500)



  
  

