library(tidyverse)
library(lubridate)
library(modelr)
library(stringr)
library(data.table)
library(openxlsx)
library(sandwich)
library(lmtest)
library(texreg)
library(rugarch)
library(RcppRoll)
library(xts)

################################################################################ Code for generating Figures 2, 3, and 5
create.ts.plot <- function(df, plotted.var) {
  df %>% ggplot(aes(x = Date)) + geom_line(aes(y = {{plotted.var}}), color = "blue", size = 1, alpha = 0.5) +
         xlab("") + ylab(deparse(substitute(plotted.var))) +  geom_vline(xintercept = monetary.policy.sur$Date, linetype = 2, alpha = 0.25) +
         theme(axis.text = element_text(size = 10), axis.title = element_text(size = 10), text = element_text(size = 8), panel.background = element_rect(fill = "white"),  axis.line = element_line(size = 1, color = "black"))
}
create.ts.plot.2.series <- function(df, plotted.var.1, plotted.var.2) {
  df %>% ggplot(aes(x = Date)) + geom_line(aes(y = {{plotted.var.1}}), color = "blue", size = 1, alpha = 0.5) +
         geom_line(aes(y = {{plotted.var.2}}), color = "red", size = 2, alpha = 0.5) +
         xlab("") + ylab( str_remove_all(deparse(substitute(plotted.var.1)), "\\dm\\.") ) +  geom_vline(xintercept = monetary.policy.sur$Date, linetype = 2, alpha = 0.25) +
         theme(axis.text = element_text(size = 10), axis.title = element_text(size = 10), text = element_text(size = 8), panel.background = element_rect(fill = "white"),  axis.line = element_line(size = 1, color = "black"))
    }
generate.ts.plots.all <- function(dta.option.quantities, page = 1) {
  main.ind.plot <- dta.option.quantities %>% create.ts.plot(Main.Spot) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  main.iv.plot <- dta.option.quantities %>% create.ts.plot.2.series(Main.1m.IV, Main.3m.IV) + ylim(0.35, 1.6)
  main.cal.plot <- dta.option.quantities %>% create.ts.plot(Main.Calendar) + ylim(-0.35, 0.15)
  cross.ind.plot <- dta.option.quantities %>% create.ts.plot(Crossover.Spot) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  cross.iv.plot <- dta.option.quantities %>% create.ts.plot.2.series(Crossover.1m.IV, Crossover.3m.IV) + ylim(0.35, 1.6)
  cross.cal.plot <- dta.option.quantities %>% create.ts.plot(Crossover.Calendar) + ylim(-0.35, 0.15)
  if(page == 1) plots.out <- (main.ind.plot / main.iv.plot) | (cross.ind.plot / cross.iv.plot)
  else plots.out <- ( main.cal.plot) | ( cross.cal.plot) 
  plots.out + plot_annotation(caption = "Data: courtesy of J.P.Morgan Chase & Co., Copyright 2022 - Morgan Markets")
}

extract.ym.for.plot <- function(date) {
  substring(paste(year(date), lubridate::month(date, label = TRUE), sep = '-'), first = 3)
}
generate.col.plot <- function(variable.str = "IV") { # event.reactions is a global variable
  aux.data <- event.reactions.3m %>% select(Date, contains(variable.str)) %>% pivot_longer(!Date, names_to = "main.cross", values_to = "val") %>%
                                     left_join(event.reactions %>% select(Date, hawkish)) %>% mutate(main.cross = str_remove_all(main.cross, "\\.[13]m"))
  event.reactions %>% select(Date, contains(variable.str)) %>% pivot_longer(!Date, names_to = "main.cross", values_to = "val") %>%
                      left_join(event.reactions %>% select(Date, hawkish)) %>% mutate(main.cross = str_remove_all(main.cross, "\\.[13]m")) %>%
                      ggplot(aes(x = as.character(Date), y = val, color = main.cross, fill = hawkish, alpha = main.cross)) + geom_col(position = "dodge") +
                      scale_color_manual(values = c("black", "black")) + scale_alpha_discrete(range = c(1,0.5)) +
                      geom_point(data = aux.data,  aes(x = as.character(Date), y = val), position = position_dodge(width = 1), size = 3, shape = 17, alpha = 1) +
                      scale_fill_manual(values = c("blue", "red", "yellow")) +
                      theme(axis.text = element_text(size = 10), axis.title = element_text(size = 12), text = element_text(size = 8), legend.position = "none",
                            panel.background = element_rect(fill = "white"),  axis.line = element_line(size = 1, color = "black")) +
                      xlab("") + ylab("bps") + scale_x_discrete(labels = extract.ym.for.plot(event.reactions$Date))
}
generate.col.plot.for.spot.or.calendar <- function(variable.str = "Spot") {
  event.reactions %>% select(Date, contains(variable.str)) %>% pivot_longer(!Date, names_to = "main.cross", values_to = "val") %>%
                      left_join(event.reactions %>% select(Date, hawkish)) %>%
                      ggplot(aes(x = as.character(Date), y = val, color = main.cross, fill = hawkish, alpha = main.cross)) + geom_col(position = "dodge") +
                      scale_color_manual(values = c("black", "black")) + scale_alpha_discrete(range = c(1,0.5)) + scale_fill_manual(values = c("blue", "red", "yellow")) +
                      theme(axis.text = element_text(size = 10), axis.title = element_text(size = 12), text = element_text(size = 8), legend.position = "none",
                            panel.background = element_rect(fill = "white"),  axis.line = element_line(size = 1, color = "black")) +
                      xlab("") + ylab("bps") + scale_x_discrete(labels = extract.ym.for.plot(event.reactions$Date))
}

################################################################################ Forecasts of vols under P
PATH.TO.INPUTS <- "" # ALL paths to data files are set to some generic locations
cds.spreads <- fread(paste0(PATH.TO.INPUTS,'/jpm/itraxx_indices_eod.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
cds.spreads <- cds.spreads %>% select(matches("(Mid)|(date)")) %>% relocate(date)
names(cds.spreads) <- c("date", "Main", "Crossover")
cds.sample <- cds.spreads %>% filter(date > ymd("2016-01-01"), date < ymd("2023-01-01"))
cds.log.R <- cds.spreads %>% mutate(across(where(is.numeric), function(x) c(NA, diff(log(x))))) %>% drop_na()
avg.var.P <- cds.log.R %>% select(date) %>% filter(date >= ymd("2016-01-01"))
avg.var.P <- avg.var.P %>% mutate(Main.1m.avg.var.P = 1, Main.3m.avg.var.P = 1, Crossover.1m.avg.var.P = 1, Crossover.3m.avg.var.P = 1)
spot.Main <- cds.log.R %>% select(date, Main) %>% rename(Spot = Main)
spot.Cross <- cds.log.R %>% select(date, Crossover) %>% rename(Spot = Crossover)
fit.GARCH <- function(df, sel.date, n.est.window = 500, spec = "sGARCH") {
  garch.spec <- ugarchspec(distribution.model = "std", variance.model = list(model = spec), mean.model = list(armaOrder = c(5,0)))
  dta.loc <- df %>% filter(date <= sel.date) %>% drop_na() %>% arrange(date) %>% tail(n.est.window)
  ugarchfit(spec = garch.spec, data = xts(dta.loc$Spot, order.by = dta.loc$date ), solver = "hybrid")   
}
forecast_avg_var <- function(fitted.model,  f.window = 22) {
  forecasts <- sigma(ugarchforecast(fitted.model, n.ahead = f.window))
  return(mean(forecasts^2)*252)
}
for(d in avg.var.P$date) {
  fitted.model <- fit.GARCH(spot.Main, as.Date(d))
  avg.var.P[avg.var.P$date == d,"Main.1m.avg.var.P"] <- forecast_avg_var(fitted.model, f.window = 22)
  avg.var.P[avg.var.P$date == d,"Main.3m.avg.var.P"] <- forecast_avg_var(fitted.model, f.window = 66)
  fitted.model <- fit.GARCH(spot.Cross, as.Date(d))
  avg.var.P[avg.var.P$date == d,"Crossover.1m.avg.var.P"] <- forecast_avg_var(fitted.model, f.window = 22)
  avg.var.P[avg.var.P$date == d,"Crossover.3m.avg.var.P"] <- forecast_avg_var(fitted.model, f.window = 66)
}
################################################################################ MFIV calculations
extrapol.iv <- function(df) {
  F <- df %>% mutate(F = 1/ (K.div.F / K)) %>% head(1) %>% pull(F)
  mdl.dta <- df %>% mutate(MN = - log(K.div.F)) %>% mutate(MN_sq = MN^2)
  mdl <- lm(vol ~ MN + MN_sq, data = mdl.dta)
  data.frame(K.div.F = c(seq(0.5, 1.50, 0.01))) %>% mutate(K = K.div.F*F, MN = - log(K.div.F), r = df$r[1]) %>% mutate(MN_sq = MN^2) %>% add_predictions(mdl) %>% rename(vol = pred) %>% mutate(vol = pmax(vol, 0.0))
}
black.price <- function(K.div.F, K, T, impl.vol, r) {
  d.1 <- (-log(K.div.F) + impl.vol^2*T/2) / (impl.vol*sqrt(T))
  d.2 <- d.1 - impl.vol*sqrt(T)
  (K.div.F > 1) * K*(1/K.div.F*pnorm(d.1) - pnorm(d.2)) + (K.div.F <= 1) * (K*(pnorm(-d.2) - 1/K.div.F*pnorm(-d.1)))
}
mfiv <- function(opt.price, K, T, r) {
  delta.K = K[2]-K[1]
  sum(opt.price * delta.K / K^2) * 2 / T
}
gen.mfiv.time.series <- function(smpl, T = 1/12) {
  smpl <- smpl %>% group_by(Date) %>% nest() %>% mutate(pred = map(data, extrapol.iv)) %>% unnest(pred) %>% select(-data)
  smpl <- smpl %>% mutate(opt.price = black.price(K.div.F, K, T, impl.vol = vol, r = r))
  smpl %>% group_by(Date) %>% summarise(mfiv = mfiv(opt.price, K, T, r))
}
dta <- fread("/cds_main_all_calls_2011-2024.csv", na.strings = "N/A", stringsAsFactors = FALSE)  %>% mutate(Date = dmy(Date)) %>% arrange(Date) %>% drop_na()
rates <- fread("/bills_mme_yld_2011-2024.csv", na.strings = "N/A", stringsAsFactors = FALSE)  %>% mutate(Date = dmy(Date)) %>% arrange(Date) %>% drop_na()
names(rates) <- c("Date", "r.1.m", "r.3.m")
rates <- rates  %>% mutate(r.1.m = case_when(abs(r.1.m) > 20 ~ NA, TRUE ~ r.1.m)) %>% fill(r.1.m)
rates %>% ggplot(aes(x=Date)) + geom_line(aes(y=r.1.m)) + geom_line(aes(y=r.3.m))
dta <- dta %>% inner_join(rates)
smpl <- dta %>% select("Date", contains("1 month"), "r.1.m") %>% rename(r = r.1.m)
names(smpl) <- str_replace_all( str_replace_all( str_replace_all(names(smpl), "% \\d month Forward", ".F"), "% \\d month Implied Vol", ".vol"), "Main ", "m")
vols <- smpl %>% select("Date", matches("vol")) %>% pivot_longer(!matches("Date"), values_to = "vol") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
strikes <- smpl %>% select("Date", matches("F$")) %>% pivot_longer(!matches("Date"), values_to = "K") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
smpl <- vols %>% inner_join(strikes) %>% inner_join(smpl %>% select("Date", "r"))
out <- gen.mfiv.time.series(smpl, T = 1/12)
mfiv.df <- out %>% rename(main.1.m = mfiv)
smpl <- dta %>% select("Date", contains("3 month"), "r.3.m") %>% rename(r = r.3.m)
names(smpl) <- str_replace_all( str_replace_all( str_replace_all(names(smpl), "% \\d month Forward", ".F"), "% \\d month Implied Vol", ".vol"), "Main ", "m")
vols <- smpl %>% select("Date", matches("vol")) %>% pivot_longer(!matches("Date"), values_to = "vol") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
strikes <- smpl %>% select("Date", matches("F$")) %>% pivot_longer(!matches("Date"), values_to = "K") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
smpl <- vols %>% inner_join(strikes) %>% inner_join(smpl %>% select("Date", "r"))
out <- gen.mfiv.time.series(smpl, T = 3/12)
mfiv.df <- mfiv.df %>% full_join(out %>% rename(main.3.m = mfiv))
dta <- fread("/cds_cross_all_calls_2011-2024.csv", na.strings = "N/A", stringsAsFactors = FALSE)  %>% mutate(Date = dmy(Date)) %>% arrange(Date) %>% drop_na()
dta <- dta %>% inner_join(rates)
smpl <- dta %>% select("Date", contains("1 month"), "r.1.m") %>% rename(r = r.1.m)
names(smpl) <- str_replace_all( str_replace_all( str_replace_all(names(smpl), "% \\d month Forward", ".F"), "% \\d month Implied Vol", ".vol"), "Main ", "m")
vols <- smpl %>% select("Date", matches("vol")) %>% pivot_longer(!matches("Date"), values_to = "vol") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
strikes <- smpl %>% select("Date", matches("F$")) %>% pivot_longer(!matches("Date"), values_to = "K") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
smpl <- vols %>% inner_join(strikes) %>% inner_join(smpl %>% select("Date", "r"))
out <- gen.mfiv.time.series(smpl, T = 1/12)
mfiv.df <- mfiv.df %>% full_join(out %>% rename(cross.1.m = mfiv))
smpl <- dta %>% select("Date", contains("3 month"), "r.3.m") %>% rename(r = r.3.m)
names(smpl) <- str_replace_all( str_replace_all( str_replace_all(names(smpl), "% \\d month Forward", ".F"), "% \\d month Implied Vol", ".vol"), "Main ", "m")
vols <- smpl %>% select("Date", matches("vol")) %>% pivot_longer(!matches("Date"), values_to = "vol") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
strikes <- smpl %>% select("Date", matches("F$")) %>% pivot_longer(!matches("Date"), values_to = "K") %>% mutate(K.div.F = as.numeric(str_extract(name, "\\d{2,3}"))/100) %>% select(-name)
smpl <- vols %>% inner_join(strikes) %>% inner_join(smpl %>% select("Date", "r"))
out <- gen.mfiv.time.series(smpl, T = 3/12)
mfiv.df <- mfiv.df %>% full_join(out %>% rename(cross.3.m = mfiv))
mfiv.df %>% fwrite("/mfiv_2011-24.csv")
################################################################################ VRP calculations
mfiv.sample <- fread("/mfiv_2011-24.csv") %>% mutate(date = ymd(Date)) %>% select(-Date)
dta.vrp <- mfiv.sample %>% inner_join(avg.var.P) %>% fill(main.1.m, main.3.m, cross.1.m, cross.3.m)
dta.vrp <- dta.vrp %>% mutate(VRP.Main.1m = main.1.m - Main.1m.avg.var.P ,
                              VRP.Main.3m = main.3.m - Main.3m.avg.var.P,
                              VRP.Crossover.1m = cross.1.m - Crossover.1m.avg.var.P,
                              VRP.Crossover.3m = cross.3.m - Crossover.3m.avg.var.P) %>% mutate(across(contains("VRP"), function(x) 10000*x))
dta.vrp %>% fwrite("/vrp_mfiv_egarch.csv")
################################################################################
################################################################################
load.data <- function(sample.start = ymd("2016-01-01"), sample.end.excl = ymd("2023-01-01")) {
  PATH.TO.INPUTS <- "/inputs/"
  out <- list()
  basis <- fread(paste0(PATH.TO.INPUTS,'/jpm/cds_bond_basis_IG.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  names(basis) <- str_replace_all( str_remove_all(names(basis), "(CDS - Bond Basis )|( Bullet To Maturity EUR)|(&amp;gt; )|(Average )"), " ", ".")
  out[['basis.ig']] <- basis %>% filter(date >= sample.start, date < sample.end.excl)
  basis.hy <- fread(paste0(PATH.TO.INPUTS,'/jpm/cds_bond_basis_HY_CRSVR.csv'), na.strings = "N/A") %>% select(matches("(date)|(Option)"), -matches("Number"))
  basis.hy <- basis.hy %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  names(basis.hy) <- str_replace_all( str_remove_all(names(basis.hy), "(CDS - Bond Basis )|(Average )|( EUR)"), " ", ".")
  names(basis.hy) <- str_remove_all(names(basis.hy), "\\.PECS$")
  out[['basis.hy']] <- basis.hy %>% filter(date >= sample.start, date < sample.end.excl)
  cds.spreads <- fread(paste0(PATH.TO.INPUTS,'/jpm/itraxx_indices_eod.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  cds.spreads <- cds.spreads %>% select(matches("(Mid)|(date)")) %>% relocate(date)
  names(cds.spreads) <- c("date", "Main", "Crossover")
  out[["cds"]] <- cds.spreads %>% filter(date >= sample.start, date < sample.end.excl)
  brp <- fread(paste0(PATH.TO.INPUTS,"brp.csv")) %>% rename(ebp.lin = eur_index_ebp_linear, ebp.quad = eur_index_ebp_quad) %>% mutate(date = mdy(date))
  brp <-brp %>% select(date, ebp.lin, ebp.quad)
  brp.v2 <- fread(paste0(PATH.TO.INPUTS,"brp_bps.csv")) %>% rename(ebp.v2 = eur_ebp_controls) %>% mutate(date = ymd(date))
  brp <- brp %>% inner_join(brp.v2)
  out[["brp"]] <- brp %>% filter(date >= sample.start, date < sample.end.excl)
  calls <- fread(paste0(PATH.TO.INPUTS,'/jpm/itraxx_atm_calls.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  calls.iv <- calls %>% select(date, matches("Implied"))
  names(calls.iv) <- str_replace_all( str_replace_all(str_remove_all(names(calls.iv), "( 100%)|(month )"), " ", "."), ".Implied.Vol", ".IV")
  out[["calls.iv"]] <- calls.iv %>% filter(date >= sample.start, date < sample.end.excl)
  out[["vrp_e"]] <- fread("/vrp_mfiv_egarch.csv") %>% mutate(date = ymd(date))
  out[["vrp_e"]] <- out[["vrp_e"]] %>% filter(date >= sample.start, date < sample.end.excl)
  blmbg.eod.df <- fread("/blmbg_eod/blmbg_eod.csv") %>% mutate(date = ymd(date)) %>% filter(date >= ymd("2010-01-01"), date < ymd("2023-01-01")) %>% drop_na()
  out[["gs_fin_cond"]] <- blmbg.eod.df %>% filter(date >= sample.start, date < sample.end.excl) %>% select(date, GS.EU.FINCOND)
  gov <- fread(paste0(PATH.TO.INPUTS,'/jpm/govt_yields_europe.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  names(gov) <- str_replace_all(str_remove_all(names(gov), " year Mid Yield"), " ", ".")
  gov <- gov %>% data.frame()
  for(tenor in c("2", "5", "10", "30")) {
    gov[paste0("It.De.",tenor)] <- gov[[paste0("Italy.",tenor)]] - gov[[paste0("Germany.",tenor)]]
  }
  out[["gov"]] <- gov %>% filter(date >= sample.start, date < sample.end.excl)
  blmbg.eod.df <- fread("/blmbg_eod/blmbg_eod.csv") %>% mutate(date = ymd(date)) %>% filter(date >= ymd("2010-01-01"), date < ymd("2023-01-01")) %>% drop_na()
  out[['usd.eur']] <- blmbg.eod.df %>% filter(date >= sample.start, date < sample.end.excl) %>% dplyr::select(date, EUR)
  fx <- fread(paste0(PATH.TO.INPUTS,'/jpm/eur_fx_atm_implied_vols.csv'), na.strings = "N/A") %>% drop_na() %>% mutate(date = dmy(Date)) %>% select(-Date)
  names(fx) <- str_replace_all( str_remove_all(names(fx), "(/USD)|(  At The Money)"), " ", ".")
  out[['eur.vols']] <- fx %>% filter(date >= sample.start, date < sample.end.excl)
  blmbg.eod.df <- fread("/blmbg_eod/blmbg_eod.csv") %>% mutate(date = ymd(date)) %>% filter(date >= ymd("2010-01-01"), date < ymd("2023-01-01")) %>% drop_na()
  out[["equities"]] <- blmbg.eod.df %>% filter(date >= sample.start, date < sample.end.excl) %>% dplyr::select(-GDBR10)
  list2env(out, .GlobalEnv)
} 
################################################################################ Functions for calculating local projections under various model specifications
local.projections.cumul <- function(dta, mps, y, transform = c("level", "log"), surpr.type = "main", sd.shock = TRUE, horizons = seq(1,21), effect.units = "ret.sd", dynamic.confound.ev.rm = FALSE) {
  transform <- match.arg(transform)
  if(transform == "log") {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(log(.x))) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
    dta <- dta %>% mutate(across(!matches("date"), ~log(.x)))
  } else {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(.x)) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
  }
  y.string <- deparse(substitute(y))
  reg.formula.main <- as.formula(paste0("cumul.diff.", y.string, " ~ mps"))
  models <- list()
  se.pval <- list()
  mps.h <- mps
  for(h in horizons) {
    dta.h <- dta %>% arrange(date) %>% transmute(date = date, "cumul.diff.{{y}}" := dplyr::lag(rollapply(data ={{y}}, width = h+1, FUN = function(x) (tail(x,1)-head(x,1)), align = "left", fill = NA),1)) %>% drop_na()
    if(dynamic.confound.ev.rm) { # removing confounding events
      if(h >= 5) mps.h <- mps.h %>% filter(date != ymd("2020-03-25"))
      if(h >= 6) mps.h <- mps.h %>% filter(date != ymd("2020-03-12"))
      if(h >= 7) mps.h <- mps.h %>% filter(date != ymd("2020-04-30"))
      if(h >= 12) mps.h <- mps.h %>% filter(date != ymd("2020-04-22"))
      if(h >= 12) mps.h <- mps.h %>% filter(date != ymd("2012-03-08"), date != ymd("2013-10-02"))
    } else {
      mps.h <- mps
    }
    if(surpr.type == "main") {
      dta.h <- dta.h %>% left_join(mps.h, by = "date") %>% mutate(mps = replace_na(mps.main, 0))
    } else {
      dta.h <- dta.h %>% left_join(mps.h, by = "date") %>% mutate(mps = replace_na(mps.cross, 0))
    }
    models[[paste0("h=",h)]] <- lm(data = dta.h, formula = reg.formula.main)
    if(h==1) {
      se.pval[[paste0("h=",h)]] <- coeftest(models[[paste0("h=",h)]], vcov = NeweyWest(models[[paste0("h=",h)]], lag = ceiling(1.5*h), prewhite = FALSE, adjust = TRUE))
    } else {
      se.pval[[paste0("h=",h)]] <- coeftest(models[[paste0("h=",h)]], vcov = NeweyWest(models[[paste0("h=",h)]], lag = ceiling(1.5*h), prewhite = FALSE, adjust = TRUE))
    }
  }
  se <- unlist(lapply(se.pval, function(x) x[,2]))
  coef <- unlist(lapply(models, function(x) x$coefficients))
  out <- data.frame(coef = coef[str_detect(names(coef), "mps")], se = se[str_detect(names(se), "mps")], expl.var = y.string, h = horizons-1)
  if(sd.shock) {
    if(surpr.type == "main") {
      sd.mps <- sd(mps %>% pull(mps.main)) 
    } else {
      sd.mps <- sd(mps %>% pull(mps.cross))
    }
  } else {
    sd.mps <- 1
  }
  if(effect.units == "ret.sd") {
    out <- out %>% inner_join(std.dev.df) %>% mutate(coef = sd.mps * coef/sd, se  = sd.mps *se/sd)
  } else {
    out <- out %>% inner_join(std.dev.df) %>% mutate(coef = sd.mps*coef, se  = sd.mps*se)
  }
  out
}

local.projections.cumul.ASYMMETRIC <- function(dta, mps, y, transform = c("level", "log"), surpr.type = "main", horizons = seq(1,21), effect.units = "ret.sd", dynamic.confound.ev.rm = FALSE) {
  transform <- match.arg(transform)
  if(transform == "log") {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(log(.x))) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
    dta <- dta %>% mutate(across(!matches("date"), ~log(.x)))
  } else {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(.x)) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
  }
  y.string <- deparse(substitute(y))
  reg.formula <- as.formula(paste0("cumul.diff.", y.string, " ~ mps.hawk + mps.dov"))
  models <- list()
  se.pval <- list()
  mps.h <- mps
  for(h in horizons) {
    dta.h <- dta %>% arrange(date) %>% transmute(date = date, "cumul.diff.{{y}}" := dplyr::lag(rollapply(data ={{y}}, width = h+1, FUN = function(x) (tail(x,1)-head(x,1)), align = "left", fill = NA),1)) %>% drop_na()
    if(dynamic.confound.ev.rm) {
      if(h >= 5) mps.h <- mps.h %>% filter(date != ymd("2020-03-25"))
      if(h >= 6) mps.h <- mps.h %>% filter(date != ymd("2020-03-12"))
      if(h >= 7) mps.h <- mps.h %>% filter(date != ymd("2020-04-30"))
      if(h >= 12) mps.h <- mps.h %>% filter(date != ymd("2020-04-22"))
      if(h >= 12) mps.h <- mps.h %>% filter(date != ymd("2012-03-08"), date != ymd("2013-10-02"))
    } else {
      mps.h <- mps
    }
    if(surpr.type == "main") {
      dta.h <- dta.h %>% left_join(mps.h, by = "date") %>% mutate(mps.main = replace_na(mps.main, 0)) %>% mutate(mps.hawk = pmax(0,mps.main), mps.dov = pmin(0,mps.main))
    } else {
      dta.h <- dta.h %>% left_join(mps.h, by = "date") %>% mutate(mps.cross = replace_na(mps.cross, 0)) %>% mutate(mps.hawk = pmax(0,mps.cross), mps.dov = pmin(0,mps.cross))
    }
    models[[paste0("h=",h)]] <- lm(data = dta.h, formula = reg.formula)
    se.pval[[paste0("h=",h)]] <- coeftest(models[[paste0("h=",h)]], vcov = NeweyWest(models[[paste0("h=",h)]], lag = ceiling(1.5*h), prewhite = FALSE, adjust = TRUE))
  }
  se <- unlist(lapply(se.pval, function(x) x[,2]))
  coef <- unlist(lapply(models, function(x) x$coefficients))
  out.h <- data.frame(coef = coef[str_detect(names(coef), "hawk")], se = se[str_detect(names(se), "hawk")], expl.var = y.string, h = horizons-1, surp = "Hawkish surprise")
  out.d <- data.frame(coef = coef[str_detect(names(coef), "dov")], se = se[str_detect(names(se), "dov")], expl.var = y.string, h = horizons-1, surp = "Dovish surprise")
  out <- bind_rows(list(out.h, out.d))
  mps.sd <- data.frame(sd.mps = c( sd(mps %>% filter(mps.main > 0) %>% pull(mps.main)), sd(mps %>% filter(mps.main < 0) %>% pull(mps.main))), surp = c("Hawkish surprise", "Dovish surprise" ))
  out <- out %>% inner_join(mps.sd)
  if(effect.units == "ret.sd") {
    out <- out %>% inner_join(std.dev.df) %>% mutate(coef = sd.mps*coef/sd, se  = sd.mps*se/sd)
  }
  out
}

local.projections.noncumul <- function(dta, mps, y, transform = c("level", "log"), surpr.type = "main", sd.shock = TRUE, horizons = seq(1,20), digits.table = 4, effect.units = "ret.sd") {
  transform <- match.arg(transform)
  if(transform == "log") {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(log(.x))) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
    dta <- dta %>% mutate(across(!matches("date"), ~log(.x)))
  } else {
    std.dev.df <- dta %>% summarise(across(where(is.numeric), ~sd(diff(.x)) )) %>% pivot_longer(everything()) %>% rename(sd = value, expl.var = name)
  }
  dta <- dta %>% left_join(mps %>% select(date, mps.main), by = "date") %>% mutate(mps.main = replace_na(mps.main, 0))
  Xy <- dta %>% mutate("lag.{{y}}" := dplyr::lag({{y}})) %>% drop_na()
  y.string <- deparse(substitute(y))
  models <- list()
  se.pval <- list()
  models[["reg.main.h=0"]] <- lm(data = Xy %>% mutate("lev.{{y}}":={{y}}), formula = as.formula(paste0("lev.", y.string, " ~ mps.main + lag.", y.string))) # special treatment of h=0
  se.pval[["reg.main.h=0"]] <- coeftest(models[["reg.main.h=0"]], vcov = vcovHC(models[["reg.main.h=0"]], type = "HC1"))
  for( h in horizons) {
    data.h <- Xy %>% mutate("lead.{{y}}" := dplyr::lead({{y}}, h)) %>% drop_na()
    model <- lm(data = data.h, formula = as.formula(paste0("lead.", y.string, " ~ mps.main + lag.", y.string)))
    models[[paste0("reg.main.h=",as.character(h))]] <- model 
    se.pval[[paste0("reg.main.h=",as.character(h))]] <- coeftest(model, vcov = NeweyWest(model, lag = ceiling(1.5*h), prewhite = FALSE, adjust = TRUE))
  }
  se <- unlist(lapply(se.pval, function(x) x[,2]))
  coef <- unlist(lapply(models, function(x) x$coefficients))
  out <- data.frame(coef = coef[str_detect(names(coef), "mps")], se = se[str_detect(names(se), "mps")], expl.var = y.string, h = c(0, horizons))
  if(sd.shock) {
    if(surpr.type == "main") {
      sd.mps <- sd(mps %>% pull(mps.main))
    } else {
      sd.mps <- sd(mps %>% pull(mps.cross))
    }
  } else {
    sd.mps <- 1
  }
  if(effect.units == "ret.sd") {
    out <- out %>% inner_join(std.dev.df) %>% mutate(coef = sd.mps * coef/sd, se  = sd.mps *se/sd)
  } else {
    out <- out %>% inner_join(std.dev.df) %>% mutate(coef = sd.mps*coef, se  = sd.mps*se)
  }
  out
}
################################################################################ Calculating local projections for various series
gen.local.proj <- function(surpr.type = "main", sd.shock = TRUE) {
  results <- list()
  results[["iTraxx CDS ind., log"]] <- bind_rows(list(cds %>% local.projections.cumul(mps, Main, "log", surpr.type, sd.shock = sd.shock),
                                                 cds %>% local.projections.cumul(mps, Crossover, "log", surpr.type, sd.shock = sd.shock)))
  results[[ "CDS-bond basis" ]] <- bind_rows(list(basis.ig %>% rename(IG = IG.Maturity.5y.Basis) %>% local.projections.cumul(mps, IG, "level", surpr.type, sd.shock = sd.shock),
                                                  basis.hy %>% rename(Crossover = iTraxx.Crossover.Basis.Option.Free) %>% local.projections.cumul(mps, Crossover, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "PECS, log" ]] <- bind_rows(list(basis.ig %>% select(date, IG.Maturity.5y.PECS) %>% rename(IG = IG.Maturity.5y.PECS) %>% local.projections.cumul(mps, IG, "log",surpr.type, sd.shock = sd.shock),
                                             basis.hy %>% select(date, iTraxx.Crossover.PECS.Option.Free) %>% rename(Crossover = iTraxx.Crossover.PECS.Option.Free) %>% local.projections.cumul(mps, Crossover, "log", surpr.type, sd.shock = sd.shock)))
  results[[ "Excess bond premium" ]] <- brp %>% local.projections.cumul(mps, ebp.v2, "level", surpr.type, sd.shock = sd.shock)
  results[[ "Swaption ATM vols, Main" ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Main.1.IV) %>% local.projections.cumul(mps, `1.month`, "level", surpr.type, sd.shock = sd.shock),
                                                           calls.iv %>% rename(`3.month` = Main.3.IV) %>% local.projections.cumul(mps, `3.month`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "Swaption ATM vols, Cross" ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Crossover.1.IV) %>% local.projections.cumul(mps, `1.month`, "level", surpr.type, sd.shock = sd.shock),
                                                             calls.iv %>% rename(`3.month` = Crossover.3.IV) %>% local.projections.cumul(mps, `3.month`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "Variance risk premium, Main" ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Main.1m) %>% local.projections.cumul(mps, `1.month`, "level", surpr.type, sd.shock = sd.shock),
                                                               vrp_e %>% rename(`3.month` = VRP.Main.3m) %>% local.projections.cumul(mps, `3.month`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "Variance risk premium, Cross." ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Crossover.1m) %>% local.projections.cumul(mps, `1.month`, "level", surpr.type, sd.shock = sd.shock),
                                                                 vrp_e %>% rename(`3.month` = VRP.Crossover.3m) %>% local.projections.cumul(mps, `3.month`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "GS Europe fin. cond., log" ]] <- gs_fin_cond %>% local.projections.cumul(mps, GS.EU.FINCOND, "log", surpr.type, sd.shock = sd.shock)
  results[[ "Italy-Germany sovereign spread" ]] <- bind_rows(list(gov %>% rename(`2.year` = It.De.2) %>% local.projections.cumul(mps, `2.year`, "level", surpr.type, sd.shock = sd.shock),
                                                                  gov %>% rename(`5.year` = It.De.5) %>% local.projections.cumul(mps, `5.year`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "USD/EUR rate, log" ]] <- usd.eur %>% local.projections.cumul(mps, EUR, "log", surpr.type, sd.shock = sd.shock)
  results[[ "EUR ATM vol." ]] <- bind_rows(list(eur.vols %>% rename(`1.month` = EUR.1.Month.At.The.Money) %>% local.projections.cumul(mps, `1.month`, "level", surpr.type, sd.shock = sd.shock),
                                                eur.vols %>% rename(`6.month` = EUR.6.Month.At.The.Money) %>% local.projections.cumul(mps, `6.month`, "level", surpr.type, sd.shock = sd.shock)))
  results[[ "Equity indices, log" ]] <- bind_rows(list(equities %>% local.projections.cumul(mps, DAX, "log", surpr.type, sd.shock = sd.shock),
                                                       equities %>% rename(Euro_Stoxx_50 = SX5E) %>% local.projections.cumul(mps, Euro_Stoxx_50, "log", surpr.type, sd.shock = sd.shock)))
  results[[ "Implied vol. of equity indices" ]] <- bind_rows(list(equities %>% mutate(DAX = DAX.VIX) %>% local.projections.cumul(mps, DAX, "level", surpr.type, sd.shock = sd.shock),
                                                                  equities %>% rename(Euro_Stoxx_50 = EU.VIX) %>% local.projections.cumul(mps, Euro_Stoxx_50, "level", surpr.type, sd.shock = sd.shock)))
  results
}

gen.local.proj.ASYM <- function(surpr.type = "main") {
  results <- list()
  results[["iTraxx CDS ind., log"]] <- bind_rows(list(cds %>% local.projections.cumul.ASYMMETRIC(mps, Main, "log", surpr.type),
                                                      cds %>% local.projections.cumul.ASYMMETRIC(mps, Crossover, "log", surpr.type)))
  results[[ "CDS-bond basis" ]] <- bind_rows(list(basis.ig %>% rename(IG = IG.Maturity.5y.Basis) %>% local.projections.cumul.ASYMMETRIC(mps, IG, "level", surpr.type),
                                                  basis.hy %>% rename(Crossover = iTraxx.Crossover.Basis.Option.Free) %>% local.projections.cumul.ASYMMETRIC(mps, Crossover, "level", surpr.type)))
  results[[ "PECS, log" ]] <- bind_rows(list(basis.ig %>% select(date, IG.Maturity.5y.PECS) %>% rename(IG = IG.Maturity.5y.PECS) %>% local.projections.cumul.ASYMMETRIC(mps, IG, "log",surpr.type),
                                             basis.hy %>% select(date, iTraxx.Crossover.PECS.Option.Free) %>% rename(Crossover = iTraxx.Crossover.PECS.Option.Free) %>% local.projections.cumul.ASYMMETRIC(mps, Crossover, "log", surpr.type)))
  results[[ "Excess bond premium" ]] <- brp %>% local.projections.cumul.ASYMMETRIC(mps, ebp.v2, "level", surpr.type)
  results[[ "Swaption ATM vols, Main" ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Main.1.IV) %>% local.projections.cumul.ASYMMETRIC(mps, `1.month`, "level", surpr.type),
                                                           calls.iv %>% rename(`3.month` = Main.3.IV) %>% local.projections.cumul.ASYMMETRIC(mps, `3.month`, "level", surpr.type)))
  results[[ "Swaption ATM vols, Cross." ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Crossover.1.IV) %>% local.projections.cumul.ASYMMETRIC(mps, `1.month`, "level", surpr.type),
                                                             calls.iv %>% rename(`3.month` = Crossover.3.IV) %>% local.projections.cumul.ASYMMETRIC(mps, `3.month`, "level", surpr.type)))
  results[[ "Variance risk premium, Main" ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Main.1m) %>% local.projections.cumul.ASYMMETRIC(mps, `1.month`, "level", surpr.type),
                                                               vrp_e %>% rename(`3.month` = VRP.Main.3m) %>% local.projections.cumul.ASYMMETRIC(mps, `3.month`, "level", surpr.type)))
  results[[ "Variance risk premium, Cross." ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Crossover.1m) %>% local.projections.cumul.ASYMMETRIC(mps, `1.month`, "level", surpr.type),
                                                                 vrp_e %>% rename(`3.month` = VRP.Crossover.3m) %>% local.projections.cumul.ASYMMETRIC(mps, `3.month`, "level", surpr.type)))
  results[[ "GS Europe fin. cond., log" ]] <- gs_fin_cond %>% local.projections.cumul.ASYMMETRIC(mps, GS.EU.FINCOND, "log", surpr.type)
  results[[ "Italy-Germany sovereign spread" ]] <- bind_rows(list(gov %>% rename(`2.year` = It.De.2) %>% local.projections.cumul.ASYMMETRIC(mps, `2.year`, "level", surpr.type),
                                                                  gov %>% rename(`5.year` = It.De.5) %>% local.projections.cumul.ASYMMETRIC(mps, `5.year`, "level", surpr.type)))
  results[[ "USD/EUR rate, log" ]] <- usd.eur %>% local.projections.cumul.ASYMMETRIC(mps, EUR, "log", surpr.type)
  results[[ "EUR ATM vol." ]] <- bind_rows(list(eur.vols %>% rename(`1.month` = EUR.1.Month.At.The.Money) %>% local.projections.cumul.ASYMMETRIC(mps, `1.month`, "level", surpr.type),
                                                eur.vols %>% rename(`6.month` = EUR.6.Month.At.The.Money) %>% local.projections.cumul.ASYMMETRIC(mps, `6.month`, "level", surpr.type)))
  results[[ "Equity indices, log" ]] <- bind_rows(list(equities %>% local.projections.cumul.ASYMMETRIC(mps, DAX, "log", surpr.type),
                                                       equities %>% rename(Euro_Stoxx_50 = SX5E) %>% local.projections.cumul.ASYMMETRIC(mps, Euro_Stoxx_50, "log", surpr.type)))
  results[[ "Implied vol. of equity indices" ]] <- bind_rows(list(equities %>% mutate(DAX = DAX.VIX) %>% local.projections.cumul.ASYMMETRIC(mps, DAX, "level", surpr.type),
                                                                  equities %>% rename(Euro_Stoxx_50 = EU.VIX) %>% local.projections.cumul.ASYMMETRIC(mps, Euro_Stoxx_50, "level", surpr.type)))
  results
}

gen.local.proj.alt.spec <- function(surpr.type = "main") {
  results <- list()
  results[["iTraxx CDS ind., log"]] <- bind_rows(list(cds %>% local.projections.noncumul(mps, Main, "log", surpr.type),
                                                      cds %>% local.projections.noncumul(mps, Crossover, "log", surpr.type)))
  results[[ "CDS-bond basis" ]] <- bind_rows(list(basis.ig %>% rename(IG = IG.Maturity.5y.Basis) %>% local.projections.noncumul(mps, IG, "level", surpr.type),
                                                  basis.hy %>% rename(Crossover = iTraxx.Crossover.Basis.Option.Free) %>% local.projections.noncumul(mps, Crossover, "level", surpr.type)))
  results[[ "PECS, log" ]] <- bind_rows(list(basis.ig %>% select(date, IG.Maturity.5y.PECS) %>% rename(IG = IG.Maturity.5y.PECS) %>% local.projections.noncumul(mps, IG, "log",surpr.type),
                                             basis.hy %>% select(date, iTraxx.Crossover.PECS.Option.Free) %>% rename(Crossover = iTraxx.Crossover.PECS.Option.Free) %>% local.projections.noncumul(mps, Crossover, "log", surpr.type)))
  results[[ "Excess bond premium" ]] <- brp %>% local.projections.noncumul(mps, ebp.v2, "level", surpr.type)
  results[[ "Swaption ATM vols, Main" ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Main.1.IV) %>% local.projections.noncumul(mps, `1.month`, "level", surpr.type),
                                                           calls.iv %>% rename(`3.month` = Main.3.IV) %>% local.projections.noncumul(mps, `3.month`, "level", surpr.type)))
  results[[ "Swaption ATM vols, Cross." ]] <- bind_rows(list(calls.iv %>% rename(`1.month` = Crossover.1.IV) %>% local.projections.noncumul(mps, `1.month`, "level", surpr.type),
                                                             calls.iv %>% rename(`3.month` = Crossover.3.IV) %>% local.projections.noncumul(mps, `3.month`, "level", surpr.type)))
  results[[ "Variance risk premium, Main" ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Main.1m) %>% local.projections.noncumul(mps, `1.month`, "level", surpr.type),
                                                               vrp_e %>% rename(`3.month` = VRP.Main.3m) %>% local.projections.noncumul(mps, `3.month`, "level", surpr.type)))
  results[[ "Variance risk premium, Cross." ]] <- bind_rows(list(vrp_e %>% rename(`1.month` = VRP.Crossover.1m) %>% local.projections.noncumul(mps, `1.month`, "level", surpr.type),
                                                                 vrp_e %>% rename(`3.month` = VRP.Crossover.3m) %>% local.projections.noncumul(mps, `3.month`, "level", surpr.type)))
  results[[ "GS Europe fin. cond., log" ]] <- gs_fin_cond %>% local.projections.noncumul(mps, GS.EU.FINCOND, "log", surpr.type)
  results[[ "Italy-Germany sovereign spread" ]] <- bind_rows(list(gov %>% rename(`2.year` = It.De.2) %>% local.projections.noncumul(mps, `2.year`, "level", surpr.type),
                                                                  gov %>% rename(`5.year` = It.De.5) %>% local.projections.noncumul(mps, `5.year`, "level", surpr.type)))
  results[[ "USD/EUR rate, log" ]] <- usd.eur %>% local.projections.noncumul(mps, EUR, "log", surpr.type)
  results[[ "EUR ATM vol." ]] <- bind_rows(list(eur.vols %>% rename(`1.month` = EUR.1.Month.At.The.Money) %>% local.projections.noncumul(mps, `1.month`, "level", surpr.type),
                                                eur.vols %>% rename(`6.month` = EUR.6.Month.At.The.Money) %>% local.projections.noncumul(mps, `6.month`, "level", surpr.type)))
  results[[ "Equity indices, log" ]] <- bind_rows(list(equities %>% local.projections.noncumul(mps, DAX, "log", surpr.type),
                                                       equities %>% rename(Euro_Stoxx_50 = SX5E) %>% local.projections.noncumul(mps, Euro_Stoxx_50, "log", surpr.type)))
  results[[ "Implied vol. of equity indices" ]] <- bind_rows(list(equities %>% mutate(DAX = DAX.VIX) %>% local.projections.noncumul(mps, DAX, "level", surpr.type),
                                                                  equities %>% rename(Euro_Stoxx_50 = EU.VIX) %>% local.projections.noncumul(mps, Euro_Stoxx_50, "level", surpr.type)))
  
  results
}
################################################################################ Function for generating plots of local projections
generate.plots <-function(results, type.desc.str, asym = FALSE, text.size = 20, width = 10, height = 5) {
  custom.theme <- theme(text = element_text(size = text.size), axis.text.x = element_text(size = text.size), axis.text.y.left = element_text(size = text.size), 
                        panel.background = element_rect(fill = "white"), axis.line = element_line(linewidth = 1, color = "black"),
                        axis.ticks.length.x = unit(-0.25, "cm"), axis.ticks.length.y = unit(-0.25, "cm"), axis.ticks = element_line(linewidth = 1, colour = "black"),
                        axis.title.y = element_text(margin = margin(r = text.size)), plot.title = element_text(vjust = 2), plot.subtitle = element_text(vjust = 2),
                        legend.position = "bottom", legend.title = element_blank(), panel.spacing = unit(4, "cm"), strip.background = element_rect(fill = "white"))
  CI.width <- qnorm(0.95)
  for(res in names(results)) {
    fig <- results[[res]] %>% mutate(expl.var = str_replace_all(expl.var, c("\\." = "-", "_" = " ")))  %>% ggplot() + geom_line(aes(x = h, y = coef, color = expl.var)) +
                       geom_ribbon(aes(x=h, ymin = (coef-CI.width*se), ymax = (coef+CI.width*se), fill = expl.var), alpha = 0.25) +
                       geom_hline(yintercept = 0) + custom.theme + ylab(res) + xlab("Cumulative change horizon, h") + scale_x_continuous(breaks = c(1, 5, 10, 15, 20))
    if( length(unique(results[[res]]$expl.var)) == 2 ) {
      fig <- fig + scale_color_manual(values = c("blue", "red")) + scale_fill_manual(values = c("blue", "red"))
    }
    if( length(unique(results[[res]]$expl.var)) == 3 ) {
      fig <- fig + scale_color_manual(values = c("blue", "red", "yellow")) + scale_fill_manual(values = c("blue", "red", "yellow"))
    }
    if(length(unique(results[[res]]$expl.var)) == 1 ) {
      fig <- fig + scale_color_manual(values = c("blue")) + scale_fill_manual(values = c("blue")) + guides(color = "none", fill = "none")
    }
    if(asym) {
      fig <- fig + facet_wrap(vars(surp))
    }
    ggsave(paste0(PATH.TO.OUTPUTS,  type.desc.str, "/", str_remove_all(res, "[/. ,//]"), ".png"), width = width, height = height)
  }
}

################################################################################ Script for generating the graphs of local projections 
PATH.TO.INPUTS <- "/inputs/"
PATH.TO.OUTPUTS <- "/outputs/figs_/"
mps <- read.xlsx(paste0(PATH.TO.INPUTS, "mps.xlsx"), sheet = "mps", detectDates = TRUE)
mps <- mps %>% rename(mps.main = main_res_f, mps.cross = crossover_res_f) %>% select(date, mps.main, mps.cross)
mps <- mps %>% mutate(date = ymd(date)) %>% arrange(date)
load.data(sample.start = ymd("2016-01-01"), sample.end.excl = ymd("2023-01-01"))
results <- gen.local.proj(surpr.type = "main") # CORE results
generate.plots(results, "core")
results.cross <- gen.local.proj(surpr.type = "cross") # Crossover surprise
generate.plots(results.cross, "mps-cross")
results.asym <- gen.local.proj.ASYM(surpr.type = "main") # Asymmetry
generate.plots(results.asym, "asym", asym = TRUE, width = 10, height = 3, text.size = 10)
results.alt.spec <- gen.local.proj.alt.spec(surpr.type = "main")
generate.plots(results.alt.spec, "alt-spec")
load.data(sample.start = ymd("2011-11-01"), sample.end.excl = ymd("2016-01-01")) # Robustness - 2011-11 to 2016
results.pre2016 <- gen.local.proj(surpr.type = "main")
generate.plots(results.pre2016, "pre2016")
############################################################################### LOADING MPS DATA ---- PCA based
mps <- read.xlsx(paste0(PATH.TO.INPUTS, "mps.xlsx"), sheet = "mps", detectDates = TRUE)
mps <- mps %>% rename(mps.main = main_res_pc, mps.cross = crossover_res_pc) %>% select(date, mps.main, mps.cross)
mps <- mps %>% mutate(date = ymd(date)) %>% arrange(date)
load.data(sample.start = ymd("2016-01-01"), sample.end.excl = ymd("2023-01-01"))
results <- gen.local.proj(surpr.type = "main")
generate.plots(results, "PCA-mps-main")

################################################################################ This script contains code used to generate results related to bond issuance
library(tidyverse)
library(lubridate)
library(modelr)
library(stringr)
library(data.table)
library(openxlsx)
library(sandwich)
library(lmtest)
library(texreg)
local.projections.noncumul <- function(dta, y, horizons = seq(1,3), digits.table = 2, file) {
  Xy <- dta %>% mutate("lag.{{y}}" := dplyr::lag({{y}})) %>% drop_na()
  y.string <- deparse(substitute(y))
  models <- list()
  se.pval <- list()
  models[["reg.main.h=0"]] <- lm(data = Xy %>% mutate("lev.{{y}}":={{y}}), formula = as.formula(paste0("lev.", y.string, " ~ mps.main + lag.", y.string))) # special treatment of h=0
  se.pval[["reg.main.h=0"]] <- coeftest(models[["reg.main.h=0"]], vcov = vcovHC(models[["reg.main.h=0"]], type = "HC1"))
  for( h in horizons) {
    data.h <- Xy %>% mutate("lead.{{y}}" := dplyr::lead({{y}}, h)) %>% drop_na()
    model <- lm(data = data.h, formula = as.formula(paste0("lead.", y.string, " ~ mps.main + lag.", y.string)))
    models[[paste0("reg.main.h=",as.character(h))]] <- model 
    se.pval[[paste0("reg.main.h=",as.character(h))]] <- coeftest(model, vcov = NeweyWest(model, lag = ceiling(1.5*h), prewhite = FALSE, adjust = TRUE))
  }
  texreg(models, override.se = lapply(se.pval, function(x) x[,2]), override.pvalues = lapply(se.pval, function(x) x[,4]),
         digits = digits.table, include.rsquared = TRUE, include.adjrs = TRUE, include.nobs = FALSE, stars = c(0.01, 0.05, 0.10))
}
PATH.TO.INPUTS <- "/inputs/"
PATH.TO.OUTPUTS <- "/outputs/"
mps <- read.xlsx(paste0(PATH.TO.INPUTS, "mps.xlsx"), sheet = "mps", detectDates = TRUE)
mps <- mps %>% rename(mps.main = main_res_f, mps.cross = crossover_res_f) %>% select(date, mps.main, mps.cross) 
mps <- mps %>% mutate(date = ymd(date)) %>% arrange(date)
mps.completed <- mps %>% mutate(Month = month(date), Year = year(date))
mps.completed <- mps.completed %>% complete(Month = seq(1,12), Year = seq(2011,2022)) %>% arrange(Year, Month) 
mps.agg <- mps.completed %>% group_by(Year, Month) %>% summarise(mps.main = sum(mps.main, na.rm = TRUE),
                                                                 mps.cross = sum(mps.cross, na.rm = TRUE),
                                                                 Year = Year[1], Month = Month[1]) %>% arrange(Year, Month)
mps.agg <- mps.agg %>% ungroup()
net.bond.issuance <- read.csv(paste0(PATH.TO.INPUTS,"/issuance/net_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) %>% rename(net.issuance.non.fin = issuanceTotal) # Net non-financial corporate bond issuance in euro zone - SOURCE: ECB - column K
net.bond.issuance <- net.bond.issuance %>% mutate(date = mdy(Date)) %>% mutate(Month = month(date), Year = year(date)) %>% arrange(date) %>% select(-Date)
dta <-  net.bond.issuance %>% full_join(mps.agg, by = c("Year" = "Year", "Month" = "Month")) %>% arrange(Year, Month) %>% mutate(date = ym(paste(Year, Month))) %>% filter(date >= ym("2016/01")) %>% drop_na()
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(net.issuance.non.fin, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/net_issuance.doc'))
dta <-  net.bond.issuance %>% full_join(mps.agg, by = c("Year" = "Year", "Month" = "Month")) %>% arrange(Year, Month) %>% mutate(date = ym(paste(Year, Month))) %>% drop_na() # Alternative sample
dta %>% filter(date >= ym("2011-11"), date < ym("2016-01")) %>% local.projections.noncumul(net.issuance.non.fin, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/net_issuance_pre2016.doc'))
bond.Issuance <- read.csv(paste0(PATH.TO.INPUTS,"issuance/gross_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) # GROSS non-financial corporate bond issuance in euro zone
bond.Issuance <- bond.Issuance %>% mutate(year.month = paste(year, month, sep = '-')) %>% mutate(date = ym(year.month)) %>% arrange(date)
bond.Issuance <- bond.Issuance %>% rename(issuance.TOTAL = amtEuro1, issuance.IG = amtEuroIG1, issuance.NEW = amtEuroNew1, Year = year, Month = month) %>% mutate(issuance.HY = issuance.TOTAL - issuance.IG)
dta <- inner_join(mps.agg, bond.Issuance) %>% arrange(Year, Month)
dta <- dta %>% select(-contains("amt"), -time)
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.IG, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_IG_issuance.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.TOTAL, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_TOTAL_issuance.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.HY, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_HY_issuance.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.NEW, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_NEW_issuance.doc'))
dta %>% filter(date >= ym("2011-11"), date < ym("2016-01")) %>% local.projections.noncumul(issuance.IG, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_IG_issuance.doc'))
dta %>% filter(date >= ym("2011-11"), date < ym("2016-01")) %>% local.projections.noncumul(issuance.TOTAL, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_TOTAL_issuance.doc'))
dta %>% filter(date >= ym("2011-11"), date < ym("2016-01")) %>% local.projections.noncumul(issuance.HY, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_HY_issuance.doc'))
################################################################################ PCA based surprises
mps <- read.xlsx(paste0(PATH.TO.INPUTS, "mps.xlsx"), sheet = "mps", detectDates = TRUE)
mps <- mps %>% rename(mps.main = main_res_pc, mps.cross = crossover_res_pc) %>% select(date, mps.main, mps.cross) 
mps <- mps %>% mutate(date = ymd(date)) %>% arrange(date)
mps.completed <- mps %>% mutate(Month = month(date), Year = year(date))
mps.completed <- mps.completed %>% complete(Month = seq(1,12), Year = seq(2011,2022)) %>% arrange(Year, Month) 
mps.agg <- mps.completed %>% group_by(Year, Month) %>% summarise(mps.main = sum(mps.main, na.rm = TRUE),
                                                                 mps.cross = sum(mps.cross, na.rm = TRUE),
                                                                 Year = Year[1], Month = Month[1]) %>% arrange(Year, Month)
mps.agg <- mps.agg %>% ungroup()
net.bond.issuance <- read.csv(paste0(PATH.TO.INPUTS,"/issuance/net_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) %>% rename(net.issuance.non.fin = issuanceTotal) # Net non-financial corporate bond issuance in euro zone - SOURCE: ECB - column K
net.bond.issuance <- net.bond.issuance %>% mutate(date = mdy(Date)) %>% mutate(Month = month(date), Year = year(date)) %>% arrange(date) %>% select(-Date)
dta <- net.bond.issuance %>% full_join(mps.agg, by = c("Year" = "Year", "Month" = "Month")) %>% arrange(Year, Month) %>% mutate(date = ym(paste(Year, Month))) %>% filter(date >= ym("2016/01")) %>% drop_na()
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(net.issuance.non.fin, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/net_issuance_PCA.doc'))
bond.Issuance <- read.csv(paste0(PATH.TO.INPUTS,"issuance/gross_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) #GROSS non-financial corporate bond issuance in euro zone
bond.Issuance <- bond.Issuance %>% mutate(year.month = paste(year, month, sep = '-')) %>% mutate(date = ym(year.month)) %>% arrange(date)
bond.Issuance <- bond.Issuance %>% rename(issuance.TOTAL = amtEuro1, issuance.IG = amtEuroIG1, issuance.NEW = amtEuroNew1, Year = year, Month = month) %>% mutate(issuance.HY = issuance.TOTAL - issuance.IG)
dta <- inner_join(mps.agg, bond.Issuance) %>% arrange(Year, Month)
dta <- dta %>% select(-contains("amt"), -time)
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.TOTAL, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_TOTAL_issuance_PCA.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.IG, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_IG_issuance_PCA.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.HY, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_HY_issuance_PCA.doc'))
################################################################################ Crossover surprises
mps <- read.xlsx(paste0(PATH.TO.INPUTS, "mps.xlsx"), sheet = "mps", detectDates = TRUE)
mps <- mps %>% rename(mps.main = main_res_pc, mps.cross = crossover_res_pc) %>% select(date, mps.main, mps.cross) %>% mutate(mps.main = mps.cross)# OVERRIDING the name
mps <- mps %>% mutate(date = ymd(date)) %>% arrange(date)
mps.completed <- mps %>% mutate(Month = month(date), Year = year(date))
mps.completed <- mps.completed %>% complete(Month = seq(1,12), Year = seq(2011,2022)) %>% arrange(Year, Month) 
mps.agg <- mps.completed %>% group_by(Year, Month) %>% summarise(mps.main = sum(mps.main, na.rm = TRUE),
                                                                 mps.cross = sum(mps.cross, na.rm = TRUE),
                                                                 Year = Year[1], Month = Month[1]) %>% arrange(Year, Month)
mps.agg <- mps.agg %>% ungroup()
net.bond.issuance <- read.csv(paste0(PATH.TO.INPUTS,"/issuance/net_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) %>% rename(net.issuance.non.fin = issuanceTotal) # Net non-financial corporate bond issuance in euro zone - SOURCE: ECB - column K
net.bond.issuance <- net.bond.issuance %>% mutate(date = mdy(Date)) %>% mutate(Month = month(date), Year = year(date)) %>% arrange(date) %>% select(-Date)
dta <- net.bond.issuance %>% full_join(mps.agg, by = c("Year" = "Year", "Month" = "Month")) %>% arrange(Year, Month) %>% mutate(date = ym(paste(Year, Month))) %>% filter(date >= ym("2016/01")) %>% drop_na()
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(net.issuance.non.fin, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/net_issuance.doc'))
bond.Issuance <- read.csv(paste0(PATH.TO.INPUTS,"issuance/gross_issuance.csv"), na.strings = "N/A", stringsAsFactors = FALSE) # GROSS non-financial corporate bond issuance in euro zone
bond.Issuance <- bond.Issuance %>% mutate(year.month = paste(year, month, sep = '-')) %>% mutate(date = ym(year.month)) %>% arrange(date)
bond.Issuance <- bond.Issuance %>% rename(issuance.TOTAL = amtEuro1, issuance.IG = amtEuroIG1, issuance.NEW = amtEuroNew1, Year = year, Month = month) %>% mutate(issuance.HY = issuance.TOTAL - issuance.IG)
dta <- inner_join(mps.agg, bond.Issuance) %>% arrange(Year, Month)
dta <- dta %>% select(-contains("amt"), -time)
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.TOTAL, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_TOTAL_issuance_CO.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.IG, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_IG_issuance_CO.doc'))
dta %>% filter(date >= ym("2016-01")) %>% local.projections.noncumul(issuance.HY, horizons = seq(1,3), file = paste0(PATH.TO.OUTPUTS, 'tbls/gross_HY_issuance_CO.doc'))