diff --git a/blackScholes.py b/blackScholes.py new file mode 100644 index 0000000..1462cd5 --- /dev/null +++ b/blackScholes.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*- +""" + +""" + + +from math import ceil +#from mpl_toolkits.mplot3d import axes3d +import matplotlib.pyplot as plt +import numpy as np +import implementSOR +import processIO + +def put_Black_Scholes(X=10, N=100, M=150, T=30, std=0.3, r=0.02, X_to_Smax = 0.5): + """ + European put options + T must be measured in days + r must be an annual rate of return + X_to_Smax is the ratio of X to SMax, must obey 0 < x_to_smax <= 1 + X_to_Smax determines how many zero values are in b initially + """ + k = T / M + A = set_A(N, k, std, r) #A stays constant throughout + alignment = set_alignment(X, k, r, std) #add this constant to b[0] at every timestep + b = init_b(X, N, X_to_Smax, alignment) #first value of b at timestep M-1 + # print(b) + + #Prepare stock price, timestep axes (x and y axes in the wireframe) + S_axis = [] + T_axis = [] + S_current_row = [] + for j in range(1, N): + S_current_row.append((X * j) / (X_to_Smax * N)) + for k in range(1, M): + T_range = [] + for l in range(1, N): + T_range.append(M - k) + T_axis.append(np.array(T_range)) + S_axis.append(np.array(S_current_row)) +# print(S_axis) +# print(T_axis) + + """ + Calculate option values using SOR + The option value at this timestep is set to b for next timestep + This is the z-axis in the wireframe + """ + Value_axis = [] + for i in range(0, M-1): + b = implementSOR.sparse_sor(A, b) + Value_axis.append(b) #add calculated value pre-alignment + b[0]+=alignment #this is vector b for the next timestep +# print(Value_axis) + +#plot this in matplotlib as a 3d wireframe + fig1 = plt.figure() + fig1.suptitle('Put Option for X='+str(X)+', N='+str(N)+', M='+str(M)+', T='+str(T)+', std='+str(std)+', r='+str(r)) + wire = fig1.add_subplot(111, projection='3d') + wire.plot_wireframe(S_axis, T_axis, Value_axis) + + wire.set_xlabel('Share price') + wire.set_ylabel('Time steps up until maturity') + wire.set_zlabel('Value of option') + plt.savefig('black_scholes_merton.png') + plt.show() + +def set_A(N, k, std, r): + A=[] + for i in range(0, N-1): + A_row = [] + for j in range(0, N-1): + if i==j: + A_row.append(1 + (k * r) + (k * (std**2) * ((i+1)**2))) + elif i==j+1: + A_row.append(-0.5 * (i + 1) * k * (((i + 1) * (std**2 )) - r)) + elif i==j-1: + A_row.append(-0.5 * (i + 1) * k * (((i + 1) * (std**2 )) + r)) + else: + A_row.append(0) + A.append(A_row) +# print(to_csr(A)) + return processIO.to_csr_format(A) + + +def set_alignment(X, k, r, std): + return 0.5 * k * X * ((std**2) - r) + + +def init_b(X, N, X_to_Smax, alignment): + b = [] + N_threshold = ceil(N * X_to_Smax) #value of N where S = X + for i in range(0, N_threshold): + b.append(X - ((i + 1) * X / N_threshold)) + b[0]+=alignment + for j in range(N_threshold, N-1): + b.append(0) + return b + + + +#put_Black_Scholes() #low value of k = T/M +put_Black_Scholes(M=30) #high value of k = T/M \ No newline at end of file diff --git a/implementSOR.py b/implementSOR.py index f504c07..f5b345f 100644 --- a/implementSOR.py +++ b/implementSOR.py @@ -65,25 +65,25 @@ def chk_converge(A, xnew, b, xold, x_seq_tol, res_tol, flag): -def sparse_sor(matrix_a, matrix_b, matrix_x, dimension_n, max_it, \ - x_tol, res_tol, w ): +def sparse_sor(matrix_a, vector_b, matrix_x, dimension_n, max_it=50, \ + x_tol=1e-13, res_tol=1e-13, w=1.25 ): num_it = 1 stop_reason = 6 #something has gone wrong if this does not get overwritten later matrix_x_last_it = np.array([0.0]) matrix_x_new = np.array(matrix_x) matrix_a_np = np.array(matrix_a) - matrix_b_np = np.array(matrix_b) + vector_b_np = np.array(vector_b) flag = True #required to enter while loop first time only while num_it <= max_it and \ - chk_converge(matrix_a_np, matrix_x_new, matrix_b_np, matrix_x_last_it, + chk_converge(matrix_a_np, matrix_x_new, vector_b_np, matrix_x_last_it, x_tol, res_tol, flag) == -1: flag = False matrix_x_last_it = np.array(matrix_x_new[:]) - for i in range(0,len(matrix_b)): + for i in range(0,len(vector_b)): sum = 0 first_nz = matrix_a_np[1][(matrix_a_np[2][i]-1)] - 1 #pos of first non zero on row for j in range(matrix_a_np[2][i]-1, matrix_a_np[2][i+1]-1): @@ -91,9 +91,9 @@ def sparse_sor(matrix_a, matrix_b, matrix_x, dimension_n, max_it, \ (matrix_a_np[2][i]-1) + first_nz] if matrix_a_np[1][j] == i+1: d = matrix_a[0][j] - matrix_x_new[i] = matrix_x_new[i] + w * (matrix_b_np[i] - sum) / d + matrix_x_new[i] = matrix_x_new[i] + w * (vector_b_np[i] - sum) / d num_it+=1 - conv = chk_converge(matrix_a_np, matrix_x_new, matrix_b_np, matrix_x_last_it, \ + conv = chk_converge(matrix_a_np, matrix_x_new, vector_b_np, matrix_x_last_it, \ x_tol, res_tol, False) if num_it-1 == max_it: stop_reason = 3 @@ -103,7 +103,7 @@ def sparse_sor(matrix_a, matrix_b, matrix_x, dimension_n, max_it, \ # set_output(stop_reason, max_it, num_it-1, x_tol, res_tol,matrix_x_new, ) return (stop_reason, max_it, num_it-1, matrix_x_new ) -def dense_SOR(matrix_a,matrix_b,dimension_n,max_iteration,w,matrix_x): +def dense_SOR(matrix_a,vector_b,dimension_n,max_iteration,w,matrix_x): print("Dense SOR Calculation") it_counter=0 while(it_counter<=max_iteration): @@ -111,7 +111,7 @@ def dense_SOR(matrix_a,matrix_b,dimension_n,max_iteration,w,matrix_x): sum=0.0 for col_counter in range(dimension_n): sum = sum + matrix_a[row_counter,col_counter]*matrix_x[col_counter] - matrix_x[row_counter]=matrix_x[row_counter]+w*(matrix_b[row_counter]-sum) \ + matrix_x[row_counter]=matrix_x[row_counter]+w*(vector_b[row_counter]-sum) \ /matrix_a[row_counter,row_counter] # print("Iteration: ",it_counter,"\n","X:",matrix_x) it_counter+=1 diff --git a/mainSOR.py b/mainSOR.py index 32c107d..9436722 100644 --- a/mainSOR.py +++ b/mainSOR.py @@ -21,9 +21,9 @@ def main(): try: #Input Processing from input file if(len(sys.argv) >= 2): - [dimension_n,matrix_a,matrix_b] = processIO.input(sys.argv[1]) + [dimension_n,matrix_a,vector_b] = processIO.input(sys.argv[1]) else: - [dimension_n,matrix_a,matrix_b]=processIO.input() + [dimension_n,matrix_a,vector_b]=processIO.input() #Constants max_it = 50 @@ -50,15 +50,19 @@ def main(): if((not validateMatrix.col_diagonally_dominant(matrix_a,dimension_n)) and \ (not validateMatrix.row_diagonally_dominant(matrix_a,dimension_n))): raise Exception("Matrix A wont converge as it is not row / column diagonally dominant.",6) - + + if(not validateMatrix.spectral_radius_convergence_check(matrix_a)): + raise Exception("Spectral radius check failed for Matrix A.",6) + + #Dense SOR Calculation -# implementSOR.dense_SOR(matrix_a,matrix_b,dimension_n,max_it,w,matrix_x) +# implementSOR.dense_SOR(matrix_a,vector_b,dimension_n,max_it,w,matrix_x) #Sparse SOR Calculation csr_matrix_a = processIO.to_csr_format(matrix_a) [stop_reason, max_it, num_it, matrix_x ] = implementSOR.sparse_sor(csr_matrix_a, \ - matrix_b, matrix_x, dimension_n, max_it, x_tol, res_tol, w ) + vector_b, matrix_x, dimension_n, max_it, x_tol, res_tol, w ) processIO.output(stop_reason, max_it, num_it, x_tol, res_tol, matrix_x,\ output_filename) diff --git a/processIO.py b/processIO.py index e02bb60..8f4dd1e 100644 --- a/processIO.py +++ b/processIO.py @@ -17,9 +17,9 @@ def input(file_name="nas_Sor.in"): if(matrix_a.shape[0]!=dimension_n or matrix_a.shape[1]!=dimension_n): raise Exception('Incorrect File Format for A Matrix dimensions') - matrix_b=np.genfromtxt(file_name,skip_header=dimension_n+1) -# print("Matrix B :",matrix_b, "Shape:",matrix_b.shape) - if(matrix_b.shape[0]!=dimension_n): + vector_b=np.genfromtxt(file_name,skip_header=dimension_n+1) +# print("Matrix B :",vector_b, "Shape:",vector_b.shape) + if(vector_b.shape[0]!=dimension_n): raise Exception('Incorrect File Format for B Matrix Dimensions') except OSError: @@ -34,7 +34,7 @@ def input(file_name="nas_Sor.in"): "Last Line should contain matrix B with n rows") sys.exit() - return (dimension_n,matrix_a,matrix_b) + return (dimension_n,matrix_a,vector_b) def output(stop_reason, max_it, num_it, x_tol, res_tol, matrix_x, file_name, \ stop_cause='something has gone wrong'): diff --git a/validateMatrix.py b/validateMatrix.py index 6d9700e..7092b55 100644 --- a/validateMatrix.py +++ b/validateMatrix.py @@ -4,6 +4,7 @@ """ import numpy as np + def non_zero_diagonal_check(matrix_a,dimension_n): for counter in range(dimension_n): if(matrix_a[counter,counter]==0.0): @@ -36,4 +37,25 @@ def col_diagonally_dominant(matrix_a,dimension_n): # print("Col Sum: ",sum,"Diag Value:",matrix_a[col_counter,col_counter]) if (abs(matrix_a[col_counter,col_counter])> sum): return True - return False \ No newline at end of file + return False + +def spectral_radius_convergence_check(matrix_a): + + matrix_l=np.tril(matrix_a, k=-1) + matrix_u=np.triu(matrix_a, k=-1) + matrix_d=np.diag(np.diag(matrix_a)) + +# print(matrix_l, matrix_u, matrix_d) + + +# print(np.linalg.inv(matrix_d+matrix_l) * -1) + + matrix_c = np.dot(-np.linalg.inv(matrix_d+matrix_l), matrix_u) +# print(matrix_c) + + radius=np.linalg.eigvals(matrix_c).max() +# print(radius) + + if radius < 1: + return False + return True \ No newline at end of file