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:
- 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)
results.append(final_spots)
df = pd.DataFrame(results)
df['payoff'] = ((df.sum(axis=1) / 5) - strike).clip(0)
df['payoff'].mean()
gives price $\sim 2.09$
- 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...
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$
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):
