Financial Econometrics: Nelson-Seigel 3-factor Model & Forecasting

In class, we discussed the Nelson-Seigel 3-factor model of the yield curve. This page estimates the model using daily data from 1/02/2000 to 10/31/2023 for the following maturities of US Treasuries:3 months, 6 months, 1 year, 2 years, 3 years, 5 years, 10 years, 20 years, and 30 years. (The data is available at: https://fred.stlouisfed.org/). The results are then used to forecast the yield curve for November 2023 through an AR(1) model.

1. Data Fetching & Pre-processing

We fetch all the necessary data into a list of lists indexed by maturity dates, then combine into a final processed dataframe.

#devtools::install_github("rstudio/tensorflow")
#devtools::install_github("rstudio/keras")
library('quantmod')
library(forecast)
library(dplyr)
library(lubridate)
library(YieldCurve)
library(xts)
library(plotly)

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

treasuries <- list()

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

##            DGS3MO DGS6MO DGS1 DGS2 DGS3 DGS5 DGS10 DGS20 DGS30
## 2000-01-03   5.48   5.81 6.09 6.38 6.42 6.50  6.58  6.94  6.61
## 2000-01-04   5.43   5.75 6.00 6.30 6.34 6.40  6.49  6.84  6.53
## 2000-01-05   5.44   5.74 6.05 6.38 6.43 6.51  6.62  6.95  6.64
## 2000-01-06   5.41   5.69 6.03 6.35 6.39 6.46  6.57  6.86  6.58
## 2000-01-07   5.38   5.66 6.00 6.31 6.35 6.42  6.52  6.82  6.55
## 2000-01-10   5.42   5.64 6.07 6.38 6.42 6.49  6.57  6.86  6.59

2.1 Nelson-Siegel Model

Nelson and Siegel fit the yield at any maturity as follows: \[\begin{align} &f_t(\tau) = \beta_{1,t} + \beta_{2,t} \left( \frac{1 - e^{-\lambda_t \tau}}{\lambda_t \tau} \right) + \beta_3 \left( \frac{1 - e^{-\lambda_t \tau}}{\lambda_t \tau} - e^{-\lambda_t \tau} \right)\\ &-\tau \text{ is the time to maturity}.\\ &-\beta_{1,t}, \beta_{2,t} \beta_{3,t} \text{ are the level(long term), slope(short term),} \\ &\text{ and curvature(medium term) parameters of the yield curve}.\\ & -\lambda \text{ is the exponential decay rate}. \end{align}\] And we can use \(R\)’s \(YieldCurve\) package to estimate the parameters:

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

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

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

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

## plotting
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")
lines(maturity.Fed,y, col=2)
legend("topleft", legend=c("observed yield curve","fitted yield curve"),
col=c(1,2),lty=1)
grid()

2.2 Nelson-Siegel Fitted & Actual（3D surface)

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

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

# Plot the actual yields surface
p <- plot_ly(x = ~rev(maturity.Fed), y = ~yield$DATE, z = ~actual_yield_matrix, type = 'surface',
             colorscale = 'Viridis',
             name = 'Actual Yields') %>%
  
      add_surface(x = ~rev(maturity.Fed), y = ~yield$DATE, z = ~fitted_yield_matrix,
                  colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(128,0,64)")),
                  name = 'Fitted Yields') %>%
      layout(title = "Yield Surface Over Time and Maturity",
             scene = list(xaxis = list(title = "Maturity (Months)"),
                          yaxis = list(title = "Date"),
                          zaxis = list(title = "Yield (%)")),
             width = '70%')

# Render the plot
p

3.1 Forecasting

Following Diebold & Li (2006), we forecast the yield curve with the following specifications: \[\begin{align} \hat{y}_{t+h|t}(\tau) &= \hat{\beta}_{1,t+h|t} + \hat{\beta}_{2,t+h|t} \left[ \frac{1-e^{-\lambda \tau}}{\lambda \tau} \right] + \hat{\beta}_{3,t+h|t} \left[ \frac{1-e^{-\lambda \tau}}{\lambda \tau} - e^{-\lambda \tau} \right]\\ \hat{\beta}_{i,t+1} &= c_i + \gamma_i \hat{\beta}_{i,t} + \epsilon_{i,t+1} \end{align}\] Where we forecast the NS factors (\(\beta\)s) as an \(ARIMA(p, d, q)\) process to handle potential non-stationarity:

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

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

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

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

# str(forecast_beta_0) # checking the structure of output

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

## Beta 1 ARIMA(p, d, q): 5, 1, 0

# For arima_beta_1
non_seasonal_order_beta_1 <- arima_beta_1$arma[1:3]
cat("Beta 2 ARIMA(p, d, q):", paste(non_seasonal_order_beta_1, collapse = ", "), "\n")

## Beta 2 ARIMA(p, d, q): 4, 1, 0

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

## Beta 3 ARIMA(p, d, q): 2, 1, 0

3.2 Forecasting the yield using the forecasted factors

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

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

# Calculate the yield curve for each day in Nov
for (i in 1:nrow(yield_forecasts)) {
  for (j in 1:ncol(yield_forecasts)) {
    yield_forecasts[i, j] <- calculate_yield(
      beta_0 = forecast_beta_0$mean[i],
      beta_1 = forecast_beta_1$mean[i],
      beta_2 = forecast_beta_2$mean[i],
      lambda = lambda,
      maturity = maturity[j]
    )
  }
}

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

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