cir.R

vasicek <- function(T, t=0, r=0.05, a=0.2, b=0.1, sigma=0.02) {
  sharpe = 0
  rbar = b + sigma*sharpe/a - 0.5*sigma^2/a^2
  
  A <- function(t, T) {
    exp(rbar*(B(t,T) + t - T) - B(t,T)^2*sigma^2/(4*a))
  }
  
  B <- function(t, T) {
    (1 - exp(-a*(T-t)))/a
  }
  
  P <- function(t, T, r) {
    A(t,T) * exp(-B(t,T) * r)
  }
  
  -log(P(t, T, r))/(T-t) # implied yield
}

cir <- function(T, t=0, r=0.05, a=0.2, b=0.1, sigma=0.02) {
  sigma = sigma/sqrt(r)
  barSharpe = 0
  gamma = sqrt((a + barSharpe)^2 + 2*sigma^2)
  
  A <- function(t, T) {
    ((2*gamma*exp( (a + barSharpe + gamma) * (T- t) / 2) )/
      ( (a + barSharpe + gamma) * (exp(gamma*(T-t))- 1) + 2*gamma ))^(2*a*b/(sigma^2))
  }
  
  B <- function(t, T) {
    2*(exp(gamma*(T-t)) - 1)/
      ((a + barSharpe + gamma) * (exp(gamma*(T-t)) - 1) + 2*gamma)
  }

  P <- function(t, T, r) {
    A(t,T) * exp(-B(t,T) * r) 
  }
  
  -log(P(t, T, r))/(T-t) # implied yield
  #P(t,T,r)
}

t = seq(0, 25, by=0.1)
sigmas = c(0.02, 0.10)

# Vasicek Model vs CIR Model
for(i in 1:length(sigmas)) {
  r1 <- numeric(0)
  r2 <- numeric(0)
  for(num in t) {
    r1 = c(r1, vasicek(num, sigma=sigmas[i])) # vasicek
    r2 = c(r2, cir(num, sigma=sigmas[i])) # cir
  }
  if(i==1) {
    plot(t, r1, type='l', lty=3,col='black', ylim=c(0.045, 0.10), ylab="r", xlab="T", main=paste0("Vasicek vs CIR, sigma = ", sigmas[i]))
    lines(t, r2, type='l', col='black')
    legend(1, 0.095, c("Vasicek", "CIR"), lty=1:2)
  } else {
    plot(t, r1, type='l', lty=3,col='black', ylim=c(0,0.09), ylab="r", xlab="T", main=paste0("Vasicek vs CIR, sigma = ", sigmas[i]))
    lines(t, r2, type='l', col='black')
    legend(1, 0.08, c("Vasicek", "CIR"), lty=1:2)
  }
}