Risk Management Techniques: Principal Component Analysis of Interest Rates Using R

1 Purpose

The purpose of this project is demonstrate principal component analysis of interest rates using R. This includes accessing data from public sources and performing the necessary calculations.

Principal component analysis is frequently used to capture the variability in the movement of interest rates along the term structure. Understanding the varability allows for creation of stressed interest rate term structures that can be applied as part of a risk management program wherever the uncertainty of interest rates is a concern. Varying degrees of stress can be applied in the fashion of analytical Value-at-Risk (VaR) techniques. This exercise focuses on the 99.5th percentile as this degree is commonly the focus for many insurance companies and is aligned with benchmarks set by industry standards such as Solvency II Standard Formula. This technique is described in many places but this exercise will follow the process outlined in How to Estimate and Calibrate Analytical VaR for Interest Rate Risk. Please refer to this helpful document as efforts have been made to avoid repeating what has already been described there.

2 Interest Rate Data Sources

LIBOR interest rate data is obtained for the following terms: 1-month, 3-month, 6-month, and 1-year. The source of the data is the Federal Reserve Bank of St. Louis and is pulled via API’s the source makes available. To obtain the data a user must register through the website and request an API key. Data can then be queried and retrieved in either JSON or XML format. More information related to usage of the API can found on the sources website.

Swap rate data is obtained for the following terms: 1-year, 3-year, 5-year, 7-year, 10-year, and 30-year. The source of the data is the Federal Reserve. The data is available in the form of downloadable comma separated value (csv) files.

3 R Libraries

This project utilizes xts, jsonlite, ggplot2, and reshape2 packages. The xts package is useful for storing time series data while jsonlite allows for reading data stored in JSON format. The latter is especially helpful for reading LIBOR data from the source.The ggplot2 package is a well known plotting package used to generate the graphs. The reshape2 package is used to structure data to allow for easier computation or plotting.

require(xts)
require(jsonlite)
require(ggplot2)
require(reshape2)

4 Load Raw Interest Rate Data

Interest rate data is loaded into an xts time-series object. This object is useful for data where values correspond to dates.

4.1 Load Libor Data

json_file <- "https://api.stlouisfed.org/fred/series/observations?series_id=USD1MTD156N&frequency=m&aggregation_method=eop&file_type=json&api_key=29fb098284eda4e97e4678353c8f25f2"
json_data <- fromJSON(txt=json_file)
libor_1m <- xts(as.numeric(json_data$observations$value), as.yearmon(json_data$observations$date, format='%Y-%m-%d'))

json_file <- "https://api.stlouisfed.org/fred/series/observations?series_id=USD3MTD156N&frequency=m&aggregation_method=eop&file_type=json&api_key=29fb098284eda4e97e4678353c8f25f2"
json_data <- fromJSON(txt=json_file)
libor_3m <- xts(as.numeric(json_data$observations$value), as.yearmon(json_data$observations$date, format='%Y-%m-%d'))
libor <- merge(libor_1m, libor_3m, all=FALSE)

json_file <- "https://api.stlouisfed.org/fred/series/observations?series_id=USD6MTD156N&frequency=m&aggregation_method=eop&file_type=json&api_key=29fb098284eda4e97e4678353c8f25f2"
json_data <- fromJSON(txt=json_file)
libor_6m <- xts(as.numeric(json_data$observations$value), as.yearmon(json_data$observations$date, format='%Y-%m-%d'))
libor <- merge(libor, libor_6m, all=FALSE)

4.2 Load Swap Data

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=86f09fbc8cfcc804b71a3d827a5fddc8&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_1yr.csv")
data <- read.csv("swap_1yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_1yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=81223b01cb8c943e81dddfbd03cd2199&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_2yr.csv")
data <- read.csv("swap_2yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_2yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap_1yr, swap_2yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=23bf93b6527f59eff7f778462706d86d&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_3yr.csv")
data <- read.csv("swap_3yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_3yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_3yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=025cf7daa62487392306aa9fb6f3d830&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_4yr.csv")
data <- read.csv("swap_4yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_4yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_4yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=05629c455ea42adce05a1c4861e0ceaa&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_5yr.csv")
data <- read.csv("swap_5yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_5yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_5yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=1f1346dca868bbf046c10ad5b912eae1&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_7yr.csv")
data <- read.csv("swap_7yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_7yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_7yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=8f401407ccab6d90aeb1db13ffe746b5&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_10yr.csv")
data <- read.csv("swap_10yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_10yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_10yr)

download.file("http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=4f2a3bbeb30b5b4fbcb019a11a513fec&lastObs=&from=&to=&filetype=csv&label=include&layout=seriescolumn", destfile="swap_30yr.csv")
data <- read.csv("swap_30yr.csv", skip=5, col.names=c("Date", "Rate"))
swap_30yr <- xts(data$Rate, as.yearmon(data$Date, format='%Y-%m'))
swap <- merge(swap, swap_30yr)

4.3 Combine Libor and Swap Data

The available interest rate data is now combined into a single xts object

rates <- merge(libor, swap, all=FALSE)

5 Graph the Rates

zoo.rates <- as.zoo(rates)
tsRainbow <- rainbow(ncol(zoo.rates))
plot(x=zoo.rates,col=tsRainbow, main="Historical Interest Rates", ylab="Rates", xlab="Time", plot.type="single", xaxt ="n", yaxt="n")
axis(1)
pts<- pretty(zoo.rates)
axis(2, at=pts, labels=paste(pts, "%", sep=""))
legend(x="topright", legend=colnames(rates), lty=1, col=tsRainbow, cex=.64,bty="n")

6 Interpolate for Half-Year Periods

swap <- interpolateSwapsPeriods(1,30,.5, swap)
rates<- merge(libor, swap, all=FALSE)
zoo.rates <- as.zoo(rates)
tsRainbow <- rainbow(ncol(zoo.rates))
plot(x=zoo.rates,col=tsRainbow, main="Historical Interest Rates", ylab="Rates", xlab="Time", plot.type="single", xaxt ="n", yaxt="n")
axis(1)
pts<- pretty(zoo.rates)
axis(2, at=pts, labels=paste(pts, "%", sep=""))
legend(x="topright", legend=colnames(rates), lty=1, col=tsRainbow, cex=.64,bty="n")

#head(rates)
#tail(rates)

7 Bootstrap to Obtain Spot Rates

In this section spot rates are obtained via the commonly used bootstrapping technique.

spot <- bootstrap(.5, libor, swap)
libor <- merge(libor_1m/100, libor_3m/100)
spot <- merge(libor, spot)
spot <- na.omit(spot)
write.zoo(spot, file="spot.csv")
#tail(spot)

8 Calculate Principal Components

8.1 Eigenvectors and Eigenvalues

spotPctDiff <- spot/lag(spot,1) - 1
spotPctDiff <- unclass(spotPctDiff)[2:length(spotPctDiff[,1]),1:length(colnames(spotPctDiff))]
scaledSpotPctDiff <- scale(spotPctDiff)
pca <- eigen(cor(scaledSpotPctDiff))
write(pca$values, file="eigenvalues.csv", sep=",", ncolumns=1)
write(t(pca$vectors), file="eigenvectors.csv", sep=",", ncolumns=length(colnames(spot)))
vol.df <- data.frame(attr(scaledSpotPctDiff, "scaled:scale"))
colnames(vol.df)<-c("vol")
write(vol.df$vol, file="vol.csv", sep=",", ncolumns=1)

8.2 Examine Amount of Variation Explained

The plot shows that three principal components describe the bulk of the variation.

8.3 Principal Component Interpretation

In this section the sign of the key principal components are changed to align with interpretations in the paper by Folpmers and Eweg referenced above. There the authors say:

1. the first principal component is highly correlated with the rates of all maturities and the correlations
2. the second principal component is negatively correlated with short-maturity rates and positively correlated with long-maturity series. This allows for an interpretation of the second principal component as a factor that tilts the term structure (“steepness”);
3. the third principal component is positively correlated with short- and long-maturity rates and negatively correlated with intermediate-maturity rates, and is therefore a factor that effects the curvature.

8.4 Graph the 3 Main Principal Components

9 PCA Shocks and Shocked Rates

Principal component analysis based shocks are developed by considering the product of the principal component vectors (eigenvectors) and their standard deviations (square root of the eigenvalues). The degree of the stress is obtained by applying the desired number of standard deviation. The normal distribution is assumed for this purpose in combination with the 99.5th percentile.

Shocked rates are obtaind by applying the principal component analysis based shocks described above to the current interest rate term structure. Finally, the square root of 12 is applied to convert the monthly volatility to an annual figure.

# Develop shocks using the principal component vectors and variances
# Need to multiply eigenvectors by eigenvalues. First must get eigevectors in form to multiply easily.
eigenGoodForm <- NULL
stdevRateGoodForm <- NULL
lastRateGoodForm <- NULL
curvePoints <- ncol(scaledSpotPctDiff)
for(i in 1:curvePoints)
{
  eigenGoodForm <- c(eigenGoodForm, c(rep(pca$values[i], curvePoints)))  
  stdevRateGoodForm <- c(stdevRateGoodForm, c(attr(scaledSpotPctDiff, "scaled:scale")))
  lastRateGoodForm <- c(lastRateGoodForm, c((unclass(last(spot))[1,])))
}
# Develop the shocked pc's and the rateShocks
pcaShockUp <- pca$vectors * eigenGoodForm^.5 * qnorm(.995)
rateShockUp <- (1+pcaShockUp*stdevRateGoodForm *sqrt(12))* lastRateGoodForm

pcaShockUp <- subset(data.frame(pcaShockUp), select=-c(4:length(colnames(data.frame(pcaShockUp)))))
colnames(pcaShockUp) <- c("PCA_1_Up", "PCA_2_Up", "PCA_3_Up")
rateShockUp <-  subset(data.frame(rateShockUp), select=-c(4:length(colnames(data.frame(rateShockUp)))))
colnames(rateShockUp) <- c("rate_PCA_1_Up", "rate_PCA_2_Up", "rate_PCA_3_Up")

pcaShockDown <- pca$vectors * eigenGoodForm^.5 * -qnorm(.995)
rateShockDown <- (1+pcaShockDown*stdevRateGoodForm * sqrt(12))* lastRateGoodForm 

pcaShockDown <- subset(data.frame(pcaShockDown), select=-c(4:length(colnames(data.frame(pcaShockDown)))))
colnames(pcaShockDown) <- c("PCA_1_Down", "PCA_2_Down", "PCA_3_Down")
rateShockDown <-  subset(data.frame(rateShockDown), select=-c(4:length(colnames(data.frame(rateShockDown)))))
colnames(rateShockDown) <- c("rate_PCA_1_Down", "rate_PCA_2_Down", "rate_PCA_3_Down")

pcaShock <- cbind(pcaShockUp, pcaShockDown)
pcaShock$Term <- c(1/12, 3/12, seq(from=.5, t=30, by=.5))
col_id <- grep(c("Term"), names(pcaShock))
pcaShock <- pcaShock[,c(col_id, (1:ncol(pcaShock))[-col_id])]

rateShock <- cbind(rateShockUp, rateShockDown)
rateShock$Term <- c(1/12, 3/12, seq(from=.5, t=30, by=.5))
rateShock$Base <- unclass(last(spot))[1,]
col_id <- grep(c("Base"), names(rateShock))
rateShock <- rateShock[,c(col_id, (1:ncol(rateShock))[-col_id])]
col_id <- grep(c("Term"), names(rateShock))
rateShock <- rateShock[,c(col_id, (1:ncol(rateShock))[-col_id])]

9.1 Graph the Shocked Rates

10 Appendix

10.1 PCA Shocks

Term	PCA_1_Up	PCA_2_Up	PCA_3_Up	PCA_1_Down	PCA_2_Down	PCA_3_Down
0.0833333	0.5226558	-1.2805887	1.9347736	-0.5226558	1.2805887	-1.9347736
0.2500000	0.2131756	-1.8095608	1.7256696	-0.2131756	1.8095608	-1.7256696
0.5000000	0.4654727	-1.8177872	1.6238846	-0.4654727	1.8177872	-1.6238846
1.0000000	1.0797353	-2.0917086	0.4787185	-1.0797353	2.0917086	-0.4787185
1.5000000	1.5978651	-1.9186782	-0.0421185	-1.5978651	1.9186782	0.0421185
2.0000000	1.8307288	-1.6672682	-0.3634506	-1.8307288	1.6672682	0.3634506
2.5000000	2.0188586	-1.4550278	-0.5149783	-2.0188586	1.4550278	0.5149783
3.0000000	2.1235027	-1.2663640	-0.6178181	-2.1235027	1.2663640	0.6178181
3.5000000	2.2120751	-1.1228490	-0.6177914	-2.2120751	1.1228490	0.6177914
4.0000000	2.2723573	-0.9927522	-0.6182653	-2.2723573	0.9927522	0.6182653
4.5000000	2.3343312	-0.8829901	-0.5737572	-2.3343312	0.8829901	0.5737572
5.0000000	2.3813653	-0.7768722	-0.5358476	-2.3813653	0.7768722	0.5358476
5.5000000	2.4215763	-0.6905226	-0.4866376	-2.4215763	0.6905226	0.4866376
6.0000000	2.4543602	-0.6043465	-0.4412225	-2.4543602	0.6043465	0.4412225
6.5000000	2.4807173	-0.5188463	-0.3992134	-2.4807173	0.5188463	0.3992134
7.0000000	2.5014039	-0.4343785	-0.3602801	-2.5014039	0.4343785	0.3602801
7.5000000	2.5168052	-0.3778303	-0.3255626	-2.5168052	0.3778303	0.3255626
8.0000000	2.5298002	-0.3205099	-0.2917583	-2.5298002	0.3205099	0.2917583
8.5000000	2.5404996	-0.2625528	-0.2588334	-2.5404996	0.2625528	0.2588334
9.0000000	2.5489887	-0.2040787	-0.2267599	-2.5489887	0.2040787	0.2267599
9.5000000	2.5553348	-0.1451964	-0.1955149	-2.5553348	0.1451964	0.1955149
10.0000000	2.5595913	-0.0860065	-0.1650807	-2.5595913	0.0860065	0.1650807
10.5000000	2.5617060	-0.0697487	-0.1513998	-2.5617060	0.0697487	0.1513998
11.0000000	2.5636223	-0.0530663	-0.1376344	-2.5636223	0.0530663	0.1376344
11.5000000	2.5653324	-0.0359803	-0.1237816	-2.5653324	0.0359803	0.1237816
12.0000000	2.5668282	-0.0185080	-0.1098392	-2.5668282	0.0185080	0.1098392
12.5000000	2.5681015	-0.0006635	-0.0958050	-2.5681015	0.0006635	0.0958050
13.0000000	2.5691438	0.0175415	-0.0816776	-2.5691438	-0.0175415	0.0816776
13.5000000	2.5699462	0.0360971	-0.0674555	-2.5699462	-0.0360971	0.0674555
14.0000000	2.5704996	0.0549954	-0.0531379	-2.5704996	-0.0549954	0.0531379
14.5000000	2.5707947	0.0742295	-0.0387238	-2.5707947	-0.0742295	0.0387238
15.0000000	2.5708217	0.0937937	-0.0242129	-2.5708217	-0.0937937	0.0242129
15.5000000	2.5705705	0.1136833	-0.0096048	-2.5705705	-0.1136833	0.0096048
16.0000000	2.5700306	0.1338942	0.0051003	-2.5700306	-0.1338942	-0.0051003
16.5000000	2.5691914	0.1544228	0.0199021	-2.5691914	-0.1544228	-0.0199021
17.0000000	2.5680418	0.1752663	0.0348003	-2.5680418	-0.1752663	-0.0348003
17.5000000	2.5665702	0.1964220	0.0497939	-2.5665702	-0.1964220	-0.0497939
18.0000000	2.5647648	0.2178877	0.0648819	-2.5647648	-0.2178877	-0.0648819
18.5000000	2.5626135	0.2396612	0.0800630	-2.5626135	-0.2396612	-0.0800630
19.0000000	2.5601037	0.2617407	0.0953355	-2.5601037	-0.2617407	-0.0953355
19.5000000	2.5572224	0.2841244	0.1106976	-2.5572224	-0.2841244	-0.1106976
20.0000000	2.5539563	0.3068105	0.1261470	-2.5539563	-0.3068105	-0.1261470
20.5000000	2.5502917	0.3297973	0.1416813	-2.5502917	-0.3297973	-0.1416813
21.0000000	2.5462145	0.3530831	0.1572975	-2.5462145	-0.3530831	-0.1572975
21.5000000	2.5417102	0.3766661	0.1729926	-2.5417102	-0.3766661	-0.1729926
22.0000000	2.5367641	0.4005443	0.1887630	-2.5367641	-0.4005443	-0.1887630
22.5000000	2.5313607	0.4247159	0.2046048	-2.5313607	-0.4247159	-0.2046048
23.0000000	2.5254846	0.4491786	0.2205137	-2.5254846	-0.4491786	-0.2205137
23.5000000	2.5191197	0.4739301	0.2364852	-2.5191197	-0.4739301	-0.2364852
24.0000000	2.5122496	0.4989679	0.2525141	-2.5122496	-0.4989679	-0.2525141
24.5000000	2.5048576	0.5242893	0.2685950	-2.5048576	-0.5242893	-0.2685950
25.0000000	2.4969264	0.5498913	0.2847220	-2.4969264	-0.5498913	-0.2847220
25.5000000	2.4884387	0.5757708	0.3008887	-2.4884387	-0.5757708	-0.3008887
26.0000000	2.4793763	0.6019242	0.3170883	-2.4793763	-0.6019242	-0.3170883
26.5000000	2.4697211	0.6283479	0.3333134	-2.4697211	-0.6283479	-0.3333134
27.0000000	2.4594542	0.6550379	0.3495562	-2.4594542	-0.6550379	-0.3495562
27.5000000	2.4485566	0.6819900	0.3658085	-2.4485566	-0.6819900	-0.3658085
28.0000000	2.4370086	0.7091996	0.3820613	-2.4370086	-0.7091996	-0.3820613
28.5000000	2.4247902	0.7366620	0.3983051	-2.4247902	-0.7366620	-0.3983051
29.0000000	2.4118809	0.7643719	0.4145300	-2.4118809	-0.7643719	-0.4145300
29.5000000	2.3982597	0.7923242	0.4307254	-2.3982597	-0.7923242	-0.4307254
30.0000000	2.3839049	0.8205131	0.4468799	-2.3839049	-0.8205131	-0.4468799

10.2 Shocked Interest Rates

A factor of square root of 12 was applied to the monthly volatility. This makes these shifts appropriate for a 1 year horizon.

Term	Base	rate_PCA_1_Up	rate_PCA_2_Up	rate_PCA_3_Up	rate_PCA_1_Down	rate_PCA_2_Down	rate_PCA_3_Down
0.0833333	0.0018650	0.0022294	0.0009721	0.0032141	0.0015006	0.0027579	0.0005159
0.2500000	0.0028320	0.0030589	0.0009058	0.0046689	0.0026051	0.0047582	0.0009951
0.5000000	0.0044485	0.0051192	0.0018292	0.0067884	0.0037778	0.0070678	0.0021086
1.0000000	0.0053082	0.0072759	0.0014961	0.0061806	0.0033404	0.0091202	0.0044357
1.5000000	0.0073217	0.0112623	0.0025899	0.0072178	0.0033811	0.0120535	0.0074256
2.0000000	0.0093438	0.0154626	0.0037713	0.0081290	0.0032249	0.0149163	0.0105585
2.5000000	0.0111723	0.0192242	0.0053692	0.0091184	0.0031204	0.0169755	0.0132262
3.0000000	0.0130116	0.0230306	0.0070367	0.0100966	0.0029925	0.0189865	0.0159265
3.5000000	0.0145005	0.0257672	0.0087815	0.0113539	0.0032337	0.0202195	0.0176471
4.0000000	0.0159996	0.0285221	0.0105287	0.0125924	0.0034770	0.0214704	0.0194067
4.5000000	0.0171936	0.0303095	0.0122323	0.0139698	0.0040777	0.0221548	0.0204174
5.0000000	0.0183966	0.0321412	0.0139127	0.0153038	0.0046520	0.0228805	0.0214894
5.5000000	0.0192871	0.0332679	0.0153004	0.0164775	0.0053063	0.0232738	0.0220967
6.0000000	0.0201847	0.0344330	0.0166763	0.0176233	0.0059364	0.0236931	0.0227461
6.5000000	0.0210897	0.0356307	0.0180485	0.0187497	0.0065488	0.0241310	0.0234298
7.0000000	0.0220027	0.0368574	0.0194231	0.0198631	0.0071479	0.0245823	0.0241422
7.5000000	0.0225332	0.0374031	0.0203009	0.0206097	0.0076634	0.0247655	0.0244567
8.0000000	0.0230687	0.0379681	0.0211810	0.0213504	0.0081693	0.0249564	0.0247870
8.5000000	0.0236091	0.0385509	0.0220649	0.0220868	0.0086674	0.0251533	0.0251314
9.0000000	0.0241546	0.0391501	0.0229540	0.0228206	0.0091592	0.0253552	0.0254886
9.5000000	0.0247052	0.0397647	0.0238495	0.0235530	0.0096458	0.0255609	0.0258575
10.0000000	0.0252612	0.0403940	0.0247527	0.0242852	0.0101283	0.0257697	0.0262372
10.5000000	0.0253617	0.0404387	0.0249512	0.0244706	0.0102846	0.0257722	0.0262527
11.0000000	0.0254649	0.0404892	0.0251539	0.0246582	0.0104405	0.0257759	0.0262715
11.5000000	0.0255705	0.0405448	0.0253605	0.0248480	0.0105962	0.0257805	0.0262930
12.0000000	0.0256784	0.0406051	0.0255707	0.0250396	0.0107516	0.0257860	0.0263171
12.5000000	0.0257883	0.0406697	0.0257845	0.0252331	0.0109069	0.0257921	0.0263435
13.0000000	0.0259002	0.0407382	0.0260015	0.0254284	0.0110622	0.0257989	0.0263719
13.5000000	0.0260139	0.0408103	0.0262217	0.0256255	0.0112174	0.0258060	0.0264022
14.0000000	0.0261293	0.0408859	0.0264450	0.0258243	0.0113727	0.0258136	0.0264344
14.5000000	0.0262464	0.0409647	0.0266714	0.0260247	0.0115281	0.0258214	0.0264681
15.0000000	0.0263651	0.0410465	0.0269007	0.0262268	0.0116836	0.0258294	0.0265033
15.5000000	0.0264853	0.0411311	0.0271330	0.0264305	0.0118394	0.0258376	0.0265400
16.0000000	0.0266070	0.0412186	0.0273682	0.0266360	0.0119954	0.0258457	0.0265780
16.5000000	0.0267301	0.0413086	0.0276064	0.0268430	0.0121516	0.0258539	0.0266172
17.0000000	0.0268547	0.0414012	0.0278475	0.0270518	0.0123082	0.0258619	0.0266576
17.5000000	0.0269807	0.0414962	0.0280916	0.0272623	0.0124652	0.0258698	0.0266991
18.0000000	0.0271081	0.0415936	0.0283387	0.0274745	0.0126225	0.0258775	0.0267416
18.5000000	0.0272368	0.0416933	0.0285888	0.0276885	0.0127803	0.0258848	0.0267852
19.0000000	0.0273670	0.0417953	0.0288421	0.0279043	0.0129386	0.0258918	0.0268297
19.5000000	0.0274985	0.0418995	0.0290985	0.0281219	0.0130975	0.0258984	0.0268751
20.0000000	0.0276314	0.0420059	0.0293582	0.0283414	0.0132568	0.0259045	0.0269214
20.5000000	0.0277656	0.0421145	0.0296212	0.0285628	0.0134168	0.0259101	0.0269685
21.0000000	0.0279013	0.0422252	0.0298876	0.0287862	0.0135774	0.0259150	0.0270164
21.5000000	0.0280384	0.0423380	0.0301575	0.0290116	0.0137387	0.0259192	0.0270651
22.0000000	0.0281768	0.0424529	0.0304310	0.0292391	0.0139007	0.0259227	0.0271145
22.5000000	0.0283167	0.0425700	0.0307082	0.0294688	0.0140635	0.0259253	0.0271647
23.0000000	0.0284581	0.0426891	0.0309892	0.0297007	0.0142271	0.0259270	0.0272155
23.5000000	0.0286009	0.0428104	0.0312742	0.0299349	0.0143915	0.0259277	0.0272670
24.0000000	0.0287453	0.0429337	0.0315633	0.0301714	0.0145568	0.0259272	0.0273191
24.5000000	0.0288911	0.0430592	0.0318566	0.0304103	0.0147230	0.0259256	0.0273719
25.0000000	0.0290385	0.0431868	0.0321543	0.0306518	0.0148902	0.0259226	0.0274252
25.5000000	0.0291875	0.0433166	0.0324566	0.0308959	0.0150583	0.0259183	0.0274790
26.0000000	0.0293381	0.0434485	0.0327637	0.0311426	0.0152276	0.0259124	0.0275335
26.5000000	0.0294903	0.0435826	0.0330757	0.0313922	0.0153979	0.0259049	0.0275884
27.0000000	0.0296442	0.0437190	0.0333928	0.0316446	0.0155694	0.0258956	0.0276438
27.5000000	0.0297998	0.0438575	0.0337153	0.0319000	0.0157421	0.0258844	0.0276997
28.0000000	0.0299573	0.0439984	0.0340434	0.0321585	0.0159161	0.0258711	0.0277560
28.5000000	0.0301165	0.0441416	0.0343774	0.0324203	0.0160914	0.0258556	0.0278127
29.0000000	0.0302775	0.0442871	0.0347174	0.0326854	0.0162680	0.0258376	0.0278697
29.5000000	0.0304405	0.0444349	0.0350639	0.0329539	0.0164461	0.0258171	0.0279271
30.0000000	0.0306054	0.0445852	0.0354171	0.0332261	0.0166256	0.0257938	0.0279848

11 Helper Functions

bootstrap<-function(freq,libor, swap){
  at <- availTenors(colnames(swap))
  max <- at$tenor[length(colnames(swap))]
  result <- libor[,match("libor_6m", colnames(libor))]/100
  colnames(result)[length(colnames(result))] <- "spot_.5yr"
  df <- (1+result[,length(colnames(result))])^-.5
  colnames(df)[length(colnames(df))] <- "df_.5yr"
  for(i in seq(from=1, to=max, by=freq)){
    txt <- paste0("swap_", i, "yr")
    coupon <- swap[,match(txt, colnames(swap))]/100/2
    result$newcol <- ((1+coupon)/(1-annuityDF(df)*coupon))^(1/i)-1
    txt <- paste0("spot_", i, "yr")
    colnames(result)[length(colnames(result))] <- txt
    txt <- paste0("df_", i, "yr")
    df$newcol <- (1+result[,length(colnames(result))])^-i
    colnames(df)[length(colnames(df))] <- txt
  }
  result <- na.omit(result)
  return(result)
}

annuityDF<-function(df)
{
  max <- length(colnames(df))
  for(i in seq(from=1, to=max,by=1)){
    if(i==1){
      result <-  df[,1]
      colnames(result)[length(colnames(result))] <- "sum"
    }else{
      result <- result + df[,i]
    }
  }
  return(result)
}

# This function combines other functions to loop from the lowest desired swap
# rate to the highest and interpolates along the way. The intention of this function
# is to develop a swap rate for every half-year point from 1 - 30 years.
interpolateSwapsPeriods<-function(yrLow, yrHigh, freq, rates){
  result <- rates
  at <- availTenors(colnames(rates))
  for(i in seq(from=yrLow, to=yrHigh, by=freq)){
    if(i %in% at$tenor){
      #do nothing
    }
    else{
      txt <- paste0("swap_", i, "yr")
      result$newcol <- interpolateSwaps(i, rates)
      colnames(result)[length(colnames(result))] = txt
      result <- sortSwapCols(result)
    }
  }
  return(result)
}

# This function takes a year and a swap rate table and linearly interpolates across
# the whole time series
interpolateSwaps<-function(yr, rates){
  hl <- highlow(yr, colnames(rates))
  
  if(is.null(hl)){
    return(NULL)
  }
  else if(hl[1]==hl[2]){
    return(rates[,hl[3]])
  }
  else{
    return((rates[,hl[3]]*(hl[2]-yr)+rates[,hl[4]]*(yr-hl[1]))/(hl[2]-hl[1]))
  }
}

# Given a year and header of a swap rate table, this function returns a vector
# of length four where the first element is the known tenor just below the given
# year, the second elements is the know tenor just above the given year, and the
# third and fourth are the location of columns related to the values respectively.
highlow <- function(yr, availSwaps){
  max <- length(availSwaps)
  
  txtStart <- regexpr("_", availSwaps[1])
  txtEnd <- regexpr("yr", availSwaps[1])
  minVal <- as.numeric(substr(availSwaps[1], txtStart[1]+1, txtEnd[1]-1))
  
  txtStart <- regexpr("_", availSwaps[max])
  txtEnd <- regexpr("yr", availSwaps[max])
  maxVal <- as.numeric(substr(availSwaps[max], txtStart[1]+1, txtEnd[1]-1))
  
  if(yr<minVal || yr>maxVal){
    return(NULL)
  }
  
  for(i in 1:max){
    txtStart <- regexpr("_", availSwaps[i])
    txtEnd <- regexpr("yr", availSwaps[i])
    l <- as.numeric(substr(availSwaps[i], txtStart[1]+1, txtEnd[1]-1))
    txtStart <- regexpr("_", availSwaps[i+1])
    txtEnd <- regexpr("yr", availSwaps[i+1])
    h <- as.numeric(substr(availSwaps[i+1], txtStart[1]+1, txtEnd[1]-1))
    if(yr == l){
      return(c(c(l,l), c(i,i)))
    }
    else if(yr == h ){
      return(c(c(h,h), c(i+l, i+1)))
    }
    else if(yr > l && yr < h){
      return(c(c(l,h),c(i,i+1)))
    }
  }
}

# This function sorts the columns of a swap rate table by tenor as identified in the
# text of the column header
sortSwapCols<- function(swap){
  # Convert the time-series to data.frame (seems that sorting is easier)
  df <- data.frame(date=index(swap), coredata(swap))
  max <- length(colnames(swap))
  
  sorter <- availTenors(colnames(swap))
  sorter <- sorter[order(sorter$tenor),]
  
  for(i in 1:max){
    df <- movetolast(df, c(colnames(swap)[sorter$col[i]]))
  }
  result <- xts(subset(df, select=-date), order.by=df$date)
  return(result)
}

# This function is used by the sorting algorithm to move a column to the end
# of a data.frame
movetolast <- function(data, move) {
  data[c(setdiff(names(data), move), move)]
}

# This function reads the header of a swap rate table and returns an object
# with elements col and tenor both of which are vectors. They identify the
# known tenors and the column location of them in the swap rate table.
availTenors<- function(header){
  max <- length(header)
  tenors <- numeric(max)
  col <- c(1:max)
  for(i in 1:max)
  {
    txtStart <- regexpr("_", header[i])
    txtEnd <- regexpr("yr", header[i]) 
    tenors[i]<- as.numeric(substr(header[i], txtStart[1]+1, txtEnd[1]-1))
  }
  result <- data.frame(col=col)
  result$tenor <- tenors
  return(result)
}