Suppose I'm simulating a European call option on a basket consisting of N stocks with slightly varying volatilities but all other parameters remain the same. From the perspective of an estimate, it seems simple to aggregate the values at expiry from a Monte Carlo simulation and subtract the strike of the basket. However, I'm wondering if there's a closed form equation to value basket options, ideally based off of the Black Scholes PDE, that could be used to compute the exact value as a comparison.

Edit: I edited to clarify what I meant by an option basket

  • 2
    $\begingroup$ What do you mean with a basket of N European call options? I think you mean you want to price an European option on a basket of N stocks/indices, but just making certain. $\endgroup$
    – user34971
    Commented Aug 10, 2020 at 15:25
  • 1
    $\begingroup$ That was my assumption too! $\endgroup$
    – StackG
    Commented Aug 10, 2020 at 15:45
  • 3
    $\begingroup$ @StackG has an excellent answer below (he/she provides good ways to deal with the lack of closed-form solutions, I upvoted). Since you do not necessarily require the Black-Scholes framework, you could assume your forward distributions are gaussian rather than lognormal (multivariate Bachelier model). You would then obtain a closed-form pricing formula, and it would work pretty well for short tenors. I've never tried it because moment matching (as StackG cites) is really the best approach. $\endgroup$
    – Brian B
    Commented Aug 10, 2020 at 19:41

3 Answers 3


I'm not completely certain from your question, but I'm going to assume you have a basket of $n$ stocks with prices $S_0(t)$ to $S_n(t)$, and you want to price an option with payoff at $C(\tau)$ at time $\tau$ equal to

\begin{align} C(\tau) = \max\Bigl({\frac 1 n}\sum^n_{i=1} S_i - K, 0\Bigr) \end{align} where $K$ is the strike of the option

I'm also going to make BS assumptions that each of these evolves according to geometric brownian motion in the RN measure, so

\begin{align} dS_i = S_i \bigl(r dt + \sigma_i dW_t) \end{align} where the brownian motions are potentially correlated with $n \times n$ correlation matrix $\tilde{\Sigma}$

This means that at time $\tau$, each of the stocks has a price \begin{align} S_i(\tau) &= S_i(0)\exp \Bigl( (r -{\frac 1 2} \sigma_i^2)\tau + \sigma_i \sqrt{\tau} x_i \Bigr)\\ &= F_i(0) \exp \Bigl( {\frac 1 2} \sigma_i^2\tau + \sigma_i \sqrt{\tau} x_i \Bigr) \end{align} where I have absorbed the $r$ term into the forward to simplify algebra, and the $x_i$ are variables drawn from an $n$-dimensional multivariate normal with mean $0$, variance $1$ and correlation matrix $\tilde{\Sigma}$ from above

Now the problem we face is that assuming these dynamics, we know how to price an option on a single stock, whose price is lognormally distributed, using the BS formula. But unfortunately, the $\sum^n_{i=0} S_i$ term in the payoff is not lognormally distributed because it's a sum of lognormals, not a product.

We have two choices:

  1. Price numerically using Monte-Carlo

Here is python to do that (here, for 5 stocks with an random correlation matrix I just made up)

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.stats import multivariate_normal

means = np.zeros(5)
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])

cov_mat = np.diag(vols).dot(corr_mat).dot(np.diag(vols))

initial_spots = np.array([100., 100., 100., 100., 100.])
tte = 1.0
strike = 100
seed = 43
num_paths = 50000

results = []
rng = multivariate_normal(means, cov_mat).rvs(size=num_paths, random_state=seed)

for i in range(num_paths):
    rns = rng[i]
    final_spots = initial_spots * np.exp(-0.5*vols*vols*tte) * np.exp(tte * rns)

df = pd.DataFrame(results)
df['payoff'] = ((df.sum(axis=1) / 5) - strike).clip(0)


gives price $\sim 2.09$

  1. Price APPROXIMATELY, using analytical techniques

We can use a trick here. The price of the sum of the options is not lognormally distributed, but the product of the prices is, so we CAN analytically price the contract with payoff \begin{align} C(\tau) = \max\Bigl(\bigl(\prod^n_{i=1} S_i\bigr)^{\frac 1 n} - K, 0\Bigr) \end{align}

The algebra is a bit involved (see bottom of answer), but it turns out that this simplifies to a vanilla option pricing problem, so we can price the option using the regular BS equations: \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} \sigma^2 &= {\frac 1 n}\sum_{i=1}^n \sigma_i^2\\ \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} \cdot \exp\Bigl(-{\frac 1 2} \bigl(\sigma^2 - \tilde{\sigma}^2 \bigr)\tau\Bigr) \end{align}

I've implemented that in scruffy python here as well:

mod_vol_1 = (vols ** 2).mean()
mod_vol_2 = vols.dot(corr).dot(vols) / len(vols)**2

mod_fwd = np.product(initial_spots)**(1/len(vols)) * np.exp(-0.5*tte*(mod_vol_1 - mod_vol_2))

d_plus = (np.log(mod_fwd / strike) + 0.5 * mod_vol_2 * tte) / np.sqrt(mod_vol_2 * tte)
d_minus = d_plus - np.sqrt(mod_vol_2 * tte)

mod_fwd * norm.cdf(d_plus) - strike * norm.cdf(d_minus)

price is $1.87$

How does this help us? Actually in two ways...

  • Moment Matching

This was developed for the geometric averaging basket, but it turns out that we can use a technique called Moment Matching to improve the approximation

To first order, this gives us the same equations as above except that $F$ is instead equal to simply \begin{align} F &= \Bigl(\prod_{i=1}^n F_i\Bigr)^{\frac 1 n} \end{align}

If we make this adjustment to our scruffy python, we match the Monte-Carlo price above almost exactly...

mod_vol_1 = (vols ** 2).mean()
mod_vol_2 = vols.dot(corr).dot(vols) / len(vols)**2

mod_fwd = np.product(initial_spots)**(1/len(vols))

d_plus = (np.log(mod_fwd / strike) + 0.5 * mod_vol_2 * tte) / np.sqrt(mod_vol_2 * tte)
d_minus = d_plus - np.sqrt(mod_vol_2 * tte)

mod_fwd * norm.cdf(d_plus) - strike * norm.cdf(d_minus)

price is $2.10$

  • Control Variates

We can also use the geometric option to improve our MC calculation, using the technique of Control Variates, which relies on the fact that because the price of the two types of option are highly correlated, MC paths that over-price one will tend to over-price the other, and vice-versa, which allows us to improve the convergence of the MC greatly.

And they are indeed highly correlated... here is a scatter plot of the two prices along the same paths (note the geometric basket is always cheaper than the arithmetic basket along a given path):

Geometric vs Arithmetic Baskets


In addition to StackG's answer, here is a good introductory overview of several (approximate and semi-analytical) methods to price baskets in a Black-Scholes framework:

Krekel et al - An analysis of pricing methods for basket options


There are many approximations for basket options, they all make different assumptions. The most accurate one cited for European baskets that have only positive weights is Ju 2002, which uses a Taylor expansion. If you want to take a look at a few implementations you can find one in MATLAB (type: open basketsensbyju.m) and one in Python (PyFENG package, multiasset_Ju2002.py), pip install pyfeng. https://github.com/PyFE/PyFENG - this package also has quite a few other models for Euro baskets implemented in multiasset.py, including a couple using ABM that can handle negative weights.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.