Monte Carlo Pricing of a Three-Asset Spread Option
Article Quant Q&A · Author: akasolace
Summary
The document examines a Monte Carlo simulation for pricing a European spread option whose payoff depends on three correlated assets. The initial attempt treats simulated normal returns as simple price multipliers, which does not model lognormal asset prices. The accepted answer instead simulates correlated standard normal shocks and applies each asset’s volatility and the lognormal drift adjustment to its forward price before calculating the payoff and discounting it.
Key ideas
- Simulate correlated normal shocks using the assets’ correlation matrix.
- Convert each shock into a lognormal forward price using volatility and time to maturity.
- Include the negative half-variance term so the simulated forward price has the intended expectation.
- Discount the average option payoff to obtain the present value.
- Monte Carlo estimates vary with sampling error, so reported values should be interpreted alongside their standard errors.
Tags
Full text
# Multivariate MC: what am I doing wrong?
# Multivariate MC: what am I doing wrong?
I am trying to generate multivariate MC results presented in this paper A Simple Generalisation of Kirk’s Approximation for Multi-Asset Spread Options by the Lie-Trotter Operator Splitting Method, by Chi-Fai Lo https://file.scirp.org/pdf/JMF_2014050615380663.pdf.
I made different attemps using different approaches but I always end up producing the same incorrect result. For example my code for generating the first element in table 5.1 produces `13.7424 +/- 0.0045` instead of `13.5763 ± 0.0089`.
Here below is the Python code I made
If someone is kind enough to tell me what I am doing wrong ....
(for info I can confirm the EK values in the paper are correct because I could find the values provided by Kirk)
```
import numpy as np
import scipy.stats
nb_simuls = 10_000_000
# parameters from https://file.scirp.org/pdf/JMF_2014050615380663.pdf Table 1 / cell 1
S1, S2, S3 = 50, 60, 150
v1, v2, v3 = 0.3, 0.3, 0.3
rho_12, rho_23, rho_13 = 0.40, 0.20, 0.80
K = 30.0
T = 0.25
r = 0.05
# from spot to forward values
F1 = S1 * np.exp(r*T)
F2 = S2 * np.exp(r*T)
F3 = S3 * np.exp(r*T)
# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3]) * np.sqrt(T)
# determine covariance matrix from correls and volatilities
correl = np.asarray([[1, rho_12, rho_13],
[rho_12, 1, rho_23],
[rho_13, rho_23, 1]])
cov = np.diag(vols).dot(correl).dot(np.diag(vols))
# simulate prices
returns = 1 + np.random.default_rng().multivariate_normal((0, 0, 0), cov, nb_simuls)
# determine exercise and value
values = []
for i in range(nb_simuls):
v = max(0, F3 * returns[i, 2] - F1 * returns[i, 0] - F2 * returns[i, 1] - K)
values.append(v)
v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"The option value is: {v} +/- {se}")
# The option value is: 13.742449851577431 +/- 0.004475778801668813
```
## Answer by akasolace (score 2, accepted)
https://quant.stackexchange.com/a/55652
Based on Quantuple comments (thank you), I fixed many mistakes and I came up with the following code:
```
import numpy as np
import scipy.stats
nb_simuls = 5_000_000
# parameters from https://file.scirp.org/pdf/JMF_2014050615380663.pdf Table 1 / cell 1
S1, S2, S3 = 50, 60, 150
v1, v2, v3 = 0.3, 0.3, 0.3
rho_12, rho_23, rho_13 = 0.40, 0.20, 0.80
K = 30.0
T = 0.25
r = 0.05
# from spot to forward values
F1 = S1 * np.exp(r*T)
F2 = S2 * np.exp(r*T)
F3 = S3 * np.exp(r*T)
# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3])
# determine covariance matrix from correls and volatilities
correl = np.asarray([[1, rho_12, rho_13],
[rho_12, 1, rho_23],
[rho_13, rho_23, 1]])
# simulate prices
returns = np.random.default_rng().multivariate_normal((0, 0, 0), correl, nb_simuls)
# determine exercise and value
values = []
for i in range(nb_simuls):
F3_exp = F3 * np.exp(-0.5 * v3 ** 2 * T + v3 * np.sqrt(T) * returns[i, 2])
F2_exp = F2 * np.exp(-0.5 * v2 ** 2 * T + v2 * np.sqrt(T) * returns[i, 1])
F1_exp = F1 * np.exp(-0.5 * v1 ** 2 * T + v1 * np.sqrt(T) * returns[i, 0])
v = max(0, F3_exp - F1_exp - F2_exp - K)
values.append(v)
v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"Method 1 (using correlation matrix): The option value is: {v} +/- {se}")
# The option value is: 13.742449851577431 +/- 0.004475778801668813
# derive volatilities from yearly volatility
vols = np. array([v1, v2, v3]) * np.sqrt(T)
# determine covariance matrix from correls and volatilities
correl = np.asarray([[1, rho_12, rho_13],
[rho_12, 1, rho_23],
[rho_13, rho_23, 1]])
cov = np.diag(vols).dot(correl).dot(np.diag(vols))
# simulate prices
returns = np.random.default_rng().multivariate_normal((0, 0, 0), cov, nb_simuls)
# determine exercise and value
values = []
for i in range(nb_simuls):
F3_exp = F3 * np.exp(-0.5 * np.diag(cov)[2] + returns[i, 2])
F2_exp = F2 * np.exp(-0.5 * np.diag(cov)[1] + returns[i, 1])
F1_exp = F1 * np.exp(-0.5 * np.diag(cov)[0] + returns[i, 0])
v = max(0, F3_exp - F1_exp - F2_exp - K)
values.append(v)
v, se = np.mean(values) * np.exp(-r*T), scipy.stats.sem(values) * np.exp(-r*T)
print(f"Method 2 (using covariance matrix): The option value is: {v} +/- {se}")
# Method 1 (using correlation matrix): The option value is: 13.5738 +/- 0.0066
# Method 2 (using covariance matrix): The option value is: 13.5803 +/- 0.0066
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.