## cross-sectional policy rule regressions using SEP forecasts

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

sep <- read_xlsx("data/SEP.xlsx", col_names = c("date", "id", "year", "rgdp", "unemp", "pce", "pcex", "ff"), skip = 1) %>%
    filter(year != "LR") %>% ## leave out LR forecasts
    mutate(date = as.Date(date),
           year = as.numeric(year),
           h = year - year(date)) %>%
    group_by(date, id) %>%
    arrange(h) %>%
    mutate(rgdpcum = cumprod(1 + rgdp/100)) %>%  # cumulative real GDP growth %>%
    filter(h<=1) %>%
    ungroup

## output gap forecasts
## (1) real-time GDP from ALFRED
gdp <- read.delim("data/alfred/GDPC1.txt", na.strings=".") %>% as_tibble %>%
    rename(date = observation_date) %>%
    pivot_longer(-date, names_to = "vin", values_to = "gdp", names_prefix = "GDPC1_") %>%
    mutate(date = as.Date(date),
           vin = as.Date(vin, format="%Y%m%d")) %>%
    na.omit %>%
    filter(quarter(date) == 4) %>% # only keep most recent end-of-year GDP level for each vintage
    group_by(vin) %>% arrange(date) %>% slice_tail(n=1) %>% ungroup %>%
    filter(year(vin) >= 2011) %>%
    mutate(year = year(date)) %>%
    select(-date)

## determine units for each vintage
gdpunits <- read.fwf("data/alfred/GDPC1_UNITS.txt", width = c(72, 13, 13), col.names = c("units", "start", "end")) %>% as_tibble %>%
    mutate(across(everything(), stringr::str_trim),  # trim whitespaces
           units = readr::parse_number(units),   # simplify units description (only number)
           across(c(start, end), as.Date)) %>%
    na.omit # removes vintages from 2018 onwards
gdpvintages <- gdp %>% select(vin) %>% distinct %>% arrange(vin) %>%
    mutate(units = gdpunits$units[findInterval(vin, gdpunits$start)]) %>%
    rename(gdpvin = vin, gdpunits = units)

## (2) real-time potential GDP projections from ALFRED
pot <- read.delim("data/alfred/GDPPOT.txt", na.strings=".") %>% as_tibble %>%
    rename(date = observation_date) %>%
    pivot_longer(-date, names_to = "vin", values_to = "pot", names_prefix = "GDPPOT_") %>%
    mutate(date = as.Date(date),
           vin = as.Date(vin, format="%Y%m%d")) %>%
    na.omit %>%
    arrange(vin) %>%
    mutate(year = year(date)) %>%
    filter(quarter(date) == 4,  # only keep end-of-year projections for each vintage
           year(vin) >= 2011,
           year >= 2010) %>%
    select(-date)
## determine units for each vintage
potunits <- read.fwf("data/alfred/GDPPOT_UNITS.txt", width = c(72, 13, 13), col.names = c("units", "start", "end")) %>%
    as_tibble %>%
    mutate(across(everything(), stringr::str_trim),
           units = readr::parse_number(units),   # simplify units description (only number)
           across(c(start, end), as.Date))
potvintages <- pot %>% select(vin) %>% distinct %>%
    mutate(units = potunits$units[findInterval(vin, potunits$start)]) %>%
    rename(potvin = vin, potunits = units)

## for each survey, find most recent end-of-year GDP (from latest real-time vintage)
fomc <- sep %>% select(date) %>% distinct %>%
    mutate(gdpvin = gdpvintages$gdpvin[findInterval(date, gdpvintages$gdpvin)])
## Jan-2012 SEP - pretend newest GDP release (with 2011-Q4 numbers) is available, released three days later
if (fomc$gdpvin[1] == make_date(2011, 12, 22))
    fomc$gdpvin[1] <- make_date(2012, 1, 27)
fomc <- left_join(fomc, gdp %>% rename(gdpvin = vin, gdpyear = year))
stopifnot(all(fomc$gdpyear == year(fomc$date) - 1))
## for each FOMC, find most recent vintage of potential GDP projections
fomc <- fomc %>%
    mutate(potvin = potvintages$potvin[findInterval(date, potvintages$potvin)]) %>%
    left_join(potvintages) %>%
    left_join(gdpvintages)
## deal with different units: when most recent potential GDP does not have same units as real GDP, use first potential GDP vintage with same units (pretend it was available in real time)
for (j in 1:nrow(fomc)) {
    if (fomc$gdpunits[j] != fomc$potunits[j]) {
        potvinnew <- potvintages %>% filter(potunits == fomc$gdpunits[j]) %>% pull(potvin) %>% min
        fomc$potvin[j] <- potvinnew
        fomc$potunits[j] <- fomc$gdpunits[j]
    }
}
## add potential GDP projections
fomc <- fomc %>%
    expand_grid(h = 0:5) %>%
    mutate(year = year(date) + h) %>%
    left_join(pot %>% rename(potvin = vin))

## merge current GDP and pot. GDP forecasts into surveys
sep <- sep %>%
    inner_join(fomc %>% select(date, h, gdp, pot)) %>%
    # construct output gap
    mutate(gdp = gdp*rgdpcum,  # gdp level forecast
           gap = 100*(gdp - pot)/pot) %>%
    arrange(date, h, id)     # nice sorting

## estimation
d <- sep %>%
    group_by(date) %>%
    nest %>%
    mutate(
        model = map(data, function(df) plm(ff ~ gap + pce, data=df, index=c("id", "h"), model = "within")),
        gamma_fe = map_dbl(model, function(m) m$coef["gap"]),
        beta_fe = map_dbl(model, function(m) m$coef["pce"]),
        V = map(model, function(m) vcovHC(m)),
        gamma_se = map_dbl(V, ~sqrt(diag(.x))["gap"]),
        beta_se = map_dbl(V, ~sqrt(diag(.x))["pce"]),
        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


load("output/bluechip_rule_regressions.RData")
bcest <- data %>%
    select(date, gamma_fe, beta_fe)

## Figure B.5
cols <- c("FOMC" = "black", "BCFF" = "steelblue")
p1 <- ggplot() +
    geom_ribbon(data = d, mapping = aes(x=date, ymax=gamma_ub, ymin=gamma_lb), fill="gray92", alpha=0.7) +
    geom_line(data = d, mapping = aes(x=date, y=gamma_fe, colour="FOMC"), linewidth=1) +
    geom_line(data = bcest, mapping = aes(x=date, y=gamma_fe, colour="BCFF"), linewidth=1) +
    geom_hline(yintercept=0) +
    geom_vline(xintercept = make_date(2015, 12, 17), linetype=2) +
    xlim(min(d$date), max(d$date)) +
    scale_colour_manual(name = "Estimate", values = cols) +
    scale_y_continuous(expand = c(0,0)) +
    ggtitle("Output gap coefficient") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0.1,0.1,-0.1,0), "cm"),
          legend.position = "inside",
          legend.position.inside = c(0.1, 0.8),
          legend.key.height = unit(.4, 'cm'),
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "gray", fill= "white"),
          legend.text = element_text(size = 8))
p2 <- ggplot() +
    geom_ribbon(data = d, mapping = aes(x=date, ymax=beta_ub, ymin=beta_lb), fill="gray92", alpha=0.7) +
    geom_line(data = d, mapping = aes(x=date, y=beta_fe, colour="FOMC"), linewidth=1) +
    geom_line(data = bcest, mapping = aes(x=date, y=beta_fe, colour="BCFF"), linewidth=1) +
    xlim(min(d$date), max(d$date)) +
    coord_cartesian(ylim=c(-1.5, 3)) +
    geom_vline(xintercept = make_date(2015, 12, 17), linetype=2) +
    scale_colour_manual(name = "Estimate", values = cols) +
    scale_y_continuous(expand = c(0,0)) +
    geom_hline(yintercept=0) +
    ggtitle("Inflation coefficient") + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0.1,0.1,-0.1,0), "cm"),
          legend.position = "none")
grid.arrange(p1, p2)
g <- arrangeGrob(p1, p2)
ggsave("figures/sep.pdf", g, width=6, height=4)
