Monte Carlo approximation of call option on the maximum of two assets

I want to compute the price of the option with payoff $$$$\max \big\{\max\{S^1_T, S^2_T\} - K, 0\big\},$$$$ where $$S^{1,2}$$ have the same dynamics with 0 correlation. So, \begin{align} dS^1_t &= r S_t^1 dt + \sigma S^1_t dW^1_t \\ dS^2_t &= r S_t^2 dt + \sigma S^2_t dW^2_t, \end{align} where $$W^1$$ and $$W^2$$ are independent Wiener processes under the pricing measure $$Q$$. This option has analytical pricing formula (e.g. in The Complete Guide to Option Pricing Formulas, p.211). However, when I try to compute the value of this option using MC method, I get values that are consistently incorrect.

Below is my code for the MC simulation. First a function to do the numerical integration of the SDEs:

# Euler scheme for two GBMs (no correlation) with same drift and volatility
# Returns the terminal value (prices at last time step)
gbm <- function(mu, sigma, max_time, num_steps, init_value){
h = max_time / num_steps
paths <- matrix(NA, num_steps+1, 2)
paths[1, ] = init_value
normals = matrix(rnorm(num_steps*2, sd=sqrt(h)), num_steps, 2)

for (i in 1:num_steps){
paths[i+1, ] = paths[i, ] + (mu * paths[i, ] * h) + (sigma * paths[i, ] * normals[i, ])
}
return(paths[num_steps, ])
}


Then the Monte Carlo method. Note that I compute the price for a call option on the maximum AND the price for just a vanilla call:

trials <- 10000
maxes <- array(NA, trials)
max_payoffs <- array(NA, trials)
vanilla_payoffs <- array(NA, trials)

for(i in 1:trials){
# Compute terminal values of the SDEs
terminal_values <- gbm(mu=0.02, sigma=0.2, max_time=3, num_steps=1000, init_value=c(1, 1))
# Vanilla call payoff just on one of the GBM - for assuring my numerical integration correct
vanilla_payoffs[i] <- max(terminal_values[1] - 1, 0)

# Call on the maximum of the two assets - strike 1
maxes[i] = max(terminal_values)
max_payoffs[i] = max(maxes[i] - 1, 0)
}

# Mean of the payoffs + 95% confidence interval
mean(max_payoffs) * exp(-0.02 * 3)
sd(max_payoffs * exp(-0.02 * 3)) * 2 / sqrt(trials)

# Mean of the vanilla call payoffs
mean(vanilla_payoffs) * exp(-0.02 * 3)


For the call on the max of two assets, my sample mean is $$0.2839 \pm 0.0064$$ which is very far off the correct value of $$0.2235$$. However my vanilla call option is almost exactly right $$0.1656$$ compared to the true value $$0.1646$$.

Just to be clear, the parameters are $$\sigma=0.2$$, $$r=0.02$$, $$S^{1,2}_0 = 1$$, $$K=1$$, $$T=3$$, $$\rho=0$$.

I'd be very grateful if anyone could explain where I am going wrong.

EDIT: I added python code which uses no numerical integration as per @Yoda And Friends's answer. It still gives incorrect price though:

def terminal_spots(trials, r, sigma, t, spot):
normals = np.random.normal(size = (trials, 2))
return spot * np.exp(t * (r - 0.5 * sigma * sigma) + sigma * np.sqrt(t) * normals)


and

def mc_call_max_two_assets(trials, r, sigma, t, spot, strike):
terminals = terminal_spots(trials, r, sigma, t, spot)
max_terminal = terminals.max(1)
payoffs = np.maximum(max_terminal - strike, 0)
mn = payoffs.mean() * np.exp(-r*t)
conf_interval = (payoffs * np.exp(-r*t)).std() * 2 / np.sqrt(trials)
return mn, conf_interval

• Are dynamics without $dt$? Or is it a typo? May 28, 2021 at 8:26
• @YodaAndFriends Sorry, yes it was a typo, edited now May 28, 2021 at 8:36

I image you want to calculate the following payoff: $$\pi_T = \max\left[ \max(S_T^1, S_T^2) - K, 0 \right]$$ If dynamics are expressed with the following dynamic (from your code, it should be the case): $$dS_t^i = rS_t^idt + \sigma S_t^idW_t^i$$ where $$W = (W^1, W^2)$$ is two dimensional SBM, then, I guess you could simplify your code a bit!

From your definition of the problem, the contingent claim is path-independent. This means you are only interested in finding $$S_T^i$$ which is a just a random variable. Let us see if we can find its distribution! Take the following transformation: $$d(\log(S_t^i)) = (r - \frac{1}{2}\sigma^2)dt + \sigma dW_t^i$$ which is just short hand notation for the following integral representation: $$\log(S_t^i) - \log(S_0^i) = \int_0^t (r - \frac{1}{2}\sigma^2)ds + \int_0^t \sigma dW_s^i$$ Being $$r$$ and $$\sigma$$ constant, and since $$\int_0^t \sigma dW_s^I \sim N(0, \int_0^t \sigma^2) = N(0, \sigma^2 t)$$ we have that: $$\log(S_t^i) \sim N\left(\log(S_0^i) + (r - \frac{1}{2}\sigma^2)t, \ \ \sigma^2 t\right)$$ We can then "simulate" a standard normal distribution $$Z \sim N(0, 1)$$ using our laptop: $$\log(S_t^i) \sim \log(S_0^i) + (r - \frac{1}{2}\sigma^2)t + \sigma \sqrt{t} Z$$ With this in mind: $$S_t = \exp(\log(S_t))$$. So that our Monte Carlo method boils down to this:

def generateStock(initialValue, sigma, rfr, t, randomGenerator):
z = RandomGenerator.nextGaussian(0, 1)
return initialValue*Math.exp((rfr - 0.5*sigma*sigma)*t + sigma*Math.sqrt(t)*z)

def monteCarloMaxOption(initialValue, sigma, rfr, strike, maturity, numberOfSimulations):
sample = [0]*numberOfSimulations #Just create an array of length nos
for w in range(numberOfSimulations){
stock1 = generateStock(initialValue, sigma, rfr, maturity)
stock2 = generateStock(initialValue, sigma, rfr, maturity)
maximumStock = Math.max(stock1, stock2)
sample[w] = Math.max(maximumStock - strike, 0.0)
}
discountFactor = Math.exp(-rfr*maturity)
priceToday = mean(sample)*discountFactor
sampleError = errorFunction(sample)
return priceToday, error


I am sorry to have used pseudo code (Python style with Java influences... :D) but I find R really difficult!

If the payoff is not what you intended, please let me know!!!

• Thanks so much! No problem with the pseudo code. Can you look at my edits to the question? I followed your approach (in Python), which gives answer consistent with my R code. However, it still does not match the true price? Any ideas what I'm doing wrong? May 28, 2021 at 10:48
• I actually tried myself with a tested library I have built. I get that the option value is 0.285188 with an error of 9.5E-4. Are you completely sure the theoretical one is 0.22? If so, are the parameters you provided here correct? May 28, 2021 at 11:32
• Yes you're right, that was the mistake. (I had set $b_1=b_2=0$ in the pricing formula). Thanks a lot for your time explaining the answers and writing the code :) May 28, 2021 at 12:55
• No worries! I actually had fun!!! :D May 28, 2021 at 13:42

I believe something is wrong with your analytical pricing formula:

I have provided an R script for the analytical pricing formula specified on (p. 211) and it gives a call price of 0.2853.

library(pbivnorm)

maxassets_analytical <- function(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho){

y1 <- (log(S1/K) + (b1 + sigma1^2/2)*T)/ (sigma1*sqrt(T))

y2 <- (log(S2/K) + (b2 + sigma2^2/2)*T)/ (sigma2*sqrt(T))

sig <- sqrt(sigma1^2 + sigma2^2 - 2 * sigma1 * sigma2 * rho)

d <- (log(S1/S2) + (b1 - b2 + sig^2/2)*T)/ (sig*sqrt(T))

rho1 <- (sigma1 - rho * sigma2)/sig

rho2 <- (sigma2 - rho * sigma1)/sig

call <- S1 * exp((b1 - r) * T)  * pbivnorm(y1, d, rho1) +
S2 * exp((b2 - r) * T)  * pbivnorm(y2, -d + sig * sqrt(T), rho2) -
K * exp(-r * T) * (1 - pbivnorm(-y1+sigma1*sqrt(T), -y2 + sigma2 * sqrt(T), rho))

return(list("call price" = call, "y1" = y1, "y2" = y2))

}


Where $$b_1$$ and $$b_2$$ is the cost-of-carry for respectively, $$S_t^1$$ and $$S_t^2$$. For non-dividend paying stocks the cost of carry equals the risk-free rate. Inserting your specified parameters as-well as $$b_1=b_2=0.02$$ gives 0.2853:

For brevity, I have provided a picture of the formula below.

Here, $$M(a,b,\rho)$$ is the cumulative bivariate normal distribution function (as seen in Chapter 13).

Verification:

I verified my code by implementing the analytical pricing formula for the corresponding put option, since he provides an example calculation specified on (p. 212). Not only that, but the put option nests the analytical formula for the call option, which gives a way to verify it:

I implemented the put pricing function in R and the code can be found below. Following his example and thus inserting the parameters, $$S_1=100$$, $$S_2=105$$, $$X=98$$, $$T=0.5$$, $$r=0.05$$, $$b_1=-0.01$$, $$b_2=-0.04$$, $$\sigma_1 = 0.11$$, $$\sigma_2=0.16$$, $$\rho = 0.63$$, I get the following results:

Which agrees with the result in his example:

I hope this helps!

Code for the Put option:

maxassets_analyticalput <- function(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho){

sig <- sqrt(sigma1^2 + sigma2^2 - 2 * sigma1 * sigma2 * rho)

d <- (log(S1/S2) + ((b1 - b2) + (sig^2)/2)*T)/ (sig*sqrt(T))

cmax0 <- S2 * exp( (b2 - r) * T) + S1 * exp( (b1 - r) * T) * pnorm(d) - S2 * exp((b2
- r)*T) * pnorm(d-sig*sqrt(T))

cmax <- maxassets_analytical(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho)

pmax <- K * exp(-r * T) -
cmax0 + cmax[[1]]

return(list("put price" = pmax, "cmax, K = 0" = cmax0, "cmax" = cmax[[1]],
"d" = d, "sigma" = sig, "y1" = cmax[[2]], "y2" = cmax[[3]]))
}

• I was worried the analytic value indicated was not the right one when I tried the problem myself! Thanks for having clarified that! May 28, 2021 at 12:00
• This is incredibly helpful @Pleb, thank you! My misunderstanding came from the cost-of-carry rate. I set $b_1 = b_2 = 0$ in the pricing formula (gives a price of 0.22...), I did not realise that it should equal the constant short rate. Thank you for taking the time to explain it to me :) May 28, 2021 at 12:53