Efficient Asian option pricing under regime switching jump diffusions and stochastic volatility models

Abstract

Utilizing frame duality and a FFT-based implementation of density projection we develop a novel and efficient transform method to price Asian options for very general asset dynamics, including regime switching Lévy processes and other jump diffusions as well as stochastic volatility models with jumps. The method combines continuous-time Markov chain approximation, with Fourier pricing techniques. In particular, our method encompasses Heston, Hull-White, Stein-Stein, 3/2 model as well as recently proposed Jacobi, $\alpha$-Hypergeometric, and 4/2 models, for virtually any type of jump amplitude distribution in the return process. This framework thus provides a ‘unified’ approach to pricing Asian options in stochastic jump diffusion models and is readily extended to alternative exotic contracts. We also derive a characteristic function recursion by generalizing the Carverhill-Clewlow factorization which enables the application of transform methods in general. Numerical results are provided to illustrate the effectiveness of the method. Various extensions of this method have since been developed, including the pricing of barrier, American, and realized variance derivatives.

Appendices

Implementation details

1.1 Nonuniform variance grid

Table 15 Stable coefficient approximations derived for the cubic basis using a five point Newton-Cotes quadrature (Boole’s rule)

To apply the Markov chain approximation, we determine the variance grid $\{v_j\}_{j=1}^{m_0}$ as follows. First, we fix $t = T/2$ and center the grid about the mean of the variance process $v_t$ by: $v_1: = \max \{{{\bar{v}}}, {\bar{\mu }}(t) - \gamma {\bar{\sigma }}(t)\}$ if the domain of $v_t$ is positive; otherwise $v_1= {\bar{\mu }}(t) - \gamma {\bar{\sigma }}(t)$. We next choose $v_{m_0}:= {\bar{\mu }}(t) + \gamma \bar{\sigma }(t)$. Here we have defined ${\bar{\mu }}(t) = {\mathbb {E}}[v_t |v_0]$ and ${\bar{\sigma }}(t)$ the standard deviation conditional on $v_0$. Further, we define the constant⁷$\gamma = 4.5$ and ${{\bar{v}}} = 0.00001$.⁸ Finally, we generate $v_2,v_3,\ldots ,v_{m_0-1}$ using the procedure

$$\begin{aligned} v_i=v_0+{\bar{\alpha }}\sinh \left( c_2\frac{i}{m_0}+c_1\Big (1-\frac{i}{m_0}\Big ) \right) \end{aligned}$$

where

$$\begin{aligned} c_1={{\,\mathrm{arcsinh}\,}}\left( \frac{v_1-v_0}{{\bar{\alpha }}} \right) , \quad c_2= {{\,\mathrm{arcsinh}\,}}\left( \frac{v_{m_0}-v_0}{{\bar{\alpha }}} \right) \end{aligned}$$

for ${\bar{\alpha }}<(v_{m_0}-v_1)$. This transformation concentrates more grid points near the critical point $v_0$, where the non-uniformity of the grid is determined by the parameter ${\bar{\alpha }}$: smaller ${\bar{\alpha }}$ results in a more nonuniform grid. For computations in this paper we choose ${\bar{\alpha }}= (v_{m_0} - v_1)/5$. Since $v_0$ is not likely a member of the variance grid, find the bracketing index $j_0$ such that $v_{j_0}\le v_0 < v_{j_0+1}$. Holding the points $v_1, v_2$ constant,⁹ we shift the remaining points $v_j, j\ge 2$ by $v_0 -v_{j_0}$ so that $v_{j_0}= v_0$ is a member of the adjusted grid.

Remark 3

We note that for Heston, 3/2 (with reformulated parameters), and 4/2 models

$$\begin{aligned} {\bar{\mu }}(t)= & {} e^{-\eta t}v_0 + \theta (1- e^{-\eta t}), \qquad {\bar{\sigma }}^2(t):=\frac{\sigma _v^2}{\eta }v_0(e^{-\eta t} - e^{-2\eta t}) \\&+ \frac{\theta \sigma ^2_v}{2\eta }(1- e^{-\eta t}+e^{-2\eta t}). \end{aligned}$$

In the Stein-Stein’s model, we have

$$\begin{aligned} {\bar{\mu }}(t) = e^{-\eta t}v_0 + \theta (1- e^{-\eta t}), \qquad {\bar{\sigma }}^2(t):=\frac{\sigma _v^2}{2\eta }(1- e^{-2\eta t}). \end{aligned}$$

For Hull-White’s model, the variance process is described by the moments

$$\begin{aligned} {\bar{\mu }}(t) = v_0 e^{a_vt}, \qquad {\bar{\sigma }}^2(t) = v_0^2 e^{2a_vt}(e^{\sigma _v^2t} -1). \end{aligned}$$

1.2 Approximation of $\Psi$

To approximate ${\bar{\Psi }} =\{{\bar{\Psi }}\}_{n,k}$ for $n =1,\ldots , N$, $k = 1,\ldots , N +N_{M-1}$, we apply Gaussian quadrature over each interval $I_k^m =[x^*_{k}-2\Delta ,x^*_{k}+2\Delta ]$, where $x_k^*$ is defined in Sect. 4.2. After a change of variables, this is accomplished by applying an $N_q$-point quadrature applied to each subinterval of $[-2,2]=[-2,-1]\cup [-1,0]\cup [0,1]\cup [1,2]$ as follows. Specifically,

$$\begin{aligned} \Psi (n,k)&: =a^{1/2} \int _{I^m_k}\left( e^y + 1\right) ^{i\xi _n}a^{1/2}\varphi (a(y-x_k^*))dy \\&= \int _{[-2,2]}\exp \left( \mathrm {i}\xi _n \ln \left( 1+\exp \left( x^*_k +\frac{y}{a}\right) \right) \right) \varphi (y)dy \\&\approx \sum _{l=1}^{4\cdot N_q} \omega _l\cdot \exp \left( \mathrm {i}\xi _n \ln \left( 1+\exp \left( x^*_k +\gamma _l\right) \right) \right) \varphi (\gamma _l), \end{aligned}$$

where $\{(\gamma _l, \omega _l)\}_{l=1}^{4\cdot N_q}$ are the nodes and weights. If we define the sample grid

$$\begin{aligned} \eta _j, \quad j=1,...,N_\eta , \qquad N_\eta :=(N +N_{M-1} -1)\cdot N_q + 4 \cdot N_q, \end{aligned}$$

(49)

the Gaussian approximation of $\Psi (n,k)$, for $n=1,\ldots ,N$ and $k=1,\ldots ,N+N_{M-1}$, is given by

$$\begin{aligned} {\bar{\Psi }} (n,k)&:= \sum _{l=1}^{4\cdot N_q} \theta ^n_{N_q(k-1)+l}\cdot \sigma _l = \sum _{l=1}^{2\cdot N_q} \left( \theta ^n_{N_q(k-1)+l} + \theta ^n_{N_q(k+3)+1-l}\right) \cdot \sigma _l , \end{aligned}$$

(50)

where $\sigma _l := \varphi ^{[3]}(\gamma _l)\cdot \omega _l,$$l = 1,..., 4 \cdot N_q$, and

$$\begin{aligned} \theta _j^n := \exp \left( \mathrm {i}\xi _n \ln \left( 1+\exp \left( \eta _j\right) \right) \right) , \qquad n = 1,...,N, \quad j=1,...,N_\eta . \end{aligned}$$

To populate ${\bar{\Psi }} (n,k)$ efficiently, note that $\theta _j^n = \theta _j^{n-1}\cdot \exp (\mathrm {i}\Delta _\xi \ln (1+\exp (\eta _j)))$. In this work, we utilize a five-point Gaussian quadrature, detailed in Appendix A.2.1.

1.2.1 Five-point Gaussian quadrature

A composite five point Gaussian quadrature is implemented by applying a five point rule over each interval $[r,r+1]$, for $r = -2,-1,0,1$. By symmetry we only evaluate for $r=-2,-1$. On $[-2,-1]$, we obtain the nodes and weights

$$\begin{aligned} \{\gamma _l\}_{l=1}^5= & {} \Big \{-\frac{3}{2}-g_3,-\frac{3}{2}- g_2, -\frac{3}{2},-\frac{3}{2}+g_2,-\frac{3}{2}+ g_3\Big \}, \\ \{\omega _l\}_{l=1}^5= & {} \frac{1}{2}\{ {{\hat{v}}}_3,{{\hat{v}}}_2,{{\hat{v}}}_1,{{\hat{v}}}_2,{{\hat{v}}}_3\}, \end{aligned}$$

while on $[-1,0]$ we have

$$\begin{aligned} \{\gamma _l\}_{l=6}^{10}= & {} \Big \{-\frac{1}{2}-g_3,-\frac{1}{2}- g_2, -\frac{1}{2},-\frac{1}{2}+g_2,-\frac{1}{2}+ g_3\Big \}, \\ \{\omega _l\}_{l=6}^{10}= & {} \frac{1}{2}\{{{\hat{v}}}_3,{{\hat{v}}}_2,{{\hat{v}}}_1, {{\hat{v}}}_2,{{\hat{v}}}_3\}, \end{aligned}$$

where we define the constants

$$\begin{aligned}&g_2: = \frac{1}{6}\sqrt{5 - 2 \sqrt{10/7}}, \quad g_3 : =\frac{1}{6}\sqrt{5 + 2 \sqrt{10/7}}\\&\quad {{\hat{v}}}_1: = 128/225, \quad {{\hat{v}}}_2: = (322 + 13\sqrt{70})/900, \quad {{\hat{v}}}_3:=(322-13\sqrt{70})/900. \end{aligned}$$

The final weights $\sigma _l = \varphi ^{[3]}(\gamma _l)\cdot w_l$ are found by evaluating the cubic generator defined in equation (38) at each of $\{\gamma _l\}_{l=1}^{10}$, and $\sigma _l$ is stored for repeated use.

1.3 Continuous monitoring

Rather than set M extremely large to estimate the option value with continuous monitoring, we can utilize the now standard technique of Richardson extrapolation from values computed with modest levels of M. Let ${\mathcal {V}}_N(M)$ denote the discretely monitored value approximation with M monitoring dates, and with N fixed. By choosing $d\in {\mathbb {N}}_+$, we can approximate the continuously monitored value with a four-point Richardson extrapolation:

$$\begin{aligned} {\mathcal {V}}_N^\infty (d) := \frac{1}{21}\left( 64{\mathcal {V}}_N(2^{d+3}) - 56 {\mathcal {V}}_N(2^{d+2})+14{\mathcal {V}}_N(2^{d+1})-{\mathcal {V}}_N(2^d) \right) . \end{aligned}$$

as applied in Zhang and Oosterlee (2013). An efficient APROJ extrapolation algorithm can be devised by reusing the matrix ${\bar{\Psi }}$ for each d. See Table 16 for an example of the continuous monitoring approximation

Table 16 Continuously monitored Asian option values by Richardson Extrapolation

Error analysis

This section demonstrates stability of the characteristic function recursion, and provides a bound on the rate of convergence of the APROJ value error.¹⁰ Returning shortly to the case of $m=2$, we have for $m\ge 3$

$$\begin{aligned} \epsilon ({\bar{\phi }}^j_{Z_{m-1}}(\xi ))&:= \phi ^j_{Z_{m-1}}(\xi )-\bar{\phi }^j_{Z_{m-1}}(\xi ) = \int _{\mathbb {R}}(e^y+1)^{\mathrm {i}\xi }(f^j_{Y_{m-1}}(y) -\bar{f}^j_{Y_{m-1}}(y))dy\\&= \int _{{\mathbb {R}}/G_m} (e^y+1)^{\mathrm {i}\xi }f^j_{Y_{m-1}}(y)dy \\&\quad + \left( \int _{G_m} (e^y+1)^{\mathrm {i}\xi }f^j_{Y_{m-1}}(y)dy - \Upsilon _{a,N}\sum _{k=1}^N\beta ^j_{Y_{m-1},k}\Psi _{m-1}(\xi ,k)\right) \\&\quad + \Upsilon _{a,N}\sum _{k=1}^N\beta ^j_{Y_{m-1},k}( \Psi _{m-1}(\xi ,k) - {\bar{\Psi }}_{m-1}(\xi ,k)) \\&\quad + \Upsilon _{a,N}\sum _{k=1}^N\bar{\Psi }_{m-1}(\xi ,k)(\beta ^j_{Y_{m-1},k}-{\bar{\beta }}^j_{Y_{m-1},k})\\&:=\left( \tau (G_{m-1}) + {\mathcal {E}}^j_{m-1,1}(\xi ) + {\mathcal {E}}^j_{m-1,2}(\xi )\right) + {\mathcal {E}}^j_{m-1}(\xi ), \end{aligned}$$

where $G_m = [{\bar{\mu }}_m - {{\bar{a}}}/2, {\bar{\mu }}_m +{{\bar{a}}}/2]$, and the final term will be further split into two components.¹¹ The term $\tau (G_{m-1})$ captures density truncation error and can be bounded by $\tau (G)$ which converges exponentially in ${{\bar{a}}} = 2^{{{\bar{P}}}}$ for most processes of interest. If we define $P_{{\mathcal {M}}_a} f^j_{Y_{m-1}}$ to be the true (untruncated) orthogonal projection of $f^j_{Y_{m-1}}$ onto the space ${\mathcal {M}}_a$, the second term satisfies

$$\begin{aligned} {\mathcal {E}}^j_{m-1,1}(\xi )&:= \int _{G_m} (e^y+1)^{\mathrm {i}\xi }f^j_{Y_{m-1}}(y)dy - \Upsilon _{a,N}\sum _{k=1}^N\beta ^j_{Y_{m-1},k}\Psi _{m-1}(\xi ,k) \\&\le \left\Vert (e^y+1)^{\mathrm {i}\xi }\right\Vert ^{G_{m-1}}_2 \cdot \left\Vert f^j_{Y_{m-1}} - P_{{\mathcal {M}}_a} f^j_{Y_{m-1}}\right\Vert ^{{\mathbb {R}}}_2 \\&\le \sqrt{{{\bar{a}}}}\cdot C_1(f^j_{Y_{m-1}})\cdot \Delta ^4, \end{aligned}$$

where $\Delta ^4$ is the theoretical convergence rate of cubic projection (in practice, the rate is much faster). The constant $C_1(f^j_{Y_{m-1}})$ is bounded by a constant multiple of $\left\Vert (-\mathrm {i}\xi )^4 \phi _{Y^j_{m-1}}(\xi )\right\Vert _2 \le \left\Vert \xi ^4\phi _{X_{\Delta _t}}^{{{\bar{k}}}}(\xi )\right\Vert _2$ for some $1\le {{\bar{k}}} \le m_0$, where the inequality follows by Lemma 1. Hence, for a constant $C_1$ we have the bound

$$\begin{aligned} {\mathcal {E}}^j_{m-1,1}(\xi ) \le C_1\cdot \sqrt{{{\bar{a}}}}\cdot \Delta ^4 \end{aligned}$$

(51)

which is independent of $\xi , j, m$. This will govern the overall rate of convergence. Next define the numerical integration error

$$\begin{aligned} \epsilon ({\bar{\Psi }}):=\sup \{|{\bar{\Psi }}(\xi _n,k) - \Psi (\xi _n,k)|: 1\le n \le N, 1\le k\le N+N_{M-1}\}. \end{aligned}$$

For a constant $C(\varphi )$, it follows from Lemma 5.2 of Kirkby (2016) that

$$\begin{aligned} {\mathcal {E}}^j_{m-1,2}(\xi )&= \Upsilon _{a,N}\sum _{k=1}^N\beta ^j_{Y_{m-1},k}(\Psi _{m-1}(\xi ,k) - {\bar{\Psi }}_{m-1}(\xi ,k)) \\&\le \sqrt{{{\bar{a}}}} \cdot \epsilon ({\bar{\Psi }})\cdot C(\varphi )\cdot \left\Vert f^j_{Y_{m-1}}\right\Vert _2 \le C_2 \cdot \sqrt{{{\bar{a}}}} \cdot \epsilon ({\bar{\Psi }}) \end{aligned}$$

for some $C_2$, where the second inequality again follows from Lemma 1.

The remaining term, which provides the link between errors though time, is

$$\begin{aligned} {\mathcal {E}}^j_{m-1}(\xi )&:= \Upsilon _{a,N}\sum _{k=1}^N\bar{\Psi }_{m-1}(\xi ,k)(\beta ^j_{Y_{m-1},k}-{\bar{\beta }}^j_{Y_{m-1},k})\\&= a^{-1/2}\sum _{k=1}^N{\bar{\Psi }}_{m-1}(\xi ,k)\cdot \epsilon ({\bar{\beta }}^j_{Y_{m-1},k}), \end{aligned}$$

where $\epsilon ({\bar{\beta }}^j_{Y_{m-1},k}):=a^{1/2} \Upsilon _{a,N}(\beta ^j_{Y_{m-1},k}-{\bar{\beta }}^j_{Y_{m-1},k})$. Since ${\bar{\beta }}^j_{Y_{m-1},k}$ contains several sources of error, we define $\breve{\beta }^j_{Y_{m-1},k}$ to be the approximation using the true $\phi ^j_{Y_{m-1}}$ rather than ${\bar{\phi }}^j_{Y_{m-1}}$. Hence

$$\begin{aligned} \epsilon ({\bar{\beta }}^j_{Y_{m-1},k})&=\left( \langle f^j_{Y_{m-1}},{\widetilde{\varphi }}_{a,k} \rangle - a^{1/2}\Upsilon _{a,N}\breve{\beta }^j_{Y_{m-1},k}\right) + a^{1/2}\Upsilon _{a,N}\left( \breve{\beta }^j_{Y_{m-1},k} - \bar{\beta }^j_{Y_{m-1},k}\right) \\&:=\epsilon _1({\bar{\beta }}^j_{Y_{m-1},k}) + \epsilon _2(\bar{\beta }^j_{Y_{m-1},k}), \end{aligned}$$

which yields ${\mathcal {E}}^j_{m-1}(\xi ) = {\mathcal {E}}^j_{m-1,3}(\xi ) + {\mathcal {E}}^j_{m-1,4}(\xi )$.

The first term can be captured by the Corollary 3.2 of Kirkby (2016), modified slightly to the present case. We assume that $\widetilde{\phi }^j_{X_{\Delta _t}}$, $j=1,\ldots ,m_0$ are analytic within a strip ${\mathcal {D}}_d:=\{z\in {\mathbb {C}}: \mathfrak {I}(z)\in (-d,d)\}$, for some $d>0$, and satisfy equation (18). We further define a constant $C_M$, which is bounded by a multiple of $\max _{1\le m\le M, 1\le j\le m_0} \left\Vert \phi ^j_{Y_m}\right\Vert ^{\mathcal H_d}$, where $\left\Vert f\right\Vert ^{{\mathcal {H}}_d}$ is the Hardy norm of f on ${\mathcal {D}}_d$.

Corollary 3

Fix $a = 2^P$ and $N = a \cdot {{\bar{a}}}$, where ${{\bar{a}}} = 2^{{{\bar{P}}}}$ for ${{\bar{P}}} >1 + \log _2|{\bar{\mu }}_M|$. Assume Lemma 1 holds for some $c_{{{\bar{k}}}},\kappa _{{{\bar{k}}}} >0$ and $\nu _{ {{\bar{k}}}} \in (0,2]$. Then for some $0<\gamma \le d$

$$\begin{aligned} \sup _{1\le k\le N}\left| a^{1/2}\Upsilon _{a,N}\cdot \breve{\beta }^j_{Y_{m},k} - \langle f^j_{Y_{m}},\widetilde{\varphi }_{a,k} \rangle \right|&\le \frac{a^{-1/2}}{\pi }\left( C_M\frac{e^{-({{\bar{a}}} - 2|\bar{\mu }_M|)\gamma /2}}{1-e^{-{{\bar{a}}} \gamma }} + \tau _{a}\left( \widetilde{\phi }_{X_{\Delta _t}}^{{{\bar{k}}}}\right) \right) \nonumber \\&:= \frac{a^{-1/2}}{\pi }\epsilon ^M(a,{{\bar{a}}}) \end{aligned}$$

(52)

independently of $1\le m \le M$ and $j=1,\ldots , m_0$ where $\tau _{a}({\widetilde{\phi }}_{X_{\Delta _t}}^{{{\bar{k}}}}) =\mathcal O(a\exp (-\Delta t c_{{{\bar{k}}}} \cdot (2\pi a)^{\nu _{{{\bar{k}}}}}))$. For large enough $a>0$, and $d<\infty$, $\gamma$ will approach d.

By Corollary 3 and $|{\bar{\Psi }}_{m-1}(\xi ,k)|\le 1$ we have

$$\begin{aligned} |{\mathcal {E}}^j_{m-1,3}(\xi )| \le a^{-1/2}\sum _{k=1}^N|\bar{\Psi }_{m-1}(\xi ,k)| \cdot |\epsilon _1({\bar{\beta }}^j_{Y_{m-1},k})| \le \frac{{{\bar{a}}}}{\pi } \epsilon ^M(a, {{\bar{a}}}). \end{aligned}$$

The discretization error represented by the first term in $\epsilon ^M(a,{{\bar{a}}})$ decays exponentially in ${{\bar{a}}}$, while the truncation error decays exponentially in a (for fixed ${{\bar{a}}}$).

The final term to estimate is

$$\begin{aligned} {\mathcal {E}}^j_{m-1,4}(\xi ) :=a^{-1/2}\sum _{k=1}^N\bar{\Psi }_{m-1}(\xi ,k) \cdot \epsilon _2({\bar{\beta }}^j_{Y_{m-1},k}) \end{aligned}$$

where, with $h^{m-1}_{a,k}(\xi ):=\widehat{\widetilde{\varphi }}(\xi /a)\exp (\mathrm {i}\xi x_k^{m-1})$,

$$\begin{aligned} \epsilon _2({\bar{\beta }}^j_{Y_{m-1},k})&= \frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}\left( \phi ^j_{Y_{m-1}}(\xi _n) - {\bar{\phi }}^j_{Y_{m-1}}(\xi _n)\right) h^{m-1}_{a,k}(\xi _n)\right) \\&=\frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}\left( \sum _{l=1}^{m_0}\left( \phi ^l_{Z_{m-2}}(\xi _n) - {\bar{\phi }}^l_{Z_{m-2}}(\xi _n)\right) {\mathcal {E}}_{l,j}(\xi _n)\right) h^{m-1}_{a,k}(\xi _n)\right) \\&=\frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}\sum _{l=1}^{m_0}\epsilon ({\bar{\phi }}^l_{Z_{m-2}}(\xi _n))\mathcal E_{l,j}(\xi _n)\cdot h^{m-1}_{a,k}(\xi _n)\right) \end{aligned}$$

where $\sum {}^{'}$ indicates that the first and last terms in the sum are halved.

Denoting $\epsilon ({\bar{\phi }}_{Z_{m-2}}):= \max _{1\le n\le N, 1\le k\le m_0}| \epsilon ({\bar{\phi }}^{k}_{Z_{m-2}}(\xi _n))|$, and $\tilde{\epsilon }({\bar{\phi }}^k_{Z_{m-2}}(\xi _n))\cdot \epsilon (\bar{\phi }_{Z_{m-2}})= \epsilon ({\bar{\phi }}^k_{Z_{m-2}}(\xi _n)),$

$$\begin{aligned} \epsilon _2({\bar{\beta }}^j_{Y_{m-1},k})&=\epsilon ({\bar{\phi }}_{Z_{m-2}}) \frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}h^{m-1}_{a,k}(\xi _n)\sum _{l=1}^{m_0}{{\tilde{\epsilon }}}({\bar{\phi }}^l_{Z_{m-2}}(\xi _n))\mathcal E_{l,j}(\xi _n) \right) \\&\le C_\epsilon \cdot \epsilon ({\bar{\phi }}_{Z_{m-2}}) \frac{a^{-1/2}}{\pi }, \end{aligned}$$

for N sufficiently large and some $0<C_\epsilon <\infty$, since we can majorize the error term by an $L^2$ function which admits an upper frame bound. Thus,

$$\begin{aligned} {\mathcal {E}}^j_{m-1,4}(\xi )&:=a^{-1/2}\sum _{k=1}^N\bar{\Psi }_{m-1}(\xi ,k) \cdot \epsilon _2({\bar{\beta }}^j_{Y_{m-1},k}) \\&=a^{-1/2}\epsilon ({\bar{\phi }}_{Z_{m-2}})\sum _{k=1}^N\bar{\Psi }_{m-1}(\xi ,k) \frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}h^{m-1}_{a,k}(\xi _n)\sum _{l=1}^{m_0}{{\tilde{\epsilon }}}({\bar{\phi }}^l_{Z_{m-2}}(\xi _n))\mathcal E_{l,j}(\xi _n) \right) \\&= {\mathcal {O}}\left( \frac{\epsilon ({\bar{\phi }}_{Z_{m-2}})}{a^{1/2}} \left| \sum _{k=1}^N{\bar{\Psi }}_{m-1}(\xi ,k) \frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}h^{m-1}_{a,k}(\xi _n)\sum _{l=1}^{m_0}\mathcal E_{l,j}(\xi _n) \right) \right| \right) \\&= {\mathcal {O}}\left( \frac{\epsilon ({\bar{\phi }}_{Z_{m-2}})}{a^{1/2}} \left| \sum _{k=1}^N{\bar{\Psi }}_{m-1}(\xi ,k) \frac{a^{-1/2}}{\pi } \mathfrak {R}\left( \Delta _\xi \sum _{n=1}^N{}^{'}h^{m-1}_{a,k}(\xi _n)\phi _{Y_1}^j(\xi _n) \right) \right| \right) \\&= {\mathcal {O}}\left( \frac{\epsilon ({\bar{\phi }}_{Z_{m-2}})}{a^{1/2}} \left| \sum _{k=1}^N{\bar{\Psi }}_{m-1}(\xi ,k) \cdot a^{1/2}\Upsilon _{a,N}{\bar{\beta }}_{Y_1,N_{m-1}+k}^j\right| \right) . \end{aligned}$$

Given the exponential decay of ${\bar{\beta }}_{Y_1,n}^j$ for processes satisfying the assumptions above, we have for some constants $C^j$:

$$\begin{aligned} {\mathcal {E}}^j_{m-1,4}(\xi )&= {\mathcal {O}}\left( \frac{\epsilon (\bar{\phi }_{Z_{m-2}})}{a^{1/2}} \left| \sum _{k=1}^N{\bar{\Psi }}_{m-1}(\xi ,k) \cdot a^{1/2}\Upsilon _{a,N}{\bar{\beta }}_{Y_1,k}^j\right| \right) \\&\le C^j \epsilon ({\bar{\phi }}_{Z_{m-2}}) a^{-1/2}|{\bar{\phi }}_{Z_1}^j(\xi )|. \end{aligned}$$

Summarizing,

$$\begin{aligned} \epsilon ({\bar{\phi }}^j_{Z_{m-1}}(\xi ))&\le \tau (G) + |\mathcal E^j_{m-1,1}(\xi )| + |{\mathcal {E}}^j_{m-1,2}(\xi )| + |\mathcal E^j_{m-1,3}(\xi )|+ |{\mathcal {E}}^j_{m-1,4}(\xi )|\\&\le \left( \tau (G) + C_1\cdot \sqrt{{{\bar{a}}}}\cdot \Delta ^4 + C_2 \cdot \sqrt{{{\bar{a}}}} \cdot \epsilon ({\bar{\Psi }}) + \frac{{{\bar{a}}}}{\pi } \epsilon ^M(a, {{\bar{a}}})\right) \\&\quad + a^{-1/2}C^j|{\bar{\phi }}_{Z_1}^j(\xi )| \epsilon ({\bar{\phi }}_{Z_{m-2}})\\&\le \gamma ^M(a,{{\bar{a}}}) + \Omega (a,\xi )\epsilon (\bar{\phi }_{Z_{m-2}}), \end{aligned}$$

where $\gamma ^M(a,{{\bar{a}}})$ is the term in parentheses, and $\Omega (a,\xi ):=a^{-1/2}\max _{1\le j\le m_0}C^j|{\bar{\phi }}_{Z_1}^j(\xi )|$. Using the same logic as above, we can show that $\epsilon ({\bar{\phi }}_{Z_{1}})\le \gamma ^M(a,{{\bar{a}}})$. Iterating from $M-1$ we obtain

$$\begin{aligned} \epsilon ({\bar{\phi }}_{Z_{M-1}}(\xi ))&\le \gamma ^M(a,{{\bar{a}}}) \sum _{m=0}^{M-3} \Omega (a,\xi )^m + \Omega (a,\xi )^{M-2} \epsilon (\bar{\phi }_{Z_{1}})\\&= \gamma ^M(a,{{\bar{a}}}) \frac{1 - \Omega (a,\xi )^{M-2}}{1-\Omega (a,\xi )} + \Omega (a,\xi )^{M-2} \epsilon ({\bar{\phi }}_{Z_{1}})\\&\le \gamma ^M(a,{{\bar{a}}})\frac{1 - \Omega (a,\xi )^{M-1}}{1-\Omega (a,\xi )} \le 2\gamma ^M(a,{{\bar{a}}}), \end{aligned}$$

where the final inequality holds for a sufficiently large. By the definition of $\phi _{Y_{M}}^j(\xi )$, and $\epsilon (\bar{\phi }^j_{Y_{M}}(\xi )):=\phi ^j_{Y_{M}}(\xi ) - \bar{\phi }^j_{Y_{M}}(\xi )$ we have with $\epsilon (\bar{\phi }^k_{Z_{M-1}}(\xi )) = \epsilon ({\bar{\phi }}_{Z_{M-1}}(\xi ))\cdot {{\tilde{\epsilon }}}({\bar{\phi }}^k_{Z_{M-1}}(\xi ))$

$$\begin{aligned} |\epsilon ({\bar{\phi }}^j_{Y_{M}}(\xi ))|&= \left| \sum _{k =1,..., m_0}\epsilon ({\bar{\phi }}^k_{Z_{M-1}}(\xi )){\mathcal {E}}_{k,j}(\xi ) \right| \\&=|\epsilon ({\bar{\phi }}_{Z_{M-1}}(\xi ))| \left| \sum _{k =1,..., m_0}{{\tilde{\epsilon }}}({\bar{\phi }}^k_{Z_{M-1}}(\xi )){\mathcal {E}}_{k,j}(\xi ) \right| \\&\le {{\bar{C}}} \gamma ^M(a,{{\bar{a}}}), \end{aligned}$$

for some ${{\bar{C}}}$ independent of $\xi$. Hence, the characteristic function recursion is stable, and the error in ${\bar{\phi }}_M^j$ can be made arbitrarily small. In particular, with the truncation error $\tau (G)$ (which converges exponentially for processes of interest) is controlled by the choice of ${{\bar{a}}}$, $\epsilon ({\bar{\Psi }})$ sufficiently small by Gaussian quadrature, and $\epsilon ^M(a,\bar{a})$ exponentially convergent by Corollary 3, the error is dominated by that of cubic projection, which is $\mathcal O(\Delta ^4)$. We now state the main result, whose proof is analogous to Section 5.2 in Kirkby (2016), which characterizes the final valuation error.

Proposition 4

Given a European-style payoff $g(A_M)$ on the arithmetic average $A_M$, suppose that the assumptions of Lemma 1 hold, and that ${{\bar{a}}}$ has been fixed sufficiently large to control truncation error. Define the value approximation ${\mathcal {V}}_N\circ g(S_0, \alpha _0)$ as in equation (46), where $N:=a{{\bar{a}}}$. Then the value error, ${\mathcal {E}}({\mathcal {V}}_N) = |{\mathcal {V}}\circ g(S_0, \alpha _0) - {\mathcal {V}}_N\circ g(S_0, \alpha _0)|$ decays as $\mathcal E({\mathcal {V}}_N) = {\mathcal {O}}(a^{-4})$ as the resolution parameter a is increased.

In practice, we typically observe a faster rate of error decay than this result would suggest, but it serves as a conservative estimate.

Additional stochastic volatility models

1.1 Heston and Bate’s models

Consider the Heston stochastic volatility model Heston (1993) with jumps (also known as the Bates model Bates (1996) in the literature):

$$\begin{aligned} \mathbf{Hes}: \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r- q-\lambda \kappa )dt+\sqrt{v_t}dW_t^1 + d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) ,\\ dv_t=\eta (\theta -v_t)dt+\sigma _v\sqrt{v_t}dW^2_t, \end{array} \right. \end{aligned}$$

(53)

where $\eta$ is the mean reversion rate, $\theta$ is the equilibrium level, and $\sigma _v>0$ is the volatility of volatility. Note that to ensure $v_t>0$, the Feller’s condition $2\eta \theta >\sigma _v^2$ is imposed (see Heston (1993)). Applying the transform in equation (21) yields

$$\begin{aligned} {\widetilde{X}}_t=\log \left( \frac{S_t}{S_0}\right) - \rho \int _{v_0}^{v_t}\frac{\varkappa (u)}{{{\hat{\sigma }}}(u)}du =\log \left( \frac{S_t}{S_0}\right) - \frac{\rho }{\sigma _v}(v_t-v_0), \end{aligned}$$

(54)

with $f(v_t,v_0) = \frac{\rho }{\sigma _v}(v_t-v_0)$. It follows that

$$\begin{aligned} \left\{ \begin{array}{l} d{\widetilde{X}}_t=\left[ (\frac{\rho \eta }{\sigma _v}-\frac{1}{2})v_t + (r-q-\frac{\rho \eta \theta }{\sigma _v}-\lambda \kappa )\right] dt+\sqrt{(1-\rho ^2)v_t}dW^*_t +d \left( \sum _{i=1}^{N(t)}J_i\right) ,\\ dv_t=\eta (\theta -v_t)dt+\sigma _v\sqrt{v_t}dW^2_t. \end{array}\right. \end{aligned}$$

(55)

1.2 3/2 Model with jumps

Next consider the dynamics of the 3/2 stochastic volatility model with jumps

$$\begin{aligned} \mathbf 3/2: \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r-q-\lambda \kappa )dt+\sqrt{v_t}dW_t^1 +d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) , \\ dv_t=v_t[\eta (\theta -v_t)dt+\sigma _v \sqrt{v_t}dW^2_t],\ v(0)=v_0, \end{array} \right. \end{aligned}$$

(56)

where in (56) $v_t$ is the variance of the asset $S_t$, r is the risk-free interest rate, $\sigma _v>0$, $\theta \in {\mathbb {R}}$ is the mean reversion level, $\eta$ is given such that $\eta v_t$-a stochastic volatility quantity-is the speed of mean reversion.

While this form can be applied directly, we prefer an alternative formulation which nests the 3/2 model within the 4/2 model introduced in Sect. 3.1.1. Applying Ito’s formula we have

$$\begin{aligned}&d\left( \frac{1}{v_t}\right) =\eta \theta \left( \frac{\eta +\sigma _v^2}{\eta \theta }-\frac{1}{v_t} \right) dt-\displaystyle \frac{\sigma _v}{\sqrt{v_t}}dW_t^2 \end{aligned}$$

(57)

$$\begin{aligned}&\widehat{v_t}:=\frac{1}{v_t},\quad {\widehat{\eta }}:=\eta \theta ,\quad {\widehat{\theta }}:=\frac{\eta +\sigma _v^2}{\eta \theta },\quad {\widehat{\sigma }}_v:=-\sigma _v, \end{aligned}$$

(58)

from which (56) is reduced to

$$\begin{aligned} \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r- q -\lambda \kappa )dt+\displaystyle \frac{1}{\sqrt{{\widehat{v}}_t}}dW_t^1 +d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) , \\ d{\widehat{v}}_t={\widehat{\eta }}[{\widehat{\theta }}-\widehat{v_t}]dt+{\widehat{\sigma }}_v \sqrt{{\widehat{v}}_t}dW_t^2,\quad {\widehat{v}}_0=1/v_0. \end{array} \right. \end{aligned}$$

(59)

Equation equation (21) prescribes the change of variables

$$\begin{aligned} \widetilde{X}_t=\log \Big (\frac{S_t}{S_0}\Big )-\frac{\rho }{{\widehat{\sigma }}_v}\log \Big (\frac{{\widehat{v}}_t}{{\widehat{v}}_0}\Big ), \end{aligned}$$

(60)

which results in the decorrelated dynamics

$$\begin{aligned} \left\{ \begin{array}{l} d{\widetilde{X}}_t=\displaystyle \Big [\frac{1}{2}\Big (\rho {\widehat{\sigma }}_v-2\frac{\rho {\widehat{\eta }}{\widehat{\theta }}}{{\widehat{\sigma }}_v}-1\Big )\frac{1}{{\widehat{v}}_t}+\frac{\rho {\widehat{\eta }}}{{\widehat{\sigma }}_v}+(r-q -\lambda \kappa )\Big ]dt+\sqrt{\frac{(1-\rho ^2)}{{\widehat{v}}_t}}dW^*_t +d \left( \sum _{i=1}^{N(t)}J_i\right) ,\\ d{\widehat{v}}_t={\widehat{\eta }}[{\widehat{\theta }}-\widehat{v_t}]dt+{\widehat{\sigma }}_v\sqrt{{\widehat{v}}_t} dW_t^2,\quad {\widehat{v}}_0=1/v_0. \end{array}\right. \end{aligned}$$

(61)

We see that the 4/2 model in (27) indeed contains the 3/2 model as a special case where $a= 0$, $b=1$, and the parameters for 3/2 are selected using the re-parameterization in equation (58).

1.3 Hull-White’s model with jumps

Augmenting the traditional Hull and White (1990) model with jumps produces

$$\begin{aligned} \mathbf{HW}: \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r-q -\lambda \kappa )dt+\sqrt{v_t}dW_t^1 +d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) ,\\ dv_t=a_vv_tdt+\sigma _vv_tdW^2_t. \end{array} \right. \end{aligned}$$

(62)

The change of variable that will help us to remove the correlation, $\rho$, between the two stochastic processes $W_t^1$, $W_t^2$ in (62) is given by

$$\begin{aligned} \widetilde{X}_t=\log \Big (\frac{S_t}{S_0}\Big )-\frac{2\rho }{\sigma _v}(\sqrt{v_t}-\sqrt{v_0}). \end{aligned}$$

(63)

The decorrelated dynamics satsify

$$\begin{aligned} \left\{ \begin{array}{l} d\widetilde{X}_t=\Big [\Big (\displaystyle \frac{\rho \sigma _v}{4}-\frac{a_v\rho }{\sigma _v}\Big )\sqrt{v_t}-\frac{1}{2}v_t +(r-q -\lambda \kappa )\Big ]dt+\sqrt{(1-\rho ^2)v_t}dW^*_t +d \left( \sum _{i=1}^{N(t)}J_i\right) ,\\ dv_t=a_v v_tdt+\sigma _v v_tdW^2_t. \end{array}\right. \end{aligned}$$

(64)

1.4 Stein-Stein’s model with jumps

Stein and Stein (1991) consider a stochastic volatility model where the two Brownian motions are independent. In this section, we extend their model by allowing for correlation and add a jump component. More specifically, the dynamics of the model is specified as follow:

$$\begin{aligned} \mathbf{SS}: \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r-q -\lambda \kappa )dt+v_tdW_t^1 +d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) ,\\ dv_t=\eta (\theta -v_t)dt+\sigma _vdW^2_t, \end{array} \right. \end{aligned}$$

(65)

The change of variable that will help us to remove the correlation between the two stochastic processes $W_t^1$, $W_t^2$ in (65) is given by

$$\begin{aligned} \widetilde{X}_t=\log \Big (\frac{S_t}{S_0}\Big )-\frac{1}{2}\frac{\rho }{\sigma _v}(v_t^2-v_0^2). \end{aligned}$$

(66)

Then, we have

$$\begin{aligned} \left\{ \begin{array}{l} d{\widetilde{X}}_t=\Big [\Big (\displaystyle \frac{\rho \eta }{\sigma _v}-\frac{1}{2}\Big )v_t^2-\frac{\rho \eta \theta }{\sigma _v}v_t+r-q -\lambda \kappa -\frac{\rho \sigma _v}{2}\Big ]dt+v_t\sqrt{(1-\rho ^2)}dW^*_t +d \left( \sum _{i=1}^{N(t)}J_i\right) ,\\ dv_t=\eta (\theta -v_t)dt+\sigma _vdW^2_t. \end{array}\right. \end{aligned}$$

(67)

1.5 $\alpha$-Hypergeometric model

The $\alpha$-Hypergeometric model was recently proposed by Da Fonseca and Martini (2016). Unlike Heston model, for $\alpha$-Hypergeometric model the strict positivity of volatility is guaranteed. The dynamics of the stock price is given by

$$\begin{aligned} \left\{ \begin{array}{l} \frac{dS_t}{S_{t^-}}=(r-q -\lambda \kappa )dt+ e^{v_t}dW_t^1+d\left( \sum _{i=1}^{N(t)} (e^{J_i}-1)\right) , \\ dv_t=(\eta -\theta e^{a_v v_t})dt+\sigma _v dW^2_t,\ v(0)=v_0, \end{array} \right. \end{aligned}$$

(68)

where $\eta ,v_0\in (-\infty ,+\infty ), \theta>0,\sigma _v>0,a_v>0$. Let

$$\begin{aligned} X_t=\log (\frac{S_t}{S_0})-\frac{\rho }{\sigma _v}(e^{v_t}-e^{v_0})-(r-q -\lambda \kappa )t, \end{aligned}$$

(69)

then we have

$$\begin{aligned} \left\{ \begin{array}{l} dX_t=\Big [\frac{\rho \theta }{\sigma _v}e^{(1+a_v)v_t}-\rho (\frac{\eta }{\sigma _v}+\frac{\sigma _v}{2})e^{v_t}-\frac{e^{2v_t}}{2}\Big ]dt +e^{v_t}\sqrt{1-\rho ^2}dW^*_t+d \left( \sum _{i=1}^{N(t)}J_i\right) , \\ dv_t=(\eta -\theta e^{a_v v_t})dt+\sigma _v dW^2_t. \end{array} \right. \end{aligned}$$

(70)

1.6 Jacobi model

An interesting recent SV model in the literature is the Jacobi model (without jump) of Ackerer et al. (2016), which specifies a bounded variance process Q(v), where $v\in [v_{min},v_{max}]$ for $0\le v_{min}<v_{max}$, defined by the quadratic function

$$\begin{aligned} Q(v)= \frac{(v-v_{min})(v_{max}-v)}{(\sqrt{v_{max}} - \sqrt{v_{min}})^2}, \quad 0\le Q(v) \le v, \quad v\in [v_{min},v_{max}]. \end{aligned}$$

With $Z_t:=\log (S_t)$, the dynamics under the Jacobi model with jumps are

$$\begin{aligned} \left\{ \begin{array}{ll} dZ_t=(r-q -\lambda \kappa - v_t/2)dt + \sqrt{v_t -\rho ^2Q(v_t)}dW_t^{*} + \rho \sqrt{Q(v_t)}dW_t^{(2)} +d \left( \sum _{i=1}^{N(t)}J_i\right) ,\\ dv_t =\eta (\theta -v_t)dt+\alpha \sqrt{Q(v_t)} dW_t^{(2)}. \end{array} \right. \end{aligned}$$

(71)

where $\eta \ge 0$, $\theta \in [v_{min},v_{max}]$, and $\alpha >0$. Note here that ${\mathbb {E}}[dW_t^{(2)}dW_t^{*}] = 0$, as the correlation structure is already incorporated with correlation parameter $\rho$. From the equality

$$\begin{aligned} \rho \int _0^t \sqrt{Q(v_s)}dW_s^{(2)} = \frac{\rho }{\alpha }(v_t - v_0) -\frac{ \rho }{\alpha }\int _{0}^t\eta (\theta -v_s)ds, \end{aligned}$$

we can derive the auxiliary process ${\tilde{X}}_t := Z_t - \frac{ \rho }{\alpha }v_{t}$ with

$$\begin{aligned} d{\tilde{X}}_t&= \left( r-q -\lambda \kappa - \frac{v_{t}}{2} -\frac{\rho }{\alpha }\eta (\theta -v_{t})\right) dt + \sqrt{v_{t} -\rho ^2Q(v_{t})}dW_t^{*}+d \left( \sum _{i=1}^{N(t)}J_i\right) \\&= \left( \left( r-q -\lambda \kappa - \frac{ \rho }{\alpha }\eta \theta \right) + v_{t}\left( \frac{ \rho }{\alpha }\eta - \frac{1}{2} \right) \right) dt + \sqrt{v_{t} -\rho ^2Q(v_{t})}dW_t^{*}\\&\quad +d \left( \sum _{i=1}^{N(t)}J_i\right) . \end{aligned}$$