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)


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
  • $\begingroup$ Are dynamics without $dt$? Or is it a typo? $\endgroup$ Commented May 28, 2021 at 8:26
  • $\begingroup$ @YodaAndFriends Sorry, yes it was a typo, edited now $\endgroup$
    – R. Rayl
    Commented May 28, 2021 at 8:36

2 Answers 2


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

  • $\begingroup$ 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? $\endgroup$
    – R. Rayl
    Commented May 28, 2021 at 10:48
  • 1
    $\begingroup$ 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? $\endgroup$ Commented May 28, 2021 at 11:32
  • $\begingroup$ 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 :) $\endgroup$
    – R. Rayl
    Commented May 28, 2021 at 12:55
  • $\begingroup$ No worries! I actually had fun!!! :D $\endgroup$ Commented 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.


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


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]]))
  • 1
    $\begingroup$ I was worried the analytic value indicated was not the right one when I tried the problem myself! Thanks for having clarified that! $\endgroup$ Commented May 28, 2021 at 12:00
  • 1
    $\begingroup$ 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 :) $\endgroup$
    – R. Rayl
    Commented May 28, 2021 at 12:53

Your Answer

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

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