## perceived rule with lagged interest rate

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

load("data/bc_May23.RData")

## lagged interest rate -- previous horizon for each survey/forecasters
bc <- bc %>%
    group_by(date, id) %>%
    arrange(h) %>%
    mutate(lff = dplyr::lag(ff)) %>%
    ungroup

## add lagged fed funds rate for h=0
ff <- read_csv("data/FEDFUNDS.csv") %>%
    mutate(date = as.Date(DATE, format="%m/%d/%Y")) %>%
    rename(ff = FEDFUNDS) %>%
    mutate(yq = year(date)*10 + quarter(date)) %>%
    group_by(yq) %>%
    summarise(ffobs = mean(ff)) %>%
    filter(yq > 19840)

bc <- bc %>%
    mutate(ldate = date %m-% months(3),
           yq = year(ldate)*10 + quarter(ldate)) %>%
    left_join(ff) %>%
    mutate(lff = ifelse(h == 0, ffobs, lff)) %>%
    select(-ldate, -yq, -ffobs)

data <- bc %>%
    select(ff, lff, gap, cpi, date, id, h) %>%
    group_by(date) %>%
    nest %>%
    mutate(
        ## forecaster FE
        model = map(data, function(df) plm(ff ~ lff + gap + cpi, data = df, index=c("id", "h"), model = "within")),
        rho = map_dbl(model, function(m) m$coef["lff"]),
        gamma = map_dbl(model, function(m) m$coef["gap"]),
        beta = map_dbl(model, function(m) m$coef["cpi"]),
        SE = map(model, function(m) {
            v <- diag(vcovDC(m))
            if (any(v < 0))
                v <- diag(vcovHC(m))
            return(sqrt(v)) }),
        rho_se = map_dbl(SE, function(s) s["lff"]),
        gamma_se = map_dbl(SE, function(s) s["gap"]),
        beta_se = map_dbl(SE, function(s) s["cpi"]),
        rho_lb = rho - 1.96*rho_se,
        rho_ub = rho + 1.96*rho_se,
        gamma_lb = gamma - 1.96*gamma_se,
        gamma_ub = gamma + 1.96*gamma_se,
        beta_lb = beta - 1.96*beta_se,
        beta_ub = beta + 1.96*beta_se) %>%
    ungroup

data <- data %>%
    select(date,
           rho, rho_lb, rho_ub,
           gamma, gamma_lb, gamma_ub,
           beta, beta_lb, beta_ub)

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

## mean of rho-hat pre/post 2000
data %>% group_by(year(date) < 2000) %>% summarize(mean(rho))

