#Code for: Designing an efficient reference rate: Lessons from SOFIA
#Date April 2026
#
#This code creates variables used in the time series analysis



packagesNeeded <- c("tidyverse",
                    "data.table",
                    "mpi", 
                    "lubridate", 
                    "bankRun",
                    "ggrba",
                    "arphit",
                    "rbatools",
                    "mpi", 
                    "tseries", 
                    "spatstat.univar",
                    "janitor", 
                    "lmtest", 
                    "sandwich",
                    "minpack.lm", 
                    "msm", 
                    "ivreg", 
                    "fixest", 
                    "ggtext",
                    "purrr",
                    "writexl")

if(any(!packagesNeeded %in% dplyr::pull(tibble::as_tibble(installed.packages()), Package))){
  
  install.packages(packagesNeeded[!packagesNeeded %in% dplyr::pull(tibble::as_tibble(installed.packages()), Package)])
  
}

for(i in packagesNeeded){library(i, character.only = TRUE)}
rm(list=ls())


startDate <- "2022-01-01"
endDate <- Sys.Date()

startDate <- asDate(startDate)
endDate <- asDate(endDate)

#-----------------------------------------------------------------
#load AX transaction level data (MPI package)
rawASXData <- fetchASXRepoData(startDate = "2022-01-01", 
                               endDate = Sys.Date(), 
                               ERUser = TRUE)


#-------------------------------------------------------------------------------  


#import official SOFIA from ASX webpage 

directory <- tempdir()

saveDir <- paste0(directory, "/sofiaBeta.xlsx")

rbatools::set_proxy()

utils::download.file("https://www.asx.com.au/content/dam/asx/connectivity-and-data/sofia-beta-version.xlsx",
                     saveDir,
                     mode = "wb")

SOFIAData <- readxl::read_excel(saveDir, skip = 3) %>%
  # dplyr::select(1:3) %>%
  tidyr::drop_na() %>%
  dplyr::mutate(Date = as.Date(as.numeric(.data$Date), origin = "1899-12-30")) %>%
  dplyr::select("Date" = 'Date', "SOFIA" = 'SOFIA Overnight (%) VWAP',
                "SOFIAVolume" = 'Eligible Volume ($)',
                "SOFIAMedian" = "SOFIA Overnight  (%) Volume Weighted Median",
                "SOFIATrades" = "Eligible Trades",
                "SOFIAMin" = "Minimum Yield",
                "SOFIAMax" = "Maximum Yield") %>%
  dplyr::mutate(SOFIAVolume = .data$SOFIAVolume/1e9)

if(is.null(endDate)){endDate <- max(SOFIAData$Date)}
if(is.null(startDate)){startDate <- endDate}


#-------------------------------------------------------------------------------
#import some files that would otherwise be pulled from system 

cashRateDaily <- read.csv("C:/repo_market/repo_market/cashRateDaily.csv") %>% 
  as_tibble() %>% 
  mutate(date = as.Date(date)) %>% 
  select(-X)

businessDays <- read.csv("C:/repo_market/repo_market/businessDaySchedule.csv") %>% 
  as_tibble() %>% 
  mutate(days = as.Date(days), 
         shifted = as.Date(shifted))


SOFIAFiltered <- SOFIAData %>%
  dplyr::filter(dplyr::between(.data$Date, .env$startDate, .env$endDate))

#----------------------------------------------------------------------------------

list_buyer <- rawASXData %>%
  dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
  mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
  mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                     TRUE ~ 0))%>%
  filter(collateralDummy == 1,                                                        
                  repoTermBusDays == 1,                                                        
                  tradeDate == settleDate,                                                     
                  originalRepoType != "Call",                                                  
                  actualRepoType != "Call"                                                   
           ) %>%
  filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode)) %>%                #exclude RBA as counterparty
  group_by(buyerCode, parentBuyerCode, buyerName) %>%
    dplyr::summarise(volume = sum(leg1FV)) %>%
    dplyr::ungroup()%>% 
  dplyr::rename("code" = "buyerCode",
                "parentCode" = "parentBuyerCode",
                "name" = "buyerName")%>%
  subset(select = -volume)

  
list_seller <- rawASXData %>%
  dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
  mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
  mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                     TRUE ~ 0))%>%
  filter(collateralDummy == 1,                                                        
         repoTermBusDays == 1,                                                        
         tradeDate == settleDate,                                                     
         originalRepoType != "Call",                                                  
         actualRepoType != "Call"                                                   
  ) %>%
  filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode)) %>%                #exclude RBA as counterparty
  group_by(sellerCode, parentSellerCode, sellerName) %>%
  dplyr::summarise(volume = sum(leg1FV)) %>%
  dplyr::ungroup()%>% 
  dplyr::rename("code" = "sellerCode",
                "parentCode" = "parentSellerCode",
                "name" = "sellerName")%>%
  subset(select = -volume)


participants<-list_buyer%>%
  full_join(list_seller)%>%
  arrange(code)%>%
  group_by(parentCode)

 
  
  #concentration variable defined as concentration after data triming (bottom 25 percetile)
  #so first we need to find cutoff repo rate and proportion of the volume in cutoff

  cutoff <- rawASXData %>%
  dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
  mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
  mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                     TRUE ~ 0))%>%
  filter(collateralDummy == 1,                                                        
         repoTermBusDays == 1,                                                        
         tradeDate == settleDate,                                                     
         originalRepoType != "Call",                                                  
         actualRepoType != "Call"                                                   
  ) %>%
  filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode))%>%
  mutate(parentBuyerCode2=ifelse(parentBuyerCode == "BNPT", "BPSS", parentBuyerCode),
         parentBuyerCode2=ifelse(parentBuyerCode == "CHAM", "JPMA", parentBuyerCode),
         parentBuyerCode2=ifelse(parentBuyerCode == "CINA", "SALS", parentBuyerCode),
         parentBuyerCode2=ifelse(parentBuyerCode == "DBAL", "DEUT", parentBuyerCode),
         parentBuyerCode2=ifelse(parentBuyerCode == "HKBA", "HKSB", parentBuyerCode),
         
         parentSellerCode2=ifelse(parentSellerCode == "BNPT", "BPSS", parentSellerCode),
         parentSellerCode2=ifelse(parentSellerCode == "CHAM", "JPMA", parentSellerCode),
         parentSellerCode2=ifelse(parentSellerCode == "CINA", "SALS", parentSellerCode),
         parentSellerCode2=ifelse(parentSellerCode == "DBAL", "DEUT", parentSellerCode),
         parentSellerCode2=ifelse(parentSellerCode == "HKBA", "HKSB", parentSellerCode))%>%
  mutate(relPartyDummy = case_when(parentBuyerCode2 == parentSellerCode2 ~ 1,      
                                   TRUE ~ 0),
         relPartyVolume=relPartyDummy*leg1FV,
         dailyVol=ave(leg1FV,tradeDate, FUN = sum ))%>%
  mutate(marketshare= leg1FV/dailyVol,
         mssquare=marketshare^2)%>%
  group_by(tradeDate, repoRate) %>%                                                  
  dplyr::summarise(volume = sum(leg1FV), 
                   trades = n(),
                   relTrades = sum(relPartyDummy),
                   relVolume = sum(relPartyVolume),
                   sumMS= sum(marketshare),
                   sumMSS = sum(mssquare))  %>% 
  group_by(tradeDate) %>% 
  arrange(tradeDate, repoRate) %>%                                                
  mutate( dailyVol=ave(volume,tradeDate, FUN = sum ),
          share = volume/sum(volume)*100,                                       
          share2= volume/dailyVol*100,
          shareTrades = trades/sum(trades)) %>%                                      
  mutate(cumVolume = cumsum(share)) %>% 
  mutate(diff25 = cumVolume - 25,
         atmin=diff25/share) %>% 
  filter(diff25 >= 0)%>%   
  mutate(cutoffRate = ave(repoRate, tradeDate, FUN = min) )%>%
  filter(repoRate==cutoffRate)%>%
    subset(select = c(tradeDate, atmin, cutoffRate) )
  
#herfindahl 
  HHI <- rawASXData %>%
    dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
    mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
    mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                       TRUE ~ 0))%>%
    filter(collateralDummy == 1,                                                        
           repoTermBusDays == 1,                                                        
           tradeDate == settleDate,                                                     
           originalRepoType != "Call",                                                  
           actualRepoType != "Call"                                                   
    ) %>%
    arrange(tradeDate, repoRate) %>% 
    dplyr::left_join(cutoff, by = "tradeDate" ) %>%
    filter(repoRate >= cutoffRate)%>%
    mutate(newVolume= case_when(repoRate == cutoffRate ~ leg1FV*atmin,       #volume after cutting data for w. median
                                TRUE ~ leg1FV),
           dailyVol=ave(newVolume,tradeDate, FUN = sum ),
           marketShare=newVolume/dailyVol,
           step1=marketShare^2)%>%
    group_by(tradeDate) %>% 
    summarise(HHI=sum(step1))
  
  #herfindahl in full data
  
  HHIfull <- rawASXData %>%
    dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
    mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
    mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                       TRUE ~ 0))%>%
    filter(collateralDummy == 1,                                                        
           repoTermBusDays == 1,                                                        
           tradeDate == settleDate,                                                     
           originalRepoType != "Call",                                                  
           actualRepoType != "Call"                                                   
    ) %>%
    filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode))%>%
    mutate( dailyVol=ave(leg1FV,tradeDate, FUN = sum ),
            marketshare = leg1FV/dailyVol,
            mssquare=marketshare^2)%>%
    group_by(tradeDate) %>%                                                  
    dplyr::summarise(HHI = sum(mssquare))  
  
  
  

  
  
  
  
  
  #data for File 2, Figure 1, Figure 2

  timeseriesdata <- rawASXData %>%
    dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
    mutate(repoTermBusDays = originalTerm - adjustment) %>%                             
    mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                       TRUE ~ 0))%>%
    filter(collateralDummy == 1,                                                        
           repoTermBusDays == 1,                                                        
           tradeDate == settleDate,                                                     
           originalRepoType != "Call",                                                  
           actualRepoType != "Call"                                                   
    ) %>%
    filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode))%>%
    mutate(parentBuyerCode2=ifelse(parentBuyerCode == "BNPT", "BPSS", parentBuyerCode),
           parentBuyerCode2=ifelse(parentBuyerCode == "CHAM", "JPMA", parentBuyerCode),
           parentBuyerCode2=ifelse(parentBuyerCode == "CINA", "SALS", parentBuyerCode),
           parentBuyerCode2=ifelse(parentBuyerCode == "DBAL", "DEUT", parentBuyerCode),
           parentBuyerCode2=ifelse(parentBuyerCode == "HKBA", "HKSB", parentBuyerCode),
           
           parentSellerCode2=ifelse(parentSellerCode == "BNPT", "BPSS", parentSellerCode),
           parentSellerCode2=ifelse(parentSellerCode == "CHAM", "JPMA", parentSellerCode),
           parentSellerCode2=ifelse(parentSellerCode == "CINA", "SALS", parentSellerCode),
           parentSellerCode2=ifelse(parentSellerCode == "DBAL", "DEUT", parentSellerCode),
           parentSellerCode2=ifelse(parentSellerCode == "HKBA", "HKSB", parentSellerCode))%>%
    mutate(relPartyDummy = case_when(parentBuyerCode2 == parentSellerCode2 ~ 1,      
                                     TRUE ~ 0),
           relPartyVolume=relPartyDummy*leg1FV,
           dailyVol=ave(leg1FV,tradeDate, FUN = sum ))%>%
    mutate(marketshare= leg1FV/dailyVol,
           mssquare=marketshare^2)%>%
    group_by(tradeDate, repoRate) %>%                                                  
    dplyr::summarise(volume = sum(leg1FV), 
                     trades = n(),
                     relTrades = sum(relPartyDummy),
                     relVolume = sum(relPartyVolume),
                     sumMS= sum(marketshare),
                     sumMSS = sum(mssquare))  %>% 
    group_by(tradeDate) %>% 
    arrange(tradeDate, repoRate) %>%                                                
    mutate( dailyVol=ave(volume,tradeDate, FUN = sum ),
            share = volume/sum(volume)*100,                                       
            share2= volume/dailyVol*100,
            shareTrades = trades/sum(trades)) %>%                                      
    mutate(cumVolume = cumsum(share)) %>% 
    mutate(diff25 = cumVolume - 25) %>% 
    filter(diff25 >= 0) %>%                                                           #keep only transactions with market share above 25%
    mutate(newVolume = case_when(diff25 == min(diff25) ~ (diff25/share)*volume,       #volume after cutting data for w. median
                                 TRUE ~ volume))%>% 
    dplyr::summarise(SOFIAMedian = spatstat.univar::weighted.quantile(x = repoRate,   
                                                                      w = newVolume,
                                                                      probs = 0.5,
                                                                      type = 1),
                     SOFIA = round(weighted.mean(x = repoRate, w = newVolume), 6),
                     SOFIAMin = min(repoRate), 
                     SOFIAMax = max(repoRate), 
                     SOFIATrades = sum(trades),
                     SOFIAVolume = round(sum(newVolume, na.rm = TRUE), 0)/1e9,
                     RelatedTrades = sum(relTrades),
                     RelatedVolume = sum(relVolume),
                     MarketShare = sum(sumMS),
                     HuiHebel = (SOFIAMax-SOFIAMin)/SOFIAVolume,
                     dailyVolume=mean(dailyVol),
                     shareRelVol = (RelatedVolume/1e9)/SOFIAVolume)%>%
    dplyr::ungroup() %>%
    dplyr::mutate(Type = "Proposed change")%>%
    left_join(HHIfull, by = "tradeDate" )%>% 
    dplyr::rename("Date" = "tradeDate")  
      
      
      
    write_xlsx(timeseriesdata, "timeseriesdata.xlsx") 
