
# Diebold and Li ----------------------------------------------------------

# Calculate diebold and li coefficients using least squares
mp_shock_diebold_li <- function(df, lambda = 0.597761, tau_variable = "resid_maturity", yield_variable = "yield" ){
  
  # check tau_variable and yield_variable exist
  
  # Run regressiona and tidy
  reg_data <- df %>%
    mutate(tau = case_when(tau_variable == "duration" ~ duration,
                           TRUE ~ !!rlang::parse_expr(tau_variable)), # tau in years
           b1 = (1- exp(-lambda*tau))/(lambda*tau),
           b2 = b1 - exp(-lambda*tau)) %>%
    nest_by(date_time, name) %>%
    rowwise() %>% 
    mutate(reg_output = list(lm(paste0(yield_variable," ~ b1 + b2") , data = data) ),
           reg_tidy = list(broom::tidy(reg_output) %>% 
                             mutate(term = case_when(str_detect(term, "Intercept") ~ "b0",
                                                                            TRUE ~ term)  ) )) %>%
    unnest(reg_tidy)
  
  # Calculate shock
  diebold_li_shock  <- reg_data %>%
    pivot_wider(id_cols = c("date_time",  "term", "name"), names_from = "name", values_from = "estimate") %>%
    mutate(mp_shock =  post - pre) # surprise caclulation
  
  # reformat
  diebold_li_shock_wide <- diebold_li_shock %>%
    pivot_wider(id_cols = c("date_time", "term"), names_from = "term", values_from = "mp_shock") %>%
    ungroup() 
  
  return( diebold_li_shock_wide)
}


# Target and Path - Simple Approach ---------------------------------------
mp_shock_target_path <- function( short_tenor = 2 , long_tenor = 10, short_type = "ags", long_type = "ags", yield_changes = NULL, ois_changes = NULL, ...){
  
  # Import necessary data
  if ("ags" %in% c(short_type, long_type)  & is.null(yield_changes) ){
    yield_changes <- extract_ags_yield(...)
  } 
  
  if ("ois" %in% c(short_type, long_type) & is.null(ois_changes)){
    ois_changes <- extract_ois(...)
  }
  

# Calculate target for AGS -------------------------------------------------
  if (short_type == "ags"){
    yield_short <- extract_target_yield(yield_changes, target = short_tenor)
  } else if (short_type == "ois"){
    yield_short <- ois_changes %>% 
      transmute(date_time, name, yield = !!sym(paste0(short_tenor, "m")))
  }
  
  if (long_type == "ags"){
    yield_long <- extract_target_yield(yield_changes, target = long_tenor)
  } else if ( long_type == "ois"){
    yield_long <- ois_changes %>% 
      transmute(date_time, name, yield = !!sym(paste0(long_tenor, "m")))
  }
  
  
  
  yield_target <- yield_short %>% 
    pivot_wider(names_from = name, values_from = yield) %>% 
    transmute(date_time, target = post - pre)
  
  yield_slope <- left_join(yield_long %>% rename(yield_long = yield),
                         yield_short %>% rename(yield_short = yield), by = c("date_time", "name")) %>% 
    transmute(date_time, name, slope = yield_long - yield_short) %>% 
    pivot_wider(names_from = name, values_from = slope) %>% 
    transmute(date_time, slope = post - pre)
  
  
  
  yield_all <- full_join(yield_target, yield_slope, by = "date_time") %>%
    arrange(date_time) %>% 
    rbatools::remove_empty_rows(strict = T) %>% 
    {mutate(. ,  slope_orth_target =  lm(slope ~ target, data =  .)$residuals)} %>% 
    transmute(date_time,
              !!paste0("target_", short_tenor, short_type , "_", long_tenor, long_type ) := target,
              !!paste0("slope_", short_tenor, short_type , "_", long_tenor, long_type  ) := slope,
              !!paste0("slope_orth_target_", short_tenor, short_type , "_", long_tenor, long_type  ):= slope_orth_target)
  
  
  return(yield_all)
}


# PCA approach ------------------------------------------------------------

mp_shock_pca <- function(ois_changes){

    # Make data -----------------------------------------------------------------
    pca_data <-  ois_changes %>%
      mp_shock_ois() %>%
      select(date_time, starts_with("ois")) # remove lowly traded ois tenors

    pca_data <- pca_data[ rowSums(is.na(pca_data)) == 0 , ] # remove row if it contains an NA

    # Perform PCA
    pca <- prcomp(pca_data %>% select(-date_time), center = T, scale. = T )

    return(list(pca_output = pca, pca_factors = bind_cols( date_time = pca_data$date_time, as_tibble(pca$x)  )  ) )

}


mp_shock_pca_scale <- function(df, pc = c("PC1", "PC2"), scale = c("ois_1m", "ois_1y") ){

  scale  <- map2_dbl(pc, scale, ~ lm( as.formula(paste0(.y, "~ ",  .x  ) ), data = df )$coefficients[[.x]] ) %>%
    set_names(pc)

  scale_factors <- map_dfc(pc, ~ df[.]*scale[.] ) %>% set_names(paste0(stringr::str_to_lower( pc), "_scaled"))

  output <- bind_cols(df, scale_factors)

  return(output)

}




# GSS Rotation and Scaling ------------------------------------------------

mp_shock_pca_gss_rotate <- function(pca){
 # THE IMPACT OF ECB MONETARY POLICY DECISIONS AND COMMUNICATION ON THE YIELD CURVE (Appendix C for details)
# By  Claus Brand, Daniel Buncic and Jarkko Turunen

  pca_output <- pca$pca_output
  # first factor
  a1 <- pca_output$rotation["ois_1m", "PC1"]/(pca_output$rotation["ois_1m", "PC1"] + pca_output$rotation["ois_1m", "PC2"])
  a2 <- pca_output$rotation["ois_1m", "PC2"]/(pca_output$rotation["ois_1m", "PC1"] + pca_output$rotation["ois_1m", "PC2"])

  pc1_rotated <- a1*pca_output$x[ , "PC1"] + a2 * pca_output$x[ , "PC2"]

  # second factor

  factor_one_variance <- var(pca_output$x[ , "PC1"])
  factor_two_variance <- var(pca_output$x[ , "PC2"])

  b1 <- (-1 * a2 * factor_two_variance) /(a1*factor_one_variance - a2*factor_two_variance )
  b2 <- (a1*factor_one_variance)/(a1*factor_one_variance - a2*factor_two_variance )

  pc2_rotated <- b1*pca_output$x[ , "PC1"] + b2*pca_output$x[ , "PC2"]

  # combine

 return( bind_cols( date_time = pca$pca_factor$date_time,
                    tibble(pc1_rotated, pc2_rotated)  ))

}


mp_shock_scale <- function(df, pc = c("pc1_rotated", "pc1_rotated + pc2_rotated"), scale = c("ois_1m", "ois_12m")){

  names <- c(pc[1], str_extract(pc[2], "\\+.*") %>% str_remove("\\+") %>% str_squish())

  scale  <- map2(pc, scale, ~ lm( as.formula(paste0(.y, "~ ",  .x  ) ), data = df )$coefficients) %>%
    set_names(names)

  scale[[names[1]]] <- scale[[names[1]]] [[names[1] ]]
  scale[[names[2]]] <- scale[[names[2]]] [[ names[2] ]]/scale[[names[2] ]] [[ names[1] ]]

  scale_factors <- map_dfc(names(scale), ~ df[.]*scale[.] ) %>% set_names( str_to_lower( paste0(names(scale), "_scaled")))

  output <- bind_cols(df, scale_factors)

  return(output)

}


# OIS ---------------------------------------------------------------------


mp_shock_ois <- function(ois_changes){
  
  
  # Make wide and take difference to create mp Shock ----------------------------------------------------------------
  mp_shock <- ois_changes %>%
    select(-value) %>%
    pivot_longer( -c(date_time, name),  names_to = "tenor") %>%
    pivot_wider( names_from = name, values_from = value) %>%
    mutate(mp_shock = post - pre)
  
  # Make wide
  mp_shock_wide <- mp_shock %>%
    select(date_time, tenor,  mp_shock) %>%
    pivot_wider(names_from = "tenor", values_from = "mp_shock")
  
  colnames(mp_shock_wide) <- c("date_time", paste0("ois_", colnames(mp_shock_wide)[2:ncol(mp_shock_wide)]))
  
  return( mp_shock_wide)
  
}


mp_shock_ois_mean <- function(ois_changes){
  
  
  ois <- ois_changes %>%
    filter(name %in% c("pre", "post"))
  
  # MP Shock ----------------------------------------------------------------
  mp_shock <- ois %>%
    select(-value) %>%
    rbatools::remove_empty_rows(strict = T) %>%
    pivot_longer( -c(date_time, name) ,  names_to = "tenor") %>%
    group_by(date_time, name) %>%
    summarise(value = mean(value)) %>%
    ungroup %>%
    filter(value != 0) %>%
    pivot_wider(id_cols = c("date_time")) %>%
    transmute(date_time,
              ois_mean = post - pre)
  
  return(mp_shock)
  
  
}

# Cash rate changes -------------------------------------------------------
mp_shock_cashrate_change <- function(){
 output <-  rbatools::board_decisions() %>%
    arrange(effective_date) %>%
    mutate(change = str_remove(change, "\\+") ) %>%
    transmute( decision_date= case_when(lubridate::year(effective_date) >= 2008 ~ effective_date - 1 ,
                                      TRUE  ~ effective_date),
               date_time = case_when(lubridate::year(decision_date) >= 2008 ~ ymd_hm(paste0(decision_date, " 14:30")),
                                              TRUE ~ ymd_hm(paste0(decision_date, " 09:30"))),
               cashrate_change =  as.numeric(change))

  # manual inputs for two months in early 1990s where changes were around a range;
  # take mid-point of range
  output$cashrate_change[output$decision_date == "1990-01-23"] <- -0.75
  output$cashrate_change[output$decision_date == "1990-04-04"] <- -1.25
  
  output <- output %>% transmute(date_time, cashrate_change)

  return(output)
}

# Cash rate level ---------------------------------------------------------

mp_shock_cashrate_level <- function(){
  output <-   rbatools::board_decisions() %>%
    arrange(effective_date) %>%
    mutate(target = as.numeric(str_remove(target, " to .*"))) %>%
    transmute( decision_date = case_when(lubridate::year(effective_date) >= 2008 ~ effective_date - 1 ,
                                      TRUE  ~ effective_date),
               decision_date_time = case_when(lubridate::year(decision_date) >= 2008 ~ ymd_hm(paste0(decision_date, " 14:30")),
                                              TRUE ~ ymd_hm(paste0(decision_date, " 09:30"))),
               cashrate_level =  case_when( decision_date <= "1990-07-04" ~ as.numeric(target) + 0.25,
                                            TRUE ~ as.numeric(target)))

  
  output <- output %>% transmute(date_time = decision_date_time, cashrate_level)
  
  return(output)
}



# Combine all HF shocks ---------------------------------------------------
mp_shock_combine <- function(pre = 30, post = 90, open_window = 60){

  ois_changes <- extract_ois(pre = pre, post = pre, open_window = open_window) %>%
    select(-`1w`, -`18m`, -`24m`, -`36m`) 
  yield_changes <- extract_ags_yield(pre = pre, post = pre, open_window = open_window)
  
  
  # Augment data with OIS AND Yields and do PCA
  yield_long <- map_dfr(2:10, function(x) extract_target_yield(yield_changes, target = x) %>% mutate(tenor = paste0(x * 12, "m" )))
  yield_ois_changes <- left_join( ois_changes, 
                                  yield_long %>% pivot_wider(names_from = tenor, values_from = yield), 
                                  by = c("date_time", "name"))  
  pca_augmented <- mp_shock_pca(yield_ois_changes)$pca_factors
  colnames(pca_augmented) <- c("date_time", paste0(colnames(pca_augmented)[2:ncol(pca_augmented)], "_augmented"))
  
  
  
  all_shocks <- list(
    cash_rate_change = mp_shock_cashrate_change(),
    cash_rate_level = mp_shock_cashrate_level(),
    ois = mp_shock_ois(yield_ois_changes ),
    ois_mean = mp_shock_ois_mean(ois_changes),
    pca = mp_shock_pca(ois_changes)$pca_factors,
    pca_augmented = pca_augmented,
    target_slope_2ags_10ags = mp_shock_target_path(short_tenor = 2, long_tenor = 10, short_type = "ags", long_type = "ags", yield_changes = yield_changes),
    target_slope_3ois_2ags = mp_shock_target_path(short_tenor = 3, long_tenor = 2, short_type = "ois", long_type = "ags", yield_changes = yield_changes, ois_changes = ois_changes),
    dl = mp_shock_diebold_li(yield_changes))

  # add gss rotated shocks
  all_shocks$gss  <- mp_shock_pca(ois_changes) %>% mp_shock_pca_gss_rotate()
  all_shocks$gss_augmented <- mp_shock_pca(yield_ois_changes) %>% mp_shock_pca_gss_rotate()
  colnames( all_shocks$gss_augmented) <- c("date_time", "pc1_augmented_rotated", "pc2_augmented_rotated")

  # add scaled ois shocks
  all_shocks <- purrr::reduce(all_shocks, left_join) %>%
    mp_shock_scale(pc = c("PC1", "PC1 + PC2"), scale = c("ois_1m", "ois_12m")) %>%
    mp_shock_scale(pc = c("PC1_augmented", "PC1_augmented + PC2_augmented"), scale = c("ois_1m", "ois_12m")) %>% 
    mp_shock_scale(pc =  c("pc1_rotated", "pc1_rotated + pc2_rotated"), scale = c("ois_1m", "ois_12m")) %>% 
    mp_shock_scale(pc =  c("pc1_augmented_rotated", "pc1_augmented_rotated + pc2_augmented_rotated"), scale = c("ois_1m", "ois_12m"))

  return(all_shocks)
}


# Save --------------------------------------------------------------------
mp_shock_save <- function(){
  output <- mp_shock_combine()
  saveRDS(output, file = "./output/high-frequency-surprises.RDS")
  cat("High frequency surprises saved: ./output/high-frequency-surprises.RDS\n")
  return(output)
}
