Parallel Option Pricing on GPU: Barrier Options and Realized Variance Options

1. Introduction

We propose parallel algorithms (named Application 1 and Application 2) implemented in CUDA subroutines ready to run on Graphics Processing Units (GPUs) to price two kinds of financial derivatives, that is: continuous barrier options and realized variance options. The really satisfactory parallel performances of these subroutines are mainly due to two reasons: the mathematical treatment of the problems considered that leads to pricing formulae very well suited for GPU parallel computing and the flexibility of CUDA language that makes possible to implement parallel algorithms to exploit these formulae. The analysis of the parallel performance on a GPU of the CUDA implementation of these algorithms shows that the use of GPU computing can be of great benefit in financial applications. The use of parallel computing has gained importance in finance in the last decades due to the fact that the speed of the decision making process is a key success factor in the financial markets operations. The decision making process in the financial institutions (banks, insurance companies, hedge funds, ...) often involves the evaluation of a huge number of contracts of the same type varying the values assigned to some of the parameters that define the contracts and/or the models used to price them (such as, for example, strike price, maturity time, volatility, drift). This is a feature of the decision making process in the financial institutions very well suited for parallel and/or distributed computing. Moreover often the evaluation of an individual contract is done using Monte Carlo methods or reduces to the evaluation of an integral via numerical quadrature. These kinds of computation are also very well suited for parallel and /or distributed computing.

The work presented here is the CUDA implementation of two algorithms that exploit pricing formulae derived in [3] (Application 1, barrier options) and in [2] (Application 2, realized variance options). A more detailed version of the work summarized here is presented in [1].

A C++ implementation of numerical methods that evaluate the pricing formulae derived in [3], [2] is available as two subroutines of the QuantLib software library (http://quantlib.org/index.shtml) respectively in the websites:

http://quantlib.org/reference/class_quant_lib_1_1_perturbative_barrier_option_engine.html. .

and

http://quantlib.org/reference/class_quant_lib_1_1_integral_heston_variance_option_engine.html

where respectively the problem of pricing barrier options and realized variance options is considered. These C++ subroutines do not make use of parallel computing. A grid enabled version of these C++ subroutines and of some other QuantLib subroutines [5] based on the use of a middleware software product of Platform Computing Inc., Toronto, (http://www.platform.com/) known as Symphony DE application can be downloaded from the website:

http://www.ceri.uniroma1.it/ceri/zirilli/w5/.

A more detailed version of the work summarized here is presented in [1]. A general reference to the work of some of the authors and of their coauthors in mathematical finance is the website: http://www.econ.univpm.it/recchioni/finance/ .

Experimental setting

In the numerical experiments discussed below the GPU used is the NVIDIA "GeForce GTX 285" having 30 multiprocessors for a total of 240 cores (each core is a GPU processing unit) and the CPU used is an Intel core 2 duo E8500 3.18 GHz. The GPU speed up factors presented are computed as the ratio T_CPU/T_GPU,p where p is the number of processing units (i.e. cores) of the GPU used, T_GPU,p is the (wall clock) execution time of the parallel algorithm on p processing units of the GPU and T_CPU is the (wall clock) execution time of the corresponding sequential algorithm on the CPU. The subroutines executed on the CPU are the QuantLib C++ subroutines mentioned above. The results presented are obtained using p=240, that is in the numerical experiments all the 240 cores of the GPU are used in the computation. Here we limit our attention to the parallel performance of the codes on the GPU, a more complete analysis that includes the study of the accuracy of the results obtained can be found in [1], [2], [3].

Application 1: Barrier options

A continuously monitored up and out single barrier option is an option that is terminated when the price of the underlying asset goes above a fixed trigger (i.e. the barrier) within the lifetime of the option contract [6].

Application 1 prices a barrier option in the Black and Scholes framework with time dependent parameters evaluating numerically a truncated series expansion of the option price. The first three terms of this series expansion have been written as expressions containing integrals of explicitly known integrands that involve some elementary and non elementary transcendental functions (see [1], [3], [5] and this website http://www.econ.univpm.it/recchioni/finance/w3 for further details). Application 1 evaluates numerically these three terms in a way well suited for GPU computing [1].

The CUDA implementation of Application 1 is rather complex. In fact in Application 1 the option value is approximated using a truncation of a perturbative series expansion, that we denote with P₀+P₁+P₂ where P₀, P₁, P₂ are respectively the zeroth, first and second order terms of the series expansion considered (see [1], [3] for details). Let N₁ and N₂ be positive integers. The computation of the zero-th order term P₀ involves only the evaluation of elementary functions (such as, for example, univariate normal distributions). The computation of the first order term P₁ involves the evaluation of a one dimensional integral of bivariate normal distributions. This integral is approximated using a numerical quadrature formula with N₁ quadrature nodes. The computation of the second order term P₂ involves the evaluation of multidimensional integrals of bivariate and trivariate normal distributions that are approximated using a numerical quadrature formula with N₁·N₂ quadrature nodes. The evaluation of the term P₂ is the most demanding part of the computation and for this reason has been implemented in CUDA and is executed on the GPU using parallel computing. The sum P₀+P₁ is always computed on the CPU. In fact when N₁, N₂ are large enough to guarantee accurate approximations of the quantities computed the numerical experience shows that the (wall clock) execution time required to compute the approximation of P₀+P₁ on the CPU is always negligible with respect to the (wall clock) execution time required to compute the approximation of P₂ on the GPU. Hence in the CUDA implementation of Application 1 the use of the GPU is limited to the computation of the second order term P₂ . In fact the GPU implementation of the computation of P₀+P₁ will not change significantly the time required to compute the approximation of the barrier option price given by P₀+P₁+P₂. This way of doing is called "hybrid CUDA implementation" of the computation P₀+P₁+P₂.

A test problem consisting in pricing a barrier option considered in [1], [3] is studied. The results obtained are shown in Table 1 where the (wall clock) execution times are shown. Table 1 shows that the "hybrid CUDA implementation" (t^*_GPU,240) of the computation of P₀+P₁+P₂ (i.e. P₀+P₁ computed on the CPU and P₂ computed on the GPU) outperforms the CPU implementation (t^*_CPU) of the same computation in term of (wall clock) execution time. Table 1 shows a speed up factor belonging to the interval [15,20] that slightly increases when the number N₁·N₂ increases. Moreover Table 1 shows that using parallel computing on the GPU we can compute the value of the barrier option considered in an affordable time (a few minutes) even when the number of quadrature nodes employed is large while the same computation on the CPU will require a few hours. Note that when N₁, N₂ are chosen as in Table 1 and N₁·N₂ is of order 10⁶ in the example considered in [1], [3] at least six, seven correct significant digits are guaranteed in the numerical quadratures.

Let n₁, n₂ be positive integers, we must mention the fact that when one tries to compute at once on the GPU n₁·n₂ function values of the integrand function of the two dimensional integrals that define P₂ and n₁, n₂ are large enough the GPU encounters a so called "watchdog timer" error. This error consists in the fact that the operating system stops the GPU functioning when the GPU remains out of the control of the operating system (lack of activity for the operating system) for a too long time. To avoid the GPU crash induced by the "watchdog timer" error we have divided the computation of the sum of the addenda contained in the numerical approximation of the integrals that define P₂ in blocks of size N_b, that is we have split the computation partitioning the addenda using a block matrix of dimension N₁·N₂/N_b and each block M_k,m, k = 1,2,¼,N₁, m = 1,2,¼N₂/N_b, is made of 1·N_b addenda. These blocks are distributed between the GPU cores to be processed. Note that, for simplicity, we are assuming that N₂/N_b is an integer. Roughly speaking, for k = 1,2,¼,N₁, m = 1,2,¼N₂/N_b the computation of the sum of the addenda contained in each block M_k,m is assigned to a GPU core and each core to perform its duty activates and executes a certain number of threads (in our numerical experience up to 512 threads). In our implementation N_b is the number of threads activated by each core. The optimal number of threads that each core must activate and execute to carry out the computation efficiently depends on the total numbers of blocks N₁·N₂/N_b. In fact the number of threads activated and executed by each core must be large enough to take the GPU busy but not so large to clog (to block) it due to the "watchdog timer" error. Figure 2 shows the (wall clock) execution time as a function of N₁·N₂ when N₁, N₂ are chosen as in Table 1 for several choices of N_b. We note that the choices N_b =64 and N_b =256 perform well.

It is known that the most effective GPU computations make use of the so called "inline" compiler directives. However Application 1 cannot be implemented using the "inline" compiler directives on the GPU that we are using in our experiments since in the case of Application 1 the execution of the subroutine calls in the GPU core from where they originate that is required by the "inline" compiler directives saturates the RAM memory available to the individual GPU core. Therefore the CUDA implementation of Application 1 has been realized using the "noinline" compiler directives. This choice limits the parallel performance of the code. However, as shown in Table 1, the performance of the hybrid CUDA implementation of the computation of P₀+P₁+P₂ with respect to the CPU implementation of the same computation in terms of (wall clock) execution time is satisfactory. In fact the speed up factors obtained are in the range fifteen, twenty. In Table 1 t^*_GPU,240 = t_P0+P1(CPU)+t_P2(GPU) is the time needed to compute P₀+P₁+P₂ using the hybrid CUDA implementation and t^*_CPU is the time needed to compute P₀+P₁+P₂ using the CPU.

The CUDA implementation of Application 1 can be downloaded here:

Software Download click here to download the CUDA code to price continuous up and out single barrier put options.

The CUDA implementation of Application 2 can be downloaded here:

Software Download click here to download the CUDA code to price realized variance options.