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