I am reading the article, where different approximations for the pricing of basket options are presented. I have tried to reproduce the result obtained by the Gentle's method in Python.

We define the price of a basket of stocks as the weighted average of the prices of $$n$$ stocks at maturity $$T$$ $$B(T) =\sum^n_{i=1} w_iS_i(T).$$

Our task is to determine the price of a call ($$\theta = 1$$) or a put ($$\theta = -1$$) with strike $$K$$ and maturity $$T$$ on the basket, i.,e. to value the payoff $$P_{Basket} (B(T), K, \theta) = [\theta (B(T) - K)]^+.$$

From the article I have used the Gentle’s approximation by geometric average and the input data.

The method.

The fact that a geometric average of log-normal random variables is again log-normally distributed allows for a Black-Scholes type valuation formula for pricing the approximating payoff. More precisely after rewriting the payoff of the basket option as $$P_{basket}=\left(\theta\left( \sum_{i=1}^n w_i S_i(T) - K \right)\right)^+ = \left(\theta\left( \left( \sum_{i=1}^n w_i F_i^T \right) \sum_{i=1}^n a_i S^*_i(T) - K \right)\right)^+,$$ where the $$T$$-forward price of stock $$i$$ $$F_i^T = S_i(0) \exp\left( \int_0^T (r(s) - d_i(s))ds \right),$$ $$r(\cdot)$$ and $$d_i(\cdot)$$ are deterministic interest rates and dividend yields, $$a_i=\frac{w_i F_i^T}{\sum_{i=1}^n w_i F_i^T}, \quad S_i^*=\frac{S_i(T)}{F_i^T}.$$ We approximate $$\sum_{i=1}^n a_i S^*_i(T)$$ by the geometric average: $$\overset{\sim}{B}(T) = \left( \sum_{i=1}^n w_i F^T_i \right) \prod_{i=1}^n (S_i^*(T))^{a_i}.$$

To correct for the mean, $$K^* = K - (E(B(T)) - E(\overset{\sim}B(T)))$$ is introduced. As approximation for $$(B(T) - K)^+$$, $$(\overset{\sim}B(T) - K^*)^+$$ is used, which -- as $$\overset{\sim}B(T)$$ is log-normally distributed -- can be valued by the Black-Scholes formula resulting in $$V_{Basket} (T) = e^{-rT}\theta\left( e^{\overset{\sim}m + \frac{1}{2} \overset{\sim}v^2} N(\theta d_1) - K^*N(\theta d_2)\right) ,$$

$$V_{basket}^{call}(T) = e^{-r \cdot T} \cdot (e^{\tilde{m} + 0.5 \cdot \tilde{v}^2} \cdot N(d_1) - K^* \cdot N(d_2))$$ and $$V_{basket}^{put}(T) = e^{-r \cdot T} \cdot (K^* \cdot N(-d_2) - e^{\tilde{m} + 0.5 \cdot \tilde{v}^2} \cdot N(-d_1)),$$ $$N(\cdot)$$ the distribution function of a standard normal random variable and $$d_1 = \frac{\overset{\sim}m - \ln K^* + \overset{\sim}v^2}{\overset{\sim}v}, \quad d_2 = d_1 - \overset{\sim}v,$$ $$\overset{\sim}m = E(\ln \overset{\sim}B(T)) = \ln \left(\sum^n_{i=1} w_i F^T_i \right) - \frac{1}{2} \sum^n_{i=1} a_i \sigma^2_i T$$ and $$\overset{\sim}v^2 = Var(\ln \overset{\sim} B(T)) = \sum^{n}_{i=1}\sum^{n}_{j=1} a_ia_j \sigma_i \sigma_j \rho_{ij}T.$$

The input data from the article. The Black-Scholes framework holds. Compute the the prices of call option on a basket, with four stocks and parameters given by $$T = 5.0$$, $$r = 0.0$$, $$\rho_{ij} = 0.5, \forall i \neq j$$, $$K = 100$$, $$F^T_i = 100$$, $$\sigma_i = 40\%$$ and $$w_i = 0.25$$, $$i, j = 1,2,3,4$$.

From the article one can see the result is $$23.78$$:

The Python code is below.

import numpy as np
from scipy.stats import norm

exp = np.exp
log = np.log
sqrt = np.sqrt
N = norm.cdf

T = 5
r = 0
rho = np.ones((4, 4)) * 0.5 + np.eye(4) * 0.5
K = 100
F = np.array([100., 100., 100., 100.])
sigma = np.array([0.4, 0.4, 0.4, 0.4])
w = np.array([0.25, 0.25, 0.25, 0.25])

a = w * F / np.sum(w * F)
K_star = K
m_tilde = log(np.sum(w * F)) - 0.5 * np.sum(a * sigma ** 2 * T)
v_tilde = sqrt((a * sigma).T @ rho @ (a * sigma) * T)

d1 = (m_tilde - log(K_star) + v_tilde ** 2) / v_tilde
d2 = d1 - v_tilde

basket_call = exp(-r * T) * (exp(m_tilde + 0.5 * v_tilde ** 2) * N(d1) - K_star * N(d2))


Basket call price: 19.294903017034283. Unfortunately, I cannot reproduce the article result.

Question. Could someone to verify the Python code? I think the weak place in the code is the command: K_star = K. How does one specify the K_star correctly?

Essentially, by approximating the algebraic average in $$B(T)$$ with the geometric average as $$\tilde B(T)$$, the payoff function should remain invariant w.r.t to this transformation

$$(B(T) - K)^+ = (\tilde B(T) - K^*)^+$$

As suggested in the method, the adjustment has the form

$$K^* = K - \mathrm{E}[B(T)] + E[\tilde B(T)]$$

where

$$\mathrm{E}[B(T)] = \Sigma_i w_i F_i(T)$$

and

$$\mathrm{E}[\tilde B(T)] = e^{ \tilde{m} + 0.5 \tilde\nu^2 }$$

Eventually this is what works for me

K_star = K - (w @ F) + exp(m_tilde + 0.5 * v_tilde ** 2)


to reproduce prices in all 4 tables.