Quantlib: Greeks of FX option in Python

Question

I'm using Quantlib in Python to price an FX option. I'm comparing the result to Bloomberg, to make sure the code is working correct.

I also want to calculate all the Greeks, and eventually use those in a Taylor expansion of the P&L (as in for example: P&L of delta hedged call option)

The option I'm trying to price, is priced in Bloomberg as follows:

It is a stylized example.

The code I use is as follows:

import QuantLib as ql


Spot = 1.1
Strike = 1.101
Sigma = 10/100
Ccy1Rate = 5/100
Ccy2Rate = 10/100
OptionType = ql.Option.Call

#Option dates in quantlib objects
EvaluationDate = ql.Date(3, 1,2022)
SettlementDate = ql.Date(5, 1, 2022) #Evaluation +2
ExpiryDate = ql.Date(10, 1, 2022) #Evaluation + term which is 1 week
DeliveryDate = ql.Date(12, 1, 2022) #Expiry +2
NumberOfDaysBetween = ExpiryDate - EvaluationDate
#print(NumberOfDaysBetween)

#Generate continuous interest rates
EurRate = Ccy1Rate
UsdRate = Ccy2Rate

#Create QuoteHandle objects. Easily to adapt later on.
#You can only access SimpleQuote objects. When you use setvalue, you can change it.
#These global variables will then be used in pricing the option.
#Everything will be adaptable except for the strike.
SpotGlobal = ql.SimpleQuote(Spot)
SpotHandle = ql.QuoteHandle(SpotGlobal)
VolGlobal = ql.SimpleQuote(Sigma)
VolHandle = ql.QuoteHandle(VolGlobal)
UsdRateGlobal = ql.SimpleQuote(UsdRate)
UsdRateHandle = ql.QuoteHandle(UsdRateGlobal)
EurRateGlobal = ql.SimpleQuote(EurRate)
EurRateHandle = ql.QuoteHandle(EurRateGlobal)

#Settings such as calendar, evaluationdate; daycount
Calendar = ql.UnitedStates()
ql.Settings.instance().evaluationDate = EvaluationDate
DayCountRate = ql.Actual360()
DayCountVolatility = ql.ActualActual()

#Create rate curves, vol surface and GK process
RiskFreeRateEUR = ql.YieldTermStructureHandle(ql.FlatForward(0, Calendar, EurRateHandle, DayCountRate))
RiskFreeRateUSD = ql.YieldTermStructureHandle(ql.FlatForward(0, Calendar, UsdRate, DayCountRate))
Volatility = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(0, Calendar, VolHandle, DayCountVolatility))
GKProcess = ql.GarmanKohlagenProcess(SpotHandle, RiskFreeRateEUR, RiskFreeRateUSD, Volatility)

#Generate option
Payoff = ql.PlainVanillaPayoff(OptionType, Strike)
Exercise = ql.EuropeanExercise(ExpiryDate)
Option = ql.VanillaOption(Payoff, Exercise)
Option.setPricingEngine(ql.AnalyticEuropeanEngine(GKProcess))
BsPrice = Option.NPV()



ql.Settings.instance().evaluationDate = EvaluationDate
print("Premium is:", Option.NPV()*1000000/Spot)
print("Gamma is:", Option.gamma()*1000000*Spot/100)
print("Vega is:", Option.vega()*1000000*(1/100)/Spot)
print("Theta is:", Option.theta()*1000000*(1/365)/Spot)
print("Delta is:", Option.delta()*1000000)

Which gives the next output:

Premium is: 5550.960519027888
Gamma is: 287777.2550015351
Vega is: 551.9015849344515
Theta is: -462.68771985750703
Delta is: 504102.4957777005

It matches Bloomberg very well, except for Theta. I tried to divide by 255 (workdays) instead of 365, but that's also wrong.

I'm wondering what the correct answer is, since it's necessary for finding the Taylor expansion P&L.

Answer 1

As mentioned in the comments, theta in BBG is bump and reprice. By definition, Greeks measure sensitivities keeping all else equal. In this case, it means keeping all inputs constant, but moving the date one day forward, and compute the value difference, which is theta.

In quantlib, it's model theta (closed form Garman Kohlhagen).

The language below is not Python, but Julia because I already had this code. However, the syntax is sufficiently similar to Python so it should be possible to follow the logic. I manually replicate Bloomberg's as well as Quantlib's computation to illustrate the differences.

Import all packages and define the cdf.

using Distributions
N(x) = cdf(Normal(0,1),x)

Define inputs and compute continuous rates

spot = 1.1
f = 1.101070
strike = 1.101
ccy1 = 0.05 # EUR
ccy2 = 0.1 # USD
vol = 0.1
days = 7
t = days/365

r1_cont = log(1+ccy1*days/360)/(days/365)
r2_cont = log(1+ccy2*days/360)/(days/365)

Define Garman Kohlhagen with forward (technically Black76) => same result as can be seen here.

function GKF(F,K,t,ccy2,σ)
    d1 = ( log(F/K) +  0.5*σ^2*t ) / (σ*sqrt(t))
    d2 = d1 - σ*sqrt(t)
    c  = exp(-ccy2*t)*(F*N(d1) - K*N(d2))
    p  =  exp(-ccy2*t)*(-F*N(-d1) + K*N(-d2))
  return c, p, d1, d2
end

Define todays option value and compute tomorrows, with exact same inputs but one day less to expiry. Theta is the numerical difference between the two NPVs.

t1 = GKF(f,strike, days/365, r2_cont, vol)[1]*1000000/spot
t2 = GKF(f,strike, (days-1)/365, r2_cont, vol)[1]*1000000/spot
theta = t2-t1

There are minor rounding differences to BBG because I only used the fwd points visible in the screenshot, which lacks the exact decimal precision.

In terms of quantlib, instead of

print("Theta is:", Option.theta()*1000000*(1/365)/Spot)

you could use

print("Theta is:", Option.thetaPerDay()*1000000/Spot)

to save the manual computation of (/365).

QL uses the closed form Garman Kohlhagen formula (the actual derivative of the option price with respect to time). Theta is defined here for example.

If we were to define this manually, we could write the following Julia code.

function θ(S,K,t,ccy2,ccy1, σ)
    d1 = (log(S/K) +  (ccy2-ccy1+0.5*σ^2)*t ) / (σ*sqrt(t))
    d2 = d1 - σ*sqrt(t)
    thetaGK = (-(S*exp(-ccy2*t)*n(d1)*σ)/(2*sqrt(t)) + ccy1*exp(-ccy1*t)*S*N(d1) - ccy2*exp(-ccy2*t)*K*N(d2))/365
    return thetaGK
end

This yields -462 like ql:

Since I do not use quantlib, I cannot comment if ql can compute finite difference theta. Quickly skipping through the c++ code makes me think it does not. The BS implementation can be seen here

It is the standard closed form theta, although I am not sure why it was implemented with this specific syntax. In any case, the derivative of the price vs time refers to the change per unit time (the change after one year). In other words, mathematically the result of the formula for theta is expressed in value per year. Since the model is continuous, you need to use 365 days (not workdays, or sqrt or whatever).

Delta as displayed in BBG (based on your setting) is computed in ql like this:

EURUSD is one of the exceptions where usually delta is not premium included. You can cross check the conventions on OVDV, for any fx pair, by clicking on Settings->Conventions.