################################################################################
###
### AUTHOR:       Carlo Favero 
### DATE:         Gerzensee 2024
### DESCRIPTION:  This script gives the basics of forecasting with VAR for ex 1
### OUTPUT:           
################################################################################

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

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

## -----------------------------------------------------------------------------
listofpackages <- c("tidyverse","ellipse","reshape2","xts","xlsx","readxl",     # Install and load the required packages
"quantmod","vars","systemfit","stargazer","tseries", "dynlm","sandwich","urca")   
source("preparePackages.R")
preparePackages(listofpackages)

######################
### read data
######################
data<- read_excel('quarterly2025.xlsx', sheet=2, col_names=TRUE)
dates <-seq(as.Date("1978-01-01"),length=nrow(data), by="quarters")
rawdata.ts <- xts(data[,2:15], order.by=dates)
typeof(rawdata.ts$cpi)

######################
### transform data
######################

# Inflation: π (pi)
rawdata.ts$LP     <- log(rawdata.ts$cpi)
rawdata.ts$pi     <- (rawdata.ts$LP - lag.xts(rawdata.ts$LP, 4))*100
rawdata.ts$pi_gap <- rawdata.ts$pi - rawdata.ts$ptr

# Output gap: always check in FRED base year GDP and GDPPOT
rawdata.ts$pgdp=rawdata.ts$gdp/rawdata.ts$gdpc1
rawdata.ts$pgdppot=rawdata.ts$ngdppot/rawdata.ts$gdppot
rawdata.ts$y      <- log(rawdata.ts$gdpc1)
rawdata.ts$y_pot  <- log(rawdata.ts$gdppot)
rawdata.ts$y_gap  <- (rawdata.ts$y - rawdata.ts$y_pot)*100
rawdata.ts$dy_pot <- (rawdata.ts$y_pot - lag.xts(rawdata.ts$y_pot, 4))*100
plot(rawdata.ts$y_gap)
#summary(dynlm(y_gap ~ L(y_gap)))

# Interest rates
rawdata.ts$y3m_cc <-100*log(1+0.01*rawdata.ts$tb3ms)
rawdata.ts$y1_cc <-100*log(1+0.01*rawdata.ts$gs1)
rawdata.ts$y10_cc <-100*log(1+0.01*rawdata.ts$gs10)
rawdata.ts$Dy3m_cc<-(rawdata.ts$y3m_cc - lag.xts(rawdata.ts$y3m_cc, 1))
rawdata.ts$spread<-rawdata.ts$y10_cc-rawdata.ts$y3m_cc
######################
### create a df to use ggplot 
######################

rawdata.df <- data.frame(dates = index(rawdata.ts), Data = coredata(rawdata.ts))
col_names <- c("date", colnames(rawdata.ts))
colnames(rawdata.df) <- col_names

######################
### plots 
######################
# Subset based on a date range
subset_df <- rawdata.df[
  rawdata.df$date >= as.Date("1980-01-01") &
    rawdata.df$date <= as.Date("2025-07-01"),
]

# --- Compute y-limits robustly (ensures 2% line always visible + margin) ---
y_vars <- c("y3m_cc", "y1_cc", "y10_cc", "pi", "ptr")
y_all  <- unlist(subset_df[, y_vars])
y_all  <- y_all[is.finite(y_all)]

y_min <- min(c(y_all, 2), na.rm = TRUE)
y_max <- max(c(y_all, 2), na.rm = TRUE)
pad   <- 0.05 * (y_max - y_min)  # 5% padding

plot <- ggplot(subset_df, aes(x = date)) +
  geom_line(aes(y = y3m_cc, color = "y3m_cc"), linewidth = 1) +
  geom_line(aes(y = y1_cc,  color = "y1_cc"),  linewidth = 1) +
  geom_line(aes(y = y10_cc, color = "y10_cc"), linewidth = 1) +
  geom_line(aes(y = pi,     color = "infl"),   linewidth = 1) +
  geom_line(aes(y = ptr,    color = "long-run"),
            linetype = "dotted", linewidth = 1) +
  
  # --- 2% line: in legend + subtle styling ---
  geom_hline(aes(yintercept = 2, color = "Target 2%"),
             linetype = "dashed", linewidth = 0.9) +
  
  # --- Label the 2% line near the left edge (clean + readable) ---
  annotate(
    "label",
    x = min(subset_df$date),   # left edge
    y = 2,
    label = "Target 2%",
    hjust = -0.05,             # slightly inside the plot
    vjust = 1.4,               # BELOW the line
    label.size = 0.2,
    size = 3.5,
    fill = "white"
  )+

  
  labs(
    title = "Trends in Yields and Inflation",
    x = "Time",
    y = "%"
  ) +
  
  scale_color_manual(
    values = c(
      "y3m_cc"   = "lightblue",
      "y1_cc"    = "blue",
      "y10_cc"   = "darkblue",
      "infl"     = "red",
      "long-run" = "red",
      "Target 2%" = "grey55"   # subtle grey target line
    ),
    name = NULL
  ) +
  
  scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
  
  # --- Ensure visibility + neat margins without distorting data ---
  coord_cartesian(ylim = c(y_min - pad, y_max + pad), expand = FALSE) +
  
  theme_minimal() +
  theme(
    axis.line = element_line(color = "black"),
    panel.grid.major = element_line(color = "grey85", linewidth = 0.35),
    panel.grid.minor = element_line(color = "grey92", linewidth = 0.25),
    legend.position = c(0.9, 0.9),
    legend.justification = c(1, 1),
    axis.text.x = element_text(size = 8)
  )

print(plot)

ggsave(
  filename = "Figure_1.pdf",
  plot = plot,
  device = "pdf",
  width = 5.72,
  height = 3.12
)



plot <- ggplot(subset_df, aes(x = date)) +
  geom_line(aes(y = spread, color = "spread"), size = 1) +
  geom_line(aes(y = Dy3m_cc, color = "change in short-term rates"),linetype = "dotted", size = 1) +
  labs(title = "Spread and change in rates", x = "Time", y = "%") +
  scale_color_manual(values = c("spread" = "blue", "change in short-term rates" = "black"), name = NULL) +
  theme_minimal() +
  theme(axis.line = element_line(color = "black")) +
  scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
  theme(
    legend.position = c(0.9, 0.4),  # Specify custom coordinates (x, y)
    legend.justification = c(1, 1),  # Specify justification for the coordinates
    axis.text.x = element_text(size = 8)  # Adjust the font size here
  )
print(plot)
ggsave(filename = "Figure_2.pdf", plot = plot, device = "pdf",width = 5.72,height=3.12)

######################
### forecasting with a VAR 
######################
params <- c("Dy3m_cc","spread")
X <- rawdata.ts[, c(params)]
X_F<-X["1980-01-01/2006-10-01"]
var_forc <- VAR(X_F, p = 1, type = "const") 
summary(var_forc)
X_FOR <- predict(var_forc, n.ahead = 64)

spread_f<-as.matrix(X_FOR$fcst$spread[,1])
Dy3m_cc_f<-as.matrix(X_FOR$fcst$Dy3m_cc[,1])

dates_for <-seq(as.Date("2007-01-01"),length=64, by="quarters")
spread_f<- xts(spread_f, order.by=dates_for)
Dy3m_cc_f<- xts(Dy3m_cc_f, order.by=dates_for)

rawdata.ts <-merge(rawdata.ts, spread_f, join = 'outer')
rawdata.ts <-merge(rawdata.ts, Dy3m_cc_f, join = 'outer')

rawdata.ts$y3m_cc_f <- NA
rawdata.ts$y10_cc_f <- NA
target_date <- as.Date("2006-10-01")
position <- which(index(rawdata.ts) == target_date)
rawdata.ts[position, "y3m_cc_f"][[1]]=rawdata.ts[position, "y3m_cc"][[1]]
rawdata.ts[position, "y10_cc_f"][[1]]=rawdata.ts[position, "y10_cc"][[1]]

for (i in (position+1):(position + 64)) {
  rawdata.ts[i, "y3m_cc_f"][[1]]=rawdata.ts[i-1, "y3m_cc_f"][[1]]+rawdata.ts[i, "Dy3m_cc_f"][[1]]
  rawdata.ts[i, "y10_cc_f"][[1]]=rawdata.ts[i, "y3m_cc_f"][[1]]+rawdata.ts[i, "spread_f"][[1]]
}


plotforc3m <- data.frame(
  Date = index(rawdata.ts),
  y3m_cc = coredata(rawdata.ts$y3m_cc),
  y3m_cc_f =coredata(rawdata.ts$y3m_cc_f) )

plotforc10y <- data.frame(
  Date = index(rawdata.ts),
  y10_cc = coredata(rawdata.ts$y10_cc),
  y10_cc_f =coredata(rawdata.ts$y10_cc_f) )

plot<-ggplot(plotforc3m, aes(x = Date)) +
  geom_line(aes(y = y3m_cc, color = "Actual"), size = 1) +
  geom_line(aes(y = y3m_cc_f, color = "Forecast"), linetype = "dashed", size = 1) +
  labs(x = "Date", y = "Actual and Forecast", title = "3-month rates") +
  scale_color_manual(values = c("Actual" = "blue", "Forecast" = "red")) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5)  # Center the title
  ) +
  guides(color = guide_legend(title = ""))
plot <- plot +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "black")

print(plot)
ggsave(filename = "Figure_VAR3m.pdf", plot = plot, device = "pdf",width = 5.72,height=3.12)

plot<-ggplot(plotforc10y, aes(x = Date)) +
  geom_line(aes(y = y10_cc, color = "Actual"), size = 1) +
  geom_line(aes(y = y10_cc_f, color = "Forecast"), linetype = "dashed", size = 1) +
  labs(x = "Date", y = "Actual and Forecast", title = "10-year Yields") +
  scale_color_manual(values = c("Actual" = "blue", "Forecast" = "red")) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5)  # Center the title
  ) +
  guides(color = guide_legend(title = ""))
plot <- plot +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "black")

print(plot)
ggsave(filename = "Figure_VAR10y.pdf", plot = plot, device = "pdf",width = 5.72,height=3.12)