Added blackScholas code

and spectral radius convergence check with some code refactoring

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

implementSOR.py

@@ -65,35 +65,35 @@ 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):
                 sum = sum + matrix_a_np[0][j] * matrix_x_new[j - \
                             (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,15 +103,15 @@ 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):
         for row_counter in range(dimension_n):
             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

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)           

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'):

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
+    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