So I coded up the solution from here: Do basket options have a closed form valuation formula?
Which provides a good solution for equally-weighted underlyings under a Black model. The simplified Python code is here after modifying it quite a bit. I tried to add a weight vector to the code which seems to diverge when ATM from the referenced solution, so I'm pretty sure my implementation is wrong (the analytical one, not MC). If I set all the weights back to equal (sum=1), then the solution seems to converge when ATM and OTM pretty well, so something is definitely wrong with my weight formula in the analytical solution. Hoping someone can spot my error here:
import numpy as np
from scipy.stats import multivariate_normal, norm
def Basket_Euro_MC(F, vols, corr, T, strike, paths, weights, callput):
seed = 1
means = np.zeros(F.shape[0])
cov_mat = np.diag(vols).dot(corr).dot(np.diag(vols))
results = np.zeros((paths, F.shape[0]))
rng = multivariate_normal(means, cov_mat).rvs(size=paths, random_state=seed)
for i in range(paths):
results[i] = F * np.exp(-0.5*vols**2*T) * np.exp(T * rng[i])
if callput == 1:
return max(np.mean(np.sum(results*weights, axis=1)-strike),0)
else:
return max(np.mean(strike-np.sum(results*weights, axis=1)),0)
def Basket_Euro_Analytic(F, vols, corr, T, strike, weights, callput):
mod_vol = vols.dot(corr).dot(vols) / len(vols)**2
mod_fwd = np.product(F)**(1/len(vols))
d_plus = (np.log(mod_fwd / max(strike,0.0001)) + 0.5 * mod_vol * T) / np.sqrt(mod_vol * T)
d_minus = d_plus - np.sqrt(mod_vol * T)
if callput == 1: # call formula is fine
return np.sum(weights.dot(mod_fwd * norm.cdf(d_plus) - strike * norm.cdf(d_minus)))
else: # put formula looks pretty suspect...
return np.sum(weights.dot(max(-(norm.cdf(-d_minus) * mod_fwd - strike * norm.cdf(-d_plus)),0)))
if __name__ == '__main__':
F = np.array([80., 85., 82., 81., 84.])
corr = np.array(([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])
paths = 100000
T=1
strike = 78# note BS makes strike >0 in divide!
weights = np.array([0.25,0.1,0.2,0.2,0.25])
callput = 1
MC_result = Basket_Euro_MC(F, vols, corr, T, strike, paths, weights, callput)
analytic_result = Basket_Euro_Analytic(F, vols, corr, T, strike, weights, callput)
print('MC result:', MC_result, 'Analytic Result', analytic_result)