## cross-sectional policy rule regressions using SPF data

library(readxl)
library(tidyr)
library(dplyr)
library(purrr) # map_int
library(readr)
require(readxl)
library(lubridate)
library(sandwich)
library(plm)
library(ggplot2)
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(gridExtra)

## read SPF individual forecasts from Excel, pivot to long format
tb <- read_excel(path = "data/SPFmicrodata.xlsx",
                 sheet = "TBILL",
                 na = "#N/A",
                 col_type = "numeric") %>%
    select(-(TBILLA:TBILLD), -INDUSTRY) %>%
    pivot_longer(TBILL1:TBILL6,
                 names_to = "h",
                 names_prefix = "TBILL",
                 values_to = "tb") %>%
    na.omit
cpi <- read_excel(path = "data/SPFmicrodata.xlsx",
                  sheet = "CPI",
                  na = "#N/A",
                  col_type = "numeric") %>%
    select(-(CPIA:CPIC), -INDUSTRY) %>%
    pivot_longer(CPI1:CPI6,
                 names_to = "h",
                 names_prefix = "CPI",
                 values_to = "cpi") %>%
    na.omit
ur <- read_excel(path = "data/SPFmicrodata.xlsx",
                 sheet = "UNEMP",
                 na = "#N/A",
                 col_type = "numeric") %>%
    select(-(UNEMPA:UNEMPD), -INDUSTRY) %>%
    pivot_longer(UNEMP1:UNEMP6,
                 names_to = "h",
                 names_prefix = "UNEMP",
                 values_to = "ur") %>%
    na.omit

## join forecast data
spf <- tb %>%
    left_join(cpi) %>%
    left_join(ur) %>%
    filter(!is.na(cpi), !is.na(ur)) %>%
    rename(id = ID) %>%
    mutate(h = as.numeric(h),
           yq = YEAR*10 + QUARTER,
           date = make_date(YEAR, QUARTER*3, 1)) %>%
    select(yq, date, id, h, tb, cpi, ur) %>%
    ## add lagged interest rate forecast
    group_by(yq, date, id) %>%
    arrange(h) %>%
    mutate(ltb = dplyr::lag(tb)) %>%
    ungroup %>%
    filter(!is.na(cpi), h > 1)  ## exclude backcasts -- h=2 is nowcast

spf %>% select(yq, id) %>% distinct %>% group_by(yq) %>% summarise(n=n()) %>% pull(n) %>% summary

## estimate baseline policy rules
d1 <- spf %>%
    group_by(yq, date) %>%
    nest %>%
    mutate(n = map_int(data, function(df) nrow(na.omit(df))),
           model_fe = map(data, function(df) plm(tb ~ cpi + ur, data = df, index=c("id", "h"), model = "within")),
           beta = map_dbl(model_fe, function(m) m$coef[1]),
           gamma = map_dbl(model_fe, function(m) m$coef[2]),
           V = map(model_fe, function(m) vcovHC(m)),
           beta_se = map_dbl(V, ~sqrt(diag(.x))[1]),
           gamma_se = map_dbl(V, ~sqrt(diag(.x))[2]),
           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

## compare to BCFF estimates
load("output/bluechip_rule_regressions.RData")
bc1 <- data %>% rename(gamma_bc=gamma_fe)

##################################################
## estimate inertial policy rules
d2 <- spf %>%
    group_by(yq, date) %>%
    nest %>%
    mutate(n = map_int(data, function(df) nrow(na.omit(df))),
           model = map(data, function(df) plm(tb ~ ltb + cpi + ur, data = df, index=c("id", "h"), model = "within")),
           rho = map_dbl(model, function(m) m$coef["ltb"]),
           beta = map_dbl(model, function(m) m$coef["cpi"]),
           gamma = map_dbl(model, function(m) m$coef["ur"]),
           V = map(model, function(m) vcovHC(m)), # forecaster clustering
           rho_se = map_dbl(V, ~sqrt(diag(.x))["ltb"]),
           beta_se = map_dbl(V, ~sqrt(diag(.x))["cpi"]),
           gamma_se = map_dbl(V, ~sqrt(diag(.x))["ur"]),
           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


## compare to BCFF estimates
load("output/bluechip_rule_inertial.RData")
bc2 <- data %>% rename(gamma_bc = gamma)

## Figure B.4

d1 <- d1 %>%
    mutate(gamma_spf = -0.5*gamma,
           gamma_spf_lb = -0.5*gamma_lb,
           gamma_spf_ub = -0.5*gamma_ub)
d2 <- d2 %>%
    mutate(gamma_spf = -0.5*gamma,
           gamma_spf_lb = -0.5*gamma_lb,
           gamma_spf_ub = -0.5*gamma_ub)

cols <- c("steelblue", "black")
p1 <- ggplot(d1) +
    geom_ribbon(aes(x=date, ymax=gamma_spf_ub, ymin=gamma_spf_lb), fill="lightgray") +
    geom_line(aes(x=date, y=gamma_spf, colour="SPF"), linewidth=0.75) +
    geom_line(data=bc1, mapping=aes(x=date, y=gamma_bc, colour="BCFF"), linewidth=0.75) +
    geom_line(aes(x=date, y=gamma_spf, colour="SPF"), linewidth=0.75) +
    geom_hline(yintercept=0) +
    scale_colour_manual(name="",values=cols) +
    ggtitle("Baseline rule") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=11)) +
    theme(legend.position = "inside",
          legend.position.inside = c(0.2, 0.85),
          legend.direction="horizontal",
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "black", fill= "white"),
          plot.margin = unit(c(0,0.2,0,0), "cm"))
p2 <- ggplot(d2) +
    geom_ribbon(aes(x=date, ymax=gamma_spf_ub, ymin=gamma_spf_lb), fill="lightgray") +
    geom_line(aes(x=date, y=gamma_spf, colour="SPF"), linewidth=0.75) +
    geom_line(data=bc2, mapping=aes(x=date, y=gamma_bc, colour="BCFF"), linewidth=0.75) +
    geom_line(aes(x=date, y=gamma_spf, colour="SPF"), linewidth=0.75) +
    geom_hline(yintercept=0) +
    scale_colour_manual(name="",values=cols) +
    ggtitle("Inertial rule") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=11),
          legend.position = "none",
          plot.margin = unit(c(0,0.2,0,-0.25), "cm"))
dev.new()
grid.arrange(p1, p2)
g <- arrangeGrob(p1, p2)
ggsave("figures/spf.pdf", g, width=6, height=5)
