""" Created on Wed Oct 01 21:58:13 2014 @author: mylurian """ from matplotlib import pyplot as PLT 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 = PLT.figure() ax = fig.add_subplot(211) #Defines the figure and name of the plots to follow (otherwise it draws new graphs) ax1 = fig.add_subplot(212) 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 ax1.plot(H, I, '-') #Plots vectors H and I with a blue line ax1.set_xlim(0, 400) #Sets the range of the x-axis as 0 to 400 ax1.set_ylim(-150, 300) #sets the range of the y-axis as -150 to 300 ax.set_xlim(0,n) ax.set_ylim(0,500) plt.xlabel('Underlying Stock Price') #Labels the x-axis plt.ylabel('Present Value') #Labels the y-axis ax1.plot(np.arange(400),np.zeros(400), '-r') #Graphs the present value = 0 with a red line ax1.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, = ax1.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(chocula,sum(N)/chocula, '.m') #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: ax1.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