## cross-sectional policy rule regressions using Bluechip data
## -> rule specification with credit spreads

library(dplyr)
library(tidyr)
library(purrr)
library(modelr)
library(lubridate)
library(ggplot2)
library(sandwich)
library(plm)
library(gridExtra)

load("data/bc_May23.RData")

d <- bc %>%
    filter(year(date) >= 2001) %>%
    mutate(spread = baa - y10) %>%
    select(ff, cpi, gap, spread, date, id, h) %>%
    group_by(date) %>%
    nest %>%
     mutate(
        ## baseline estimates
        model1 = map(data, function(df) plm(ff ~ gap + cpi, data=df, index=c("id", "h"), model = "within")),
        gamma1 = map_dbl(model1, function(m) m$coef["gap"]),
        SE1 = map(model1, function(m) {
            v <- diag(vcovSCC(m, maxlag=2))
            return(sqrt(v))}),
        gamma1_se = map_dbl(SE1, function(s) s["gap"]),
        gamma1_lb = gamma1 - 1.96*gamma1_se,
        gamma1_ub = gamma1 + 1.96*gamma1_se,
        ## with credit spreads
        model2 = map(data, function(df) plm(ff ~ gap + cpi + spread, data=df, index=c("id", "h"), model = "within")),
        gamma2 = map_dbl(model2, function(m) m$coef["gap"]),
        delta = map_dbl(model2, function(m) m$coef["spread"]),
        SE2 = map(model2, function(m) {
            v <- diag(vcovSCC(m, maxlag=2))
            return(sqrt(v))}),
        gamma2_se = map_dbl(SE2, function(s) s["gap"]),
        delta_se = map_dbl(SE2, function(s) s["spread"]),
        gamma2_lb = gamma2 - 1.96*gamma2_se,
        gamma2_ub = gamma2 + 1.96*gamma2_se,
        delta_lb = delta - 1.96*delta_se,
        delta_ub = delta + 1.96*delta_se
     ) %>%
    ungroup

cols <- c("Credit spreads" = "black",  "Baseline" = "steelblue")
p1 <- ggplot(d) +
    geom_ribbon(aes(x=date, ymax=gamma1_ub, ymin=gamma1_lb), fill="lightgray") +
    geom_line(aes(x=date, y=gamma1, colour="Baseline"), linewidth=1) +
    geom_line(aes(x=date, y=gamma2, colour="Credit spreads"), linewidth=1) +
    geom_hline(yintercept=0) +
    scale_colour_manual(name="Estimate",values=cols) +
    ggtitle("Output gap coefficient") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=11)) +
    theme(legend.position = "inside",
          legend.position.inside = c(0.2, 0.15),
          legend.direction="horizontal",
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "black", fill= "white"))
p2 <- ggplot(d) +
    geom_ribbon(aes(x=date, ymax=delta_ub, ymin=delta_lb), fill="lightgray") +
    geom_line(aes(x=date, y=delta, colour="Credit spreads"), linewidth=1) +
    scale_colour_manual(name="Estimate",values=cols) +
    geom_hline(yintercept=0) +
    ggtitle("Spread coefficient") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=11)) +
    theme(legend.position = "none")
grid.arrange(p1, p2)

## compare to baseline estimates
d %>% select(gamma1, gamma2) %>% cor
