A volatility-of-volatility expansion of the option prices in the SABR stochastic volatility model

TL;DR

The paper introduces a rapid, explicit series expansion method for approximating solutions to the SABR PDE, including derivatives and implied volatility, with improved simplicity and applicability over existing approaches.

Contribution

It develops a small-parameter series expansion for the SABR PDE that yields explicit, fast-to-evaluate approximations for option prices and implied volatility, including mean reversion effects.

Findings

The method provides accurate approximations compared to Hagan's formula.

It efficiently computes derivatives of the solution.

Numerical tests confirm good performance with market data.

Abstract

We propose a general, very fast method to quickly approximate the solution of a parabolic Partial Differential Equation (PDEs) with explicit formulas. Our method also provides equaly fast approximations of the derivatives of the solution, which is a challenge for many other methods. Our approach is based on a computable series expansion in terms of a "small" parameter. As an example, we treat in detail the important case of the SABR PDE for $\beta = 1$, namely $\partial_{\tau}u = \sigma^2 \big [ \frac{1}{2} (\partial^2_xu - \partial_xu) + \nu \rho \partial_x\partial_\sigma u + \frac{1}{2} \nu^2 \partial^2_\sigma u \, \big ] + \kappa (\theta - \sigma) \partial_\sigma$, by choosing $\nu$ as small parameter. This yields $u = u_0 + \nu u_1 + \nu^2 u_2 + \ldots$, with $u_j$ independent of $\nu$. The terms $u_j$ are explicitly computable, which is also a challenge for many other, related methods. Truncating this expansion leads to computable approximations of $u$ that are in "closed form," and hence can be evaluated very quickly. Most of the other related methods use the "time" $\tau$ as a small parameter. The advantage of our method is that it leads to shorter and hence easier to determine and to generalize formulas. We obtain also an explicit expansion for the implied volatility in the SABR model in terms of $\nu$, similar to Hagan's formula, but including also the {\em mean reverting term.} We provide several numerical tests that show the performance of our method. In particular, we compare our formula to the one due to Hagan. Our results also behave well when used for actual market data and show the mean reverting property of the volatility.

Tables

Table 1. Table 1. Comparison of the implied volatilities of the models with the market data

Obj.	ISE	OSE	$\overline{\nu }$	$std\left(\nu \right)$	$\overline{\sigma }$	$std\left(\sigma \right)$	$\overline{\rho }$	$std\left(\rho \right)$
${\sigma }_D$	.0152	.0179	1.335	0.3	.1889	.078	-0.54	.07
${\sigma }_H$	.0154	.0181	1.332	0.29	.1902	.079	-0.58	.069

Table 2. Table 2. Comparison of the actual prices of the models with the market data

Obj.	ISE	OSE	$\overline{\nu }$	$std\left(\nu \right)$	$\overline{\sigma }$	$std\left(\sigma \right)$	$\overline{\rho }$	$std\left(\rho \right)$
$C_H$	.0038	.0043	1.0463	.24	.197	.07	-0.56	.07
${\tilde{C}}_H$	.0036	.0042	1.0687	.28	.192	.07	-0.39	.15
$C_D$	.0037	.0042	1.0502	.25	.196	.07	-0.53	.08
$C_{SA,2}$	.0037	.0043	1.0566	.25	.196	.07	-0.55	.07

Table 3. Table 3. Comparison of the log prices of the models with the market data

Obj Funct	ISE	OSE	$\overline{\nu }$	$std\left(\nu \right)$	$\overline{\sigma }$	$std\left(\sigma \right)$	$\overline{\rho }$	$std\left(\rho \right)$
$\ln \left(C_H\right)$	.069	.077	1.2	.2	.19	.08	-0.57	.06
$\ln \left({\tilde{C}}_H\right)$	.056	.066	1.4	.4	.18	.08	-0.28	.12
$\ln \left(C_D\right)$	.059	.068	1.5	.3	.19	.08	-0.54	.04
$\ln \left(C_{SA,2}\right)$	.067	.076	1.3	.3	.19	.08	-0.55	.05
$\ln \left(C_{\kappa }\right)$	.064	.074	1.2	.3	.19	.08	-0.55	.04
$\ln \left(C_{BS}\right)$	.119	.123			.18	.06

Table 4. Table 4. Estimate of the residuals R := ‖ ( ∂ t − L ) ​ C ‖ 2 assign 𝑅 subscript norm subscript 𝑡 𝐿 𝐶 2 R:=\|(\partial_{t}-L)C\|_{2}

$C=$	$C_H$	$C_D$	$C_{SA,2}$	$C_{BS}$
${10}^{3}R$	0.489	0.181	0.163	16.436

Table 5. Table 5. Estimates of residuals (continued)

$C=$	$C_H$	$C_D$	$C_{SA,2}$	$C_{BS}$	$\nu$	$T\in$	$y\in$
${10}^{2}R$	.33	.07	.072	1.8	.1	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}30}]$	$[-0.3,.3]$
${10}^{2}R$	.38	.28	.37	1.8	.1	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}30}]$	$[-1.5,\mathrm{\hspace{0.17em}1.5}]$
${10}^{2}R$	.02	.016	.018	1.4	.25	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.2}]$	$[-1.5,\mathrm{\hspace{0.17em}1.5}]$
${10}^{2}R$	3.2	2.4	5.2	18.8	1	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}1}]$	$[-0.2,\mathrm{\hspace{0.17em}0.2}]$
${10}^{2}R$	4.3	26.	14.8	15.5	1	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}1}]$	$[-1,\mathrm{\hspace{0.17em}1}]$
${10}^{2}R$	6.1	4.	7.	17.5	1	$[\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}2}]$	$[-0.2,\mathrm{\hspace{0.17em}0.2}]$

Table 6. Table 6. Comparison of Dyson series method with Hagan’s formula

Obj Funct	${\ell }_1$	${\ell }_2$	${\ell }_{\infty }$
${\sigma }_H-{\sigma }_D$	.0057	.0139	.167
$C_H-C_D$	.0010	.0026	.0266
$\ln \left(C_H\right)-\ln \left(C_D\right)$	.0326	.0938	3.598

Table 7. Table 7. Comparison of the closed form approximate solutions with the iterative FD approximations and ρ = − .2 𝜌 .2 \rho=-.2 . All norms were multiplied by 100.

T	$\nu$	$\Vert \delta C_H\Vert$	${\Vert \delta C_H\Vert }_{\infty }$	log $C_H$	$\Vert \delta C_{SA,2}\Vert$	${\Vert \delta C_{SA,2}\Vert }_{\infty }$	log $C_{SA,2}$	$\Vert {\delta }_{HD}\Vert$	est error
5	1	10.5	22.7	74.3	7.35	16.8	55.3	5.32	.036
2	1.5	6.32	14.2	78.	2.91	7.72	56.7	5.9	.0102
2	1	1.45	3.28	35.3	0.939	2.26	38.6	1.72	.0091
2	.5	.1	0.23	4.03	.136	.398	6.68	.179	.0022
1*	1.5	1.25	2.81	43.9	.741	1.6	56.4	1.7	.0442
1	1.5	1.24	2.9	43.3	.732	1.61	56.4	1.69	.011
1	1	.241	.53	14.3	.2379	.608	22.2	0.403	.0101
.5	1	.029	.088	2.13	.051	.179	4.83	.0665	.003

Table 8. Table 8. Comparison of the closed form approximate solutions with the iterative FD approximations and ρ = − .5 𝜌 .5 \rho=-.5 . All norms were multiplied by 100.

T	$\nu$	$\Vert \delta C_H\Vert$	${\Vert \delta C_H\Vert }_{\infty }$	log $C_H$	$\Vert \delta C_{SA,2}\Vert$	${\Vert \delta C_{SA,2}\Vert }_{\infty }$	log $C_{SA,2}$	${\delta }_{HD}$	FD-err	cut-off
2	1	1.24	2.88	16.8	1.03	2.47	39.4	1.54	.0185	.0418
1	1	.253	.543	2.4	.314	.752	25.9	.462	.0212	.0016
0.5	1	.0497	.137	1.41	.0671	.177	13.7	.0936	.0058	.0012

Equations

$$\begin{cases}\ \partial_{\tau}u\,=\,\sigma^{2}\big{[}\frac{1}{2}(\partial^{2}_{x}u-\partial_{x}u)+\nu\rho\partial_{x}\partial_{\sigma}u+\frac{1}{2}\nu^{2}\partial^{2}_{\sigma}u\,\big{]}\ +\ \kappa(\theta-\sigma)\partial_{\sigma}&\\ \ u(x,\sigma,0)=v(x,\sigma)\,.&\\ \end{cases}$$

$${\partial }_{\tau }u-\frac{1}{2}{\sigma }^2({\partial }_x^2-{\partial }_x)u=0.$$

$$u=u_0+\nu u_1+{\nu }^2{u}_2+\dots +{\nu }^k{u}_k+O({\nu }^{k+1}).$$

$$h(x,\sigma )=|e^x-K{|}_{+}:=\max \{e^x-K,0\}.$$

$$F_{SA}(S,K,0,\sigma ,\rho ,\tau )=F_{BS}(S,K,\sigma ,\tau ),$$

$$F_{SA}(S,K,\nu ,\sigma ,\rho ,\tau )=F_{BS}(S,K,\sigma ,\tau )+\nu F_1+{\nu }^2{F}_2+\dots +{\nu }^k{F}_k+O({\nu }^{k+1}).$$

$${\partial }_{t}u=Lu,u(0)=v.$$

$$e^{tL}v=I_0+I_1+\dots +I_{n-1}+{\mathcal{E}}_n,$$

$$I_k:={\int }_0^t{\int }_0^{{\tau }_1}\dots {\int }_0^{{\tau }_{k-1}}{e}^{(t-{\tau }_1)L_0}V{e}^{({\tau }_1-{\tau }_2)L_0}\dots V{e}^{({\tau }_{k-1}-{\tau }_k)L_0}V{e}^{{\tau }_k{L}_0}vd\tau ,$$

$$e^{(t-{\tau }_1)L_0}V{e}^{({\tau }_1-{\tau }_2)L_0}V\dots e^{({\tau }_{k-1}-{\tau }_k)L_0}V{e}^{{\tau }_k{L}_0}=e^{t{L}_0}{P}_k(t,\tau ),$$

$$d_{\pm }:=\frac{\ln (S{e}^{rt}/K)}{\sigma \sqrt{t}}\pm \frac{\sigma \sqrt{t}}{2}=\frac{x-\ln K}{\sigma \sqrt{t}}\pm \frac{\sigma \sqrt{t}}{2}$$

$$F_{1}(S,K,\sigma,t)\,=\,\frac{Kt}{2}\big{(}\kappa_{0}(\theta-\sigma)\sqrt{t}-\rho\sigma d_{-}\big{)}N^{\prime}(d_{-})\,.$$

$${\sigma }_{imp}=\sigma +\nu e_1+{\nu }^2{e}_2+O({\nu }^3).$$

$$e_1=-\rho d_{-}\sigma \sqrt{t}/2$$

$$F_{SA,2}(S,K,\nu ,\sigma ,\rho ,t)=F_{BS}(S,K,\sigma ,t)+\nu F_1(S,K,\sigma ,t)+{\nu }^2{F}_2(S,K,\sigma ,t),$$

$$F_D(S,K,\nu ,\sigma ,\rho ,t)=F_{BS}(S,K,\sigma +\nu e_1+{\nu }^2{e}_2,t),$$

$$C_M(S,K,t)=C_{BS}(S,K,{\sigma }_{imp},t),$$

$$e^{tL}v=e^{t{L}_0}v+{\int }_0^t{e}^{(t-\tau )L_0}V{e}^{\tau L}vd\tau .$$

$${\mathcal{E}}_n:={\int }_0^t{\int }_0^{{\tau }_1}\dots {\int }_0^{{\tau }_{n-1}}{e}^{(t-{\tau }_1)L_0}V{e}^{({\tau }_1-{\tau }_2)L_0}\dots e^{({\tau }_{n-1}-{\tau }_n)L_0}V{e}^{{\tau }_{n}L}vd\tau .$$

$$P_m(L_0,V;\theta ):=\sum _{k=0}^m\frac{{\theta }^k}{k!}{\mathrm{ad}}_{L_0}^k(V)=L+\theta [L_0,V]+\frac{{\theta }^2}{2}[L_0,[L_0,V]]+\dots .$$

$$e^{\theta L_0}V=P_m(L_0,V;\theta )e^{\theta L_0},$$

$$V{e}^{\theta L_0}=e^{\theta L_0}{P}_m(L_0,V;-\theta ).$$

$$\begin{gathered} e^{(t-{\tau }_1)L_0}V{e}^{({\tau }_1-{\tau }_2)L_0}\dots e^{({\tau }_{k-1}-{\tau }_k)L_0}V{e}^{{\tau }_{k}L} \\ =e^{t{L}_0}{P}_m(L_0,V;-{\tau }_1)\dots P_m(L_0,V;-{\tau }_{k-1})P_m(L_0,V;-{\tau }_k) \end{gathered}$$

$$\begin{array}{c}{L}_0=\frac{1}{2}{\sigma }^2({\partial }_x^2-{\partial }_x),L_1=\rho {\sigma }^2{\partial }_x{\partial }_{\sigma }+{\kappa }_0(\theta -\sigma ){\partial }_{\sigma },L_2=\frac{1}{2}{\sigma }^2{\partial }_{\sigma }^2,\text{mbox}and\\ L:=L_0+\nu L_1+{\nu }^2{L}_2.\end{array}$$

$$e^{tL}=e^{t{L}_0}+\nu B_1+{\nu }^2{B}_2+\dots {\nu }^{n-1}{B}_{n-1}+{\nu }^n{\mathcal{R}}_n,$$

$$\begin{array}{c}{J}_1:={\int }_0^t{e}^{(t-{\tau }_1)L_0}{L}_1{e}^{{\tau }_1{L}_0}d{\tau }_1,\\ J_2:={\int }_0^t{e}^{(t-{\tau }_1)L_0}{L}_2{e}^{{\tau }_1{L}_0}d{\tau }_1,\text{mbox}and\\ I(T_1,T_2):={\int }_0^t{\int }_0^{{\tau }_1}{e}^{(t-{\tau }_1)L_0}{T}_1{e}^{({\tau }_1-{\tau }_2)L_0}{T}_2{e}^{{\tau }_2{L}_0}d{\tau }_{2}d{\tau }_1\end{array}$$

$$\begin{gathered} {\mathcal{R}}_2=I(L_1,L_2)+I(L_2,L_1)+\nu I(L_2,L_2) \\ +{\int }_0^t{\int }_0^{{\tau }_1}{\int }_0^{{\tau }_2}{e}^{(t-{\tau }_1)L_0}V{e}^{({\tau }_1-{\tau }_2)L_0}V{e}^{({\tau }_2-{\tau }_3)L_0}V{e}^{{\tau }_{3}L}d{\tau }_{3}d{\tau }_{2}d{\tau }_1. \end{gathered}$$

$$e^{tL}=e^{t{L}_0}+\nu J_1+{\nu }^2(J_2+I(L_1,L_1))+{\nu }^3{\mathcal{R}}_2.$$

$$B_1=J_1\text{mbox}and{B}_2=J_2+I(L_1,L_1).$$

$$e^{tL}v=e^{t{L}_0}v+\nu J_{1}v+{\nu }^2(J_{2}v+I(L_1,L_1))+{\nu }^3{\mathcal{R}}_{2}v,$$

Full text

The opinions expressed in this paper are those of the authors and do not necessarily reflect the views of the Board Governors of the Federal Reserve System or any other individuals within the Federal Reserve System. Manuscripts available from **http://www.math.psu.edu/nistor/. **

1. 1 The Duhamel-Dyson perturbative series expansion

2. 1.1 General results: The Duhamel and Campbell-Hausdorff-Backer formulas

3. 1.2 Commutator calculations for SABR

4. 1.3 Expansion in the powers of $\nu$

5. 1.4 Symbolic formulas for $B_{j}$, $j=1,2$.

6. 1.5 Heat kernels and convolutions

7. 1.6 Calculation of the differential operators $\tilde{G}_{j}$

8. 1.7 Remarks on the implementation

9. 2 Applications: approximations of the derivatives and of the implied volatility

10. 2.1 Approximation of derivatives

11. 2.2 Implied Volatility

12. 3 Model calibration and market tests

13. 3.1 Description of the data

14. 3.2 General description of the method

15. 3.3 The first type of market data tests: implied volatility

16. 3.3.1 Objective function: $\sigma_{D}$

17. 3.3.2 Objective function: $\sigma_{H}$

18. 3.3.3 Conclusion

19. 3.4 The second type of market data tests: actual prices

20. 3.4.1 Objective functions: $C_{D}$ and $C_{SA,2}$

21. 3.4.2 Objective functions: $C_{H}$ and $\tilde{C}_{H}$

22. 3.4.3 Conclusion

23. 3.5 The third type of market data tests: log-prices

24. 4 Numerical tests

25. 4.1 The residual of the approximations: substituting in the PDE

26. 4.2 Comparison of our implied volatility formula with Hagan’s formula for maket data

27. 4.3 A Monte Carlo simulation

28. 4.4 Comparison with the Finite Difference approximate solution

29. 4.4.1 Outline of the method and of its challenges

30. 4.4.2 Cut-off errors

31. 4.4.3 Discretization errors

32. 4.4.4 Total error estimation

33. 5 Extensions of the method

34. 5.1 The log-normal model with mean-reverting volatility

35. 5.2 The SABR model (with general $\beta$)

Introduction

We propose a new, general, non-iterative method to construct fast approximations to suitable parabolic Partial Differential Equations by constructing closed form approximate formulas for their solutions. In addition to being very fast, our method has also the advantage that it can be used to approximate also the derivatives of the solution. In a nut-shell, our method is to first use a Dyson perturbative series to expand the solution $u$ of our parabolic PDE in terms of a “small” parameter and then to explicitly compute the terms of the series expansion. Dyson series expansions have long been used in calculations; our achievement is to realize such a series with computable terms.

Our method expands the scope of the general perturbative-type method devised and used in [10, 14, 35]. We illustrate our method by providing closed form approximations to the solution $u$ of the SABR partial differential equation (PDE) with mean reversion (the $\lambda$SABR model) that arises in option pricing. Being able to treat the mean reversion term in the SABR PDE is another advantage of our method. Whereas in the aforementioned papers the small parameter is the time, in this paper, the small parameter is the “volatility-of-volatility” parameter $\nu$ (see Equation (1)). This sets apart our approach from the previous papers using series expansions, but see also [27, 31, 32, 33, 49, 44, 54]. The Dyson perturbative series expansion is obtained by iterating Duhamel’s formula, so we will call it the Duhamel-Dyson perturbative series expansion. See below for more on the earlier results, in general, and on the Duhamel-Dyson perturbative series expansions, in particular.

We are interested in this paper in obtaining concrete results, so we tried to keep the theoretical considerations (including proofs) to a minimum, striving instead to carefully explain our method and to obtain as quickly as possible our formulas, which we then tested numerically.

The mathematical problem: the $\lambda$SABR PDE

While our method is quite general, since we want to illustrate and test our method through explicit results, we shall concentrate mainly in this paper on the PDE

[TABLE]

This PDE is obtained from the usual SABR PDE [31, 32] by using the change of variables $x=\ln(Se^{rt})$ and $\tau=T-t$, by fixing $\beta=1$, and by including the standard mean reverting term $\kappa(\theta-\sigma)\partial_{\sigma}$. We shall use the term $\lambda$SABR PDE for Equation (1). When $\kappa=0$, we shall also use the term SABR PDE.

In order to set up our perturbative method, we write $\kappa=\nu\kappa_{0}$, for some fixed parameter $\kappa_{0}$. That is, we consider $\kappa$ of the same order as $\nu$. Therefore, when $\nu=0$, the SABR PDE reduces to the forward Black-Scholes PDE (for the forward price and in $x=\ln(Se^{r\tau})$)

[TABLE]

(here $\tau=T-t$). If $u_{0}$ is the solution of the Black-Scholes equation, our method provides an expansion of the form

[TABLE]

For general initial data $v$ in Equation (1), the terms $u_{j}$ can be computed using an integration of $v$ with respect to the Green functions, but this is not fully explicit (see Remarks 1.7 and 1.10). To obtain fully explicit (or “closed form”) results, we specialize our initial condition $v(x,\sigma)=u(x,\sigma,0)$ to the one of interest in practice, that is, to $v=h$ given by

[TABLE]

For this initial data, we also use $u=F_{SA}(S,K,\nu,\sigma,\rho,\tau)$ to denote the solution of the $\lambda$SABR PDE, Equation (1). Let $F_{BS}$ be the Black-Scholes formula (or function), which is the solution of the Black-Scholes PDE with initial data $h=|e^{x}-K|_{+}$. Therefore,

[TABLE]

with $F_{BS}(S,K,\sigma,\tau)$ explicitly computable in terms of the cumulative normal distribution function $N$. For the initial condition $v=h$, the expansion (3) becomes

[TABLE]

One of our main results is to provide a method to explicitly compute the coefficients $F_{j}$ in this expansion. We thus obtain, in principle, closed form approximations of the solution $F_{SA}$. We carry through the calculations to provide explicit, closed-form formulas for $F_{1}$ and $F_{2}$. This leads to an explicit, closed-form approximation of $u=F_{SA}$ with error of the order $O(\nu^{3})$. See also the next subsection for a more precise and complete description of our results.

Unlike the usual iterative numerical methods, our method based on closed form solutions has a limited precision and hence a limited range of applications (small $\nu$ and not too large $\tau$). However, the type of explicit, closed form formulas that we obtained using our method are favored in financial applications, where a great precision is not needed, but it is important to have very fast, easy to implement methods. In those applications, $F_{SA}(S,K,\nu,\sigma,\rho,\tau)$ represents the no-arbitrage forward price of a European Call Option with strike $K$, volatility $\sigma$, and time to expiry $\tau$ in the SABR model. Nevertheless, a suitable modification of our method can also be used to obtain approximations of the solution with an arbitrary predetermined precision and for any initial data. This uses an approximation of the Green function of the $\lambda$SABR PDE that is of high order in time in combination with the bootstrap procedure used in [10, 14]. The bootstrap procedure is, however, iterative, it relies on numerical integration and time discretization, which increase its cost, and is no longer explicit. See Remarks 1.7 and 1.10. The numerical integration leads us to Fredholm integral operators, see also [29].

While, from a mathematical point of view, the SABR model without mean reversion is well understood since it fits into standard theories, this is not the case once one includes the mean-reverting term. For this reason, our results are less complete in the mean-reverting case. For instance, when $\kappa=0$, we do not need boundary conditions at $\sigma=0$, but when $\kappa>0$, this is not so clear. This paper was first circulated in 2014 as an SSRN Preprint and is based on the second named author’s 2012 Master Thesis at Pennsylvania State University (under the supervision of A. Mazzucato and V. Nistor) [34]. In fact, the whole project has started while all the authors were still at Pennsylvania State University, whose continued support we are glad to acknowledge.

The method, history, and the main results

Let us now explain our method, state our main results, and briefly discuss how they compare with earlier works. As explained above, we shall consider the $\lambda$SABR PDE, defined in Equation (1), and our goal is to approximate its solution $u$, in general, and the solution $u=F_{SA}(S,K,\nu,\sigma,\rho,t$ for the initial data $u(0)=h=|e^{x}-K|_{+}$, in particular. (From now on we shall use $t$ instead of $\tau$ in the $\lambda$SABR PDE.)

To set up some notation and concepts, let us consider the general evolution equation

[TABLE]

We shall formally write the solution of this equation as $u(t)=e^{tL}v$. Assume that $L=L_{0}+V$. A typical result, expressing $e^{tL}v$ in terms of $L_{0}$, $e^{sL_{0}}$, and $V$, is the following Duhamel-Dyson perturbative series expansion:

[TABLE]

where $I_{0}=e^{tL_{0}}v$, $I_{1}=\int_{0}^{t}e^{(t-\tau)L_{0}}Ve^{(t-\tau)L_{0}}v\,d\tau$, and, in general,

[TABLE]

$d\tau=d\tau_{1}\ldots d\tau_{k}$, $k<n$. A similar formula gives the error term $\mathcal{E}_{n}$, except that the last $L_{0}$ is replaced by $L$, see Equation (19). This expansion was used also in [32] to study the distribution kernel of the SABR PDE.

The terms of the expansion in Equation (8) are thus integrals of expressions of the form $e^{s_{0}L_{0}}Ve^{s_{1}L_{0}}Ve^{s_{2}L_{0}}\ldots e^{s_{n-1}L_{0}}Ve^{s_{n}L_{0}}v$, where $s_{0}+s_{1}+\ldots+s_{n}=t$, $s_{j}\geq 0$ (the exponentials and the $V$s alternate). In general, we do not know any method to explicitly compute integral of this form. However, if $L_{0}$ and $V$ generate a finite dimensional Lie algebra, we can “move” all the exponentials to one side of the formula such that the resulting integral is of a form that can easily be computed. Concretely, this is achieved through several applications of the Campbell-Hausdorff-Backer formula provided by Theorem 1.1, as in [10, 14]. See [54] for a form of this procedure that leads to solvable Lie algebras (instead of nilpotent ones). In the specific case of the $\lambda$SABR PDE, we obtain

[TABLE]

where $t\geq\tau_{1}\ldots\geq\tau_{k}\geq 0$, $V$ and the exponentials alternate, and $P_{k}(t,\tau)$ is a polynomial in $t$ and $\tau=(\tau_{1},\ldots,\tau_{k})$ with coefficients differential operators, similar results, but taking $t$ as a small parameter were obtained in [10, 11, 14, 44]. The integration over $\tau$ can then be carried out in a straightforward manner since the “exponential” $e^{tL_{0}}$ (defined using the semi-group generated by $L_{0}$) does not depend on $\tau$ anymore, and thus can be factored out of the integral. The connection between the Duhamel-Dyson perturbative series expansion and Lie algebra techniques was first noticed in [14].

Perturbative series expansions, in general, and the Duhamel-Dyson perturbative series expansion, in particular, were used before in [32] for similar purposes in conjunction to the explicit formulas for the heat kernel of the Laplacian on the hyperbolic plane to study the SABR PDE. In general, perturbative expansions were very extensively studied in physical applications. They were used in [10, 11, 14, 44], using the “time” $t$ as a small parameter, to compute in general integrals of the form $I_{k}$ up to an error of order $t^{\infty}$ (that is, to an arbitrary order in $t$). An important progress in this direction was achieved in [44], where complete explicit calculations for the SABR PDE were obtained (taking time as a small parameter). As that paper shows, this is a very tedious, albeit elementary calculation. It also undescores the need to find alternative perturbations using other “small parameters,” in particular, our calculations for the $\lambda$SABR PDE are less tedious. See also [49, 43, 41] for related calculations for the small time asymptotics. Going further back in time, one needs to mention, among many other contributions, the works of Henry-Labordere [35, 37, 36] who used Riemannian geometry heat kernel approximations, the works of Gatheral and his collaborators who used heat kernel asymptotics to study the implied volatility, and the works of Lesniewski and his collaborators [31, 32, 33], who introduced and studied the SABR model, the work of Fouque, Papanicolaou, Sircar, and Solna [24, 23, 25], who studied stochastic volatility models and singular perturbation techniques in option pricing, see also [22]. For many of these authors, the motivation was to study stochastic volatility models, which are discussed below in more detail, and whose importance is underscored also by our results.

The explicit calculation of the terms $e^{tL_{0}}P_{k}(t,\tau)h$, albeit elementary, becomes more and more tedious as $k$ growths. Because of this, we only compute explicitly the second order approximation (in $\nu$) of $F_{SA}$. That amounts to express $V$ as a multiple of $\nu$ and to collect the terms with the same powers of $\nu$ up to order two, the remaining terms being included in the error. For the coefficients of $\nu$ and $\nu^{2}$, we obtain complete, closed form (i.e. fully explicit) expression that we subsequently test numerically. (The first term, more precisely, the coefficient of $\nu^{0}$, is given by the Black-Scholes formula, so there is no additional work to be done.) Moreover, for the initial value $v=h=|e^{x}-K|_{+}$, the calculation of $e^{tL_{0}}P_{k}(t,\tau)h$ greatly simplifies, since $\partial_{\sigma}h=0$.

Here is now an example of the type of results that we obtain. Let

[TABLE]

be the usual terms appearing in the Black-Scholes formula and $N$ be the normal cumulative distribution function. Hence $N^{\prime}(x)=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}$. Then, for the term $F_{1}$ appearing in Equation (6), we obtain

[TABLE]

The formula for $F_{2}$ is significantly more complicated, as seen in Theorem 1.11.

In practice, one is often interested in the “implied volatility,” or, more precisely, in the “SABR model implied volatility.” The SABR model implied volatility $\sigma_{imp}$ is defined by $F_{SA}(S,K,\nu,\sigma,\rho,t)\,=\,F_{BS}(S,K,\sigma_{imp},t)$. (See also Equation (17) for the “Black-Scholes market implied volatility.”) The expansion (6) then yields an expansion in $\nu$ for the SABR model implied volatility, namely,

[TABLE]

The coefficients $e_{1}$ and $e_{2}$ are obtained by a tedious, but elementary calculation. For example,

[TABLE]

if $\kappa=0$. See Theorem 2.2 for the more complicated formula for $e_{2}$. We also extend our approach to approximate $\partial_{S}F_{SA}(S,K,\nu,\sigma,\rho,\tau)$, which is important in practice, see Theorem 2.1.

We thus obtain two second order approximations of $F_{SA}(S,K,\nu,\sigma,\rho,t)$, namely,

[TABLE]

obtained by truncating Equation (6) and

[TABLE]

obtained by truncating the formula for $\sigma_{imp}$ to second order. To numerically test our results, we compare these two approximations with approximations obtained via other methods. Thus, we first compare our two approximation formulas with equation (2.17) in [31], which yields a well known approximation of $F_{SA}$, often used in practice. Our approximation agrees almost exactly with Hagan’s approximation when $t$ or $\nu$ are small. In fact, the power series approximation in $y=\ln(F/K)$ of the first order coefficient of $\nu$ in Hagan’s formula coincides with our $e_{1}$ term. When $t$ is taken one year, both approximation are very accurate, the largest difference between the simulated price is 0.8 cents, or around 0.1–0.2% of the option price, which is well within the bid-ask spread. As $t$ increases to 10 years, the largest error increases to 1.5%. However, note that the option also is more expensive for a longer time to maturity, so the percentage (or relative) error is still reasonable. We have also compare our two approximate formulas with Hagan’s formula and with Monte-Carlo and Finite Difference approximations of the option price across a range of strikes and maturities. Again, the results are good, at least for $t\leq 1$. See Sections 3 and 4 for details. We note that in these numerical tests, it is a challenge to perform an accurate enough test using either Monte Carlo or Finite Differences, which shows the importance of having alternative, faster methods. It is possible that a modification of our method, combined with the results of [10], will allow us to extend the range of $\beta$.

Practical motivation

The methods and results in this paper can be properly understood only when put into the perspective of their applications. Let us say a few words about this, without entering into unnecessary details.

Options and other derivatives have long been used in financial applications. They were mathematically rigorously priced for the first time by the famous Black-Scholes model [9], under the assumption that the underlying asset (for example, a stock) follows a log-normal distribution, or geometric Brownian motion. While the use of the Black–Scholes model is widespread in practice, the Black-Scholes model is known to produce option values significantly different from the ones in the market. The reasons for these differences are studied, for instance, in Rubinstein’s seminal paper [51]. For us, the most important reason is that the volatility of the underlying asset is non constant.

It is therefore increasingly common these days to consider models that relax an assumption of the constant volatility of the underlying asset. In practice, the time-varying volatility has lead to the concept of Black-Scholes market implied volatility (sometimes, simply, implied volatility) $\sigma_{imp}$, defined as the volatility that would have to be used in the Black-Scholes formula to recover (at a given time) the market price. That is,

[TABLE]

where $C_{M}(S,K,\tau)$ is given either by market prices or by another model (such as the SABR model in our case). The forward price $F$ is related to the actual price by the formula $F=e^{rt}C$, where from now on, $t$ denotes the time to expiry. If the geometric Brownian motion provided a perfect description of the behavior of the underlying, then the Black-Scholes implied volatility would be constant as a function of time to expiry $\tau$ and strike $K$. However, the practice shows immediately that this is not the case. See [26, 28, 52] for more information.

The fact that the volatility is not a constant in practice has motivated the introduction of several other models. Among them, we are interested in the stochastic volatility models (see [18, 31, 38, 40], for example). In stochastic volatility models, the volatility is not only non-constant in time, but is in fact a random variable with its own volatility, denoted here $\nu$ and called volatility of volatility or vol-of-vol. This explains why the SABR PDE has two space variables ($x$ and $\sigma$ in our convention), unlike the Black-Scholes PDE, which has only one space variable ($x$ in our convention). One of the disadvantages of stochastic volatility models, however, is their greater computational cost. It is important therefore to find good closed form approximations of the solutions to option pricing in stochastic volatility models. It is the purpose of this paper to provide such an approximation, for small values of the vol-of-vol parameter $\nu$, by a Taylor-type expansion in $\nu$. See [22, 44] for more on stochastic volatility models and for further references to this topic.

We note that our results provide realistic estimates of the volatility of volatility parameter $\nu$. They also show that that using stochastic volatility improves the fit to the data, and, by further including the mean reverting term, we get yet even better results. See Subsection 3.5 and, especially, Table 3.

Contents of the paper

The paper is organized as follows. Section 1 contains the main results of the paper. In that section, we explain our method in general and then we perform in detail the calculations for the particular case of the $\lambda$SABR PDE, thus obtaining a third order accurate in $\nu$ expansion of the solution. That amounts to finding the exact formulas for $F_{1}$ and $F_{2}$ in the expansion $F_{SA}=F_{BS}+\nu F_{1}+\nu^{2}F_{2}+\ldots$. In this section, we also provide an approximation for the distribution kernel of the $\lambda$SABR PDE, which allows, in principle, to approximate its solution to an arbitrary predetermined precision for any initial data. One of the advantages of our method is that it allows us to include the mean reverting term, so we keep it in our formulas for most of our calculations. In Section 2, we use extend our results by providing an asymptotic expansion for the implied volatility $\sigma_{imp}=\sigma+\nu e_{1}+\nu^{2}e_{2}+\ldots$ obtained from the equation $F_{SA}=F_{BS}(\sigma_{imp})$. We obtain a similar expansion for the $\Delta:=\partial_{S}C$ hedging parameter. In Section 3 we use our two approximations of the price, namely $F_{SA,2}:=F_{BS}+\nu F_{1}+\nu^{2}F_{2}$ and for $F_{D}:=e^{rt}C_{BS}(\sigma+\nu e1+\nu^{2}e_{2})$. We use these formulas, as well as Hagan’s formula to calibrate our model using market data and to compare these models. Our approximations are seen to be quite competitive in this regard. The results of the market calibration also provide us with a realistic range of the parameters for the numerical tests that we perform in Section 4. In that section, we perform several tests to see how well our results approximate the exact solution of the $\lambda$SABR PDE. In particular, we also compare our solutions to the one obtained using a Finite Difference (FD) test. We discuss in detail the results of the FD implementation and its challenges. In the last section, we outline some possible extensions of our methods and results: a more precise method to deal with the mean reverting factor and a method to treat $\beta\neq 1$ in the original SABR PDE.

We thank Wen Cheng, Anna Mazzucato, Dan Pirjol, Camelia Pop, and Siyan Zhang for useful discussions.

1. The Duhamel-Dyson perturbative series expansion

In this section, we introduce and explain in detail our method. To illustrate it, we perform some general calculations leading to the second order approximation of the Green function of the $\lambda$SABR PDE with initial condition $h$. This then leads to a determination of the coefficients $F_{1}$ and $F_{2}$ in Equation (6), which is our main approximation formula of the solution $F_{SA}=F_{BS}+\nu F_{1}+\nu^{2}F_{2}+\ldots$ of the $\lambda$SABR PDE.

As we are interested in this paper mainly in numerical methods for the $\lambda$SABR PDE, we keep the theoretical aspects (including proofs), to a minimum, trying to obtain as quickly as possible our formulas, which we then test numerically.

1.1. General results: The Duhamel and Campbell-Hausdorff-Backer formulas

As in the Introduction, we cast our study of the SABR PDE as a particular case of the problem of approximating general evolution equations (7) (i.e. of the form $\partial_{t}u-Lu=0$, $u(0)=v$). In this subsection, we recall some general results pertaining to the Equation (7), including Duhamel’s formula and its generalization, the Duhamel-Dyson perturbative series expansion, as well as the Campbell-Hausdorff-Backer formula. Beginning with the next subsection, we specialize to the case of the $\lambda$SABR PDE.

The operator $L$ appearing in Equation (7) is called the parabolic generator (of the given PDE or of the semi-group $e^{tL}$). Whenever the solution of the evolution problem (7) exists, is unique, and depends continuously on the initial data (in a suitable functions space $V$), we shall say that our problem is well posed (in Hadamard’s sense). In that case, we shall write $u(t)=e^{tL}h$, where $e^{tL}:V\to V$ is a continuous, linear map such that $e^{tL}v$ depends continuously on $v\in V$ (i.e. $e^{tL}$ is a $c_{0}$-semi-group [2, 50, 53]). We shall often use this notation in a formal way, that is, even if it was not fully justified mathematically.

Assume that $L=L_{0}+V$ in our evolution equation, Equation (7), and that $u(0)=v$, as before. Assume that both $L_{0}$ and $L$ generate $c_{0}$-semi-groups. We then use $u(t)=e^{tL}v$ and $\partial_{t}u-L_{0}u=Vu$ to write Duhamel’s formula in the form

[TABLE]

By substituting Duhamel’s formula for $u(\tau)$, that is, for $t$ replaced by $\tau$, back in the last integral of the original Duhamel’s formula and then iterating this procedure $(n-1)$-times, we obtain Equations (8) and (9) of the Introduction. The error term $\mathcal{E}_{n}$ in Equation (8) is the last integral of the resulting formula and is similarly given by

[TABLE]

Notice that in the last exponential we have $L$ instead of $L_{0}$. If we were to iterate one more time to obtain the Duhamel-Dyson formula for $n+1$, we would substitute Duhamels’s formula (18) applied to $e^{\tau_{n}L}v$ in the integral defining $\mathcal{E}_{n}$.

In general, integrals of the kind defining $I_{k}$ are notoriously hard to compute, and a lot of work has been devoted to understanding them. However, if $L_{0}$ and $V$ generate a finite dimensional Lie algebra, one can compute these integrals by establishing a “Campbell-Hausdorff-Backer” formula (CHB) [10, 14, 54]. The idea is to use the CHB formula to shift all the exponential operators $e^{tL_{0}}$ to the left in our cases of interest. This is what we are going to do in this paper as well. We note that in our paper, the CHB formula reduces to a finite sum, so there are no convergence issues. It is possible to use the CHB formula also in certain situations when the sum is not finite; see [54] for an example.

Recall that $[T_{1},T_{2}]:=T_{1}T_{2}-T_{2}T_{1}$ denotes the commutator of $T_{1}$ and $T_{2}$ and $ad_{T}(T_{1})=[T,T_{1}]$. Clearly, We are ready to state now one of our main technical ingredients, the Campbell-Hausdorff-Backer (CHB) formula (see [10, 14, 54] for a proof in a greater generality).

Theorem 1.1**.**

Let $L_{0}$ and $V$ be operators such that $\operatorname{ad}_{L_{0}}^{m+1}(V)=0$ for some $m$. For any $\theta>0$, define

[TABLE]

Then

[TABLE]

and, similarly,

[TABLE]

Let $\mathbb{R}[\sigma,\partial_{\sigma},\partial_{x}]$ be the algebra of differential operators in $\partial_{\sigma}$ and $\partial_{x}$ with coefficients polynomials in $\sigma$. The CHB formula allows us to compute the terms $I_{k}$ of Equation (9) in the Duhamel-Dyson perturbative series expansion, Equation (8). Indeed, a direct calculation gives the following corollary.

Corollary 1.2**.**

With the assumptions of Theorem 1.1, we have

[TABLE]

and hence $I_{k}=e^{tL_{0}}p_{k}(t)$, where $p_{k}(t)\in\mathbb{R}[\sigma,\partial_{\sigma},\partial_{x}]$ depends polynomialy on $t$, $L_{0}$, and $V$.

Similarly, by moving the exponentials to the right, we obtain that $I_{k}=q_{k}(t)e^{tL_{0}}$, where $q_{k}(t)\in\mathbb{R}[\sigma,\partial_{\sigma},\partial_{x}]$ depends polynomialy on $t$, $L_{0}$, and $V$.

1.2. Commutator calculations for SABR

We would like to use the CHB formula (Theorem 1.1) and its Corollary 1.2 to identify explicitly the terms $I_{k}$ in the Duhamel-Dyson perturbative series expansion for the $\lambda$SABR PDE, namely in Equation (8). To this end, we need to introduce the decomposition $L=L_{0}+V$ and compute the relevant commutators in order to see that we are in position to use the CHB formula.

Let us introduce the differential operators:

[TABLE]

In this way, the $\lambda$SABR PDE (Equation (1)) becomes a particular case of our general evolution equation, Equation (7), for $L\,:=\,L_{0}+\nu L_{1}+\nu^{2}L_{2}$. We can then use the general results of the previous subsection, for $L=L_{0}+V$, where $V:=\nu L_{1}+\nu^{2}L_{2}$. The notation for the differential operators $L_{0}$, $L_{1}$, $L_{2}$, $L$, and $V$ will remain fixed throughout the paper.

The following lemma shows that, in the case of the $\lambda$SABR PDE, we are in position to use the CHB formula an hence in position to compute explicitly the terms in Duhamel-Dyson perturbative series expansion for $L_{0}$ and $V$:

Lemma 1.3**.**

Let $b:=\partial_{x}(\partial_{x}-1)$. Then $\operatorname{ad}_{L_{0}}(L_{1})=-\sigma[\rho\sigma^{2}\partial_{x}+\kappa_{0}(\theta-\sigma)]b$ and $\operatorname{ad}^{j}_{L_{0}}(L_{1})=0$, if $j>1$. Similarly, $\operatorname{ad}_{L_{0}}(L_{2})=-\frac{1}{2}(\sigma^{2}+2\sigma^{3}\partial_{\sigma})b$,  $\operatorname{ad}^{2}_{L_{0}}(L_{2})=\sigma^{4}b^{2},$ and $\operatorname{ad}^{j}_{L_{0}}(L_{2})=0$ if $j>2$.

In particular, $\operatorname{ad}^{j}_{L_{0}}(V)=0$ if $j>2$, and hence we can use the CHB formula.

1.3. Expansion in the powers of $\nu$

Let us now expand the terms $I_{k}$ in the Duhamel-Dyson perturbative series expansion in powers of $\nu$, see Equation (8). Then we use the CHB formula of Theorem 1.1 and its corollary, Corollary 1.2 (whose use was justified by Lemma 1.3), to obtain

Proposition 1.4**.**

We have

[TABLE]

where $B_{j}=P_{j}e^{tL_{0}}=e^{tL_{0}}Q_{j}$, $j<n$, with $P_{j},Q_{j}\in\mathbb{R}[\sigma,\partial_{\sigma},\partial_{x}]$it. Moreover, $\mathcal{R}_{n}=Pe^{tL_{0}}+\mathcal{E}_{n}=e^{tL_{0}}Q+\mathcal{E}_{n}$, with $\mathcal{E}_{n}$ as in Equation (19).

We now turn to the calculation of the terms $B_{j}$, $j=1,2$. The goal is to ultimately obtain exact formulas for the terms $F_{1}=B_{1}h$ and $F_{2}=B_{2}h$ in the expansion $F_{SA}=F_{BS}+\nu F_{1}+\nu^{2}F_{2}+\ldots$ of Equation (6).

In view of Equations (8) and (9) and of the relation $V=\nu L_{1}+\nu L_{2}$, let us introduce the notation

[TABLE]

The integrals above are defined as in [14, 54]. Also, let us consider also the “error term”

[TABLE]

(We notice again that the last exponential in the last formula is different from the other ones, with $L$ being used instead of $L_{0}$.) By collecting the terms with like powers of $\nu$ in the Duhamel-Dyson perturbative series expansion for $L=L_{0}+V$, where $L_{0}$, $V$ and all the other differential operators are as introduced in Equation (20), we see that the expansion of Propositon 1.4 becomes

[TABLE]

That is,

[TABLE]

Note that Equation (23) is an exact formula – not just an approximation – whenever the exponentials (or semi-groups) are defined. Also, it is an equality of operators and hence

[TABLE]

which we shall mainly use for $v(x)=h(x):=|e^{x}-K|_{+}$.

1.4. Symbolic formulas for $B_{j}$, $j=1,2$.

We now turn to a symbolic calculation of the terms $B_{j}$, $j=1,2$, of Proposition 1.4. In view of Equation (24), in order to compute $B_{1}$ and $B_{2}$, we need to compute $J_{1}$, $J_{2}$, and $I(L_{1},L_{1})$. To this end, we use the CHB formula to successively move the exponentials to the left in the integrals defining $J_{1}$, $J_{2}$, and $I(L_{1},L_{1})$. First, the definitions of the integral $J_{1}$ in Equation (22) gives

[TABLE]

Similarly,

[TABLE]

Finally, the definition of the integral $I$ in Equation (22) gives

[TABLE]

This completes the symbolic calculations of $B_{1}$ and $B_{2}$, that is, in terms of the exponentials $e^{tL_{0}}$ and the commutators of $L_{j}$.

Remark 1.5**.**

Note that we obtained formulas for which all the exponentials of the form $e^{tL_{0}}$ have been shifted to the right and formulas for which they have all been shifted to the left. The two forms are both needed, but they serve different purposes. The formulas that we have obtained are explicit enough to allow for exact calculations. In particular, the coefficients $B_{1}=J_{1}$ and $B_{2}=J_{2}+I(L_{1},L_{1})$ of Equation (23) can be written in the form $B_{j}=P_{j}e^{tL_{0}}$ and $B_{j}=e^{tL_{0}}Q_{j}$, with $P_{j},Q_{j}\in\mathbb{R}[\sigma,\partial_{\sigma},\partial_{x}]$, $j=1,2$, a fact that will be exploited shortly. The same is true for the other terms in the expansion ($j\geq 3$).

Remark 1.6**.**

Let us notice that $e^{tL_{0}}$ is defined (since it corresponds to a family of heat equations). In particular,

[TABLE]

For $e^{tL}$, the situation is more complicated. Let $M:=\mathbb{R}\times[0,\infty)\ni(x,\sigma)$. In case $\kappa_{0}=0$, that is, if there is no mean reverting term, then $L$ is generated by the derivatives $\sigma\partial_{x}$ and $\sigma\partial_{\sigma}$ together with smooth, totally bounded functions (or coefficients) in $\sigma$ and $x$. This type of differential operators were studied in [5, 6] in the framework of differential operators on “Lie manifolds.” In our case, $M$ is a Lie manifold for with the structure generated by the standard compactification of $M$ to a disk using hyperbolic coordinates. The geometry is that of the hyperbolic plane, something noticed also in [32]. Since Lie manifolds are of bounded geometry [5], the theory developed in [46] and, more recently by Herbert Amann [3, 4] shows that the equation $\eqref{eq.IVgen}$ is well-posed in the Sobolev spaces associated to this geometry, that is, in spaces of the form

[TABLE]

This well-posedness result is, unfortunately, not known in the mean reverting case. See however [33, 54]. Because of this, our treatment of the mean reverting case will be shorter and somewhat formal. Analysis on this type of manifolds is related to that on hyperbolic spaces [32]. It is also related to the analysis of edge singularities for PDEs on polyhedral domains [7, 13, 16, 20, 42].

Remark 1.7**.**

We can, in principle, compute exactly the distribution kernels (or Green functions) of the operators $B_{j}$. They are very closely related to the kernels of the operators $e^{tL_{0}}$. We will do that for $B_{1}$, see Remark 1.10. These kernels can then be used to integrate against any initial data $v$ to obtain, by truncating in the expansion in the small parameter, a semi-discretization (i.e. a discretization only in time) of our PDE. This semi-discretization can be turned into a full discretization via an approximation of $v$ by functions in a finite dimensional space and then by using numerical integration. This approximation is very good for small time. It can be turned out into an arbitrary precision approximation for any time using the bootstrap method in [10, 14]. See Remark 1.10 for more details.

1.5. Heat kernels and convolutions

The Green functions of the operators $B_{j}$ discussed in Remark 1.7, while explicit in principle, are given by pretty complicated formulas. The calculation of $B_{j}h$, however, can further be simplified using the specific properties of $h$ and $e^{tL_{0}}$. Let us define

[TABLE]

Then it is well known that

[TABLE]

in the sense that $e^{tL_{0}}f(x)=\int_{\mathbb{R}}\,e^{tL_{0}}(x,y)f(y)dy$ and $u(t,x,\sigma):=e^{tL_{0}}f(x)$ satisfies the partial differential equation $\partial_{t}u(t)-L_{0}u(t)=0$ with initial condition $u(0)=f$. (That is, $e^{tL_{0}}(x,y)$ is the Green function of $\partial_{t}-L_{0}$.) In particular, this gives $e^{tL_{0}}h$, which will be needed in the final formula, as follows. Let $N(x):=(2\pi)^{-1/2}\int_{-\infty}^{x}e^{-t^{2}/2}dt$ be the normal cumulative function. Then

[TABLE]

where we used the notation of Equation (11).

This suggests to consider integral kernel operators $T_{k}$, where $k(x,y)$ is a suitable measurable function and $T_{k}f(x):=\int_{\mathbb{R}}k(x,y)f(y)dy$. Then $k$ is called the distribution kernel of $T_{k}$. If $k(x,y)=\phi(x-y)$, then we shall write also $C_{\phi}=T_{k}$ and we shall refer to $C_{\phi}$ as the operator of convolution with $\phi$. So

[TABLE]

In particular, $e^{tL_{0}}(x,y)$ is the operator of convolution with $\phi_{t}$:

[TABLE]

All our convolution operators will be convolutions in the $x$ variable, but they will often depend on $\sigma$ as a parameter. The parameter $\sigma$ in the notation of the convolution will sometimes be omitted, as in the last equation.

Let also $h(x)=|e^{x}-K|_{+}=(e^{x}-K)^{+}$ be the usual pay-off of a European call option, Equation (4). Recall that a function $\phi:\mathbb{R}\to\mathbb{C}$ is said to be of Schwartz class if $x^{k}\partial_{x}^{j}\phi(x)$ is bounded for all $k,j\geq 0$. If $\phi$ is of Schwartz class, then $\phi$ and all its derivatives are integrable. We have the following lemma, which underscores the important role played by the polynomial $b=\partial_{x}^{2}-\partial_{x}$.

Lemma 1.8**.**

Assume $\phi$ is function of Schwartz class. Then

[TABLE]

in the sense that $(C_{\phi}\psi)^{\prime}=C_{\phi}(\psi^{\prime})=C_{\phi^{\prime}}(\psi)$, for $\psi$ of Schwartz class. If, moreover, there exist two constants $C>0$ and $\epsilon>1$ such that $|\phi^{(k)}(x)|\leq Ce^{-\epsilon|x|}$ for $k=0,1,2$, then

[TABLE]

In particular, if $P\in\mathbb{R}[\sigma,\partial_{x}]$, the operators of the form $Pe^{tL_{0}}$ are convolution operators in the $x$ variable.

Note that $Pe^{tL_{0}}=e^{tL_{0}}P$ if $P\in\mathbb{R}[\sigma,\partial_{x}]$. More precisely, we obtain that $Pe^{tL_{0}}=e^{tL_{0}}P$ is a family parametrized by $\sigma$ of convolution operators in the $x$ variable.

Proof.

The proof is by direct calculation. For the convenience of the reader, we note that the second relation follows from the first relation and from $bh=(\partial_{x}^{2}-\partial_{x})h=K\delta_{\ln K}$, where $\delta_{y}$ is the Dirac distribution concentrated at $y$. ∎

The following remark explains the idea of the last step of the calculation yielding the terms $F_{1}=B_{1}h$ and $F_{2}=B_{2}h$ in the expansion $F_{SA}=F_{BS}+\nu F_{1}+\nu^{2}F_{2}+\ldots$ of Equation (6).

Remark 1.9**.**

Let us first notice that

[TABLE]

for some polynomials $\tilde{H}_{n}(x)$, $g\geq 0$, whose coefficients are functions of $\sigma$. Let $H_{n}$ be the Hermite polynomials and

[TABLE]

Then $\tilde{H}_{n}=(-\frac{1}{\sigma\sqrt{t}})^{n}H_{n}(\ell(x))$. See Remark 1.14 for more details. Let $B_{j}=e^{tL_{0}}Q_{j}$, as in Proposition 1.4, and let $bf=\partial_{x}^{2}-\partial_{x}$, as before. Then $Q_{j}=\sum_{i}G_{ji}\partial_{\sigma}^{i}$ (a finite sum), with $G_{ji}\in\mathbb{R}[\sigma,\partial_{x}]$. We will see that $G_{j0}=\tilde{G}_{j}b$, with $\tilde{G}_{j}=\sum_{i}a_{ji}\partial_{x}^{i}$, where $a_{ji}=a_{ji}(\sigma)$ are polynomials in $\sigma$. Recall also that $\partial_{x}$, $\sigma$, and $e^{tL_{0}}$ commute, and hence $e^{tL_{0}}$ and $G_{j0}$ commute. Lemma 1.8 then simplifies the calculation of $B_{j}h$, $j=1,2$, as follows:

[TABLE]

In order to complete our calculation, all that is thus left to do is to find the differential operators $\tilde{G}_{j}$, $j=1,2$, or, equivalently, the polynomials $a_{ji}(\sigma)$. We also notice that $\partial_{x}(\psi(x-a))=(\partial\psi)(x-a)$, so there is no danger of confusion when writing $\tilde{G}_{j}\phi_{t}(x-\ln K,\sigma)$.

The last remark explains also why, if we are interested in evaluating terms of the form $B_{j}h$, it is more convenient to keep the exponential to the left of the formula (to get rid of the $\partial_{\sigma}$).

This discussion allows us to identify the distribution kernel of $B_{1}$.

Remark 1.10**.**

Recall that B_{1}=J_{1}=e^{tL_{0}}\big{(}tL_{1}-\frac{t^{2}}{2}\operatorname{ad}_{L_{0}}(L_{1})\big{)}, which combines the formula for $J_{1}$, Equation (22) with the identification (24). Then

[TABLE]

This gives $B_{1}v=C_{\psi_{1}}\partial_{\sigma}v+C_{\psi_{2}}v$. Then the solution of any initial value problem with $u(0)=v$ can be approximated by $u\approx D_{t}v$,

[TABLE]

Of course, a better estimate will be obtained if one includes also the $B_{2}$ term, which is, however, much more complicated. One can perform a bootstrap, that is, divide the interval $[0,T]$ into $N$ subintervals, and solve on each interval the corresponding initial problem approximately. That is

[TABLE]

This method works well once one establishes an error $\|u(t)-D_{t}v\|=O(t^{a})$ with $a>1$. See [10, 14].

1.6. Calculation of the differential operators $\tilde{G}_{j}$

We now pursue in detail the method outlined in Remark 1.9 to compute explicitly closed form expressions for $F_{1}=B_{1}h$ and $F_{2}=B_{2}h$, where, we recall, $B_{1}=J_{1}$ and $B_{2}=J_{2}+I(L_{1},L_{1})$. The differential operators $\tilde{G}_{j}\in\mathbb{R}[\sigma,\partial_{x}]$, $j=1,2$, considered in this subsection are as defined in that remark and thus satisfy $B_{j}h(x)=K\tilde{G}_{j}\phi_{t}(x-\ln K,\sigma)$.

We continue to use the notation $b=\partial_{x}^{2}-\partial_{x}$, for simplicity. Recall also from Lemma 1.3 that $\operatorname{ad}_{L_{0}}(L_{1})=-\sigma[\rho\sigma^{2}\partial_{x}+\kappa_{0}(\theta-\sigma)]b$. Lemma 1.8 and Equation 26 then give (since $t=0$)

[TABLE]

Hence $\tilde{G}_{1}=:\frac{t^{2}\sigma}{2}[\rho\sigma^{2}\partial_{x}+\kappa_{0}(\theta-\sigma)]=a_{10}+a_{11}\partial_{x}$ and

[TABLE]

which yields right away the formula (12) in the Introduction.

The second term is similar, but the calculation requires more work. We shall use that $L_{1}h=L_{2}h=0$. We first calculate the intermediate term

[TABLE]

We then use this calculation together with $[L_{0},L_{1}]^{2}h\,=\,\sigma^{2}[\rho\sigma^{2}\partial_{x}+\kappa_{0}(\theta-\sigma)]^{2}b^{2}h$, the formulas (27) and (28), the relation $\partial_{\sigma}h=0$, and the fact that $\partial_{\sigma}$ commutes with $\partial_{x}$, to obtain

[TABLE]

Let

[TABLE]

Then $\tilde{G}_{2}=\sum a_{2i}\partial_{x}^{i}$ and, finally,

[TABLE]

Putting together our calculations, we obtain the following result. (Recall the notation of Equation (6), but see also Equation (15).)

Theorem 1.11**.**

The formal second order approximation of solution $F$ of the $\lambda$SABR PDE is

[TABLE]

where $e^{tL_{0}}h$ is given by Equation (31), $B_{1}h=F_{1}$ and $B_{2}h=F_{2}$ are given explicitly by Equations (38) and (40), and $x=\ln(Se^{rt})$, as before.

Remark 1.12**.**

Since the initial data $h(x)=|e^{x}-K|_{+}=(e^{x}-K)^{+}$ is very smooth in $\sigma$ – in fact, even independent of $\sigma$ in our stochastic volatility models – the methods of [14] will give that $\mathcal{R}h$ is bounded (and hence that $\nu^{3}\mathcal{R}h$ is very small for $\nu$ small). That is,

[TABLE]

A rigorous proof of this fact is beyond the scope of this paper, however.

Similar formulas for $\kappa=0$ but for general $\beta$ were obtained in [44]. They use a small time asymptotic method, very similar to the one in [14], which leads to much longer formulas.

1.7. Remarks on the implementation

When implementing the formula for $F_{SA,2}(S,K,\nu,\sigma,\rho,t)$, it will be convenient to take into account the following remarks. For simplicity, we restrict to the case $\kappa_{0}=0$ from now on.

Remark 1.13**.**

First of all, for our calculations, it was convenient to work with the forward prices, usually denoted $F$ and decorated with various indices. In practice, however, one may need to use the actual prices, denoted $C$ and decorated with the corresponding indices. They are related by the formula $F=e^{rt}C$, which corresponds to the fact that the price $C$ is given in “today’s” currency, whereas the forward price is quoted in using the value of the currency at the expiration. Here, of course, $r$ is the interest rate and $t$ is the time to expiration, as before. This is in agreement with our calculations. Indeed, let $d_{\pm}:=\frac{\ln S-\ln K+rt}{\sigma\sqrt{t}}\pm\frac{\sigma\sqrt{t}}{2}:=\frac{\ln(F/K)}{\sigma\sqrt{t}}\pm\frac{\sigma\sqrt{t}}{2}$, as before, Equation (11). In particular, $\ell(x-\ln K)=d_{-}$. We recall then that the Black-Scholes formula for the price $C_{BS}(S,K,\sigma,t)$ of a call option with strike $K$, underlying $S$, and volatility $\sigma$ is given by

[TABLE]

Recall that $F_{BS}=e^{tL_{0}}h$, and hence the formula $F_{BS}=e^{rt}C_{BS}$ is verified, by Equation (31), since $x-\ln K=\ln(Se^{rt}/K)$.

The term $\ln(Se^{rt}/K)$ is called the log-moneyness.

Remark 1.14**.**

Recall that the Hermite polynomials $H_{n}(x)$ (the probabilist’s version), are given by

[TABLE]

Of these, only $H_{0}(x)=1$, $H_{1}(x)=x$, $H_{2}(x)=x^{2}-1$, $H_{3}(x)=x^{3}-3x$, and $H_{4}(x)=x^{4}-6x^{2}+3$ are needed for the evaluation of $B_{1}$ and $B_{2}$. If $\ell(x)$ is any linear function of $x$ (so $\ell^{\prime}$ is a constant), then

[TABLE]

Let us fix $\ell(x):=\frac{x}{\sigma\sqrt{t}}-\frac{\sigma\sqrt{t}}{2}$ from now on, as in (36), which gives $\phi_{t}(x,\sigma)=ce^{-\ell(x)^{2}/2}$, $c\in\mathbb{R}$, and hence we have $\tilde{H}_{n}(x):=(-\ell^{\prime})^{n}H_{n}(\ell(x))$. Recall that $\phi_{t}$ is given by Equation (29) (but see also Equation (33)) and the polynomials $\tilde{H}_{n}$ are as defined in Equation (35) (that is, $\partial_{x}^{n}\phi_{t}(x)=\tilde{H}_{n}(x)\,\phi_{t}(x)$). Furthermore, in implementation, it will be useful to notice that $\ell(x-\ln K)=d_{-}$, that \tilde{H}_{n}(x-\ln K)=\big{(}-\frac{1}{\sigma\sqrt{t}}\big{)}^{n}H_{n}(d_{-}), and that $\phi_{t}(x-lnK,\sigma)=\frac{1}{\sigma\sqrt{2\pi t}}e^{-d_{-}^{2}/2}.$ Recall that $d_{-}$ was defined in Equation (11).

Remark 1.15**.**

Let $f(S,K)$ be an arbitrary function. We shall say that it is homogeneous of degree one in $(S,K)$ if $f(\lambda S,\lambda K)=\lambda f(S,K)$. If that is the case, the function $g=K^{-1}f$ can be written solely in terms of $y=\ln(Se^{rt}/K)$ (we assume $r$ and $t$ to be parameters in this discussion). For instance, the Black-Scholes pricing formula is homogeneous of degree one in $(S,K)$. The form of the Equation (1) tells us that the forward price $F_{SA}$ of a European call in the SABR model is also homogeneous of degree one in $(S,K)$. Since $C_{SA}:=e^{-rt}F_{SA}$, we have that $C_{SA}$ is also homogeneous of degree one in $(S,K)$. We also see that $B_{1}h$ and $B_{2}h$ are homogeneous of degree one in $(S,K)$ as well, which is consistent with the properties of $F_{SA}$. In particular, $F_{SA,2}$ and $C_{SA}:=e^{-rt}F_{SA}$ are also homogeneous in $(S,K)$ of degree one.

Remark 1.16**.**

Let $y:=x-\ln K$ denote the log-moneyness. We have already noticed that $\ell(y)=d_{-}$. Let

[TABLE]

so that $C_{rel}$ depends directly only on $y$ (and only indirectly on $S$ or $K$, as explained in Remark 1.15). This will be useful in implementation, where often we will be given $y$, but not $S$ and $K$, and we will want to estimate relative errors of the form

[TABLE]

with $y_{j}=ln(S_{j}e^{rt}/K)$. Also for the purpose of implementation, let us notice that we have

[TABLE]

where

[TABLE]

which follows from Equations (38) and (40) and Theorem 1.11.

2. Applications: approximations of the derivatives

and of the implied volatility

We now present two applications of the methods that we have developed.

2.1. Approximation of derivatives

Our method is especially well suited for approximating the derivatives of $C_{SA}$ (some of these derivatives are called “Greeks”). Indeed, the derivatives with respect to parameters (other than $\nu$) of the solution $C_{SA}$ of Equation (1) and Taylor expansions in $\nu$ are seen to commute. We shall look in what follows at $\partial_{S}C_{SA}$. The asymptotic expansion in $\nu$ for $\partial_{S}C_{SA}$ thus can be obtained from the corresponding expansion of $C_{SA}$ by taking the derivative term by term with respect to $S$ in Equation (6). Let us, for example, compute the approximation $\Delta_{2}:=\partial_{S}C_{SA,2}(S,K,\nu,\sigma,t)$. This approximation can be obtained directly by differentiating in Equation (41) with respect to $S=e^{x-rt}$ and $\sigma$. When differentiating with respect to $S$, we also take into account the fact that $\partial_{S}=S^{-1}\partial_{x}$, $\partial_{S}^{2}=S^{-2}(\partial_{x}^{2}-\partial_{x})$, and hence that $\partial_{S}$ and $\partial_{S}^{2}$ commute with all the other differential operators (this is because all the other differential operators are polynomials in $\sigma$, $\partial_{\sigma}$, and $\partial_{x}$). We continue to assume that $\kappa_{0}=0$.

For instance, Equation (38) gives

[TABLE]

Similarly,

[TABLE]

and hence

[TABLE]

Using also that $\partial_{S}C_{BS}=N(d_{+})$, we then obtain the following theorem, which is an analog for the “hedging parameter” $\Delta:=\partial_{S}C_{SA}$ of Theorem 1.11.

Theorem 2.1**.**

Let $t$ be the time to expiry and $e^{x}=Se^{rt}$. Then the formal second order approximation of the $\partial_{S}C_{SA}(S,K,\nu,\sigma,t)$, denoted $\partial_{S}C_{SA,2}(S,K,\nu,\sigma,t)$, is given by

[TABLE]

where $\partial_{x}B_{1}h$ and $\partial_{x}B_{2}h$ are given explicitly by Equations (49) and (51).

In other words, we can just compute the derivatives in the Duhamel-Dyson series term by term.

2.2. Implied Volatility

To gain a qualitative insight into the behavior of the implied volatility in our approximation of the SABR model, it is necessary to convert (41) into an asymptotic expansion of the implied volatility. We thus assume that the second order approximation for the implied volatility is of the form

[TABLE]

with the coefficients $e_{1}$ and $e_{2}$ to be calculated. We then let $P(\sigma):=C_{BS}(S,K,\sigma,t)$ and define $\sigma_{imp}$ by the equation

[TABLE]

We then substitute the assumed form of the $\sigma_{imp}$ function into the Black Scholes formula $C_{BS}$ and Taylor expand in $\nu$ this composite function around the initial volatility $\sigma$. See also [27, 30, 49, 44]. That is, we have

[TABLE]

Recall that $F=e^{x}=Se^{rt}$ and $y=\ln(F/K)=x-\ln K$ is the log-moneyness. Let $z=d_{-}=\ell(y)$ and $v=\sigma\sqrt{t}$. We shall need also the first two derivatives of $P:=C_{BS}$ with respect to $\sigma$, so we record them here (note that $FN^{\prime}(d_{+})=KN^{\prime}(d_{-})$)

[TABLE]

which can be made more explicit by substituting $N^{\prime}(d_{-})=\frac{1}{\sqrt{2\pi}}e^{-d_{-}^{2}/2}$. Comparing Equations (41) and (55) and matching the $\nu$-coefficients according to Equation (54), we obtain

[TABLE]

Similarly, let us denote

[TABLE]

the long factor appearing in the formula for $B_{2}h$, Equation (40). Then, using $\partial_{\sigma}P(\sigma)=Ke^{-rt}\sqrt{t}N^{\prime}(d_{-})$ and $\sqrt{2\pi}N^{\prime}(d_{-})=e^{-\frac{(y-\sigma^{2}t/2)^{2}}{2\sigma^{2}t}}$, we obtain

[TABLE]

We have the following result:

Theorem 2.2**.**

Let $t$ be the time to expiry and $y=\ln(F/K)=\ln(Se^{rt}/K)$. Then the implied volatility $\sigma_{imp}$ has an asymptotic expansion of the form

[TABLE]

The coefficients $e_{1}$ and $e_{2}$ are given by $e_{1}=-\rho\sigma\sqrt{t}d_{-}/2$ and

[TABLE]

3. Model calibration and market tests

Before we discuss our numerical tests, it is useful to see how our formulas perform when using market data. The primary reason for looking at market data is not to show that our method is good (although we do achieve this, at least partially), but rather to find the most interesting set of data for which to test numerically our results in the next section. We thus apply our method to calibrate our model on some specific data described next and compare the results to those obtained by using Hagan’s approximation recalled below.

3.1. Description of the data

The data that we use are options on the SP500 index for 2517 consecutive trading days. Each day has $260=2\times 10\times 13$ data points: 2 for the choice Call/Put, 10 for the choice of expiry time $T$:

[TABLE]

and 13 for the choices of $\Delta$, namely,

[TABLE]

(the value $y=\ln(Se^{rt}/K)$ is determined from $\Delta$ as explained below, Remark 3.1).

3.2. General description of the method

We run three types of tests using the market data. Each type of tests is distinguished by its objective function. More precisely, for the first type of tests we use the implied volatility as an objective function, for the second type of tests we use the actual price (obtained from the implied volatility and the Black-Scholes formula) as an objective function, and, finally, for the third type of tests we use the logarithm of the actual price. This is explained in more detail next.

For each of the three objective functions and each of the 2517 trading days, we use the least squares optimization to estimate the model parameters $(\nu,\sigma,\rho)$ to best fit the market data. (See Equation (65) for an explicit formula.) We list the fitting error (which is an “in-sample” error) in the first column of Table 1. We also use the parameters $(\nu,\sigma,\rho)$ to predict the market data for the next day, and then we list the corresponding (“out-of-sample”) error. We looked at different types of “errors” (more precisely, norms): $\ell_{1}$-, $\ell_{2}$-, and $\ell_{\infty}$-norms, but for simplicity we will mostly discuss the $\ell_{2}$-norm error. All these norms are normalized, in the sense that we always use averages (or probability measures). More specifically, for a sequence $a:=(a_{i})_{i\in I}$, where $I$ is a finite set with $|I|$ elements, its $\ell_{2}$-norm $\|a\|_{2}$ is given by

[TABLE]

(When complex numbers are used, $a_{i}^{2}$ is replaced with $|a_{i}|^{2}$, but that will not be necessary in what follows.) Similarly, the $\ell_{1}$- and the $\ell_{\infty}$-norms of $a:=(a_{i})_{i\in I}$ are given by $\|a\|_{1}:=|I|^{-1}\sum_{i\in I}|a_{i}|$ and by $\|a\|_{\infty}:=\max_{i\in I}|a_{i}|$, respectively. For the $\ell_{\infty}$-norm, the normalization plays, of course, no role.

Thus, for each trading day $\tau=1,\ldots,N$ we have a set of $J=260$ triples $\mathcal{M}_{\tau}=\{(y_{\tau,j},\Sigma_{\tau,j},T_{\tau,j})\}$, (so $j=1,\ldots,J$). Each triple in the above set represents the moneyness, the quoted volatility, and the time to expiry, respectively. Assuming that the strike is always $K=1$ and $r=0$, which does not change our calculations in view of Remarks 1.15 and 1.16, the price of the option for this trade is then obtained using the function $C_{rel}$ of Equation (45)

[TABLE]

More precisely, the procedure that we have outlined in the beginning of this subsection amounts to the following. For each objective function and each trading day $\tau$, we use the least squares optimization to find the the parameters $(\nu_{\tau},\sigma_{\tau},\rho_{\tau})$ for which the resulting objective function best matches (in an $\ell_{2}$ sense) the quoted data.

Before explaining this in even more detail for each of the three types of tests, let us make first the following remark about the structure of our data.

Remark 3.1**.**

In our data, the market data does not provide the log-moneyness value $y_{\tau,j}$, but rather the hedging parameter $\Delta_{\tau,j}$. The value $y=ln(Se^{rt}/K)$ is then estimated from $\Delta$ using the formula

[TABLE]

where $\sigma_{\tau-1}$ is the implied volatility determined in the previous trading day and $T=t$. For the first trading day we took for the parameters as an initial value the average value from a previous run of the program. We also estimated $y_{\tau,j}$ using $\Sigma_{\tau,j}$ instead of $\sigma_{\tau-1}$, but did not find any significant differences in the error estimates reported in Table 1.

3.3. The first type of market data tests: implied volatility

For the first type of tests, we have used two types of implied volatility approximations as an objective function: $\sigma_{D}$ and $\sigma_{H}$ defined below.

3.3.1. Objective function: $\sigma_{D}$

The first choice of implied volatility function, namely $\sigma_{D}$, is the one provided by Theorem 2.2 by truncating to the second order approximation in powers of $\nu$, that is,

[TABLE]

We used this objective function as follows.

Recall that $J_{\tau}$ is the set of market data for the day $\tau$. For each trading day $\tau$, we have computed the parameters $(\nu_{\tau},\sigma_{\tau},\rho_{\tau})$ that minimize the (square of the) $\ell_{2}$-error between our objective function ($\sigma_{D}$) and the market provided analog function (the implied volatility, denoted $\Sigma$ decorated with various indices in this case). That is, we have chosen $(\nu_{\tau},\sigma_{\tau},\rho_{\tau})$ to minimize

[TABLE]

The resulting minimum value for the $\ell_{2}$ norm is the fitting error of our model, or the “in-sample-error” (ISE), obtained by replacing $(\nu,\sigma,\rho)$ with $(\nu_{\tau},\sigma_{\tau},\rho_{\tau})$:

[TABLE]

We have also computed the “out-of-sample” (OSE) error obtained by replacing $(\nu,\sigma,\rho)$ with $(\nu_{\tau-1},\sigma_{\tau-1},\rho_{\tau-1})$, that is, with the parameters obtained in the previous day.

[TABLE]

The second line of the following table (Table 1) summarizes the results of our test, by providing the ISE and OSE, the average values $\overline{\nu}$, $\overline{\sigma}$, $\overline{\rho}$, as well as their standard deviations $std(\nu)$, $std(\sigma)$, and $std(\rho)$ for the specified objective function ($\sigma_{D}$ in this case).

3.3.2. Objective function: $\sigma_{H}$

Recall that in [31], Hagan, Kumar, Lesniewski, and Woodward applied a singular perturbation technique to derive a closed form implied volatility approximation under SABR dynamics, denoted $\sigma_{H}$. Let $z=\frac{\nu}{\sigma}\ln(F/K)$ and

[TABLE]

For $\beta=1$, the only case considered in this paper, their implied volatility approximation is given by the formula:

[TABLE]

(See also [1].) The last data line of at table, Table 1, summarizes the same information as in the previous paragraph, but for $\sigma_{H}$ replacing $\sigma_{D}$.

3.3.3. Conclusion

The calculations summarized in Table 1 show that the two formulas, $\sigma_{H}$ and $\sigma_{D}$, yield very similar results for the parameters appearing in the data set (this is not true in general!).

3.4. The second type of market data tests: actual prices

In this subsection we report the results of the tests that compared the actual prices The choice of an approximation of the implied volatility yields an approximation of the price through the Black-Scholes formula (43). For instance, if we approximate $\sigma_{imp}$ with $\sigma_{D}$ (as in the first set of tests in the previous subsection), the resulting approximation for the price is given by

[TABLE]

(In particular, $C_{D}=e^{-rt}F_{D}$, see Equation (16).) Analogously, if we approximate $\sigma_{imp}$ with $\sigma_{H}$, the resulting formula will be denoted

[TABLE]

The approximation $C_{H}$ of the solution of the SABR PDE is widely considered as very accurate for short dated options. We have compared the resulting values $C_{H}$ and $C_{D}$ with the prices $C_{M}$ provided by the market data. In addition to these approximate price functions, we have also tested a regularization $\tilde{C}_{H}$ of Hagan’s formula, as well as our second order approximation of $C_{SA}$, that is $C_{SA,2}=e^{-rt}F_{SA,2}$ of Theorem 1.11. As explained in Remarks 1.15 and 1.16, we can assume $K=1$ and $r=0$, up to multiplying all the results with a global factor.

3.4.1. Objective functions: $C_{D}$ and $C_{SA,2}$

We have then proceeded as in the previous subsection, but we have used $C_{D}$ instead of $\sigma_{D}$ and the actual prices $C_{M}=C_{BS}(S,K,\Sigma,t)$ instead of $\Sigma$, where $\Sigma$ is the implied volatility and is provided by the data set. We have also included a test for the formula for the approximate price including the mean reverting term for $\kappa=.25$. In our tests, we have assumed $Ke^{-rt}=1$. Instead of (65), we have thus chosen $(\nu,\sigma,\rho)$ to minimize

[TABLE]

where $C_{rel}$ is as defined in Equation (45) (we have thus replaced $\sigma_{D}(y_{\tau,j},\nu,\sigma,\rho,T_{\tau,j})-\Sigma_{\tau,j}$ with $C_{rel}(y_{\tau,j},\sigma_{D}(y_{\tau,j},\nu,\sigma,\rho,T_{\tau,j}),T_{\tau,j})-C_{rel}(y_{\tau,j},\Sigma_{\tau,j},T_{\tau,j})$ in Equation (65)). To obtain the In Sample Error (ISE) and the Out of Sample Error (OSE), we have performed some similar modifications to Equations (66) and (67). That has amounted to replacing $(\nu,\sigma,\rho)$ in Equation (72) with $(\nu_{\tau},\sigma_{\tau},\rho_{\tau})$ in order to obtain the ISE and with $(\nu_{\tau-1},\sigma_{\tau-1},\rho_{\tau-1})$ in order to obtain the OSE. We proceeded similarly for $C_{SA,2}=e^{-rt}F_{SA,2}$ (see Equation (15)). The results are provided in the last two columns of Table 2.

3.4.2. Objective functions: $C_{H}$ and $\tilde{C}_{H}$

The similar results for $C_{D}$ replaced with $C_{H}$ are provided in the second column of that table. It turns out that there are some numerical instabilities in the formula for $\sigma_{H}$ due to the quotient $z/\xi(z)$, since $\xi(0)=0$, so for very small $z$ we get into issues of machine precision. For this reason, for $z$ very small, we interpolate between $z/\xi(z)$ and $1-\frac{\sigma z}{2}\approx z/\xi(z)$ to obtain a new formula $\tilde{C}_{H}$, for which we also summarize the (slightly better) results. When dealing with market data, $z$ is, in fact, never zero, but when dealing with numerical tests, we do get $z=0$, so this interpolation becomes indispensable.

3.4.3. Conclusion

Replacing $C_{H}$ with $\tilde{C}_{H}$ improves the performance of this formula by a small, but significant amount. In any case, all these models (except the Black-Scholes model) perform in a very similar way, as seen by examining the In Sample Errors (ISE) and the Out of Sample Errors (OSE).

3.5. The third type of market data tests: log-prices

The reader may have already observed the fallacy of the method used in the previous section: by looking at the differences $C_{M}-C_{D}$, for example, we do not take account the fact that we want a smaller error when $C_{M}$ is small. Because of this, we now perform the same tests, but for the logarithm of the prices. We thus replace $C_{D}$ with $\ln(C_{D})$ (and similarly for the other $C$s). Thus, for the objective function $\ln(C_{D})$, we minimize

[TABLE]

All the other formulas for the other objective functions change in a similar way. We have also tested the mean reverting term. We obtain the results summarized in Table 3, with the mean reverting term in the row corresponding to $\ln(C_{\kappa})$.

Here are some comments on these results. The best model is the modified Hagan model, followed closely by our implied volatility model (so it is an advantage in this case to consider rather the implied volatility formula of Theorem 2.2 instead of the approximate price formula of Theorem 1.11, although they differ by a term of order $O(\nu^{3})$). For comparison, on the last line, we also show the performance of the Black-Scholes model, which is seen to be significantly worse than all the other ones. Using an expansion in $\nu$ of the implied volatility improves the results compared to the similar approximation of the price. Introducing the mean reverting term also improves (for small $\kappa$) the results. This is seen by comparing the row corresponding to $\ln(C_{\kappa})$ with the row corresponding to $\ln(C_{SA,2})$, with the formula for $C_{\kappa}$ being obtained from the formula for $C_{SA,2}$ by including the mean reverting term.

Remark 3.2**.**

We have included the information on the average values and the standard deviations of the parameters $(\nu,\sigma,\rho)$ since we will use them to select the most relevant parameters for our numerical tests. In this sense, we notice that $y$ has mean $-.01$ and standard deviation $0.28$.

4. Numerical tests

We have tested numerically our approximation formulas $C_{SA,2}$ and $C_{D}$ in several ways. The numerical tests confirm the efficiency of our method in a suitable range of the parameters. In all our numerical tests below, we have chosen $K=1$ and $\kappa=0$. We have also chosen $r=0$, so the distinction between the forward prices $F=e^{rt}C$ and the actual prices $C$ disappears.

4.1. The residual of the approximations: substituting in the

PDE

The simplest numerical test to estimate the performance of the various methods considered in this paper is to check if the resulting price function $C$ satisfy the SABR PDE (Equation (1)). Recall that $L$ is the generator of the SABR PDE and $\kappa=0$ in our numerical tests. That is, we compute the residual $\partial_{t}C-LC$, where $C=C_{H}$ for Hagan’s model, $C=C_{D}$ for our model using the approximate volatility $\sigma_{D}$, and $C=C_{SA,2}$ for our second order approximation. If $C$ was an exact solution of the SABR PDE, the residual would be zero. So the norm $R=\|\partial_{t}C-LC\|$ of the residual gives us an idea how far $C$ is from the actual solution. For comparison, we compute also the norm of the residual when $C=C_{BS}$, that is, when $C$ is given simply by the Black-Scholes formula. One can say that if the norm of the residual is small, then $C$ is close to the actual solution, since the PDE is well-posed for $\kappa=0$. The converse is not true, however, in that a large residual does not imply that $C$ is far from the solution (this is the same phenomenon as the one that gives that $f$ is small if $f^{\prime}$ is small, but not the other way around).

We computed the $\ell_{2}$-norm of the residual over various regions. For instance, let us consider the set of data points for which $\nu=0.125$, $T$ belongs to a set of equally spaced nodes between $0.1$ and $1$, $\rho=-0.4$, $\sigma$ belongs to a set of equally spaced nodes between $0.1$ and $0.3$, and $y$ belongs to a set of equally spaced nodes between $-0.5$ and $0.5$. The norms of the residuals for our choices of $C$ are summarized in Table 4. (In that table, we showed the norms times 1000, since the numbers were otherwise very small.)

We performed the same test for various other regions, including for regions in which the values of the parameters were close to the market data. The results, as expected, get worse with the increase of $\nu$, of $T$, and of the size of the interval for $y$. In these tests, our model behaves better for small $\nu$, even for large $T$. However, for larger values of $\nu$, Hagan’s model behaves better (for this type of tests).

In all these tests, we have used only values of $T\geq 0.1$, since the numerical differentiation used in the program becomes less reliable for $T<0.1$.

4.2. Comparison of our implied volatility formula with Hagan’s

formula for maket data

For all the parameters $q_{\tau,j}:=(\nu_{\tau},\sigma_{\tau},\rho_{\tau},y_{\tau,j},T_{\tau,j})$ determined as explained in the previous section, we have also computed the differences between the implied volatilities $\sigma_{H}-\sigma_{D}$, the prices

[TABLE]

We have similarly computed the differences between the corresponding log-prices:

[TABLE]

For all these options, we have computed the $\ell_{1}$, $\ell_{2}$, and $\ell_{\infty}$ norm, results that are summarized in Table 6

The conclusion is that, on average, $\sigma_{H}$ and $\sigma_{D}$ predict some very close values. However, occasionally, these values may be very different. The “best results” are obtained for the differences in price, but these are not too relevant, since they are not dimensionless (they do not take into account the magnitude of the prices).

4.3. A Monte Carlo simulation

None of the formulas $C_{H}$, $C_{D}$, or $C_{SA,2}$ is an exact solution of the SABR PDE. In order to compare them to the true solution, we need a method to approximate very well the true solution, so to see which of the methods approximates it best. Recall that $\kappa=0$ in our numerical test.

In a first set of tests, we have used the Monte Carlo Method to approximate the true solution of the SABR PDE. To carry out the tests, we have first computed prices of call options through Monte Carlo simulation (denoted below by $C_{MC}$) using 30000 paths and time step $10^{-4}$. We took thise values as the benchmark. We have then calculated the differences $E_{H}=C_{H}-C_{MC}$ and $E_{D}=C_{D}-C_{MC}$ for a range of $K$ and $t$ and ploted them as a function of the moneyness $y=\ln{F/K}$. We choose the following parameters throughout the test: $F=10$, $\sigma=0.2$, $\nu=0.2$, $\rho=-0.3$.

The approximations $C_{H}$ and $C_{D}$ agree almost exactly when $t=1$, and are both fairly accurate. The largest error is 0.8%, or around 0.1% to 0.2% of the benchmark price $C_{MC}$. As $t$ increase to 3 and 10, the error starts to increase and the two approximations gradually diverge. When both $C_{H}$ and $C_{D}$ overestimate $C_{MC}$, Hagan’s approximation tends to give a better approximation, whereas if they both underestimate $C_{MC}$, Duhamel-Dyson perturbative series approximation has a smaller error. So overall it is not possible to tell which of the two approximations of the implied volatility (Hagan’s and ours) is better for small $t$. Finally, when $t$ is as large as 30 years. The difference between $C_{H}$ and $C_{D}$ becomes significant. It is interesting to see that both $C_{H}$ and $C_{D}$ overestimate $C_{MC}$ most of the time and, moreover, $C_{D}$ is almost systematically below $C_{H}$. As a result, Duhamel-Dyson perturbative series method in this case gives a better approximation for a majority of strikes. Indeed, the largest relative error for Hagan’s approximation is 22%, however for Duhamel-Dyson perturbative series expansion method is only 12%.

For $\nu=1.5$, our Monte-Carlo simulation, however, does not converge even for $10^{6}$ paths. For this reason we tried also a Finite Difference simulation.

4.4. Comparison with the Finite Difference approximate solution

In a second set of tests, we took as benchmark a Finite Difference approximation of the solution to the SABR PDE. This is the most relevant method (so we saved the best for the last), since the FD method is more precise than the Monte Carlo method and, in our case, leads to a rather good approximation of the solution (but not perfect).

In fact, we have computed a sequence $w_{k}$ of Finite Difference (FD) approximations of the solution $F_{SA}$ of Equation (1) with $\kappa_{0}=0$ obtained by successively refining our set of nodes. For this sequence $w_{k}$ of approximations, we have tested the convergence of the FD approximations and we have compared them to the various approximations $C$ that we have considered before: $C_{SA,2}$, $C_{D}$, $C_{H}$, and $C_{BS}$. The results of these comparisons for the finest discretization (largest $k$) are shown in Tables 7 and 8. Let us now explain the results included in those talbles. We continue to assume $r=0$, so $F_{SA}=C_{SA}$, $F_{D}=C_{D}$, and so on.

4.4.1. Outline of the method and of its challenges

The tests in Section 3 (using market data) have given for $y$ a range approximately contained in $[-0.3,0.3]$. Therefore, in our FD tests, we have taken $I:=[-1,1]$ as an interval of interest for $y$. That means that we have compared our FD approximations $w_{k}$ with the values $C$ predicted by the other methods only for $y\in[-1,1]$. We similarly took as an interval of interest for $\sigma$ the interval

[TABLE]

(that is, $J$ is an interval symmetric with respect to $0.18$ on a logarithmic scale; but note that we have rounded off our values to two significant digits). In view of our results involving the market data (see Section 3), the choice of $0.18$ for the average volatility seems to be a reasonable choice. We have then compared our FD approximations $w_{k}$ with the the other approximations $C$ (namely $C_{H}$, $C_{V}$, $C_{SA,2}$, and $C_{BS}$) at the nodal points of

[TABLE]

The parameters $\nu,\rho,T$ have also been chosen to be close to the ones provided by the market data (as we will explain below).

In the implementation and error analysis, we had to deal with the fact that the domain is non-compact and that the coefficients of the SABR PDE are very far from being constant, in fact, they are very small in some regions and very large in others. Let us explain next how we specifically dealt with these issues.

4.4.2. Cut-off errors

Recall that the equation that we want to solve is Equation (1) (for $\kappa=0$, in these tests) on the domain $(x,\sigma)\in\mathbb{R}\times(0,\infty)$. Since the domain is non-compact, in order to discretize this equation using Finite Differences (FD), we have performed a cut-off of the domain by restricting to a rectangular domain $\Omega:=[-x_{MAX},x_{MAX}]\times[\sigma_{min},\sigma_{MAX}]$ with $\sigma_{min}\sigma_{MAX}=0.18^{2}$ and $\Omega$ significantly larger than $I\times J$. This required us then to specify the values for the solution on the boundary $\partial\Omega$ of the domain $\Omega$, something that had not been needed before the cut-off (since, instead of prescribing boundary conditions, we had assumed that our solution $u(t)$ is in a suitable weighted Sobolev space). We have thus replaced our problem (1) with

[TABLE]

By classical results, see [8, 15, 17, 39, 45, 47, 48] for instance, the FD approximations of this equation converge to the solution $w$ of Equation (75), something that we have noticed very clearly in our implementation.

The solution $w$ of Equation (75) is, however, different from the solution $C_{SA}$ of Equation (1) for $v=h$ and with $\kappa=0$. Nevertheless, in view of Remark 1.6 and of the exponential decay of the Green function of the equation $\partial_{t}-L$, one can prove for $g$ not too far from the actual solution on the boundary (even $g=0$ will do) that $\|u-w\|_{I\times J}\to 0$ as $\Omega$ approaches the total domain $\mathbb{R}\times(0,\infty)$ (see [12, 46] and the references therein. See also [19]. Here the norm is the normalized $\ell_{2}$ norm defined using the nodal points in the domain of interest $I\times J$. Therefore, a large cut-off will not affect too much the desired values on the small rectangle of interest $I\times J$. The question, however, is how large this cut-off should be chosen in implementation. To answer this question, we would need precise estimates on the cut-off error. Rigorous estimates for this cut-off error are very difficult to obtain and the simplest estimates seem to be, in any case, much larger than what we have seen in our implementation. Instead, we have chosen to numerically estimate the size of $\Omega$ by repeatedly increasing it until our FD solution changed very little on the domain of interest $I\times J$.

Specifically, we have then taken $\Omega:=[-3,3]\times[0.0193,1.6803]$ for the largest value of $\nu$ that we have tested ($\nu=1.5$) and for the largest value of $T$ that we have tested ($T=5$). Some smaller intervals have been used for some of the smaller value of $T$ and $\nu$ when this changed only marginally the results. Notice that, on a logarithmic scale, the interval in $\sigma$ has been chosen significantly larger than the corresponding interval in $x$ (when compared to the intervals of interest), so the discretization in $\sigma$ required more nodes than the one in $x$.

We found that the method converged faster if the boundary condition $g$ was given by the solution of the Black-Scholes equation on the boundary of $\Omega$, which corresponds to the solution of our PDE for $\nu=0$.

4.4.3. Discretization errors

For the discretization, be have begun with a uniform grid on the interval $[-3,3]$ with mesh size $h=1/2$. However, on the interval $[0.0193,1.6803]$, we have chosen a geometric progression grid with $19$ nodes (thus the end points of our interval of interest, whose precise values are $0.1404$ and $0.2307$, are among the chosen nodes). For $T=5$, we have used $190$ time-steps for the first iteration. We have used a direct method, because it is easy to implement and reasonably efficient in our case, at least after choosing the right grid in the $\sigma$ variable (i.e. a geometric progression grid). We have then successively refined our grid by dividing each interval into two smaller intervals, while preserving the nature of the grids in each dimension. (That is, in the $x$ direction, we have chosen the arithmetic mean of the end points to divide an old interval, whereas, in the $\sigma$ direction, we have used the geometric mean.) Correspondingly, we had to multiply by 4 the number of time steps each time when we performed a refinement in order to satisfy the stability condition for the FD implementation. This has thus lead us to a sequence $w_{k}$ of approximations of the solution $w$ of Equation (75). In particular, each approximate solution $w_{k}$ had required (essentially) 16 times more time to compute than the previous one. This has imposed some stringent limits on the precision of our approximations due to the running time of the code (which had been written in C++ and run on a simple laptop). The expected rate of convergence $\|w_{k+1}-w_{k}\|\approx 4^{-1}\|w_{k}-w_{k-1}\|$ has been observed almost immediately and that has lead to good error estimates of the form $\|w-w_{k}\|\approx 3^{-1}\|w_{k}-w_{k-1}\|$ (unlike the cut-off error, which was estimated empirically, although we knew that it decayed faster than any exponential due to the decay of the Green function of $\partial_{t}-L$). This rate of convergence was observed both on the domain of interest $I\times J$ and on the total domain $\Omega$.

To see the kind of precision obtained, we included the entry that begins with $1*$ in table 7. It corresponds not to the last term $w_{k}$ in our sequence of FD approximations, but rather to $w_{k-1}$. The entry corresponding to $w_{k}$ is right underneath it, and we see that the numbers are very close (except the estimated error $\|w-w_{k}\|_{I\times J}$, which is, as expected, about four times smaller for the $k$th term). Notice that the values recorded in that table represent 100 times the values of the corresponding norms, since we were dealing with small numbers.

4.4.4. Total error estimation

Using the notation of the previous paragraphs, we see that the total error (on the domain of interest) is

[TABLE]

In this estimate, $\|C_{SA}-w\|_{I\times J}$ is the cut-off error, due to our restricting the FD test to a bounded domain $\Omega$. The cut-off error is independent of $k$ in our tests, so the domain $\Omega$ had to be chosen rather large to start with. The term $\|w-w_{k}\|_{I\times J}$ represents the error due to the FD discretization and goes to 0 as expected as $c4^{-k}$. The resulting estimates for $\|C_{SA}-w_{k}\|_{I\times J}$ are included in the last columns of Tables 7 and 8. Let us explain the meaning of the other columns in those tables.

Let $\|\,\cdot\,\|:=\|\,\cdot\,\|_{\ell_{2}(I\times J)}$ and $\|\,\cdot\,\|_{\infty}:=\|\,\cdot\,\|_{\ell_{\infty}(I\times J)}$. Let $C_{H}$ be the solution predicted by Hagan’s model and $C_{SA,2}$ be the solution predicted by our Duhamel-Dyson series model, as before. The “log” columns show the $\|\ln(C_{H})-\ln(w_{k})\|$ and $\|\ln(C_{SA,2})-\ln(w_{k})\|$ norms of the differences at the nodal points, where $k$ corresponds to the last discretization. We also write $\delta C_{H}:=C_{H}-w_{k}$ and, similarly, $\delta C_{SA,2}:=C_{SA,2}-w_{k}$ and $\delta_{HD}:=C_{H}-C_{SA,2}$. We have performed our tests for $\rho=-.2$ and for different values of the of the other parameters. The results are included in Table 7. Similar results, but for $\rho=-.5$ are included in Table 8, where we show separately the estimates for the FD error and for the cut-off error. Because we the numbers were very small, all norms in these tables were multiplied by 100.

The above tests (taking into account also the discretization and cut-off errors) seem to suggest that the Duhamel-Dyson method is at least as competitive as Hagan’s method (if not better for large $T$ and parameters compatible with market data). The quantities that seem the most relevant to support our conslusions are thus $\|\delta C_{H}\|$ and $\|\delta C_{SA,2}\|$.

The results of the tests in Tables 7 and 8 sometimes allow us to compare the different other approximation methods as follows. Let us do that for Hagan’s formula and for our second order approximation $C_{SA,2}$. We have (all norms are on the domain of interest $I\times J$)

[TABLE]

This gives that $C_{SA,2}$ will be closer to the actual solution $C_{SA}$ (on the grid points of $I\times J$) whenever we have $\|\delta C_{H}\|-\|\delta C_{SA,2}\|>2\|C_{SA}-w_{k}\|$. On the other hand, if the norm $\|C_{H}-C_{SA,2}\|$ is small, it will be hard to decide which of $\|C_{H}-C_{SA}\|$ and $\|C_{SA,2}-C_{SA}\|$ is smaller. This is the case when $T$ and $\nu$ are small. In general, the difference $\|C_{H}-C_{SA,2}\|$ goes faster to 0 than $\|C_{SA}-w_{k}\|$, so for these values, it is more difficult to tell which of $C_{H}$ or $C_{SA,2}$ is closer to $C_{SA}$. Thus, for small values of $\nu$ or $T$, we need more precision (i.e. a larger $k$) in order to distinguish between models, both of which seem to approximate very well the actual solution (for parameters compatible with the market data, in particular, for $T\leq 2$). Since increasing $k$ by 1 increases the time for running the program approximately by a factor of 15, distinguishing $C_{H}$ and $C_{SA,2}$ for small $T$ or $\nu$ is a true challenge in our tests.

Another challenge in the FD implementation is that the cost of the method increases very fast with $T$. If we could keep the cut-off domain fixed (with $T$), the cost would grow linearly in $T$ (which is already a challenge). However, as $T$ growth, we need to take a larger $\Omega$, since the Green function decreases roughly like $e^{-d^{2}/(2T)}$, where $d$ is the hyperbolic distance on $\mathbb{R}\times(0,\infty)$. To maintain the same precision for a larger $T$, we would expect then the distance to $\partial\Omega$ to grow at least as fast as $\sqrt{T}$. So if replace $T$ by $4T$, then we would have to double the (hyperbolic) distance to $\partial M$. That would mean to multiply $\sigma_{MAX}$ by roughly $e^{2\nu}$. Since we had chosen a geometric mesh in the $\sigma$ direction (exactly for this reason), that would not increase too much the number of nodes in the $\sigma$ direction. (This number growth like $\ln(\Omega_{MAX})$.) However, it would affect the stability of the matrix in the (direct) FD method by a factor of $e^{4\nu}$ (because of the coefficient $\sigma^{2}$ in front of $\partial_{x}^{2}$) and thus would require a proportional decrease in the size of the time steps. The total number of time steps would hence need to be thus increased by a factor of $4e^{4\nu}$. When $\nu$ is large, this is very expensive.

See [21, 33, 54] for some recent results on the mean-reverting case. We stress, however, than many theoretical results that are true in the non-mean-reverting case (i.e. $\kappa_{0}=0$) are not known and even may not true in the case $\kappa_{0}\neq 0$.

5. Extensions of the method

So far, we have completely solved a special case of the much more general $\lambda$-SABR model. That is, we have taken $\beta=1$ and $\kappa=0$. However, the applicability of the commutator approach does not require $F(t)$ and $\sigma(t)$ to be log-normal. As mentioned before, the key condition for this method to work is that the model coefficients are all polynomial functions of the state variables. At first glance, the SABR model with a general $\beta$ violates this condition because the coefficient, $F(t)^{\beta}$, is not a polynomial function of $F(t)$. Fortunately, we will show that, if the volatility process is log-normal (i.e. $\kappa=0$), a simple change of variable transforms the SABR equations into one with polynomial coefficients.

The other assumption, that $\kappa=0$, is not crucial either. When, $\kappa\neq 0$, the mean-reverting drift term can be eliminated by considering instead the forward volatility $\mathbb{E}_{t}[\sigma(T)]$ as our state variable. This creates slight complication in our model, as the transformed SDE’s will not be time-homogeneous. We will show that this is not an essential issue for the application of the commutator method.

With the above being said, the most general case, that $\kappa\neq 0$ and $\beta\neq 1$, is substantially more difficult and cannot be dealt with by the commutator method alone. Therefore, in the rest of the section, we consider only generalizations along the two directions, 1) $\kappa\neq 0$ and $\beta=1$, 2) $\beta\neq 1$ and $\kappa=0$.

5.1. The log-normal model with mean-reverting volatility

We consider first the simpler case that $\beta=1$, while the volatility process, instead of being log-normal, is driven by a process similar to an Ornstein-Uhlenbeck process. Then our price-volatility dynamics under the risk-neutral measure is assumed to be

[TABLE]

To eliminate the volatility drift term, we consider the forward volatility, defined as the time-$t$ conditional expectation of the terminal volatility $\sigma(T)$,

[TABLE]

Then $z(t)$ is a martingale by construction, and hence there is no drift in its dynamics. To derive the exact relationship between $z(t)$ and $\sigma(t)$, we calculate

[TABLE]

Ingrate the above from $t$ to $T$, and take the conditional expectation $\mathbb{E}_{t}[\cdot]$ on both sides, we obtain

[TABLE]

The above allows us to write the variable $\sigma(t)$, as a function of our new state variable $z(t)$ and time, as

[TABLE]

It is easy to verify that

[TABLE]

Therefore, we arrive at the transformed model,

[TABLE]

We are now back at the previous case where the “volatility” process has no drift, although the coefficients are now time-dependent functions of the state variables. After replacing the price process $F(t)$ with $X(t)=\ln{(F(t))}$, the option price $u(x,z,t)$ under the transformed model solves the following initial value problem

[TABLE]

The corresponding differential operators are

[TABLE]

and

[TABLE]

Keep in mind that $\sigma$ here is a function of $z$ and $\tau$ given by (79), and the $L_{i}$’s are no longer time-homogeneous. Therefore the fundamental solution of the system is no longer a semi-group, but an evolution system, in the form of $e^{\int_{0}^{\tau}L(s)ds}$. Repeated application of Duhamel’s principle yields the following Duhamel-Dyson perturbative series expansion of $e^{\int_{0}^{\tau}L(s)ds}$

[TABLE]

Applying again the CHB identity, we can easily arrive at the following lemma

Lemma 5.1**.**

[TABLE]

with

[TABLE]

Proof.

The proof is a straightforward calculation using CHB identity and the fact that $u(x,z,0)$ depends on $x$ only. ∎

Lemma 5.2**.**

The commutators are given by:

[TABLE]

Proof.

Let us first calculate the time-indexed commutator $[L_{1}(t_{1}),\,L_{0}(t_{0})]$. For brevity, we use the notation $v_{i}$ to denote $\sigma(z,t_{i})^{2}$, and use $v_{i}^{\prime}$ to denote $\frac{\partial\,\sigma(z,t_{i})^{2}}{\partial\,z}$. Recall that $\sigma(z,\tau)=e^{\kappa\tau}z-\theta(e^{\kappa\tau}-1)$. As a result,

[TABLE]

[TABLE]

Hence

[TABLE]

Rearranging the terms yields the first identity. Calculations for the second and third identities are similar. ∎

Combing the above two lemmas, $J_{1}$ and $J_{2}$ can be obtained through straightforward (but lengthy) calculations. Then, same as before, the operator $1+\nu J_{1}+\nu^{2}J_{2}$ can be written as a linear combination of $\partial_{x}$, $\partial_{x}^{2}$, $\ldots$, $\partial^{6}_{x}$, with each coefficient being elementary functions of parameters. The only question left is, how do we calculate the transition density $e^{\int_{t}^{T}L_{0}(s)ds}(x,t;y,T)$, in the presence of volatility mean-reversion?

Note that when $\nu=0$, we arrive at the Black-Scholes model with deterministic volatility, given by

[TABLE]

Hence the terminal log-price, $X(T)$ can be solved explicitly as,

[TABLE]

Conditional on $\mathcal{F}_{t}$, $X(T)$ is normally distributed with mean $X(t)-\frac{1}{2}\overline{V}(t)$ and variance $\overline{V}(t)$, with $\overline{V}(t)$ denote the time aggregated variance. Then the transition density in the case of $\nu=0$ is given by

[TABLE]

Here $\overline{V}(t)$ can be obtained by solving the ODE (85), as

[TABLE]

$e^{\int_{t}^{T}L(s)ds}(x,t;y,T)$ as well as the option price follow from lemmas 5.1 and 5.2 and equation 86.

5.2. The SABR model (with general $\beta$)

We now consider the original SABR model with the volatility following a pure log-normal process, i.e.

[TABLE]

With $\beta\neq 1$, the model can no longer be thought of as the Black-Scholes model perturbed by the vol-of-vol coefficient. Indeed, when $\nu=0$, (88) degenerates to the CEV model which is itself hard to solve. Therefore, our first objective is to apply a change of variable technique such that the transformed model becomes again a perturbed Black-Scholes model.

First notice that the “instantaneous volatility” is not $\sigma(t)$, but $\sigma(t)F(t)^{\beta-1}$. We call it the local volatility process and denote it by $M(t)$. The good news that, when $\sigma(t)$ follows a driftless log-normal process as in (88), the Ito-differential of $M(t)$ does not depend on $\sigma(t)$ and $F(t)$ individually, but on $M(t)$ alone. To see this, apply Ito’s lemma

[TABLE]

Let $\epsilon$ denote $\beta-1$, the SABR model is now rewritten as

[TABLE]

Therefore, we can now consider the SABR model as the Black-Scholes model perturbed by two “vol-of-vol” parameters, $\epsilon$ and $\nu$. Moreover, the coefficients in the above equation are all polynomial functions of $F(t)$ and $M(t)$. As a result, if we perform a Taylor series expansion for the option price, $C_{SA}(f,m)$ on both $\epsilon$ and $\nu$, the CHB theorem guarantees that the coefficients can be calculated exactly with finitely many terms.

Note that the above result no longer holds when the volatility is mean-reverting, as the local volatility process $M(t)$ is no longer autonomous. We also point out that, unlike in the previous case with mean-reverting volatility, we do not try to get rid of the “volatility drift” term in the $M(t)$-equation. The reason being that we are now performing Taylor expansion on both $\epsilon$ and $\nu$, hence the drift of $M(t)$ does not enter our zeroth order operator and complicate the rest of the calculations.

Equation (90) implies that the option price, $u(x,m,\tau)$ solves the following initial value problem

[TABLE]

The second order differential operator on the right hand side can be written as,

[TABLE]

where

[TABLE]