Introduction
A common practice in order to obtain the discount factor for different maturities, given a spot term structure is to use linear interpolation. Another common practice, often used by central banks, is to calibrate a model, such as the Nelson-Siegel model, on the term structure, in this way they ensure to obtain a function with they can obtain the interest rates and so the discount factor for every maturities. In this article is presented a modified version of the original Nelson-Siegel model, in this version were introduced 2 parameters \(\beta_3\) and \(\tau_2\), in order to have more flexibility. This modified version is called Nelson-Siegel-Svensson model and the equation, with parameters \((\beta_0,\beta_1, \beta_2, \beta_3, \tau_1, \tau_2)\), is given by:
\[\hat{h}(0,T) = \beta_0 + \beta_1 f^1(T, \tau_1) + \beta_2 f^2(T, \tau_1) + \beta_3 f^2(T, \tau_2)\] where:
\[f^1(T,\tau_1) = \frac{1 - e^{-\frac{T}{\tau_1}}}{\frac{T}{\tau_1}}\] and
\[f^2(T,\tau) = \frac{1 - e^{-\frac{T}{\tau}}}{\frac{T}{\tau}} - e^{-\frac{T}{\tau}}\]
A benefit of the Nelson-Siegel-Svensson model is the interpretations of \(\beta_0\) and \(\beta_0 + \beta_1\). In fact \(\beta_0\) can be interpreted as the long term interest rate, while \(\beta_0 + \beta_1\) as the instantaneous short term interest rate (Yallup, 2011). In order to have this contraints valid we will need to perform a contraint optimization problem in which the function to optimize is given by the mean square error of the fitted value with respect to the real values (over discrete maturities \((T_a, ....T_b)\)):
\[\underset{\beta_0, \beta_1, \beta_2, \beta_3, \tau_1, \tau_2}{argmin}\biggl\{ \sum_{i = a}^{i = b} \bigl[ h(0,T_i) - \hat{h}(0,T_i) \bigl]^2 \biggl\}\] Subject to the contraint: \[\underset{\beta_0, \beta_1}{argmin}\bigl\{ | \beta_0 + beta_1 - h(0,T)| \bigl\} \underline{<} \alpha\] where \(\alpha = 0.01\). is a threshold level for the accuracy of the contraint.
In order to find the optimal parameters, is needed to solve a minimization problem starting from a set of 6 parameters, that in general, are initializated randomly or a priori. In this case, the following alghorithm was implemented:
- fix \(\beta_0\) and \(\beta_1\) both equal to \(\frac{h(0,t_1)}{2}\), in such a way the initial parameters respect the contraint.
- generate randomly the other 4 parameters in a reasonable bound \([-5, 5]\).
- perform the minimization problem subject the contraint.
Since the other 4 parameters are randomly generated, in order to find the best set of six parameters that minimizes the mean square error between the real values and the fitted ones, is common to repeat this operation many times. In this way the minimization problem is repeated many times and in the end one is able to choose the best parameters.
Nelson-Siegel base function
Takes as input a vector of maturities and a set of 6 parameters in this order \((\beta_0,\beta_1, \beta_2, \beta_3, \tau_1, \tau_2)\), and gives in output a value for the estimates spot interest rates for that maturities.
nelson_siegel <- function(maturity, params){
# parameters
beta0 = params[1]
beta1 = params[2]
beta2 = params[3]
beta3 = params[4]
tau1 = params[5]
tau2 = params[6]
# function
beta0 +
beta1 * (1 - exp(-maturity/tau1)) / (maturity/tau1) +
beta2 * ((1 - exp(-maturity/tau1)) / (maturity/tau1) - exp(-maturity/tau1)) +
beta3 * ((1 - exp(-maturity/tau2)) / (maturity/tau2) - exp(-maturity/tau2))
}
# example
params = c(0.78, -1.32, -1.72, 3.61, 0.21, 7.16)
maturity = c(0.5, 1, 2, 3, 4, 5)Minimization Problem
Given a set of initial parameters our aim is to find the best set of parameters such that the loss function is minumum and the contraint on \(\beta_0 + \beta_1\) is respected with an error given by the parameter contraint_threshold, set by default equal to 0.01. The function will return a tibble with 1 row and 5 columns:
- term_structure: a list, containing the initial term structure.
- optim_params: a list, containing the optimal parameters.
- start_params: a list, containing the starting parameters.
- fit: a list, containing the fitted_values.
- mse: a number, the mean square error of the fit.
optim_nelson_siegel <- function(term_structure = NULL, maturity = NULL, params = rep(0.1, 6), contraint_threshold = 0.01){
loss_function <- function(params){
# parameters
beta0 = params[1]
beta1 = params[2]
beta2 = params[3]
beta3 = params[4]
tau1 = params[5]
tau2 = params[6]
# contraints
contraint = abs(beta0 + beta1 - term_structure[1])
if(contraint > contraint_threshold ){
return(NA)
}
sum(
(term_structure - beta0 -
beta1 * (1 - exp(-maturity/tau1)) / (maturity/tau1) -
beta2 * ((1 - exp(-maturity/tau1)) / (maturity/tau1) - exp(-maturity/tau1)) -
beta3 * ((1 - exp(-maturity/tau2)) / (maturity/tau2) - exp(-maturity/tau2))
)^2)
}
# set initial values for to beta0 and beta1 that respect the contraint
params = c(term_structure[1]/2, term_structure[1]/2, params[3:6] )
# minimization of the loss function
optimization = optim(params, loss_function)
# optimal parameters
optim_params = c(beta0 = optimization$par[1],
beta1 = optimization$par[2],
beta2 = optimization$par[3],
beta3 = optimization$par[4],
tau1 = optimization$par[5],
tau2 = optimization$par[6])
# fitted value from Nelson-Siegel base function
fitted_values = nelson_siegel(maturity, optim_params)
# compute the mean square error
mse_fit = sd(fitted_values - term_structure, na.rm = TRUE)
# output
dplyr::tibble(
maturity = list(maturity),
term_structure = list(term_structure),
start_params = list(params),
optim_params = list(optim_params),
fit = list(fitted_values),
mse = mse_fit
)
} Simulation for Optimal Parameters
Since the starting parameters are arbitrary and we do not have particulary contraints except the one specified before (except the fact that we do not want extreme huge values for the parameters), in order to find the best vector of 6 parameters, such that the mean square error of the model is minimum we can:
- generate randomly in an interval \((\beta_2, \beta_3, \tau_1, \tau_2)\) and set \(\beta_0 = \beta_1 = \frac{h(t, T_1)}{2}\).
- perform the minimization problem with the function optim.nelson_siegel and compute the mse.
- repeat again for n-times and choose the parameters that gives a mean square error minimum.
# n.params: number of parameter to generate
# params.min: for random generation, minimum parameter
# params.max: for random generation, maximum parameter
# seed: to control randomness
random_params <- function(n.params = 1, params.min = 0, params.max = 5, seed = 1){
set.seed(seed)
runif(n.params, params.min, params.max)
}
# object: an output of the function "optim.nelson_siegel"
# n: number of simulations
# params.min: for random generation, minimum parameter
# params.max: for random generation, maximum parameter
# verbose: disply progress in the importation.
calibrate_nelson_siegel <- function(object, n = 100, params.min = 0, params.max = 5, contraint_threshold = 0.01, verbose = TRUE ){
# initialize the parameters
term_structure = object$term_structure[[1]]
maturity = object$maturity[[1]]
# list containing all the simulations
simulations = list()
# safe version to avoid errors if we made many simulations
safe_optim = purrr::safely(optim_nelson_siegel)
for(i in 1:n){
# generate a random seed
random_seed = mean(random_params(10, 0, 100000, seed = i))
# generate random parameters
random_params = random_params(6, params.min, params.max, seed = random_seed)
simulations[[i]] = safe_optim(term_structure = term_structure, maturity = maturity, params = random_params, contraint_threshold = contraint_threshold)$result
if( verbose &(i%%50 == 0)){message("Simulations: ", i, "/", n)}
}
# unique dataset for all the simulation
simulations_df = dplyr::bind_rows(simulations)
# add the initial object
simulations_df = dplyr::bind_rows(object, simulations)
# index for the simulations
simulations_df = dplyr::mutate(simulations_df, n = 1:nrow(simulations_df))
return(simulations_df)
}A unique function
We can put the functions before constructed in an unique function that easily perform n simulations for the parameters and gives as output the best result, a plot of the fitted
# 4) optimal_params_ns
#
# term_structure = NULL
# maturity = NULL
# n = 1000
# params.init = rep(0.1, 6)
# params.min = -5
# params.max = 5
# contraint_threshold = 0.01
# label = NULL (label for the plot)
# verbose = TRUE
optimal_params_ns <- function(term_structure = NULL, maturity = NULL,
n = 1000, params.init = rep(0.1, 6), params.min = -5, params.max = 5,
contraint_threshold = 0.01, label = NULL, verbose = TRUE){
# first fit
first_fit_ns = optim_nelson_siegel(term_structure = term_structure, maturity = maturity, params = params.init)
# simulations
sim_fit_ns = calibrate_nelson_siegel(first_fit_ns, n = n, params.min = params.min, params.max = params.max, contraint_threshold = contraint_threshold)
# best parameters
df_optim_params = sim_fit_ns[which(sim_fit_ns$mse == min(sim_fit_ns$mse, na.rm = TRUE)),]
# setting the title of the plot
if(!is.null(label) & is.character(label)){
plot_title = paste0("Fitted Nelson-Siegel-Svensonn vs Real Value ", "(", label, ")" )
} else {
plot_title = "Fitted Nelson-Siegel-Svensonn vs Real Value"
}
# plot of fitted vs real values
plot_df = dplyr::inner_join(
dplyr::tibble(
t = maturity,
pred = df_optim_params$fit[[1]]
),
dplyr::tibble(
t = maturity,
real = term_structure
),
by = "t"
)
# Plot Real value vs Fitted Values
plot_ns = plot_df %>%
mutate(label = paste0("T = ", round(t, 3)))%>%
ggplot()+
geom_point(aes(t, pred), color = "red", size = 2, alpha = 0.8) +
geom_point(aes(t, real), color = "black", alpha = 0.5)+
geom_line(aes(t, pred), color = "red", size = 1) +
geom_line(aes(t, real), color = "black", linetype = "dashed") +
geom_label(aes(t+1, real-0.001, label = label), size = 1.5)+
ggthemes::theme_solarized()+
theme(axis.text.x = element_text(angle = 25, face = "bold", size = 7),
axis.text.y = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic"),
plot.caption = element_text(face = "italic"),
panel.grid.major.x = element_line(colour="grey60", linetype="dotted"),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(colour="grey60", linetype="dotted"),
legend.text = element_text(face = "italic", size = 5),
legend.title = element_text(face = "bold"),
legend.position = "top" ) +
scale_x_continuous(breaks=c(0, 0.5, 1,2,3,4,5,7,10,15, 20, 25, 30))+
ggtitle(plot_title, subtitle = "Fitted Value in Red and Real Values in Black" )+
xlab("Maturities")+
ylab("") +
labs(caption = paste0("Mean Squared Error of the Fit: ", round(df_optim_params$mse, 6)))
# output
structure(
list(
optim_params = tibble( beta0=df_optim_params$optim_params[[1]][1],
beta1=df_optim_params$optim_params[[1]][2],
beta2=df_optim_params$optim_params[[1]][3],
beta3=df_optim_params$optim_params[[1]][4],
tau1=df_optim_params$optim_params[[1]][5],
tau2=df_optim_params$optim_params[[1]][6]),
plot_ns = plot_ns,
simulations = sim_fit_ns,
df_optim = df_optim_params
)
)
}Example of Fit
# calibration interest rates
interest_rate = ts_discount$h0T # term structure interest rates
maturity = ts_discount$t # maturities
# calibration using monte carlo simulation
interest_rate_ns = optim_nelson_siegel(term_structure = interest_rate, maturity = maturity, params = rep(0.1, 6))
optim_interest_rate_ns = calibrate_nelson_siegel(interest_rate_ns, n = 1000, params.min = -5, params.max = 5, verbose = FALSE)
# find the best parameters
optim_params_interest_rates = optim_interest_rate_ns[which(optim_interest_rate_ns$mse == min(optim_interest_rate_ns$mse, na.rm = TRUE)),]
# fitted interest
fitted_interest_rate = nelson_siegel(maturity, optim_params_interest_rates$optim_params[[1]])
# compare the real values with the fitted ones
plot_NSS = inner_join(
tibble(
t = maturity,
pred = fitted_interest_rate
),
tibble(
t = maturity,
real = ts_discount$h0T
), by = "t"
) %>%
mutate(label = paste0("T = ", round(t, 3)))%>%
ggplot()+
geom_point(aes(t, pred), color = "red", size = 2, alpha = 0.8) +
geom_point(aes(t, real), color = "black", alpha = 0.5)+
geom_line(aes(t, pred), color = "red", size = 1) +
geom_line(aes(t, real), color = "black", linetype = "dashed") +
geom_label(aes(t+1, real-0.001, label = label), size = 1.5)+
ggthemes::theme_solarized()+
theme(axis.text.x = element_text(angle = 25, face = "bold", size = 7),
axis.text.y = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic"),
plot.caption = element_text(face = "italic"),
panel.grid.major.x = element_line(colour="grey60", linetype="dotted"),
panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_line(colour="grey60", linetype="dotted"),
legend.text = element_text(face = "italic", size = 5),
legend.title = element_text(face = "bold"),
legend.position = "top" ) +
scale_x_continuous(breaks=c(0, 0.5, 1,2,3,4,5,7,10,15, 20, 25, 30))+
ggtitle("Nelson-Siegel-Svensonn model")+
xlab("Maturities")+
ylab("Interest Rates") +
labs(caption = "Fitted Values vs Real Values")
plot_NSS