Numerical Greeks

In this notebook, I’ll build on the facilities provided by the Instrument class (that is, its ability to detect changes in its inputs and recalculate accordingly) to show how to calculate numerical Greeks when a given engine doesn’t provide them.

As usual, we import the QuantLib module and set the evaluation date:

import QuantLib as ql

today = ql.Date(8, ql.October, 2014)
ql.Settings.instance().evaluationDate = today

A somewhat exotic option

As an example, we’ll use a knock-in barrier option:

strike = 100.0
barrier = 120.0
rebate = 30.0

option = ql.BarrierOption(
    ql.Barrier.UpIn,
    barrier,
    rebate,
    ql.PlainVanillaPayoff(ql.Option.Call, strike),
    ql.EuropeanExercise(ql.Date(8, ql.January, 2015)),
)

For the purpose of this example, the market data are the underlying value, the risk-free rate and the volatility. We wrap them in quotes, so that the instrument will be notified of any changes…

u = ql.SimpleQuote(100.0)
r = ql.SimpleQuote(0.01)
σ = ql.SimpleQuote(0.20)

…and from the quotes we build the flat curves and the process that the engine requires. As explained in another notebook, we build the term structures so that they move with the evaluation date; this will be useful further on.

riskFreeCurve = ql.FlatForward(
    0, ql.TARGET(), ql.QuoteHandle(r), ql.Actual360()
)
volatility = ql.BlackConstantVol(
    0, ql.TARGET(), ql.QuoteHandle(σ), ql.Actual360()
)

process = ql.BlackScholesProcess(
    ql.QuoteHandle(u),
    ql.YieldTermStructureHandle(riskFreeCurve),
    ql.BlackVolTermStructureHandle(volatility),
)

Finally, we build the engine (the library provides one based on an analytic formula) and set it to the option.

option.setPricingEngine(ql.AnalyticBarrierEngine(process))

Now we can ask the option for its value…

print(option.NPV())

29.24999528163375

…but we’re not so lucky when it comes to Greeks:

try:
    print(option.delta())
except Exception as e:
    print(f"{type(e).__name__}: {e}")

RuntimeError: delta not provided

The engine doesn’t provide the delta, so asking for it raises an error.

Numerical calculation

What does a quant have to do? We can use numerical differentiation to approximate the Greeks: that is, we can approximate the derivative by calculating the option value for two slightly different values of the underlying and by taking the slope between the resulting points.

The relevant formulas are:

\[ \Delta = \frac{P(u_0+h)-P(u_0-h)}{2h} \; \; \; \; \; \; \Gamma = \frac{P(u_0+h)-2P(u_0)+P(u_0-h)}{h^2} \]

where \(P(u)\) is the price of the option for a given value of the underlying \(u\), \(u_0\) is the current value of the underlying and \(h\) is a tiny increment.

Thanks to the framework we set in place, getting the perturbed prices is easy enough: we can set the relevant quote to the new value and ask the option for its price again. Thus, we choose a value for \(h\) and start. First, we save the current value of the option…

u0 = u.value()
h = 0.01

P0 = option.NPV()
print(P0)

29.24999528163375

…then we increase the underlying value and get the new option value…

u.setValue(u0 + h)
P_plus = option.NPV()
print(P_plus)

29.24851583696518

…then we do the same after decreasing the underlying value.

u.setValue(u0 - h)
P_minus = option.NPV()
print(P_minus)

29.25147220877793

Finally, we set the underlying value back to its current value.

u.setValue(u0)

Applying the formulas above give us the desired Greeks:

Δ = (P_plus - P_minus) / (2 * h)
Γ = (P_plus - 2 * P0 + P_minus) / (h * h)
print(Δ)
print(Γ)

-0.14781859063752734
-0.025175243933972524

The approach works for almost any Greeks. We can use the two-sided formula above, or the one-sided formula below if we want to minimize the number of evaluations:

\[ \frac{\partial P}{\partial x} = \frac{P(x_0+h)-P(x_0)}{h} \]

For instance, here we calculate Rho and Vega:

r0 = r.value()
h = 0.0001

r.setValue(r0 + h)
P_plus = option.NPV()
r.setValue(r0)

ρ = (P_plus - P0) / h
print(ρ)

-9.984440333461464

σ0 = σ.value()
h = 0.0001

σ.setValue(σ0 + h)
P_plus = option.NPV()
σ.setValue(σ0)

Vega = (P_plus - P0) / h
print(Vega)

-12.997818541258255

The approach for the Theta is a bit different, although it still relies on the fact that the option reacts to the change in the market data. The difference comes from the fact that we don’t have the time to maturity available as a quote, as was the case for the other quantities. Instead, since we set up the term structures so that they move with the evaluation date, we can set it to tomorrow’s date to get the corresponding option value. In this case, the value of the time increment \(h\) should be equivalent to one day:

ql.Settings.instance().evaluationDate = today + 1
P1 = option.NPV()
ql.Settings.instance().evaluationDate = today

h = 1.0 / 365
Θ = (P1 - P0) / h
print(Θ)

5.548998890566246

Why not in the engine?

We just saw that it’s possible to calculate the missing ones numerically. So why don’t engines use this method when Greeks are not available otherwise?

It’s a good question; there are a few layers to it, and it can lead us to some insights into the architecture of the library.

The first answer is that numerical Greeks are costly: each one needs to reprice the instrument once or twice, depending on the accuracy we want. Therefore, we don’t want them to be calculated by default. You should incur the cost only if you need them.

That doesn’t need to stop us, though. We could configure the engine by passing argument to its constructor that specify which Greeks to calculate. It would add a bit of housekeeping, but it’s doable, so let’s assume we went ahead with the plan. (I won’t bother showing you the code keeping track of the configuration. That’s boring stuff, and moreover, I don’t have it.)

How would we write the calculation of the Greeks, then? As an example, let’s take the MCEuropeanEngine class. It takes as its first input a shared_ptr to a GeneralizedBlackScholesProcess (I’m talking of the C++ class, of course, since that’s where we would implement the change); in turn, this contains a Handle<Quote> containing the value of the underlying, two Handle<YieldTermStructure> holding the risk-free rate and dividend yield curves, and a Handle<BlackVolTermStructure> holding the volatility. Those are the inputs we want to bump before repricing the option.

It turns out that it’s fairly difficult to modify them. Even in the simplest case, that is, the Handle<Quote> containing the underlying, we have no guarantee that the contained Quote can be cast to a SimpleQuote and set a new value; it couls be another user-defined class. The same holds for the other handles; without knowledge of the actual derived class of the contained YieldTermStructure or BlackVolTermStructure instances, there’s no way we can modify them.

We don’t actually need to do it, though, right? Let’s take, e.g., the risk-free curve; we can wrap the current instance in, say, a ZeroSpreadedTermStructure instance to shift it upwards or backwards a given increment. Something like:

Handle<Quote> dr(make_shared<SimpleQuote>(1e-4)); // 1 bp

auto new_curve =
    make_shared<ZeroSpreadedTermStructure>(process_->riskFreeRate(), dr);

The same goes for the underlying quote: we can ask for the current value, perturb it, and store it in a new SimpleQuote instance, as in:

Real du = 0.001;

auto new_u =
    make_shared<SimpleQuote>(process_->x0() + du);

However, the problem now is that we can’t put the new inputs into the process. It contains Handle instances, and those can’t be relinked so that they store the new curve and quote; only RelinkableHandle instances can. Hmm. It seems like we went out of our way to make this operation difficult, doesn’t it?

Well, in fact, we did—we wrote the code explicitly so that the inputs couldn’t be changed. We had a reason for that, though: the instrument might not be the only one currently using them. If we had a bunch of options on the same underlying, for instance, they could be sharing the same process instance; and instruments in the same currency would probably share the same risk-free curve. There’s no point creating and maintaining an instance for each instrument.

This means that any one instrument isn’t free to modify its inputs: it could affect other, unrelated instruments. That’s why inputs are only ever stored and made available through handles. (As a historical side note, we were once bitten by this before we started coding this way. In fact, we were bitten so badly we had to retire a release because it was too dangerous to use: version 0.3.5 was removed from the download page and version 0.3.6 replaced it as soon as we found out. You can check the fix here if you’re curious.)

If you’re determined enough and really want to add numerical Greeks to your engine, you can still find a way. Instead of modifying the process, you’ll have to clone it. For instance, if you want to shift the risk-free rate, you can create a shifted curve as above and write:

auto new_process =
    make_shared<GeneralizedBlackScholesProcess>(
        process_->stateVariable(),
        process_->dividendYield(),
        Handle<YieldTermStructure>(new_curve),
        process_->blackVolatility());

after which the engine can perform the calculations using the cloned process.

You might also want to think what should happen if the Greeks calculations were to raise an exception. Do you throw away the calculated price, or keep it and just bail on the Greeks, and maybe try the next one? Granted, this requires some thought even for ordinary Greeks; but those are usually calculated together with the price, and it’s more likely that they succeed or fail together.

In the end, though, we preferred to leave it to users to calculate missing Greeks. On the one hand, they can do it much more easily by perturbing the inputs that they set up and probably still hold on to; and on the other hand, unlike us, they’re also able to optimize the calculation when multiple instruments depend on the same process. There’s no need for each of them to create a new process; users can modify the relevant input once, read the new prices from all the instruments, and restore the input to its former value.