Pricing Futures Spread Options with Margrabe, Integration, and Approximations

# How to price a futures spread option?


# How to price a futures spread option?












Let's say I have two futures contract $F_1(0,T)$ and $F_2(0,T)$ on two different correlated underlyings.

If I assume that both underlying follow a GBM with volatility $\sigma_1$ and $\sigma_2$ respectively, and their return are correlated by $\rho$, is there a closed-form formula I can use to price an European option on the spread $X_{1,2} = F_1(S,T) - F_2(S,T)$ at some expiry date $S \leq T$?

Otherwise, how do we usually proceed to price these things? Monte-Carlo using risk-neutral dynamics and discounting payoffs?

The payoff of such an option at maturity reads $\max(X_{1,2}-K,0)$ where $K \in \mathbb{R}$ is the strike price. In particular we have $K=0$ for an at-the-money option.

## Answer by jherek (score 1)

https://quant.stackexchange.com/a/77678

Surprised to see so many bad answers here. The comments on the question are good however.

With strike=0, the pricing formula is known as the Margrabe formula, it is Black-Scholes using one of the futures price as numeraire.

As mentioned by @AFK, with a non-zero strike, the correct way to do this is a one-dimensional numerical integration: the double integral on the payoff $max(F1-F2-K,0)$ over the joint density reduces to a single integral on the Black-Scholes formula. The reference is:

Neil D Pearson. An efficient approach for pricing spread options. The Journal of Derivatives, 3(1):76–91, 1995.

It is exact, fast. The Kirk approximation may be faster, but being an approximation, will behave badly on some cases.

More generally, you may want to model the spread directly by a normal (Bachelier) model. The issue is then to go back to exposures on F1 and F2 (and not just F1-F2).

## Answer by jaehyukchoi49 (score 0)

https://quant.stackexchange.com/a/32492

The best approximation that I know is Li, Deng and Zhou, 2006. It's an analytic approximation where the price is expressed as a direct formula, so easy to implement.

If you want to be VERY accurate, here's my paper, J Choi (2018) (Arxiv). It handles the options on any linear combination of assets such as basket and Asian options as well as spread option. There is some discussion on the performances of other analytic approximations.

## Answer by Rusan Kax (score -1)

https://quant.stackexchange.com/a/15965

You could take a look at this, if you have not already Fast Fourier transform for pricing spread options

The classic reference on this is a paper by Carr and Madan, reference number 4 in that document.

It is not a closed form solution, per se, but the method is quick and easy to implement.

## Answer by Matt (score -1)

https://quant.stackexchange.com/a/70482

As Jaehyuk Choi mentioned, he has a very accurate spread option + basket + Asian option pricer. What he didn't mention is that he already coded it up in R for you here: https://github.com/PyFE/SumBSM-R/blob/master/SumBSM/blksmd.R

He noted the Li, Deng, Zhou 2006 as an accurate easy one to code up. He coded that up already (for calls), and I just added the put-call parity to calculate put values (note it gives erroneous results with negative strikes, maybe someone can code the actual formula as I can't access the paper). You pass spot like: `spot <- c(spot1, spot2)` and vol like `vol <- c(vol1, vol2)` - I'm not an R user at all, but still can follow the code a bit. There are some boundary "break" conditions that aren't handled by the code, but at least this gives you a head start:

```
blksmd_CalcSpreadOptLDZ <- function( strk, spot, t.exp, vol, rho, callput, r = 0, d = 0 ){
  call <- rep(NA, length(strk))
  fwd <- spot * exp((r-d)*t.exp)
  std <- vol * sqrt(t.exp)
  
  rho.comp <- sqrt(1-rho*rho)*std[1]
  mu1 <- log(fwd[1]) - 0.5*std[1]*std[1]
  mu2.exp <- fwd[2] * exp(-0.5*std[2]*std[2]) # R = exp(mu2) with y0 = 0

  # vector for the rest
  r.plus.k <- mu2.exp + strk # (R+K) in DLZ (vectorized)

  epsilon <- -1/(2*rho.comp) * (std[2]*std[2]*mu2.exp*strk)/(r.plus.k*r.plus.k)
  
  C3 <- ( mu1 - log(r.plus.k) )/rho.comp
  D3 <- ( rho*std[1] - std[2]*mu2.exp/r.plus.k )/rho.comp
  
  C2 <- C3 + std[2]*( D3 + epsilon*std[2] )
  D2 <- D3 + 2*std[2]*epsilon
  
  C1 <- C3 + rho*std[1]*(D3 + epsilon*rho*std[1]) + rho.comp
  D1 <- D3 + 2*rho*std[1]*epsilon

  I_S1 <- blksmd_CalcSpreadOptLDZ_I( C1, D1, epsilon )
  I_S2 <- blksmd_CalcSpreadOptLDZ_I( C2, D2, epsilon )
  I_K  <- blksmd_CalcSpreadOptLDZ_I( C3, D3, epsilon )

  price.fwd <- fwd[1]*I_S1 - fwd[2]*I_S2 - strk*I_K

  if (callput) {
      return(exp(-r*t.exp)*price.fwd)
      } else {
      return((exp(-r*t.exp)*price.fwd) + exp(-r*t.exp)*strk - (spot[1]-spot[2]))
      }  
}

#' An auxilary function used by blksmd_CalcSpreadOptLDZ
#' 
blksmd_CalcSpreadOptLDZ_I <- function( u, v, eps ) {
  u2 <- u*u
  u4 <- u2*u2
  v2 <- v*v
  v4 <- v2*v2
  v6 <- v4*v2
  v2.sqrt <- sqrt(1+v2)
  
  arg.norm <- u/v2.sqrt
  J0 <- pnorm(arg.norm)
  J1 <- (1+(1+u2)*v2)/v2.sqrt^5L * dnorm(arg.norm)
  J2 <- (6*(1-u2)*v2 + (21-2*u2-u4)*v4 + 4*(3+u2)*v6 - 3)*u/v2.sqrt^11L * dnorm(arg.norm)
  
  return( J0 + eps*(J1 + 0.5*eps*J2))
}
```

So this is his code (with one addition for puts), you should upvote his answer if you like it. He also has his very complicated procedure coded up as well at the above link, thank you Jaehyuk!!!