A Fourier transform method for spread option pricing^††thanks: Research supported by the Natural Sciences and Engineering Research Council of Canada and MITACS, Mathematics of Information Technology and Complex Systems Canada

When $S_{jt},j=1,2$ are two asset price procsses, the basic spread option with maturity $T$ and strike $K\geq 0$ is the contract that pays $(S_{1T}-S_{2T}-K)^{+}$ at time $T$. From the risk-neutral expectation formula, the time $0$ price of this option, assuming a constant interest rate $r$, must be

The literature on applications of spread options is extensive, and is reviewed by Carmona and Durrleman [1] who explore further applications of spread options beyond the case of equities modelled by geometric Brownian motion, notably energy trading. For example, the “crack spread” is the difference between the prices of crude oil and natural gas. “Spark spreads” refer to differences between the price of electricity and the price of fuel: options on spark spreads are widely used by power plant operators to optimize their revenue streams. Energy pricing requires models with mean reversion and jumps very different from geometric Brownian motion, and pricing spread options in such situations can be challenging.

Closed formulas for (1) are known only for a limited set of cases. In the Bachelier stock model, $S_{t}=(S_{1t},S_{2t})$ is an arithmetic Brownian motion, and in this case (1) has a Black-Scholes type formula for any $T,K$. In the special case $K=0$ when $S_{t}$ is geometric Brownian motion, (1) is given by the Margrabe formula [8].

In the important case where $S_{t}$ is geometric Brownian motion and $K>0$, no explicit pricing formula is known. Therefore there is a long history of approximation methods. Numerical integration methods, typically Monte Carlo based, are often employed. Dempster and Hong [4] introduced a numerical integration method based on double fast Fourier transforms (FFT) that is efficient for geometric Brownian motion and certain more general asset price processes.

Another stream of research develops analytical methods applicable to log normal models that involve linear approximations of the nonlinear exercise boundary. Such methods are often very fast, but their accuracy is usually not easy to determine. Kirk [5] presented an analytical approximation that performs well in practice. Carmona and Durrleman [2] used a family of lower bounds for the spread option price to obtain an analytical approximation formula that apparently leads to acceptable accuracy. Moreover these results extend to approximate values for the Greeks. Deng, Li and Zhou [6] have developed a systematic analytical approximation scheme which is perhaps even more reliable.

Next to having analytical option pricing formulas, the fastest option pricing engines by numerical integration are usually those based on the fast Fourier transform methods introduced by Carr and Madan [3]. Their interest was in option pricing for geometric Lévy process models like the variance gamma (VG) model, but their basic framework has been adapted to a host of other models where the characteristic function of the returns is known.

The main purpose of the present paper is to give a new numerical integration method for computing spread options in two or higher dimensions via FFTs, applicable to any model for which the characteristic function of the joint return process is given analytically. Our method is completely different from prior approaches, and consequently has distinct advantages, notably the flexibility to include a wide range of desirable asset return models, with a competitive computational expense.

The results we describe all stem from the following fundamental new result that gives a Fourier representation of the basic spread option payoff function $P(x_{1},x_{2})=(e^{x_{1}}-e^{x_{2}}-1)^{+}$. Note that this spread option payoff is without loss of generality: We can reduce the general case $K\neq 0$ to $K=1$ by using scaling and interchange of $S_{1}$ and $S_{2}$.

Using this result, whose proof is given in the appendix, we will find we can follow the logic of Carr and Madan to derive numerical algorithms for efficient computation of a variety of spread options and their Greeks. The basic strategy to compute (1) is to combine (2) with an explicit formula for the characteristic function of the bivariate random variable $X_{t}=\log S_{t}$. For the remainder of this paper, we make a simplifying assumption that the dynamics of $X$ is homogeneous, which is equivalent to the following factorization of the characteristic function

Then, if $S_{t}=e^{X_{t}}$ is the price process for two traded assets, the spread option formula can be written

The Greeks are handled in exactly the same way. For example, the Delta $\Delta^{1}:=\partial\mbox{Spr}/\partial S_{10}$ is obtained as a function of $S_{0}$ by replacing $\Phi$ in (4) by $\partial\Phi/\partial S_{10}$.

Explicit double Fourier integrals like this can be approximated numerically by a two dimensional FFT. Such approximations involve both a truncation and discretization of the integral, and the two properties that determine their accuracy are the decay of the integrand of (4) in $u$-space, and the decay of the function Spr in $x$-space. The remaining issue of computing the gamma function is not a real difficulty. Fast and accurate computation of the complex gamma function in for example, Matlab, is based on the Lanczos approximation popularized by [9]¹¹1According to these authors, computing the gamma function becomes “not much more difficult than other built-in functions that we take for granted, such as $\sin x$ or $e^{x}$”..

In this paper, we demonstrate how our method performs in three different models for spread options on two stocks, namely the geometric Brownian motion (GBM), a three factor stochastic volatility (SV) model and the variance gamma (VG) model. Section 2 provides the essential definitions of the three types of models, including explicit formulas for their bivariate characteristic functions. Section 3 discusses how the two dimensional FFT can be implemented for our problem, and gives a heuristic picture of how the accuracy and speed will depend on the choices made in the implementation. Section 4 gives the detailed results of the performance of the method in the three models. In that section, the accuracy of each model is compared to benchmark values computed by an independent method for a reference set of option prices. We also demonstrate that the computation of the spread option Greeks in such models is equally feasible. Section 5 extends all the above results to spread options on three or more assets. While in such cases the formulation is simple, the resulting higher dimensional FFTs are in practice much slower to compute.

To compute (4) in these models we approximate the double integral by a double sum over the lattice

for appropriate choices of $N,\eta,\bar{u}:=N\eta/2.$ For the FFT it is convenient to take $N$ to be a power of $2$ and lattice spacing $\eta$ such that truncation of the $u$-integrals to $[-\bar{u},\bar{u}]$ and discretization leads to an acceptable error. Finally, we choose initial values $X_{0}=\log S_{0}$ to lie on the reciprocal lattice with spacing $\eta^{*}=2\pi/N\eta=\pi/\bar{u}$:

For any $S_{0}=e^{X_{0}}$ with $X_{0}=x(\ell)\in\Gamma^{*}$ we then have the approximation

Now, as usual for the discrete FFT, as long as $N$ is even,

This leads to the double inverse discrete Fourier transform

where

The selection of suitable values for $\epsilon,N$ and $\eta$ is a somewhat subtle issue when implementing the above FFT approximation of (4). We now briefly discuss separately the truncation error and discretization error. The pure truncation error, defined by taking $\eta\to 0,N\to\infty$ keeping $\bar{u}=N\eta/2$ fixed, can be made small if the integrand of (4) is small outside the square $[-\bar{u}+i\epsilon_{1},\bar{u}+i\epsilon_{1}]\times[-\bar{u}+i\epsilon_{2},\bar{u}+i\epsilon_{2}]$. Corollary 4, proved in the Appendix, gives an $O(|u|^{-2})$ upper bound on $\hat{P}$, while $\Phi(u)$ can generally been seen directly to have some $u$-decay. Typically, with some caveats, if one picks $\bar{u}$ large enough so that $|\Phi|<\delta_{1}\ll 1$ and has decay outside the square, then the truncation error will be less than $O(\delta_{1})$.

The pure discretization error, defined by taking $\bar{u}\to\infty,N\to\infty$ while keeping $\eta$ fixed, can be made small if $e^{\epsilon X_{0}^{\prime}}\mbox{Spr}$, taken as a function of $X_{0}\in\mathbb{R}^{2}$, has rapid decay in $X_{0}$. This is related to the smoothness of $\Phi(u)$ and the choice of $\epsilon$. The first two models are not very sensitive to $\epsilon$, but in the VG model the following conditions are needed to ensure that singularities in $u$ space are avoided:

The error in the approximation (11) with $N=\infty$ and $\eta$ finite, is given by the formula

Typically, with some caveats, the norm of this can be made smaller than $O(\delta_{2})\ll 1$ for $X_{0}\in[-c\bar{x},c\bar{x}]^{2}$ if $0<c\ll 1$ and $e^{\epsilon X_{0}^{\prime}}\mbox{Spr}(X_{0})$ is $O(\delta_{2})$ and rapidly decaying outside the square $[-\bar{x},\bar{x}]^{2}$.

In summary, one expects that the combined truncation and discretization error will be $O(\delta_{1}+\delta_{2})$ if $\bar{u}$ and $\eta=\pi/\bar{x}$ are each chosen as above.

The FFT method can also be applied to the Greeks, enabling us to tackle hedging and other interesting problems. It is particularly efficient for the GBM model, where differentiation under the integral sign is always permissible. For instance, the FFT formula for vega (the sensitivity to $\sigma$) takes the form:

where $\frac{\partial{\sigma^{2}}}{\partial{\sigma_{1}}}=[2\sigma_{1},0]$ and $\frac{\partial{\Sigma}}{\partial{\sigma_{1}}}=[2\sigma_{1},\rho\sigma_{2};\rho\sigma_{2},0]$. Other Greeks including those of higher orders can be computed in similar fashion. This method needs to be used with care for the SV and VG models, since it is possible that differentiation leads to an integrand that decays slowly.

Our numerical experiments were coded and implemented in Matlab version 7.6.0 on an Intel 2.80 GHz machine running under Linux with 1 GB physical memory. If they were coded in C++ with similar algorithms, we should expect to see much faster performance. All the following experiments were run with $\epsilon_{1}=-3,\epsilon_{2}=1$.

The strength of the FFT method is demonstrated by comparison with benchmarks. Based on a selection of $\{S_{10}^{i},S_{20}^{i},K^{i}\},i\in 1,2,3\dots n$, the objective function we measure is defined as

where $M^{i}$ and $B^{i}$ are FFT computed prices and benchmark prices for the $ith$ combination. We chose $\log{S_{10}}=\frac{i\pi}{10},\log{S_{20}}=-\frac{\pi}{5}+\frac{j\pi}{10},i,j\in{1,2,3,4,5,6}$ and fixed the strike $K=1$, leading to 36 prices contributing to the objective function. These choices cover a wide range of moneyness, from deep out-of-the-money to deep in-the-money. Since these combinations all lie on lattices $\Gamma^{*}$ corresponding to $N=2^{n}$ and $\bar{u}/10=2^{m}$ for integers $n,m$, all 36 prices $M^{i}$ can be computed simultaneously with one FFT.

Figure 1 shows how the FFT method performs in the 2-dimensional Geometric Brownian motion model for different choices of $N$ and $\bar{u}$. Since the two factors are bivariate normal, benchmark prices can be calculated to high accuracy by one dimensional integrations. In Figure 1 we can clearly see the effects of both truncation errors and discretization errors. For a fixed $\bar{u}$, the objective function decreases when $N$ increases. The $\bar{u}=20$ curve flattens out near $10^{-5}$ due to its truncation error of that magnitude. In turn, we can quantify its discretization errors with respect to $N$ by subtracting the truncation error from the total error. The flattening of the other three curves near $10^{-14}$ should be attributed to Matlab round-off errors: because of the rapid decrease of the characteristic function $\Phi$, their truncation error is negligible. For a fixed $N$, increasing $\bar{u}$ brings two effects: reducing truncation error and enlarging discretization error. These effects are well demonstrated in figure 1.

For the stochastic volatility model, no analytical or numerical method we know is consistently more accurate to serve as benchmark. Instead, we used Monte Carlo simulation to achieve a moderate degree of accuracy. Each benchmark price was created by $1,000,000$ simulations, each consisting of 2000 time steps. No variance reduction was employed. The resulting objective function is shown in Figure 2. The objective function only decreases to near $4\times 10^{-4}$, reflecting the inaccuracy of the benchmark. A comparable graph (not shown), using benchmark prices computed with the FFT method with $N=2^{12}$ and $\bar{u}=80$, showed similar behaviour to Figure 1, and is consistent with accuracies at the level of round-off. The computational cost to reduce the Monte Carlo simulation error to a comparable level would be prohibitive.

The VG model also lacks a reliable benchmark, so again we used Monte Carlo simulation. Since the VG process can be simulated as a Brownian motion time-changed by a gamma process, MC is here much more efficient than for the SV model, and we were able to run $10^{9}$ simulations with antithetic variates for each price. The objective function is shown in Figure 3. The truncation error for $\bar{u}=20$ is about $2\times 10^{-5}$. The other three curves flatten out near $7\times 10^{-6}$, which is presumably the accuracy of the benchmark. A comparable graph (not shown), using benchmark prices computed with the FFT method with $N=2^{12}$ and $\bar{u}=80$, showed similar behaviour to Figure 1, and is consistent with the FFT method being capable of producing accuracies at the level of round-off.

The strength of the FFT method is further illustrated by the computation of individual prices shown in Tables 1 to 3. The columns labeled “Analytic” and “Simulation” refer to the type of benchmark used. One can observe that an FFT with $N=256$ is capable of producing very high accuracy in all three models. In the GBM model, the errors are simply due to round-off. The discrepancies in the VG model are small too, less than 1 basis point. The apparent discrepancies in the SV model we attribute to the inaccurate benchmark, rather than the FFT method itself.

The FFT computes in a single iteration an $N\times N$ panel of prices $spread$ corresponding to initial values $S_{10}=e^{x_{10}+\ell_{1}\eta^{*}}$, $S_{20}=e^{x_{20}+\ell_{2}\eta^{*}}$, $K=1$, ($\ell_{1},\ell_{2})\in\{0,\dots,N-1\}^{2}$. If the desired selection of $\{S_{10},S_{20},K\}$ fits into this panel of prices, or its scaling, a single FFT suffices. If not, then one has to match ($x_{10},x_{20}$) with each combination, and run several FFTs, with a consequent increase in computation time. However, we have found that an interpolation technique is very accurate for practical purposes. For instance, prices for multiple strikes with the same $S_{10}$ and $S_{20}$ are approximated by a polynomial fit along the diagonal of the price panel: $\mbox{Spr}(S_{0};K_{1})=K_{1}\cdot spread(1,1)$, $\mbox{Spr}(S_{0};K_{1}e^{-\eta^{*}})=K_{1}e^{-\eta^{*}}\cdot spread(2,2)$, $\mbox{Spr}(S_{0};K_{1}e^{-2\eta^{*}})=K_{1}e^{-2\eta^{*}}\cdot spread(3,3)\dots$. The results of this technique are recorded in Table 2 and Table 3 in the column “Interpolation”. We can see this technique generates very competitive results and moreover, saves computational resource.

Finally, we computed first order Greeks using the method described at the beginning of Section 3 and compared them with finite differences. As seen in Table 4, the two methods come up with very consistent results. The Greeks of our at-the-money spread option exhibit some resemblance to the at-the-money European put/call option. The delta of $S_{1}$ is close to the delta of the call option, which is about 0.5. And the delta of $S_{2}$ is close to the delta of the put option, which is also about 0.5. The time premium of the spread option is positive. The option price is much more sensitive to $S_{1}$ volatility than to $S_{2}$ volatility. It is an important feature that the option price is negatively correlated with the underlying correlation: Intuitively speaking, if the two underlyings are strongly correlated, their co-movements diminish the probability that $S_{1T}$ develops a wide spread over $S_{2T}$. This result is consistent with observations made by [6].

Since the FFT method naturally generates a panel of prices, and interpolation can be implemented accurately with negligible additional computational cost, it is appropriate to measure the efficiency of the method by timing the computation of a panel of prices. Such computing times are shown in Table 5. For the FFT method, the main computational cost comes from the calculation of the matrix $H$ in (11) and the subsequent FFT of $H$. We see that the GBM model is noticeably faster than the SV and VG models: This is due to a recursive method used to calculate the $H$ matrix entries of the GBM model, which is not applicable for the SV and VG models. The number of calculations for $H$ is of order$N^{2}$ which for large $N$ exceeds the $N\log{N}$ of the FFT of $H$, and thus the advantage of this efficient algorithm for GBM is magnified as $N$ increases. However, our FFT method is still very fast for the SV and VG models and is able to generate a large panel of prices within a couple of seconds.

The ideas of section 2 turn out to extend naturally to basket options with payoffs $(\tilde{S}_{T}-S_{1T}-\dots-S_{MT}-1)^{+}$ for $M\geq 2$. The analysis is based on a simple result proved in the Appendix:

As before we write $x_{m}=\log S_{m},\tilde{x}=\log\tilde{S}$ and seek to compute for $u\in\mathbb{C}^{M}$, $\tilde{u}\in\mathbb{C}$

We introduce the factor $1=\int_{\mathbb{R}}\delta(e^{z}-\sum_{m=1}^{M}e^{x_{m}})e^{z}dz$ and interchange the $z$ integral with the $x$ integrals. Then using Proposition 2 one finds

When we can apply Proposition 1, we obtain the result:

This paper presents a fundamentally new approach to the valuation of spread options, an important class of financial contracts. The method is based on a newly discovered explicit formula for the Fourier transform of the spread option payoff in terms of the gamma function.

This mathematical result leads to simple and transparent algorithms for computing spread options in all dimensions. The powerful tool of the Fast Fourier Transform provides an accurate and efficient implementation of the fundamental result. The difficulties and pitfalls of the FFT, of which there are admittedly several, are by now well understood, and thus the reliability and stability properties of our method are rather clear.

Many important processes in finance, particularly affine models and Lévy jump models, have well known explicit characteristic functions, and can be included in the method with little difficulty. Thus the method can be easily applied to important problems arising in energy and commodity markets.

Finally, the Greeks can be systematically evaluated for such models, with similar performance and little extra work.