Yield curve forecasting

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.

 
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:

 
# 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: