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} = \max( F_1(S,T) - F_2(S,T), 0 )$ at some expiry date $S \leq T$?

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

  • $\begingroup$ Might be simpler to model the spread directly as a normal process, then price with a Bachelier model? Otherwise you could check Haug's book, I don't have it handy... $\endgroup$ Dec 22, 2014 at 18:53
  • $\begingroup$ Do you mean that you want to price payoffs $\max( F_1(S,T) - F_2(S,T), 0 )$ at $S$ or European payoffs $g(X_{1,2})$ (for example, if $g(x) = (x-K)_+$, that would be $\max( F_1(S,T) - F_2(S,T) - K, 0)$). In the first case you have an exchange option and this reduces to the BS formula by using the second underlying as a numeraire. In the second case, there is no closed form. $\endgroup$
    – AFK
    Dec 24, 2014 at 13:52

3 Answers 3


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.

  • $\begingroup$ Actually the links in your post are reversed... The first one is Choi and the second is Li Deng Zhou. $\endgroup$
    – nbbo2
    Feb 17, 2017 at 17:10
  • 1
    $\begingroup$ Fixed and updated with better links. $\endgroup$ Aug 4, 2018 at 0:18

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) {
      } 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!!!


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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