## policy rule panel regressions using Bluechip data

library(dplyr)
library(tidyr) # nest
library(purrr) # map
library(sandwich)
library(plm)

load("data/bc_May23.RData")

data <- bc %>%
    select(ff, gap, cpi, date, id, h) %>%
    group_by(date) %>%
    nest %>%
    mutate(
        ## Pooled OLS
        model = map(data, function(df) plm(ff ~ gap + cpi, data=df, index=c("id", "h"), model = "pooling")),
        alpha_ols = map_dbl(model, function(m) m$coef["(Intercept)"]),
        gamma_ols = map_dbl(model, function(m) m$coef["gap"]),
        beta_ols = map_dbl(model, function(m) m$coef["cpi"]),
        ## forecaster FE
        model_fe = map(data, function(df) plm(ff ~ gap + cpi, data=df, index=c("id", "h"), model = "within")),
        gamma_fe = map_dbl(model_fe, function(m) m$coef["gap"]),
        beta_fe = map_dbl(model_fe, function(m) m$coef["cpi"]),
        SE = map(model_fe, function(m) {
            ## two-way clustering
            v <- diag(vcovDC(m))
            if (any(v < 0))
                v <- diag(vcovHC(m)) # if N/A, one-way clustering
            return(sqrt(v))}),
        gamma_se = map_dbl(SE, function(s) s["gap"]),
        beta_se = map_dbl(SE, function(s) s["cpi"]),
        gamma_lb = gamma_fe - 1.96*gamma_se,
        gamma_ub = gamma_fe + 1.96*gamma_se,
        beta_lb = beta_fe - 1.96*beta_se,
        beta_ub = beta_fe + 1.96*beta_se
    ) %>%
    ungroup

## save estimates
data <- data %>%
    select(date,
           gamma_ols, gamma_fe, gamma_lb, gamma_ub,
           beta_ols, beta_fe, beta_lb, beta_ub,
           alpha_ols)

write.csv(data, file="output/bluechip_rule_regressions.csv", row.names=FALSE, quote=FALSE)
save(data, file="output/bluechip_rule_regressions.RData")

