## How does gamma change in response to business cycle and FOMC surprises?

library(dplyr)
library(readr)
library(readxl)
library(lubridate)
library(sandwich)
library(zoo)
library(latex2exp)
library(stargazer)
library(ggplot2)
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(gridExtra)

load("output/bluechip_rule_regressions.RData")
d <- data %>%
    select(date, gamma_fe) # simple/baseline gamma
load("output/bluechip_rule_inertial.RData")
d2 <- data %>%
    rename(gamma_sm = gamma) %>% # for SMoothing
    select(date, gamma_sm) # inertial gamma
d <- inner_join(d, d2) %>%
     mutate(ym = year(date)*100 + month(date))

## output gap
ogap <- read_csv("data/output_gap.csv", col_names = c("date", "ogap"), skip = 1)
d <- left_join(d, ogap) %>%
    mutate(ogap = na.locf(ogap))

## FOMC surprises
d <- mutate(d, ym = year(date)*100 + month(date))
fomc <- read_xlsx("data/fomc_bauer_swanson.xlsx", sheet="FOMC Announcements", na="NA") %>%
    rename(date = Date, mps = MPS) %>%
    ## rename(date = Date, mps = MPS_ORT) %>%  ## use orthogonalized surprises
    mutate(date = as.Date(date)) %>%
    select(date, mps) %>%
    mutate(ym = year(date)*100 + month(date)) %>%
    group_by(ym) %>%
    summarise(mps = sum(mps))

d <- left_join(d, fomc) %>%
    mutate(mps = ifelse(is.na(mps), 0, mps),
           D = ogap < median(ogap), # recession indicator
           i1 = D*1, i2 = 1-i1,
           z1 = mps*i1, z2 = mps*i2)

##################################################
## local projections

get_lp <- function(dat, N = 12, confint = 1.96) {
    res <- tibble(j = 0:N, R2 = NA, b1 = NA, se1 = NA, b2 = NA, se2 = NA)
    dat <- mutate(dat, lgamma = dplyr::lag(gamma, 1))
    for (h in 0:N) {
        dat$y <- dplyr::lead(dat$gamma, h)
        mod <- lm(y ~ z1 + z2 + D + lgamma, dat)
        res$R2[h+1] <- summary(mod)$r.squared
        nw_lags <- ceiling(h*1.5)
        V <- NeweyWest(mod, lag = nw_lags, prewhite=FALSE)
        res$b1[h+1] <- mod$coef[2]
        res$se1[h+1] <- sqrt(diag(V))[2]
        res$b2[h+1] <- mod$coef[3]
        res$se2[h+1] <- sqrt(diag(V))[3]
    }
    res %>% mutate(t1 = b1/se1, t2 = b2/se2,
                   lb1 = b1 - confint*se1, lb2 = b2 - confint*se2,
                   ub1 = b1 + confint*se1, ub2 = b2 + confint*se2)
}

res1 <- get_lp(d %>% mutate(gamma = gamma_fe))
res2 <- get_lp(d %>% mutate(gamma = gamma_sm))
ymin <- min(res1$lb1, res1$lb2, res2$lb1, res2$lb2)
ymax <- max(res1$ub1, res1$ub2, res2$ub1, res2$ub2)
theme_update(plot.margin = unit(c(5,5,2,0), "pt"),
             text = element_text(size=10),
             plot.title = element_text(size=10))
p1 <- ggplot(res1) +
    geom_ribbon(aes(x=j, ymax=ub2, ymin=lb2), fill="steelblue", alpha=.5) +
    geom_line(aes(x=j, y=lb2)) +
    geom_line(aes(x=j, y=ub2)) +
    geom_line(aes(x=j, y=b2), linewidth=1.5) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(-0.5, ymax)) +
    ylab("") +
    xlab("")+
    ggtitle(TeX("Response of baseline $\\hat{\\gamma}$ in strong economy")) +
    scale_x_continuous(expand = c(0,0)) +
    theme(plot.margin = unit(c(0.1,0.2,0,-0.2), "cm"))
p2 <- ggplot(res1) +
    geom_ribbon(aes(x=j, ymax=ub1, ymin=lb1), fill="steelblue", alpha=.5) +
    geom_line(aes(x=j, y=lb1)) +
    geom_line(aes(x=j, y=ub1)) +
    geom_line(aes(x=j, y=b1), linewidth=1.5) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(ymin, 0.75)) +
    ylab("") +
    ggtitle(TeX("Response of baseline $\\hat{\\gamma}$ in weak economy")) +
    xlab("Horizon (months)") +
    scale_x_continuous(expand = c(0,0)) +
    theme(plot.margin = unit(c(0.1,0.2,0,0), "cm"))
p3 <- ggplot(res2) +
    geom_ribbon(aes(x=j, ymax=ub2, ymin=lb2), fill="steelblue", alpha=.5) +
    geom_line(aes(x=j, y=lb2)) +
    geom_line(aes(x=j, y=ub2)) +
    geom_line(aes(x=j, y=b2), linewidth=1.5) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(-0.5, ymax)) +
    ylab("") +
    xlab("")+
    ggtitle(TeX("Response of inertial $\\hat{\\gamma}$ in strong economy")) +
    scale_x_continuous(expand = c(0,0)) +
    theme(plot.margin = unit(c(0.1,0.2,0,-0.2), "cm"))
p4 <- ggplot(res2) +
    geom_ribbon(aes(x=j, ymax=ub1, ymin=lb1), fill="steelblue", alpha=.5) +
    geom_line(aes(x=j, y=lb1)) +
    geom_line(aes(x=j, y=ub1)) +
    geom_line(aes(x=j, y=b1), linewidth=1.5) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(ymin, 0.75)) +
    ylab("") +
    ggtitle(TeX("Response of inertial $\\hat{\\gamma}$ in weak economy")) +
    xlab("Horizon (months)") +
    scale_x_continuous(expand = c(0,0)) +
    theme(plot.margin = unit(c(0.1,0.2,0,0), "cm"))
dev.new()
grid.arrange(p1, p3, p2, p4)
g <- arrangeGrob(p1, p3, p2, p4)
ggsave("figures/gamma_lp.pdf", g, width=6, height=4.2)

cat("Sample period:", range(d$ym), "\n")
cat("Maximum absolute values:\n")
cat("baseline gamma, strong economy:", max(res1$b2), "\n")
cat("baseline gamma, weak economy:", min(res1$b1), "\n")
cat("inertial gamma, strong economy:", max(res2$b2), "\n")
cat("inertial gamma, weak economy:", min(res2$b1), "\n")

## regressions with interaction effect---testing for difference
mods <- list()
d <- d %>% mutate(gamma = gamma_fe, lgamma = dplyr::lag(gamma, 1))
i <- 1
for (h in c(3, 6, 9, 12)) {
    d$y <- dplyr::lead(d$gamma, h)
    mod <- lm(y ~ mps + mps:D + D + lgamma, d)
    mods[[i]] <- mod
    b <- mod$coef
    V <- NeweyWest(mod, lag = ceiling(h*1.5), prewhite=FALSE)
    SEs <- sqrt(diag(V))
    tstats <- b / SEs
    all.equal(summary(mod)$r.squared, res1$R2[h+1])
    all.equal(mod$coef[2], res1$b2[h+1], check.attributes=FALSE)
    all.equal(mod$coef[2] + mod$coef[5], res1$b1[h+1], check.attributes=FALSE)
    all.equal(SEs[2], res1$se2[h+1], check.attributes=FALSE)
    i <- i + 1
}
d <- d %>% mutate(gamma = gamma_sm, lgamma = dplyr::lag(gamma, 1))
for (h in c(3, 6, 9, 12)) {
    d$y <- dplyr::lead(d$gamma, h)
    mod <- lm(y ~ mps + mps:D + D + lgamma, d)
    mods[[i]] <- mod
    b <- mod$coef
    V <- NeweyWest(mod, lag = ceiling(h*1.5), prewhite=FALSE)
    SEs <- sqrt(diag(V))
    tstats <- b / SEs
    all.equal(summary(mod)$r.squared, res2$R2[h+1])
    all.equal(mod$coef[2], res2$b2[h+1], check.attributes=FALSE)
    all.equal(mod$coef[2] + mod$coef[5], res2$b1[h+1], check.attributes=FALSE)
    all.equal(SEs[2], res2$se2[h+1], check.attributes=FALSE)
    i <- i + 1
}
rob_se <- mapply(function(mod, h)
    sqrt(diag(NeweyWest(mod, lag = ceiling(h*1.5), prewhite = FALSE))),
    mods, c(3,6,9,12), SIMPLIFY=FALSE)
all.equal(rob_se[[1]][2], res1$se2[4], check.attributes=FALSE)
stargazer(mods, digits=2, header=FALSE, type="text", se=rob_se, report="vc*t", no.space=TRUE, omit.stat=c("LL", "ser", "f", "adj.rsq"))
