"""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)