Analyzing the yield curve for seasonality

################################################################################
### Script to analyze the yield curve for seasonality
###    by <Mathew.Binkley@Vanderbilt.edu>
################################################################################

### Set your FRED API key here.  You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
api_key_fred <- "WHATEVER_YOUR_FRED_API_KEY_IS"

### Load necessary libraries, installing them if necessary
packages <- c("fredr", "forecast", "lubridate", "tidyverse", "tsibble")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
invisible(lapply(packages, "library", quietly = TRUE, character.only = TRUE))

### Now that the "fredr" library is loaded, set the FRED API key
fredr_set_key(api_key_fred)

### Select the yield curve series to graph.
### * 10 yr/3 mo [T10Y3M] curve is considered the most accurate
### * 10 yr/2 y3 [T10Y2Y] curve was used in the original research
series <- "T10Y3M"

### Pull yield curve data from FRED
data    <- fredr(series_id = series,
                 observation_start = as.Date("1982-01-01"),
                 frequency = "d") %>%
                 select(date, value) %>%
                 na.omit()

### There are 250 trading days per year, so use that for our frequency
### when constructing the time series
data_ts <- ts(data$value,
              frequency = 250,
              start = c(year(min(data$date)), month(min(data$date))))

### Now break our data down into trend, seasonal, and remainder
data_decomposed <- mstl(data_ts)

actual    <- ts(data$value)
trend     <- trendcycle(data_decomposed)
seasonal  <- seasonal(data_decomposed)
remainder <- remainder(data_decomposed)

################################################################################
### Graph 1:  Yield Curve using raw data
################################################################################

plot(data, type = "l", xlab = "year", ylab = "%",
     main = paste(series, "Yield Curve (Not Seasonally Adjusted)"))
abline(h = 0)

readline(prompt = "Press [ENTER] to continue...")

################################################################################
### Graph 2:  Decomposed yield curve into trend, seasonal, remainder
################################################################################

plot(data_decomposed,
     main = paste(series, "Yield Curve Decomposition"), xlab = "Year")

readline(prompt = "Press [ENTER] to continue...")

################################################################################
### Graph 3:  The seasonal component of the yield curve data
################################################################################

### Find the amplitude over the past 4 years
subset <- tail(seasonal, n = 4 * 250)
amplitude <- round(max(subset) - min(subset), digits = 4)

title <- paste(series, " Yield Curve Seasonal Component\n[Amplitude = ",
           amplitude, "%]", sep = "")

plot(decimal_date(data$date), seasonal, type = "l", lty = 1,
     main = title, xlab = "Year", ylab = "%",
     xlim = c(2018.0, 2019.0),
     ylim = c(min(subset), max(subset)))

readline(prompt = "Press [ENTER] to continue...")

################################################################################
### Graph 4:  Deseasonalized yield curve in our region of interest (2019-2020)
################################################################################
### Near yield curve inversions, the yield curve value is always near 0%.   So
### the Seasonal component can dominate the data and distort our view as to
### what's actually going on.   So let's remove it and plot trend + remainder

### Find the amplitude over the past year
subset_deseasonalized <- tail(trend + remainder, n = 1 * 250)
subset_actual         <- tail(actual,            n = 1 * 250)
subset_date           <- tail(data$date,         n = 1 * 250)

### Find appropriate y limits for our graph
min <- min(subset_deseasonalized, subset_actual)
max <- max(subset_deseasonalized, subset_actual)

t1 <- "Yield Curve\n[Original in dotted blue, Deseasonalized in black]\n"
title <- paste(series, t1)

plot(decimal_date(subset_date), subset_deseasonalized, type = "l",
     main = title, xlab = "Year", ylab = "%",
     ylim = c(min, max))

abline(h = 0, col = "red") # Draw a line at 0% to highlight inversions

lines(decimal_date(subset_date), subset_actual,
      type = "l", col = "blue", lty = 3)

readline(prompt = "Press [ENTER] to continue...")

################################################################################
### Graph 5:  3 month forecast of yield curve based on its history
################################################################################

forecast <- data_ts %>%
            mstl() %>%
            forecast(h = round(250 * 3 / 12), level = c(68.27, 95.45))

### Find appropriate y limits for our graph
ymin <- forecast$lower %>% data.frame() %>% pull("X95.45.") %>% min()
ymax <- forecast$upper %>% data.frame() %>% pull("X95.45.") %>% max()

plot(forecast,
     main = paste(series, "Yield Curve Three Month Forecast"),
     xlab = "Year",
     ylab = "Percent",
     xlim = c(2019, 2020.25),
     ylim = c(ymin, ymax))

abline(h = 0)