ABSTRACT

We develop two new pricing formulae for European options. The purpose of these formulae is to better understand the impact of each term of the model, as well as improve the speed of the calculations. We consider the SABR model (with $\beta =1$ ) of stochastic volatility, which we analyze by tools from Malliavin Calculus. We follow the approach of Alòs et al. (2006) who showed that under stochastic volatility framework, the option prices can be written as the sum of the classic Hull-White (1987) term and a correction due to correlation. We derive the Hull-White term, by using the conditional density of the average volatility, and write it as a two-dimensional integral. For the correction part, we use two different approaches. Both approaches rely on the pairing of the exponential formula developed by Jin, Peng, and Schellhorn (2016) with analytical calculations. The first approach, which we call “Dyson series on the return’s idiosyncratic noise” yields a complete series expansion but necessitates the calculation of a 7-dimensional integral. Two of these dimensions come from the use of Yor’s (1992) formula for the joint density of a Brownian motion and the time-integral of geometric Brownian motion. The second approach, which we call “Dyson series on the common noise” necessitates the calculation of only a one-dimensional integral, but the formula is more complex. This research consisted of both analytical derivations and numerical calculations. The latter show that our formulae are in general more exact, yet more time-consuming to calculate, than the first order expansion of Hagan et al. (2002).

Keywords:

SABR Model, Stochastic Volatility, Malliavin Calculus, Exponential Formula, Option Pricing

1. Introduction

European options are traditionally priced and hedged by Black-Scholes [1] (1973) model, one of the natural extensions of the Black-Scholes model to make volatility stochastic. The simplest stochastic volatility models assume that the volatility and the noise driving stock prices are uncorrelated. Moreover, the Hull-White formula [2] (1987) establishes that the European option price is the expectation of the Black-Scholes option pricing formula with a time-dependent volatility. An important success of this model is that it calculates European prices which implied volatilities smile. The development of local volatility models by Dupire and Derman (1994) was a major development in handling smiles and skews. However its predictions contradict empirical findings. Thus the SABR (stochastic alpha beta rho) model, a stochastic volatility model in which the asset price is correlated with its volatility was derived by Hagan et al. [3] (2002) to resolve this problem. Alòs [4] (2006) extended the classical Hull-White formula to the correlated case by means of Malliavin calculus. The new generalization decomposes option prices as the sum of the same derivative price if there was no correlation and a correction due by correlation. Another popular model is the Heston (1993) model. In that model, the volatility is mean-reverting. The general asymptotic method presented by Fouque, Papanicolau and Sircar (2000) [5] can be used to analyze Heston’s model. For more information on stochastic volatility models, we refer the reader to Gatheral [6] (2006).

Nevertheless, there are still terms of conditional expectation of functions of non-adapted processes in the new generalization of Hull-White formula. Jin, Peng and Schellhorn [7] (2016) showed that under certain smoothness conditions, a Brownian martingale can be represented via an exponential formula when evaluated at a fixed time. It is a powerful tool similar to Clark-Ocone formula that allows us to work with the conditional expectation of a random variable instead of the random variable itself.

The main goal of this research was to obtain an option pricing formula for the special case of the SABR model with $\beta =1$. We used two different approaches. Both approaches rely on the pairing of the exponential formula developed by Jin, Peng, and Schellhorn (2016) with analytical calculations, and start by conditioning on the path of the common noise term W. In the first approach, which we call “Dyson series in the return’s idiosyncratic noise’’, we first apply a Dyson series in the idiosyncratic noise term Z and then apply Yor’s [8] formula (1992) for the joint density of a Brownian motion and the time-integral of geometric Brownian motion to integrate with respect to the common noise term W. We note that Yor’s formula is used for pricing Asian options, but it is ideally suited to analyze realized volatility in the SABR model with $\beta =1$, since volatility is a geometric Brownian motion. Faà di Bruno’s formula is used for analytical differentiation. The first approach yields a complete series expansion but necessitates the calculation of a 7-dimensional integral. Two of these dimensions come from the analytical expression of the joint density of a Brownian motion and the time-integral of geometric Brownian motion. In the second approach, which we call “Dyson series in the common noise’’, we first integrate away the idiosyncratic noise term Z and then apply a Dyson series in the common noise term W. This results in a formula which necessitates the calculation of only a one-dimensional integral, but the formula is more complex, and we carried the calculation only of the first term of the series.

The organization of this paper is as follows. In Section 2 we present a brief introduction to Malliavin Calculus as well as a representation theorem for smooth Brownian Martingales. Section 3 is a review of basic option pricing theory and an extension to stochastic volatility models. In Section 4, we present several Hull-White formulas for European call option prices with different model assumptions. In Section 5, we derive the Dyson series in the return’s idiosyncratic noise for the call price. In Section 6, we derive the Dyson series in the common noise for the call price, and compare numerically all approaches.

2. Preliminaries on Malliavin Calculus

The following section briefly reviews some basic facts of Malliavin Calculus required along the paper. For a complete exposition we refer to Nualart [9] (1995) and Øksendal [10] (2008). Let $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}\mathrm{,}\mathbb{P}\right)$ be a complete filtered probability space where $\left\{{\mathcal{F}}_t\right\}$ is generated by a standard Brownian motion ${\left\{W_t\right\}}_{t\ge 0}$. In Section 2.4, we will enlarge our probability space to consider two standard Brownian motions.

2.1. Malliavin Derivative

Let $L^2\left({\left[\mathrm{0,}T\right]}^n\right)$ be the standard space of square integrable Borel real functions on ${\left[\mathrm{0,}T\right]}^n$ and let ${\tilde{L}}^2\left({\left[\mathrm{0,}T\right]}^n\right)\in L^2\left({\left[\mathrm{0,}T\right]}^n\right)$ be the space of symmetric square integrable Borel real functions on ${\left[\mathrm{0,}T\right]}^n$, consider the set

$S_n=\left\{\left(t_1\mathrm{,}\cdots \mathrm{,}{t}_n\right)\in {\left[\mathrm{0,}T\right]}^n\mathrm{<mo>\; :\; </mo>0}\le t_1\le \cdots \le t_n\le T\right\}$.

Definition 2.1 If f is a deterministic function defined on $S_n\left(n\ge 1\right)$ such that ${\Vert f\Vert }_{L^2\left(S_n\right)}^2:={\displaystyle {\int }_{S_n}}{f}^2\left(t_1,\cdots ,t_n\right)\text{d}{t}_1\cdots \text{d}{t}_n<\infty$, then the n-fold iterated It integral is defined as

$J_n\left(f\right):={\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_n}}\cdots {\displaystyle {\int }_0^{t_3}}{\displaystyle {\int }_0^{t_2}}f\left(t_1,\cdots ,t_n\right)\text{d}{W}_{t_1}\cdots \text{d}{W}_{t_{n-1}}\text{d}{W}_{t_n},$ (1)

and if $g\in {\tilde{L}}^2\left({\left[\mathrm{0,}T\right]}^n\right)$ we define

$I_n\left(g\right)={\displaystyle {\int }_{{\left[0,T\right]}^n}}g\left(t_1,\cdots ,t_n\right)\text{d}{W}_{t_1}\cdots \text{d}{W}_{t_n}:=n!J_n\left(g\right).$ (2)

Theorem 2.2 The Wiener-It Chaos Expansion. Let F be an ${\mathcal{F}}_T$ -measurable random variable in $L^2\left(P\right)$. Then there exists a unique sequence ${\left\{f_n\right\}}_0^{\infty }$ of functions $f_n\in {\tilde{L}}^2\left({\left[\mathrm{0,}T\right]}^n\right)$ such that $F={\displaystyle {\sum }_0^{\infty }}{I}_n\left(f_n\right)$.

Definition 2.3 Let $u\left(t\right)$,$t\in \left[\mathrm{0,}T\right]$, be a measurable stochastic process such that for all $t\in \left[\mathrm{0,}T\right]$ the random variable $u\left(t\right)$ is ${\mathcal{F}}_T$ -measurable and $E\left[{\displaystyle {\int }_0^T}{u}^2\left(t\right)\text{d}t\right]<\infty$. Let its Wiener-It chaos expansion be

$u\left(t\right)={\displaystyle \underset{0}{\overset{\infty }{\sum }}}{I}_n\left(f_{n,t}\right)={\displaystyle \underset{0}{\overset{\infty }{\sum }}}{I}_n\left(f_n\left(\cdot ,t\right)\right).$ (3)

Then we define the Skorohod integral of u by

$\delta \left(u\right):={\displaystyle {\int }_0^T}u\left(t\right)\delta W_t:={\displaystyle \underset{0}{\overset{\infty }{\sum }}}{I}_{n+1}\left({\tilde{f}}_n\right),$ (4)

when converge in $L^2\left(P\right)$, we say that u is Skorohod integrable and we write $u\in Dom\left(\delta \right)$ if the series in (4) converges in $L^2\left(P\right)$.

The operator $\delta$ is an extension of the It integral, in the sense that the set $L^2\left(P\right)$ of square integrable and adapted processes is included in $Dom\left(\delta \right)$ and the operator $\delta$ restricted to $L^2\left(P\right)$ coincides with the It stochastic integral.

Theorem 2.4 Let $u=u\left(t\right),t\in \left[0,T\right]$, be a measurable $\mathbb{F}$ -adapted stochastic process such that $E\left[{\displaystyle {\int }_0^T}{u}^2\left(t\right)\text{d}t\right]<\infty$. Then $u\in Dom\left(\delta \right)$ and its Skorohod integral coincides with the It integral

${\displaystyle {\int }_0^T}u\left(t\right)\delta W_t={\displaystyle {\int }_0^T}u\left(t\right)\text{d}{W}_t\mathrm{.}$ (5)

Definition 2.5 Let $F\in L^2\left(P\right)$ be ${\mathcal{F}}_T$ -measurable with chaos expansion $F={\displaystyle {\sum }_0^{\infty }}{I}_n\left(f_n\right)$, where $f_{\mathrm{.}}\in {\tilde{L}}^2\left({\left[\mathrm{0,}T\right]}^n\right)$, for $n=1,2,\cdots$, we say that $F\in {\mathbb{D}}_{\mathrm{1,2}}$ if ${\Vert F\Vert }_{{\mathbb{D}}_{1,2}}^2:={\displaystyle {\sum }_0^{\infty }}nn!{\Vert f_n\Vert }_{L^2\left({\left[0,T\right]}^n\right)}^2<\infty$. If $F\in {\mathbb{D}}_{\mathrm{1,2}}$ we define the Malliavin derivative $D_{t}F$ of F at time t as the expansion

$D_{t}F={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}}n{I}_{n-1}\left(f_n\left(\cdot ,t\right)\right),t\in \left[0,T\right].$ (6)

We will need the following results on the Malliavin derivative.

Theorem 2.6 Product rule for Malliavin derivative. Suppose $F_1\mathrm{,}{F}_2\in {\mathbb{D}}_{\mathrm{1,2}}^0$. Then $F_1\mathrm{,}{F}_2\in {\mathbb{D}}_{\mathrm{1,2}}$ and also $F_1{F}_2\in {\mathbb{D}}_{\mathrm{1,2}}$ with

$D_t\left(F_1{F}_2\right)=F_1{D}_t{F}_2+F_2{D}_t{F}_1\mathrm{.}$ (7)

Theorem 2.7 Chain rule. Let $G\in {\mathbb{D}}_{\mathrm{1,2}}$ and $g\in {\mathcal{C}}^1\left(ℝ\right)$ with bounded derivative. Then $g\left(G\right)\in {\mathbb{D}}_{\mathrm{1,2}}$ and

$D_{t}g\left(G\right)=g^{\prime }\left(G\right)D_{t}G\mathrm{.}$ (8)

Example 2.1

$D_t{\left({\displaystyle {\int }_0^T}f\left(s\right)\text{d}{W}_s\right)}^n=n{\left({\displaystyle {\int }_0^T}f\left(s\right)\text{d}{W}_s\right)}^{n-1}{D}_t\left({\displaystyle {\int }_0^T}f\left(s\right)\text{d}{W}_s\right)=n{\left({\displaystyle {\int }_0^T}f\left(s\right)\text{d}{W}_s\right)}^{n-1}f\left(t\right)$.

Theorem 2.8 The fundamental theorem of calculus. Let $u=u\left(s\right),s\in \left[0,T\right]$, be a stochastic process such that $E\left[{\displaystyle {\int }_0^T}{u}^2\left(s\right)\text{d}s\right]<\infty$ and assume that, for all $s\in \left[\mathrm{0,}T\right]\mathrm{,}u\left(s\right)\in {\mathbb{D}}_{\mathrm{1,2}}$ and that, for all $t\in \left[\mathrm{0,}T\right]\mathrm{,}{D}_{t}u\in Dom\left(\delta \right)$. Assume also that $E\left[{\displaystyle {\int }_0^T}{\left(\delta \left(D_{t}u\right)\right)}^2\text{d}t\right]<\infty$. Then ${\displaystyle {\int }_0^T}u\left(s\right)\delta W_s$ is well-defined and belongs to ${\mathbb{D}}_{\mathrm{1,2}}$ and

$D_t\left({\displaystyle {\int }_0^T}u\left(s\right)\delta W_s\right)={\displaystyle {\int }_0^T}{D}_{t}u\left(s\right)\delta W_s+u\left(t\right)\mathrm{.}$ (9)

2.2. Exponential Formula

A Brownian motion martingale can be represented via an exponential formula when evaluated at a fixed time under certain smoothness conditions.

Definition 2.9 Given $\omega \in \Omega$, a freezing operator ${\omega }^t$ is defined as:

$W\left(s\mathrm{,}{\omega }^t\left(\omega \right)\right)=(\begin{array}{ll}W\left(s\mathrm{,}\omega \right)\mathrm{,}\hfill & \text{if}s\le t\mathrm{<mo>\; ;\; </mo>}\hfill \\ W\left(t\mathrm{,}\omega \right)\mathrm{,}\hfill & \text{if}t<s\le T\mathrm{.}\hfill \end{array}$ (10)

The freezing operator ${\omega }^t$ is a mapping from $\Omega$ to $\Omega$. The following equations show some properties of the freezing operator:

Proposition 2.10

1) For $p\in \mathcal{P}$, space of polynomials, suppose $F=p\left(W_{s_1},\cdots ,W_{s_n}\right)$, then $F\left({\omega }^t\right)=p\left(W_{s_1\wedge t}\mathrm{,}\cdots \mathrm{,}{W}_{s_n\wedge t}\right)$ ;

2) $\left({\displaystyle {\int }_0^T}f\left(s\right)\text{d}{W}_s\right)\left({\omega }^t\right)={\displaystyle {\int }_0^t}f\left(s\right)\text{d}{W}_s$ ;

3) $\left({\displaystyle {\int }_0^T}{W}_s\text{d}s\right)\left({\omega }^t\right)={\displaystyle {\int }_0^t}{W}_s\text{d}s+W_t\left(T-t\right)$ ;

4) $\left({\displaystyle {\int }_0^T}{W}_s\text{d}{W}_s\right)\left({\omega }^t\right)=\left(\frac{W_T^2-T}{2}\right)\left({\omega }^t\right)=\frac{W_t^2-T}{2}$.

We denote the Malliavin derivative of order l of F at time t by $D_t^{l}F$, as a shorthand notation for $D_t\cdots D_{t}F$. We call ${\mathbb{D}}_{\infty }\left(\left[\mathrm{0,}T\right]\right)$ the set of random variables which are ${\mathcal{F}}_t$ -measurable and infinitely Malliavin differentiable.

Definition 2.11 A random variable F is said to be infinitely Malliavin differentiable if for any integer n:

$E\left[{\left(\underset{s_1,\cdots ,s_n\in \left(t,T\right)}{\sup }\left|\left(D_{s_n}^2\cdots D_{s_1}^{2}F\right)\right|\right)}^2\right]<+\infty .$ (11)

In particular, we denote by ${\mathbb{D}}^N\left(\left[\mathrm{0,}T\right]\right)$ the space of all random variables F which satisfy (14) for all $n\le N$.

The next theorem, or exponential formula, was obtained by Jin et al. (2016). The resulting series (12) is called a Dyson series.

Theorem 2.12 Suppose $F\in {\mathbb{D}}_{\infty }\left(\left[\mathrm{0,}T\right]\right)$ satisfies the following condition:

$\frac{{\left(T-t\right)}^{2n}}{{\left(2^{n}n!\right)}^2}E\left[{\left(\underset{u_1,,u_n\in \left(t,T\right)}{\sup }\left|\left(D_{u_n}^2\cdots D_{u_1}^{2}F\right)\left({\omega }^t\right)\right|\right)}^2\right]\underset{n\to \infty }{\to }0,$

for fixed $t\in \left[\mathrm{0,}T\right]$, then

$E\left[F|{\mathcal{F}}_t\right]={\displaystyle \underset{i=0}{\overset{\infty }{\sum }}}\frac{1}{2^{i}i!}{\displaystyle {\int }_{{\left[t,T\right]}^i}}\left(D_{s_i}^2\cdots D_{s_1}^{2}F\right)\left({\omega }^t\right)\text{d}{s}_i\cdots \text{d}{s}_1.$ (12)

Example 2.2 An example of applying the Exponential formula: Let $F=W_T^2$, then for $t\le s\le T$ : $F\left({\omega }^t\right)=W_t^2$ and $\left(D_s^{2}F\right)\left({\omega }^t\right)=2$, then by Theorem 2.2 we have

$E\left[F|{\mathcal{F}}_t\right]=F\left({\omega }^t\right)+\frac{1}{2}{\displaystyle {\int }_t^T}\left(D_s^{2}F\right)\left({\omega }^t\right)\text{d}s=W_t^2+T-t.$ (13)

2.3. Faà di Bruno’s Formula

Lemma 2.13 Faà di Bruno’s formula. If f and g are functions with a sufficient number of derivatives, then

$\frac{{\text{d}}^n}{\text{d}{x}^n}f\left(g\left(x\right)\right)={\displaystyle \sum }\frac{n!}{{\displaystyle {\prod }_{i=1}^n{m}_i!}}{f}^{\left({\displaystyle {\sum }_{k=1}^n{m}_k}\right)}\left(g\left(x\right)\right)\cdot {\displaystyle {\prod }_{j=1}^n{\left(\frac{g^{\left(j\right)}\left(x\right)}{j!}\right)}^{m_j}},$ (14)

where the sum is over all n-tuples of non-negative integers $\left(m_1\mathrm{,}\cdots \mathrm{,}{m}_n\right)$ satisfying the constraint ${\displaystyle {\sum }_{k=1}^n}k{m}_k=n$. Combining the terms with the same value of ${\displaystyle {\sum }_{i=1}^n}{m}_i=k$ leads to a simpler formula expressed in terms of Bell polynomials $B_{n\mathrm{,}k}\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_{n-k+1}\right)$ :

$\frac{{\text{d}}^n}{\text{d}{x}^n}f\left(g\left(x\right)\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{f}^{\left(k\right)}\left(g\left(x\right)\right)\cdot B_{n,k}\left(g^{\prime }\left(x\right),g^{″}\left(x\right),\cdots ,g^{n-k+1}\left(x\right)\right).$ (15)

Definition 2.14 Exponential Bell polynomials. The partial or incomplete exponential Bell polynomials are a triangular array of polynomials given by

$B_{n,k}\left(x_1,x_2,\cdots ,x_{n-k+1}\right)={\displaystyle \sum }\frac{n!}{{\displaystyle {\prod }_{i=1}^{n-k+1}{j}_i!}}{\displaystyle {\prod }_{i=1}^{n-k+1}{\left(\frac{x_i}{i!}\right)}^{j_i}},$ (16)

where the sum is taken over all sequences $j_1\mathrm{,}{j}_2\mathrm{,}\cdots \mathrm{,}{j}_{n-k+1}$ non-negative integers such that these two conditions are satisfied: ${\displaystyle {\sum }_{i=1}^{n-k+1}}{j}_i=k$ and ${\displaystyle {\sum }_{i=1}^{n-k+1}}i\cdot j_i=n$. The sum

$B_n\left(x_1,\cdots ,x_n\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{B}_{n,k}\left(x_1,x_2,\cdots ,x_{n-k+1}\right),$ (17)

is called the $n$ th complete exponential Bell polynomials.

The Faà di Bruno’s formula can be generalized to Malliavin derivative in the following way:

Lemma 2.15 Faà di Bruno’s formula for Malliavin derivative. If f and g are functions with a sufficient number of derivatives, then for a random variable $F\in {\mathbb{D}}^N\left(\left[\mathrm{0,}T\right]\right)$ and $\forall n\le N$, by theorem 2.7 and lemma 15 we have

$D_t^{n}f\left(g\left(F\right)\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{f}^{\left(k\right)}\left(g\left(F\right)\right)\cdot B_{n,k}\left(g^{\prime }\left(F\right),g^{″}\left(F\right),\cdots ,g^{n-k+1}\left(F\right)\right)D_t^{n}F,$ (18)

where $B_{n\mathrm{,}k}\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_{n-k+1}\right)$ are the incomplete exponential Bell polynomials.

2.4. Extension to Two Brownian Motions

In what follows, we work with two independent Brownian motions ${\left\{W_t\right\}}_{t\ge 0}$ and ${\left\{Z_t\right\}}_{t\ge 0}$ defined in a probability space $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}\mathrm{,}\mathbb{P}\right)$, let $\left\{{\mathcal{F}}_t^W\right\}$ and $\left\{{\mathcal{F}}_t^Z\right\}$ be the filtrations generated by the Brownian motion $W_t$ and $Z_t$ respectively. Let ${\mathcal{F}}_{t_1}^W\vee {\mathcal{F}}_{t_2}^Z\mathrm{<mo>\; :\; </mo>}=\sigma \left\{W_{s_1}\mathrm{,}{Z}_{s_2}\mathrm{,}{s}_1\le t_1\mathrm{,}{s}_2\le t_2\right\}$ be the filtration generated by two Brownian motions $W_t$ and $Z_t$. When $t_1=t_2=t$, we keep the symbol ${\mathcal{F}}_t:={\mathcal{F}}_t^W\vee {\mathcal{F}}_t^Z$ for the sigma-algebra generated by both Brownian motions.

Let $D^W$ and $D^Z$ be the Malliavin derivation operator w.r.t the Brownian motion $W_t$ and $Z_t$, this implies that for a ${\mathcal{F}}_T$ measurable random variable $F\left(\omega \right)$, the 2-dimensional directional derivative of F at the point $\omega \in \Omega$ in the direction $\gamma \left({\gamma }_1\mathrm{,}{\gamma }_2\right)\in \Omega$ by:

$\begin{array}{c}{D}_{\gamma }F\left(\omega \right):=\underset{ϵ\to 0}{\lim }\frac{F\left(\omega +ϵ\gamma \right)-F\left(w\right)}{ϵ}\\ ={\displaystyle {\int }_0^T}{D}_s^{W}F\left(\omega \right)\frac{\text{d}{\gamma }_1}{\text{d}s}\text{d}s+{\displaystyle {\int }_0^T}{D}_s^{Z}F\left(\omega \right)\frac{\text{d}{\gamma }_2}{\text{d}s}\text{d}s.\end{array}$ (19)

The freezing operators ${\omega }_W^t$ and ${\omega }_Z^t$ follow the same definition 2.9 as the one dimension case. However, each random variable are depend on the the path of single Brownian motion indicated by its subscript.

3. Preliminaries on Option Pricing

Throughout this paper we shall operate in the context of a complete financial market. Options are an example of a broader class of assets called contingent claims. We will study European call option pricing under stochastic framework. The aim of this section is to review the basic objects, ideas and results of the classical Black-Scholes theory, stochastic volatility models of derivative pricing [11] .

Definition 3.1

1) A contingent claim is any asset whose future payoff is contingent on the outcome of some uncertain event.

2) A European call option is a contract that gives its holder the right, but not the obligation, to buy one unit of an underlying asset for predetermined strike price K on the maturity date T.

3.1. The Black-Scholes Theory

The Black-Scholes model is widely used for the dynamics of a financial market containing derivative investment instruments. From the Black-Scholes equation, one can deduce the Black-Scholes formula, which gives a theoretical estimate of the price of European-style options. The Black-Scholes model with constant volatility under risk-neutral probability measure is that the stock price $S_t$ satisfies the following stochastic differential equation:

$\text{d}{S}_t=r{S}_t\text{d}t+\sigma S_t\text{d}{W}_t,$ (20)

where r and $\sigma$ are constants. For reasons of convenience, we make the change of variable in the following sections, let $X_t=\ln S_t$ denote the logarithm of stock price, then

$\text{d}{X}_t=r-\frac{1}{2}{\sigma }^2\text{d}t+\sigma \text{d}{W}_t,$ (21)

the price $V_t$ of an European call option with payoff ${\left(X_T-K\right)}_{+}$ at time t for this model with constant volatility $\sigma$, current stock price ${\text{e}}^x$, maturity time T and interest rate r, satisfy the risk-neutral pricing formula [12] :

$V_t={\text{e}}^{-r\left(T-t\right)}E\left[{\left(S_T-K\right)}^{+}|{\mathcal{F}}_t\right].$ (22)

And the closed-form solution of Black-Scholes PDE is the Black-Scholes-Merton formula:

$V_t=BS\left(t,x,\sigma \right):={\text{e}}^{x}N\left(d_{+}\right)-K{\text{e}}^{-r\left(T-t\right)}N\left(d_{-}\right),$ (23)

where

$d_{\pm }\left(t,\sigma \right)=\frac{X_t-\ln K+\left(r\pm \frac{{\sigma }^2}{2}\right)\left(T-t\right)}{\sigma \sqrt{T-t}},$ (24)

and

$N\left(x\right)=\frac{1}{\sqrt{2\text{π}}}{\displaystyle {\int }_{-\infty }^x}{\text{e}}^{\frac{-y^2}{2}}\text{d}y=\frac{1}{\sqrt{2\text{π}}}{\displaystyle {\int }_{-x}^{\infty }}{\text{e}}^{\frac{-y^2}{2}}\text{d}y.$ (25)

is the standard normal cumulative distribution function. The derivation consists of finding a self-financing investment strategy, that replicates the call option payoff structure and assume that one continuously adjusts the replicating portfolio over time.

3.2. Stochastic Volatility Models

That it might make sense to model volatility as a random variable should be clear to the most casual observer of equity markets. Nevertheless, given the success of the Black-Scholes model in parsimoniously describing market option prices, it’s not immediately obvious what the benefit of making such a modeling choice might be.

3.2.1. SABR (Stochastic Alpha Beta Rho) Model with $\beta =1$

Stochastic volatility models are useful because they explain in a self-consistent way why options with different strike and expiration have different Black-Scholes implied volatility. And moreover, stochastic volatility models assume realistic dynamics for the underlying. Specifically, the SABR model is an extension of the Black Scholes model in which the volatility parameter follows a stochastic process:

$\text{d}{S}_t=r{S}_t\text{d}t+{\sigma }_t{S}_t^{\beta }\left(\rho \text{d}{W}_t+\sqrt{1-{\rho }^2}\text{d}{Z}_t\right),$ (26)

$\text{d}{\sigma }_t=\alpha {\sigma }_t\text{d}{W}_t.$ (27)

The two Brownian motions, $W_t$ and $Z_t$ are independent. It can be shown by Lévy’s Theorem that $M_t:=\rho \text{d}{W}_t+\sqrt{1-{\rho }^2}\text{d}{Z}_t$ is a Brownian motion, thus $\text{d}{M}_t\text{d}{W}_t=\rho \text{d}t$. Volatility does note mean revert in the SABR model, so it is only good for short expirations. Nevertheless the model has the virtue of having an exact expression for the implied volatility smile in the short-expiration limit $\tau :=T-t\to 0$. The resulting functional form can be used to fit observed short-dated implied volatilities and the model parameters $\alpha \mathrm{,}\beta$ and $\rho$ thereby extracted.

Hagan et al. derived, with perturbation techniques, an approximating direct formula for this implied volatility under the SABR model in [3] :

$\begin{array}{l}{\sigma }_{BS}\left(S_0\mathrm{,}K\right)\\ =\frac{{\sigma }_0}{{\left(S_{0}K\right)}^{\left(1-\beta \right)/2}\left[1+\frac{{\left(1-\beta \right)}^2}{24}{{\displaystyle \ln }}^2\frac{S_0}{K}+\frac{{\left(1-\beta \right)}^4}{1920}{{\displaystyle \ln }}^4\frac{S_0}{K}+\cdots \right]}\frac{z}{x\left(z\right)}\\ \cdot \left[1+\left(\frac{{\left(1-\beta \right)}^2}{24}\frac{{\sigma }_0^2}{{\left(S_{0}K\right)}^{1-\beta }}+\frac{1}{4}\frac{\rho \beta \alpha {\sigma }_0}{{\left(S_{0}K\right)}^{\left(1-\beta \right)/2}}+\frac{2-3{\rho }^2}{24}{\alpha }^2\right)\tau +O\left({\tau }^2\right)\right]\mathrm{,}\end{array}$ (28)

where $z:=-\frac{\alpha }{{\sigma }_0}{\left(S_{0}K\right)}^{\left(1-\beta \right)/2}\log \left(\frac{S_0}{K}\right)$ and $x\left(z\right)=\ln \left(\frac{\sqrt{1-2\rho z+z^2}+z-\rho }{1-\rho }\right)$.

For the case of at-the money options, i.e. when $S_0=K$, this formula reduces to

${\sigma }_{BS}\left(S_0\mathrm{,}{S}_0\right)=\frac{{\sigma }_0}{S_0^{1-\beta }}\cdot \left[1+\left(\frac{{\left(1-\beta \right)}^2{\sigma }_0^2}{24S_0^{2-2\beta }}+\frac{\rho \beta \alpha {\sigma }_0}{4S_0^{1-\beta }}+\frac{2-3{\rho }^2}{24}{\alpha }^2\right)\tau +O\left({\tau }^2\right)\right]\mathrm{.}$ (29)

In the special case $\beta =1$, the SABR implied volatility formula reduces to

${\sigma }_{BS}\left(S_0\mathrm{,}K\right)={\sigma }_0\frac{y}{f\left(y\right)}\left[1+\left(\frac{1}{4}\rho \alpha {\sigma }_0+\frac{2-3{\rho }^2}{24}{\alpha }^2\right)\tau +O\left({\tau }^2\right)\right]\mathrm{,}$ (30)

where $y:=-\frac{\alpha }{{\sigma }_0}\log \left(\frac{S_0}{K}\right)$ and $f\left(y\right)=\ln \left(\frac{\sqrt{1-2\rho y+y^2}+y-\rho }{1-\rho }\right)$.

3.2.2. Exponential Functions of Brownian Motion

Marc Yor’s discovery (1992) of an integral formula for joint density of the distribution of a Brownian motion and the integral of exponential Brownian motion taken over a finite time interval has been computed in the case $\sigma =2$.

Proposition 3.2 Marc Yor’s formula. Applying Brownian motion rescaling [13] , this joint density of $\left({\displaystyle {\int }_0^t}{\text{e}}^{\sigma W_s}\text{d}s,W_t\right)$,$\sigma >0$ can be written for an arbitrary volatility parameter $\sigma$ as

$\begin{array}{c}{\varphi }_{t\mathrm{,}\sigma }\left(x\mathrm{,}y\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{\text{d}x\text{d}y}\mathbb{P}\left({\displaystyle {\int }_0^t}{\text{e}}^{\sigma W_s}\text{d}s\in \text{d}x\mathrm{,}{W}_t\in \text{d}y\right)\\ =\frac{\sigma }{2x}{\text{e}}^{-\frac{2}{{\sigma }^{2}x}\left(1+{\text{e}}^{\sigma y}\right)}\cdot \theta \left(\frac{4{\text{e}}^{\sigma y/2}}{{\sigma }^{2}x}\mathrm{,}\frac{{\sigma }^{2}t}{4}\right)\mathrm{,}\end{array}$ (31)

for $x>0,y\in R,t>0$, where

$\theta \left(r,t\right)=\frac{r}{\sqrt{2{\text{π}}^{3}t}}{\text{e}}^{\frac{{\text{π}}^2}{2t}}{\displaystyle {\int }_0^{\infty }}{\text{e}}^{-\frac{{\xi }^2}{2t}}\cdot {\text{e}}^{-r\cosh \xi }\sinh \xi \sin \frac{\text{π}\xi }{t}\text{d}\xi ,r,t>0.$ (32)

By Lyasoff [14] , (32) is equivalent to the following:

$\theta \left(r,t\right)=\frac{r}{\sqrt{2{\text{π}}^{3}t}}{\text{e}}^{\frac{{\text{π}}^2}{2t}}{\displaystyle {\int }_0^{\infty }}{\text{e}}^{-\frac{{\xi }^2}{2t}}\cdot \cosh \xi \cos \left(r\sinh \xi -\frac{\text{π}\xi }{2t}\right)\text{d}\xi ,r,t>0.$ (33)

From computational point of view, the $\frac{\text{π}}{2}$ -formula: (31) with $\theta (\cdot )$ defined as (33), may be preferable to the $\text{π}$ -formula: (31) with $\theta (\cdot )$ defined as (32).

Proposition 3.3 A straightforward application of the Cameron-Martin-Girsanov theorem implies that the joint density of $\left({\displaystyle {\int }_0^t}{\text{e}}^{\sigma W_s-\mu s}\text{d}s\mathrm{,}{W}_t\right)$,$\sigma >0,\mu \in ℝ$, which we denote by ${\varphi }_{t\mathrm{,}\sigma \mathrm{,}\mu }\left(x\mathrm{,}y\right)$,$x>0,y\in ℝ$, can be connected with the density ${\varphi }_{t,\sigma }\left(x,y\right)={\varphi }_{t,\sigma ,0}\left(x,y\right)$ through the formula

${\varphi }_{t,\sigma ,\mu }\left(x,y\right)={\text{e}}^{-\frac{\mu }{\sigma }y+\frac{{\mu }^{2}t}{2{\sigma }^2}}{\varphi }_{t,\sigma ,0}\left(x,y-\frac{\mu }{\sigma }t\right).$ (34)

4. Hull and White Formula and Extension

The no-arbitrage price at time t using the risk-neutral theory for any derivatives with terminal time T and payoff function $h\left(x\right)$ is given by the risk-neutral formula below:

$V_t=E\left[{\text{e}}^{-r\left(T-t\right)}h\left(X_T\right)\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]\mathrm{.}$ (35)

Thus $V_t$ is a no-arbitrage price for the contingent claim. In what follows, we consider the pricing of a call option, i.e.:

$h\left(X_T\right)={\left({\text{e}}^{X_T}-K\right)}^{+}.$ (36)

4.1. Hull-White Formula: Uncorrelated Volatility

Under the assumption that the volatility ${\sigma }_t$ is uncorrelated with the asset price driven by another Brownian motion $Z_t$, i.e. when $\rho =0$, the pricing formula (35) can be further simplified. By conditioning on the path of the volatility process and using the iterated conditioning property, the European call option price is given by

$V\left(t\mathrm{,}{X}_t\mathrm{,}{\sigma }_t\right)={\text{e}}^{-r\left(T-t\right)}E\left[E\left[{\left({\text{e}}^{X_T}-K\right)}^{+}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\vee {\mathcal{F}}_T^W\right]\right].$ (37)

The inner expectation is the Black-Scholes computation with a time-dependent volatility. Since ${\sigma }_t$ is a Markov process, we can apply the Black-Scholes formula, and obtain:

$V\left(t\mathrm{,}x\mathrm{,}y\right)=E\left[BS\left(t\mathrm{,}x\mathrm{<mo>\; ;\; </mo>}K\mathrm{,}T\mathrm{<mo>\; ;\; </mo>}{v}_t\right)\mathrm{<mo>\; |\; </mo>}{Y}_t=y\right]\mathrm{,}$ (38)

where

$v_t^2=\frac{1}{T-t}{\displaystyle {\int }_t^T}{\sigma }_s^2\text{d}s\mathrm{,}$ (39)

is the root-mean-square time future average volatility.

4.2. Hull-White Formula: Correlated Volatility

In general, the situation is more complicated when volatility is correlated with the Brownian motion $W_t$ driving the stock price. Again we can use iterated expectation to price a European call option.

$V\left(t,x,y\right)=E\left[{\xi }_{t}BS\left(t,x;K{\xi }_t^{-1},T;{\bar{\sigma }}_{\rho }\right)|Y_t=y\right],$ (40)

where

${\xi }_t=\exp \left(\rho {\displaystyle {\int }_t^T}{\sigma }_s\text{d}{\hat{Z}}_s-\frac{1}{2}{\rho }^2{\displaystyle {\int }_t^T}{\sigma }_s^2\text{d}s\right),$ (41)

${\bar{\sigma }}_{\rho }^2=\frac{1}{T-t}{\displaystyle {\int }_t^T}\left(1-{\rho }^2\right){\sigma }_s^2\text{d}s.$ (42)

The Hull-White formula is of practical use for Monte Carlo simulation of prices in a correlated stochastic volatility model since only one Brownian motion path has to be generated. However, it does not directly reveal any information about the implied volatility curve like the uncorrelated case.

4.3. A Generalization of Hull-White Formula

The classical Hull-White formula for option pricing can be extended, by means of Malliavin Calculus, to the correlated case. The main problem is that average future volatility is not adapted, however, this issue can be resolved by anticipating stochastic calculus. And this method decomposes option prices as the sum of the same derivative price if there is no correlation and a correction due by correlation. The following theorem is due to Alòs et al. (2006).

Theorem 4.1 Consider model (26)-(27) with $\beta =1$, and assume the following hypotheses hold:

1) The payoff function $h\mathrm{<mo>\; :\; </mo>}ℝ\to {ℝ}_{+}$ is continuous and piecewise ${\mathcal{C}}^1$ ;

2) There exists a positive real constant a such that $a\le {\sigma }_t^2$ for all $t\in \left[\mathrm{0,}T\right]$ ;

3) ${\sigma }^2\in {\mathcal{L}}_W^{\mathrm{1,2}}\left(\left[\mathrm{0,}T\right]\right)$ ;

4) For all $t\in \left[\mathrm{0,}T\right]$ there exists a positive constant C such that for all $s\in \left[\mathrm{0,}T\right]$,

$\left|E\left[\left({\displaystyle {\int }_s^T}{D}_s^W{\sigma }_r^2\text{d}r\right){\sigma }_s\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]\right|\le C\mathrm{.}$ (43)

Then, for all $t\in \left[\mathrm{0,}T\right]$,

$V_t=E\left[BS\left(t\mathrm{,}{X}_t\mathrm{,}{v}_t\right)\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]+\frac{\rho }{2}E\left[{\displaystyle {\int }_t^T}{\text{e}}^{-r\left(s-t\right)}H\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right){\Lambda }_s\text{d}s\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]\mathrm{,}$ (44)

where $v_t^2$ is the future average volatility defined (39) in Subsection 4.1 and

$H\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)\mathrm{<mo>\; :\; </mo>}=\left(\frac{{\partial }^3}{\partial x^3}-\frac{{\partial }^2}{\partial x^2}\right)BS\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)\mathrm{,}$ (45)

${\Lambda }_s\mathrm{<mo>\; :\; </mo>}=\left({\displaystyle {\int }_s^T}{D}_s^W{\sigma }_r^2\text{d}r\right){\sigma }_s\mathrm{.}$ (46)

Notice that formula (44) does not reduce the dimensionality of the problem but identifies the impact of correlation. When $\rho =0$, it is the same as (38).

5. Dyson Series in the Return’s Idiosyncratic Noise

5.1. Application of Marc Yor’s Formula

Throughout this and next section we denote by $V_{t_1,t_2}:={\displaystyle {\int }_{t_1}^{t_2}}{\sigma }_u^2\text{d}u$ a cumulative time integral, from $t_1$ to $t_2$, of future volatility, i.e. $V_{s,T}=v_s^2\left(T-s\right)$. The aim of this paper is to extend Theorem 4.1 to a deterministic form by specifically assuming the underlying asset and volatility process follow (26), then ${\sigma }_t$ is a square integrable process adapted to $\left\{{\mathcal{F}}_t^W\right\}$.

Lemma 5.1 The conditional probability density function of $V_{s\mathrm{,}T}$ is $\frac{1}{{\sigma }_s^2}{\psi }_{V_{s\mathrm{,}T}}\left(\frac{v}{{\sigma }_s^2}\right)$.

where

$y_{V_{s\mathrm{,}T}}\left(v\right)={\displaystyle {\int }_{-\infty }^{\infty }}{\varphi }_{T-s\mathrm{,2}\alpha \mathrm{,}{\alpha }^2}\left(v\mathrm{,}z\right)\text{d}z\mathrm{.}$ (47)

One straightforward application of (47) is using the conditional density of $V_{t\mathrm{,}T}$ to obtain the first conditional expectation in (44):

Theorem 5.2 The conditional expectation of $BS\left(t\mathrm{,}{X}_t\mathrm{,}{v}_t\right)$ is

$E\left[BS\left(t\mathrm{,}{X}_t\mathrm{,}{v}_t\right)\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]={\displaystyle {\int }_0^{\infty }}BS\left(t\mathrm{,}{X}_t\mathrm{,}\sqrt{\frac{v}{T-t}}\right)\frac{1}{{\sigma }_t^2}{\psi }_{V_{t\mathrm{,}T}}\left(\frac{v}{{\sigma }_t^2}\right)\text{d}v\mathrm{.}$ (48)

5.2. Application of Exponential Formula

Theorem 5.3 For $t\le s$, define $G\left(s,X_s,v_s\right)={\displaystyle {\sum }_{n=0}^{\infty }}\frac{1}{2^{n}n!}{g}_n\left(s,X_t,v_t\right)$, where $g_n\left(s,X_t,v_t\right)={\omega }_Z^t\circ {\displaystyle {\int }_{{\left[t,T\right]}^n}}{D}_{\tau \otimes n}^{2n,Z}H\left(s,X_s,v_s\right)\text{d}{\tau }^{\otimes n}$. Let $G_s\mathrm{,}{H}_s$ be the short notation for $G\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$ and $H\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$, then the option price (44) can be further simplified as the following:

$V_t=E\left[BS\left(t,X_t,v_t\right)|{\mathcal{F}}_t\right]+\frac{\rho }{2}{\displaystyle {\int }_t^T}{\text{e}}^{-r\left(s-t\right)}E\left[{\Lambda }_s{G}_s|{\mathcal{F}}_t\right]\text{d}s.$ (49)

here $\otimes$ denotes the tensor power. In general, the tensor product $f\otimes g$ of two functions $f\mathrm{,}g$ is defined by $f\otimes g\left(x_1\mathrm{,}{x}_2\right)=f\left(x_1\right)g\left(x_2\right)$. See [15] for more detail.

Now we use Malliavin calculus to deduce a full asymptotic series for $G\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$ and use it to obtain $E\left[{\Lambda }_s{G}_s\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]$ in (49), which gives us a deterministic formula for European option price. By (21) and (23),

$\begin{array}{c}{X}_s=X_t+{\displaystyle {\int }_t^s}r-\frac{1}{2}{\sigma }_u^2\text{d}u+{\displaystyle {\int }_t^s}{\sigma }_u\left(\rho \text{d}{W}_u+\sqrt{1-{\rho }^2}\text{d}{Z}_u\right)\\ =X_t+r\left(s-t\right)-\frac{1}{2}{V}_{t\mathrm{,}s}+\frac{\rho }{\alpha }\left({\sigma }_s-{\sigma }_t\right)+\sqrt{1-{\rho }^2}{\displaystyle {\int }_t^s}{\sigma }_u\text{d}{Z}_u\mathrm{.}\end{array}$ (50)

Note that the volatility process in (26) is the differential notation for ${\sigma }_t-{\sigma }_0={\displaystyle {\int }_0^t}\alpha {\sigma }_u\text{d}{W}_u$, and obviously for $0<t<s$,${\sigma }_s-{\sigma }_t={\displaystyle {\int }_t^s}\alpha {\sigma }_u\text{d}{W}_u$. And accordingly by (50) and (24), $d_{\pm }\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)=\frac{X_s-\ln K+\left(r\pm \frac{v_s^2}{2}\right)\left(T-s\right)}{\sqrt{V_{s\mathrm{,}T}}}$. Therefore by Theorem 2.7,

$D_{\tau }^Z{X}_s=D_{\tau }^Z{\displaystyle {\int }_t^s}{\sigma }_u\sqrt{1-{\rho }^2}\text{d}{Z}_u={\sigma }_{\tau }\sqrt{1-{\rho }^2}1{\text{l}}_{\left\{\tau \le s\right\}}\mathrm{,}$ (51)

$D_{\tau }^Z{d}_{+}=D_{\tau }^Z{d}_{-}=\frac{D_{\tau }^Z{X}_s}{v_s\sqrt{T-s}}=\frac{{\sigma }_{\tau }\sqrt{1-{\rho }^2}1{\text{l}}_{\left\{\tau \le s\right\}}}{v_s\sqrt{T-s}}\mathrm{.}$ (52)

Lemma 5.4 Let two real-valued functions $p\left(t\mathrm{,}x\mathrm{,}\sigma \right)$ and $q\left(t\mathrm{,}x\mathrm{,}\sigma \right)$ be defined as following:

$p\left(t\mathrm{,}x\mathrm{,}\sigma \right)=x-\frac{d_{+}^2\left(x\mathrm{,}t\right)}{2}+\ln \left(-d_{-}\left(x\mathrm{,}t\right)\right)\mathrm{,}$ (53)

$q\left(t,x,\sigma \right)=\frac{1}{\sqrt{2\text{π}}{\sigma }^2\left(T-t\right)}{\text{e}}^x.$ (54)

Then the 2n^th order Malliavin derivative of $H_s$ can be expressed as:

$D_{\tau \otimes n}^{2n\mathrm{,}Z}{H}_s={\left(1-{\rho }^2\right)}^n{H}_s{B}_{2n}\left(p^{\prime }(\cdot )\mathrm{,}{p}^{″}(\cdot )\mathrm{,}\cdots \mathrm{,}{p}^{\left(2n\right)}(\cdot )\right){\displaystyle {\prod }_{i=1}^n{\sigma }_{{\tau }_i}^2{\text{1l}}_{\left\{{\tau }_i\le s\right\}},}$ (55)

where

$p^{\left(j\right)}\left(s,X_s,v_s\right)=(\begin{array}{l}\frac{{\left(-1\right)}^{j+1}-d_{-}^2}{{\left(\sqrt{V_{s,T}}{d}_{-}\right)}^j}\text{when}j=1,2;\\ \frac{{\left(-1\right)}^{j+1}\left(j-1\right)!}{{\left(\sqrt{V_{s,T}}{d}_{-}\right)}^j}\text{for}j\ge 3.\end{array}$ (56)

for $d_{-}$ evaluated at $\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$.

The second step to calculate $G\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$ is to apply freezing operator ${\omega }_Z^t$ to $D_{\tau \otimes n}^{2n\mathrm{,}Z}{H}_s$ for $n=0,1,2,\cdots$. Let $\mathcal{X}$ be any random variable depend on Brownian motion ${\left\{Z_t\right\}}_{t\ge 0}$, denote ${\mathcal{X}}_Z^{\omega }\mathrm{<mo>\; :\; </mo>}={\omega }_Z^t\circ \mathcal{X}$ be the random variable $\mathcal{X}$ applied by the freezing operator ${\omega }_Z^t$, by Proposition 2.10,

$\begin{array}{c}{X}_s^{\omega }\mathrm{<mo>\; :\; </mo>}=X_t+r\left(s-t\right)-\frac{1}{2}{\displaystyle {\int }_t^s}{\sigma }_u^2\text{d}u+\frac{\rho }{\alpha }\left({\sigma }_s-{\sigma }_t\right)+{\omega }_Z^t\circ {\displaystyle {\int }_t^s}{\sigma }_u\sqrt{1-{\rho }^2}\text{d}{Z}_u\\ =X_t+r\left(s-t\right)-\frac{1}{2}{V}_{t\mathrm{,}s}+\frac{\rho }{\alpha }\left({\sigma }_s-{\sigma }_t\right)\mathrm{,}\end{array}$ (57)

and accordingly we have

$d_{\pm }^{\omega }\left(s,X_s,v_s\right)=d_{\pm }\left(s,X_s^{\omega },v_s\right)=\frac{X_s^{\omega }-\ln K+\left(r\pm \frac{v_s^2}{2}\right)\left(T-s\right)}{\sqrt{V_{s,T}}},$ (58)

$H_s^{\omega }=\frac{\sqrt{V_{s\mathrm{,}T}}-d_{+}^{\omega }}{\sqrt{2\text{π}}{V}_{s\mathrm{,}T}}{\text{e}}^{X_s^{\omega }-\frac{d_{+}^{\omega }{}^2}{2}}=\frac{-d_{-}^{\omega }}{\sqrt{2\text{π}}{V}_{s\mathrm{,}T}}{\text{e}}^{X_s^{\omega }-\frac{d_{+}^{\omega }{}^2}{2}}\mathrm{.}$ (59)

Therefore in general, for $n=0,1,2,\cdots$,

$\begin{array}{l}{D}_{\tau \otimes n}^{2n\mathrm{,}Z}{H}_s^{\omega }\\ ={\omega }_Z^t\circ \left[{\left(1-{\rho }^2\right)}^n{H}_s{B}_{2n}\left(p^{\prime }\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)\mathrm{,}\cdots \mathrm{,}{p}^{\left(2n\right)}\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)\right){\displaystyle {\prod }_{i=1}^n{\sigma }_{{\tau }_i}^2}\right]\\ ={\left(1-{\rho }^2\right)}^n{H}_s^{\omega }{B}_{2n}\left(p^{\prime }\left(s\mathrm{,}{X}_s^{\omega }\mathrm{,}{v}_s\right)\mathrm{,}\cdots \mathrm{,}{p}^{\left(2n\right)}\left(s\mathrm{,}{X}_s^{\omega }\mathrm{,}{v}_s\right)\right){\displaystyle {\prod }_{i=1}^n{\sigma }_{{\tau }_i}^2{\text{1l}}_{\left\{{\tau }_i\le s\right\}}\mathrm{,}}\end{array}$ (60)

where $p^{\left(j\right)}\left(s\mathrm{,}{X}_s^{\omega }\mathrm{,}{v}_s\right)\mathrm{<mo>\; :\; </mo>}={\omega }_Z^t\circ p^{\left(j\right)}\left(s\mathrm{,}{X}_s\mathrm{,}{v}_s\right)$ is given by (56) except that $d_{-}$ is now evaluated at $\left(s\mathrm{,}{X}_s^{\omega }\mathrm{,}{v}_s\right)$. Thus, by (60), we are able to compute $G_s$ in the following:

$\begin{array}{c}{G}_s={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}\frac{1}{2^{n}n!}{\omega }_Z^t\circ {\displaystyle {\int }_{{\left[t,T\right]}^n}}{D}_{\tau \otimes n}^{2n,Z}H\left(s,X_s,v_s\right)\text{d}{\tau }^{\otimes n}\\ ={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}\frac{{\left(1-{\rho }^2\right)}^n}{2^{n}n!}{H}_s^{\omega }{B}_{2n}\left(p^{\prime }\left(s,X_s,v_s\right),\cdots ,p^{\left(2n\right)}\left(s,X_s,v_s\right)\right)\\ \cdot {\displaystyle {\int }_{{\left[t,T\right]}^n}}{\displaystyle {\prod }_{i=1}^n{\sigma }_{{\tau }_i}^2{\text{1l}}_{\left\{{\tau }_i\le s\right\}}\text{d}{\tau }^{\otimes n}}\\ =H_s^{\omega }{\displaystyle \underset{n=0}{\overset{\infty }{\sum }}}\frac{{\left(1-{\rho }^2\right)}^n}{2^{n}n!}{\left(V_{t,s}\right)}^n{B}_{2n}\left(p^{\prime }\left(s,X_s^{\omega },v_s\right),\cdots ,p^{\left(2n\right)}\left(s,X_s^{\omega },v_s\right)\right).\end{array}$ (61)

5.3. Option Pricing Formula for SABR Model

Lemma 5.5 Let $L_s^{\omega }=V_{s,T}{H}_s^{\omega }=\frac{-d_{-}^{\omega }}{\sqrt{2\text{π}}}{\text{e}}^{X_s^{\omega }-\frac{d_{+}^{\omega }{}^2}{2}}$ and