1 Introduction
An implied volatility is the value of the volatility parameter for which an observed market price fits the Black–Scholes option pricing formula [Reference Black and Scholes3]. Maturity and strike-dependent computed implied volatilities, contrary to the assumption of constant volatility in the Black–Scholes model, resulted in the search for models capable of fitting market-observed volatilities. Dupire [Reference Dupire11] introduced a local volatility model, where, for a given maturity, the asset price process matches the strike-dependent volatility profile, usually referred to as a smile dynamics. One drawback of the Dupire model is that the smile shifts are opposite to the observed market behaviour of asset prices. This observation led to the introduction of the stochastic alpha–beta–rho (SABR) model [Reference Hagan, Kumar, Lesniewski and Woodward16]. The widespread use of SABR in financial markets stems from the facts that, in addition to being a local stochastic volatility model and having the capability of fitting various market volatility structures, this option pricing model also eliminates the problem of asset prices and smiles moving in opposite directions as in Dupire’s model.
The SABR model incorporates a constant elasticity of variance (CEV)-type diffusion process for the forward price whose volatility follows a Black–Scholes-type diffusion with zero drift. The probability of the forward price process hitting zero is positive, but it exponentially decays to zero with shrinking time horizon with faster convergence rates for large initial asset prices or elasticity parameter of the CEV dynamics in the model or small values of the forward price’s initial volatility or the volatility of its volatility process [Reference Chen and Yang6].
Computations of prices which are arbitrage-free require the imposition of an absorbing boundary condition at zero. Ignoring this boundary condition in a numerical method runs the risk of producing negative probability densities for low strikes [Reference Hagan, Kumar, Lesniewski and Woodward17]. The quantification of not having an absorbing boundary condition carried out by Chen and Yang [Reference Chen and Yang6] shows that the approximation error in prices is negligible only for small maturities.
Specifying an absorbing boundary condition at zero and using an expansion around a one-dimensional Bessel process, the closed-form approximation for European options [Reference Yang, Chen, Liu and Wan33] has been shown to be accurate for small maturity problems [Reference Thakoor, Tangman and Bhuruth31]. This implies that the search for a numerical method capable of accurately pricing options with large maturities needs to incorporate the absorbing boundary condition in order to ensure that computed prices are arbitrage-free.
Observing that the boundary layer next to zero forward price has a significant influence on pricing, an effective one-dimensional equation for the probability density with an absorbing boundary condition was derived [Reference Hagan, Kumar, Lesniewski and Woodward17]. By numerical solution of this equation using a Crank–Nicolson time-stepping scheme [Reference Crank and Nicolson8], call and put values can be obtained by numerical integration and this procedure leads to arbitrage-free prices. Other methods proposed for the numerical determination of the probability density include the trapezoidal rule with the second-order backward difference formula [Reference Le Floc’h and Kennedy21] and an exponential time integration scheme [Reference Thakoor29]. This one-dimensional approach is not exact, and the solution of the two-dimensional SABR partial differential equation (PDE) provides a more accurate means for price determination.
This work addresses the issue of accurate price computations by proposing a localized radial basis functions (RBFs) method for solving the two-dimensional SABR pricing equation. Applications of RBFs in a local mode have gained much importance in the numerical solutions of PDEs in the field of computational fluid dynamics. Their applications to PDEs in finance have been much less considered, and the present work develops a new local RBF method for the solution of the arbitrage-free SABR PDE. Weighting coefficients for generalized multiquadrics (GMQ) presented here are new. The numerical discretization in space is then carried out using multiquadrics and inverse multiquadrics, which are special cases of GMQ. One clear advantage of using RBFs is that due to their meshless nature, unstructured grids or randomly scattered grid points can be employed to achieve accurate prices. An extensive set of numerical examples using structured and randomly scattered points are considered to demonstrate that accurate no-arbitrage prices are computed.
The presentation of the work carried out is as follows. The SABR PDE is described in Section 2 and a brief account of some recent studies related to the SABR model is mentioned. Local generalized multiquadrics approximations for the first and second derivatives are derived in Section 3. The numerical method for the valuation of options under SABR is developed in Section 4 and the time-stepping scheme for the solution of the system of differential equations is described. Numerical examples are considered in Section 5 and a summary of the work, conclusions reached and possible future work are given in Section 6.
2 No-arbitrage SABR PDE
Consider a financial market with a risky asset with forward price process
$F_t$
and volatility
$\alpha _t$
over a time horizon
. Let
$W_t$
and
$\widehat {W}_t$
be Brownian motions on the filtered probability space
$(\Omega ,\mathcal {F},\mathcal {F}_t,\mathbb {Q})$
, where
$\Omega $
is the sample space,
$\mathcal {F}$
is a
$\sigma $
-algebra and
$\mathbb {Q}$
is a martingale measure. The filtration
$\mathcal {F}_t$
is defined below. Let
$\rho \in (-1,1)$
be the correlation parameter between the Brownian motion
$W_t$
driving the constant elasticity of variance dynamics with elasticity parameter
$\beta $
for
$F_t$
and the Brownian motion
$\widehat {W}_t$
driving the geometric Brownian motion process with zero drift for
$\alpha _t$
[Reference Hagan, Kumar, Lesniewski and Woodward16]. Let
$$ \begin{align*} \overline{W}_t=(1-\rho^2)^{-1/2}(W_t-\rho \widehat{W}_t).\end{align*} $$
Consider the filtration
$\mathcal {F}_t=\mathcal {F}_t^{\overline {W}}\otimes \mathcal {F}_t^{\widehat {W}}$
, where
$\mathcal {F}_t^{\overline {W}}$
and
$\mathcal {F}_t^{\widehat {W}}$
are the filtrations generated by the independent Brownian motions
$\overline {W}_t$
and
$\widehat {W}_t$
, respectively. Denoting by
$\nu $
the volatility of the volatility process
$\alpha _t$
, the SABR dynamics of
$(F_t,\alpha _t)$
can be specified by the system of two stochastic differential equations given by
$$ \begin{align} \begin{aligned} dF_t & =\alpha_t\bigg(\sqrt{1-\rho^2}d\overline{W}_t+\rho d\widehat{W}_t\bigg)F_t^{\beta},\quad F_0=f, \\ d\alpha_t & =\nu\alpha_t d\widehat{W}_t,\quad \alpha_0=\hat{\alpha}. \end{aligned} \end{align} $$
Let
$\phi (F_t,\alpha _t)$
be the payoff at time t of a contingent claim on the underlying asset. The infinitesimal generator
$\mathcal {L}$
of the payoff
$\phi (F,\alpha )$
defined as
is given by Cui et al. [Reference Cui, Kirkby and Nguyen9] as
$$ \begin{align*} \mathcal{L}\phi=\frac{1}{2}\alpha^2F^{2\beta}\frac{\partial^2\phi}{\partial F^2}+\rho\nu\alpha^2\frac{\partial^2\phi}{\partial F\partial \alpha}+\frac{1}{2}\nu^2\alpha^2 \frac{\partial^2\phi}{\partial \alpha^2}.\end{align*} $$
The payoff for a call option with strike K is
$\phi (F,\alpha )=(F-K)^{+}=\max (F-K,0)$
, and for a put option
$\phi (F,\,\alpha )=(K-F)^+$
. Let
$\eta _t$
be the first time that the forward price process hits zero, and let the time-t price of a contingent claim on the underlying asset be given by
$$ \begin{align*} V(F,\alpha,t)=\mathbb{E}[\phi (F_t,\alpha_t)\mid F_t=F,\alpha_t=\alpha].\end{align*} $$
Let the corresponding call price be
$V_c (F,\alpha ,t)$
and the put price be
$V_p(F,\alpha ,\,t)$
. Then, denoting the indicator function for a set S by
$\mathbf {1}_S$
, the call and put prices are given by [Reference Yang, Chen, Liu and Wan33]
$$ \begin{align*} V_c(F,\alpha,t)=\mathbb{E} [(F_T-K)^+\mathbf{1}_{\{\eta_t>T\}}\mid F_t=f,\,\alpha_t=\alpha] \end{align*} $$
and
Pricing a European call option with strike K is the same as pricing a down-and-out call option with a knockout boundary at zero. The call price
$V_c(F,\,\alpha ,\,t)$
is the unique solution [Reference Yang, Chen, Liu and Wan33] of the PDE given by
$$ \begin{align} \frac{\partial V(F,\alpha,t)}{\partial t} +\mathcal{L}V(F,\alpha,t)=0,\quad (F,\alpha,t)\in \mathbb{R}^+\times \mathbb{R}^+\times [0,T] \end{align} $$
with terminal condition
$V(F,\alpha ,T) =(F-K)^+$
and boundary condition
$V(0,\alpha ,t)=0$
.
Pricing a European put option with strike K with an absorbing boundary condition at zero is the same as pricing a down-and-out put option with a rebate payment of K if the knock-out boundary at zero is reached. An analytical approximation for the probability that the forward price process hits zero was derived by Yang and Wan [Reference Yang and Wan35]. Let
$$ \begin{align*} P_0(F,\alpha,t)&=\mathbb{E}[\mathbf{1}_{\{\eta_t>T\}}\mid F_t=f,\,\alpha_t=\alpha]\\ &=\mathbb{P}(\eta_t>T\mid F_t=f,\,\alpha_t=\alpha). \end{align*} $$
Then
$P_0(F,\alpha ,t)$
is the unique solution of the PDE
$$ \begin{align*} \frac{\partial P_0(F,\alpha,t)}{\partial t} + \mathcal{L}P_0(F,\alpha,t)=0,\quad (F,\alpha,t)\in \mathbb{R}^+\times \mathbb{R}^+\times [0,T] \end{align*} $$
with condition
$P_0(0,\alpha ,t) =0$
and boundary condition
$P_0(F,\alpha ,T)=1$
.
It then follows that the European put price is given by
$$ \begin{align*} V_p(F,\alpha,t)=P_1(F,\alpha,t)+K(1-P_0(F,\alpha,t)), \end{align*} $$
where
$P_1(F,\alpha ,t)$
is the unique solution of the PDE
$$ \begin{align*} \frac{\partial P_1(F,\alpha,t)}{\partial t} +\mathcal{L}P_1(F,\alpha,t)=0,\quad (F,\alpha,t)\in \mathbb{R}^+\times \mathbb{R}^+\times [0,T] \end{align*} $$
with terminal condition
$P_1(F,\alpha ,T) =(K-F)^+$
and boundary condition
${P_1 (0,\alpha ,t)=0}$
.
3 General multiquadrics finite-difference approximations
A numerical scheme for the solution of (2.2) requires approximations of the space derivatives
$\partial ^2 V/\partial F^2$
,
$\partial ^2 V/\partial F\partial \alpha $
and
$\partial ^2V/\partial \alpha ^2$
. The method proposed here is based on obtaining three-point finite-difference (FD) approximations of the derivatives using general multiquadrics (GMQ) radial basis functions.
Let
$u(x)$
be a one-dimensional function whose derivatives
$u'(x_0)$
and
$u''(x_0)$
are to be approximated using values of u at the nodes
$x_0-h$
,
$x_0$
and
$x_0+h$
. Thus, denoting by
$u^{(k)}(x_0)$
the kth derivative of u at the node
$x_0$
, we seek approximations of the form
$$ \begin{align*} u^{(k)}(x)\approx a_{-1}^{(k)}u(x_0-h)+a_0^{(k)}u(x_0)+a_1^{(k)}u(x_0+h).\end{align*} $$
The weighting coefficients
are to be found using GMQs given by
such that
Note that
$$ \begin{align*} \varphi_i'(x)=\frac{m(x-x_0-ih)\varphi_i(x)}{c^2+(x-x_0-ih)^2},\quad \varphi_i''(x)=\frac{m(c^2+(m-1)(x-x_0-ih)^2)\varphi_i(x)}{(c^2+(x-x_0-ih)^2)^2}.\end{align*} $$
Letting
the weighting coefficients for the approximation of the first derivative then satisfy the equation
$$ \begin{align*} \begin{pmatrix} w_0 & w_1 & w_2 \\ w_1 & w_0 & w_1 \\ w_2 & w_1 & w_0 \\ \end{pmatrix}\begin{pmatrix} a_{-1}^{(1)} \\[.1cm] a_0^{(1)} \\[.1cm] a_1^{(1)} \\ \end{pmatrix}=q_0\begin{pmatrix} 1 \\ 0 \\ -1 \\ \end{pmatrix}.\end{align*} $$
Solving the above linear system gives the weighting coefficients as
$$ \begin{align*} a_{-1}^{(1)}=\frac{q_0}{c^m-(c^2+4 h^2)^{m/2}},\quad a_0=0,\quad a_1^{(1)}=-\frac{q_0}{c^m-(c^2+4 h^2)^{m/2}}.\end{align*} $$
For finding the weighting coefficients for the second derivative, let
$$ \begin{align*} q_1=\frac{m(c^2+(m-1)h^2)(c^2+h^2)^{m/2}}{(c^2+h^2)^2}.\end{align*} $$
The weighting coefficients are solutions of the linear system
$$ \begin{align*} \begin{pmatrix} w_0 & w_1 & w_2 \\ w_1 & w_0 & w_1 \\ w_2 & w_1 & w_0 \\ \end{pmatrix}\begin{pmatrix} a_{-1}^{(2)} \\[.1cm] a_0^{(2)} \\[.1cm] a_1^{(2)} \\ \end{pmatrix}=\begin{pmatrix} q_1 \\[.1cm] mc^{m-2} \\ q_1 \\ \end{pmatrix}.\end{align*} $$
Solving the above linear system gives
$$ \begin{align*} a_{-1}^{(2)}=a_1^{(2)}=\frac{-mc^{m-2}h^2(c^2+h^2)^{{m/2-2}}(h^2-(m-3)c^2)}{ c^m (c^m+(c^2+4h^2)^{m/2})-2(c^2+h^2)^m}\end{align*} $$
and
$$ \begin{align*} a_0^{(2)}=\frac{m}{c^2}\bigg(1+\frac{2h^2(c^2+h^2)^{m-2}(h^2-(m-3)c^2)}{ c^m(c^m+(c^2+4h^2)^{m/2})-2(c^2+h^2)^m}\bigg).\end{align*} $$
In the limit when the shape parameter c is much larger than the grid size h, that is,
$c\gg h$
, it can be shown that
$$ \begin{align*} a_1^{(1)}=\frac{1}{2h}\bigg(1+\bigg(1-\frac{m}{2}\bigg)\frac{h^2}{c^2} -\frac{m^2-30m+56}{24}\frac{h^4}{c^4}\bigg),\end{align*} $$
$$ \begin{align*} a_{1}^{(2)}=\frac{1}{h^2}\bigg(1-\frac{(m-2)(m-5)}{2(m-3)}\frac{h^2}{c^2}+\frac{(m-2)(m^3 +8m^2-111m+282)}{24(m-3)^2}\frac{h^4}{c^4}\bigg)\end{align*} $$
and
$$ \begin{align*} a_0^{(2)}=-\frac{2}{h^2}\bigg(1-\frac{(m-5)(m-2)}{2(m-3)}\frac{h^2}{c^2}-\frac{(m-2)(m^3-19m^2+87m-141)}{12 (m-3)^2}\frac{h^4}{c^4}\bigg).\end{align*} $$
The formulas presented in this work are new and these recover the special cases of RBF-FDs for multiquadrics [Reference Bayona, Moscoso, Carretero and Kindelan1] and inverse multiquadrics [Reference Soleymani, Barfeie and Haghani26].
4 Development of RBF-FD method for SABR PDE
This section develops the numerical method for the solution of the SABR PDE (2.2) using the RBF-FD approximations derived in the previous section. For the pricing of options under the Heston model, there exist a few alternating direction implicit (ADI) schemes [Reference Düring and Miles13, Reference in’t Hout and Foulon18, Reference Zhu and Chen36]. Although it is possible to develop an ADI-RBF-FD method in a similar way, a different approach is developed in this work.
Assume that
$\beta>0$
and
$F\geq 0$
. Using the substitution
$\tau =T-t$
, the pricing equation (2.2) becomes
$$ \begin{align} \frac{\partial V}{\partial \tau}=\frac{1}{2}\alpha^2F^{2\beta}\frac{\partial^2 V}{\partial F^2}+\rho\nu\alpha^2 F^{\beta}\frac{\partial^2 V}{\partial F\partial \alpha}+\frac{1}{2}\nu^2\alpha^2\frac{\partial^2 V}{\partial \alpha^2},\quad (F,\alpha,\tau)\in \mathbb{R}^+\times \mathbb{R}^+\times [0,T] \end{align} $$
with initial condition
$V(F,\alpha ,0)=\phi (F,\alpha )$
and absorbing boundary condition
$V(0,\alpha ,\tau )=0$
.
The pricing equation (4.1) is defined on an unbounded domain. For development of the RBF-FD method, the problem is localized to a finite domain
$[0,F_{\max }]\times [0,\alpha _{\max }]\times [0,T]$
. Localizing the problem requires that additional boundary conditions be imposed. For a call option, these conditions are given by Kienitz et al. [Reference Kienitz, McWalter and Sheppard20] as
$$ \begin{align*} \frac{\partial V}{\partial F}(F_{\max},\alpha,\tau)=1,\quad \frac{\partial V}{\partial \tau}(F,0,\tau)=0\end{align*} $$
and
$$ \begin{align} \begin{aligned} \frac{\partial V}{\partial \tau}(F,\alpha_{\max},\tau)&=\frac{1}{2}\alpha_{\max}^2 F^{2\beta}\frac{\partial^2 V}{\partial F^2}(F,\alpha_{\max},\tau). \end{aligned} \end{align} $$
For a put option, the additional boundary conditions are
$$ \begin{align*} \frac{\partial V}{\partial F} (F_{\max},\alpha,\tau)=0,\quad \frac{\partial V}{\partial \tau} (F,0,\tau)=0\end{align*} $$
and
$$ \begin{align*} \frac{\partial V}{\partial \tau} (F,\alpha_{\max},\tau)=\frac{1}{2}\alpha_{\max}^2 F^{2\beta}\frac{\partial^2 V}{\partial F^2} (F,\alpha_{\max},\tau).\end{align*} $$
Consider uniformly spaced grid nodes in the F- and
$\alpha $
-directions with spacings
$h_F=F_{\max }/M$
and
$h_{\alpha }=\alpha _{\max }/J$
, respectively. Denote grid nodes in the
$(F,\alpha )$
-space by
$(F_i,\alpha _j)$
, where
and
$\alpha _j=jh_{\alpha }$
for
.
Let
$V_{ij}(\tau )$
denote the option price at the node
$(F_i,\alpha _j,\tau )$
and consider the following RBF-FD approximations:
$$ \begin{align*} D_{FF}V_{ij}(\tau)\approx \frac{\partial^2 V}{\partial F^2}(F_i,\alpha_j,\tau),\quad D_{\alpha \alpha}V_{ij}(\tau)\approx \frac{\partial^2 V}{\partial \alpha^2}(F_i,\alpha_j,\tau)\end{align*} $$
given by
$$ \begin{align*} D_{FF}V_{ij}(\tau)&=a_{-1,F}^{(2)}V_{i-1,j}(\tau)+a_{0,F}^{(2)} V_{ij}(\tau)+a_{1,F}^{(2)}V_{i+1,j}(\tau), \\ D_{\alpha\alpha}V_{ij}(\tau)&=a_{-1,\alpha}^{(2)}V_{i,j-1}(\tau)+a_{0,\alpha}^{(2)} V_{ij}(\tau)+a_{1,\alpha}^{(2)}V_{i,j+1}(\tau). \end{align*} $$
The algorithm development is presented for multiquadrics (MQ), which corresponds to the case
$m=1$
in (3.1). For
$c \gg h_F$
, retaining only terms in
$h_F^2/c^2$
, the weighting coefficients
$a_{i,F}^{(2)}$
for
are given by
$$ \begin{align*} a_{i,F}^{(2)}=\frac{1}{h_F^2}\bigg(1+\frac{h_F^2}{c^2}\bigg),\quad i=-1,\,1,\quad a_{0,F}^{(2)} =-\frac{2}{h_F^2}\bigg(1+\frac{h_F^2}{c^2}\bigg).\end{align*} $$
The coefficients
$a_{i,\alpha }^{(2)}$
,
, are given by similar expressions. Since the SABR PDE does not contain any convection term, no upwinding discretizations need to be employed as for the Heston PDE where the presence of convection terms due to the drift terms in the asset and volatility dynamics may lead to oscillations in the numerical solutions when convection is dominant.
The RBF-FD approximation
$D_{F\alpha }V_{ij}(\tau )$
for the cross-derivative term
$({\partial ^2 V}/{\partial F \partial \alpha }) (F_i,\alpha _j,\tau )$
is given as
$$ \begin{align*} D_{F\alpha} V_{ij}(\tau)=a_{12}( V_{i+1,j+1}(\tau)- V_{i-1,j+1}(\tau)- V_{i+1,j-1}(\tau)+V_{i-1,j-1}(\tau)),\end{align*} $$
where
$$ \begin{align*} a_{12}=\frac{1}{4h_Fh_{\alpha}}\bigg(1+\frac{3h_F^2}{2c^2} +\frac{3h_{\alpha}^2}{2c^2}-\frac{9h_F^2h_{\alpha}^2}{4c^4}\bigg).\end{align*} $$
This discretization of the cross-derivative term leads to a nine-point computational stencil at an interior point of the computational domain. Other possibilities leading to seven-point computational stencils are given by Thakoor et al. [Reference Thakoor, Tangman and Bhuruth30].
Let
$\hat {a}_{ij}=\alpha _j^2F_i^{2\beta }/2$
,
$\hat {b}_{ij}=\rho \nu \alpha _j^2F_i^{\beta }$
and
$\hat {c}_j=\nu ^2\alpha _j^2/2$
. Then, at the interior grid points
$(F_i,\alpha _j)$
of the
$(F,\alpha )$
-computational domain, the discretization of (4.1) is given by
Consider the vector
$\mathbf {U}(\tau )\in \mathbf {R}^{(M-1)(J-1)}$
of call prices at time level
$\tau $
given by
where
For
and
, denote
$$ \begin{align*} \hat{d}_{i-1,j}&=\hat{a}_{ij}a_{-1,F}^{(2)}, \qquad\hat{d}_{ij}=\hat{a}_{ij}a_{0,F}^{(2)}+\hat{c}_ja_{0,\alpha}^{(2)},\quad \hat{d}_{i+1,j}=\hat{a}_{ij}a_{1,F}^{(2)},\\ \hat{d}_{i-1,j-1}&=\hat{b}_{ij}a_{12}, \quad\ \ \,\,\hat{d}_{i,j-1}=\hat{c}_ja_{-1,\alpha}^{(2)}, \qquad\quad\,\, \hat{d}_{i+1,j-1}=-\hat{b}_{ij}a_{12},\\ \hat{d}_{i-1,j+1}&=-\hat{b}_{ij}a_{12}, \quad\,\hat{d}_{i,j+1}=a_{1,\alpha}^{(2)}\hat{c}_j, \qquad\qquad\, \hat{d}_{i+1,j+1}=\hat{b}_{ij}a_{12}. \end{align*} $$
For
, consider the tridiagonal matrices
$A_{1j}$
,
$A_{2j}$
and
$A_{3j}$
in
$\mathbb {R}^{(M-1)\times (M-1)}$
given by
$$ \begin{align*} A_{1j}=\begin{pmatrix} \hat{d}_{1,j-1} & \hat{d}_{2,j-1} & & & \\ \hat{d}_{1,j-1} & \hat{d}_{2,j-1} & \hat{d}_{3,j-1}& & \\ & \ddots & \ddots & \ddots & \\ & & \hat{d}_{M-3,j-1} & \hat{d}_{M-2,j-1} & \hat{d}_{M-1,j-1} \\ & & & \hat{d}_{M-2,j-1} & \hat{d}_{M-1,j-1} \\ \end{pmatrix},\end{align*} $$
$$ \begin{align*} A_{2j}=\begin{pmatrix} \hat{d}_{1j} & \hat{d}_{2j} & & & \\ \hat{d}_{1j} & \hat{d}_{2,j} & \hat{d}_{3,j}& & \\ & \ddots & \ddots & \ddots & \\ & & \hat{d}_{M-3,j} & \hat{d}_{M-2,j} & \hat{d}_{M-1,j} \\ & & & \hat{d}_{M-2,j} & \hat{d}_{M-1,j} \\ \end{pmatrix}\end{align*} $$
and
$$ \begin{align*} A_{3j}=\begin{pmatrix} \hat{d}_{1,j+1} & \hat{d}_{2,j+1} & & & \\ \hat{d}_{1,j+1} & \hat{d}_{2,j+1} & \hat{d}_{3,j+1}& & \\ & \ddots & \ddots & \ddots & \\ & & \hat{d}_{M-3,j+1} & \hat{d}_{M-2,j+1} & \hat{d}_{M-1,j+1} \\ & & & \hat{d}_{M-2,j+1} & \hat{d}_{M-1,j+1} \\ \end{pmatrix}.\end{align*} $$
Let the block-tridiagonal matrix
$\textbf {A}\in \mathbb {R}^{(M-1)(J-1)\times (M-1)(J-1)}$
be given by
$$ \begin{align*} \textbf{A}=\begin{pmatrix} A_{21} & A_{31} & & & \\ A_{12} & A_{22} & A_{32}& & \\ & \ddots & \ddots & \ddots & \\ & & A_{1,J-2} & A_{2,J-2} & A_{3,J-2} \\ & & & A_{1,J-1} & A_{2,J-1} \\ \end{pmatrix}.\end{align*} $$
Not accounting for the Neumann boundary conditions (4.2) (we show later in this work how these are incorporated to obtain the full set of ordinary differential equations), we obtain the initial value problem given by
with initial condition
, where
$\mathbf {u}_i=(F_i-K)^+\mathbf {\bar {\unicode{x0142}}}$
and
$\mathbf {\bar {\unicode{x0142}}}\in \mathbb {R}^{M-1}$
is a vector of ones. The vector
consists of option values not in the vector
$\mathbf {U}(\tau )$
. The vector
$ \tilde {\mathbf {b}}_j$
is given by
where its components are
For
, the vectors
$\tilde {\mathbf {b}}_j(\tau )\in \mathbb {R}^{M-1}$
have entries
and the entries of the vector
$\tilde {\mathbf {b}}_{J-1}\in \mathbb {R}^{M-1}$
are given by
and
$$ \begin{align*} \tilde{b}_{M-1,J-1}(\tau)&=\hat{d}_{M-2,J}V_{M-2,J}(\tau)+\hat{d}_{M-1,J}V_{M-1,J}(\tau)+\hat{d}_{M,J-1}V_{M,J-2}(\tau) \\ &\quad +\hat{d}_{M,J-1}V_{M,J-1}(\tau)+\hat{d}_{M,J}V_{M,J}(\tau). \end{align*} $$
The following shows how to obtain the full linear system of ordinary differential equations for time stepping. First, since the boundary condition
${V(0,\,\alpha ,\tau )=0}$
, we have
$V_{0,j}(\tau )=0$
for
. The other boundary conditions which are of Neumann type are accounted for in the following way. The condition
${(\partial V/\partial \tau )(F,0,\tau )=0}$
means that
$V_{i,0}(\tau )=0$
for
. At the boundary
$F=F_{\max }$
, the boundary condition
$(\partial V/\partial F) (F_{\max },\alpha ,\tau )=1$
is discretized by a left-sided derivative
, where
$$ \begin{align*} D^{-}_FV_{M,j}(\tau)=\frac{1}{2h_F}\bigg(1-\frac{h_F^2}{c^2}\bigg)[V_{M-2,j}(\tau)-4V_{M-1,j} (\tau)+3V_{M,j}(\tau) ].\end{align*} $$
Then, on the boundary
$F=F_{\max }$
,
$$ \begin{align} V_{M,j}(\tau)=\frac{1}{3}\bigg(\frac{2h_F}{1-h_F^2/c^2}-V_{M-2,J}(\tau)+4V_{M-1,J}(\tau)\bigg). \end{align} $$
At
$\alpha =\alpha _{\max }$
, the condition (4.2) is discretized as
Then
Augmenting
$\mathbf {U}(\tau )$
by the vector of option values when
$\alpha =\alpha _{\max }$
and denoting
we obtain the full system of differential equations
$$ \begin{align} \hat{\mathbf{U}}'(\tau)=\hat{\mathbf{A}}\hat{\mathbf{U}}(\tau)+\hat{\mathbf{b}}(\tau), \end{align} $$
where
$\hat {\mathbf {A}}\in \mathbb {R}^{(M-1)J\times (M-1)J}$
and
$\hat {\mathbf {b}}(\tau ) \in \mathbb {R}^{(M-1)J}$
are modifications brought to the matrix
$\mathbf {A}$
and vector
$\bar {\mathbf {b}}(\tau )$
in (4.4) to account for the boundary conditions at
$F=F_{\max }$
and
$\alpha =\alpha _{\max }$
given by (4.5) and (4.6), respectively.
4.1 Time-integration scheme
There are different possibilities for the time integration of (4.7). A one-step exponential time integration [Reference Rambeerich, Tangman, Lollchund and Bhuruth23] would require the computation of an exponential matrix. This can be efficiently carried out via the Carathéodory–Fejer approximation and contour integrals [Reference Schmelzer and Trefethen25]. For simplicity, an implicit Euler time stepping is described below. Using N time steps in
$[0,T]$
, let
$k=T/N$
denote the step size with time nodes
$\tau _n=nk$
,
. Let
and
$\hat {\mathbf {b}}^n=\hat {\mathbf {b}}(\tau _n)$
. An implicit Euler discretization of (4.7) gives the time-stepping procedure
with initial condition
$\hat {\mathbf {U}}^0$
determined by the payoff function as previously described.
4.2 Numerical stability
Problem (4.1) is a two-dimensional diffusion equation with a mixed derivative term. For multidimensional problems, but with constant coefficients and periodic boundary conditions, in’t Hout and Welfert [Reference in’t Hout and Welfert19] have presented a stability analysis for second-order ADI schemes. The numerical method in this work uses an implicit Euler time stepping, but variable coefficients and nonperiodic boundary conditions for the problem considered here will bring difficulties in a matrix method for carrying out the stability analysis.
Another approach for the stability analysis would be to use the principle of frozen coefficient problems. Von Neumann analyses [Reference Charney, Fjørtoft and Von Neumann5] can then be carried out for constant-coefficient problems which are obtained by fixing the coefficients at their values at each point in the computational domain [Reference Mishra and Svard22, Reference Strikwerda28]. The variable-coefficient problem is stable if each frozen coefficient problem is stable. This would require finding the amplification factor for the fully discrete scheme.
The implicit Euler discretization of (4.3) can be written as
$$ \begin{align} V_{ij}^{n+1}=V_{ij}^n+k\hat{a}_{ij}D_{FF}V_{ij}^{n+1}+k\hat{b}_{ij}D_{F\alpha} V_{ij}^{n+1}+k\hat{c}_jD_{\alpha\alpha}V_{ij}^{n+1},\quad n\ge 0. \end{align} $$
Let
$$ \begin{align} V_{ij}^n=g^n (\xi_1h_F,\,\xi_2h_{\alpha})e^{\imath (i\xi h_F+j \xi_2h_{\alpha})}, \end{align} $$
where
$\imath $
is the unit imaginary complex number,
$\xi _1$
and
$\xi _2$
are wave numbers and
$g(\xi _1h_F,\,\xi _2h_{\alpha })$
is the amplification factor. The stability condition would require that the amplification factor satisfies
. Substituting (4.9) in (4.8),
where
$$ \begin{align*} \gamma_F=\frac{k}{h_F^2}\bigg(1+\frac{h_F^2}{c^2}\bigg),\quad \gamma_{\alpha}=\frac{k}{h_{\alpha}^2}\bigg(1+\frac{h_{\alpha}^2}{c^2}\bigg),\quad \gamma_{F\alpha}=\frac{k}{4h_{\alpha}h_F}\bigg(1+\frac{3h_F^2}{2c^2} +\frac{3h_{\alpha}^2}{2c^2}-\frac{9h_F^2h_{\alpha}^2}{4c^4}\bigg).\end{align*} $$
We assume that
$\beta>0$
and
$F\ge 0$
. Since
$\hat {a}_{ij}=\alpha _j^2F_i^{2\beta }/2$
and
$\hat {c}_j=\nu ^2\alpha _j^2/2$
, we find that both
$\hat {a}_{ij}$
and
$\hat {c}_j$
are nonnegative. When the correlation coefficient
$\rho>0$
, the coefficient
$\hat {b}_{ij}=\rho \nu \alpha _j^2F_i^{\beta }$
is also nonnegative. Therefore, for a positive correlation coefficient
$\rho $
, the amplification factor (4.10) satisfies the stability condition
.
When the correlation
$\rho <0$
, establishing that
is not straightforward. In this case, in order to validate the stability of the numerical scheme, an experimental study is described in the next section. The results show that the errors in the numerical solution remain small with larger time steps. Numerical investigation of the magnitude of the amplification factor for different parabolic mesh ratios also did not reveal any violation of the inequality
.
5 Numerical results
In this section, some numerical examples are considered to illustrate the performance of the proposed scheme in achieving arbitrage-free and accurate prices. All numerical experiments have been performed using Matlab R2017a on a Core i5 laptop with 16 GB RAM and speed 3.60 GHz. In all the numerical experiments, we choose
$F_{\max }=4E$
and
$\alpha _{\max }=4 \hat {\alpha }$
unless stated otherwise.
5.1 No-arbitrage pricing
The first numerical example draws a comparison of computed implied volatilities obtained using the RBF-FD method, the no-arbitrage closed-form (No-Arb CF) approximation method [Reference Yang, Chen, Liu and Wan33] and the standard approximation (Std SABR) [Reference Hagan, Kumar, Lesniewski and Woodward16] to demonstrate that the proposed method is arbitrage-free. Figure 1 shows the implied volatility curve for a low-maturity problem with
$T=1$
. The other parameters are given by an initial volatility of
$\hat \alpha =0.1$
, a CEV elasticity of
$\beta =0.1$
, with zero correlation
$\rho =0$
, a volatility-of-volatility (vol-vol) of
$\nu =0.1$
and the initial forward price is chosen as
$f=0.05$
[Reference Yang, Chen, Liu and Wan33]. The results are benchmarked with a Monte Carlo (MC) simulation technique with 100 000 time steps. Note that the implied volatilities computed using the no-arbitrage closed-form approximation method and the RBF-FD method agree closely with those computed using the MC simulation method, while the computed implied volatilities of the standard SABR model [Reference Hagan, Kumar, Lesniewski and Woodward16] deviate from the MC results for low strikes. Yang et al. [Reference Yang, Chen, Liu and Wan33] mentioned that their approximations work well when the total vol-vol
$(\nu \sqrt {T})$
parameter is small, but, as T becomes large, the formula leads to biased results. We therefore consider a longer maturity problem in the next example.
Figure 1 Implied volatility for low-strikes problem (colour available online).
Doust [Reference Doust10] showed that the implied volatility approximation [Reference Hagan, Kumar, Lesniewski and Woodward16] can result in a negative probability density function for long-dated options. Using the same set of parameters,
$f=0.0488$
with the model parameters given by
$(\hat \alpha ,\beta ,\rho ,\nu )=(0.026,0.5,-0.1,0.4)$
, we give the no-arbitrage density function obtained by the RBF-FD multiquadric (MQ) method in Figure 2 and the RBF-FD inverse multiquadric (IMQ) method, which corresponds to the case
$m=-1$
(3.1), in Figure 3. Note that both methods are always positive, whereas that of the standard SABR [Reference Hagan, Kumar, Lesniewski and Woodward16] goes negative for low strikes. Comparison is carried out against the no-arbitrage one-dimensional PDE (No-Arb 1D-PDE) method [Reference Hagan, Kumar, Lesniewski and Woodward17] and the no-arbitrage closed-form (No-Arb CF) approximation by Yang et al. [Reference Yang, Chen, Liu and Wan33]. The corresponding implied volatilities over the range of strike prices
$[0, 0.14]$
are given in Figure 4. We observe that while the standard SABR [Reference Hagan, Kumar, Lesniewski and Woodward16] deviates as the strike price decreases, there is a close agreement between the RBF-FD and the No-Arb 1D-PDE methods.
Figure 2 Positive probability density using the RBF-FD MQ method.
Figure 3 Positive probability density using the RBF-FD IMQ method.
Figure 4 Implied volatilities for different strike prices.
5.2 Option prices: uncorrelated case
Table 1 gives the results for pricing a call option with a strike price of
$K=1.1$
and a very short maturity of one month
$(T=1/12)$
in the uncorrelated case. The initial forward price is taken as
$f=1.1$
. In this setting, the expansion formula [Reference Hagan, Kumar, Lesniewski and Woodward16] and the analytical price approximation [Reference Yang, Chen, Liu and Wan33] are expected to be highly accurate, since T,
$\hat \alpha f^{\beta -1}\nu $
and
are all small. The results indicate that the RBF-FD method is in close agreement with the methods using the expansion formula and the analytical price approximation. For the option with a longer maturity of
$T=5$
years, it can be seen that the numerical method has a smooth second-order convergence. The reference price comes from computations using a much refined mesh in the F-direction.
Table 1 Uncorrelated case.
The next numerical example tests the accuracy of the RBF-FD scheme for the special case when the expansion formula [Reference Hagan, Kumar, Lesniewski and Woodward16] is not accurate enough. Parameters in the model (2.1) are given by
$\rho =0$
,
$\beta =0.8$
and
$\nu =0.4$
and initial forward price and volatility are
$f=1.1$
and
$\hat \alpha =0.3$
, respectively. The maturity T is allowed to vary in the range
$[1,10]$
for a strike price of
$K=1.1$
. The results are benchmarked against the exact simulation method [Reference Cai, Song and Chen4]. It is observed that as T increases, the expansion formula of Hagan suffers from an increasing upward bias while the proposed no-arbitrage RBF-FD PDE approach prices coincide with the results of the exact simulation method throughout the range of values of the maturity T.
5.3 Option prices: correlated case
Computed RBF-FD prices are reported for the nonzero correlation case with
$\rho =-0.25$
in Table 2. Other parameters are the same as those computed using the continuous-time Markov chain (CTMC) approximation [Reference Cui, Kirkby and Nguyen9]. The RBF-FD prices are also compared against the analytical approximation [Reference Yang, Chen, Liu and Wan33] and the exact simulation method [Reference Song27]. The results demonstrate close agreement of RBF-FD prices with those given by CTMC and the exact simulation methods. Also note that the analytical approximation method [Reference Yang, Chen, Liu and Wan33] yields less accurate prices.
Table 2 Comparison of option prices in the correlated case.
5.4 Optimal shape parameter
The selection of the shape parameter in a RBF method plays an important influence on the efficiency and accuracy of the numerical method. An algorithm for selection of the shape parameter based on the leave-one-out cross-validation method was proposed by Rippa [Reference Rippa24] and extensions of this approach were described by Fasshauer and Zhang [Reference Fasshauer and Zhang14]. The choice of the optimal shape parameter for PDE problems has been recently considered by Chen et al. [Reference Chen, Hong and Lin7] and Wang et al. [Reference Wang, Chen, Zhang and Hua32]. An algorithm based on minimization of local approximations to the RBF truncation error has been developed by Bayona et al. [Reference Bayona, Moscoso and Kindelan2]. However, we believe that adapting this minimization algorithm for determining the optimal shape parameter for the problem considered in this work would require a separate study and we have not pursued it further.
For studying the influence of the shape parameter on the accuracy of the proposed method, we have studied the dependence of the
$\ell ^2$
error on the shape parameter and the results of this experimental study for different values of c are displayed in Figure 6. It is observed that the error decreases until a critical value of
$c\approx 1.5$
is reached and, at this approximate value, the error attains a least value. Further increase in c results in a flattening of the
$\ell ^2$
error curve. It is therefore observed that a value of
$c=1.5$
in the numerical method is a good choice.
Figure 5 Accuracy comparison of the different methods with the exact simulation method when the expansion formula [Reference Hagan, Kumar, Lesniewski and Woodward16] is unreliable.
Figure 6 Dependence of the
$\ell ^2$
error on the shape parameter c (colour available online).
5.5 Numerical tests for stability
A numerical approach for the detection of stability restrictions for option pricing PDEs under stochastic volatility models was employed by Düring and Fournié [Reference Düring and Fournié12]. A similar approach is used here and the
$\ell ^2$
errors of the numerical solutions are computed for varying values of the parabolic mesh ratios
$\gamma _F$
and
$\gamma _{\alpha }$
.
European call option prices for the uncorrelated and correlated cases are calculated using the parameters in Table 1 with
$T=0.5$
and in Table 2. The corresponding stability plots for the
$\ell ^2$
error are given in Figures 7(a) and 7(b), respectively. The stability plots are obtained for
$\gamma _F$
varying in the range
and a sequence of spatial grid sizes
for the uncorrelated case and
for the correlated case. The results displayed in Figures 7–9 show that there is no stability restriction on the time-step size and that accurate numerical solutions are computed in all the cases as evidenced by the small
$\ell ^2$
errors. From Figures 8(a) and 8(b), the same conclusion is reached when the parameter
$\gamma _\alpha $
is varied.
Figure 7 Contour plots for the
$\ell ^2$
error norm against the parabolic mesh ratio
$\gamma _F$
against the spatial mesh
$h_F$
in the F-direction.
Figure 8 Contour plots for the
$\ell ^2$
error norm against the parabolic mesh ratio
$\gamma _\alpha $
against the spatial mesh
$h_\alpha $
in the
$\alpha $
-direction.
Figure 9 Surface plot for the
$\ell ^2$
error against parabolic mesh ratios
$\gamma _F$
and
$\gamma _\alpha $
.
Figures 9(a) and 9(b) show the
$\ell ^2$
error as a function of
$(\gamma _F,\,\gamma _\alpha )$
. We observe that the errors remain small and this validates our earlier analysis that the numerical scheme is unconditionally stable.
5.6 Barrier options
Yang et al. [Reference Yang, Liu and Cui34] derived closed-form approximation formulas for pricing continuously monitored down-and-out and up-and-out options. The prices given in their work are used to demonstrate that the proposed RBF-FD method is also an accurate and efficient method for the pricing of barrier options under the SABR dynamics. The parameters in the SABR model are
$\beta =0.5$
and
$\nu =0.1$
, the initial forward price is taken as
$f=100$
and the initial volatility is
$\hat {\alpha }=0.1$
. A down-and-out call option with strike price
$K=100$
is priced with a lower barrier at
$L=98$
using RBF-FD and the results are given in Table 3. Reference prices are from a MC method via Euler discretization with 1 000 000 simulation samples [Reference Yang, Liu and Cui34]. The table also shows relative errors (RE) computed using the formula
$$ \begin{align*}\text{RE}=\frac{|V_{\text{Approx}}-V_{\text{MC}}|}{V_{\text{MC}}}.\end{align*} $$
Table 3 Down-and-out call option prices for varying values of
$\rho $
.
Table 3 shows that the relative errors are all less than
$1\%$
and that the RBF-FD method slightly outperforms the results obtained by the closed-form approximation method in the case of small maturities of
$T=0.5$
and
$T=1$
year.
Figure 10(a) shows that there is a close agreement between the RBF-FD method and the closed-form approximation for
$T=1$
. However, for a longer maturity of
$T=5$
years, the authors reported a relative error of
$5.0\%$
with MC prices and that the relative error increases as maturity increases. They further mentioned that their formula is not valid for relatively large maturities, since the approximation error is at first order for a positive barrier level. This is illustrated in Figure 10(b), where we observe that as the maturity increases, the closed-form analytical approximation price deviates from the RBF-FD price.
Figure 10 Comparison of down-and-out call options for small- and large-maturity problems.
In Table 4, we report the values of double knock-out call options for different values of
$\beta $
with knock-out boundaries at
$L=90$
and
$U=120$
. The analytical approximation approach of Yang et al. [Reference Yang, Liu and Cui34] works only for single-barrier options and the authors mentioned that double-knock-out options must be tackled differently. We therefore do not have any benchmark results for comparison of the RBF-FD computed prices, but, given the high accuracy of the proposed method for pricing single-knock-out options, it is expected that the results reported here are also highly accurate.
Table 4 Double knock-out call options prices using the RBF-FD method.
5.7 Randomly scattered nodes
We finally consider the pricing of options using the RBF-FD method on scattered nodes by generating Halton points constructed from the Van der Corput sequences [Reference Fasshauer15]. Figure 11 shows the Halton points constructed using 100 interior points in a unit square in
$\mathbb {R}^2$
, and along the boundaries the nodes are partitioned uniformly.
Figure 11 Halton grid [Reference Fasshauer15] with
$M=100$
interior points for
$F\in [0,\,2]$
and
$\alpha \in [0,\,1]$
.
For the parameters
$f=1,\,E=1,\,\hat \alpha =0.2,\,\beta =0.5,\,\rho =-0.06,\,\nu =0.2$
and
$T = 0.5$
with
$F_{\max }=2$
and
$\alpha _{\max }=1$
, the interpolated option prices on the scattered nodes and the resulting fitted surface are shown in Figure 12.
Figure 12 Option price fitted surface for
$F\in [0,\,2]$
and
$\alpha \in [0,\,1]$
.
Finally, a comparison is made between the magnitude of the error norms obtained by the following three methods: a RBF-FD method using randomly scattered grid points, a RBF-FD method on a uniform grid and a standard finite-difference (FD) scheme using central difference approximations. The FD and RBF-FD results are based on the same grid construction. From the results given in Table 5, it is observed that choosing the same range for F for both the FD and RBF-FD schemes gives numerical errors of similar magnitudes. It is also observed that increasing the number of support nodes
$(ns)$
for the random grid method leads to more accurate prices.
Table 5 Option prices obtained on scattered nodes and structured nodes.
6 Conclusion
Valuation of options under the popular SABR model using asymptotic implied volatility expansions or analytical approximations using expansions of Bessel operators are unlikely to yield accurate prices for large maturities. Our work addresses this lack of accuracy problem by proposing a novel approach based on the direct numerical solution of the two-dimensional pricing equation with an absorbing boundary condition at zero. The proposed method employed radial basis functions approximations in a finite-difference mode of the spatial derivatives.
New formulas for approximating the first- and second-order derivatives using generalized multiquadrics were derived, and these were shown to recover the approximations for multiquadrics and inverse multiquadrics derived in other studies. The system of ordinary differential equations arising in the computational method was formulated and time integrated using an implicit Euler scheme. Various examples of option pricing problems under SABR were numerically solved. Confirmation of no-arbitrage prices was established by the validation of nonnegative probability densities computed using multiquadrics and inverse multiquadrics.
Comparisons with existing MC simulation methods showed that accurate European and barrier prices were computed by the new RBF-FD method. One clear advantage of the RBF-FD method over standard finite-difference approaches is that the method can be applied in a meshless approach. Using randomly scattered points in the computational domain, the method was shown to yield accurate prices. As such, the method proposed in this work provides an efficient approach for option pricing under SABR. In a follow-up work, the pricing of American options using the RBF-FD approach will be studied.
Acknowledgements
The author would like to thank the two anonymous reviewers for their valuable suggestions, which have brought significant and important enhancements to this work.
