# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Evaluating option prices using the FFT
# ### Calculating call and put options for the Carr & Madan method

import numpy as np
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(12,10))
import pandas as pd
from scipy import fftpack
import itertools

import timeit

import numpy as np
import math
import scipy.stats as si
from math import log, sqrt, exp
from scipy import stats
import matplotlib.pyplot as plt


# ## Probability density function function

def ft_pdf(rf,sg,t,g):
    mu = rf-sg*sg/2
    #res = exp(i*t*mu*g-0.5*t*sg*sg.*g.^2)
    res=[np.e**(complex(0,t*mu*i)-0.5*t*sg*sg*i**(2)) for i in g]
    return res


# ## Characteristic function

def ft_option(rf,sg,t,a,g):
    #g_=[g - i*(a+1)]
    g_=[i - complex(0,(a+1)) for i in g]
    
    #y1 = ft_pdf(rf,sg,t,g_)
    y1 = ft_pdf(rf,sg,t,g_)
    
    #y2 = a^2 + a - g.^2 + i*(2*a+1)*g;
    y2=[a**(2)+a-i**(2)+complex(0,2*a+1)*i for i in g]
    
    # = y1./ y2;
    y=[ i/j for i,j in zip(y1,y2)]
    return y


# +

def FFT_Pricing_CallPut(S,rf,sg,t,alpha,npower):
    #find the number of points
    
    N = 2**(npower)
    
    #spacing in the strike
    lambda_ = 0.01;
    
    #parameters for the characteristic function grid (0,N*eta)
    eta = 2*np.pi/(N*lambda_)
    
    #create grid for the cf 
    u=list(range(N))
    g=[i * eta for i in u]
     
    #create grid for the strike
    k0 = (-N*lambda_)/2
    k = [k0+(i * lambda_) for i in u]
    
    #compute the characteristic function
    fc=ft_option(rf,sg,t,alpha,g)
    
    
    #trapezoidal rule
    w = [0.5]+ list(np.ones(N-2)) +[0.5]; 
    #h2=[a * b for a, b in zip ([ (np.e**(i*complex(0,-k0))) for i in g],[j*i * eta for i,j in zip(fc,w)])]
    
    h2=[(np.e**(complex(0,-k0*i)))*j*k*eta for i,j,k in itertools.zip_longest(g,fc,w)]
        
    #Fourier inversion
    g=fftpack.fft(h2)
    g2=[((np.e**(-alpha*a))/np.pi*b).real for a, b in zip (k,g)]
    
    K = [S*np.e**(i) for i in k] 
    K =np.transpose(K)


    Y=[np.e**(-rf*t)*S*i for i in g2]
    Y =np.transpose(Y)
    Dict = {'Strike Price':K, 'Option Price':Y}

    return Dict 


# -

W = FFT_Pricing_CallPut(100,0.05,0.2,1,7,8) #List of call prices from FFT
Strk = W['Strike Price']
Y = W['Option Price']






WPa = FFT_Pricing_CallPut(100,0.05,0.2,1,2.1,8)
StrkP = W['Strike Price']
YPa = WPa['Option Price']

# +


plt.plot(Strk,Y,label='BS-FFT')
#plt.plot(Strk,Y,label='BS-FFt')
plt.xlabel('Strike Price')
plt.ylabel('European call option')
plt.title('Black-Scholes (FFT) call price vs Strike price')
plt.legend()
plt.grid(True) 
#return df.head(10)
# -

# ## Changing the sign of $\alpha$ gives us the put options.



WP = FFT_Pricing_CallPut(100,0.05,0.2,1,-3.1,8)
StrkP = WP['Strike Price']
YP = WP['Option Price']


# +
plt.plot(StrkP,YP,label='BS-FFT')

plt.xlabel('Strike Price')
plt.ylabel('European Put option')
plt.title('Black-Scholes(FFT) put price vs Strike price')
plt.legend()
plt.grid(True) 
#return df.head(10)


# -



# # Analytical Black-Scholes Call option code

# +

# List of strikes generated from the FFT code
Strk = W['Strike Price']

"""
S = 100 
T = 1
r = 0.05
alpha = 2
sigma = 0.2
"""

call_list = [] # Call prices calculated using Black-Scholes analytical formula are appended to this list

def BScall(S, K, T, r, sigma):

    float(S)
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    value = (S * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))
    #si.norm.cdf() computes the cumulative distribution function of a normal random variable
    return value


for K in Strk:
    call = BScall(100, K, 1, 0.05, 0.2) #Calculates call price for each Strike K given
    call_list.append(call) 


#plt.plot(Strk,call_list,label='BS')
#plt.plot(Strk,Y,label='BS-FFt')
#plt.xlabel('Strike Price')
#plt.ylabel('European call option')
#plt.title('Black-Scholes call price vs Strike price')
#plt.legend()
#plt.grid(True) 


# -

# # Analytical Black-Scholes put option code

# +

# List of strikes generated from the FFT code
Strk = W['Strike Price']


put_list = [] # Call prices calculated using Black-Scholes analytical formula are appended to this list
#Btime = 
def BSput(S, K, T, r, sigma):

    float(S)
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    value = (-S * si.norm.cdf(-d1, 0.0, 1.0) + K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0))
    #si.norm.cdf() computes the cumulative distribution function of a normal random variable
    return value
  
for K in Strk:
    put = BSput(100, K, 1, 0.05, 0.2) #Calculates call price for each Strike K given
    put_list.append(put) 




#for K in Strk:
 #   start = timeit.timeit()
 #   put = BSput(100, K, 1, 0.05, 0.2) #Calculates call price for each Strike K given
  #  put_list.append(put) 
  #  stop = timeit.time()
   # duration = stop-start
   # print(duration)


#plt.plot(Strk,put_list,label='BS')
#plt.plot(Strk,Y,label='BS-FFt')
##plt.xlabel('Strike Price')
#plt.ylabel('European put option')
#plt.title('Black-Scholes put price vs Strike price')
#plt.legend()
#plt.grid(True) 



# +
#put_list
# -

 BSput(100, 150, 1, 0.05, 0.2)

plt.plot(StrkP,YP,label='BS-FFT')
plt.plot(Strk,put_list,label='BS')
#plt.plot(Strk,Y,label='BS-FFt')
plt.xlabel('Strike Price')
plt.ylabel('European put option')
plt.title('Black-Scholes put price vs Strike price')
plt.legend()
plt.grid(True) 

# ## Plotting the call price of the analytical solution against the FFT solution

plt.plot(Strk,call_list,label='BS')
plt.plot(Strk,Y,label='BS-FFT')
plt.xlabel('Strike Price')
plt.ylabel('European call value')
plt.title('Black-Scholes call price vs Strike price')
plt.legend()
plt.grid(True)



# ## Testing for the run-time

# +

elapsed_time = timeit.timeit(Btime, number=1)
print(elapsed_time)

elapsed_time = timeit.timeit(Btime, number=1)
print(elapsed_time)

elapsed_time = timeit.timeit(Btime, number=1)
print(elapsed_time)

elapsed_time = timeit.timeit(Btime, number=1)
print(elapsed_time)

elapsed_time = timeit.timeit(Btime, number=1)
print(elapsed_time)

elapsed_timefft = timeit.timeit(qfft, number=1)
print(elapsed_time)
# -

# ## True absolute error for calls

# +
Error = np.zeros(len(Strk))

for i in range(len(Strk)):
    Error[i] = np.abs(call_list[i]-Y[i])
    
# -

# ## True absolute error for puts

# +
Errorput = np.zeros(len(Strk))

for i in range(len(Strk)):
    Errorput[i] = np.abs(put_list[i]-YP[i])

# +
#Error
# -

Errp=Errorput[125:256]

Strkp=Strk[125:256]

# +
plt.plot(Strk, Error)
#plt.plot(Strk, Error)
#plt.plot(Strk,YPa,label='BS-FFt')
plt.xlabel('Strike Price')
plt.ylabel('Error')
#plt.title('Black-Scholes(FFT) put price vs Strike price')
#plt.legend()
plt.grid(True) 
#plt.ylim(0, 0.0010)
#plt.xlim(20,100)
#return df.head(10)


# -

plt.plot(Strk, Errorput)
#plt.plot(Strk,YPa,label='BS-FFt')
plt.xlabel('Strike Price')
plt.ylabel('Error')
#plt.title('Black-Scholes(FFT) put price vs Strike price')
#plt.legend()
plt.grid(True) 
#plt.ylim(0, 5)
#plt.xlim(0,10)
#return df.head(10)

# #Different values of $\alpha$.

# +
Strk = W['Strike Price']
Y = W['Option Price']
AlpCall = []


alpha1 = FFT_Pricing_CallPut(100,0.05,0.2,1,0.5,8)['Option Price']
alpha2 = FFT_Pricing_CallPut(100,0.05,0.2,1,1,8)['Option Price']
alpha3 = FFT_Pricing_CallPut(100,0.05,0.2,1,2,8)['Option Price']
alpha4 = FFT_Pricing_CallPut(100,0.05,0.2,1,5,8)['Option Price']
alpha5 = FFT_Pricing_CallPut(100,0.05,0.2,1,6,8)['Option Price']

# +
plt.plot(Strk,alpha1,label='alpha=0.5')
plt.plot(Strk,alpha2,label='alpha=1')
plt.plot(Strk,alpha3,label='alpha=2')
plt.plot(Strk,alpha4,label='alpha=5')
plt.plot(Strk,alpha5,label='alpha=6')


plt.xlabel('Strike Price')
plt.ylabel('European call value')
plt.title('Black-Scholes call price vs Strike price')
plt.legend()
plt.grid(True)



# +
#ND = np.linspace(0,Strk)

# +
#math.log(2,2)
# -



# +
Strk = W['Strike Price']
Y = W['Option Price']
AlpCall = []


alpha1 = FFT_Pricing_CallPut(100,0.05,0.2,1,-1.5,8)['Option Price']
alpha2 = FFT_Pricing_CallPut(100,0.05,0.2,1,-2,8)['Option Price']
alpha3 = FFT_Pricing_CallPut(100,0.05,0.2,1,-3,8)['Option Price']
alpha4 = FFT_Pricing_CallPut(100,0.05,0.2,1,-5,8)['Option Price']
alpha5 = FFT_Pricing_CallPut(100,0.05,0.2,1,-6,8)['Option Price']

# +
plt.plot(Strk,alpha1,label='alpha=0.5')
plt.plot(Strk,alpha2,label='alpha=1')
plt.plot(Strk,alpha3,label='alpha=2')
plt.plot(Strk,alpha4,label='alpha=5')
plt.plot(Strk,alpha5,label='alpha=6')


plt.xlabel('Strike Price')
plt.ylabel('European put value')
plt.title('Black-Scholes put price vs Strike price')
plt.legend()
plt.grid(True)
# -



StarcoderdataJupyterScriptsDedupFiltered