A boundary element approach to barrier option pricing in Black–Scholes framework

Abstract

We treat the way to achieve great computational savings and accuracy in the evaluation of barrier options through boundary element method (BEM). The proposed method applies to quite general pricing models. The only requirement is the knowledge of the characteristic function for the underlying asset distribution, usually available under general asset models. This paper serves as an introductory work to illustrate the implementation of BEM using numerical Fourier inverse transform of the characteristic function and to numerically show its stability and efficiency under simple frameworks such as the Black–Scholes model.

Keywords:

• boundary element method
• Black–Scholes equation
• barrier options

AMS Subject Classifications:

• 65R20
• 65N38
• 91G60
• 91G20
• 91G80 in MSC2010

Acknowledgments

We thank the two anonymous referees for their helpful comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

1 The single layer potential $\mathcal{V}$ related to a function φ is defined as the operator (9) $\mathcal{V}\left[\phi \right](x,\tau ):={\int }_0^{\tau }{\int }_{\mathrm{\Gamma }}\phi (y,s)G(x-y,\tau -s)\mathrm{d}y\mathrm{d}s.$(9) It is related to the Green theorem, where the full representation formula for the solution of an interior problem defined in a domain Ω with boundary Γ is defined as (10) $u(x,\tau ):={\int }_{\mathrm{\Omega }}u(y,0)G(x-y,\tau )\mathrm{d}y+\mathcal{V}\left[\frac{\mathrm{\partial }u}{\mathrm{\partial }{\mathrm{n}}_y}\right](x,\tau )-\mathcal{K}\left[u\right](x,\tau )$(10) and, in general, it involves also the double layer potential $\mathcal{K}$ defined as (11) $\mathcal{K}\left[\phi \right](x,\tau ):={\int }_0^{\tau }{\int }_{\mathrm{\Gamma }}\phi (y,s)\frac{\mathrm{\partial }G(x-y,\tau -s)}{\mathrm{\partial }{\mathrm{n}}_y}\mathrm{d}y\mathrm{d}s.$(11) where ${\mathrm{n}}_y$ is the unitary outward normal at $y\in \mathrm{\Gamma }$. For a deep examination of the properties of these integral operators, look at [5].

2 We define the Fourier transform of a function $f\left(x\right)$ as $\hat{f}\left(\omega \right):=\mathfrak{F}\left[f\right]\left(\omega \right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f\left(x\right){\mathrm{e}}^{\mathrm{i}\omega x}\mathrm{d}x$ and the Fourier inverse transform as $f\left(x\right):={\mathfrak{F}}^{-1}\left[\hat{f}\right]\left(x\right)=\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\hat{f}\left(\omega \right){\mathrm{e}}^{-\mathrm{i}\omega x}\mathrm{d}\omega .$

3 Although referring to other differential problems, similar difficulties were noted by the author in [10, p. 309]:

Due to the oscillating character of the kernels in the Fourier-BEM and due to the infinite integration intervals, the computing time is not in the range of the CPU for standard BEM integrations. Hence, further work is needed to develop particular numerical integration routines to take into account the oscillating character of the kernels and the infinite bounds of the integrals.

4 When considering N nodes $x_i$ and weights $w_i$ of the Gauss–Hermite quadrature rule, the following approximation of the integral: ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f\left(x\right){\mathrm{e}}^{-b(x-a{)}^2}\mathrm{d}x\approx \sum _{i=1}^N{w}_{i}f\left(x_i\right)$ has order of accuracy equal to $2N-1$, i.e. the approximation is exact if $f\left(x\right)$ is a polynomial of order $2N-1$. Convergence of Gauss–Hermite formula has been proved for smooth functions f but also for functions with singularities. This integral has been evaluated using the algorithm in [12] for Gauss–Hermite integration over $(-\mathrm{\infty },\mathrm{\infty })$.

5 Consider the following integral: $I={\int }_a^{b}f\left(x\right)\mathrm{d}x={\int }_0^{1}f\left(\right(b-a)s+a)(b-a)\mathrm{d}s$ and then the transformation $\varphi :[0,1]\to [0,1]$ s.t. its derivatives ${\varphi }^{\left(i\right)}\left(0\right)={\varphi }^{\left(j\right)}\left(1\right)=0$ for $i=1\cdots p-1$ and $j=1\cdots q-1$: $\varphi \left(t\right)=\frac{(p+q-1)!}{(p-1)!(q-1)!}{\int }_0^t{u}^{p-1}(1-u{)}^{q-1}\mathrm{d}u.$ We apply the Gauss–Legendre quadrature rule to the integral $I={\int }_0^{1}f\left(\right(b-a\left)\varphi \right(t)+a)(b-a){\varphi }^{\prime }\left(t\right)\mathrm{d}t,$ it turns out that nodes move to the lower bound a if p increases and they move to the upper bound b if q increases.

6 The transformation is ${\int }_0^{\mathrm{\infty }}F\left(y\right)\mathrm{d}y={\int }_0^{1}F\left(\frac{\tilde{y}}{1-\tilde{y}}\right)\frac{1}{(1-\tilde{y}{)}^2}\mathrm{d}\tilde{y}$ omitting ∼ over y to simplify the notation.

7 To exploit this technique, nonlinear transformations like $\omega =L\cot \tilde{\omega }$ or $\omega =\tilde{\omega }/(1-{\tilde{\omega }}^2)$ are not suitable.

8 $\mathrm{ceil}[\cdot ]$:=function that rounds its argument to the nearest integers towards infinity.

9 A numerical method has order of convergence equal to p if the absolute error between the exact solution and the numerical approximation obtained with discretization step h is proportional to h, i.e. ${\mathrm{Err}}_h=O\left(h^p\right).$ Following this definition, when we halve the discretization step we expect that, for a method of order p, $\mathrm{Ratio}=\frac{{\mathrm{Err}}_h}{{\mathrm{Err}}_{h/2}}\sim 2^p.$

10 All the numerical simulations have been performed with a laptop computer: CPU Intel i5, 4Gb RAM. From a theoretical standpoint, we could also compare the computational costs evaluated according to relations (38(38) $C_{\mathrm{BEM}}\sim N_{\mathrm{\Delta }t}\cdot N_Q+N_{\mathrm{\Delta }t}^2+N\cdot (N_{\mathrm{\Delta }t}+1)\cdot N_Q.$(38) ), (39(39) $C_{\mathrm{BEM}}^{\mathfrak{F}}\sim N_{\mathrm{\Delta }t}\cdot N_Q\cdot N_{\mathfrak{F}}\log \left(N_{\mathfrak{F}}\right)+N_{\mathrm{\Delta }t}^2+N\cdot (N_{\mathrm{\Delta }t}+1)\cdot N_Q\cdot N_{\mathfrak{F}}\log \left(N_{\mathfrak{F}}\right).$(39) ), (48(48) $C_{\mathrm{FD}}\sim N_{\mathrm{\Delta }t}\cdot (N_{\mathrm{\Delta }x}-1).$(48) ), assuming for easiness that for all BEM and Fourier integrals the number of nodes is equal to 48. Although in this way the effective costs of BEM are abundantly increased, we will come to the same considerations about BEM efficiency even when considering numerical transform.

11 Alternatively, when considering time-dependent volatility (i.e. the diffusion coefficient in the partial differential equation), we need to better specify the integral formulation (12(12) $u(x,\tau )={\int }_{\mathrm{\Omega }}u(y,0)G(x-y,\tau )\mathrm{d}y+{\int }_0^{\tau }{\int }_{\mathrm{\Gamma }}p(y,s)G(x-y,\tau -s)\mathrm{d}y\mathrm{d}s,$(12) ). In fact, the Green theorem yields the following boundary integral representation of the solution u: (51) $u(x,\tau )={\int }_{\mathrm{\Omega }}u(y,0)G(x-y,\tau )\mathrm{d}y+{\int }_0^{\tau }{\int }_{\mathrm{\Gamma }}p(y,s)G(x-y,\tau -s)\mathrm{d}y\mathrm{d}s,$(51) where the unknown term in this case is in the form $p(y,s)={\sigma }^2\left(s\right)/2\mathrm{\partial }u(y,s)/\mathrm{\partial }y$, i.e. it contains the time-varying volatility as well, and the fundamental solution (available from [24]) takes the form (52) $G(x-y,\tau -s)=\frac{1}{\sqrt{2\pi {\int }_0^{\tau -s}{\sigma }^2\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }}}\exp \left\{-\frac{[x-y+{\int }_0^{\tau -s}(r-d-{\sigma }^2\left(t^{\prime }\right)/2)\mathrm{d}{t}^{\prime }{]}^2}{2{\int }_0^{\tau -s}{\sigma }^2\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }}\right\},\tau >s.$(52) Something similar happens in Heston model too, where the volatility is a stochastic process, and it will be analysed in detail in a further paper [16].

Additional information

Funding

This work was supported by ‘National Group of Computing Science (GNCS-INdAM)’.