## How is gamma related to yield curve, cyclical variables, monetary policy cycle?
## VIX and uncertainty?

library(dplyr)
library(lubridate)
library(sandwich)
library(readr)
library(stargazer)
library(zoo)

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
dat <- inner_join(d, d2) %>%
    mutate(date = dplyr::lag(date), # lead gamma by one month - publication lag
           ym = year(date)*100 + month(date)) %>%
    na.omit

## yields
gsw <- read_csv("data/feds200628.csv", skip=9) %>%
    rename(date = Date) %>%
    filter(year(date) >= 1984) %>%
    select(date, SVENY01:SVENY10) %>%
    na.omit %>%
    mutate(ym = year(date)*100 + month(date)) %>%
    group_by(ym) %>%
    summarize(across(SVENY01:SVENY10, mean)) # monthly averages
Y <- gsw %>% select(SVENY01:SVENY10) %>% data.matrix
W <- eigen(cov(Y))$vectors[,1:3]
W[,1] <- W[,1]/sum(W[,1])
gsw[c("L", "S", "C")] <- Y %*% W
gsw <- gsw %>%
    mutate(lagS = dplyr::lag(S, 12)) %>% # 12m-lagged slope for regressions
    select(ym, SVENY01, L, S, lagS)
dat <- left_join(dat, gsw)

## unemployment rate
ur <- read_csv("data/UNRATE.csv") %>%
    rename(date = DATE) %>%
    mutate(ym = year(date)*100 + month(date)) %>%
    select(ym, UNRATE) %>%
    mutate(lagUNRATE = dplyr::lag(UNRATE, 12)) # 12m-lagged unempl. rate

dat <- left_join(dat, ur)

## VIX
vix <- read_csv("data/VIXCLS_monthly.csv", na=".") %>%
    rename(date = DATE, vix = VIXCLS) %>% na.omit
dat <- left_join(dat, vix)
vxo <- read_csv("data/VXOCLS_monthly.csv", na=".") %>%
    rename(date = DATE, vxo = VXOCLS) %>% na.omit
dat <- left_join(dat, vxo)
dat <- mutate(dat, vix = ifelse(is.na(vix), vxo, vix)) %>% select(-vxo)

## Fed easing and tightening cycles, NBER recessions
easing <- data.frame(
    xstart = as.Date(c("1989-06-01", "1990-07-13", "1995-07-06", "1998-09-29", "2001-01-03", "2007-09-18", "2019-08-01")),
    xend = as.Date(c("1990-01-31", "1992-09-04", "1996-01-31", "1998-11-17", "2003-06-25", "2008-12-16", "2020-06-30")))
tightening <- data.frame(
    xstart = as.Date(c("1988-04-01", "1994-02-04", "1997-02-01", "1999-06-30", "2004-06-30", "2015-12-17", "2022-03-17")),
    xend = as.Date(c("1989-04-30", "1995-02-01", "1997-04-30", "2000-05-16", "2006-06-29", "2018-12-20", "2023-07-27")))
zlbperiods <- data.frame(xstart = as.Date(c("2008-12-16", "2020-03-16")),
                         xend = c(as.Date(c("2015-12-17", "2022-03-16"))))
make_dummy <- function(datemat, ddates) {
    x <- rep(0, length(ddates))
    for (i in 1:nrow(datemat)) {
        if (max(ddates) >= datemat[i,1]) {
            ind1 <- min(which(ddates >= datemat[i,1]))
            ind2 <- max(which(ddates <= datemat[i,2]))
            x[ind1:ind2] <- 1
        }
    }
    x
}

dat <- dat %>%
    mutate(easing = make_dummy(easing, date),
           tightening = make_dummy(tightening, date),
           zlb = make_dummy(zlbperiods, date))

##################################################
## CHOOSE gamma simple or inertial HERE

## Panel A
dat <- mutate(dat, gamma = gamma_fe)   # simple/baseline rule

## Panel B
##dat <- mutate(dat, gamma = gamma_sm)   # inertial rule

##################################################

mods <- list(
    lm(gamma ~ lagS, dat),
    lm(gamma ~ tightening, dat),
    lm(gamma ~ UNRATE, dat),
    lm(gamma ~ zlb, dat),
    lm(gamma ~ vix, dat))
rob_se <- lapply(mods, function(m) sqrt(diag(NeweyWest(m, lag = 12, prewhite=FALSE))))
stargazer(mods, digits=2, header=FALSE, se=rob_se, title="Determinants of gamma", type="text",  no.space=TRUE, omit.stat=c("LL", "ser", "f", "adj.rsq"))

