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.

  • 1
    $\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
  • 1
    $\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

4 Answers 4


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

  • 1
    $\begingroup$ Regarding "bad answers". What do you think of JChoi's answer? And is Carr and Madan worth considering or not? $\endgroup$
    – nbbo2
    Dec 11, 2023 at 15:28
  • 1
    $\begingroup$ Carr and Madan becomes interesting when the assets do not follow a GBM but some process with a known characteristic function (affine stochastic volatility typically). There is no point to use it for GBMs. JChoi approximation is for a basket spread, it may be relevant in the case where the number of futures (or assets in general) is larger or equal to 3. This is not the context of the question however. $\endgroup$
    – jherek
    Dec 11, 2023 at 20:17

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$
    – jChoi
    Aug 4, 2018 at 0:18

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.


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


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.