Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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 below.

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}
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://kiandlee.blogspot.com#——————————————————–## Nelson-Siegel fitted yield curve using rjags#========================================================# graphics.off()  # clear all graphsrm(list = ls()) # remove all files from your workspace library(rjags)library(coda) # Prepare datamat <– c(3, 6, 9, 12, 24, 36, 48, 60, 72,         84, 96, 108, 120, 144, 180, 240) yield <– c(5.638781, 5.857775, 6.062339, 6.232739,            6.533613, 6.586592, 6.604274, 6.672092,            6.643833, 6.636791, 6.621937, 6.608467,            6.608286, 6.631803, 6.640675, 6.569184) vnm <– c(“beta0”,“beta1”,“beta2”,“la”,“sigma”) Data <– list(x = mat, y = yield) # for use inside modelData$y_short <– yieldData$y_long  <– yieldData$y_mid <– yield # Model stringjags.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 parameterx11(); plot(mod1.samples[][,c(1,2,3)])x11(); plot(mod1.samples[][,c(4,5)]) # descriptive statistics sm <– summary(mod1.samples) # correlationcor(mod1.samples[]) # means of posteriors of parametersb0 <– sm$statistics[vnm, “Mean”]b1 <– sm$statistics[vnm, “Mean”]b2 <– sm$statistics[vnm, “Mean”]la <– sm\$statistics[vnm, “Mean”]c(b0, b1, b2, la) # in-sample fitypred <– b0 +         b1*((1–exp(–la*mat))/(la*mat)) +         b2*((1–exp(–la*mat))/(la*mat) – exp(–la*mat)) # graphx11(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) Colored by Color Scripter cs