The prompt and base code for this project can be found on Datacamp. The basis of this is to explore the volatility of Zero Coupon Bond Yields (from year maturities 1-30) and see how that has changed since the 1960’s. Below, I loaded the necessary librares and did some data exploration of the time series data obtained from Quandl. This was obtained making an API call.

Necessary packages

options(warn = -1) 

library(dygraphs)
library(Quandl)
library(dplyr)
library(rugarch)
library(plotly)
library(RColorBrewer)
#get dataset, filter for coronavirus dates
df <- Quandl("FED/SVENY", api_key="15cxBvvbCzucYDswsDfJ")
covid_dates <- subset(df, Date > '2020-01-01' & Date < '2020-09-07')

str(df)
## 'data.frame':    14765 obs. of  31 variables:
##  $ Date   : Date, format: "2020-08-28" "2020-08-27" ...
##  $ SVENY01: num  0.148 0.156 0.159 0.157 0.154 ...
##  $ SVENY02: num  0.137 0.163 0.16 0.158 0.157 ...
##  $ SVENY03: num  0.173 0.205 0.197 0.197 0.193 ...
##  $ SVENY04: num  0.237 0.269 0.255 0.256 0.247 ...
##  $ SVENY05: num  0.317 0.344 0.324 0.326 0.311 ...
##  $ SVENY06: num  0.404 0.426 0.399 0.4 0.38 ...
##  $ SVENY07: num  0.494 0.51 0.476 0.475 0.45 ...
##  $ SVENY08: num  0.583 0.594 0.553 0.55 0.52 ...
##  $ SVENY09: num  0.67 0.675 0.627 0.623 0.589 ...
##  $ SVENY10: num  0.753 0.754 0.7 0.692 0.656 ...
##  $ SVENY11: num  0.832 0.83 0.769 0.759 0.72 ...
##  $ SVENY12: num  0.907 0.901 0.835 0.823 0.782 ...
##  $ SVENY13: num  0.978 0.969 0.898 0.885 0.842 ...
##  $ SVENY14: num  1.044 1.033 0.959 0.943 0.9 ...
##  $ SVENY15: num  1.107 1.094 1.016 0.998 0.954 ...
##  $ SVENY16: num  1.17 1.15 1.07 1.05 1.01 ...
##  $ SVENY17: num  1.22 1.21 1.12 1.1 1.06 ...
##  $ SVENY18: num  1.27 1.26 1.17 1.15 1.11 ...
##  $ SVENY19: num  1.32 1.3 1.22 1.2 1.15 ...
##  $ SVENY20: num  1.37 1.35 1.26 1.24 1.2 ...
##  $ SVENY21: num  1.41 1.39 1.3 1.28 1.24 ...
##  $ SVENY22: num  1.45 1.43 1.34 1.32 1.28 ...
##  $ SVENY23: num  1.49 1.47 1.38 1.36 1.31 ...
##  $ SVENY24: num  1.53 1.51 1.42 1.39 1.35 ...
##  $ SVENY25: num  1.56 1.54 1.45 1.43 1.39 ...
##  $ SVENY26: num  1.6 1.57 1.48 1.46 1.42 ...
##  $ SVENY27: num  1.63 1.6 1.51 1.49 1.45 ...
##  $ SVENY28: num  1.66 1.63 1.54 1.52 1.48 ...
##  $ SVENY29: num  1.68 1.66 1.57 1.55 1.51 ...
##  $ SVENY30: num  1.71 1.69 1.6 1.57 1.53 ...
##  - attr(*, "freq")= chr "daily"

Plotting daily estimates for all zero coupon yields in 2020

dfasxts <- as.xts(x = df[, -1], order.by = df$Date)
covid_dates <- as.xts(x = covid_dates[, -1], order.by = covid_dates$Date)
dygraph(covid_dates, main = "All Zero Coupon Yields (1-30) 2020", ylab = "Value") %>%
            dyAxis('x', axisLabelFontSize = 12) %>%
            dyRangeSelector()
df$Date <- as.Date(df$Date)
df$year <- format(as.Date(df$Date, format="%m/%d/%Y"),"%Y")
df <- select(df, -Date)
df <- na.omit(df)

df <- df %>%
  group_by(year) %>%
  summarise_all(mean)
SVENY01 <- df$SVENY01
SVENY10 <- df$SVENY10
SVENY30 <- df$SVENY30


data <- data.frame(df, SVENY01, SVENY10, SVENY30)
data <- select(data, year, SVENY01, SVENY10, SVENY30)
fig <- plot_ly(data, x = data$year, y = ~SVENY01, name = 'SVENY01', type = 'scatter', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY10, name = 'SVENY10', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY30, name = 'SVENY30', mode = 'lines + markers')
fig <- fig %>% layout(title = "",
         xaxis = list(title = "year"),
         yaxis = list (title = ""))

1, 10, and 30 year maturity bond yield estimate averages across time

fig
row.names(df) <- df$year
df <- select(df, -year)
df_matrix <- data.matrix(df)

df_heatmap <- heatmap(df_matrix, Rowv=NA, Colv=NA, col = brewer.pal(6, "Greys"), scale="column", margins=c(4,1), main = "Bond Yield Estimate Averages Heatmap")

# plotting the evaluation of bond yields
library(viridisLite)
yields <- dfasxts
plot.type <- "single"
plot.palette <- magma(n = 30)
asset.names <- colnames(dfasxts)
plot.zoo(x = dfasxts, plot.type = "single", col = plot.palette, ylab = "", xlab = "")
legend(x = "topleft", legend = asset.names, col = plot.palette, cex = 0.45, lwd = 3)

#plotting yield changes from differentiated time series
dfasxts_d <- diff(dfasxts)

plot.zoo(x = dfasxts_d, plot.type = "multiple", ylim = c(-0.5, 0.5), cex.axis = 0.7, ylab = 1:30, col = plot.palette, main = "", xlab = "")

dfasxts <- dfasxts_d["2000/",]
x_1 <- dfasxts[,"SVENY01"]
x_20 <- dfasxts[, "SVENY20"]

# Plot the autocorrelations of the yield changes)
par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,2))
acf_1 <- acf(x_1)
acf_20 <- acf(x_20)

# Plot the autocorrelations of the absolute changes of yields
acf_abs_1 <- acf(abs(x_1))
acf_abs_20 <- acf(abs(x_20))

# if autocorrelation is close to one, then past data value will be very close to next day value
# if autocorrelation is close to zero, then past data value will not have an effect on next day value
#GARCH modeling ; Fitting one year and twenty year maturities
spec <- ugarchspec(distribution.model = "sstd")


fit_1 <- ugarchfit(x_1, spec = spec)


vol_1 <- sigma(fit_1)
res_1 <- scale(residuals(fit_1, standardize = TRUE)) * sd(x_1) + mean(x_1)


merge_1 <- merge.xts(x_1, vol_1, res_1)
plot.zoo(merge_1, xlab = "Year")

####################

fit_20 <- ugarchfit(x_20, spec = spec)


vol_20 <- sigma(fit_20)
res_20 <- scale(residuals(fit_20, standardize = TRUE)) * sd(x_20) + mean(x_20)


merge_20 <- merge.xts(x_20, vol_20, res_20)
plot.zoo(merge_20, xlab = "Year")

par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,1))
hist(res_1)
hist(res_20)

ugarchspec()
## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : sGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(1,0,1)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE
#specs of the GARCH Model, and histogram if year one and year twenty residuals
#year one model has greater and larger volatility than the twenty year model. Lets plot a density diagram and examine the distributionsx
density_x_1 <- density(x_1)
density_res_1 <- density(res_1)


plot(density_x_1)
lines(density_res_1, col = "red")


norm_dist <- dnorm(seq(-0.4, 0.4, by = .01), mean = mean(x_1), sd = sd(x_1))
lines(seq(-0.4, 0.4, by = .01), 
      norm_dist, 
      col = "darkblue"
     )

# Add legend
legend <- c("Before GARCH", "After GARCH", "Normal distribution")
legend("topleft", legend = legend, 
       col = c("black", "red", "darkblue"), lty=c(1,1))

distribution <- qnorm


qqnorm(x_1, ylim = c(-0.5, 0.5))
qqline(x_1, distribution = distribution, col = "darkgreen")


par(new=TRUE)
qqnorm(res_1 * 0.614256270265139, col = "red", ylim = c(-0.5, 0.5))
qqline(res_1 * 0.614256270265139, distribution = distribution, col = "darkgreen")
legend("topleft", c("Before GARCH", "After GARCH"), col = c("black", "red"), pch=c(1,1))

Conclusion

To begin the time series data was pulled from Quandl using an API call. Out of curiosity, I plotted the daily estimates for all bond yields in the year 2020 using a dygraph. From there I averaged the bond yield estimates for all years to see how they have changed over time. Using plotly I displayed the averages of SVENY01 SVENY10, and SVENY30 per year. A heatmap was also generated to display the average estimates for all zero coupon bond yields. Moving on, I plotted the evaluation of bond yields. We can see the level of bond yields for some maturities but to recognize volatility we need to see changes in the yield levels. We need to differentiate the time series and make the time series independent of time. After plotting the differentiated series, we can dive into statistics. Testing for autocorrelation on the previous graph reveals two things; If autocorrelation is close to one, then past data value will be very close to next day value, if autocorrelation is close to zero, then past data value will not have an effect on next day value. After this, I utilized Garch (Generalized AutoRegressive Conditional Heteroskedasticity) which specializes in handeling changing volatility in financial time series data. I fit Garch modeling on the one year and twenty year maturities. The first year model is showing more erratic behavior than the twenty year model, this is because the year one model has greater volatility. Then we wanted to see what the fitted model did with the distribution. Garch modeling brought the residuals closer to normal distribution. Overall, Garch modeling also revealed how volatility has changed over time.