#Code for: Designing an efficient reference rate: Lessons from SOFIA
#Date April 2026
#
#This code creates variables SOFIA with 2 std deviation cut

#step 1. set cutoff rate
#step 2. create 2 std deviations of this post cutoff value
#step 3. delete the extreme values of sofia
#step 4. how many observations deleted? ->222
#step 5. create new synthetic sofia
#step 6. compare new synthetic sofia to ASX sofia
#step 7. compare to 95th trim







# Packages ----------------------------------------------------------------

dir <- "C:/repo_market/repo_market"

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

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)}

#-------------------------------------------------------------------------------
#set date

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

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


#-------------------------------------------------------------------------------  
#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(.data$Date)) %>%
  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 (no)

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))

#-------------------------------------------------------------------------------  
#get raw repo data


rawASXData <- fetchASXRepoData(startDate = "2022-01-01", 
                               endDate = Sys.Date(), 
                               ERUser = TRUE)

#  businessDays <- tibble::tibble(days = generateBusinessDays(startDate = min(SOFIAData$Date),
#                                                             endDate = max(SOFIAData$Date)))   %>%
#    dplyr::mutate(shifted = shiftBusinessDays(vecDate = days, shiftBy = 1),
#                  adjustment = as.numeric(shifted - days - 1))

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


#step 1. set cutoff rate



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) )



#step 2. create 2 std deviations of this post cutoff value


stddev<-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))%>%
  group_by(tradeDate) %>%  
  mutate(mean_value_by_day = mean(repoRate, na.rm = TRUE),
    sd_value_by_day = sd(repoRate, na.rm = TRUE)) %>%
  ungroup()%>%
  filter(repoRate>=mean_value_by_day+2*sd_value_by_day)%>%
  mutate(std2=1)%>%
  select(uniqueID, std2)


#step 5. create new synthetic sofia



SOFIANew <- rawASXData %>%
  dplyr::left_join(stddev, by = c("uniqueID"))%>%
  mutate(
    std2 = ifelse(is.na(std2), 0, std2))%>%
  filter(std2 != 1)%>%
  dplyr::left_join(businessDays, by = c("tradeDate" = "days")) %>%
  mutate(repoTermBusDays = originalTerm - adjustment) %>%                             #need to check that
  mutate(collateralDummy = case_when(issuerType %in% c("Govt", "Semi-Govt") ~ 1,      #create collateral dummy (AGS and Semi)
                                     TRUE ~ 0)) %>%
  filter(collateralDummy == 1,                                                        #keep only GC1
         repoTermBusDays == 1,                                                        #check that
         tradeDate == settleDate,                                                     # same day settlement
         originalRepoType != "Call",                                                  #exclude call repos
         actualRepoType != "Call"                                                    #exclude call repos
  ) %>% 
  filter(., !grepl("RBAA", buyerCode), !grepl("RBAA", sellerCode)) %>%                #exclude RBA as counterparty
  group_by(tradeDate, repoRate) %>%                                                   #question2
  summarise(volume = sum(leg1FV), 
            trades = n()) %>% 
  group_by(tradeDate) %>% 
  arrange(tradeDate, repoRate) %>%                                                   #order by repo rate
  mutate(share = volume/sum(volume)*100,                                            #volume market share
         shareTrades = trades/sum(trades)) %>%                                      #transaction market share
  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) %>%
  dplyr::ungroup() %>% 
  dplyr::rename("Date" = "tradeDate") %>% 
  dplyr::mutate(Type = "excl. >2 std dev")

write_xlsx(SOFIANew, "SOFIA_minus_outliers.xlsx")


#step 6. compare new synthetic sofia to ASX sofia

#add ASX version of SOFIA to compare
SOFIASplitCalc <- SOFIANew %>% 
  bind_rows(SOFIAData %>% mutate(Type = 'ASX')) %>% 
  dplyr::left_join(cashRateDaily, by = c("Date" = "date")) %>%
  dplyr::mutate(SOFIAspread = round((.data$SOFIA - round(ESRate, 4))*100, 2),
                SOFIAspreadMedian = (.data$SOFIAMedian - .data$ESRate)*100,
                SOFIAspreadMin = (.data$SOFIAMin - .data$ESRate)*100,
                SOFIAspreadMax = (.data$SOFIAMax - .data$ESRate)*100) %>% 
  filter(Date >= as.Date("2022-01-01"), 
         Date <= as.Date("2025-03-30")) %>% 
  filter(!Date %in% as.Date(c("2022-07-29", "2022-09-30", "2023-08-04", "2023-10-02", "2023-08-07", "2023-09-29")))



#step 7. compare to 95th trim

#plot graphs of both SOFIA
SOFIANewGraph111 <- arphit::arphitgg(layout = "3h") +
  arphit::agg_line(data = SOFIASplitCalc,
                   aes = arphit::agg_aes(x = .data$Date, y = .data$SOFIAspread, group = Type),
                   colour = c(arphit::RBA["Red2"], arphit::RBA["Blue9"], arphit::RBA["Grey10"]),
                   panel = "1") +
  arphit::agg_line(data = SOFIASplitCalc,
                   aes = arphit::agg_aes(x = .data$Date, y = .data$SOFIAVolume, group = Type),
                   colour = c(arphit::RBA["Red2"], arphit::RBA["Blue9"], arphit::RBA["Grey10"]),
                   panel = "3") +
  arphit::agg_legend(x= "topleft", ncol = 1) + 
  
  arphit::agg_line(data = SOFIASplitCalc,
                   aes = arphit::agg_aes(x = .data$Date, y = .data$SOFIATrades, group = Type),
                   colour = c(arphit::RBA["Red2"], arphit::RBA["Blue9"], arphit::RBA["Grey10"]),
                   panel = "5") +
  arphit::agg_units("bps", panel = "1") +
  arphit::agg_units("$b", panel = "3") +
  arphit::agg_units("no", panel = "5") +
  arphit::agg_title("") +
  arphit::agg_title("spread to ES rate", panel = "1") +
  arphit::agg_title("eligible volume", panel = "3") +
  arphit::agg_title("eligible trades", panel = "5") +
  arphit::agg_xlim(yearFrac("2021-12-01"), NA, panel = "1") +
  arphit::agg_xlim(yearFrac("2021-12-01"), NA, panel = "3") +
  arphit::agg_xlim(yearFrac("2021-12-01"), NA, panel = "5") +
  arphit::agg_source(c("ASX", "RBA")) 
# arphit::agg_footnote(paste0("ASX's official beta estimate for SOFIA. ", "SOFIA estimates up until ", max(SOFIAData$Date), "."))

agg_draw(SOFIANewGraph111,  "C:\\repo_market\\repo_market\\sofia_excl_outliers.png")

