**Abstract**

This paper introduces parallel computation for spread options using twodimensional 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.

**Keywords:** spread option, single core, parallel computing

### **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:

$$\frac{\text{GPU norm} - \text{single core norm}}{\text{single core norm}} \times 100\tag{1}$$

S*t*ð Þ¼ *S*1, *S*2, *K e*

*DOI: http://dx.doi.org/10.5772/intechopen.93430*

S*t*ð Þ¼ *S*1, *S*2, *K e*

where *<sup>X</sup>*<sup>0</sup> <sup>¼</sup> log *<sup>S</sup>*, *<sup>Φ</sup>*ð Þ¼ *<sup>u</sup>*, *<sup>T</sup> E e*<sup>i</sup>*u X*<sup>0</sup>

down the speed of price computation.

**2.2 Model characteristic function**

the log return.

<sup>¼</sup> <sup>1</sup> ð Þ <sup>2</sup>*<sup>π</sup>* <sup>2</sup> *<sup>e</sup>*

*Spread Option Pricing on Single-Core and Parallel Computing Architectures*

Γð Þ¼ *z*

�*r T*ð Þ �*<sup>t</sup> <sup>E</sup><sup>Q</sup>* ð Þ *<sup>S</sup>*1ð Þ� *<sup>T</sup> <sup>S</sup>*2ð Þ� *<sup>T</sup> <sup>K</sup>* þjF*<sup>t</sup>*

�*r T*ð Þ �*<sup>t</sup> EQ* ð Þ *<sup>S</sup>*1ð Þ� *<sup>T</sup> <sup>S</sup>*2ð Þ� *<sup>T</sup> <sup>K</sup>* þjF*<sup>t</sup>*

*R*þiϵ *e* i*uX*0

� �

<sup>0</sup>Φð Þ *<sup>u</sup>*, *<sup>T</sup> P u* ^ð Þ*d*<sup>2</sup>

h i, is the joint characteristic function of

<sup>Γ</sup>ð Þ <sup>i</sup>*u*<sup>1</sup> <sup>þ</sup> <sup>1</sup> (5)

*<sup>z</sup>*�<sup>1</sup> *dt*, (6)

*u*

(4)

Under normality assumption, Eq. (3) can be analytically approximated as done in [1]. 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, *S*<sup>1</sup> and *S*<sup>2</sup> 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 [3] that gives the price for spread option

�*rT* ð ð

*P u* ^ð Þ¼ <sup>Γ</sup>ð Þ <sup>i</sup>ð Þ� *<sup>u</sup>*<sup>1</sup> <sup>þ</sup> *<sup>u</sup>*<sup>2</sup> <sup>1</sup> <sup>Γ</sup>ð Þ �i*u*<sup>2</sup>

ð<sup>∞</sup> 0 *e* �*t t*

*<sup>T</sup>*�*X*<sup>0</sup> ð Þ<sup>0</sup>

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.

The characteristic function for a spread option comprised of two assets, each of

<sup>1</sup> *σ*1*σ*2*ρ σ*1*σ*2*ρ σ*<sup>2</sup>

!

*r*, *σ<sup>i</sup>* : *i* ¼ 1, 2 denote the risk-free rate, volatilities, respectively, and *ρ* is the

2

2 � �<sup>0</sup>

!

� *<sup>u</sup>*Σ*GBMu*<sup>0</sup>

2

*T*

, and *<sup>σ</sup>*<sup>2</sup> <sup>¼</sup> diagð Þ <sup>Σ</sup>*GBM* , i <sup>¼</sup> ffiffiffiffiffiffi

(7)

�<sup>1</sup> <sup>p</sup> , .

*2.2.1 Two-dimensional geometric Brownian motion (GBM)*

*<sup>Φ</sup>*GBMð Þ¼ *<sup>u</sup>*; *<sup>T</sup>* exp i*u rTe* � *<sup>σ</sup>*<sup>2</sup>*<sup>T</sup>*

which is modeled as a correlated GBM, is given by

where *<sup>e</sup>* <sup>¼</sup> ½ � 1, 1 , <sup>Σ</sup>*GBM* <sup>¼</sup> *<sup>σ</sup>*<sup>2</sup>

**95**

Fast and accurate pricing is often the most desirable feature of the model. In [4], 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

� �*:* (3)

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.1 Spread option valuation**

We fix the trading time horizon *T* and consider a filtered probability space *Ω*, F, *P*,ð Þ F*<sup>t</sup> <sup>t</sup>*<*<sup>T</sup>* <sup>&</sup>lt;*<sup>T</sup>* defined in the usual way.

The goal is to compute the value of an European spread option between stock price processes *S*<sup>1</sup> ¼ f g *S*1ð Þ*t* <sup>0</sup>≤*t*≤*<sup>T</sup>* and *S*<sup>2</sup> ¼ f g *S*2ð Þ*t* <sup>0</sup>≤*t*≤*<sup>T</sup>* with maturity time *T* and exercise price *K*.

At expiration time *T*, the payoff is given by

$$\mathcal{S}\_T(\mathbb{S}\_1, \mathbb{S}\_2, K) = \max\left\{ \mathbb{S}\_1(T) - \mathbb{S}\_2(T) - K, 0 \right\},\tag{2}$$

and under a risk-neutral conditional measure<sup>1</sup> *Q* its value at time *t* is given by

<sup>1</sup> 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.

*Spread Option Pricing on Single-Core and Parallel Computing Architectures DOI: http://dx.doi.org/10.5772/intechopen.93430*

$$\mathcal{S}\_t(\mathbb{S}\_1, \mathbb{S}\_2, K) = e^{-r(T-t)} E^Q \left[ (\mathbb{S}\_1(T) - \mathbb{S}\_2(T) - K)^+ | \mathcal{F}\_t \right]. \tag{3}$$

Under normality assumption, Eq. (3) can be analytically approximated as done in [1]. 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, *S*<sup>1</sup> and *S*<sup>2</sup> 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 [3] that gives the price for spread option

$$\begin{split} \mathcal{S}\_{t}(\mathcal{S}\_{1}, \mathcal{S}\_{2}, K) &= e^{-r(T-t)} E^{Q} \left[ \left( \mathcal{S}\_{1}(T) - \mathcal{S}\_{2}(T) - K \right)^{+} | \mathcal{F}\_{t} \right] \\ &= \frac{1}{\left( 2\pi \right)^{2}} e^{-rT} \int \int\_{R + i\epsilon} e^{i\omega \mathcal{X}\_{0}'} \Phi(u, T) \hat{P}(u) d^{2}u \end{split} \tag{4}$$

$$\hat{P}(\boldsymbol{\mu}) = \frac{\Gamma(\mathbf{i}(\boldsymbol{\mu}\_1 + \boldsymbol{\mu}\_2) - \mathbf{1})\Gamma(-\mathbf{i}\boldsymbol{\mu}\_2)}{\Gamma(\mathbf{i}\boldsymbol{\mu}\_1 + \mathbf{1})} \tag{5}$$

$$
\Gamma(x) = \int\_0^\infty e^{-t} t^{x-1} dt,\tag{6}
$$

where *<sup>X</sup>*<sup>0</sup> <sup>¼</sup> log *<sup>S</sup>*, *<sup>Φ</sup>*ð Þ¼ *<sup>u</sup>*, *<sup>T</sup> E e*<sup>i</sup>*u X*<sup>0</sup> *<sup>T</sup>*�*X*<sup>0</sup> ð Þ<sup>0</sup> h i, is the joint characteristic function of

the log return.

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

*Real Perspective of Fourier Transforms and Current Developments in Superconductivity*

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:

GPU norm � single core norm

price spread option accurately and in a fastest possible way.

**2. Model description**

*Ω*, F, *P*,ð Þ F*<sup>t</sup> <sup>t</sup>*<*<sup>T</sup>* <sup>&</sup>lt;*<sup>T</sup>*

exercise price *K*.

**94**

**2.1 Spread option valuation**

defined in the usual way.

At expiration time *T*, the payoff is given by

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

We fix the trading time horizon *T* and consider a filtered probability space

The goal is to compute the value of an European spread option between stock price processes *S*<sup>1</sup> ¼ f g *S*1ð Þ*t* <sup>0</sup>≤*t*≤*<sup>T</sup>* and *S*<sup>2</sup> ¼ f g *S*2ð Þ*t* <sup>0</sup>≤*t*≤*<sup>T</sup>* with maturity time *T* and

and under a risk-neutral conditional measure<sup>1</sup> *Q* its value at time *t* is given by

<sup>1</sup> 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.

S*T*ð Þ¼ *S*1, *S*2, *K* max f g *S*1ð Þ� *T S*2ð Þ� *T K*, 0 , (2)

single core norm � <sup>100</sup> (1)

Fourier domain.

Fast and accurate pricing is often the most desirable feature of the model. In [4], 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

$$\phi\_{\rm GBM}(u;T) = \exp\left(\mathrm{i}u\left(rTe - \frac{\sigma^2 T}{2}\right)' - \frac{u\Sigma\_{\rm GBM} u'T}{2}\right) \tag{7}$$

where *<sup>e</sup>* <sup>¼</sup> ½ � 1, 1 , <sup>Σ</sup>*GBM* <sup>¼</sup> *<sup>σ</sup>*<sup>2</sup> <sup>1</sup> *σ*1*σ*2*ρ σ*1*σ*2*ρ σ*<sup>2</sup> 2 !, and *<sup>σ</sup>*<sup>2</sup> <sup>¼</sup> diagð Þ <sup>Σ</sup>*GBM* , i <sup>¼</sup> ffiffiffiffiffiffi �<sup>1</sup> <sup>p</sup> , . correlation parameter between two asset prices processes *S*<sup>1</sup> and *S*2. Here, Σ*GBM* is the covariance matrix.

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

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

*d*<sup>1</sup> and *d*<sup>2</sup> represent the dividend rate for *S*<sup>1</sup> and *S*2, respectively.

*Spread Option Pricing on Single-Core and Parallel Computing Architectures*

*DOI: http://dx.doi.org/10.5772/intechopen.93430*

rmse <sup>¼</sup> <sup>X</sup> *NO*

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

two-dimensional FFT.

*Market pricing parameters.*

**Table 1.**

**Table 2.**

**97**

*j*¼1

*P j FFT* � *<sup>P</sup> <sup>j</sup>*

where *PMonte* represents the benchmark price computed using Monte Carlo method with 1000000 simulations and 1000 time step and *PFFT* is the price from

*S***<sup>1</sup>** *S***<sup>2</sup>** *rTd***<sup>1</sup>** *d***<sup>2</sup>** 100 96 0.1 1.0 0.05 0.05

**Grid points: 2***<sup>N</sup>* **Pricing accuracy Single core (s) GPU parallel (s) GPU speed factor (a) 2d-GBM model:** *σ***<sup>1</sup>** ¼ **0***:***2,** *σ***<sup>2</sup>** ¼ **0***:***1,** *ρ* ¼ **0***:***5** 1.10E-04 2.671126 0.213765 12.49561902 4.29E-05 10.155822 0.254708 39.87241076 1.97E-05 39.884087 0.650358 61.32635718 1.06E-05 157.833762 2.422815 65.14478489 6.19E-06 317.858516 9.474695 33.54815284 3.80E-06 1316.938937 38.519949 34.18849119 **(b) 2d-NIG model:** *μ***<sup>1</sup>** ¼ **0,** *μ***<sup>2</sup>** ¼ **0,** *α* ¼ **6***:***20,** *δ* ¼ **0***:***150** *β***<sup>1</sup>** ¼ �**3***:***80,** *β***<sup>2</sup>** ¼ �**2***:***50,** *ρ* ¼ **0,** *μ***<sup>1</sup>** ¼ **0,** *μ***<sup>2</sup>** ¼ **0 and Δ** ¼ I 1.20E-04 2.528521 0.21573 11.72087536157 5.21E-05 9.768831 0.31099 31.41244621944 3.00E-05 39.372496 1.03477 38.04969824066 2.12E-05 162.069306 4.57217 35.44691939427 1.72E-05 651.850506 19.92393 32.71695934733 1.50E-05 2592.307476 83.77870 30.94232156861

*P j Monte* !<sup>2</sup>

*Monte*

, (10)

### *2.2.2 Two-dimensional normal inverse Gaussian (NIG) Levy process*

Let *S* denote a two-dimension NIG random variable. The characteristic function of *S* ¼ ð Þ *S*1, *S*<sup>2</sup> is given by

$$\Phi\_{\rm NIG}(u;T) = \exp\left(\mathrm{i}u'\mu T + \delta T \left[\sqrt{a^2 - \beta\Delta\beta'} - \sqrt{a^2 - (\beta + \mathrm{i}u)\Delta(\beta + \mathrm{i}u)'}\right]\right), \tag{8}$$

where *<sup>α</sup>*, *<sup>δ</sup>*∈*R*þ, *<sup>β</sup>*, *<sup>μ</sup>*∈*R*<sup>2</sup> , and *Δ*∈*R*2�<sup>2</sup> is a symmetric, positive-definite matrix. Moreover, the structural matrix *Δ* is assumed to have determinant Detð Þ¼ *Δ* 1.

The covariance matrix corresponding to the two-dimensional NIG-distributed random variable *S* is

$$
\Sigma\_{\rm NIG} = \delta \left( a^2 - \beta \Delta \beta' \right)^{-\frac{1}{2}} \left( \Delta + \left( a^2 - \beta \Delta \beta' \right)^{-1} \Delta \beta \beta' \Delta \right) \tag{9}
$$

### **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

*Spread Option Pricing on Single-Core and Parallel Computing Architectures DOI: http://dx.doi.org/10.5772/intechopen.93430*

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 *d*<sup>1</sup> and *d*<sup>2</sup> represent the dividend rate for *S*<sup>1</sup> and *S*2, respectively.

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

$$\text{rmsse} = \sum\_{j=1}^{No} \left( \frac{P\_{\text{FFT}}^j - P\_{\text{Moutte}}^j}{P\_{\text{Moutte}}^j} \right)^2,\tag{10}$$

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.


**Table 1.** *Market pricing parameters.*

correlation parameter between two asset prices processes *S*<sup>1</sup> and *S*2. Here, Σ*GBM* is

*Real Perspective of Fourier Transforms and Current Developments in Superconductivity*

Let *S* denote a two-dimension NIG random variable. The characteristic function

� � � � q

<sup>2</sup> *<sup>Δ</sup>* <sup>þ</sup> *<sup>α</sup>*<sup>2</sup> � *βΔβ*<sup>0</sup> � ��<sup>1</sup>

� �

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *α*<sup>2</sup> � ð Þ *β* þ i*u* Δð Þ *β* þ i*u* <sup>0</sup>

> Δ*ββ*<sup>0</sup> *Δ*

, and *Δ*∈*R*2�<sup>2</sup> is a symmetric, positive-definite matrix.

, (8)

(9)

*<sup>α</sup>*<sup>2</sup> � *<sup>β</sup>*Δ*β*<sup>0</sup> <sup>p</sup> �

Moreover, the structural matrix *Δ* is assumed to have determinant Detð Þ¼ *Δ* 1. The covariance matrix corresponding to the two-dimensional NIG-distributed

The FFT algorithm for spread option pricing along the line of [3] can be

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

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

array operations in parallel, both on multi-core CPUs as well as GPUs.

*2.2.2 Two-dimensional normal inverse Gaussian (NIG) Levy process*

*<sup>μ</sup><sup>T</sup>* <sup>þ</sup> *<sup>δ</sup><sup>T</sup>* ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

<sup>Σ</sup>*NIG* <sup>¼</sup> *δ α*<sup>2</sup> � *βΔβ*<sup>0</sup> � ��<sup>1</sup>

the covariance matrix.

of *S* ¼ ð Þ *S*1, *S*<sup>2</sup> is given by

ΦNIGð Þ¼ *u*; *T* exp i*u*<sup>0</sup>

random variable *S* is

**2.3 FFT algorithm**

described as follows

**3. Numerical results**

**96**

**3.1 Implementation outlook**

where *<sup>α</sup>*, *<sup>δ</sup>*∈*R*þ, *<sup>β</sup>*, *<sup>μ</sup>*∈*R*<sup>2</sup>


**Table 2.**

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

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.

*Spread Option Pricing on Single-Core and Parallel Computing Architectures*

of GPUs through the tools available in MATLAB.

*DOI: http://dx.doi.org/10.5772/intechopen.93430*

**Author details**

**99**

Shiam Kannan<sup>1</sup> and Mesias Alfeus2,3\*

3 University of Stellenbosch, South Africa

provided the original work is properly cited.

2 University of Wollongong, Wollongong, Australia

\*Address all correspondence to: malfeus@uow.edu.au; mesias@sun.ac.za

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

1 Ridge High School, NJ, USA

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

### **Figure 1.** *Numerical results for spread option pricing.*

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 *Φ*NIGð Þ *u*; *T* , 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

### *Spread Option Pricing on Single-Core and Parallel Computing Architectures DOI: http://dx.doi.org/10.5772/intechopen.93430*

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.
