"""ir.py: Script that compares various Interest Rate term structure models.
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 .
"""
__author__ = "Alexander Baker"
__copyright__ = "Copyright 2012, The Hatch Project"
__credits__ = ["Alexander Baker"]
__license__ = "GPL"
__version__ = "1.0.1"
__maintainer__ = "Alexander Baker"
__email__ = "baker.alexander@gmail.com"
__status__ = "Production"
import numpy as np
from numpy.random import standard_normal
from numpy import array, zeros, sqrt, shape
from pylab import plot, title, xlabel, ylabel, show, figure, grid
r0 = 0.5 # current UK funding rate
theta = 1.00 # 1 % long term interest rate
k = 0.3
beta = 0.03
## simulate short rate paths
n = 10000 # MC simulation trials
T = 24. # total time
m = 100. # subintervals
dt = T/m # difference in time each subinterval
r = np.zeros(shape=(n, m), dtype=float) # matrix to hold short rate paths
fig1 = figure(1)
fig1.clf()
ax1a = fig1.add_subplot(121)
for j in np.arange(0,n): # number of MC simulations
r[j,0] = r0
for i in np.arange(1,m): #trials per simulation
r[j,i] = r[j,i-1] + k*(theta-r[j,i-1])*dt + beta*sqrt(dt)*standard_normal();
#print j,i, r[j,i]
ax1a.plot(np.arange(0, T, dt), r[j])
## plot paths
t = np.arange(0, T, dt)
rT_expected = theta + (r0-theta)*pow(np.e,-k*t)
rT_stdev = sqrt( pow(beta,2)/(2*k)*(1-pow(np.e,-2*k*t)))
print 'expected', rT_expected, 'std', rT_stdev
plot(np.arange(0, T, dt), rT_expected, '-or')
plot(np.arange(0, T, dt), rT_expected+2*rT_stdev, '-ob')
plot(np.arange(0, T, dt), rT_expected-2*rT_stdev, '-ob')
print shape(t), shape(r)
title('Simulations %d Steps %d r0 %.2f alpha %.2f beta %.2f sigma %.2f' % (int(n), int(m), r0, k, theta, beta))
xlabel('steps')
ylabel('short term interest rate')
grid()
stock_ending_price = r[:,-1]
print 'average ending rate ', np.mean(stock_ending_price)
ax2a = fig1.add_subplot(122)
ax2a.hist(stock_ending_price, 100)
ax2a.grid()
show()
#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)
