
# -----------------------------------------------------------------------------
# This code imports the data and plots the line maps in section 4 of the paper.
#
# version: 26/12/2025
#
# -----------------------------------------------------------------------------
# Bubble size methodology:
# The bubble for each country represents the total liquidity line size in bn USD,
# averaged across borrowing and lending across all counterparties.
#
# For lines with a stated cap, the deal amount is converted to USD using
# contemporaneous exchange rates.
#
# For unlimited lines (notably the Fed's standing arrangements with the
# Bank of Canada, Bank of England, Bank of Japan, ECB, and Swiss National Bank),
# the cap is replaced with the maximum observed drawing in bn USD, sourced from
# central bank operational data and:
#   Bahaj, Reis (2020) "Central bank swap lines during the Covid-19 pandemic"
#   Bahaj, Reis (2022) "Central bank swap lines: Evidence on the Effect of
#                        the Lender of Last Resort"
#
# The pre-computed deal sizes in Input/Map_data/deal_size_YYYY.xlsx are provided
# as inputs so that the replication code can reproduce the figures without
# requiring access to the proprietary exchange rate data used for currency
# conversion.
# -----------------------------------------------------------------------------


# -----------------------------------------------------------------------------
#   STEP 1: prepare environment 

# Set directory from command line argument
args <- commandArgs(trailingOnly = TRUE)
if (length(args) > 0) {
  current_dir <- args[1]
} else {
  current_dir <- dirname(rstudioapi::getSourceEditorContext()$path)
}
parent_dir    <- dirname(current_dir)
new_path      <- file.path(parent_dir, "Data")
output_path   <- file.path(parent_dir, "Output", "figures")
setwd(new_path)

rm(current_dir, new_path, parent_dir)

# packages
packages <- c("tidyverse", "knitr", "haven", "readxl", "writexl", "dplyr",
              "lubridate", "tinytex", "colorspace", "magrittr", "moments",
              "maps", "geosphere", "assertthat", "igraph", "ggraph", "ggmap",
              "remotes", "extrafont", "ggnewscale", "scales")

installed_packages <- installed.packages()
new_packages <- packages[!(packages %in% installed_packages[, "Package"])] 
if (length(new_packages) > 0) {
  install.packages(new_packages, dependencies = TRUE, repos = "https://cloud.r-project.org")
}
lapply(packages, require, character.only = TRUE)
rm(new_packages, packages, installed_packages)


# Load inputs
ccc <- read_xlsx("Input/country_code_conversion.xlsx", sheet = "Sheet1") %>%
  select(ID, longitude, latitude) %>%
  rename(lon = longitude, lat = latitude)
network_data <- read_excel("liquidity_lines.xlsx", sheet = "LiquidityLines")
network_data$start_date <- as.Date(network_data$start_date)
network_data$end_date <- as.Date(network_data$end_date)

#-------------------------------------------------------------------------------
# Function to generate node and edge data for a given year
# Reads pre-computed deal sizes from Input/Map_data/ (produced by graph_replication_v6.R)
generate_year_data <- function(target_year, network_data, ccc) {

  year_start <- as.Date(paste0(target_year, "-01-01"))
  year_end <- as.Date(paste0(target_year, "-12-31"))

  active_network <- filter(network_data, start_date <= year_end & end_date >= year_start)

  country_data <- select(active_network, c(start_date, end_date, ISO_source, ISO_recipient))
  country_data <- country_data %>%
    mutate(ISO_source = ifelse(ISO_source == "EA", "DEU", ISO_source),
           ISO_recipient = ifelse(ISO_recipient == "EA", "DEU", ISO_recipient))

  country_data <- select(country_data, c(ISO_recipient, ISO_source))
  country_data <- unique(country_data)
  names(country_data)[names(country_data) == "ISO_recipient"] <- "target_country"
  names(country_data)[names(country_data) == "ISO_source"] <- "source_country"

  country_data$country_pair <- ifelse(country_data$source_country < country_data$target_country,
                                      paste(country_data$source_country, country_data$target_country, sep = "-"),
                                      paste(country_data$target_country, country_data$source_country, sep = "-"))
  country_data_unique <- unique(country_data[, c("country_pair", "source_country", "target_country")])
  country_data_unique <- country_data_unique[, -1]

  # Read pre-computed deal sizes
  deal_size <- read_xlsx(paste0("Input/Map_data/deal_size_", target_year, ".xlsx"))
  deal_size <- deal_size %>%
    mutate(ISO_source = ifelse(ISO_source == "EA", "DEU", ISO_source))

  dtf_edge <- select(country_data_unique, c(source_country, target_country))
  dtf_node <- union(unique(country_data$source_country), unique(country_data$target_country))
  dtf_node <- data.frame(ID = unique(dtf_node))

  dtf_node <- dtf_node %>%
    left_join(deal_size, by = c("ID" = "ISO_source")) %>%
    rename(amount = total_USD_amount) %>%
    mutate(amount = replace(amount, is.na(amount), 0)) %>%
    left_join(ccc, by = "ID") %>%
    rename(longitude = lon, latitude = lat)

  dtf_edge <- dtf_edge %>%
    rename(Source = source_country, Target = target_country)

  return(list(node = dtf_node, edge = dtf_edge))
}

#-------------------------------------------------------------------------------
# Generate all countries reference and colors
# (Need to do one pass to get all unique countries across all years)
all_nodes <- data.frame()
for (yr in c(2000, 2009, 2020, 2023, 2025)) {
  yr_data <- generate_year_data(yr, network_data, ccc)
  all_nodes <- rbind(all_nodes, yr_data$node)
}

dtf_all_countries <- all_nodes %>%
  select(ID, longitude, latitude) %>%
  unique() %>%
  filter(!is.na(longitude)) %>%
  mutate(region = case_when(
    ID == "BIS" ~ "BIS",
    longitude < -35 ~ "Americas",
    longitude < 35 & latitude > 32 ~ "Europe/Mediterranean",
    longitude < 35 & latitude <= 32 ~ "Africa",
    TRUE ~ "Asia"
  )) %>%
  mutate(region = as.factor(region))

predefined_colors <- c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
assign_region_colors <- function(dtf) {
  unique_regions <- levels(dtf$region)
  custom_order <- c("Americas", "Asia", "Europe/Mediterranean", "Africa")
  sorted_regions <- custom_order[custom_order %in% unique_regions]
  num_regions <- length(sorted_regions)
  region_colors <- setNames(predefined_colors[1:num_regions], sorted_regions)
  return(region_colors)
}
colours <- assign_region_colors(dtf_all_countries)

rm(all_nodes, yr_data)

#-------------------------------------------------------------------------------
# make_plot function
make_plot <- function(year, data_node, data_edge, suffix = "", breaks = c(500, 250, 5)) {
  
  data_node %<>% left_join(dtf_all_countries %>% select(ID, region), by = "ID")
  
  data_edgelist <- left_join(data_edge, data_node %>% select(ID, longitude, latitude, region),
                             by = c("Source" = "ID")) %>%
    rename(lon_start = longitude, lat_start = latitude, region_start = region)
  data_edgelist %<>% left_join(data_node %>% select(ID, longitude, latitude, region),
                               by = c("Target" = "ID")) %>%
    rename(lon_end = longitude, lat_end = latitude, region_end = region)
  
  index <- data_edgelist[, "region_start"] == data_edgelist[, "region_end"]
  data_curvelist_intra <- data_edgelist[index, ]
  data_curvelist_inter <- data_edgelist[!index, ]
  international_routes <- unique(data_curvelist_inter %>% select(region_start, region_end))
  
  linesize <- 0.25
  alpha <- 0.5
  
  maptheme <- theme(panel.grid = element_blank()) +
    theme(axis.text = element_blank()) +
    theme(axis.ticks = element_blank()) +
    theme(axis.title = element_blank()) +
    theme(legend.position = "inside") +
    theme(legend.position.inside = c(0.13, 0.26)) +
    theme(panel.background = element_rect(fill = "white")) +
    theme(plot.margin = unit(c(0, 0, 0.5, 0), 'cm')) +
    theme(legend.key = element_blank()) +
    theme(legend.key.size = unit(0.5, 'cm')) +
    theme(legend.title = element_text(size = 10)) +
    theme(legend.background = element_rect(fill = "transparent")) +
    theme(legend.box.background = element_blank())
  
  country_shapes <- geom_polygon(aes(x = long, y = lat, group = group),
                                 data = map_data('world') %>% filter(region != "Antarctica"),
                                 fill = "#E6E6E6", color = "#515151", size = 0.15)
  mapcoords <- coord_fixed(xlim = c(-150, 180), ylim = c(-80, 80))
  
  trans_size <- trans_new(name = "trans_size",
                          transform = function(x) x^0.7,
                          inverse = function(y) y^(1/0.7),
                          breaks = extended_breaks(),
                          minor_breaks = regular_minor_breaks(),
                          format = format_format(),
                          domain = c(-Inf, Inf))
  
  gg <- ggplot(data_node) +
    country_shapes +
    geom_curve(aes(x = lon_start, y = lat_start, xend = lon_end, yend = lat_end, col = region_start),
               data = data_curvelist_intra, size = linesize, curvature = 0.27, alpha = alpha) +
    mapcoords + maptheme +
    scale_colour_manual(guide = "none", values = colours)
  
  gg_inter <- gg
  if (nrow(international_routes) > 0) {
    for (i in 1:nrow(international_routes)) {
      route <- as.data.frame(international_routes[i, ])
      region_start <- route["region_start"][1, ]
      region_end <- route["region_end"][1, ]
      
      data_route <- data_curvelist_inter[
        data_curvelist_inter$region_start == route$region_start &
          data_curvelist_inter$region_end == route$region_end, ]
      
      curves <- NULL
      for (j in 1:nrow(data_route)) {
        current_edge <- data_route[j, ]
        geos <- c("lon_start", "lat_start", "lon_end", "lat_end")
        current_edge_geo <- as.numeric(current_edge[geos])
        names(current_edge_geo) <- geos
        
        interp <- gcIntermediate(c(current_edge_geo["lon_start"], current_edge_geo["lat_start"]),
                                 c(current_edge_geo["lon_end"], current_edge_geo["lat_end"]),
                                 n = 200, breakAtDateLine = TRUE)
        if (is.list(interp)) {
          great_circle <- greatCircle(c(current_edge_geo["lon_start"], current_edge_geo["lat_start"]),
                                      c(current_edge_geo["lon_end"], current_edge_geo["lat_end"]), n = 400)
          great_circle <- as.data.frame(great_circle)
          bEastwards <- TRUE
          if (max(interp[[1]][, "lon"]) < current_edge_geo["lon_start"]) bEastwards <- FALSE
          if (bEastwards) {
            interp <- great_circle[great_circle$lon < current_edge_geo["lon_start"] &
                                     great_circle$lon > current_edge_geo["lon_end"], ]
            onecurve <- cbind(interp, seq(1, 0, length.out = nrow(interp)), j)
          } else {
            interp <- great_circle[great_circle$lon > current_edge_geo["lon_start"] &
                                     great_circle$lon < current_edge_geo["lon_end"], ]
            onecurve <- cbind(interp, seq(0, 1, length.out = nrow(interp)), j)
          }
        } else {
          onecurve <- cbind(interp, seq(0, 1, length.out = nrow(interp)), j)
        }
        colnames(onecurve) <- c("lon", "lat", "pathing", "group")
        curves <- rbind(curves, onecurve)
      }
      assign(paste0("data_curve_inter_", i), as.data.frame(curves))
      
      gg_inter <- gg_inter +
        new_scale_color() +
        geom_path(aes(x = lon, y = lat, col = pathing, group = group),
                  data = get(paste0("data_curve_inter_", i)), lineend = "square",
                  size = linesize, alpha = alpha) +
        scale_colour_gradient(guide = "none",
                              low = colours[levels(region_start)[region_start]],
                              high = colours[levels(region_start)[region_end]])
    }
  }
  
  gg_inter <- gg_inter +
    geom_label(aes(x = longitude, y = latitude, label = ID),
               data = data_node %>% filter(ID %in% c("BIS")),
               hjust = 0, nudge_x = 4, nudge_y = 2, size = 2.5, label.size = NA,
               color = "black", fill = "white", alpha = 0, fontface = "bold")
  
  gg_inter <- gg_inter +
    new_scale_color() +
    geom_point(aes(x = longitude, y = latitude, size = amount, col = region),
               alpha = 0.5, stroke = 0.5) +
    scale_size_continuous(name = "Liquidity Line Size, in bn USD", range = c(1, 8),
                          trans = trans_size, breaks = breaks,
                          guide = guide_legend(override.aes = list(col = "darkblue", alpha = 0.5))) +
    scale_colour_manual(guide = "none", values = colours)
  
  temp_file_name <- paste0("map_temp_", year)
  if (hasArg(suffix) && length(suffix) == 1 && suffix != "") {
    temp_file_name <- paste0("map_temp_", year, "_", suffix)
  }
  
  output_pdf <- file.path(output_path, paste0(temp_file_name, ".pdf"))
  output_eps <- file.path(output_path, paste0(temp_file_name, ".eps"))
  output_png <- file.path(output_path, paste0(temp_file_name, ".png"))

  ggsave(filename = output_pdf, plot = gg_inter, device = "pdf")
  ggsave(filename = output_eps, plot = gg_inter, device = cairo_ps)
  ggsave(filename = output_png, plot = gg_inter, device = "png", dpi = 300)
  cat("Saved:", temp_file_name, ".pdf, .eps and .png\n")
}

#-------------------------------------------------------------------------------
# Generate plots for all years
target_years <- c(2000, 2009, 2020, 2023,2025)
breaks_by_year <- list(
  "2000" = c(5, 2.5, 0.5),
  "2009" = c(500, 250, 5),
  "2020" = c(500, 250, 5),
  "2023" = c(500, 250, 5),
  "2025" = c(500, 250, 5)
)

for (yr in target_years) {
  cat("Processing year:", yr, "\n")
  year_data <- generate_year_data(yr, network_data, ccc)
  make_plot(yr, year_data$node, year_data$edge, breaks = breaks_by_year[[as.character(yr)]])
  cat("Completed year:", yr, "\n\n")
}
