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.



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 .

Simulated Paths

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 parameters
r0 <- 0.03
theta <- 0.10
k <- 0.3
beta <- 0.03
 
## simulate short rate paths
n <- 10    # MC simulation trials
T <- 10    # total time
m <- 200   # subintervals
dt <- T/m  # difference in time each subinterval
 
r <- matrix(0,m+1,n)  # matrix to hold short rate paths
r[1,] <- r0
 
for(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 paths
t <- 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)

Bond Pricing and Yield Curve

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 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.



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 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 . 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 model
VasicekZCBprice <- 
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 curves
theta <- 0.10
k <- 0.5
beta <- 0.03
 
r0 <- 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)

Formula vs Monte Carlo

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.9614 MC Price: 0.9623 MC Standard Error: 0.0005

## define model parameters
r0 <- 0.03
theta <- 0.10
k <- 0.3
beta <- 0.03
 
## simulate short rate paths
n <- 1000  # MC simulation trials
T <- 1     # total time
m <- 200   # subintervals
dt <- T/m  # difference in time each subinterval
 
r <- matrix(0,m+1,n)  # matrix to hold short rate paths
r[1,] <- r0
for(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 solution
ss <- colSums(r[2:(m+1),]*dt)  # integral estimate
c <- 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')