Market pricing parameters.
This paper introduces parallel computation for spread options using two-dimensional Fourier transform. Spread options are multi-asset options whose payoffs depend on the difference of two underlying financial securities. Pricing these securities, however, cannot be done using closed-form methods; as such, we propose an algorithm which employs the fast Fourier Transform (FFT) method to numerically solve spread option prices in a reasonable amount of short time while preserving the pricing accuracy. Our results indicate a significant increase in computational performance when the algorithm is performed on multiple CPU cores and GPU. Moreover, the literature on spread option pricing using FFT methods documents that the pricing accuracy increases with FFT grid size while the computational speed has opposite effect. By using the multi-core/GPU implementation, the trade-off between pricing accuracy and speed is taken into account effectively.
- spread option
- single core
- parallel computing
Spread options have widespread uses across many industries, remarkably in the seasonal commodity market futures. One notable example is the crack spread, which is the pricing difference between a barrel of crude oil and the petroleum products refined from it. Traditionally, crack spreads involve the purchase of oil futures and the simultaneous sale of futures of the refined product, whether it be gasoline, heating oil, or other similar products. Refiners seek a positive spread between the prices of crude oil and refined products, meaning that the price of the input (oil in this case) is lower than the price of the output (gasoline, kerosene, etc.).
Beyond the oil industry, spread options are utilized by suppliers in industries as disparate as the soybean and electricity markets. Soybean spread options are known as crush spreads, and electricity spread options are known as spark spreads. Similar to crack spreads, both these options seek to maximize the spread between input costs and output prices for suppliers to maximize profit.
The use of spread options across such disparate fields is but a testament to their widespread use. Considering the popularity of spread options, it thus prompts the needs of accurate and fast pricing of price of these options. This paper dwells on the discussion of fast algorithms for these options.
As mentioned previously, pricing spread options tends to be something that cannot be easily performed using traditional closed-form solutions. Therefore, we explore the utilization of the fast Fourier transform method to price these securities. Specifically, we perform the FFT on the spread options payoff, assuming knowledge of the model joint characteristic function, which we represent as a pointwise multiplication of the characteristic function and the complex gamma function in the Fourier domain.
Regarding our implementation, we adapt the parallel computing Toolbox in MATLAB to take advantage of the multi-core capabilities of GPU processing, to substantially improve the performance and computational efficiency of the algorithm for spread options. This methodology may serve a great deal especially in model calibration and risk management approach. For measuring performance, we only time the execution of the inverse Fast Fourier Transform, as the prior steps are merely initialization of the necessary arrays. For measuring accuracy, we compute the Euclidean norms of each of the resulting option price arrays from each implementation (single-core and GPU), and find the percent error between the norms as follows:
The main contribution of this paper is the use of parallel computation for the spread option value using two-dimensional Fourier transform. Our implementation is developed both for single-core processors, as well as for parallel processing on multi-core/GPU systems. For both execution methods, we have implemented our algorithm using MATLAB built-in functions to produce a version of the algorithm compatible for multi-core systems. To ascertain the impact of the different environments as well as the different methods of execution on the computational efficiency of the algorithm, we record the times of execution for different values of FFT grid size , which is the number of grid points of discretization of the characteristic function along the two asset dimensions, for both the classical single-core implementation and the multi-core/GPU implementation. This approach completely eliminates the trade off between computational accuracy and speed, that is, we price spread option accurately and in a fastest possible way.
2. Model description
2.1 Spread option valuation
We fix the trading time horizon and consider a filtered probability space defined in the usual way.
The goal is to compute the value of an European spread option between stock price processes and with maturity time and exercise price .
At expiration time , the payoff is given by
and under a risk-neutral conditional measure1 its value at time is given by
Under normality assumption, Eq. (3) can be analytically approximated as done in . However, a departure from normality assumption ushers into numerical computational difficulties. In option pricing literature of finance, Fourier transform-based method is usually the best candidate to approximate the solution for Eq. (3), in this case whenever the joint characteristic function of the asset price processes, and are available. For spread options pricing valuation methodology using two-dimensional fast Fourier transform (FFT) techniques was coined in [2, 3]. We cite the important formula from  that gives the price for spread option
where , is the joint characteristic function of the log return.
Fast and accurate pricing is often the most desirable feature of the model. In , the authors consider spread options pricing in C++ using fast Fourier transform in the west (FFTW). They observed the trade off between fast computation and numerical accuracy; pricing accuracy is monotonic in the number of FFT grid size used in the price computation. However, using a large number of FFT grid size slow down the speed of price computation.
2.2 Model characteristic function
Fast Fourier transform (FFT) method is generically applicable in finance because it only requires the specification of the characteristic function of the random variable. In terms of spread options, one just need a characteristic function of the joint distribution of the financial variables in question. Here, we employ two characteristic functions: one based on two-dimensional normal distribution and the other one based on two-dimensional normal inverse Gaussian (NIG) distribution.
2.2.1 Two-dimensional geometric Brownian motion (GBM)
The characteristic function for a spread option comprised of two assets, each of which is modeled as a correlated GBM, is given by
where . denote the risk-free rate, volatilities, respectively, and is the correlation parameter between two asset prices processes and . Here, is the covariance matrix.
2.2.2 Two-dimensional normal inverse Gaussian (NIG) Levy process
Let denote a two-dimension NIG random variable. The characteristic function of is given by
where , , and is a symmetric, positive-definite matrix. Moreover, the structural matrix is assumed to have determinant .
The covariance matrix corresponding to the two-dimensional NIG-distributed random variable is
2.3 FFT algorithm
The FFT algorithm for spread option pricing along the line of  can be described as follows
3. Numerical results
3.1 Implementation outlook
As mentioned earlier, two versions of the algorithm were programmed in MATLAB, namely a single-core variant and a multi-core GPU variant. In MATLAB, the Parallel Processing Toolbox was used to exploit multi-core GPU capabilities to run the algorithm. Among its capabilities is the ability to run for loops and perform array operations in parallel, both on multi-core CPUs as well as GPUs.
As in the Algorithm 1 given above, the for loop in lines 3–6 was run in parallel, across six CPU cores, employing the parfor directive available in the MATLAB Parallel Processing Toolbox. We also sought to run computationally heavy functions on the GPU we had available, to improve the efficiency of our algorithm beyond what would be possible on a multi-core CPU. In that regard, we executed the inverse FFT, as described in line 7, on the GPU. We accomplished this by copying our H and A arrays onto the GPU, such that any further processing of those arrays would only occur on the GPU. To perform this operation, we utilized the MATLAB inbuilt function gpuArray() and copied the two aforementioned arrays to the GPU after the for loop. To transfer the GPU results (following execution of the inverse FFT) back to the local workspace, we used the MATLAB function gather().
The computer used to produce the following results was an ASUS ROG Strix Scar II GL704GW, with an Intel Core i7-8750H processor clocked at 2.20 Hz and comprising of 6 cores, 16GB RAM, and an NVIDIA GeForce RTX 2070 GPU with 8GB memory, running Windows 10. The computational times of the algorithm are tabulated in Table 2 and Figure 1. In a single run, we compute the price of options, that is, . The pricing accuracy is gauged using the root mean square error (rmse):
|Grid points:||Pricing accuracy||Single core (s)||GPU parallel (s)||GPU speed factor|
|(a) 2d-GBM model:|
|(b) 2d-NIG model: , and|
where represents the benchmark price computed using Monte Carlo method with simulations and time step and is the price from two-dimensional FFT.
From Table 2, we can see that the optimal values of (in terms of computation efficiency) are in the middle of the tested range, and for both the 2d-GBM and 2d-NIG models. The decline in performance for larger values of is due to the increased memory requirements. When compared to 2d-GBM, 2d-NIG seems to have less of an increase in performance when executed on GPU, which could be because it is more computationally heavy (such as calculating the characteristic function , which involves two square root and exponential calculations, as opposed to simply one in the GBM model).
In this work, we built on the literature on fast and accurate pricing of spread options based on two-dimensional FFT method using parallel computation. We examined the effectiveness of this approach by comparing the computational times of CPU and GPU implementations of the FFT Spread Option Pricing Algorithm in MATLAB. We have taken benchmarks prices from Monte Carlo simulations with paths and discretization time steps. Our results decisively conclude that the execution of the algorithm on a GPU significantly improves computational performance, decreasing the time taken to run by a factor of up to almost 60x. Considering how common spread options are in the financial market, a faster way to price these securities means increased efficiency in transactions involving spread options, and the FFT algorithm implemented for this project also vastly improves the accuracy of spread option pricing. This approach is very useful to accurate calibration of spread options which is recognized to be a challenging exercise.
As an extension to this work, one could develop a 3-asset spread option pricing algorithm using the 3D Fast Fourier Transform Algorithm. Such a scheme, while computationally heavy, could be rendered more efficient by harnessing the power of GPUs through the tools available in MATLAB.
Kirkpatrick S, Gelatt C, Vecchi M. Optiminization by simulated annealing. Science. 1983; 1(23):671-680
Dempster M, Hong S. Spread option valuation and the fast Fourier transform. In: Mathematical Finance-Bachelier Congress. Vol. 2000. 2002. pp. 203-220
Hurd T, Zhou Z. A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics. 2010; 1:142-157
Alfeus M, Schlögl E. On spread option pricing using two-dimensional Fourier transform. International Journal of Theoretical and Applied Finance. 2019; 22(5)
Hurd T, Zhou Z. A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics. 2010; 1( 1):142-157
- 3 Specifically: Under the risk–neutral measure associated with taking the continuously compounded savings account as the numeraire, and (for expositional simplicity) assuming a constant interest rate r.