# Pricing European Call Closed Form Spread Options in Python

I am currently trying to correctly price European Call Closed Form Spread Options using Python. The main problem I am currently running into is that I have nothing to compare the option price so that I know that my code is working correctly.

Does anyone know if this code is correct ? Am I missing something ? Any help / info is appreciated.

This is my code:

import numpy as np

S1 = 110  # Spot price of first asset
S2 = 100  # Spot price of second asset
K = 25    # Strike price
T = 1     # Time to maturity in years
r = 0.05  # Risk-free interest rate
sigma1 = 0.1  # Volatility of first asset
sigma2 = 0.15  # Volatility of second asset
rho = 0.8    # Correlation between the two assets
num_simulations = 10000000  # Number of Monte Carlo simulations

def monte_carlo_spread_option_price(S1, S2, K, T, r, sigma1, sigma2, rho, num_simulations):
np.random.seed(0) # for reproducible results

# Simulating correlated paths
dt = T
Z1 = np.random.standard_normal(num_simulations)
Z2 = np.random.standard_normal(num_simulations)
Z3 = rho * Z1 + np.sqrt(1 - rho**2) * Z2

# GBM formula
S1_T = S1 * np.exp((r - sigma1**2 / 2) * dt + sigma1 * np.sqrt(dt) * Z1)
S2_T = S2 * np.exp((r - sigma2**2 / 2) * dt + sigma2 * np.sqrt(dt) * Z3)

# Calculate the payoff for each path at maturity (T)
payoff_T = np.maximum(S1_T - S2_T - K, 0)

# Average the payoffs and discount back to present value
option_price = np.exp(-r * T) * np.mean(payoff_T)

return option_price

option_price = monte_carlo_spread_option_price(S1, S2, K, T, r, sigma1, sigma2, rho, num_simulations)
print(f'The price of the spread call option is: {option_price}')


I am also trying to add control variates to reduce the simulation error, something like this:

$$E[e^{-rT} \cdot C(T)] = c + E[e^{-rT}(C(T)-c(T))]$$

where

$$C(T) = \max((S_1(T) - S_2(T) - K),0)$$,

\begin{align} c(T) &= C(a,b) \\ & = E_0 e^{-rT} \left((S_1(T) - S_2(T) - K) I \left(S_1(T) \geq \frac{a (S_2(T))^b}{E[(S_2(T))^b]}\right)\right), \end{align}

$$a = F_2 + K$$,

$$b = \frac{F_2}{F_2+K}$$,

and where $$F_2$$ is the forward price of the asset $$S_2$$.

• Please try to use MathJax in your posts. It makes them much easier to read. Jul 12, 2023 at 10:59
• I tried to rewrite your formulas using MathJax but the one involving c(T) appears to be a bit confusing. The changes should be visible to you as well if and when approved. Feel free to correct them as per necessary. In addition, though it looks like the output of your code is simply a single number, it is good practice to include in your post whatever output you are getting when you run yourself exactly the code you have shared. Jul 13, 2023 at 17:54

## 2 Answers

1. Firstly, you can really use Margrabe's formula (as @Rylan said). It's exact solution, so there is good first test.
2. You can use numerical integration (scipy methods for example) for test your Monte-Carlo realization. Let $$S_1(T) = F_1(T) e^{-\frac12 \sigma_1^2 T + \sigma_1 x \sqrt{T}}, \\ S_2(T) = F_2(T) e^{-\frac12 \sigma_2^2 T + \sigma_2 y \sqrt{T}}, \\ \Sigma = \left(\begin{array}{ccc} 1 & \rho\\ \rho & 1 \end{array}\right), \\ \theta (x) = Heaviside \; step \; function$$

$$e^{-rT} E \left[ \max \left( S_1(T) - S_2(T) - K , 0 \right) \right] = \\ \theta \left( S_1(T) - S_2(T) - K \right)\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}\left( S_1(T) - S_2(T) - K \right) e^{-\frac12 (x \quad y)^{'} \Sigma^{-1} (x \quad y)} \cfrac{dx dy}{2 \pi \sqrt{\det \Sigma}}$$

from scipy.integrate import dblquad
import numpy as np
from scipy.stats import norm

S_10 = 110  # Spot price of first asset
S_20 = 100  # Spot price of second asset
K = 10 # Strike price
T = 1     # Time to maturity in years
r = 0.05  # Risk-free interest rate
sigma1 = 0.1  # Volatility of first asset
sigma2 = 0.15  # Volatility of second asset
rho = 0.8    # Correlation between the two assets
num_simulations = int(1e8)  # Number of Monte Carlo simulations

F_1 = S_10*np.exp(r*T)
F_2 = S_20*np.exp(r*T)

def monte_carlo_spread_option_price(S1, S2, K, T, r, sigma1, sigma2, rho, num_simulations):
np.random.seed(0) # for reproducible results

# Simulating correlated paths
dt = T
Z1 = np.random.standard_normal(num_simulations)
Z2 = np.random.standard_normal(num_simulations)
Z3 = rho * Z1 + np.sqrt(1 - rho**2) * Z2

# GBM formula
S1_T = S_10 * np.exp((r - sigma1**2 / 2) * dt + sigma1 * np.sqrt(dt) * Z1)
S2_T = S_20 * np.exp((r - sigma2**2 / 2) * dt + sigma2 * np.sqrt(dt) * Z3)

# Calculate the payoff for each path at maturity (T)
payoff_T = np.maximum(S1_T - S2_T - K, 0)

# Average the payoffs and discount back to present value
option_price = np.exp(-r * T) * np.mean(payoff_T)

return option_price

def numericalSpreadOptionCallPrice(x, y):
cov = np.array([[1, rho],
[rho, 1]]).reshape(2, 2)
inverse_cov = np.linalg.inv(cov)
S1 = F_1*np.exp(-1/2*sigma1**2*T + sigma1*x*np.sqrt(T))
S2 = F_2*np.exp(-1/2*sigma2**2*T + sigma2*y*np.sqrt(T))
Heaviside = np.heaviside(S1 - S2 - K, 1)

xi = np.array([x, y]).reshape(2, 1)
xi_transpose = xi.transpose()

quadratic_form = float(np.dot(xi_transpose, np.dot(inverse_cov, xi)))

Result = (1 / (2 * np.pi * ((np.sqrt(np.linalg.det(cov)))))) * Heaviside * (S1 - S2 - K) * np.exp(
-1 / 2 * quadratic_form)
return Result

def margrabeCall(F_1, F_2, sigma_1, sigma_2, rho, T, r, accuracy=10):

N = norm.cdf
def sigma_evaluation(sigma_1, sigma_2, rho):
sigma = np.sqrt(sigma_1**2 + sigma_2**2 - 2*rho*sigma_1*sigma_2)
return sigma

sigma = sigma_evaluation(sigma_1, sigma_2, rho)
d_1 = (np.log(F_1/F_2) + 0.5*sigma**2*T)/sigma/np.sqrt(T)
d_2 = d_1 - np.sqrt(T) * sigma
C = np.exp(-r*T)*(F_1*N(d_1) - F_2 * N(d_2))

return round(C, accuracy)

_round = 8
option_price = monte_carlo_spread_option_price(S_10, S_20, K, T, r, sigma1, sigma2, rho, num_simulations)
print(f'The price of the spread call option is: {round(option_price, _round)}')
NumericalPrice = round(np.exp(-r*T)*dblquad(numericalSpreadOptionCallPrice, -9, 9, lambda x: -9, lambda x:9)[0], _round)
print("Numerical integration reslut: %.8f" % round(NumericalPrice, _round))

# out:
# The price of the spread call option is: 3.86809862 (with 1e6 paths)
# The price of the spread call option is: 3.87142894 (with 1e8 paths)
# Numerical integration result: 3.87177977

print('\nMargrabe test')
K=0
NumericalPrice = round(np.exp(-r*T)*dblquad(numericalSpreadOptionCallPrice, -9, 9, lambda x: -9, lambda x:9)[0], 9)
option_price = monte_carlo_spread_option_price(S_10, S_20, K, T, r, sigma1, sigma2, rho, num_simulations)
print(f'The price of the spread call option is: {round(option_price, _round)}')
print("Numerical Integration: %.8f" % round(NumericalPrice, _round))
print("Margrabe analitical formula: %.8f" % round(margrabeCall(F_1, F_2, sigma1, sigma2, rho, T, r), _round))

# out:
# Margrabe test
# The price of the spread call option is: 10.74810401 (with 1e6 paths)
# The price of the spread call option is: 10.75422032 (with 1e8 paths)
# Numerical Integration: 10.75459894
# Margrabe analitical formula: 10.75459891


For the case of $$K=0$$, you can use Margrabe's formula as a diagnostic. This can help you figure out whether the code is working incorrectly, at least.