Early Sale of an American Call at a Fair Price

# Pricing European Call Closed Form Spread Options in Python


# 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$.

## Answer by K. Roman (score 2)

https://quant.stackexchange.com/a/76161

- Firstly, you can really use Margrabe's formula (as @Rylan said). It's exact solution, so there is good first test.

- 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
```

## Answer by Rylan (score 0)

https://quant.stackexchange.com/a/76103

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.