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

Question

I want to compute the price of the option with payoff \begin{equation} \max \big\{\max\{S^1_T, S^2_T\} - K, 0\big\}, \end{equation} 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

Answer 1

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

Answer 2

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]]))
}