[This article was first published on stotastic » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
A common model used in the financial industry for modelling the short rate (think overnight rate, but actually an infinitesimally short amount of time) is the Vasicek model. Although it is unlikely to perfectly fit the yield curve, it has some nice properties that make it a good model to work with. The dynamics of the Vasicek model are describe below.
dt+%5Cbeta%20dW_t%5E%5Cmathbb%7BQ%7D&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (1)")
In this model, the parameters 
are constants, and the random motion is generated by the Q measure Brownian motion 
. An important property of the Vasicek model is that the interest rate is mean reverting to 
, and the tendency to revert is controlled by 
. Also, this process is a diffusion process, hence Markovian, which will lead to some nice closed form formulas. Finally, the future value of the interest rate is normally distributed with the distribution e%5E%7B-%5Ckappa%20T%7D,%20%5C:%20%5Cfrac%7B%5Cbeta%5E2%7D%7B2%5Ckappa%7D(1-e%5E%7B-2%5Ckappa%20T%7D)%20%5Cright)&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (6)")
.
To give you an idea of what this process looks like, I have generate some sample paths using the Euler discretization method.
One important things you will notice is that this process starts at 0.03, but is pulled towards 0.10, which is the value of . I’ve added dashed lines to highlight this, as well as lines for the expected value of the process and confidence bands representing two standard deviations. Also, you may notice the major problem with this process is that the interest rate can drop below 0%. In the real world, this shouldn’t happen. The Cox-Ross-Ingersoll model, or CIR model, actually corrects for this, but the process is no longer normally distributed. The following is the R code used to generate the chart above.
## Simulate Sample Paths #### define model parametersr0 <- 0.03theta <- 0.10k <- 0.3beta <- 0.03## simulate short rate pathsn <- 10 # MC simulation trialsT <- 10 # total timem <- 200 # subintervalsdt <- T/m # difference in time each subintervalr <- matrix(0,m+1,n) # matrix to hold short rate pathsr[1,] <- r0for(j in 1:n){ for(i in 2:(m+1)){ dr <- k*(theta-r[i-1,j])*dt + beta*sqrt(dt)*rnorm(1,0,1) r[i,j] <- r[i-1,j] + dr }} ## plot pathst <- seq(0, T, dt)rT.expected <- theta + (r0-theta)*exp(-k*t)rT.stdev <- sqrt( beta^2/(2*k)*(1-exp(-2*k*t)))matplot(t, r[,1:10], type="l", lty=1, main="Short Rate Paths", ylab="rt") abline(h=theta, col="red", lty=2)lines(t, rT.expected, lty=2) lines(t, rT.expected + 2*rT.stdev, lty=2) lines(t, rT.expected - 2*rT.stdev, lty=2) points(0,r0)One can show that a zero coupon bond with a maturity at time T can be found by calculating the following expectation under the risk neutral measure =%5Cmathbb%7BE%7D%5E%5Cmathbb%7BQ%7D%20%5Cleft%5B%20e%5E%7B-%5Cint_0%5ET%7Br_u%20du%7D%7D%20%5Cright%5D&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (9)")
where the short rate process 
can be pretty much any process. We could estimate this expectation using Monte Carlo simulation, but the Vasicek model allows us to calculate the value of the zero coupon bond by using the Markov property and PDE techniques. Using this method, we can price the bond by the following equations.
=e%5E%7B-a(T)-b(T)r_0%7D%20%5C%5C%20%5Ctext%20%7Bwhere%7D%20%5C%5C%20b(%5Ctau)=%5Cfrac%7B1-e%5E%7B-%5Ctau%20%5Ckappa%7D%7D%7B%5Ckappa%7D%20%5C%5C%20a(%5Ctau)=(%5Ctheta%20-%20%5Cfrac%7B%5Cbeta%5E2%7D%7B2%5Ckappa%7D)%20(%5Ctau-b(%5Ctau))+%5Cfrac%7B%5Cbeta%5E2%7D%7B4%5Ckappa%7Db(%5Ctau)%5E2%20&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (11)")
This closed form solution for a zero coupon bond makes our lives much easier since we don’t need to compute the expectation under the martingale measure to find the price of a bond. Additionally, it will allow us to easily calculate the yield curve implied by the model. If we note that =e%5E%7B-y_T%20T%7D&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (12)")
where 
is the continuous time yield for a zero coupon bond with maturity T, then we can use our model estimate of the bond price to rearrange this formula and solve for the yield by )%7D%7BT%7D=y_T&bg=ffffff&%23038;fg=000000&%23038;s=0 "Fun with the Vasicek Interest Rate Model (14)")
. Finally, to demonstrate some of the possible curve shapes allowed by the Vasicek model, I produced the following chart. Notice that the yield curve can be both increasing or decreasing with maturity.
The following R code implements the closed form pricing function for a bond under the Vasicek model, and then uses the pricing function to generate the chart above.
## function to find ZCB price using Vasicek modelVasicekZCBprice <- function(r0, k, theta, beta, T){ b.vas <- (1/k)*(1-exp(-T*k)) a.vas <- (theta-beta^2/(2*k^2))*(T-b.vas)+(beta^2)/(4*k)*b.vas^2 return(exp(-a.vas-b.vas*r0))}## define model parameters for plotting yield curvestheta <- 0.10k <- 0.5beta <- 0.03r0 <- seq(0.00, 0.20, 0.05)n <- length(r0)yield <- matrix(0, 10, n)for(i in 1:n){ for(T in 1:10){ yield[T,i] <- -log(VasicekZCBprice(r0[i], k, theta, beta, T))/T }}maturity <- seq(1, 10, 1)matplot(maturity, yield, type="l", col="black", lty=1, main="Yield Curve Shapes")abline(h=theta, col="red", lty=2)Just to make sure the bond pricing formula was implemented correctly, I compared the price using the formula versus the price using Monte Carlo simulation to estimate the expectation under the martingale measure. Below are the results and R code. Seems everything is in order, although the Euler discretization method seems to be causing some error in the Monte Carlo results.
Exact Vasicek Price: 0.9614MC Price: 0.9623MC Standard Error: 0.0005
## define model parametersr0 <- 0.03theta <- 0.10k <- 0.3beta <- 0.03## simulate short rate pathsn <- 1000 # MC simulation trialsT <- 1 # total timem <- 200 # subintervalsdt <- T/m # difference in time each subintervalr <- matrix(0,m+1,n) # matrix to hold short rate pathsr[1,] <- r0for(j in 1:n){ for(i in 2:(m+1)){ dr <- k*(theta-r[i-1,j])*dt + beta*sqrt(dt)*rnorm(1,0,1) r[i,j] <- r[i-1,j] + dr }} ## calculate Monte Carlo bond price and compare to Exact Vasicek solutionss <- colSums(r[2:(m+1),]*dt) # integral estimatec <- exp(-ss)estimate <- mean(c)stdErr <- sd(c)/sqrt(n)exact <- VasicekZCBprice(r0, k, theta, beta, T)cat('Exact Vasicek Price:', round(exact,4), 'n')cat('MC Price:', round(estimate,4), 'n')cat('MC Standard Error:', round(stdErr,4), 'n')Related
To leave a comment for the author, please follow the link and comment on their blog: stotastic » R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
