"""
Created on Wed Oct 01 21:58:13 2014
@author: mylurian
"""
import pylab as plt #Imported a library
import numpy as np #Imported a library
u = input('Enter u: ') #User defined up factor
d = input('Enter d: ') #User defined down factor
R = input('Enter R: ') #User defined Discount factor
G = input('Enter S: ') #User defined Original Stock Price
S = G #Duplicates the OSP for computational purposes
K = input('Enter K: ') #User defined Strike Price
T = input('Enter T: ') #User defined time to maturity
a = input('Enter a: ') #User defined u/d factor (per 10.10)
n = input('Enter n: ') #Total number of Monte Carlo iterations
s=[] #Creates an emptyset to be populated by calculated stock values at time T
count = 0 #Sets the initial count value for the loop Process 'count'
for _ in range(n*T): #Sets the total number of loops to processed
if count == T:
S = G #Resets the underlying stock value to the OSP after T loops
if count == T:
count = 0 #Resets the count to 0 after T loops
p = (R - d*S**-a)/((u*S**a)-(d*S**-a)) #Calculates the risk neutral probability
v = np.random.random(1) #Creates a vector with T elements of random numbers between 0 and 1
x = v #Duplicates the vector v for computational purposes
x[x <= p] = 1 #Sets the value of all random numbers in v to 1 if they are within the risk neutral probability
x[x < 1] = 0 #Sets the value of all random numbers in v to 0 if they are outside the risk neutral probability
X = x[0] #Defines X to call the 'zero'th element of vector x
def f(X):
return (S**(1 - a + (2*a*X)))*(u**X)*(d**(1-X)) #Defines the Underlying Stock Value at time t = count
S = f(X) #Duplicates the Underlying Stock value at time t = count
count = count + 1 #Increases the count 'count' by 1
if count == T:
s.append(S) #Populates the vector s with the Underlying Stock Values at time t = T
k = np.ones(n) * K #Creates a vector with n elements, all = K, for computational purposes
b = s - k #Creates a vector populated by all the contingent claims
b[b <= 0] = 0 #Sets all negative values of our contingient claims to 0
L = sum(s)/n #Averages the stock values at time T
A = sum(b)/n #Averages the Call Options at time T
B = (R**-T)*sum(b)/n #Returns the Call Option Premium of our OSP (per 10.10)
print 'Your Underlying Stock Price estimate is: ', L
print 'Your Contingient Claim estimate is: ', A
print 'Your Call Option Premium estimate is: ', B
if a == 0 and T == 10: #Allows the graphic info to be called only when a = 0 AND T = 10(the graph doesn't work otherwise)
fig, ax = plt.subplots() #Defines the figure and name of the plots to follow (otherwise it draws new graphs)
H = sorted(s) #Sorts the stock price at time T vector numerically for graphing purposes
H.insert(0,0) #Inserts a 0 in the '0'th element of vector H for graphing purposes
H.append(400) #Appends a 400 in last element of vector H for graphing purposes
J = b-G #Reduces our contingient claims by the OSP for graphing purposes
I = sorted(J) #Sorts the OSP reduced contingent claims numerically for graphing purposes
I.insert(0,-K) #Inserts a -K in the '0'th element of vector I for graphing purposes
I.append(400-2*K) #Appends a 400-2*K in the last element of vector I for graphing purposes
plt.plot(H, I, '-') #Plots vectors H and I with a blue line
ax.set_xlim(0, 400) #Sets the range of the x-axis as 0 to 400
ax.set_ylim(-150, 300) #sets the range of the y-axis as -150 to 300
plt.xlabel('Underlying Stock Price') #Labels the x-axis
plt.ylabel('Present Value') #Labels the y-axis
ax.plot(np.arange(400),np.zeros(400), '-r') #Graphs the present value = 0 with a red line
ax.plot(K,0, '^g') #Plots the Strike Price on the x-axis with a green triangle
C = [] #Creates an emptyset to be populated by present values of our portfolio
N = [] #Creates an emptyset to be populated by stock prices at time = T
chocula = 0 #Sets the initial count value for the Loop Process 'chocula'
for _ in range (n): #Sets the total number of loops to be processed
if chocula == 0:
N.append(s[chocula]) #Populates the vector N with the 'chocula'th element of the vector s
C.append(b[chocula]-G) #Populates the vector H with the 'chucula'th element of the vector b minus our OSP
points, = ax.plot(N,C, '*y')#Defines our original ordered pairs to be graphed (Our contingient claims)
plt.title('Monte Carlo Simulation of a Call Option') #Prints a title on our graph
text = plt.text(50, 200, 'n = 1') #Defines the original text for our n counter
plex = plt.text(50, 275, 'Your Underlying Stock Price estimate is: ') #Defines original text
glex = plt.text(50, 250, 'Your Contingient Claim estimate is: ') #Defines original text
blex = plt.text(50, 225, 'Your Call Option Premium estimate is: ') #Defines original text
if chocula < n:
N.append(s[chocula]) #Appends the 'chocula'th element of the vector x to the vector N
C.append(b[chocula]-G) #Appends the 'chocula'th element of the vector b - our OSP to the vector C
ax.plot(sum(N)/chocula,0, 'xm') #Plots the re-estimated values of our stock price at time T on the x-axis
points.set_data(N, C) #Replaces the original ordered pairs to be graphed with 'Chocula'th + 1 points of data
chocula = chocula + 1 #Increases the count 'chocula' by 1
text.set_text('n = {0}'.format(chocula)) #Updates the n counter with each loop
plt.pause(0.0001) #Artificially defines the time interval between individual plotted points
if chocula == n-1:
ax.plot(sum(N)/n,0, 'Dk') #Plots a black diamond on our final estimate of our Underlying Stock Price
plex.set_text('Your Underlying Stock Price estimate is: {0}'.format(L)) #Updates text after the final loop
glex.set_text('Your Contingient Claim estimate is: {0}'.format(A)) #Updates text after the final loop
blex.set_text('Your Call Option Premium estimate is: {0}'.format(B)) #Updates text after the final loop
ax.plot(sum(N)/n,-G+A, 'Dk') #Plots a black diamond on our final estimate of our Contingient Claim