# Tag Archives: Yield curve

## Calibrated Hull and White short-rates with RQuantLib and ESGtoolkit

In this post, I use R packages RQuantLib and ESGtoolkit for the calibration and simulation of the famous Hull and White short-rate model.

QuantLib is an open source C++ library for quantitative analysis, modeling, trading, and risk management of financial assets. RQuantLib is built upon it, providing R users with an interface to the library .

ESGtoolkit provides tools for building Economic Scenarios Generators (ESG) for Insurance. The package is primarily built for research purpose, and comes with no warranty. For an introduction to ESGtoolkit, you can read this slideshare, or this blog post. A development version of the package is available on Github with, I must admit, only 2 or 3 commits for now.

The Hull and White (1994) model was proposed to address Vasicek’s model poor fitting of the initial term structure of interest rates. The model is defined as:

$dr(t) = \left( \theta(t) - a r(t) \right) dt + \sigma dW(t)$

Where $a$ and $\sigma$ are positive constants, and $(W(t))_{t \geq 0}$ is a standard brownian motion under a risk-neutral probability. $t \mapsto \theta(t)$, which is constant in Vasicek’s model, is a function constructed so as to correctly match the initial term structure of interest rates.

An alternative and convenient representation of the model is:

$r(t) = x(t) + \alpha(t),$

where

$dx(t) = - a x(t) dt + \sigma dW(t),$

$x(0) = 0,$

$\alpha(t) = f^M(0, t) + \frac{\sigma^2}{2 a^2}(1 - e^{-at})^2$

and $f^M(0, t)$ are market-implied instantaneous forward rates for maturities $t \geq 0$.

In insurance market consistent pricing, the model is often calibrated to swaptions, as there are no market prices for  embedded options and guarantees.

Two parameters $a$ and $\sigma$ are surely not enough to get back the whole swaptions volatility surface (or even ATM swaptions). But a perfect calibration to market-quoted swaptions isn’t vital and may lead to unnecessary overfitting. A more complex model may fit more precisely the swaptions volatility surface, but could still be proved to be more wrong for the purpose: insurance liabilities do not even match exactly market swaptions’ characteristics.

It’s worth mentioning that the yield curve bootstrapping procedure currently implemented in RQuantLib, makes the implicit assumption that LIBOR is a good proxy for risk-free rates, and collateral doesn’t matter. These assumptions are still widely used in insurance, along with simple (parallel) adjustments for credit/liquidity risks, but were abandoned by markets after the 2007 subprime crisis.

For more details on the new multiple-curve approach, see for example http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2219548. In this paper, the authors introduce a multiple-curve bootstrapping procedure, available in QuantLib.

# Cleaning the workspace
rm(list=ls())

suppressPackageStartupMessages(library(RQuantLib))
suppressPackageStartupMessages(library(ESGtoolkit))

# Frequency of simulation and interpolation
freq <- "monthly"
delta_t <- 1/12

# This data is taken from sample code shipped with QuantLib 0.3.10.
settleDate=as.Date('2002-2-19'),
payFixed=TRUE,
dt=delta_t,
strike=.06,
method="HWAnalytic",
interpWhat="zero",
interpHow= "spline")

# Market data used to construct the term structure of interest rates
# Deposits and swaps
tsQuotes <- list(d1w =0.0382,
d1m =0.0372,
d3m = 0.0363,
d6m = 0.0353,
d9m = 0.0348,
d1y = 0.0345,
s2y = 0.037125,
s3y =0.0398,
s5y =0.0443,
s10y =0.05165,
s15y =0.055175)

# Swaption volatility matrix with corresponding maturities and tenors
swaptionMaturities <- c(1,2,3,4,5)
swapTenors <- c(1,2,3,4,5)
volMatrix <- matrix(
c(0.1490, 0.1340, 0.1228, 0.1189, 0.1148,
0.1290, 0.1201, 0.1146, 0.1108, 0.1040,
0.1149, 0.1112, 0.1070, 0.1010, 0.0957,
0.1047, 0.1021, 0.0980, 0.0951, 0.1270,
0.1000, 0.0950, 0.0900, 0.1230, 0.1160),
ncol=5, byrow=TRUE)

# Pricing the Bermudan swaptions
pricing <- RQuantLib::BermudanSwaption(params, tsQuotes,
swaptionMaturities, swapTenors, volMatrix)
summary(pricing)

# Constructing the spot term structure of interest rates
# based on input market data
times <- seq(from = delta_t, to = 5, by = delta_t)
curves <- RQuantLib::DiscountCurve(params, tsQuotes, times)
maturities <- curves$times marketzerorates <- curves$zerorates
marketprices <- curves$discounts ############# Hull-White short-rates simulaton # Horizon, number of simulations, frequency horizon <- 5 # I take horizon = 5 because of swaptions maturities nb.sims <- 10000 # Calibrated Hull-White parameters from RQuantLib a <- pricing$a
sigma <- pricing$sigma # Simulation of gaussian shocks with ESGtoolkit set.seed(4) eps <- ESGtoolkit::simshocks(n = nb.sims, horizon = horizon, frequency = freq) # Simulation of the factor x with ESGtoolkit x <- ESGtoolkit::simdiff(n = nb.sims, horizon = horizon, frequency = freq, model = "OU", x0 = 0, theta1 = 0, theta2 = a, theta3 = sigma, eps = eps) # I use RQuantlib's forward rates. With the low monthly frequency, # I consider them as being instantaneous forward rates fwdrates <- ts(replicate(nb.sims, curves$forwards),
start = start(x),
deltat = deltat(x))

# alpha
t.out <- seq(from = 0, to = horizon, by = delta_t)
param.alpha <- ts(replicate(nb.sims, 0.5*(sigma^2)*(1 - exp(-a*t.out))^2/(a^2)),
start = start(x), deltat = deltat(x))
alpha <- fwdrates + param.alpha

# The short-rate
r <- x + alpha

# Stochastic discount factors (numerical integration currently is very basic)
Dt <- ESGtoolkit::esgdiscountfactor(r = r, X = 1)

# Monte Carlo prices and zero rates deduced from stochastic discount factors
montecarloprices <- rowMeans(Dt)
montecarlozerorates <- -log(montecarloprices)/maturities # RQuantLib uses continuous compounding

# Confidence interval for the difference between market and monte carlo prices
conf.int <- t(apply((Dt - marketprices)[-1, ], 1, function(x) t.test(x)\$conf.int))

# Viz
par(mfrow = c(2, 2))
# short-rate quantiles
ESGtoolkit::esgplotbands(r, xlab = "maturities", ylab = "short-rate quantiles",
main = "short-rate quantiles")
# monte carlo vs market zero rates
plot(maturities, montecarlozerorates, type='l', col = 'blue', lwd = 3,
main = "monte carlo vs market \n zero rates")
points(maturities, marketzerorates, col = 'red')
# monte carlo vs market zero-coupon prices
plot(maturities, montecarloprices, type ='l', col = 'blue', lwd = 3,
main = "monte carlo vs market \n zero-coupon prices")
points(maturities, marketprices, col = 'red')
# confidence interval for the price difference
matplot(maturities[-1], conf.int, type = 'l',
main = "confidence interval \n for the price difference")



Filed under R stuff

## Solvency II yield curve ‘fits exactly’. Too tightly to explain ?

For the QIS5 and the LTGA (Quantitative Impact Studies), the Prudential Authority (EIOPA) suggested the use of Smith-Wilson method for the purpose of discounting future cash flows. A description of this method can be found here. The technical specifications for the QIS5 argued about this method that : “the term structure is fitted exactly to all observed zero coupon bond prices”. The method was subsequently used for the LTGA.

But : is fitting exactly, a desirable feature for a yield curve ? I won’t give a definitive answer here, but some hints, as the question could be answered in a lot of different ways (well… this post is mainly an excuse to present you ycinterextra in action after this, for you to tell me how to improve the package).

I’ll fit the Nelson-Siegel (NS), Svensson (SV) and Smith-Wilson (SW) models (implemented in ycinterextra) to the data from Diebold, F. X., & Li, C. (2006). Forecasting the term structure of government bond yields. Journal of econometrics, 130(2), 337-364. And  we’ll see how accurate are the  values they produce on out-of-sample data. I used mixed ideas from Diebold and Li (2006), along with Gilli, Manfred and Schumann, Enrico, A Note on ‘Good Starting Values’ in Numerical Optimisation (June 3, 2010). Available at SSRN: http://ssrn.com/abstract=1620083.

The data can be found here : http://www.ssc.upenn.edu/~fdiebold/papers/paper49/FBFITTED.txt.

# I put the data in a text file.
data.diebold.li <- scan("FBFITTED.txt", skip = 14)


library(ycinterextra)

# Number of observation dates from January 1985 through December 2000
nbdates <- 372

# Time to maturities for each observation date
mats <- round(c(1, 3, 6, 9, 12, 15, 18, 21, 24, 30, 36, 48, 60, 72, 84, 96, 108, 120)/12, 3)
nbmats <- length(mats)
nbcol <- nbmats + 1

# Formatting the data
data.diebold.li.final <- NULL
for(i in seq_len(nbdates))
{
data.diebold.li.final <- rbind(data.diebold.li.final, data.diebold.li[(nbcol*(i -1)+1):(nbcol*(i -1)+nbcol)])
}


A 3D plot of the observed data can be obtained with the following function :

x <- mats
y <- data.diebold.li.final[,1]
z <- data.diebold.li.final[,-1]

graph3D <- function(x, y, z)
{
par(bg = "white")
nrz <- nrow(z)
ncz <- ncol(z)
# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette

nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(y, x, z, col = color[facetcol], phi = 30, theta = -30, ticktype= 'detailed', nticks = 3,
xlab = "observation date", ylab = "maturity", zlab = "zero rates in pct.")
}

graph3D(x, y, z)


Now, let’s calibrate the 3 models, NS, SV and SW to the observed yield curves (this takes a lot of time, as there are 372 curves, each with 18 maturities) :

nrz <- nrow(z)
ncz <- ncol(z)
yM.NS <- matrix(NA, nrow = nrz, ncol = ncz)
yM.SV <- matrix(NA, nrow = nrz, ncol = ncz)
yM.SW <- matrix(NA, nrow = nrz, ncol = ncz)

# Loop over the dates from January 1985 through December 2000
for (t in seq_len(nbdates))
{
# market yields
yM <- as.numeric(data.diebold.li.final[t, -1])

yM.NS[t, ] <- fitted(ycinter(yM = yM, matsin = mats, matsout = mats,
method = "NS", typeres = "rates"))

yM.SV[t, ] <- fitted(ycinter(yM = yM, matsin = mats, matsout = mats,
method = "SV", typeres = "rates"))

yM.SW[t, ] <- fitted(ycinter(yM = yM, matsin = mats, matsout = mats,
method = "SW", typeres = "rates"))
}


After a little time… SW curve seems to fit very well indeed, and the plot below looks pretty similar to the one we made before :

graph3D(x, y, yM.SW)


With a plot of the average yield curve (average of the 372 observed curves as in Diebold and Li (2006)), we can see that all of the three models seem to provide fairly good fits to the data :

# Average yield curves
plot(mats, apply(z, 2, mean), pch = 3, xlab = "time to maturity",
ylab = "yield to maturity",
main = paste("Actual and fitted (model-based)", "\n", "average yield curves"))
lines(mats, apply(yM.NS, 2, mean), col = "blue", lwd= 2)
lines(mats, apply(yM.SV, 2, mean), col = "red", lwd= 2)
lines(mats, apply(yM.SW, 2, mean), col = "green", lwd=2, lty = 2)
legend("bottomright", c("Actual", "NS", "SV", "SW"),
col=c("black", "blue", "red", "green"), text.col = c("black", "blue", "red", "green"),
lty=c(NA, 1, 1, 2), lwd = c(1, 2, 2, 2), pch=c(3, NA, NA, NA))


Here are plots of the residuals for NS and SW (for SV, the plot is pretty similar to the plot for NS) :

# residuals
graph3D(x, y, yM.NS-z)
graph3D(x, y, yM.SW-z)


Now, taking a closer look at some given cross-sections (here on 3/31/1989), we see what’s actually going on (click for a better resolution) :

t <- 231

par(mfrow=c(2,2))

plot(mats, data.diebold.li.final[t, -1], xlab = "time to maturity",
ylab = "yield to maturity", col="red", main = "Observed data")
points(mats, data.diebold.li.final[t, -1], col="red")

plot(mats, yM.NS [t,], type='l', xlab = "time to maturity",
ylab = "yield to maturity", col = "blue", main = "Nelson-Siegel fit")
points(mats, data.diebold.li.final[t, -1], col="red")

plot(mats, yM.SV[t,], type='l', xlab = "time to maturity",
ylab = "yield to maturity", col = "blue", main = "Svensson fit")
points(mats, data.diebold.li.final[t, -1], col="red")

plot(mats, yM.SW[t,], type='l', xlab = "time to maturity",
ylab = "yield to maturity", col = "blue", main = "Smith-Wilson fit")
points(mats, data.diebold.li.final[t, -1], col="red")



And for 7/31/1989 (with t == 235) :

NS fits, well. SV fits a litte bit better, due to its extended number of parameters. SW fits perfectly to each point, due to its laaarge number of parameters.

Let’s see how accurate are the values produced by the 3 methods on out-of-sample values, by the use of a simple machine learning test set method (a  method like k-fold cross-validation can be used as well, it’s computationnally more expensive, but simple though). I divide the sample data into 2 sets :

•  A learning set with the zero rates observed at time to maturities : c(1, 6, 9, 15, 18, 21, 24, 30, 48, 72, 84, 96, 120)/12
• A test set with the zero rates observed at time to maturities : c(3, 12, 36, 60, 108)/12

The models’ parameters are estimated using the learning set, and the test set is used to assess how well (with less error) the model can predict unknown values :

# Maturities for test set
matsoutsample <- c(3, 12, 36, 60, 108)/12
matchoutsample <- match(matsoutsample, mats)
m <- length(matsoutsample)

# Maturities for learning set
matsinsample <- mats[-matchoutsample]


# Residuals matrices
res.NS <- matrix(NA, nrow = nrz, ncol = length(matchoutsample))
res.SV <- matrix(NA, nrow = nrz, ncol = length(matchoutsample))
res.SW <- matrix(NA, nrow = nrz, ncol = length(matchoutsample))

# Loop over the dates from January 1985 through December 2000
for (t in seq_len(nbdates))
{
yM.t <- as.numeric(data.diebold.li.final[t, -c(1, matchoutsample)])

yM.t.bechm <- as.numeric(data.diebold.li.final[t, -1])

res.NS[t, ] <- (fitted(ycinter(yM = yM.t, matsin = matsinsample, matsout = mats,
method = "NS", typeres = "rates")) - yM.t.bechm)[matchoutsample]

res.SV[t, ] <- (fitted(ycinter(yM = yM.t, matsin = matsinsample, matsout = mats,
method = "SV", typeres = "rates")) - yM.t.bechm)[matchoutsample]

res.SW[t, ] <- (fitted(ycinter(yM = yM.t, matsin = matsinsample, matsout = mats,
method = "SW", typeres = "rates")) - yM.t.bechm)[matchoutsample]
}


For each model, we get the 372 mean squared errors as follows :

mse.ns <- apply((res.NS), 1, function(x) crossprod(x))/m
mse.sv <- apply((res.SV), 1, function(x) crossprod(x))/m
mse.sw <- apply((res.SW), 1, function(x) crossprod(x))/m


A summary of these mean squared errors shows that SW actually overfits, and fails to predict the out-of-sample data :

summary(mse.ns)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0005895 0.0054050 0.0117700 0.0227100 0.0198400 0.5957000
summary(mse.sv)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0003694 0.0075310 0.0154200 0.0329400 0.0319700 0.6079000
summary(mse.sw)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.057 0.556 1.478 431.200 10.820 29040.000