One of my Ph.D papers was published recently. It deals with yield curve forecasting.

Here is the code for applying the Nelson-Siegel model to any yield curve.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
create_NS_residuals = function(dat, horiz = 12, initobs = 108, tau = 0.0609){ n = NROW(dat) nyields <- NCOL(dat) a = c(seq(3,24,3),30,36,48,60,72,84,96,108,120) # maturity in months fc = matrix(nrow = n , ncol = 3) # recursive estimation of factors fcfar = matrix(nrow = n , ncol = 3) ; lml = list() # forecasted factors for the AR ypred = matrix(nrow = n , ncol = nyields) load2 = (1 - exp(-tau*a))/(tau*a) ; load3 = (1 - exp(-tau*a))/(tau*a) - exp(-tau*a) for (i in (1):n){ lm1 = lm(dat[i,]~load2+load3) # Estimation of the factors at time t. fc[i,] = lm1$coef } for (i in 1:(n-initobs-horiz)){ for (h in 1:3){ # get the coefficients only when we can, note we stop training before time # t: "i-horiz+initobs" lml[[h]] = lm(fc[(1+horiz):(i+initobs),h]~fc[1:(i-horiz +initobs),h]) fcfar[(initobs+i+horiz),h] = t(lml[[h]]$coef)%*%c(1,fc[(i+initobs),h]) } ypred[(i+initobs+horiz),] = fcfar[(i+initobs+horiz),]%*%t(cbind(rep(1,nyields),load2, load3)) } resmat = dat - ypred # the residual matrix list(ypred = ypred, resmat = resmat,fc=fc,fcfar=fcfar) } |

The previous function returns an out-of-sample residual matrix. The paper maintains that there is still information contained in those residuals, and that this information can be extracted to improve forecasting performance, like so:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
# Once you got the residual matrix we predict the bias using AR(1) res_pred = function(dat,residualmat,ypred, horiz = 1, initobs = 108){ # The args here are the residual matrix, the predictions and the raw data n = NROW(residualmat) nyields = NCOL(residualmat) respred0 = matrix(nrow = n, ncol = nyields) # Container ynpred = matrix(nrow = n, ncol = nyields) # Container for (j in 1:(nyields)){ for (i in 1:(n-initobs-horiz) ){ # lm1 is the AR(1) of the residuals, # note: "i-horiz+initobs" (we stop before we start forecasting) lm1 = lm(residualmat[(1+horiz):(i+initobs),j]~residualmat[1:(i-horiz+initobs),j]) # predict the new residual: respred0[(i+initobs+horiz),j] = t(lm1$coef)%*%c(1,residualmat[(i+initobs),j]) # Adjust the forecast, equation (5) in the published version ynpred[(i+initobs+horiz),j] = ypred[(i+initobs+horiz),j]+respred0[(i+initobs+horiz),j] } } nresmat = dat - ynpred # The new residuals, needed for evaluation list(ynpred = ynpred, nresmat = nresmat) } |

When used, please cite. Enjoy.

Here is a 5 minutes long presentation about the paper: