## key figures for the paper

library(dplyr)
library(ggplot2)
library(gridExtra)
library(scales)
library(lubridate)
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(latex2exp)

cols = c("Perceived" = "steelblue", "Historical" = "black")

### Figure 1: regression estimates of simple perceived rule
## estimated in bluechip_rules_regressions.R
load("output/bluechip_rule_regressions.RData")
data <- data %>%
    rename(gamma_hat = gamma_fe,
           beta_hat = beta_fe)
breaks <- seq(min(data$date), max(data$date), "years")
labels <- ifelse(year(breaks)%%10==0, year(breaks), "")
actual <- new.env()
load("output/actual_rule_simple.RData", env = actual)
data <- inner_join(data, actual$data)  # join with actual/historical rule
p1 <- ggplot(data) +
    geom_ribbon(aes(x=date, ymax=gamma_actual_ub, ymin=gamma_actual_lb), fill="gray92") +
    geom_ribbon(aes(x=date, ymax=gamma_ub, ymin=gamma_lb), fill="lightblue") +
    geom_line(aes(x=date, y=gamma_actual, colour="Historical"), linewidth=0.75) +
    geom_line(aes(x=date, y=gamma_hat, colour="Perceived"), linewidth=0.75) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(-0.75, 1.65)) +
    ggtitle(TeX("Output gap coefficient, $\\hat{\\gamma}$")) + ylab("") + xlab("") +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0), breaks = seq(-0.5, 2, 0.5)) +
    scale_colour_manual(name = "Estimate", values = cols) +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0.5,0,-0.25), "cm"),
          legend.position = "inside",
          legend.position.inside = c(0.2, 0.88),
          legend.key.height = unit(.4, 'cm'),
          legend.direction="horizontal",
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "gray", fill= "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
p2 <- ggplot(data) +
    geom_ribbon(aes(x=date, ymax=beta_actual_ub, ymin=beta_actual_lb),
                fill="gray92") +
    geom_ribbon(aes(x=date, ymax=beta_ub, ymin=beta_lb), fill="lightblue") +
    geom_line(aes(x=date, y=beta_actual, colour="Historical"), linewidth=0.75) +
    geom_line(aes(x=date, y=beta_hat, colour="Perceived"), linewidth=0.75) +
    scale_colour_manual(name = "Estimate", values = cols) +
    coord_cartesian(ylim=c(-1.25, 1.75)) +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0), breaks = seq(-1, 2, 0.5)) +
    geom_hline(yintercept=0) +
    ggtitle(TeX("Inflation coefficient, $\\hat{\\beta}$")) + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0.5,0,0), "cm"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
dev.new()
grid.arrange(p1, p2)
g <- arrangeGrob(p1, p2)
ggsave("figures/bluechip_rule_simple.pdf", g, width=6, height=5)

## correlation between perceived and actual gamma
data %>%
    group_by(period = ifelse(date < as.Date("2008-12-01"), "pre-ZLB", "after")) %>%
    summarize(mean(gamma_hat), mean(gamma_actual), cor(gamma_actual, gamma_hat))

### Figure 2: inertial rule -- bluechip_rules_inertial.R
load("output/bluechip_rule_inertial.RData")
data <- data %>%
    rename(gamma_hat = gamma,
           beta_hat = beta,
           rho_hat = rho)
actual <- new.env()
load("output/actual_rule_inertial.RData", env = actual)
data <- inner_join(data, actual$data) # join with actual/historical rule
p1 <- ggplot(data) +
    geom_ribbon(aes(x=date, ymax=gamma_actual_ub, ymin=gamma_actual_lb), fill="gray92") +
    geom_ribbon(aes(x=date, ymax=gamma_ub, ymin=gamma_lb), fill="lightblue") +
    geom_line(aes(x=date, y=gamma_actual, colour="Historical"), linewidth=0.75) +
    geom_line(aes(x=date, y=gamma_hat, colour="Perceived"), linewidth=0.75) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(-0.5, 1.0)) +
    ggtitle(TeX("Output gap coefficient, $\\hat{\\gamma}$")) + ylab("") + xlab("") +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0)) +
    scale_colour_manual(name = "Estimate", values = cols) +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0.5,0,0), "cm"),
          legend.position = "inside",
          legend.position.inside = c(0.8, 0.87),
          legend.key.height = unit(.4, 'cm'),
          legend.direction="horizontal",
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "gray", fill= "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
p2 <- ggplot(data) +
    geom_ribbon(aes(x=date, ymax=beta_actual_ub, ymin=beta_actual_lb),
                fill="gray92") +
    geom_ribbon(aes(x=date, ymax=beta_ub, ymin=beta_lb), fill="lightblue") +
    geom_line(aes(x=date, y=beta_actual, colour="Historical"), linewidth=0.75) +
    geom_line(aes(x=date, y=beta_hat, colour="Perceived"), linewidth=0.75) +
    coord_cartesian(ylim=c(-0.6, 1.2)) +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0)) +
    scale_colour_manual(name = "Estimate", values = cols) +
    geom_hline(yintercept=0) +
    ggtitle(TeX("Inflation coefficient,  $\\hat{\\beta}$")) + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0.5,0,0), "cm"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
p3 <- ggplot(data) +
    geom_ribbon(aes(x=date, ymax=rho_actual_ub, ymin=rho_actual_lb),
                fill="gray92") +
    geom_ribbon(aes(x=date, ymax=rho_ub, ymin=rho_lb), fill="lightblue") +
    geom_line(aes(x=date, y=rho_actual, colour="Historical"), linewidth=0.75) +
    geom_line(aes(x=date, y=rho_hat, colour="Perceived"), linewidth=0.75) +
    geom_hline(yintercept=0) +
    coord_cartesian(ylim=c(-0.5, 2.5)) +
    ggtitle(TeX("Lagged interest rate coefficient, $\\hat{\\rho}$")) + ylab("") + xlab("") +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0)) +
    scale_colour_manual(name = "Estimate", values = cols) +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0.5,0,0), "cm"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
dev.new()
grid.arrange(p1, p2, p3)
g <- arrangeGrob(p1, p2, p3)
ggsave("figures/bluechip_rule_inertial.pdf", g, width=6, height=7.5)

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

## Appendix Figure B.1
## -> SSM, OLS and FE  (all simple rules, not inertial)
load("output/ssm_ols.RData")
d <- data %>%
    rename(gamma_ssm_ols = gamma,
           beta_ssm_ols = beta,
           beta_ssm_ols_lb = beta_lb,
           beta_ssm_ols_ub = beta_ub,
           gamma_ssm_ols_lb = gamma_lb,
           gamma_ssm_ols_ub = gamma_ub)
load("output/ssm_fe.RData")
data <- data %>%
    rename(gamma_ssm_fe = gamma,
           beta_ssm_fe = beta,
           beta_ssm_fe_lb = beta_lb,
           beta_ssm_fe_ub = beta_ub,
           gamma_ssm_fe_lb = gamma_lb,
           gamma_ssm_fe_ub = gamma_ub)
d <- inner_join(d, data)
load("output/bluechip_rule_regressions.RData")
d <- inner_join(d, data %>% select(date, gamma_ols, beta_ols, gamma_fe, beta_fe))
breaks <- seq(min(data$date), max(data$date), "years")
labels <- ifelse(year(breaks)%%10==0, year(breaks), "")

cols <- c("SSM OLS" = "black", "SSM FE" = "steelblue", "OLS" = "red", "FE" = "green")
p1 <- ggplot(d) +
    geom_ribbon(aes(x=date, ymax=gamma_ssm_ols_ub, ymin=gamma_ssm_ols_lb), fill="lightblue") +
    geom_ribbon(aes(x=date, ymax=gamma_ssm_fe_ub, ymin=gamma_ssm_fe_lb), fill="lightblue") +
    geom_line(aes(x=date, y=gamma_ols, colour="OLS"), linewidth=0.5) +
    geom_line(aes(x=date, y=gamma_fe, colour="FE"), linewidth=0.5) +
    geom_line(aes(x=date, y=gamma_ssm_ols, colour="SSM OLS"), linewidth=0.5) +
    geom_line(aes(x=date, y=gamma_ssm_fe, colour="SSM FE"), linewidth=0.5) +
    geom_hline(yintercept=0) +
    scale_colour_manual(name = "Estimate", values = cols) +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0)) +
    ggtitle(TeX("Output gap coefficient $\\hat{\\gamma}$")) + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0,0,-0.17), "cm"),
          legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
p2 <- ggplot(d) +
    geom_ribbon(aes(x=date, ymax=beta_ssm_ols_ub, ymin=beta_ssm_ols_lb), fill="lightblue") +
    geom_ribbon(aes(x=date, ymax=beta_ssm_fe_ub, ymin=beta_ssm_fe_lb), fill="lightblue") +
    geom_line(aes(x=date, y=beta_ols, colour="OLS"), linewidth=0.5) +
    geom_line(aes(x=date, y=beta_fe, colour="FE"), linewidth=0.5) +
    geom_line(aes(x=date, y=beta_ssm_ols, colour="SSM OLS"), linewidth=0.5) +
    geom_line(aes(x=date, y=beta_ssm_fe, colour="SSM FE"), linewidth=0.5) +
    geom_hline(yintercept=0) +
    scale_colour_manual(name = "Estimate", values = cols) +
    scale_x_date(expand = c(0,0), breaks = breaks, labels = labels) +
    scale_y_continuous(expand = c(0,0)) +
    ggtitle(TeX("Inflation coefficient $\\hat{\\beta}$")) + ylab("") + xlab("") +
    theme(plot.title = element_text(size=10),
          plot.margin = unit(c(0,0,0,-0.4), "cm"),
          legend.position = "inside",
          legend.position.inside = c(0.72, 0.9),
          legend.key.height = unit(.4, 'cm'),
          legend.direction="horizontal",
          legend.title = element_blank(),
          legend.spacing.y = unit(0, "mm"),
          legend.background = element_rect(colour = "gray", fill= "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
g <- arrangeGrob(p1, p2)
ggsave("figures/ssm.pdf", g, width=6.5, height=5.5)
dev.new()
grid.arrange(p1, p2)

## correlation between SSM and panel estimates
with(d, cor(gamma_ssm_ols, gamma_ols))
with(d, cor(gamma_ssm_fe, gamma_fe))
