################################################################################
###
### AUTHOR:       Carlo Favero 
### DATE:         Milan, December 2026
### DESCRIPTION:  This script replicates Pflueger-Viceira on US data smpl 1999-2009 
### OUTPUT:           
################################################################################

## -----------------------------------------------------------------------------
rm(list=ls())                                                                   # Clear the environment 

## -----------------------------------------------------------------------------
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))                     # Set the appropriate environment to work in every machine

## -----------------------------------------------------------------------------
listofpackages <- c("readxl","dplyr","lubridate","zoo","stringr",     # Install and load the required packages
                    "ggplot2","scales","sandwich","lmtest","stargazer")   
source("preparePackages.R")
preparePackages(listofpackages)


# ============================================================
# 1) Import data from Excel 
# ============================================================

df <- read_xlsx("monthly_2025.xlsx", sheet = "Sheet1")
#options(zoo.yearmon.format = "%Y-%m")
df <- df %>%
  mutate(
    date = as.yearmon(date)   # monthly index: "1998-01"
  ) %>%
  arrange(date)

# check
head(df$date)
#options(zoo.yearmon.format = "%Y-%m")
#head(df$date)

# ============================================================
# 2) Data construction 
# ============================================================
# Helper:  lag/lead functions
L <- function(x, k = 1) dplyr::lag(x, k)
F <- function(x, k = 1) dplyr::lead(x, k)

df <- df %>%
  mutate(
    infl_q  = log(cpi / L(cpi, 3)),
    infl_qa = 4 * log(cpi / L(cpi, 3)),
    infl_a  = log(cpi / L(cpi, 12)),
    
    b_10y   = sveny10 - tipsy10,
    
    # returns ( quarterly returns annualized )
    ret10y_n = (0.25 * L(sveny10, 3) - 9.75 * (sveny10 - L(sveny10, 3))) * 4,
    ret10y_r = (0.25 * L(tipsy10, 3) - 9.75 * (tipsy10 - L(tipsy10, 3))) * 4,
    ret10y_b = ret10y_n - ret10y_r,
    
    ret3m = tb3m_cca
  )

# ============================================================
# 3) Estimation: auxiliary model for the real rate (HAC)
# ============================================================

df <- df %>%
  mutate(
    lhs = ret3m - F(infl_qa, 3),
    x1  = ret3m,
    x2  = L(ret3m, 3) - infl_qa,
    x3  = infl_a
  )

df_aux <- df %>%
  filter(date >= as.yearmon("1998-01"),
         date <= as.yearmon("2021-10")) %>%
  filter(complete.cases(lhs, x1, x2, x3))

fit_3m <- lm(lhs ~ x1 + x2 + x3, data = df_aux)
summary(fit_3m)

# HAC (Newey-West). Choose lag ~ 12 for monthly as a sensible default.
# 
vcov_hac <- sandwich::NeweyWest(fit_3m, lag = 12, prewhite = FALSE, adjust = TRUE)
eq_3m_hac <- lmtest::coeftest(fit_3m, vcov. = vcov_hac)
print(eq_3m_hac)

# Residuals and implied real rate 
df_aux <- df_aux %>%
  mutate(resids01 = residuals(fit_3m),
         rr3m = lhs - resids01)

# generate excess returns
df_aux <- df_aux %>%
    mutate(
    b_3m = ret3m - rr3m,
    
    # excess returns
    exret_n = ret10y_n - L(ret3m, 3),
    exret_r = ret10y_r - L(rr3m, 3),
    exret_b = ret10y_b - L(b_3m, 3),
    
    # spreads and their lags 
    sp_10y_3m_n = sveny10 - ret3m,
    sp_10y_3m_r = tipsy10 - rr3m,
    sp_10y_3m_b = b_10y - b_3m,
    sp_n_l3 = L(sp_10y_3m_n, 3),
    sp_r_l3 = L(sp_10y_3m_r, 3),
    sp_b_l3 = L(sp_10y_3m_b, 3)
  )

# ============================================================
# 4) Risk premia regressions 
# exret_* = c + beta * spread(-3)
# rp_* = fitted = exret - resid
# ============================================================

# 1) run regressions (df_aux is fine)
fit_rp_n <- lm(exret_n ~ sp_n_l3, data = df_aux, na.action = na.omit)
fit_rp_r <- lm(exret_r ~ sp_r_l3, data = df_aux, na.action = na.omit)
fit_rp_b <- lm(exret_b ~ sp_b_l3, data = df_aux, na.action = na.omit)

summary(fit_rp_n)
summary(fit_rp_r)
summary(fit_rp_r)
# 2) create empty columns, then fill only where the model was estimated
df_aux <- df_aux %>%
  mutate(rp_n = NA_real_, rp_r = NA_real_, rp_b = NA_real_)

df_aux$rp_n[as.integer(rownames(model.frame(fit_rp_n)))] <- fitted(fit_rp_n)
df_aux$rp_r[as.integer(rownames(model.frame(fit_rp_r)))] <- fitted(fit_rp_r)
df_aux$rp_b[as.integer(rownames(model.frame(fit_rp_b)))] <- fitted(fit_rp_b)

# ============================================================
# 6) “Super nice” ggplot figures (your three figures)
# ============================================================

theme_paper <- function() {
  theme_minimal(base_size = 12) +
    theme(
      panel.grid.minor = element_blank(),
      plot.title = element_text(face = "bold"),
      legend.title = element_blank(),
      legend.position = "right"
    )
}

# ---- Figure 1: 3-month rates (rr3m, realized real rate ret3m-infl_qa(+3), nominal ret3m) ----
fig_rr <- df_aux %>%
  filter(date >= as.yearmon("1999-04-01"),
         date <= as.yearmon("2021-10-01")) %>%
  transmute(
    date,
    `Predicted real rate` = rr3m,
    `Realized real rate`  = ret3m - F(infl_qa, 3),
    `Nominal rate`        = ret3m
  ) %>%
  tidyr::pivot_longer(-date, names_to = "series", values_to = "value")


p_rr <- ggplot(
  fig_rr,
  aes(x = date, y = value, linetype = series, color = series)
) +
  geom_hline(yintercept = 0, linewidth = 0.4, alpha = 0.6) +
  geom_line(linewidth = 1.1) +
  scale_x_yearmon(n = 5, format = "%Y") +
  scale_color_manual(
    values = c(
      "Nominal rate"        = "blue",
      "Predicted real rate" = "green4",
      "Realized real rate"  = "red"
    )
  ) +
  labs(
    title = "3-month rates",
    x = NULL,
    y = NULL
  ) +
  theme_paper()

p_rr


# ---- Figure 3: Excess returns  (rp_r, rp_n, rp_b) ----
fig_er <- df_aux %>%
  filter(date >= as.yearmon("1999-04")) %>%
  transmute(
    date,
    `TIPS`         = exret_r,
    `Nominal`      = exret_n,
    `Breakeven RP` = exret_b
  ) %>%
  tidyr::pivot_longer(-date, names_to = "series", values_to = "value")

p_er <- ggplot(fig_er, aes(x = date, y = value, linetype = series, color = series)) +
  geom_hline(yintercept = 0, linewidth = 0.4, alpha = 0.6) +
  geom_line(linewidth = 1.0) +
  scale_x_yearmon(n = 5, format = "%Y") +
  scale_color_manual(
    values = c(
      "Nominal"      = "blue",
      "TIPS"         = "red",
      "Breakeven RP" = "green4"
    )
  ) +
  labs(title = "Bond Excess Returns", x = NULL, y = NULL) +
  theme_paper()

p_er

# ---- Figure 2: Bond risk premia (rp_r, rp_n, rp_b) ----
fig_rp <- df_aux %>%
  filter(date >= as.yearmon("1999-04")) %>%
  transmute(
    date,
    `TIPS`         = rp_r,
    `Nominal`      = rp_n,
    `Breakeven RP` = rp_b
  ) %>%
  tidyr::pivot_longer(-date, names_to = "series", values_to = "value")

p_rp <- ggplot(fig_rp, aes(x = date, y = value, linetype = series, color = series)) +
  geom_hline(yintercept = 0, linewidth = 0.4, alpha = 0.6) +
  geom_line(linewidth = 1.0) +
  scale_x_yearmon(n = 5, format = "%Y") +
  scale_color_manual(
    values = c(
      "Nominal"      = "blue",
      "TIPS"         = "red",
      "Breakeven RP" = "green4"
    )
  ) +
  labs(title = "Bond Risk Premia", x = NULL, y = NULL) +
  theme_paper()

p_rp
# ============================================================
# 7) Descriptive statistics & correlations 
# ============================================================

desc_stats <- function(x) {
  x <- x[is.finite(x)]
  tibble(
    mean = mean(x),
    sd   = sd(x),
    min  = min(x),
    p25  = quantile(x, 0.25),
    med  = median(x),
    p75  = quantile(x, 0.75),
    max  = max(x),
    n    = length(x)
  )
}

grp_ret <- df_aux %>% select(infl_q, ret3m, rr3m)
grp_sp  <- df_aux %>% select(sp_10y_3m_n, sp_10y_3m_r, sp_10y_3m_b)
grp_ex  <- df_aux %>% select(exret_n, exret_r, exret_b)

Table_1 <- grp_ret %>% summarise(across(everything(), desc_stats)) %>% tidyr::pivot_longer(everything())
Table_2 <- cor(grp_ret, use = "pairwise.complete.obs")

Table_3 <- grp_sp  %>% summarise(across(everything(), desc_stats)) %>% tidyr::pivot_longer(everything())
Table_4 <- cor(grp_sp, use = "pairwise.complete.obs")

Table_5 <- grp_ex  %>% summarise(across(everything(), desc_stats)) %>% tidyr::pivot_longer(everything())
Table_6 <- cor(grp_ex, use = "pairwise.complete.obs")

# ============================================================
# 5) Extracting Inflation Expectations 
# EViews:
# smpl 1999:1 2021m10
# eq_ar1rp: rp_b = c + phi * rp_b(-1)
# then for i=1..73: rolling forecast path rp_b_f1 over 120 months,
# average expected future rp over 120 months: sum_{j=0..119} rp_b_f1(+j) / 120
# rp_av assigned at each origin; then infl_e10y = b_10y - rp_av
# ============================================================

df_ie <- df_aux %>%
  filter(date >= as.yearmon("1999-04"),
         date <= as.yearmon("2021-10")) %>%
  mutate(rp_b_l1 = L(rp_b, 1)) %>%
  filter(complete.cases(rp_b, rp_b_l1))

fit_ar1 <- lm(rp_b ~ rp_b_l1, data = df_ie)
summary(fit_ar1)
c0 <- coef(fit_ar1)[1]
phi <- coef(fit_ar1)[2]

# Forecast-horizon = 120 months (10 years)
H <- 120

# loop assigns rp_av for dates in a given range

origins <- seq(from = as.yearmon("2015-01"),
               to   = as.yearmon("2021-10"),
               by   = 1/12)
# Function: average of expected future rp over next H months given rp_t
avg_future_ar1 <- function(rp_t, c0, phi, H) {
  # E[rp_{t+h}] = mu + phi^h (rp_t - mu), mu = c0/(1-phi)
  if (is.na(rp_t)) return(NA_real_)
  if (abs(1 - phi) < 1e-8) {
    # nearly unit root: fallback to linear expectation
    exp_path <- rp_t + c0 * (0:(H-1))
  } else {
    mu <- c0 / (1 - phi)
    h <- 0:(H-1)
    exp_path <- mu + (phi^h) * (rp_t - mu)
  }
  mean(exp_path)
}

# Compute rp_av at each origin using rp_b observed at that origin
df_aux <- df_aux %>%
  mutate(
    rp_av = if_else(date %in% origins,
                    purrr::map_dbl(rp_b, ~ NA_real_),  # placeholder
                    NA_real_)
  )

# Fill rp_av only on origin dates
# (Use a join approach)
rp_av_tbl <- df_aux %>%
  filter(date %in% origins) %>%
  transmute(date, rp_b_now = rp_b) %>%
  mutate(rp_av = vapply(rp_b_now, avg_future_ar1, numeric(1), c0 = c0, phi = phi, H = H)) %>%
  select(date, rp_av)

df_aux <- df_aux %>%
  left_join(rp_av_tbl, by = "date", suffix = c("", "_new")) %>%
  mutate(rp_av = coalesce(rp_av_new, rp_av)) %>%
  select(-rp_av_new)

# Final series 
df_final <- df_aux %>%
  filter(date >= as.yearmon("2015-01"),
         date <= as.yearmon("2021-10")) %>%
  mutate(infl_e10y = b_10y - rp_av)

# ============================================================
# 6) “Super nice” ggplot figures 
# ============================================================

theme_paper <- function() {
  theme_minimal(base_size = 12) +
    theme(
      panel.grid.minor = element_blank(),
      plot.title = element_text(face = "bold"),
      legend.title = element_blank(),
      legend.position = "right"
    )
}

# ---- Figure 4: Breakeven, IRP, Expected inflation (2015-01 to 2025-01) ----
fig_final <- df_final %>%
  transmute(
    date,
    `Breakeven Inflation (10y)` = b_10y,
    `Risk Premium`    = rp_av,
    `Expected Inflation (10y)`  = infl_e10y
  ) %>%
  tidyr::pivot_longer(-date, names_to = "series", values_to = "value")

fig_final <- fig_final %>%
  mutate(date_plot = as.Date(date))

p_final <- ggplot(fig_final, aes(x = date_plot, y = value, colour = series)) +
  geom_hline(yintercept = 0, linewidth = 0.4, alpha = 0.6) +
  geom_line(linewidth = 1.0) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits = c(-0.01, 0.04)) +
  labs(
    title = "Inflation expectations, breakeven inflation and RP",
    x = NULL,
    y = NULL,
    colour = NULL
  ) +
  theme_paper()

p_final