Yield Curve, with wavelet-denoised trend and 3 mo forecast

################################################################################
### Script to analyze the yield curve for seasonality
###    by /u/MetricT
################################################################################

### 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", "WaveletComp")

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))

### Set start/end dates for our graph.   Note that this only affects the final
### graph, all computations are still done with the full dataset.
date_start <- "2019-01-01" %>% as.Date()
date_end   <- Sys.Date()   %>% as.Date()

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

### A list of possible series we can use
series_10yr <- c("DGS10",  "10 Yr")
series_2yr  <- c("DGS2",   "2 Yr")
series_3mo  <- c("DGS3MO", "3 Mo")

### Select the yield curve series to graph.
series_longmaturity  <- series_10yr
series_shortmaturity <- series_3mo

### Pull data for the long and short maturity asset
data_longmaturity  <- fredr(series_id = series_longmaturity[1],
                            frequency = "d") %>%
                            select(date, value) %>%
                            na.omit()

data_shortmaturity <- fredr(series_id = series_shortmaturity[1],
                            frequency = "d") %>%
                            select(date, value) %>%
                            na.omit()

### Find the common ground between the two datasets
data_full <- data_longmaturity %>%
             inner_join(data_shortmaturity, by = "date") %>%
             `colnames<-`(c("date", "longmaturity", "shortmaturity"))

### Create a text version of the data to use in the graph title
date_txt <- data_full$date %>% tail(n = 1) %>% format("%B %d, %Y")

### Remove any fractional year data at the beginning of the list, it confuses
### the forecast() method for some reason.
month_min <- data_full$date %>% head(n = 1) %>% month()
if (month_min != 1) {
     year_min  <- data_full$date %>% head(n = 1) %>% year()
     year_next <- year_min + 1
     date_rounded <- paste(year_next, "-01-01", sep = "")
     data_full <- data_full %>% filter(date >= date_rounded)
}

data_min <- min(data_full$date)

### Computing the yield curve.  This part is easy...
data_full$yield_curve <- data_full$longmaturity - data_full$shortmaturity

### Pull out the raw data
date               <- data_full %>% pull("date")
longmaturity       <- data_full %>% pull("longmaturity")
shortmaturity      <- data_full %>% pull("shortmaturity")

### We'll need this defined later, saves on typing
ts_date <- c(year(min(date)), month(min(date)))

### Decompose the data and extract trend/seasonal/remainder
mstl_longmaturity  <- data_full$longmaturity  %>%
                      ts(frequency = 250, start = ts_date) %>%
                      mstl()

mstl_shortmaturity <- data_full$shortmaturity %>%
                      ts(frequency = 250, start = ts_date) %>%
                      mstl()

trend_longmaturity      <- mstl_longmaturity  %>% trendcycle()
seasonal_longmaturity   <- mstl_longmaturity  %>% seasonal()
remainder_longmaturity  <- mstl_longmaturity  %>% remainder()

trend_shortmaturity     <- mstl_shortmaturity %>% trendcycle()
seasonal_shortmaturity  <- mstl_shortmaturity %>% seasonal()
remainder_shortmaturity <- mstl_shortmaturity %>% remainder()

### Let's save a yield curve based solely on trend
data_full$yield_curve_trend <- trend_longmaturity - trend_shortmaturity

### and another yield curve based on trend + remainder
data_full$yield_curve_trend_remainder <- trend_longmaturity +
                                         remainder_longmaturity -
                                         trend_shortmaturity -
                                         remainder_shortmaturity

### Denoise the non-trend component with wavelets
wavelet_longmaturity  <- (longmaturity - trend_longmaturity) %>%
                         data.frame() %>%
                         `colnames<-`(c("value")) %>%
                         analyze.wavelet("value", dt = 250)

wavelet_shortmaturity <- (shortmaturity - trend_shortmaturity) %>%
                         data.frame() %>%
                         `colnames<-`(c("value")) %>%
                         analyze.wavelet("value", dt = 250)

r_longmaturity  <- reconstruct(wavelet_longmaturity,  plot.rec = FALSE)
r_shortmaturity <- reconstruct(wavelet_shortmaturity, plot.rec = FALSE)

denoised_longmaturity_remainder  <- r_longmaturity$series$value.r +
                                    r_longmaturity$series$value.trend

denoised_shortmaturity_remainder <- r_shortmaturity$series$value.r +
                                    r_shortmaturity$series$value.trend

### Add the wavelet-denoised yield curve
data_full$yield_curve_denoised  <- trend_longmaturity +
                                   denoised_longmaturity_remainder -
                                   trend_shortmaturity -
                                   denoised_shortmaturity_remainder

### Filter the data to only graph points between our start/end dates at the
### top of the script
data_subset <- data_full %>%
               filter(date >= date_start & date <= date_end)

### Derive a three month forecast for the yield curve
forecast <- data_full$yield_curve %>%
            ts(frequency = 250, start = ts_date) %>%
            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("X68.27.") %>% min(data_subset$yield_curve)
ymax <- forecast$upper %>% data.frame() %>% pull("X68.27.") %>% max(data_subset$yield_curve)

### Labels for our graph
m1 <- paste(series_longmaturity[2], "/",
            series_shortmaturity[2], "Yield Curve for", date_txt, "\n")
m2 <- "With wavelet-denoised trend and three-month forecast"
main <- paste(m1, m2, sep = "")
xlab <- "Year"
ylab <- "Percent"

### Plot the forecast
plot(forecast, main = main, xlab = xlab, ylab = ylab,
     xlim = decimal_date(c(date_start, date_end + months(3))),
     ylim = c(ymin, ymax))

### Add our denoised trend
lines(decimal_date(data_subset$date),  data_subset$yield_curve_denoised, type = "l", col = "blue")
# lines(decimal_date(data_subset$date),  data_subset$yield_curve_trend, type = "l", col = "green")
# lines(decimal_date(data_subset$date),  data_subset$yield_curve_trend_remainder, type = "l", col = "red")

### Add a line at h = 0 denoting the inversion point
abline(h = 0)

### Add lines for the Fed rate changes
abline(v = decimal_date(as.Date("2019-08-01")), lty = 3)
abline(v = decimal_date(as.Date("2019-09-19")), lty = 3)
abline(v = decimal_date(as.Date("2019-10-31")), lty = 3)

### Add lines at 0.4 and 0.2 to aid in seeing how quickly it's declining
abline(h = 0.4, lty = 2)
abline(h = 0.2, lty = 2)