Bayesian Estimation of Nelson-Siegel model using rjags R package

This post shows how to use rjags R package to estimate Nelson-Siegel yield curve model based using a Bayesian MCMC in a compact way.


Nelson-Siegel model and rjags R package



We installed and used rjags R package and performed a Bayesian estimation of multiple linear regression as an example in the previous post

In an interesting and practical example, let's run a Bayesian estimation of the Nelson-Siegel yield curve model.


Nelson-Siegel model


The Nelson-Siegel model is widely used in practice for fitting the cross-sectional term structure of interest rates in a parsimonious way. It is so famous and lots of extended models have been emerged based on the this model. It assumes that an instantaneous forward rate is of the following form.

\[\begin{align} f(\tau) = \beta_0 + \beta_1e^{- \tau \lambda } + \beta_2 \tau \lambda e^{- \tau \lambda } \end{align}\]
Averaging through an integration of \(f(\tau)\) results in a continuously compounded spot rate.

\[\begin{align} y(\tau) &= \frac{1}{\tau} \int_{0}^{\tau} f(u) du \\ &= \beta_0 + \beta_1 \left( \frac{1-e^{- \tau \lambda }}{\tau \lambda }\right) + \beta_2 \left(\frac{1-e^{- \tau \lambda }}{\tau \lambda }-e^{- \tau \lambda }\right) \end{align}\]
Here, \(\tau\) is a maturity and \(y(\tau)\) is a continuously compounded spot rate. \(\beta_0, \beta_1, \beta_2\) are coefficients and \(\lambda\) is the shape or time decay parameter.


Likelihood


The likelihood of the Nelson-Siegel model is constructed from the following model specification.

\[\begin{align} Y_t &\sim \text{Normal}(\mu, \psi) \\ \mu &= \beta_0 + \beta_1 \left( \frac{1-e^{- \tau \lambda }}{\tau \lambda }\right) + \beta_2 \left(\frac{1-e^{- \tau \lambda }}{\tau \lambda }-e^{- \tau \lambda }\right) \\ \sigma &= 1/\sqrt{\psi} \end{align}\]
Here, \(\psi\) in \(\text{Normal}(\mu, \psi)\) is not a variance but a precision which is the reciprocal of a variance (\(\psi = \frac{1}{\sigma^2}\)).


Prior


\(\beta_0, \beta_1, \beta_2\) are assumed to have normal priors but \(\psi\) is non-negative and is assumed to have a gamma prior. \(\lambda\) is also non-negative and is assumed to have a uniform prior.

\[\begin{align} \beta_0 &\sim \text{Normal}(y(\tau_{max}), 0.01) \\ \beta_1 &\sim \text{Normal}(y(\tau_{min})-y(\tau_{max}), 0.01) \\ \beta_2 &\sim \text{Normal}(2y(\tau_{mid})-y(\tau_{min})-y(\tau_{max}), 0.01) \\ \lambda &\sim \text{Unit}(0.02,0.15) \\ \psi &\sim \text{Gamma}(0.1,0.001) \end{align}\]

R code


The R code below performs Bayesian estimation of Nelson-Siegel model using a cross-section of yields with the help of rjags R package.

#========================================================#
# Quantitative ALM, Financial Econometrics & Derivatives 
# ML/DL using R, Python, Tensorflow by Sang-Heon Lee 
#
# https://shleeai.blogspot.com
#--------------------------------------------------------#
# Nelson-Siegel fitted yield curve using rjags
#========================================================#
 
graphics.off()  # clear all graphs
rm(list = ls()) # remove all files from your workspace
 
library(rjags)
library(coda)
 
# Prepare data
mat <- c(369122436486072,
         8496108120144180240)
 
yield <- c(5.6387815.8577756.0623396.232739
           6.5336136.5865926.6042746.672092
           6.6438336.6367916.6219376.608467
           6.6082866.6318036.6406756.569184)
 
vnm <- c("beta0","beta1","beta2","la","sigma")
 
Data <- list(x = mat, y = yield)
 
# for use inside model
Data$y_short <- yield[1]
Data$y_long  <- yield[16]
Data$y_mid   <- yield[6]
 
# Model string
jags.script <- "
model{
  # priors
  beta0 ~ dnorm(y_long,0.01)
  beta1 ~ dnorm(y_short-y_long,0.01)
  beta2 ~ dnorm(2*y_mid-y_short-y_long,0.01)  
  la    ~ dunif(0.02,0.15)
  psi   ~ dgamma(0.1,0.001)
  
  # transform
  sigma <- 1 / sqrt(psi)
  
  # likelihood
  for( i in 1:length(x[])) {
  
    y[i] ~ dnorm(mu[i], psi)
    mu[i] <- beta0 + beta1*((1-exp(-la*x[i]))/(la*x[i]))+
             beta2*((1-exp(-la*x[i]))/(la*x[i])-exp(-la*x[i]))
  }
}
" 
 
mod1 <- jags.model(textConnection(jags.script), data = Data, 
                   n.chains = 4, n.adapt = 1000)
 
update(mod1, 4000)
 
mod1.samples <- coda.samples(model = mod1, 
                             variable.names = vnm, 
                             n.iter = 2000)
 
# plot trace and posterior density for each parameter
x11(); plot(mod1.samples[][,c(1,2,3)])
x11(); plot(mod1.samples[][,c(4,5)])
 
# descriptive statistics 
sm <- summary(mod1.samples)  
 
# correlation
cor(mod1.samples[[1]])
 
# means of posteriors of parameters
b0 <- sm$statistics[vnm[1], "Mean"]
b1 <- sm$statistics[vnm[2], "Mean"]
b2 <- sm$statistics[vnm[3], "Mean"]
la <- sm$statistics[vnm[4], "Mean"]
c(b0, b1, b2, la)
 
# in-sample fit
ypred <- b0 +
         b1*((1-exp(-la*mat))/(la*mat)) +
         b2*((1-exp(-la*mat))/(la*mat) - exp(-la*mat))
 
# graph
x11(width = 16/2, height = 9/2); 
plot(mat, yield, col = "red", lwd=3, pch=16
     ylab = "Yield(%)",xlab = "Maturity(month)")
lines(mat, ypred, col = "blue", lwd = 3)
 
cs


Running this R code produces the following estimation results.

Nelson-Siegel model by using rjags R Package


Estimation results also contain of the following posterior distributions for each parameters.

Nelson-Siegel model by using rjags R Package
Nelson-Siegel model by using rjags R Package

We can generate the fitted yield curve with average coefficients.

Nelson-Siegel model by using rjags R Package


Concluding Remarks


This post shows how to use rjags R package to perform Bayesian MCMC estimation of Nelson-Siegel model. Of course, as my setting for priors is not an answer but an alternative proposal, I think more elaborations on setting priors are needed. \(\blacksquare\)


No comments:

Post a Comment