Nelson-Seigel 3-factor Model & Forecasting

1
library('quantmod')
2
library(forecast)
3
library(dplyr)
4
library(lubridate)
5
library(YieldCurve)
6
library(xts)
7
library(plotly)
8

9
# query the data
10
maturities <- c("DGS3MO", "DGS6MO", "DGS1", "DGS2", "DGS3", "DGS5", "DGS10", "DGS20", "DGS30")
11
start <- as.Date('2000-01-02')
12
end <- as.Date('2023-10-31')
13

14
treasuries <- list()
15

16
for (tick in maturities){
17
  # call data
18
  getSymbols(tick, src = 'FRED', from = start, to = end)
19
  treasuries[[tick]] <- get(tick) # R lists can be indexed by name
20
}
21
# combine the indexed list to a df
22
combined_data <- do.call(merge, c(treasuries, all = TRUE))
23
print(head(combined_data))

1
# NS yield as a function of factors and maturities
2
# to be used later
3
calculate_yield <- function(beta_0, beta_1, beta_2, lambda, maturity) {
4
  term1 <- (1 - exp(-lambda * maturity)) / (lambda * maturity)
5
  term2 <- term1 - exp(-lambda * maturity)
6
  yield <- beta_0 + (beta_1 * term1) + (beta_2 * term2)
7
  return(yield)
8
}
9

10
# More cleaning & convert xts to data frame
11
# make index into a column
12
combined_data_df <- data.frame(DATE = index(combined_data), coredata(combined_data))
13
combined_data_df$DATE <- as.Date(combined_data_df$DATE)
14

15
# prepare data for YieldCurve package
16
yield <- combined_data_df # rename for easier reference
17
yield$DATE <- as.Date(yield$DATE, format = '%Y-%m-%d')
18
curve <- xts(yield[, -1], order.by = yield$DATE)
19
curve_clean <- na.omit(curve) # drop NAs, 5962 obs remaining
20
maturity.Fed <- c(3/12, 0.5, 1,2,3,5,10, 20, 30 )
21

22
# Nielson-Siegel
23
NSParameters <- Nelson.Siegel( rate=curve_clean, maturity=maturity.Fed)
24
y <- NSrates(NSParameters[5962,], maturity.Fed)
25

26
## plotting
27
plot(maturity.Fed, curve_clean[5962,], main="Fitting Nelson-Siegel yield curve (23/10/31)", xlab=c("Maturity (months)"), ylab=c('Market Yield'), type="o")
28
lines(maturity.Fed,y, col=2)
29
legend("topleft", legend=c("observed yield curve","fitted yield curve"),
30
col=c(1,2),lty=1)
31
grid()

1
# To see our result in 3D, we apply the calculate_yield function to each row in NSParameters
2
fitted_yield_matrix <- t(apply(NSParameters, 1, function(params) {
3
  sapply(maturity.Fed, calculate_yield, beta_0 = params['beta_0'], beta_1 = params['beta_1'], beta_2 = params['beta_2'], lambda = params['lambda'])
4
}))
5

6
# More cleaning
7
actual_yield_matrix <- as.matrix(curve_clean)
8

9
# Plot the actual yields surface
10
p <- plot_ly(x = ~rev(maturity.Fed), y = ~yield$DATE, z = ~actual_yield_matrix, type = 'surface',
11
             colorscale = 'Viridis',
12
             name = 'Actual Yields') %>%
13

14
      add_surface(x = ~rev(maturity.Fed), y = ~yield$DATE, z = ~fitted_yield_matrix,
15
                  colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(128,0,64)")),
16
                  name = 'Fitted Yields') %>%
17
      layout(title = "Yield Surface Over Time and Maturity",
18
             scene = list(xaxis = list(title = "Maturity (Months)"),
19
                          yaxis = list(title = "Date"),
20
                          zaxis = list(title = "Yield (%)")),
21
             width = '70%')
22

23
# Render the plot
24
p

1
# Extract the individual beta factors as separate time series
2
beta_0_ts <- ts(NSParameters$beta_0, start = c(year(start(NSParameters)), month(start(NSParameters))), frequency = 12)
3
beta_1_ts <- ts(NSParameters$beta_1, start = c(year(start(NSParameters)), month(start(NSParameters))), frequency = 12)
4
beta_2_ts <- ts(NSParameters$beta_2, start = c(year(start(NSParameters)), month(start(NSParameters))), frequency = 12)
5

6
# Fit AR(1) models to each of the factors
7
arima_beta_0 <- auto.arima(beta_0_ts)
8
arima_beta_1 <- auto.arima(beta_1_ts)
9
arima_beta_2 <- auto.arima(beta_2_ts)
10

11
# Number of days to forecast
12
forecast_horizon <- length(seq(ymd("2023-11-01"), ymd("2023-11-30"), by = "days"))
13

14
# Forecast the factors for the horizon
15
forecast_beta_0 <- forecast(arima_beta_0, h = forecast_horizon)
16
forecast_beta_1 <- forecast(arima_beta_1, h = forecast_horizon)
17
forecast_beta_2 <- forecast(arima_beta_2, h = forecast_horizon)
18
lambda <- mean(NSParameters$lambda)
19

20
# For arima_beta_0
21
non_seasonal_order_beta_0 <- arima_beta_0$arma[1:3]
22
cat("Beta 1 ARIMA(p, d, q):", paste(non_seasonal_order_beta_0, collapse = ", "), "\n")
23

24
# For arima_beta_1
25
non_seasonal_order_beta_1 <- arima_beta_1$arma[1:
26

27
3]
28
cat("Beta 2 ARIMA(p, d, q):", paste(non_seasonal_order_beta_1, collapse = ", "), "\n")
29

30
# For arima_beta_2
31
non_seasonal_order_beta_2 <- arima_beta_2$arma[1:3]
32
cat("Beta 3 ARIMA(p, d, q):", paste(non_seasonal_order_beta_2, collapse = ", "), "\n")

1
# Define the maturities for the yield curve
2
maturity <- c(3/12, 0.5, 1,2,3,5,10, 20, 30 )
3

4
# Initialize a matrix to hold the yield curve forecasts
5
forecast_dates <- seq.Date(from = as.Date("2023-11-01"), to = as.Date("2023-11-30"), by = "day")
6
yield_forecasts <- matrix(nrow = forecast_horizon, ncol = length(maturity))
7

8
# Calculate the yield curve for each day in Nov
9
for (i in 1:nrow(yield_forecasts)) {
10
  for (j in 1:ncol(yield_forecasts)) {
11
    yield_forecasts[i, j] <- calculate_yield(
12
      beta_0 = forecast_beta_0$mean[i],
13
      beta_1 = forecast_beta_1$mean[i],
14
      beta_2 = forecast_beta_2$mean[i],
15
      lambda = lambda,
16
      maturity = maturity[j]
17
    )
18
  }
19
}
20

21
# Combine the forecasted yields with the dates
22
forecast_df <- data.frame(date = forecast_dates, yield_forecasts)
23

24
# Melt the data frame to a long format for easier plotting and analysis
25
forecast_df_long <- reshape2::melt(forecast_df, id.vars = 'date', variable.name = 'maturity', value.name = 'yield')

1
library(plotly)
2
yield_matrix <- as.matrix(forecast_df[,-1])  # Excluding the date column
3

4
# Use plot_ly to create an interactive 3D plot
5
plot_ly(z = ~yield_matrix, x = ~rev(forecast_df$date), y = ~maturities, type = "surface") %>%
6
  layout(title = "Forecasted Yield Surface (November 2023)",
7
         scene = list(xaxis = list(title = "Date"),
8
                      yaxis = list(title = "Maturity (Months)"),
9
                      zaxis = list(title = "Yield (%)")),
10
         width = '70%')