3.1. The HE and DHE distributions

Given a set of $n$ exponential random variables $\varepsilon _i\sim \mathsf {Exp}(\eta _i)$, $i=1,\ldots,n$, each with a density function $f_{\varepsilon _i}(y) = \eta _i e^{-\eta _i y}$, $y\ge 0$. Suppose that $\eta _i\ne \eta _j$ for all $i, j = 1,\ldots,n, i\ne j$. The distribution followed by a sum of these nonhomogeneous exponential random variables is called a HE distribution (see, e.g. [Reference Ross12]), which is denoted as below:

$$H_{(n)} = \sum_{i=1}^n \varepsilon_i\sim\mathsf{HE}(\boldsymbol{\eta}_{(n)}),\quad \text{where}\ \boldsymbol{\eta}_{(n)} = [\eta_1, \eta_2,\ldots, \eta_n] \in \mathbb{R}^n.$$

For the above nonnegative HE random variable $H_{(n)}$, its probability density function (pdf), moment generating function (mgf), and moments are given analytically as follows.

Lemma 1. For a HE random variable $H_{(n)}\sim \mathsf {HE}(\boldsymbol {\eta }_{(n)})$, its pdf is given by

(3.4) \begin{equation} f_{H_{(n)}}(y) = \sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)}) \eta_i e^{-\eta_i y} {1}_{\{y\geq 0\}}, \end{equation}

where $\Omega (i;\boldsymbol {\eta }_{(n)}) = \prod _{j=1, j\neq i}^n {\eta _j}/{(\eta _j-\eta _i)},$ and its mgf is given by

(3.5) \begin{equation} \mathrm{E}[e^{\alpha H_{(n)}}] = \prod_{i=1}^n \frac{\eta_i}{\eta_i-\alpha}. \end{equation}

The $k$th moment of $H_{(n)}$, $k=1,2,\ldots$, can be obtained from the following formula

(3.6) \begin{equation} \mathrm{E}[H_{(n)}^k] = ({-}1)^k k! \,\boldsymbol{\pi} \mathbf{D}^{{-}k} \mathbf{1}, \end{equation}

where

$$\boldsymbol{\pi}=\underbrace{[1,0,\ldots,0]}_{1\times n}, \quad \mathbf{D} = \underbrace{\left[\begin{array}{ccccc} -\eta_1 & \eta_1 & \cdots & 0 & 0 \\ 0 & -\eta_2 & \cdots & 0 & 0 \\ \cdots & \cdots & \ddots & \cdots & \cdots \\ 0 & 0 & \cdots & -\eta_{n-1} & \eta_{n-1} \\ 0 & 0 & \cdots & 0 & -\eta_n \end{array}\right]}_{n\times n}, \quad \mathbf{1} = \underbrace{\left[\begin{array}{c} 1 \\ \vdots \\ 1 \end{array}\right]}_{n\times 1}.$$

It is clear from (3.3) that for the analysis of $X_t$, we need to work with the difference of two nonnegative random variables. For the convenience of our presentation, below we give a fundamental formula that will be used repeatedly in our subsequent derivations.

Lemma 2. Let $X_1$ and $X_2$ be two continuous nonnegative random variables with density functions $f_1(x)$ and $f_2(x)$. The density function of $X_1-X_2$ is given by

(3.7) \begin{equation} g(y) = \left\{\begin{array}{ll} \displaystyle\int_0^\infty f_1(x+y) f_2(x)\, dx, & y \ge 0, \\[9pt] \displaystyle\int_0^\infty f_1(x) f_2(x-y)\, dx, & y \lt 0. \end{array}\right. \end{equation}

Proof. Suppose that $X_1$ and $X_2$ are continuous real-valued random variables, the density function of $X_1-X_2$ takes the following (convolution) form:

$$g(y) = \int_{-\infty}^\infty f_1(x) f_2(x-y)\, dx.$$

It simplifies to (3.7) by noting that the arguments of both density functions take nonnegative values as $X_1$ and $X_2$ are nonnegative in our case.

Now, we introduce a new distribution which is followed by the difference of two HE random variables. Let $\widetilde {H}_{(m)}\sim \mathsf {HE}(\boldsymbol {\tilde {\eta }}_{(m)})$ where $\boldsymbol {\tilde {\eta }}_{(m)} = [\tilde {\eta }_1, \tilde {\eta }_2,\ldots, \tilde {\eta }_m]$ $\in \mathbb {R}^m$ be another HE random variable. The distribution of the difference of $H_{(n)}$ and $\widetilde {H}_{(m)}$ is called a differenced hypoexponential (DHE) distribution, which is denoted as below:

$$H_{(n)} - \widetilde{H}_{(m)} \sim\mathsf{DHE}(\boldsymbol{\eta}_{(n)}, \widetilde{\boldsymbol{\eta}}_{(m)}).$$

The following result shows that for this DHE distribution, the pdf, mgf, and moments can be derived in closed-form.

Lemma 3. The pdf of the DHE random variable $H_{(n)} - \widetilde {H}_{(m)}$ is given by

(3.8) \begin{equation} f_{H_{(n)} - \widetilde{H}_{(m)}}(y) = \sum_{i=1}^n \sum_{j=1}^{m} \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)}) \frac{\eta_i\tilde{\eta}_j}{\eta_i + \tilde{\eta}_j} (e^{-\eta_i y}{1}_{\{y\geq 0\}} + e^{\tilde{\eta}_j y}{1}_{\{y \lt 0\}}). \end{equation}

Its mgf is given by

(3.9) \begin{equation} \mathrm{E}[e^{\alpha (H_{(n)}-\widetilde{H}_{(m)})}] = \left(\prod_{i=1}^n \frac{\eta_i}{\eta_i-\alpha}\right)\left(\prod_{j=1}^m \frac{\tilde{\eta}_j-\alpha}{\tilde{\eta}_j}\right). \end{equation}

The $k$th moment of the DHE random variable, $k=1,2,\ldots$, can be expressed in terms of the moments of the constituent HE random variables as follows:

(3.10) \begin{equation} \mathrm{E}[(H_{(n)}-\widetilde{H}_{(m)})^k] =\sum_{j=0}^k \binom{k}{j}({-}1)^j\, \mathrm{E}[H_{(n)}^{k-j}] \,\mathrm{E}[\widetilde{H}_{(m)}^j]. \end{equation}

Proof. The pdf in (3.8) is obtained by (3.4) and (3.7) as follows. For $y\ge 0$,

\begin{align*} f_{H_{(n)} - \widetilde{H}_{(m)}}(y) & = \int_0^\infty f_{H_{(n)}}(x+y;\boldsymbol{\eta}_{(n)}) ~ f_{\widetilde{H}_{(m)}}(x;\boldsymbol{\tilde{\eta}}_{(m)}) \,dx \nonumber \\ & =\sum_{i=1}^n \sum_{j=1}^m \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)}) \frac{\eta_i\tilde{\eta}_j}{\eta_i + \tilde{\eta}_j} e^{-\eta_i y}, \end{align*}

whereas for $y \lt 0$,

\begin{align*} f_{H_{(n)} - \widetilde{H}_{(m)}}(y) & = \int_0^\infty f_{H_{(n)}}(x;\boldsymbol{\eta}_{(n)}) ~ f_{\widetilde{H}_{(m)}}(x-y;\boldsymbol{\tilde{\eta}}_{(m)}) \,dx \nonumber \\ & =\sum_{i=1}^n \sum_{j=1}^{m} \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)}) \frac{\eta_i\tilde{\eta}_j}{\eta_i + \tilde{\eta}_j} e^{\tilde{\eta}_j y}. \end{align*}

The mgf in (3.9) is obtained from the mgf of the two constituent HE random variables. The moment formula in (3.10) is obtained from the binomial theorem.

3.2. The distribution of the jump accumulation

Now that the HE and DHE distributions are well defined, the (conditional) jump accumulation process in (3.3) can be written as $X_t^{(n,m)} = H_{(m)} - \widetilde {H}_{(n-m)}$ which follows a DHE distribution with density function given in Lemma 3. The unconditional version of $X_t^{(n,m)}$ can be written as

(3.11) \begin{equation} X_t = \sum_{j=1}^{N_t} \ln Y_j = H_{(M_t)} - \widetilde{H}_{(N_t-M_t)} = \sum_{k=1}^{M_t} \varepsilon_k - \sum_{\ell=1}^{N_t-M_t} \tilde{\varepsilon}_\ell. \end{equation}

Its distribution can be obtained by unconditioning on $N_t=n$ and $M_t=m$ as follows.

Proposition 4. The pdf of $X_t$ is given by

\begin{align*} & f_{X_t}(y) = e^{-\lambda t}\delta(y) \nonumber\\ & \quad + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda t) \left[p^n \sum_{k=1}^n\Omega(k;\boldsymbol{\eta}_{(n)}) \eta_k e^{-\eta_k y} {1}_{\{y\geq 0\}} + (1-p)^n \sum_{\ell=1}^n\Omega(\ell;\widetilde{\boldsymbol{\eta}}_{(n)}) \eta_\ell e^{\eta_\ell y} {1}_{\{y \lt 0\}}\right] \nonumber \\ & \quad + \sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda t) \left[\binom{n}{m}p^m(1-p)^{n-m}\vphantom{\sum_{\ell=1}^{n-m}}\right.\nonumber\\ & \quad \times\left.\sum_{k=1}^m \sum_{\ell=1}^{n-m} \Omega (k;\boldsymbol{\eta}_{(m)})\Omega (\ell;\boldsymbol{\tilde{\eta}}_{(n-m)}) \, \frac{\eta_k\tilde{\eta}_\ell}{\eta_k + \tilde{\eta}_\ell} (e^{-\eta_k y}{1}_{\{y\geq 0\}} + e^{\tilde{\eta}_\ell y}{1}_{\{y \lt 0\}})\right], \end{align*}

where $\delta (y)$ is the Dirac delta function. The mgf of $X_t$ is given by

(3.12) \begin{align} \mathrm{E}[e^{\alpha X_t}] & = \mathrm{E}[e^{\alpha (\sum_{k=1}^{M_t} \varepsilon_k - \sum_{\ell=1}^{N_t-M_t} \tilde{\varepsilon}_\ell)}] \nonumber\\ & = e^{-\lambda t} + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda t) \left[p^n \prod_{k=1}^n \frac{\eta_k}{\eta_k-\alpha} + (1-p)^n \prod_{\ell=1}^n \frac{\tilde{\eta}_\ell-\alpha}{\tilde{\eta}_\ell} \right] \nonumber\\ & \quad + \sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda t) \left[\binom{n}{m}p^m(1-p)^{n-m} \left(\prod_{k=1}^m \frac{\eta_k}{\eta_k-\alpha}\right)\left(\prod_{\ell=1}^{n-m} \frac{\tilde{\eta}_\ell-\alpha}{\tilde{\eta}_\ell}\right)\right]. \end{align}

The $k$th moments of $X_t$, $k=1,2,\ldots$, are given by

(3.13) \begin{equation} \mathrm{E}[X_t^k] = \sum_{n=0}^{n^*}\sum_{m=0}^{n} f_{\mathsf{SP}}(n;\lambda t) \binom{n}{m}p^m(1-p)^{n-m} \left[\sum_{j=0}^k \binom{k}{j}({-}1)^j\, \mathrm{E}[H_{(m)}^{k-j}] \,\mathrm{E}[\widetilde{H}_{(n-m)}^j]\right], \end{equation}

where we define $H_{(0)} = \widetilde {H}_{(0)} = 0$ and $H_{(0)}^0 = \widetilde {H}_{(0)}^0 = 1$.

Proof. For a given $N_t = n\ge 1$, the probability of $M_t = m$ is $\binom {n}{m}p^m(1-p)^{n-m}$. Therefore, $X_t^{(n)}$ can be expressed as a mixture of HE and DHE random variables as described below:

$$X_t^{(n)} = \left\{\begin{array}{ll} X^{(n,n)} = H_{(n)}, & \text{with prob}\ p^n\ \text{(when}\ m=n), \\ X^{(n,0)} ={-}\widetilde{H}_{(n)}, & \text{with prob}\ (1-p)^n\ \text{(when}\ m=0),\\ X^{(n,m)} = H_{(m)} -\widetilde{H}_{(n-m)}, & \text{with prob}\ \binom{n}{m}p^m(1-p)^{n-m},\ m=1,\ldots,n-1. \end{array}\right.$$

Moreover, its unconditional version $X_t$ is a mixture of 0 and $X_t^{(n)}$ where $n\ge 1$:

$$X_t = \left\{\begin{array}{ll} 0, & \text{for} \ N_t=0\ \text{with probability}\ e^{-\lambda t}, \\ X_t^{(n)}, & \text{for} \ N_t = n \ge 1\ \text{with probability}\ f_{\mathsf{SP}}(n;\lambda t). \end{array}\right.$$

The claimed results are obtained by unconditioning on $N_t$ and $M_t$. Specifically, the pdf is

\begin{align*} f_{X_t}(y) & = e^{-\lambda t}\delta(y) + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda t) [p^n f_{H_{(n)}}(y;\boldsymbol{\eta}_{(n)}) + (1-p)^n f_{\widetilde{H}_{(n)}}({-}y;\widetilde{\boldsymbol{\eta}}_{(n)})]\\ & \quad +\sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda t) \left[\binom{n}{m}p^m(1-p)^{n-m} f_{H_{(m)} - \widetilde{H}_{(n-m)}}(y; \boldsymbol{\eta}_{(m)},\boldsymbol{\tilde{\eta}}_{(n-m)})\right], \end{align*}

which leads to (3.12). The mgf in (3.12) is obtained in a similar way. The $k$th moment formula (3.13) is obtained from its conditional version by using the binomial theorem:

$$\mathrm{E}[(X_t^{(n,m)})^k]= \mathrm{E}[(H_{(m)}-\widetilde{H}_{(n-m)})^k] =\sum_{j=0}^k \binom{k}{j}({-}1)^j \,\mathrm{E}[H_{(m)}^{k-j}]\, \mathrm{E}[\widetilde{H}_{(n-m)}^j].$$

The density function of $X_t$ is in fact a mixture of a discrete random variable (i.e., a constant 0 corresponding to $N_t = 0$; it is described by the Dirac delta function in (3.12)) and many continuous (HE, DHE) random variables. A typical density function is depicted in Figure 5 where $x=0$ is a point carrying nonzero probability measure.

Figure 5. The density function of $X_t$ where the delta function at $x=0$ corresponds to $N_t = 0$.

4.1. The risk-neutral stock return distribution

Let $T$ denote the maturity time of the option, and the time variable $t$ as seen in previous sections is now replaced by the constant $T$ for convenience. Let $R_T$ denote the stock return over the time period $[0,T]$. By (3.1), we have

(4.1) \begin{equation} R_T = \ln\left(\frac{S_T}{S_0}\right) = \left(\mu - \frac{\sigma^2}{2}\right) T + \sigma W_T + X_T. \end{equation}

It is clear that $R_T$ can be seen as a normal random variable (diffusion) added to the compound Poisson random variable (jump accumulation) $X_T$. The key to the analysis is the treatment of the sum and difference between a normal and the HE/DHE distributions.

When it comes to option pricing, we need to consider the risk-neutral measure under which the martingale condition must hold, that is, $\mathrm {E}[S_T] = S_0 e^{(r-q)T}$ (where $r$ and $q$ stand for the risk-free rate and stock dividend yield). The free variable $\mu$ in (4.1) can be adjusted to make this condition hold, that is, $\mu = r - q - \chi$ where

(4.2) \begin{align} \chi & =\frac{1}{T}\ln(\mathbb{E}[e^{X_T}]) \nonumber\\ & =\frac{1}{T}\ln\left[e^{-\lambda T} + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda T) \left(p^n \prod_{k=1}^n \frac{\eta_k}{\eta_k-1} + (1-p)^n \prod_{\ell=1}^n \frac{\tilde{\eta}_\ell-1}{\tilde{\eta}_\ell} \right) \right. \nonumber\\ & \quad \left. +\sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda T) \binom{n}{m}p^m(1-p)^{n-m} \left(\prod_{k=1}^m \frac{\eta_k}{\eta_k-1}\right)\left(\prod_{\ell=1}^{n-m} \frac{\tilde{\eta}_\ell-1}{\tilde{\eta}_\ell}\right)\right]. \end{align}

Note that $\chi$ defined above may have dependence on $T$. In (4.1), let $\nu _T = (\mu -{\sigma ^2}/{2})T = (r-q-\chi -{\sigma ^2}/{2})T$, then we have the following simpler expression for stock return under the risk-neutral measure

(4.3) \begin{equation} R_T = \nu_T + \sigma W_T + X_T. \end{equation}

For the convenience of the subsequent analysis, let us introduce the following three new distributions that are extended from HE and DHE distributions (where $Z \sim \mathsf {N}(0,\sigma ^2)$):

• • NHE$^+$ distribution, defined as the sum of normal and HE distributions, that is,

  $$Z + H_{(n)} \sim\mathsf{{NHE}}^+(\sigma^2,\boldsymbol{\eta}_{(n)}).$$

• • NHE$^-$ distribution, defined as the difference of normal and HE distributions, that is,

  $$Z - \widetilde{H}_{(n)} \sim\mathsf{{NHE}}^-(\sigma^2,\widetilde{\boldsymbol{\eta}}_{(n)}).$$

• • NDHE distribution, defined as the sum of normal and DHE distributions, that is,

  $$Z + H_{(n)} - \widetilde{H}_{(m)}\sim\mathsf{NDHE}(\sigma^2,\boldsymbol{\eta}_{(n)}, \widetilde{\boldsymbol{\eta}}_{(m)}).$$

The following lemma shows that the density functions for the above three distributions can be derived in closed-form.

Lemma 5. The density functions of the NHE$^+$, NHE$^-$, and NDHE distributions as defined above are given respectively by

(4.4) \begin{align} f_{Z+H_{(n)}}(y) & = \phi^{+}(y;\sigma^2,\boldsymbol{\eta}_{(n)}) = \sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)})\eta_i e^{{(\sigma\eta_i)^2}/{2}-y\eta_i}\Phi\left(\frac{y}{\sigma}-\sigma\eta_i\right), \end{align}

(4.5) \begin{align} f_{Z-\widetilde{H}_{(n)}}(y) & =\phi^{-}(y;\sigma^2,\widetilde{\boldsymbol{\eta}}_{(n)}) = \sum_{i=1}^n\Omega(i;\widetilde{\boldsymbol{\eta}}_{(n)}) \tilde{\eta}_i e^{{(\sigma\tilde{\eta}_i)^2}/{2}+y\tilde{\eta}_i} \Phi\left(\frac{-y}{\sigma} - \sigma\tilde{\eta}_i \right), \end{align}

(4.6) \begin{align} f_{Z + H_{(n)} -\widetilde{H}_{(m)}}(y) & = \varphi(y;\sigma^2,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)}) = \sum_{i=1}^n \sum_{j=1}^m \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)})\left[\frac{\eta_i\widetilde{\eta_j}}{\eta_i + \tilde{\eta}_j}\right. \nonumber\\ & \quad \left.\times\left(e^{{(\sigma\eta_i)^2}/{2}-y\eta_i} \Phi\left(\frac{y}{\sigma}-\sigma\eta_i\right) +e^{{(\sigma\tilde{\eta}_j)^2}/{2}+y\tilde{\eta}_j}\Phi\left(\frac{-y}{\sigma}-\sigma\tilde{\eta}_j\right) \right)\right]. \end{align}

Proof. Firstly, by using (3.4) together with (3.7) (Lemma 2), we may derive

\begin{align*} \phi^{+}(y;\sigma^2,\boldsymbol{\eta}_{(n)}) & = \int_{-\infty}^\infty f_{Z}(x) f_{H_{(n)}}(y-x) \,dx \\ & =\int_{-\infty}^y \tfrac{1}{\sigma\sqrt{2\pi}} e^{-{x^2}/{2\sigma^2}} \left(\sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)}) \eta_i e^{-\eta_i(y-x)}\right) dx \nonumber\\ & = \sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)}) \eta_i e^{{(\sigma\eta_i)^2}/{2}-y\eta_i} \int_{-\infty}^y \tfrac{1}{\sigma\sqrt{2\pi}} e^{-{(x-\sigma^2\eta_i)^2}/{2\sigma^2}}\, dx, \end{align*}

which leads to (4.4). The formula (4.5) for $\phi ^-(y;\sigma ^2,\boldsymbol {\eta }_{(n)})$ can be derived in a similar way. Lastly, by using (3.8), we derive

\begin{align*} \varphi(y;\sigma^2,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)}) & = \int_{-\infty}^\infty f_{Z}(x) f_{H_{(n)} - \widetilde{H}_{(m)}}(y-x) \,dx \nonumber \\ & =\sum_{i=1}^n \sum_{j=1}^m \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)}) \nonumber \\ & \quad \times \int_{-\infty}^\infty f_{Z}(x) \left(\frac{\tilde{\eta}_j}{\eta_i + \tilde{\eta}_j} ~f_{\varepsilon_i}(x-y) + \frac{\eta_i}{\eta_i + \tilde{\eta}_j} ~f_{\tilde{\varepsilon}_j}(y-x)\right) dx, \end{align*}

which leads to (4.6).

Using the above results, the density function of $R_T$ can be expressed in terms of the $\phi ^+$, $\phi ^-$, and $\varphi$ functions as given in (4.4), (4.5), and (4.6).

Proposition 6. The density function of $R_T$ is given by

(4.7) \begin{align} f_{R_T}(y) & = \frac{1}{\sqrt{2\pi\sigma^2 T}}e^{-({(y-\nu_T)^2}/{2\sigma^2 T}+\lambda T)} \nonumber\\ & \quad +\sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda T) [p^n \phi^+(y-\nu_T;\sigma^2 T, \boldsymbol{\eta}_{(n)}) + (1-p)^n \phi^-(y-\nu_T;\sigma^2 T, \widetilde{\boldsymbol{\eta}}_{(n)})] \nonumber\\ & \quad +\sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda T)\left[\binom{n}{m}p^m(1-p)^{n-m}\varphi(y-\nu_T;\sigma^2 T, \boldsymbol{\eta}_{(m)}, \widetilde{\boldsymbol{\eta}}_{(n-m)})\right], \end{align}

where $\nu _T = (r-q-\chi -{\sigma ^2}/{2})T$ with $\chi$ defined in (4.2).

Proof. Conditional on $N_T = n$ ($n\ge 1$), $M_T = m$ ($m = 0,1,\ldots,n$), we have

$$R_T^{(n,m)} = \nu_T + \sigma W_T + X_T^{(n,m)}.$$

Since $\sigma W_T\sim \mathsf {N}(0,\sigma ^2 T)$, it is clear that

• • When $m=n$ (where $n\ge 1$):

  $$\sigma W_T+ X_T^{(n,m)} = \sigma W_T + \sum_{i=1}^n \varepsilon_i\sim\mathsf{NHE}^+(\sigma^2 T, \boldsymbol{\eta}_{(n)}).$$

• • When $m=0$ (where $n\ge 1$):

  $$\sigma W_T + X_T^{(n,m)} = \sigma W_T - \sum_{i=1}^n \tilde{\varepsilon}_i\sim\mathsf{NHE}^-(\sigma^2 T, \widetilde{\boldsymbol{\eta}}_{(n)}).$$

• • When $1\le m\le n-1$ (where $n\ge 2$):

  $$\sigma W_T + X_T^{(n,m)} = \sigma W_T + \sum_{i=1}^m \varepsilon_i - \sum_{j=1}^{n-m} \tilde{\varepsilon}_j\sim\mathsf{NDHE}(\sigma^2 T, \boldsymbol{\eta}_{(m)}, \widetilde{\boldsymbol{\eta}}_{(n-m)}).$$

The density function of $R_T$ is then obtained by unconditioning on $N_t$ and $M_t$.

The mgf of $R_T$ can be easily obtained as a product of a normal mgf and (3.12). As to the $k$th moment of $R_T$, its general form can be obtained as

(4.8) \begin{equation} \mathrm{E}[R_T^k] = \sum_{j=0}^{k} \binom{k}{j}\, \mathrm{E}[(\nu_T+\sigma W_T)^j]\, \mathrm{E}[X_T^{k-j}]. \end{equation}

We are mainly interested in the first four moments of $R_T$, which can be obtained together with the moments of normal random variable $\nu _T + \sigma W_T\sim \mathsf {N}(\nu _T,\sigma ^2 T)$ as given below:

\begin{align*} & \mathrm{E}[\nu_T+\sigma W_T] = \nu_T, \quad \mathrm{E}[(\nu_T+\sigma W_T)^2] = \nu_T^2 + \sigma^2 T, \\ & \mathrm{E}[(\nu_T+\sigma W_T)^3] = \nu_T^3 + 3\nu_T\sigma^2 T, \quad \mathrm{E}[(\nu_T+\sigma W_T)^4] = \nu_T^4 + 6\nu_T^2\sigma^2 T + 3\sigma^4 T^2. \end{align*}

4.2. Option pricing formula

We are now in a position to derive the pricing formula for a European option. Consider a call option with strike price $K$ and maturity time $T$. Its price is given by the following risk-neutral pricing formula:

(4.9) \begin{equation} C = e^{{-}rT}\,\mathrm{E}[(S_T-K)^+] = e^{{-}rT}\,\mathrm{E}[S_T {1}_{\{S_T\ge K\}}] - Ke^{{-}rT}\mathrm{P}(S_T\ge K). \end{equation}

To deal with the probability $\mathrm {P}(S_T\ge K)$ in the second term, we define the complementary distribution functions of the NHE$^+$, NHE$^-$, and NDHE distributions (where $Z \sim \mathsf {N}(0,\sigma ^2)$).

\begin{align*} \bullet \quad & \mathrm{U}^+(x;\sigma^2,\boldsymbol{\eta}_{(n)}) = \mathrm{P}(Z + H_{(n)} \geq x),\\ \bullet \quad & \mathrm{U}^-(x;\sigma^2,\widetilde{\boldsymbol{\eta}}_{(m)}) =\mathrm{P}(Z - \widetilde{H}_{(m)} \geq x),\\ \bullet \quad & \mathrm{V}(x;\sigma^2,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)}) =\mathrm{P}(Z + H_{(n)} - \widetilde{H}_{(m)} \geq x). \end{align*}

As it turns out, the closed-form expressions of these functions can be derived.

Lemma 7. The closed-form expressions of the three functions $\mathrm {U}^+, \mathrm {U}^-, \mathrm {V}$ as defined above are respectively given by

(4.10) \begin{align} \mathrm{U}^+(x;\sigma^2,\boldsymbol{\eta}_{(n)}) & = \sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)})\left(\Phi\left(-\frac{x}{\sigma}\right) + e^{{(\sigma\eta_i)^2}/{2}-x\eta_i} \Phi\left(\frac{x}{\sigma} - \sigma\eta_i \right)\right), \end{align}

(4.11) \begin{align} \mathrm{U}^-(x;\sigma^2,\widetilde{\boldsymbol{\eta}}_{(m)}) & = \sum_{j=1}^n\Omega(j;\boldsymbol{\eta}_{(m)})\left(\Phi\left(-\frac{x}{\sigma}\right) - e^{{(\sigma\eta_j)^2}/{2}+x\eta_j} \Phi\left(\frac{-x}{\sigma} - \sigma\eta_j \right)\right), \end{align}

(4.12) \begin{align} \mathrm{V}(x;\sigma^2,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)}) & = \sum_{i=1}^n \sum_{j=1}^m \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)}) \nonumber\\ & \quad \times \left[\frac{\tilde{\eta}_j}{\eta_i + \tilde{\eta}_i}\left(\Phi\left(-\frac{x}{\sigma}\right) + e^{\frac{(\sigma\eta_i)^2}{2}-x\eta_i}\Phi\left(\frac{x}{\sigma} - \sigma\eta_i \right)\right)\right.\nonumber\\ & \quad \left. + \frac{\eta_i}{\eta_i + \tilde{\eta}_j}\left(\Phi\left(-\frac{x}{\sigma}\right) - e^{\frac{(\sigma\tilde{\eta}_j)^2}{2}+x\tilde{\eta}_j} \Phi\left(\frac{-x}{\sigma} - \sigma\tilde{\eta}_j \right)\right)\right]. \end{align}

Next, to deal with the expected value $\mathrm {E}[S_T {1}_{\{S_T\ge K\}}]$ in (4.9), we need to further define the following three functions related to the NHE$^+$, NHE$^-$, and NDHE distributions.

\begin{align*} \bullet \quad & \mathrm{\Pi}^+(x;\sigma^2,\boldsymbol{\eta}_{(n)})= \mathrm{E}[e^{Z + H_{(n)}}{1}_{\{Z + H_{(n)} \geq x\}}],\\ \bullet \quad & \mathrm{\Pi}^-(x;\sigma^2,\widetilde{\boldsymbol{\eta}}_{(m)})= \mathrm{E}[e^{Z - \widetilde{H}_{(m)}}{1}_{\{Z - \widetilde{H}_{(m)} \geq x\}}],\\ \bullet \quad & \mathrm{\Lambda}(x;\sigma^2,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)}) =\mathrm{E}[e^{Z + H_{(n)} - \widetilde{H}_{(m)}}{1}_{\{Z + H_{(n)} - \widetilde{H}_{(m)} \geq x\}}]. \end{align*}

Fortunately, once again we are able to derive their closed-form results, as seen subsequently.

Lemma 8. Suppose that $\eta _i \gt 1$, $i=1,\ldots,n$. The closed-form expressions of the three functions $\mathrm {\Pi }^+, \mathrm {\Pi }^-, \mathrm {\Lambda }$ as defined above are respectively given by

(4.13) \begin{align} \mathrm{\Pi}^+(x;\sigma^2,\boldsymbol{\eta}_{(n)}) & = \sum_{i=1}^n\Omega(i;\boldsymbol{\eta}_{(n)})\frac{\eta_i e^{\frac{(\sigma\eta_i)^2}{2}}}{\eta_i-1} \nonumber\\ & \quad \times \left(e^{-{(\eta_i^2-1)\sigma^2}/{2}}\Phi\left(\frac{-x}{\sigma}+\sigma\right) + e^{-(\eta_i-1)x}\Phi\left(\frac{x}{\sigma}-\sigma\eta_i\right)\right), \end{align}

(4.14) \begin{align} \mathrm{\Pi}^-(x;\sigma^2,\widetilde{\boldsymbol{\eta}}_{(m)}) & = \sum_{j=1}^m\Omega(j;\widetilde{\boldsymbol{\eta}}_{(m)})\,\frac{\tilde{\eta}_j e^{\frac{(\sigma\tilde{\eta}_j)^2}{2}}}{\tilde{\eta}_j+1} \nonumber\\ & \quad \times \left(e^{-{(\tilde{\eta}_j^2-1)\sigma^2}/{2}}\Phi\left(\frac{-x}{\sigma}+\sigma\right) - e^{(\tilde{\eta}_j+1)x}\Phi\left(\frac{-x}{\sigma}-\sigma\tilde{\eta}_j\right)\right), \end{align}

(4.15) \begin{align} \mathrm{\Lambda}(x;\sigma^2,\boldsymbol{\eta}_{(n)}, \widetilde{\boldsymbol{\eta}}_{(m)}) & = \sum_{i=1}^n \sum_{j=1}^m \Omega (i;\boldsymbol{\eta}_{(n)})\Omega (j;\boldsymbol{\tilde{\eta}}_{(m)})\frac{\eta_i\widetilde{\eta_j}}{\eta_i + \tilde{\eta}_j} \nonumber\\ & \quad \times{\times}\left[\frac{ e^{{(\sigma\eta_i)^2}/{2}}}{\eta_i-1}\left(e^{-{(\eta_i^2-1)\sigma^2}/{2}}\Phi \left(\frac{-x}{\sigma}+\sigma\right) + e^{-(\eta_i-1)x}\Phi\left(\frac{x}{\sigma}-\sigma\eta_i\right)\right) \right. \nonumber\\ & \quad + \left. \frac{e^{{(\sigma\tilde{\eta}_j)^2}/{2}}}{\tilde{\eta}_j+1} \left(e^{-{(\tilde{\eta}_j^2-1)\sigma^2}/{2}}\Phi\left(\frac{-x}{\sigma}+\sigma\right) - e^{(\tilde{\eta}_j+1)x}\Phi\left(\frac{-x}{\sigma}-\sigma\tilde{\eta}_j\right)\right)\right]. \end{align}

Armed with Lemmas 7 and 8, we are now ready to derive the European call option price under the proposed NDEJD model.

Proposition 9. The price of a European call option with strike price $K$ and maturity time $T$ under the NDEJD model is given by

(4.16) \begin{equation} C = S_0 e^{-(q + {\sigma^2}/{2} + \chi)T} \,\mathrm{\Upsilon}_1(\alpha) - K e^{{-}rT} \mathrm{\Upsilon}_2(\alpha), \end{equation}

where $\alpha = \ln (K/S_0)-\nu _T$ and $\mathrm {\Upsilon }_1(\alpha )$, $\mathrm {\Upsilon }_2(\alpha )$ are respectively given by

(4.17) \begin{align} \mathrm{\Upsilon}_1(\alpha) & = e^{{\sigma^2T}/{2}-\lambda T}\Phi\left(\frac{\sigma^2T - \alpha}{\sigma\sqrt{T}}\right) \nonumber\\ & \quad + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda T) [p^n \Pi^+(\alpha;\sigma^2 T,\boldsymbol{\eta}_{(n)}) + (1-p)^n \Pi^-(\alpha;\sigma^2 T,\widetilde{\boldsymbol{\eta}}_{(n)})] \nonumber\\ & \quad + \sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda T) \left[\binom{n}{m}p^m(1-p)^{n-m} \Lambda(\alpha;\sigma^2 T,\boldsymbol{\eta}_{(n)}, \widetilde{\boldsymbol{\eta}}_{(m)})\right], \end{align}

(4.18) \begin{align} \mathrm{\Upsilon}_2(\alpha) & = e^{-\lambda T}\Phi\left(-\frac{\alpha}{\sigma\sqrt{T}}\right)\nonumber\\ & \quad + \sum_{n=1}^{n^*} f_{\mathsf{SP}}(n;\lambda T)[p^n \mathrm{U}^+(\alpha;\sigma^2 T,\boldsymbol{\eta}_{(n)}) + (1-p)^n \mathrm{U}^-(\alpha;\sigma^2 T,\widetilde{\boldsymbol{\eta}}_{(n)})] \nonumber\\ & \quad + \sum_{n=2}^{n^*} \sum_{m=1}^{n-1} f_{\mathsf{SP}}(n;\lambda T) \left[\binom{n}{m}p^m(1-p)^{n-m} \mathrm{V}(\alpha;\sigma^2 T,\boldsymbol{\eta}_{(n)},\widetilde{\boldsymbol{\eta}}_{(m)})\right]. \end{align}

Proof. It is clear from (4.9) together with (4.3) that

$$C = S_0 e^{-(q + {\sigma^2}/{2} + \chi)T}\,\mathrm{E}[e^{\sigma W_T + X_T}{1}_{\{\sigma W_T + X_T\ge\alpha\}}] - Ke^{{-}rT}\mathrm{P}(\sigma W_T + X_T\ge \alpha),$$

where we note that the event $\{S_T\ge K\}$ is equivalent to the event $\{\sigma W_T + X_T\ge \alpha \}$. Let us define the two functions:

$$\mathrm{\Upsilon}_1(\alpha) = \mathrm{E}[e^{\sigma W_T+X_T}{1}_{\{\sigma W_T+X_T\geq \alpha\}}],\quad \mathrm{\Upsilon}_2(\alpha) = \Pr(\sigma W_T+X_T\geq \alpha).$$

Following the unconditioning technique on $N_t$ and $M_t$ as in the proof of Proposition 6, the two functions can be derived as claimed using the results in Lemmas 7 and 8.

5.1. Return distributions

We first present the numerical results for the risk-neutral return distributions under the NDEJD model. The basic parameters that are not related to jumps (about drift and diffusion) are $r = 0.05$, $q = 0.01$, and $\sigma = 0.2$. Some parameters about jumps are fixed: $\lambda = 5$, $n^\ast = 10$, $\eta _{0} = 10$, and $\tilde {\eta }_{0} = 5$ ($\tilde {\eta }_{0} \lt \eta _0$ indicates that the left tail is fatter than the right tail, as commonly seen empirically). In order to investigate the influences from the growing and diminishing effects of jumps, we consider the following parameters about jump size distributions: $p\in \{0.2, 0.5,0.8\}$, $\gamma = \{0.8, 0.9, 1.0, 1.1, 1.2\}$. Three maturity times are considered: $T\in \{0.25, 0.5, 1.0\}$.

Figure 6 shows the risk-neutral distribution of $R_T$. The effect of maturity time $T$ can be seen from top to bottom. As $T$ grows, the distribution becomes more widely spread as expected. The effect of $p$ (the probability of seeing positive jumps) can be seen from left to right. Larger $p$ indicates that positive jumps are more likely to happen than negative jumps, and this makes the right tail fatter. The differences between left and right tails become more prominent when $T$ gets longer. In each plot, five curves are provided to show the effect of $\gamma$. We observe that when jumps have growing effects ($\gamma \lt 1$), the distribution tends to be fatter-tailed and the left tail is generally (see the plots with $p=0.5$) lifted up more because the influences from negative jumps are stronger than those from positive jumps (this is because of our chosen parameters for the initial jump, i.e., $\tilde {\eta }_0 \lt \eta _0$). On the contrary, if the effects from jumps are diminishing ($\gamma \gt 1$), the distribution becomes more narrowly spread with thinner tails on both sides. When $p = 0.2$, the left tail becomes much fatter because negative jumps happen more frequently than positive jumps. As $p$ gets larger, the right tails are lifted up more. This shows that positive jumps bring about more impacts to the distribution.

Figure 6. Return distributions (pdf) under various combinations of $p$ and $T$.

The basic statistics (mean, standard deviation, skewness, kurtosis) of $R_T$ are provided in Table 1. The results for Kou's (2002) DEJD model (where jumps are driven by a Poisson process) are given below the results for the benchmark NDEJD model with $\gamma = 1$ (where jumps are driven by a stopped Poisson process with $n^\ast = 10$) in order to show that the two models are almost identical. We observe close agreements bewteen Figure 6 and Table 1, but the figures in Table 1 provide more insights. We see that as $\gamma$ gets smaller, kurtosis becomes larger but skewness becomes more negative when $p=0.2$ (left tails receive more impacts) and less negative when $p = 0.8$ (right tails receive more impacts). Moreover, it is worth pointing out that, as $T$ gets longer, the standard deviation (SD) increases as expected, but both the skewness (in absolute value) and kurtosis seem to become smaller, indicating that the distribution becomes less skewed and less fat-tailed. In contrast, however, as seen in Figure 6, the influences on both tails appear to be more prominent as $T$ gets longer. This shows that the statistical figures in Table 1 provide a different view which is not directly observable from Figure 6.

Table 1. The basic statistics of $R_T$ under various combinations of $p$ and $T$.

5.2. Option prices

Next, we present the numerical results about European call option prices. To investigate the nonlinear influences from $\gamma$ and other factors on call option prices, we conduct a regression analysis. Let $M = K/S_0$ denote the moneyness for a call option (it is at-the-money if $M=1$). We consider option contracts with different moneyness $M = 0.5, 0.6,\ldots, 1.5$ and maturity $T = 0.25, 0.5, 0.75, 1.0$. Take $\gamma = 1$ (seen as almost equivalent to the DEJD model) as the benchmark and define “Bias” for the $\gamma \ne 1$ cases as the percentage error of the call option price relative to the benchmark case, that is,

$$\text{Bias} = \frac{C_\gamma-C_{\gamma=1}}{C_{\gamma=1}}.$$

We expect that Bias is a nonlinear function of $\gamma -1$, $p, M, T$, written as $\text {Bias} = f(\gamma -1, p, M, T)$. In our analysis, we intend to use the following quadratic formulation for the function $f$, that is,

(5.1) \begin{equation} \text{Bias} = \beta_0(p,M,T) + \beta_1(p,M,T) (\gamma-1) + \beta_2(p,M,T) (\gamma-1)^2, \end{equation}

to model the dependence of Bias on the main factor $(\gamma -1)$. The contributions from $p, M, T$ are left to the coefficients $\beta _i(p,M,T), i=0,1,2$. Specifically, the linear, quadratic, and interactive influences that come purely from $p, M, T$ are captured by $\beta _0(p,M,T)$, whereas their interactions with $(\gamma -1)$ and $(\gamma -1)^2$ are respectively captured by $\beta _1(p,M,T)$ and $\beta _2(p,M,T)$.

Table 2 presents the estimation results of eight nested regression models based on (5.1). Model 1 is a purely linear model which takes the following form ($c_i$ represents the coefficients)

$$\text{Bias} = c_0 + c_1 p + c_2 M + c_3 T + c_4 (\gamma-1),\quad R^2 = 32\%.$$

Such a linear model provides a baseline level of adjusted $R^2 = 32\%$. Model 2 adds the cross terms between $p,M,T$ and $(\gamma -1)$ to model their interactions (these cross terms are part of $\beta _1(p,M,T)$ in (5.1)), that is,

$$\text{Bias} = c_0 + c_1 p + c_2 M + c_3 T + (c_4 + c_5 p + c_6 M + c_7 T)(\gamma-1), ~R^2 = 59\%,$$

lifting the adjusted $R^2$ to 59%. By adding quadratic and interaction terms between $p,M,T$ to $\beta _1(p,M,T)$, the adjusted $R^2$ is improved to 62% in Model 3 and further to 70% in Model 4. The stepwise inclusions of all higher-order interaction terms increase the adjusted $R^2$ to nearly 90% in Model 8. This is our full model where $\beta _i(p,M,T)$, $i=0,1,2$ exhibit the following structure:

\begin{align*} \beta_0 & = (c_0 + c_1 p + c_2 M + c_3 T), \\ \beta_1 & = (c_4 + c_5 p + c_6 M + c_7 T + c_8 p^2 + c_9 M^2 + c_{10} T^2 + c_{11} pM + c_{12} pT + c_{13} MT), \\ \beta_2 & = (c_{14} + c_{15} p + c_{16} M + c_{17} T + c_{18} p^2 + c_{19} M^2 + c_{20} T^2 + c_{21} pM + c_{22} pT + c_{23} MT). \end{align*}

In Table 2, the significant coefficients (marked with $\ast$) of the nonlinear terms show that the Bias depends nonlinearly on its influencing factors (i.e., $\gamma$ as well as $p,M,T$). The significance of the coefficients in $\beta _1$ (i.e., $c_4,\ldots,c_{13}$, which are almost identical from Model 4 to Model 8) reveals that the main factor $\gamma$ may cause strong Bias (percentage error) in call price if the nonhomogeneity in jumps is not considered. The significance of the coefficients in $\beta _2$ (i.e., $c_{14},\ldots,c_{23}$) further shows that the main factor $\gamma$ influences Bias in an asymmetric way. For example, from the positiveness of $c_{14} = 3.693$ (Model 5), we see that a decrease in $\gamma$ (e.g., $1\to 0.9$) would cause a greater bias than an increase in $\gamma$ (e.g., $1\to 1.1$). The fact that many of $c_{14},\ldots,c_{23}$ (Model 8) are significant indicates that $\beta _2$ depends in a complicated nonlinear way on $p, M, T$.

Table 2. Estimation results of the nested regression models for “Bias”.

$^\ast$Indicates that the coefficient is significant with $p$-value <0.05.

5.3. Implied volatility smiles

We move on to present our numerical results on implied volatiliy smiles which reveal the risk structure behind the model. Suppose that the market follows the proposed NDEJD model, the implied volatility $\hat {\sigma }$ is the particular volatility in the Black–Scholes formula that gives the same option price, that is,

$$C_{\text{NDEJD}} = C_{\text{BS}}(S_0, K, T, r, q, \hat{\sigma}),$$

where $C_{\text {NDEJD}}$ is obtained from our pricing formula (4.16). The volatility smile ($\hat {\sigma }$ as a function of moneyness, which is redefined as $M = \ln (K/S_0)$) is then obtained by solving the above nonlinear equation for all $M$. (See [Reference Cont and Fonseca4] for its implication and [Reference Dumas, Fleming and Whaley6] for an example of its empirical analysis.)

The results are presented in Figure 7 where the at-the-money ($M=0$) implied volatility is fixed at 0.4 (by setting $\lambda = 2$ with $\sigma$ left free; other parameters remain the same as in Figure 6). We see that the smile curves are pulled higher (lower) when $\gamma$ decreases (increases), indicating that the level of risk is heightened (lowered) when there are growing (diminishing) impacts from jumps. The smiles for $p=0.5$ tend to be more symmetric than $p = 0.2$ or 0.8. When $p = 0.2$ (0.8), more influences are made on the left (right) tails of the return distributions, and thus the left (right) wings of the smiles are lifted or lowered more. Moreover, the smiles exhibit the highest curvature when $T = 0.25$ (which is consistent with Table 1 where higher kurtosis is seen in smaller $T$ cases). As $T$ increases, the smiles become less curved, but the differences across all $\gamma$ become more noticeable.

Figure 7. Implied volatility smiles under various combinations of $p$ and $T$.

In contrast with the above qualitative observations, below we give a quantitative analysis on how the shapes of the IV smiles vary with $\gamma$ and other parameters. To this end, we consider a typical IV smiles as depicted in Figure 8, where ${\rm IV}_1$, ${\rm IV}_2$, and ${\rm IV}_3$ are the IV for three moneyness $M_1, M_2, M_3$ ($M_2$ refers to the lowest point and $M_1 = M_2-1$, $M_3 = M_2+1$). The following five parameters are used to capture the main features of a typical smile curve:

• • $a = M_2$: the moneyness where the smile drops to the lowest point,

• • $b = {\rm IV}_2$: the minimal level of IV at $M_2$,

• • $c = {\rm IV}_3 - 2{\rm IV}_2 + {\rm IV}_1$: the curvature of the smile,

• • $d = {\rm IV}_2 - {\rm IV}_1$: the (negative) slope of the left wing,

• • $e = {\rm IV}_3-{\rm IV}_2$: the (positive) slope of the right wing.

Figure 8. The features of a typical smile curve.

The five shape parameters of all IV smiles in Figure 7 are summarized in Table 3, from which we may take a closer inspection of the smile curves. It is seen from the results of parameter $a$ that the lowest point of smiles moves leftward as $p$ gets large but it moves rightward as $T$ gets large. From the results of parameter $b$, we observe that smaller $\gamma$ leads to higher $b$ (the overall volatility level is also pulled up, as seen in every plot). The shape parameters $c$, $d$, $e$ carry more important information about the risk-neutral distributions. We observe that as either $\gamma$ or $T$ increases, the curvature $c$ decreases, indicating that the distribution becomes less fat-tailed (this is in close agreement with Table 1). The parameters $d$ and $e$ are concerned with the slopes of the left and right wings. It is observed that smaller $\gamma$ makes the two wings steeper, and therefore, the values of $|d|$ and $e$ become greater. But the parameter $p$ tends to influence one side of the smile. Smaller $p$ makes the left wing steeper, whereas greater $p$ makes the right wing steeper. All these observations are consistent with Figure 6 and Table 1.

Table 3. The five shape parameters $a, b, c, d, e$ of IV smiles.

5.4. Empirical analysis

To further justify the proposed NDEJD model which exhibits a richer shape of the IV smile relative to the DEJD model as shown in the preceding subsection, we provide an empirical analysis. The standard S&P 500 index option (SPX) data are collected and the two competing models are calibrated to their implied volatility smiles for various times to maturity. The data for the whole year of 2021 are collected and analyzed; among them, eight data sets (see Table 4) are picked for our presentation. The chosen two observation dates are the very first trading days in the first half year (sets (a)(b)(c)(d)) and in the second half year (sets (e)(f)(g)(h)). For each observation day, four maturity days are considered in order to cover different times to maturity ($\tau$). Table 4 also shows the initial index ($S_0$), the maximal and minimal strike prices ($K_{\rm min}$, $K_{\rm max}$), and the number of strike prices ($n$) for each maturity date.

Table 4. Eight data sets for the comparison of model calibration.

With the interest rate and dividend yield chosen as $r = 0.01, q = 0$, we calibrate the two models to determine the model parameter set $\boldsymbol {x}$ by minimizing the sum of squared errors (SSE) of the implied volatility curve, that is,

$$\min_{\boldsymbol{x}} \sum_{i=1}^n (\text{IV}_{\mathrm{model}}(\boldsymbol{x}, K_i) - \text{IV}_{\mathrm{market}}(K_i) )^2,$$

where the parameter set is $\boldsymbol {x} = (\sigma,\lambda,\eta _0,\tilde {\eta }_0,p)$ for DEJD, and is $\boldsymbol {x} = (\sigma,\lambda,\eta _0,\tilde {\eta }_0,p,\gamma )$ for NDEJD. The results are given in Table 5, where we observe that the NDEJD parameter $\gamma$ may deviate from 1 (it could be greater or less than 1). The inclusion of $\gamma$ provides one more degree of freedom which reduces the SSE significantly, indicating the advantage of NDEJD over DEJD for these particular data sets.

Table 5. Results of model calibration.

Figure 9 provides a graphical illustration of the calibration results from the eight data sets. In each plot, the horizontal axis is moneyness $M$ which is defined to be $M = \ln (K/S_0)$. While DEJD provides satisfactory results, it is not difficult to see that our NDEJD model provides an even closer fit to the market IV smiles. In particular, the deviations between the two models are more prominent at the left ends (with large negative $M$), showing that the flexibility in NDEJD provides advantages for the options that are far-away-from-the-money. This gives a justification for the proposed NDEJD model.

Figure 9. Graphical illustration of the calibration results in Table 5.