QuantLib Python Option Pricing

In a previous article we presented how to build and install QuantLib so that it could be called in a Python script.

Here we introduce how to use QuantLib in a script that we already used in the past – though, a bit tweaked to get subplots instead of single charts – to produce some 3D charts of option sensitivities.

Below you will find the code of a class called QLAnalytics that returns a specific option figure (premium, delta, …).

import QuantLib as ql
import datetime
from Utils import to_date_vector, to_real_vector

def QLOption(args):
    try:
        calcdate = args[0]
        type = args[1]
        expiry = args[2]
        spot = args[3]
        strike = args[4]
        r = args[5]
        div = args[6]
        expirations = args[7]
        strikes = args[8]
        volMatrix = args[9]
        output = args[10]

        todaysDate = ql.Date(calcdate.day,calcdate.month,calcdate.year)
        ql.Settings.instance().evaluationDate = todaysDate
        calendar = ql.TARGET()
        dayCounter = ql.ActualActual()
        settlementDate = todaysDate
        riskFreeRate = ql.FlatForward(settlementDate, r, dayCounter)

        expirations = to_date_vector(expirations)
        strikes = to_real_vector(strikes)

        # surface
        volatilitySurface = ql.BlackVarianceSurface(todaysDate,calendar, expirations, strikes, volMatrix, dayCounter)
        volatilitySurface.enableExtrapolation()

        # option parameters
        expint = ql.Date(int(expiry))
        exercise = ql.EuropeanExercise(expint)
        payoff = ql.PlainVanillaPayoff(type, strike)

        # market data
        underlying = ql.SimpleQuote(spot)
        dividendYield = ql.FlatForward(settlementDate, div, dayCounter)

        process = ql.GeneralizedBlackScholesProcess(ql.QuoteHandle(underlying),
                                            ql.YieldTermStructureHandle(dividendYield),
                                            ql.YieldTermStructureHandle(riskFreeRate),
                                            ql.BlackVolTermStructureHandle(volatilitySurface))

        option = ql.EuropeanOption(payoff,exercise)

        engine = ql.AnalyticEuropeanEngine(process)
        
        # method: analytic
        option.setPricingEngine(engine)

        if (output == "PREMIUM"):
            value = option.NPV()
        elif (output == "DELTA"):
            value = option.delta()
        elif (output == "GAMMA"):
            value = option.gamma()
        elif (output == "DIVRHO"):
            value = option.dividendRho()
        elif (output == "RHO"):
            value = option.rho()
        elif (output == "VEGA"):
            value = option.vega()
        elif (output == "THETA"):
            value = option.theta()
        elif (output == "DPDSTRIKE"):
            value = option.strikeSensitivity()
        elif (output == "THETADAY"):
            value = option.thetaPerDay()

        return value
    except Exception:
        return float('NaN')

import QuantLib as ql

import datetime

from Utils import to_date_vector, to_real_vector

def QLOption(args):

try:

calcdate = args[0]

type = args[1]

expiry = args[2]

spot = args[3]

strike = args[4]

r = args[5]

div = args[6]

expirations = args[7]

strikes = args[8]

volMatrix = args[9]

output = args[10]

todaysDate = ql.Date(calcdate.day,calcdate.month,calcdate.year)

ql.Settings.instance().evaluationDate = todaysDate

calendar = ql.TARGET()

dayCounter = ql.ActualActual()

settlementDate = todaysDate

riskFreeRate = ql.FlatForward(settlementDate, r, dayCounter)

expirations = to_date_vector(expirations)

strikes = to_real_vector(strikes)

# surface

volatilitySurface = ql.BlackVarianceSurface(todaysDate,calendar, expirations, strikes, volMatrix, dayCounter)

volatilitySurface.enableExtrapolation()

# option parameters

expint = ql.Date(int(expiry))

exercise = ql.EuropeanExercise(expint)

payoff = ql.PlainVanillaPayoff(type, strike)

# market data

underlying = ql.SimpleQuote(spot)

dividendYield = ql.FlatForward(settlementDate, div, dayCounter)

process = ql.GeneralizedBlackScholesProcess(ql.QuoteHandle(underlying),

ql.YieldTermStructureHandle(dividendYield),

ql.YieldTermStructureHandle(riskFreeRate),

ql.BlackVolTermStructureHandle(volatilitySurface))

option = ql.EuropeanOption(payoff,exercise)

engine = ql.AnalyticEuropeanEngine(process)

# method: analytic

option.setPricingEngine(engine)

if (output == "PREMIUM"):

value = option.NPV()

elif (output == "DELTA"):

value = option.delta()

elif (output == "GAMMA"):

value = option.gamma()

elif (output == "DIVRHO"):

value = option.dividendRho()

elif (output == "RHO"):

value = option.rho()

elif (output == "VEGA"):

value = option.vega()

elif (output == "THETA"):

value = option.theta()

elif (output == "DPDSTRIKE"):

value = option.strikeSensitivity()

elif (output == "THETADAY"):

value = option.thetaPerDay()

return value

except Exception:

return float('NaN')

As you can see we rely on a file called Utils.py which contains the following functions:

import QuantLib as ql
import datetime

def to_date_vector(int_arr):
    tmpvec = ql.DateVector()
    i = 0
    for intdate in int_arr:
        intdate = int(int_arr[i])
        tmpvec.append(ql.Date(intdate))
        i = i + 1
    return tmpvec

def to_real_vector(real_arr):
    tmpvec = ql.DoubleVector()
    i = 0
    for real in real_arr:
        real = float(real_arr[i])
        tmpvec.append(real)
        i = i + 1
    return tmpvec

def pythondate_to_excel_date(pythondate):
    temp = datetime.date(1899, 12, 30)
    delta = pythondate - temp
    intdatetime = int(float(delta.days) + (float(delta.seconds) / 86400))
    return intdatetime

import QuantLib as ql

import datetime

def to_date_vector(int_arr):

tmpvec = ql.DateVector()

i = 0

for intdate in int_arr:

intdate = int(int_arr[i])

tmpvec.append(ql.Date(intdate))

i = i + 1

return tmpvec

def to_real_vector(real_arr):

tmpvec = ql.DoubleVector()

i = 0

for real in real_arr:

real = float(real_arr[i])

tmpvec.append(real)

i = i + 1

return tmpvec

def pythondate_to_excel_date(pythondate):

temp = datetime.date(1899, 12, 30)

delta = pythondate - temp

intdatetime = int(float(delta.days) + (float(delta.seconds) / 86400))

return intdatetime

And finally we need a main script that call the QLOption() function contained in QLAnalytics.py:

from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib import ticker
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.mlab import griddata
import numpy as np
from math import log
import datetime
import QuantLib as ql
import QLAnalytics
from QLAnalytics import QLOption

def ShowChart(xx, yy, zz,output,nbchart, color):

    print "Plotting " + output + " surface ..."

    ax = fig.add_subplot(3, 4, nbchart, projection='3d')
    ax.set_title(output)

    surf = ax.plot_surface(xx, yy, zz,rstride=1, cstride=1,alpha=0.65,cmap=color,vmin=zz.min(), vmax=zz.max())
    ax.set_xlabel('S')
    ax.set_ylabel('T')
    ax.set_zlabel(output)
    # Plot 3D contour
    zzlevels = np.linspace(zz.min(),zz.max(),num=3,endpoint=True)
    xxlevels = np.linspace(xx.min(),xx.max(),num=3,endpoint=True)
    yylevels = np.linspace(yy.min(),yy.max(),num=3,endpoint=True)
    cset = ax.contour(xx, yy, zz, zzlevels, zdir='z',offset=zz.min(),
    cmap=color,linestyles='dashed')
    cset = ax.contour(xx, yy, zz, xxlevels, zdir='x',offset=xx.min(),
    cmap=color,linestyles='dashed')
    cset = ax.contour(xx, yy, zz, yylevels, zdir='y',offset=yy.max(),
    cmap=color,linestyles='dashed')
    for c in cset.collections:
        c.set_dashes([(0, (2.0, 2.0))]) # Dash contours
    plt.clabel(cset,fontsize=8, inline=1)
    ax.set_xlim(xx.min(),xx.max())
    ax.set_ylim(yy.min(),yy.max())
    ax.set_zlim(zz.min(),zz.max())

    for tick in ax.xaxis.get_major_ticks():
                    tick.label.set_fontsize(6) 
    for tick in ax.yaxis.get_major_ticks():
                    tick.label.set_fontsize(6) 
    for tick in ax.zaxis.get_major_ticks():
                    tick.label.set_fontsize(6) 

    plt.xticks(np.arange(xx.min(), xx.max(), int((xx.max() - xx.min()) / 3)))
    plt.yticks(np.arange(yy.min(), yy.max(), int((yy.max() - yy.min()) / 3)))

    #Colorbar
    colbar = plt.colorbar(surf, shrink=1.0, extend='both', aspect = 10)
    l,b,w,h = plt.gca().get_position().bounds
    ll,bb,ww,hh = colbar.ax.get_position().bounds
    colbar.ax.set_position([ll, b+0.1*h, ww, h*0.8])
    tick_locator = ticker.MaxNLocator(nbins=4)
    colbar.locator = tick_locator
    colbar.update_ticks()

if __name__ == '__main__':

    # Pricing parameters
    todaysdate = datetime.date(2015,4,21)
    expirydate = datetime.date(2015,7,21)

    strike = 5500
    epsilon = 0.20

    shortexpiry = 42145
    longexpiry = 42846

    r = 0.02
    q = 0.01

    ################ Charting: 1st, IV Surface ##################

    data = np.genfromtxt('ImpliedVol.txt')
    #expiry dates
    y = data[:,0]
    #strikes
    x = data[:,1]
    #implied volatilities
    z = data[:,2]
 
    #expiries chart axis
    uniquemat = np.unique(y)
    end = log(np.max(uniquemat)) / log(uniquemat[0])
    #the expiries axis is arranged in log space
    yi = np.logspace(1, end,len(uniquemat),True,uniquemat[0])
 
    ##strike chart axis
    xi = np.unique(x)

    X, Y = np.meshgrid(xi, yi)
    Z = griddata(x, y, z, xi, yi,interp='linear')

    fig = plt.figure()
    #ex: cm.hsv, jet, terrain, rainbow, winter, summer, spring, cool. see http://matplotlib.org/users/colormaps.html
    color = cm.jet

    ################   Charting: IV Surface    ################## 

    ShowChart(X,Y,Z,"IV",1,color)
 
    ################  Chart: General settings  ##################

    dx = 10
    dy = 10

    surface = ql.Matrix(len(xi),len(uniquemat))

    dateiter = 0
    strikeiter = 0
    iter = 0

    for volitem in data:
        surface[strikeiter][dateiter] = data[iter,2] / 100.0
        iter=iter+1
        dateiter = dateiter + 1
        if (dateiter == len(uniquemat)):
            dateiter = 0
            strikeiter = strikeiter + 1

    xx, yy = np.meshgrid(np.linspace(strike*(1-epsilon), (1+epsilon)*strike, dx), np.linspace(shortexpiry,longexpiry,dy))

    ################ Charting: Premium Surface ##################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"PREMIUM"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Premium",2,color)

    ################  Charting: Delta Surface  ##################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"DELTA"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Delta",3, color)

    ################  Charting: Gamma Surface  ################## 

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"GAMMA"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Gamma",4,color)

    ################  Charting: Vega Surface  ################## 

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"VEGA"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Vega",5,color)

    ################  Charting: Theta Surface  ##################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"THETA"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Theta",6,color)
 
    ################  Charting: Rho Surface  ####################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"RHO"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Rho",7,color)

    ################  Charting: DivRho Surface  #################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"DIVRHO"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Div Rho",8,color)

    ################  Charting: dPdX Surface  ###################

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"DPDSTRIKE"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"dPdX",9,color)

    ################  Charting: Theta/1day Surface  #############

    zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"THETADAY"]) for
                   x,y in zip(np.ravel(xx), np.ravel(yy))])
    zz = zz.reshape(xx.shape)

    ShowChart(xx,yy,zz,"Theta 1d",10,color)

    # Show subplots
    plt.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

from matplotlib import cm

import matplotlib.pyplot as plt

from matplotlib import ticker

from mpl_toolkits.mplot3d import Axes3D

from matplotlib.mlab import griddata

import numpy as np

from math import log

import datetime

import QuantLib as ql

import QLAnalytics

from QLAnalytics import QLOption

def ShowChart(xx, yy, zz,output,nbchart, color):

print "Plotting " + output + " surface ..."

ax = fig.add_subplot(3, 4, nbchart, projection='3d')

ax.set_title(output)

surf = ax.plot_surface(xx, yy, zz,rstride=1, cstride=1,alpha=0.65,cmap=color,vmin=zz.min(), vmax=zz.max())

ax.set_xlabel('S')

ax.set_ylabel('T')

ax.set_zlabel(output)

# Plot 3D contour

zzlevels = np.linspace(zz.min(),zz.max(),num=3,endpoint=True)

xxlevels = np.linspace(xx.min(),xx.max(),num=3,endpoint=True)

yylevels = np.linspace(yy.min(),yy.max(),num=3,endpoint=True)

cset = ax.contour(xx, yy, zz, zzlevels, zdir='z',offset=zz.min(),

cmap=color,linestyles='dashed')

cset = ax.contour(xx, yy, zz, xxlevels, zdir='x',offset=xx.min(),

cmap=color,linestyles='dashed')

cset = ax.contour(xx, yy, zz, yylevels, zdir='y',offset=yy.max(),

cmap=color,linestyles='dashed')

for c in cset.collections:

c.set_dashes([(0, (2.0, 2.0))]) # Dash contours

plt.clabel(cset,fontsize=8, inline=1)

ax.set_xlim(xx.min(),xx.max())

ax.set_ylim(yy.min(),yy.max())

ax.set_zlim(zz.min(),zz.max())

for tick in ax.xaxis.get_major_ticks():

tick.label.set_fontsize(6)

for tick in ax.yaxis.get_major_ticks():

tick.label.set_fontsize(6)

for tick in ax.zaxis.get_major_ticks():

tick.label.set_fontsize(6)

plt.xticks(np.arange(xx.min(), xx.max(), int((xx.max() - xx.min()) / 3)))

plt.yticks(np.arange(yy.min(), yy.max(), int((yy.max() - yy.min()) / 3)))

#Colorbar

colbar = plt.colorbar(surf, shrink=1.0, extend='both', aspect = 10)

l,b,w,h = plt.gca().get_position().bounds

ll,bb,ww,hh = colbar.ax.get_position().bounds

colbar.ax.set_position([ll, b+0.1*h, ww, h*0.8])

tick_locator = ticker.MaxNLocator(nbins=4)

colbar.locator = tick_locator

colbar.update_ticks()

if __name__ == '__main__':

# Pricing parameters

todaysdate = datetime.date(2015,4,21)

expirydate = datetime.date(2015,7,21)

strike = 5500

epsilon = 0.20

shortexpiry = 42145

longexpiry = 42846

r = 0.02

q = 0.01

################ Charting: 1st, IV Surface ##################

data = np.genfromtxt('ImpliedVol.txt')

#expiry dates

y = data[:,0]

#strikes

x = data[:,1]

#implied volatilities

z = data[:,2]

#expiries chart axis

uniquemat = np.unique(y)

end = log(np.max(uniquemat)) / log(uniquemat[0])

#the expiries axis is arranged in log space

yi = np.logspace(1, end,len(uniquemat),True,uniquemat[0])

##strike chart axis

xi = np.unique(x)

X, Y = np.meshgrid(xi, yi)

Z = griddata(x, y, z, xi, yi,interp='linear')

fig = plt.figure()

#ex: cm.hsv, jet, terrain, rainbow, winter, summer, spring, cool. see http://matplotlib.org/users/colormaps.html

color = cm.jet

################ Charting: IV Surface ##################

ShowChart(X,Y,Z,"IV",1,color)

################ Chart: General settings ##################

dx = 10

dy = 10

surface = ql.Matrix(len(xi),len(uniquemat))

dateiter = 0

strikeiter = 0

iter = 0

for volitem in data:

surface[strikeiter][dateiter] = data[iter,2] / 100.0

iter=iter+1

dateiter = dateiter + 1

if (dateiter == len(uniquemat)):

dateiter = 0

strikeiter = strikeiter + 1

xx, yy = np.meshgrid(np.linspace(strike*(1-epsilon), (1+epsilon)*strike, dx), np.linspace(shortexpiry,longexpiry,dy))

################ Charting: Premium Surface ##################

zz = np.array([QLOption([todaysdate,ql.Option.Call,y, x,strike,r,q,uniquemat,xi,surface,"PREMIUM"]) for