Spread Option Pricing on Single-Core and Parallel Computing Architectures

1. Introduction

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 N, 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.3 FFT algorithm

The FFT algorithm for spread option pricing along the line of [3] 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().

Table 1 shows the market parameters, and these inputs are taken from [5] where d1 and d2 represent the dividend rate for S1 and S2, respectively.

S1	S2	r	T	d1	d2
100	96	0.1	1.0	0.05	0.05

Table 1.

Market pricing parameters.

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 10 options, that is, NO=10. The pricing accuracy is gauged using the root mean square error (rmse):

Grid points: 2N	Pricing accuracy	Single core (s)	GPU parallel (s)	GPU speed factor
(a) 2d-GBM model: σ1=0.2,σ2=0.1,ρ=0.5
8	1.10E-04	2.671126	0.213765	12.49561902
9	4.29E-05	10.155822	0.254708	39.87241076
10	1.97E-05	39.884087	0.650358	61.32635718
11	1.06E-05	157.833762	2.422815	65.14478489
12	6.19E-06	317.858516	9.474695	33.54815284
13	3.80E-06	1316.938937	38.519949	34.18849119

(b) 2d-NIG model: μ1=0,μ2=0,α=6.20,δ=0.150β1=−3.80,β2=−2.50,ρ=0, μ1=0,μ2=0 and Δ=I
8	1.20E-04	2.528521	0.21573	11.72087536157
9	5.21E-05	9.768831	0.31099	31.41244621944
10	3.00E-05	39.372496	1.03477	38.04969824066
11	2.12E-05	162.069306	4.57217	35.44691939427
12	1.72E-05	651.850506	19.92393	32.71695934733
13	1.50E-05	2592.307476	83.77870	30.94232156861

Table 2.

Computational accuracy and processing times (a) 2d-GBM, (b) 2d-NIG.

rmse=∑j=1NOPFFTj−PMontejPMontej2,E10

where PMonte represents the benchmark price computed using Monte Carlo method with 1000000 simulations and 1000 time step and PFFT is the price from two-dimensional FFT.

From Table 2, we can see that the optimal values of N (in terms of computation efficiency) are in the middle of the tested range, 9 and 10 for both the 2d-GBM and 2d-NIG models. The decline in performance for larger values of N 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 ΦNIGuT, which involves two square root and exponential calculations, as opposed to simply one in the GBM model).

4. Conclusion

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 1000000 paths and 100 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.