0
$\begingroup$

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

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Please try to use MathJax in your posts. It makes them much easier to read. $\endgroup$
    – amdopt
    Jul 12 at 10:59
  • $\begingroup$ 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. $\endgroup$
    – Alper
    Jul 13 at 17:54

2 Answers 2

Reset to default
2
$\begingroup$
  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
Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.