ZTPiQDeg

1################################################################################
2### Script to analyze the yield curve for seasonality
3###    by /u/MetricT
4################################################################################
5
6### Set your FRED API key here.  You may request an API key at:
7### https://research.stlouisfed.org/useraccount/apikeys
8api_key_fred <- "WHATEVER_YOUR_FRED_API_KEY_IS"
9
10### Load necessary libraries, installing them if necessary
11packages <- c("fredr", "forecast", "lubridate", "tidyverse",
12              "tsibble", "WaveletComp")
13			  
14new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
15if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
16invisible(lapply(packages, "library", quietly = TRUE, character.only = TRUE))
17
18### Set start/end dates for our graph.   Note that this only affects the final
19### graph, all computations are still done with the full dataset.
20date_start <- "2019-01-01" %>% as.Date()
21date_end   <- Sys.Date()   %>% as.Date()
22
23### Now that the "fredr" library is loaded, set the FRED API key
24fredr_set_key(api_key_fred)
25
26### A list of possible series we can use
27series_10yr <- c("DGS10",  "10 Yr")
28series_2yr  <- c("DGS2",   "2 Yr")
29series_3mo  <- c("DGS3MO", "3 Mo")
30
31### Select the yield curve series to graph.
32series_longmaturity  <- series_10yr
33series_shortmaturity <- series_3mo
34
35### Pull data for the long and short maturity asset
36data_longmaturity  <- fredr(series_id = series_longmaturity[1],
37                            frequency = "d") %>%
38                            select(date, value) %>%
39                            na.omit()
40
41data_shortmaturity <- fredr(series_id = series_shortmaturity[1],
42                            frequency = "d") %>%
43                            select(date, value) %>%
44                            na.omit()
45
46### Find the common ground between the two datasets
47data_full <- data_longmaturity %>%
48             inner_join(data_shortmaturity, by = "date") %>%
49             `colnames<-`(c("date", "longmaturity", "shortmaturity"))
50
51### Create a text version of the data to use in the graph title
52date_txt <- data_full$date %>% tail(n = 1) %>% format("%B %d, %Y")
53
54### Remove any fractional year data at the beginning of the list, it confuses
55### the forecast() method for some reason.
56month_min <- data_full$date %>% head(n = 1) %>% month()
57if (month_min != 1) {
58     year_min  <- data_full$date %>% head(n = 1) %>% year()
59     year_next <- year_min + 1
60     date_rounded <- paste(year_next, "-01-01", sep = "")
61     data_full <- data_full %>% filter(date >= date_rounded)     
62}
63
64data_min <- min(data_full$date)
65
66### Computing the yield curve.  This part is easy...
67data_full$yield_curve <- data_full$longmaturity - data_full$shortmaturity
68
69### Pull out the raw data
70date               <- data_full %>% pull("date")
71longmaturity       <- data_full %>% pull("longmaturity")
72shortmaturity      <- data_full %>% pull("shortmaturity")
73
74### We'll need this defined later, saves on typing
75ts_date <- c(year(min(date)), month(min(date)))
76
77### Decompose the data and extract trend/seasonal/remainder
78mstl_longmaturity  <- data_full$longmaturity  %>%
79                      ts(frequency = 250, start = ts_date) %>%
80                      mstl()
81
82mstl_shortmaturity <- data_full$shortmaturity %>%
83                      ts(frequency = 250, start = ts_date) %>%
84                      mstl()
85
86trend_longmaturity      <- mstl_longmaturity  %>% trendcycle()
87seasonal_longmaturity   <- mstl_longmaturity  %>% seasonal()
88remainder_longmaturity  <- mstl_longmaturity  %>% remainder()
89
90trend_shortmaturity     <- mstl_shortmaturity %>% trendcycle()
91seasonal_shortmaturity  <- mstl_shortmaturity %>% seasonal()
92remainder_shortmaturity <- mstl_shortmaturity %>% remainder()
93
94### Let's save a yield curve based solely on trend
95data_full$yield_curve_trend <- trend_longmaturity - trend_shortmaturity
96
97### and another yield curve based on trend + remainder
98data_full$yield_curve_trend_remainder <- trend_longmaturity +
99                                         remainder_longmaturity - 
100                                         trend_shortmaturity -
101                                         remainder_shortmaturity
102
103### Denoise the non-trend component with wavelets
104wavelet_longmaturity  <- (longmaturity - trend_longmaturity) %>%
105                         data.frame() %>%
106                         `colnames<-`(c("value")) %>%
107                         analyze.wavelet("value", dt = 250)
108
109wavelet_shortmaturity <- (shortmaturity - trend_shortmaturity) %>%
110                         data.frame() %>%
111                         `colnames<-`(c("value")) %>%
112                         analyze.wavelet("value", dt = 250)
113
114r_longmaturity  <- reconstruct(wavelet_longmaturity,  plot.rec = FALSE)
115r_shortmaturity <- reconstruct(wavelet_shortmaturity, plot.rec = FALSE)
116
117denoised_longmaturity_remainder  <- r_longmaturity$series$value.r +
118                                    r_longmaturity$series$value.trend
119
120denoised_shortmaturity_remainder <- r_shortmaturity$series$value.r +
121                                    r_shortmaturity$series$value.trend
122
123### Add the wavelet-denoised yield curve
124data_full$yield_curve_denoised  <- trend_longmaturity +
125                                   denoised_longmaturity_remainder -
126                                   trend_shortmaturity -
127                                   denoised_shortmaturity_remainder
128
129### Filter the data to only graph points between our start/end dates at the
130### top of the script
131data_subset <- data_full %>%
132               filter(date >= date_start & date <= date_end)
133
134### Derive a three month forecast for the yield curve
135forecast <- data_full$yield_curve %>%
136            ts(frequency = 250, start = ts_date) %>%
137            mstl() %>%
138            forecast(h = round(250 * 3 / 12), level = c(68.27)) #, 95.45))
139
140### Find appropriate y limits for our graph
141ymin <- forecast$lower %>% data.frame() %>% pull("X68.27.") %>% min(data_subset$yield_curve)
142ymax <- forecast$upper %>% data.frame() %>% pull("X68.27.") %>% max(data_subset$yield_curve)
143
144### Labels for our graph
145m1 <- paste(series_longmaturity[2], "/",
146            series_shortmaturity[2], "Yield Curve for", date_txt, "\n")
147m2 <- "With wavelet-denoised trend and three-month forecast"
148main <- paste(m1, m2, sep = "")
149xlab <- "Year"
150ylab <- "Percent"
151
152### Plot the forecast
153plot(forecast, main = main, xlab = xlab, ylab = ylab,
154     xlim = decimal_date(c(date_start, date_end + months(3))),
155     ylim = c(ymin, ymax))
156
157### Add our denoised trend
158lines(decimal_date(data_subset$date),  data_subset$yield_curve_denoised, type = "l", col = "blue")
159# lines(decimal_date(data_subset$date),  data_subset$yield_curve_trend, type = "l", col = "green")
160# lines(decimal_date(data_subset$date),  data_subset$yield_curve_trend_remainder, type = "l", col = "red")
161
162### Add a line at h = 0 denoting the inversion point
163abline(h = 0)
164
165### Add lines for the Fed rate changes
166abline(v = decimal_date(as.Date("2019-08-01")), lty = 3)
167abline(v = decimal_date(as.Date("2019-09-19")), lty = 3)
168abline(v = decimal_date(as.Date("2019-10-31")), lty = 3)
169
170### Add lines at 0.4 and 0.2 to aid in seeing how quickly it's declining
171abline(h = 0.4, lty = 2)
172abline(h = 0.2, lty = 2)