Does anyone have any suggestions on using Monte Carlo simulations to calculate Greeks of basket option?

Question

I'd ideally like to use algorithmic differentiation or finite difference methods to approximate the Greeks of a basket option. It would be a European style basket on $N$ stocks with the payoff being $\max(B(T)-K,0)$ where $B(T)$ is the value of the basket at time $T$ and $K$ is the strike of the basket.

For the delta of this basket option, I'm assuming that at least for the finite difference approximate, it is just a matter of bumping the stock prices in the monte carlo method for the overall price of the option. I'm not sure how to factor in the correlation or covariance matrix particularly in the case of rho and vega.

Answer 1

I wrote a bit about pricing basket options here: Do basket options have a closed form valuation formula?

A vanilla option only has a single underlying, but note that a basket has $N$ underlying securities, so will have $N$ deltas, vegas etc. Even in BS, there is no exact analytical expression for the price (and hence greeks) of an arithmetic basket option, so you have two choices:

1. Semi-Analytic Expression

In the link above, we saw that the price of an arithmetic basket in BS using a 0th order moment-matching approximation is given by: \begin{align} C(0) &= \delta \bigl(F\Phi(d_{+}) - K \Phi(d_{-})\bigr)\\ d_{+} &= {\frac {\ln{\frac F K} + {\frac 1 2} \tilde{\sigma}^2 \tau} {\tilde{\sigma}\sqrt{\tau}}}\\ d_{-} &= d_{+} - \tilde{\sigma}\sqrt{\tau} \end{align}

but the values that we need to insert for $F$ and $\tilde{\sigma}$ are: \begin{align} \tilde{\sigma}^2 &= {\frac 1 {n^2}} \sum_{i,j=0}^n \rho_{ij} \sigma_i \sigma_j\\ F &= \Bigl(\prod_{i=1}^n F_i\Bigr)^{\frac 1 n} \end{align}

If you implement this expression (code is in the link), it should be easy to bump the individual underlying prices/vols etc., observe the change in the price, and calculate greeks by finite differences.

2. Numerical Methods (eg. Monte Carlo)

Below is some scruffy python code to produce the price of a basket option.

You can test the sensitivity of the option to any of the five underlying prices and vols by bumping that parameter and re-pricing. HOWEVER, note that it is important to use the SAME seed for your random number generator, so that the pricing is done along the same MC paths, or else your greeks will take a long time to converge

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm
import QuantLib as ql

spot_bump = 1e-5

initial_spots_0 = np.array([100., 100., 100., 100., 100.])
initial_spots_1 = np.array([100., 100. + spot_bump, 100., 100., 100.])

corr_mat = np.matrix([[1, 0.1, -0.1, 0, 0], [0.1, 1, 0, 0, 0.2], [-0.1, 0, 1, 0, 0], [0, 0, 0, 1, 0.15], [0, 0.2, 0, 0.15, 1]])
vols = np.array([0.1, 0.12, 0.13, 0.09, 0.11])

today = ql.Date().todaysDate()
exp_date = today + ql.Period(1, ql.Years)
strike = 100
number_of_underlyings = 5

exercise = ql.EuropeanExercise(exp_date)
vanillaPayoff = ql.PlainVanillaPayoff(ql.Option.Call, strike)

payoffAverage = ql.AverageBasketPayoff(vanillaPayoff, number_of_underlyings)
basketOptionAverage = ql.BasketOption(payoffAverage, exercise)

day_count = ql.Actual365Fixed()
calendar = ql.NullCalendar()

riskFreeTS = ql.YieldTermStructureHandle(ql.FlatForward(today, 0.0, day_count))
dividendTS = ql.YieldTermStructureHandle(ql.FlatForward(today, 0.0, day_count))

prices = []
for initial_spots in [initial_spots_0, initial_spots_1]:

    processes = [ql.BlackScholesMertonProcess(ql.QuoteHandle(ql.SimpleQuote(x)), dividendTS, riskFreeTS,
                    ql.BlackVolTermStructureHandle(ql.BlackConstantVol(today, calendar, y, day_count)))
                 for x, y in zip(initial_spots, vols)]

    process = ql.StochasticProcessArray(processes, corr_mat.tolist())

    rng = "pseudorandom"

    basketOptionAverage.setPricingEngine(
        ql.MCEuropeanBasketEngine(process, rng, timeStepsPerYear=1, requiredSamples=500000, seed=42) # requiredTolerance=0.01, 
    )

    prices.append(basketOptionAverage.NPV())

print((prices[1] - prices[0]) / spot_bump)