Option Pricing under a Generalized Black–Scholes Model with Stochastic Interest Rates, Stochastic Strings, and Lévy Jumps

Abstract

We introduce a novel option pricing model that features stochastic interest rates along with an underlying price process driven by stochastic string shocks combined with pure jump Lévy processes. Substituting the Brownian motion in the Black–Scholes model with a stochastic string leads to a class of option pricing models with expiration-dependent volatility. Further extending this Generalized Black–Scholes (GBS) model by adding Lévy jumps to the returns generating processes results in a new framework generalizing all exponential Lévy models. We derive four distinct versions of the model, with each case featuring a different jump process: the finite activity lognormal and double–exponential jump diffusions, as well as the infinite activity CGMY process and generalized hyperbolic Lévy motion. In each case, we obtain closed or semi-closed form expressions for European call option prices which generalize the results obtained for the original models. Empirically, we evaluate the performance of our model against the skews of S&P 500 call options, considering three distinct volatility regimes. Our findings indicate that: (a) model performance is enhanced with the inclusion of jumps; (b) the GBS plus jumps model outperform the alternative models with the same jumps; (c) the GBS-CGMY jump model offers the best fit across volatility regimes.

1. Introduction

It is well documented that the Black–Scholes option pricing formula ([1]) is founded on several assumptions which, while simplifying, are frequently misaligned with empirical evidence. For instance, the constant risk-free interest rate assumption is not supported by the observed dynamics of market interest rates (enhancements to include stochastic interest rates in option pricing models are developed in [2,3]). Similarly, the assumption that underlying security prices follow geometric Brownian motions with constant volatility is challenged by the time- and strike-dependent volatilities implied by the Black–Scholes formula using market option prices.

However, financial practice has devised an accommodation in the form of the volatility surface that allows the Black–Scholes formula to maintain a central role as a real-world modeling tool. The volatility surface for options on a given underlying security is a three-dimensional graph of the values of implied volatility, varying in the dimensions of moneyness and the expiration of the options ([4]). Given a volatility surface, if we consider a specific expiration, we obtain a graph of implied volatility as a function of moneyness, which turns out to be a convex curve called volatility smile or skew.

Subsequent studies have produced option pricing models that avoid certain simplifying assumptions of the Black–Scholes model in order to explain the behavior of volatility skews. Some of them introduce stochastic volatility ([5,6,7]), others include jumps ([8,9,10,11,12]), and a few models present both features ([13,14]).

Several studies in the financial econometrics literature have refined methodologies to disentangle jumps from continuous volatility in price trajectories, yielding robust evidence of jump discontinuities in equity prices ([15,16]). Indeed, one explanation for a non-constant volatility surface is that the Black–Scholes-implied volatilities include not only continuous volatility but also the effects of omitted risk factors such as jumps in the prices of underlying securities.

In [14], Bakshi, Cao, and Chen (BCC) derived an option pricing formula for a model that features stochastic volatility, stochastic interest rates, and random jumps in the stock price process. The BCC model nests many of the closed-form option pricing formulas in the literature and thus serves well as a benchmark for the comparative evaluation of option pricing models. Their findings suggest that the configuration featuring both stochastic volatility and jumps (the so-called SVJ) outperforms all other specifications nested in the BCC framework, with the traditional Black–Scholes model ranking least favorably.

In a recent paper, Ref. [17] introduces a novel generalization of the Black–Scholes model incorporating stochastic interest rates and stochastic string shocks ([18,19]). Termed the generalized Black and Scholes (GBS) model, this new framework yields a Black–Scholes-like formula for European call option prices in which volatility is a function of the expiration times of underlying options. Such a specification offers an explanation for the variability of implied volatilities in the dimension of option expiry, the so-called term structure of implied volatilities.

The effect of option expiration times on underlying price dynamics is motivated in [17] by the “pin risk” phenomenon, in which underlying volatility is accentuated in proximity to the strike price as expiration looms. Recent empirical research on the nascent zero-day-to-expiry options market sheds additional light on this mechanism. In [20], evidence is presented that an increment in the trading volume of zero-day-to-expiry options results in an increase in the volatility of the underlying. This finding lends support to the existence of a feedback mechanism whereby trading in options affects the underlying volatility. The introduction of stochastic string shocks in the GBS model allows for the expiration times of options to be explicitly incorporated in the underlying dynamics.

Given this backdrop, our paper contributes to the debate on the necessity of stochastic volatility in option pricing models to account for observed option skews. Could a model combining expiration-dependent volatility with random jumps in the underlying price outperform a model with stochastic volatility and random jumps? To investigate this question, we extend the GBS model by incorporating pure jump Lévy processes into the stochastic string-driven price dynamics, thereby relaxing the three most contentious assumptions of the Black–Scholes model: constant volatility, static interest rates, and continuous sample paths.

Note that our model cannot be embedded into the framework of [14], since they use the Brownian motion as the source of continuous risk in the price processes, whereas our model is based on stochastic string processes, which are more general and, in fact, subsume Brownian motion as a particular case. Moreover, our model is not limited to compound Poisson jumps with lognormal distributions; instead, it is versatile enough to accommodate any pure jump Lévy process.

In order to determine which type of jump best captures observed discontinuities in observed security prices, we study four different examples of pure-jump Lévy processes used in the option pricing literature: compound Poisson processes with lognormally (LN) distributed jumps ([8]) and with a double-exponential (DE) distribution of jumps ([9]), generalized hyperbolic Lévy motion (GHLM) ([11]), and CGMY ([10]). First, we obtain closed or semi-closed form expressions for European call option prices under each type of jump, and verify the coherence with their established Brownian and no-jump benchmarks. Subsequently, we conduct the in-sample empirical comparisons of our extended GBS model against a suite of alternative models, both with and without jump components. The competing models are Black–Scholes and Heston, and their respective extensions include LN jumps, which correspond, respectively, to [8] and to the SVJ model in [14]. We also include in the comparison the pure jump models used. To check the robustness of our approach, we perform the comparison under three distinct historical volatility regimes.

Our results show that the GBS model augmented with CGMY jumps outperforms the alternative models (with or without jumps), except in the high volatility scenario in which the Heston model with jumps (SVJ) performs better, albeit at the cost of estimating one additional parameter. Nevertheless, when we compare the SVJ model with a GBS plus CGMY model with the same number of parameters, we find that the latter provides the best fit.

The remainder of the paper is structured as follows. Section presents the mathematical prerequisites for the rest of this paper: stochastic string and Lévy processes, and the GBS model. Section 3 explores the mathematical characteristics of the processes we are working with, which are special semi-martingales, but not Lévy processes in general. Section 4 leads to the general expression for the underlying dynamics under the GBS model with jumps. Section 5 contains the results related to option pricing with the GBS model with the four different types of jumps. Section 6 develops the main characteristics and results of the models alternative to GBS in the empirical tests. Section 7 includes the calibration results of all competing models. Section 8 summarizes and concludes.

2. Preliminaries

In this section, we recall definitions and briefly summarize relevant facts about the GBS model, stochastic strings, and Lévy processes that we will use in this paper. For stochastic string processes, the main references are [18,19]. The GBS model was introduced in [17]. For Lévy processes, there are many excellent references, including [21] (Chapter 1).

Let $(\Omega ,\mathcal{F},\mathbb{F},\mathbb{P})$ be a filtered probability space. The specific form of the filtration $\mathbb{F}={\left({\mathcal{F}}_t\right)}_{0\le t\le \mathcal{Y}}$ will be determined later, and we consider that $\mathcal{F}={\mathcal{F}}_{\mathcal{Y}}$, where $\mathcal{Y}$ is the finite time horizon. We will assume the existence of risk-free security, $B\left(t\right)$, whose evolution satisfies $dB\left(t\right)=r_{t}B\left(t\right)dt$, $B\left(0\right)=1$, where $r_t$ is the stochastic risk-free rate. In order to preclude arbitrage opportunities, we will also assume the existence of a probability measure $\mathbb{Q}$, equivalent to $\mathbb{P}$, such that asset prices discounted with $B\left(t\right)$ are $\mathbb{Q}$-martingales. We will call this measure $\mathbb{Q}$ the equivalent martingale measure.

2.1. Stochastic String Processes

A stochastic string process is an infinite-dimensional stochastic process (or random field) consisting of a continuum of adapted stochastic processes $Z\left(·,x,\omega \right)$ indexed by a (second) time parameter, x. Specifically,

$$\begin{array}{cccc}Z:\hfill & {\Delta }^2\times \Omega \hfill & \to \hfill & \mathbb{R}\hfill \\ \hfill & \left(t,x,\omega \right)\hfill & \mapsto \hfill & Z\left(t,x,\omega \right)\hfill \end{array}$$

where ${\Delta }^2=\left\{\left(t,x\right)\in {\mathbb{R}}^2:0\le t\le \mathcal{Y},x\ge 0\right\}$. Henceforth, we will drop the explicit mention of the dependence on $\omega$.

The stochastic string process Z is assumed to satisfy the following properties:

(a)

The stochastic processes $Z\left(·,x\right)$ and $Z\left(t,·\right)$ are continuous for each $x\ge 0$ and for each $t\in \left[0,\mathcal{Y}\right]$, respectively.

(b)

The process $Z\left(·,x\right)$ is a martingale for each $x\ge 0$.

(c)

The process $Z\left(t,·\right)$ is differentiable for each $t\in \left[0,\mathcal{Y}\right]$.

(d)

For each $x,y\ge 0$, it is the case that

$$d{\left[Z\left(·,x\right),Z\left(·,y\right)\right]}_t=c\left(t,x,y\right)dt$$

(1)

where $c\left(t,x,y\right)$ is an admissible, continuous, and differentiable correlation function for each t. Here, admissibility means that for each t, $c\left(t,·,·\right)$ is symmetric, positive semidefinite, and satisfies $\left|c\left(t,x,y\right)\right|\le 1$ and $c\left(t,x,x\right)=1,\forall x,y\ge 0$ ([18]).

Properties (a), (b), and (d), jointly with Lévy’s Theorem guarantee that, for each $x\ge 0$, the process $Z\left(·,x\right)$ is a Brownian motion.

2.2. The Generalized Black–Scholes Model

Stochastic string processes are used in [17] to obtain a generalization of the Black–Scholes model that accounts for volatilities dependent on option expirations and stochastic interest rates. Concretely, the no-arbitrage dynamics for the discounted underlying, ${\overline{S}}_t=B^{-1}\left(t\right)S_t$ is given by

$$\frac{d{\overline{S}}_t}{{\overline{S}}_t}={\int }_0^{\infty }\sigma \left(t,u\right)d\tilde{Z}\left(t,u\right)du$$

(2)

where $\sigma \left(t,u\right)$ is the stochastic string volatility and $\tilde{Z}\left(t,u\right)$ is a stochastic string process under $\mathbb{Q}$. The parameter u should be interpreted as representing times until the expiration of options written on the underlying stock. In the case in which there are no options written on the stock, we have $\sigma \left(t,u\right)=\sigma \left(t\right)\delta \left(u\right)$, where $\delta (·)$ is the Dirac delta, recovering Merton’s model [2]. If we go further and make $r_t=r$ and $\sigma \left(t\right)=\sigma$, Equation (2) reduces to the Black–Scholes dynamics [1].

From Equation (2), the explicit expression for the discounted stock price at time $T>t$ is obtained as

$${\overline{S}}_T={\overline{S}}_{t}exp\left\{{\int }_{s=t}^T{\int }_{u=0}^{\infty }\sigma \left(s,u\right)d\tilde{Z}\left(s,u\right)du-\frac{1}{2}{\int }_{s=t}^T{\int }_{y=0}^{\infty }{\int }_{u=0}^{\infty }{R}_s\left(u,y\right)dudyds\right\}$$

where we have defined $R_s\left(u,y\right)\equiv \sigma \left(s,u\right)\sigma \left(s,y\right)c\left(s,u,y\right)$.

Taking into account the definitions

$$\begin{array}{cc}{Z}_{t,T}^{\mathbb{Q}}\hfill & \equiv {\int }_{s=t}^T{\int }_{u=0}^{\infty }\sigma \left(s,u\right)d\tilde{Z}\left(s,u\right)du\hfill \\ R_{t,T}\hfill & \equiv {\int }_{s=t}^T{\int }_{y=0}^{\infty }{\int }_{u=0}^{\infty }{R}_s\left(u,y\right)dudyds\hfill \\ Y_{t,T}\hfill & \equiv Z_{t,T}^{\mathbb{Q}}-\frac{1}{2}{R}_{t,T}\hfill \end{array}$$

(3)

the previous expression for ${\overline{S}}_T$ can be rewritten in a more compact form as

$${\overline{S}}_T={\overline{S}}_{t}exp\left\{Y_{t,T}\right\}$$

(4)

In order to obtain a closed-form expression for a call option price, [17] introduces the additional assumptions of Gaussian interest rates, the deterministic term-structure of volatilities and correlations, and the independence of interest rates and discounted stock prices. Under these assumptions, we have that $R_{t,T}$ is deterministic, and that

$$Z_{t,T}^{\mathbb{Q}}\sim N\left(0,R_{t,T}\right)$$

(5)

and

$$Y_{t,T}\sim N\left(-\frac{1}{2}{R}_{t,T},R_{t,T}\right)$$

under $\mathbb{Q}$ on ${\mathcal{F}}_t$.

With these assumptions, the time t price of a European call option with maturity T and strike K on the underlying S, $C_{GBS}\left(t,T,K\right)$ is obtained as

$$C_{GBS}\left(t,T,K\right)=S_t\Phi \left(d_1\right)-KP\left(t,T\right)\Phi \left(d_2\right)$$

(6)

where $P\left(t,T\right)$ is the price, at time t, of a zero-coupon bond maturing at T, and

$$\begin{array}{ccc}\hfill d_1& =& \frac{1}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\left[ln\left(\frac{{\overline{S}}_t}{K}\right)+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}\right]\hfill \\ \hfill d_2& =& d_1-\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\hfill \end{array}$$

with

$$\begin{array}{ccc}\hfill {\mu }_{t,T}& =& {\mathbb{E}}^{\mathbb{Q}}\left\{{\int }_0^T{r}_{s}ds\left|{\mathcal{F}}_t\right.\right\}\hfill \\ \hfill {\sigma }_{t,T}^2& =& Va{r}^{\mathbb{Q}}\left\{{\int }_0^T{r}_{s}ds\left|{\mathcal{F}}_t\right.\right\}\hfill \end{array}$$

Under constant interest rates, $r_t=r$, we have the same Equation (6), but with $d_1$ and $d_2$ replaced by ${\overline{d}}_1$ and ${\overline{d}}_2$ given by

$$\begin{array}{ccc}\hfill {\overline{d}}_1& =& \frac{1}{\sqrt{R_{t,T}}}\left[ln\left(\frac{S_t}{K}\right)+r\left(T-t\right)+\frac{1}{2}{R}_{t,T}\right]\hfill \\ \hfill {\overline{d}}_2& =& {\overline{d}}_1-\sqrt{R_{t,T}}\hfill \end{array}$$

Moreover, assuming that $R_{t,T}$ can be written as a power series of $T-t$, the call option price can be obtained by modifying the Black–Scholes option pricing formula so that the constant volatility parameter is replaced by the following deterministic function for the term-structure of volatility

$${\sigma }_{TS}\left(T-t\right)=\sqrt{{\sigma }^2+a_2\left(T-t\right)+a_3{\left(T-t\right)}^2+\cdots }$$

(7)

where $T-t$ is the time until option expiration, and $\sigma >0$ and the $a_j$ are constants.

2.3. Lévy Processes

A real-valued stochastic process, $L_t$, that is adapted to the filtration of the probability space is a Lévy process if

• $\mathbb{P}\{L_0=0\}=1$;
• For $n\ge 1$ and $0\le t_0<t_1<\cdots <t_n$, the random variables $L_{t_0}$, $L_{t_1}-L_{t_0}$, $\dots ,L_{t_n}-L_{t_{n-1}}$ are independent (independent increments);
• For any $\epsilon >0$, ${lim}_{t\to 0}\mathbb{P}\left\{\right|L_{s+t}-L_s|>\epsilon \}=0$ for all $s\ge 0$ (stochastic continuity);
• Sample paths are right-continuous with left-limits (RCLL).

A number of interesting facts about the sample path behavior can be deduced from the definition of a Lévy process. For instance, the stochastic continuity implies that the probability of a jump at any fixed time in the future is 0. On any interval of finite length $[0,t]$, the property that paths are RCLL implies that (1) the number of jumps is at most countably infinite, and (2) for fixed $\epsilon >0$, there can be only finitely many jumps on $[0,t]$ with magnitudes greater than $\epsilon$.

A Lévy process can be decomposed as the sum of a linear drift, a Brownian motion, and a jump process. The Lévy triple $(\gamma ,\sigma ,\nu )$ corresponding to these three parts, with $\gamma$ being the drift rate, $\sigma$ being the diffusion coefficient, and $\nu$ being a measure with support on $\mathbb{R}\setminus \left\{0\right\}$ such that ${\int }_{\mathbb{R}}{\left(\right|x|}^2\wedge 1)\nu \left(dx\right)<\infty$, called a Lévy measure. Given a Lévy measure $\nu$, it follows that for all $\epsilon >0$, $\nu (\mathbb{R}\setminus (-\epsilon ,\epsilon \left)\right)<\infty$. Thus, a probability measure can be defined by restricting $\nu$ and normalizing so that total mass is 1, that is,

$${\mathbb{P}}_{\epsilon }\left(dx\right)=\frac{\nu \left(dx\right)}{\nu (\mathbb{R}\setminus (-\epsilon ,\epsilon \left)\right)}{1}_{\left\{\right|x|\ge \epsilon \}}.$$

If ${lim}_{\epsilon \to 0}\nu (\mathbb{R}\setminus (-\epsilon ,\epsilon ))=\lambda <\infty$, then $\nu \left(\mathbb{R}\right)<\infty$ (since $\nu$ has no mass at 0), and the Lévy process is said to have finite activity since sample paths of $L_t$ have only finitely many jumps on any closed interval $[0,T]$. In this case, the jump part of $L_t$ is necessarily a compound Poisson process and $\nu \left(dx\right)=\lambda F\left(dx\right)$, where $\lambda$ is the arrival rate and F is the distribution function of the jump size. Alternatively, if $\nu$ is singular in some neighborhood of 0, then $\nu \left(\mathbb{R}\right)=\infty$ and the Lévy process is said to have infinite activity since the sample paths of $L_t$ almost surely have a countably infinite number of jumps on any closed interval $[0,T]$.

By the Lévy-Khintchine Theorem, $L_t$ has characteristic exponent given by

$${\Psi }_{L_1}\left(u\right)=i\gamma u-\frac{{\sigma }^2}{2}{u}^2+{\int }_{\mathbb{R}}(e^{iux}-1-iuh\left(x\right))\nu \left(dx\right),$$

where $h\left(x\right)$ is a truncation function which is required in some cases when ${\int }_{\mathbb{R}}\nu \left(dx\right)=\infty$, depending on the behavior of $\nu$ near 0. The characteristic function of a Lévy process can be defined in terms of its characteristic exponent, $E\left(e^{iu{L}_t}\right)=e^{t{\Psi }_{L_1}\left(u\right)}$. Thus, a Lévy process is completely determined by its characteristic exponent.

Now, we specialize to the case of a Lévy process, $L_t$, with Lévy triple $(\gamma ,0,\nu )$ (no Brownian motion part). The Lévy–Itô decomposition of $L_t$ can be written as

$$L_t=\gamma t+{\int }_0^t{\int }_{\left|x\right|\ge 1}xN(ds,dx)+{\int }_0^t{\int }_{\left|x\right|<1}x(N(ds,dx)-\nu \left(dx\right)ds),$$

(8)

where N is a Poisson random measure and $\mu \otimes \nu (dt,dx)=\nu \left(dx\right)dt$ is the product measure of Lebesgue measure, $\mu$, and the Lévy measure, $\nu$. The first jump term represents large jumps (arbitrarily defined as being $\ge 1$ in magnitude) and the second jump term represents small compensated jumps (with jumps of size 0 disallowed). The representation of $L_t$ in (8) corresponds to the truncation function $h\left(x\right)=x{1}_{\left\{\right|x|<1\}}$.

The Lévy processes that we consider in this paper have finite first moments, that is, ${\int }_{\left|x\right|>1}\left|x\right|\nu \left(dx\right)<\infty$. In this case, the Lévy–Itô decomposition of $L_t$ can be written as

$$L_t=\gamma t+{\int }_0^t{\int }_{\mathbb{R}}x(N(ds,dx)-\nu \left(dx\right)ds),$$

(9)

corresponding to $h\left(x\right)=x$. Here, the jump part with the compensator is a martingale. The integral in (9) cannot be split into two terms unless the jumps have finite variation, which is implied by ${\int }_{0<\left|x\right|<1}\left|x\right|\nu \left(dx\right)<\infty$. In this case, a truncation function is not needed for the process to be well defined. Thus, we can choose $h\left(x\right)=x$, $h\left(x\right)=x{1}_{\left\{\right|x|<c\}}$, $c\in {\mathbb{R}}^{+}$, or $h\left(x\right)=0$. In the case of $h\left(x\right)=0$, the jump process is uncompensated and (9) becomes

$$L_t=\gamma t+{\int }_0^t{\int }_{\mathbb{R}}xN(ds,dx).$$

(10)

The $\sigma$ and Lévy measure $\nu$ are invariant for a given Lévy process, and once an appropriate truncation function is chosen, the Lévy triple $({\gamma }_h,\sigma ,\nu )$ is unique. In this paper, we will follow the convention of offsetting through ${\gamma }_h$ any changes to the linear drift due to compensating the jump process. That is, when we give Lévy triples for specific Lévy processes, the linear drift coefficient will be ${\gamma }_h=\gamma +{\int }_{\mathbb{R}}h\left(x\right)\nu \left(dx\right)$.

3. Combining Stochastic Strings and Jump Processes

A stochastic process $X_t$ on $(\Omega ,\mathcal{F},{\left({\mathcal{F}}_t\right)}_{t\ge 0},\mathbb{P})$ is called a semimartingale if it can be written

$$X_t=A_t+M_t,$$

(11)

where $A_t$ is a process with RCLL sample paths of finite variation and $M_t$ is a local martingale. (a property is said to hold locally for a stochastic process X if there exists a sequence of stopping times ${\left({\tau }_n\right)}_{n\in \mathbb{N}}$ such that, for all n, the stopped process $X_{t\wedge {\tau }_n}$ has the property). A locally integrable semimartingale is called a special semimartingale, in which case a decomposition as in (11) exists with $A_t$ being a predictable process of finite variation. Moreover, the local martingale $M_t$ can be written as a sum of a continuous local martingale and a purely discontinuous local martingale. Thus, a special semimartingale has the decomposition

$$X_t=X_0+A_t+X_t^c+X_t^d,$$

(12)

where $X_0^c=M_0=0$, $A_t$ is a predictable process of finite variation, $X_t^c$ is a continuous local martingale, and $X_t^d$ is a local martingale orthogonal to all continuous local martingales.

Lévy processes are semimartingales and it is proven in [19] that stochastic integrals with respect to stochastic string processes are continuous martingales. Thus, if $X_t$ is the sum of a stochastic string integral and a Lévy process with Lévy triple $(\gamma ,0,\nu )$ and finite first moment, then it is a semimartingale. Moreover, such a process is a special semimartingale since it has decomposition as in (12) with

$$A_t=\gamma t,X_t^c={\int }_0^t{\int }_{{\mathbb{R}}^{+}}\sigma (s,x)dZ(s,x)dx,X_t^d={\int }_0^t{\int }_{\mathbb{R}}x(N(ds,dx)-\nu \left(dx\right)ds).$$

(13)

Special semimartingales with decompositions as in (13) are the focus of this paper. It is important to note that $X_t$ is not a Lévy process, in general, because the term with the stochastic integral does not appear in the Lévy–Itô decomposition. It only becomes a Lévy process in the particular case in which the stochastic string integral becomes a stochastic integral with respect to Brownian motion (see Section 2.2 ). Therefore, the results of our model can be seen as generalizations of all the (exponential) Lévy models.

It will be useful to state an Itô’s lemma for such special semimartingales. We first recall Itô’s lemma for a general one-dimensional semimartingale. This result can be found in [21] (Chapter 1).

Theorem 1.

Let $X_t$ be a semimartingale and $f:\mathbb{R}\to \mathbb{R}$ a ${\mathcal{C}}^2$ function. Then, the stochastic process $f\left(X_t\right)$ is also a semimartingale and

$$\begin{array}{cc}\hfill f\left(X_t\right)& =f\left(X_0\right)+{\int }_0^t{f}^{\prime }\left(X_{s-}\right)d{X}_s+\frac{1}{2}{\int }_0^t{f}^{\prime \prime }\left(X_{s-}\right)d{\langle X^c\rangle }_s\hfill \\ \hfill & +\sum _{s\le t}\left(f\left(X_s\right)-f\left(X_{s-}\right)-f^{\prime }\left(X_{s-}\right)\Delta X_s\right).\hfill \end{array}$$

(14)

To describe the dynamics of a one-dimensional process driven by a stochastic string, the inherent multidimensionality of the stochastic string must be recognized. For the semimartingale decomposition above, with

$$X_t^c={\int }_0^t{\int }_0^{\infty }\sigma (s,x)dZ(s,x)dx,$$

we have

$$d{\langle X^c\rangle }_t={\int }_0^{\infty }{\int }_0^{\infty }\sigma (t,x)\sigma (t,y)d{[Z(·,u),Z(·,y)]}_{t}dudy.$$

In the case of the stochastic string plus a Lévy process with finite first moment, it can be shown that (14) becomes

$$\begin{array}{cc}\hfill f\left(X_t\right)& =f\left(X_0\right)+{\int }_0^t\gamma f^{\prime }\left(X_{s-}\right)ds+{\int }_0^t{f}^{\prime }\left(X_{s-}\right){\int }_0^{\infty }\sigma (s,u)Z(s,u)duds\hfill \\ \hfill & +\frac{1}{2}{\int }_0^t{f}^{\prime \prime }\left(X_{s-}\right){\int }_0^{\infty }{\int }_0^{\infty }{R}_s\left(u,y\right)dudyds\hfill \\ \hfill & +{\int }_0^t{\int }_{-\infty }^{\infty }\left(f(X_{s-}+x)-f\left(X_{s-}\right)\right)(N(ds,dx)-\nu \left(dx\right)ds)\hfill \end{array}$$

4. Option Pricing with the Generalized Black and Scholes Model with Jumps

In this section, we specify the general form for the stochastic processes we will use to model underlying price dynamics. We posit two independent random processes in the underlying security price consisting of a stochastic string process $Z(t,x)$ and a Lévy process $L_t$. The instantaneous interest rate $r_t$ is driven by a standard Brownian motion that is assumed to be independent of both the stochastic string and the Lévy process. The available information at any time $t\ge 0$ is given by the filtration

$${\mathcal{F}}_t=\sigma \left\{Z(s,x),L_s,B_s:0\le s\le t,x\ge 0\right\},t\le \mathcal{Y}.$$

We start by including the Lévy process in Expression (4) for the discounted underlying price of the GBS model in the form

$${\overline{S}}_T={\overline{S}}_{t}exp\left\{Y_{t,T}+L_{t,T}+{\alpha }_{t,T}\right\}$$

(15)

where $L_{t,T}\equiv L_T-L_t$ and ${\alpha }_{t,T}$ is a deterministic term which accounts for the martingale property of ${\overline{S}}_T$. Applying this property, we have

$${\overline{S}}_t={\mathbb{E}}^{\mathbb{Q}}\left\{{\overline{S}}_T\left|{\mathcal{F}}_t\right.\right\}={\overline{S}}_t{e}^{{\alpha }_{t,T}}{\mathbb{E}}^{\mathbb{Q}}\left\{e^{Y_{t,T}}\left|{\mathcal{F}}_t\right.\right\}{\mathbb{E}}^{\mathbb{Q}}\left\{e^{L_{t,T}}\left|{\mathcal{F}}_t\right.\right\}={\overline{S}}_t{e}^{{\alpha }_{t,T}}{e}^{\left(T-t\right){\Psi }_{L_1}\left(-i\right)}$$

from where we obtain ${\alpha }_{t,T}=-\left(T-t\right){\Psi }_{L_1}\left(-i\right)$. Replacing this value in (15), we obtain

$${\overline{S}}_T={\overline{S}}_{t}exp\left\{Y_{t,T}-\left(T-t\right){\Psi }_{L_1}(-i)+L_{t,T}\right\}$$

which can be rewritten as

$$S_T=S_{t}exp\left\{{\int }_t^T{r}_{s}ds+Y_{t,T}-\left(T-t\right){\Psi }_{L_1}(-i)+L_{t,T}\right\}$$

(16)

Equation (16) is the general expression for the price dynamics in our model.

5. Jump Processes

In this section, we present four specific cases of the Lévy jumps, $L_t$ in order to obtain different option pricing models. The first two are finite-activity Lévy processes, that is, compound Poisson processes, with lognormal ([8]) and double-exponential ([9]) distributions of jumps, respectively. The other two jump processes are Lévy processes of infinite activity: the generalized hyperbolic Lévy motion ([11]) and the CGMY process ([10]).

We will defer to Section 7 the numerical calculation of option prices through the characteristic function of Gaussian string plus jump dynamics.

5.1. Compound Poisson Process with Lognormal Jumps

Consider a compound Poisson process, $L_t={\sum }_{i=1}^{N_t}{J}_i$, where $N_t$ is a Poisson process with parameter $\lambda$ (the jump intensity) and $J_i$ are the jump sizes. The compound Poisson process has sample paths of bounded variation and its characteristic function is given by ([22], Proposition 3.4)

$${\phi }_t^{CP}\left(u\right)=exp\left\{t\lambda {\int }_{-\infty }^{+\infty }\left(e^{iux}-1\right)f\left(dx\right)\right\}$$

where f is the jump size distribution. In our case, as in [8], the jump sizes are i.i.d. with $J_i\sim N\left(m,{\delta }^2\right)$, and thus, we have

$${\Psi }_{L_1}\left(-i\right)=\lambda {\int }_{-\infty }^{+\infty }\left(e^x-1\right)e^{-\frac{{\left(x-m\right)}^2}{2{\delta }^2}}dx=\lambda \left(e^{m+\frac{{\delta }^2}{2}}-1\right)$$

Replacing these expressions for the compound Poisson process in (16), we obtain

$$S_T=S_{t}exp\left\{{\int }_t^T{r}_{s}ds+Y_{t,T}-\lambda \left(T-t\right)\left(exp\left(m+\frac{{\delta }^2}{2}\right)-1\right)+\sum _{i=N_t+1}^{N_T}{J}_i\right\}$$

where $N_t$, the Poisson process that counts the jumps of $S_t$, is independent from r and Z, and the jump sizes, $J_i$, are also independent from r, Z, and N.

Theorem 2.

The price, $C_{GBS}^{LN}$, of a European call option with payoff ${\left[S_T-K\right]}_{+}$, is given by

$$C_{GBS}^{LN}(t,S_t)=e^{-\lambda \left(T-t\right)}\sum _{n=0}^{\infty }\frac{{\left[\lambda \left(T-t\right)\right]}^n}{n!}{C}_{GBS}\left(t,T,S_t^{LN},R_{t,T}+n{\delta }^2\right)$$

where $C_{GBS}\left(t,T,S_t,R_{t,T}\right)$ is the call option price of the GBS model, Equation (6), and

$$S_t^{LN}\equiv S_{t}exp\left[nm+\frac{n{\delta }^2}{2}-\lambda \left(T-t\right)\left(exp\left(m+\frac{{\delta }^2}{2}\right)-1\right)\right]$$

Proof. 

By risk-neutral valuation, we have

$$C_{GBS}^{LN}(t,S_t)={\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\left[S_T-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}$$

Conditioning on the number of jumps between t and T, $N_T-N_t=N_{T-t}$, we can write the option value as follows

$$\begin{array}{ccc}& & C_{GBS}^{LN}(t,S_t)\hfill \\ & =& \sum _{n=0}^{\infty }\mathbb{Q}\left(N_{T-t}=n\right){\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\left[{\tilde{S}}_{t}exp\left(Y_{t,T}+\sum _{i=1}^n{J}_i\right)-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& \sum _{n=0}^{\infty }\frac{e^{-\lambda \left(T-t\right)}{\left[\lambda \left(T-t\right)\right]}^n}{n!}B\left(t\right){\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_0^T{r}_{s}ds}{\left[{\tilde{S}}_{t}exp\left(Y_{t,T}+\sum _{i=1}^n{J}_i\right)-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}\hfill \end{array}$$

where we have defined ${\tilde{S}}_t\equiv S_{t}exp\left[{\int }_t^T{r}_{s}ds-\lambda \left(T-t\right)\left(e^{m+\frac{{\delta }^2}{2}}-1\right)\right]$. Taking into account that ${\int }_0^T{r}_{s}ds\sim N\left({\mu }_{t,T},{\sigma }_{t,T}^2\right)$ and $Y_{t,T}+{\sum }_{i=1}^n{J}_i\sim N\left(-\frac{1}{2}{R}_{t,T}+nm,R_{t,T}+n{\delta }^2\right)$ on ${\mathcal{F}}_t$ under $\mathbb{Q}$, and using standard calculations with lognormal random variables, we have

$$\begin{array}{ccc}& & B\left(t\right){\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_0^T{r}_{s}ds}{\left[{\tilde{S}}_{t}exp\left(Y_{t,T}+\sum _{i=1}^n{J}_i\right)-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& S_t^{LN}\Phi \left(\frac{ln\frac{{\overline{S}}_t^{LN}}{K}+{\mu }_{t,T}+\frac{1}{2}\left(R_{t,T}+n{\delta }^2\right)}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}+n{\delta }^2}}\right)\hfill \\ & & -KP\left(t,T\right)\Phi \left(\frac{ln\frac{{\overline{S}}_t^{LN}}{K}+{\mu }_{t,T}-\frac{1}{2}(R_{t,T}+n{\delta }^2)-{\sigma }_{t,T}^2}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}+n{\delta }^2}}\right)\hfill \end{array}$$

where we have taking into account that $P\left(t,T\right)=B\left(t\right)e^{-{\mu }_{t,T}+\frac{{\sigma }_{t,T}^2}{2}}$ ([17], proof of Theorem 1). □

The following corollary provides the consistency of the previous theorem with the models in [8,17].

Corollary 1.

In the cases of no jumps and Black–Scholes price dynamics with constant interest rates, the expression of Theorem 2 reduces to the call option price in the models of [8,17], respectively.

Proof. 

The no-jumps case can be obtained by making $\lambda \to 0$, from where we obtain $\mathbb{Q}\left(N_{T-t}=n\right)=\frac{e^{-\lambda \left(T-t\right)}{\left[\lambda \left(T-t\right)\right]}^n}{n!}\to {\delta }_{n,0}$, $S_t^{LN}\to S_t$, and then $C_{GBS}^{LN}\left(t,S_t\right)\to C_{GBS}\left(t,T,S_t,R_{t,T}\right)$. The Black–Scholes case with jumps correspond to making the substitutions $r_t\to r$, $\sigma \left(t,u\right)\to \sigma \delta \left(u\right)$ and $c(t,u,y)\to 1$. From these changes, we obtain $R_{t,T}\to {\sigma }^2\left(T-t\right)$, ${\sigma }_{t,T}\to 0$, ${\mu }_{t,T}\to rT$ and $P\left(t,T\right)\to e^{-r\left(T-t\right)}$. Replacing them in the expression of Theorem 2, we arrive at the call option price in the [8] model. □

5.2. Compound Poisson Process with Double-Exponential Jumps

In [9], a finite activity jump-diffusion model is presented in which the density function of the jump size is given by

$$f_{DE}\left(x\right)=p{\theta }_1{e}^{-{\theta }_1}x{1}_{x<0}+(1-p){\theta }_2{e}^{{\theta }_{2}x}{1}_{x>0}$$

In this case, the characteristic exponent of the jump process, $L_t$, is given by

$$\begin{array}{ccc}\hfill {\Psi }_{L_1}\left(u\right)& =& \lambda {\int }_{-\infty }^{+\infty }\left(e^{iux}-1\right)f_{DE}\left(dx\right)\hfill \\ & =& \lambda \left[p{\theta }_1{\int }_0^{\infty }\left(e^{iux}-1\right)e^{-{\theta }_{1}x}dx+\left(1-p\right){\theta }_2{\int }_{-\infty }^0\left(e^{iux}-1\right)e^{{\theta }_{2}x}dx\right]\hfill \\ & =& \lambda \left[p{\theta }_1\left(\frac{{\theta }_1+iu}{{\theta }_1^2+u^2}-\frac{1}{{\theta }_1}\right)+\left(1-p\right){\theta }_2\left(\frac{{\theta }_2-iu}{{\theta }_2^2+u^2}-\frac{1}{{\theta }_2}\right)\right]\hfill \end{array}$$

from which it follows that

$${\Psi }_{L_1}\left(-i\right)=\lambda \left[\frac{p{\theta }_1}{{\theta }_1-1}+\frac{\left(1-p\right){\theta }_2}{{\theta }_2+1}-1\right]$$

Using this expression in Equation (16), we obtain

$$\begin{array}{cc}{S}_T\hfill & =S_{t}exp\left\{{\int }_t^T{r}_{s}ds+Y_{t,T}-\lambda \left(T-t\right)\left[\frac{p{\theta }_1}{{\theta }_1-1}+\frac{\left(1-p\right){\theta }_2}{{\theta }_2+1}-1\right]+L_{t,T}\right\}\hfill \\ \hfill & =S_t^{DE}exp\left[{\int }_t^T{r}_{s}ds+Y_{t,T}+L_{t,T}\right]\hfill \end{array}$$

(17)

where

$$S_t^{DE}\equiv S_{t}exp\left\{-\lambda \left(T-t\right)\left[\frac{p{\theta }_1}{{\theta }_1-1}+\frac{\left(1-p\right){\theta }_2}{{\theta }_2+1}-1\right]\right\}$$

Theorem 3.

The price, $C_{GBS}^{DE}$, of a European call option with payoff ${\left[S_T-K\right]}_{+}$, is given by

$$C_{GBS}^{DE}\left(t,S_t\right)=S_t{\Xi }^{+}-KP\left(t,T\right){\Xi }^{-}$$

with

$$\begin{array}{ccc}\hfill {\Xi }^{\pm }& =& e^{-\lambda \left(T-t\right)}\Phi \left(\frac{ln\left(\frac{{\overline{S}}_t^{DE}}{K}\right)\pm \frac{1}{2}{R}_{t,T}+{\mu }_{t,T}\pm x^{\pm }}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\right)\hfill \\ & & +\frac{exp\left(\frac{{\theta }_1^2\left({\sigma }_{t,T}^2+R_{t,T}\right)}{2}\right)}{\sqrt{2\pi \left({\sigma }_{t,T}^2+R_{t,T}\right)}}\sum _{n=1}^{\infty }\frac{e^{-\lambda \left(T-t\right)}{\left[\lambda \left(T-t\right)\right]}^n}{n!}\sum _{k=1}^n{P}_{n,k}{\left({\theta }_1\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\right)}^k\hfill \\ & & \times I_{k-1}\left(ln\left(\frac{K}{{\overline{S}}_t^{DE}}\right)\mp \frac{1}{2}{R}_{t,T}-{\mu }_{t,T}\mp x^{\mp };-{\theta }_1,-\frac{1}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}},-{\theta }_1\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\right)\hfill \\ & & +\frac{exp\left(\frac{{\theta }_2^2\left({\sigma }_{t,T}^2+R_{t,T}\right)}{2}\right)}{\sqrt{2\pi \left({\sigma }_{t,T}^2+R_{t,T}\right)}}\sum _{n=1}^{\infty }\frac{e^{-\lambda \left(T-t\right)}{\left[\lambda \left(T-t\right)\right]}^n}{n!}\sum _{k=1}^n{Q}_{n,k}{\left({\theta }_2\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\right)}^k\hfill \\ & & \times I_{k-1}\left(ln\left(\frac{K}{{\overline{S}}_t^{DE}}\right)\mp \frac{1}{2}{R}_{t,T}-{\mu }_{t,T}\mp x^{\mp };{\theta }_2,\frac{1}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}},-{\theta }_2\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\right)\hfill \end{array}$$

where $x^{+}=0$, $x^{-}={\sigma }_{t,T}^2$, and

$$\begin{array}{ccc}\hfill P_{n,k}& =& \sum _{i=k}^{n-1}\left(\genfrac{}{}{0pt}{}{n-k-1}{i-k}\right)\left(\genfrac{}{}{0pt}{}{n}{i}\right){\left(\frac{{\theta }_1}{{\theta }_1+{\theta }_2}\right)}^{i-k}{\left(\frac{{\theta }_2}{{\theta }_1+{\theta }_2}\right)}^{n-i}{p}^i{\left(1-p\right)}^{n-i}\hfill \\ \hfill Q_{n,k}& =& \sum _{i=k}^{n-1}\left(\genfrac{}{}{0pt}{}{n-k-1}{i-k}\right)\left(\genfrac{}{}{0pt}{}{n}{i}\right){\left(\frac{{\theta }_1}{{\theta }_1+{\theta }_2}\right)}^{n-i}{\left(\frac{{\theta }_2}{{\theta }_1+{\theta }_2}\right)}^{i-k}{p}^{n-i}{\left(1-p\right)}^i\hfill \end{array}$$

for $1\le k\le n-1$, $P_{n,n}=p^n$, $Q_{n,n}={\left(1-p\right)}^n$, $\left(\genfrac{}{}{0pt}{}{0}{0}\right)\equiv 1$, and for all $n\ge -1$

$$\begin{array}{ccc}\hfill I_n\left(c;\alpha ,\beta ,\delta \right)& =& -\frac{e^{\alpha c}}{\alpha }\sum _{i=0}^n{\left(\frac{\beta }{\alpha }\right)}^{n-i}H{h}_i\left(\beta c-\delta \right)\hfill \\ & & +{\left(\frac{\beta }{\alpha }\right)}^{n+1}\frac{\sqrt{2\pi }}{\beta }{e}^{\frac{\alpha \delta }{\beta }+\frac{{\alpha }^2}{2{\beta }^2}}\Phi \left(-\beta c+\delta +\frac{\alpha }{\beta }\right)\hfill \end{array}$$

if $\beta >0$ and $\alpha \ne 0$, and

$$\begin{array}{ccc}\hfill I_n\left(c;\alpha ,\beta ,\delta \right)& =& -\frac{e^{\alpha c}}{\alpha }\sum _{i=0}^n{\left(\frac{\beta }{\alpha }\right)}^{n-i}H{h}_i\left(\beta c-\delta \right)\hfill \\ & & -{\left(\frac{\beta }{\alpha }\right)}^{n+1}\frac{\sqrt{2\pi }}{\beta }{e}^{\frac{\alpha \delta }{\beta }+\frac{{\alpha }^2}{2{\beta }^2}}\Phi \left(\beta c-\delta -\frac{\alpha }{\beta }\right)\hfill \end{array}$$

if $\beta <0$ and $\alpha <0$, where $H{h}_n$ is defined for $n\ge 0$ by

$$H{h}_n=\frac{1}{n!}{\int }_x^{\infty }{\left(t-x\right)}^n{e}^{-\frac{t^2}{2}}dt$$

and Φ is the cumulative distribution function of a standard normal random variable.

Proof. 

The risk-neutral price of the option is given by

$$\begin{array}{ccc}\hfill C_{GBS}^{DE}\left(t,S_t\right)& =& {\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\left[S_T-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& {\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}-K{\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \end{array}$$

In order to obtain the first expectation, we need to work under the measure ${\mathbb{Q}}_1$, defined by the Radon–Nikodym derivative $\frac{d\mathbb{Q}}{d{\mathbb{Q}}_1}=\frac{S_0{e}^{{\int }_0^T{r}_{s}ds}}{S_T}$. By Bayes formula, we have

$$\begin{array}{cc}{\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill & =\frac{{\mathbb{E}}^{{\mathbb{Q}}_1}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\frac{d\mathbb{Q}}{d{\mathbb{Q}}_1}\left|{\mathcal{F}}_t\right.\right\}}{{\mathbb{E}}^{{\mathbb{Q}}_1}\left\{\frac{d\mathbb{Q}}{d{\mathbb{Q}}_1}\left|{\mathcal{F}}_t\right.\right\}}\hfill \\ \hfill & =\frac{e^{{\int }_0^t{r}_{s}ds}{\mathbb{Q}}_1\left\{S_T>K\left|{\mathcal{F}}_t\right.\right\}}{{\mathbb{E}}^{{\mathbb{Q}}_1}\left\{{\overline{S}}_T^{-1}\left|{\mathcal{F}}_t\right.\right\}}\hfill \end{array}$$

(18)

Applying the Bayes formula again and using the martingale property of $\mathbb{Q}$, it is not difficult to show that the inverse of the discounted price is a martingale under ${\mathbb{Q}}_1$, namely ${\mathbb{E}}^{{\mathbb{Q}}_1}\left\{{\overline{S}}_T^{-1}\left|{\mathcal{F}}_t\right.\right\}={\overline{S}}_t^{-1}$. Replacing this in (18), we obtain

$${\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}=S_t{\mathbb{Q}}_1\left\{S_T>K\left|{\mathcal{F}}_t\right.\right\}$$

Now, to obtain the conditional probability, we need to work with expression (17) under ${\mathbb{Q}}_1$. The stochastic string shock under a measure ${\mathbb{Q}}_i$ equivalent to $\mathbb{Q}$, $d{Z}^{{\mathbb{Q}}_i}\left(t,u\right)$, is given by

$$d{Z}^{{\mathbb{Q}}_i}\left(t,u\right)=d\tilde{Z}\left(t,u\right)-dt{\int }_0^{\infty }c\left(t,u,y\right)h\left(t,y\right)dy$$

with $h\left(t,y\right)$ being a function that is determined by ([19], Lemma 4.8). Imposing that ${\overline{S}}_t^{-1}$ is a ${\mathbb{Q}}_1$–martingale, it can be shown that $h\left(t,y\right)=\sigma \left(t,y\right)$ and then

$$d{Z}^{{\mathbb{Q}}_1}\left(t,u\right)=d\tilde{Z}\left(t,u\right)-dt{\int }_0^{\infty }c\left(t,u,y\right)\sigma \left(t,y\right)dy$$

With this change, we can write $Y_{t,T}=Z_{t,T}^{{\mathbb{Q}}_1}+\frac{1}{2}{R}_{t,T}$, with $Z_{t,T}^{{\mathbb{Q}}_1}\equiv {\int }_t^T{\int }_0^{\infty }\sigma \left(s,u\right)d{Z}^{{\mathbb{Q}}_1}\left(s,u\right)du$, and Equation (17) reads

$$\begin{array}{ccc}\hfill S_T& =& S_t^{DE}exp\left[{\int }_t^T{r}_{s}ds+Z_{t,T}^{{\mathbb{Q}}_1}+\frac{1}{2}{R}_{t,T}+L_{t,T}\right]\hfill \\ & =& {\overline{S}}_t^{DE}exp\left[{\int }_0^T{r}_{s}ds+Z_{t,T}^{{\mathbb{Q}}_1}+\frac{1}{2}{R}_{t,T}+L_{t,T}\right]\hfill \end{array}$$

Applying the Bayes rule and the independence assumption, it can be shown that the conditional mean and variance of ${\int }_0^T{r}_{s}ds$ do not change under the measure ${\mathbb{Q}}_1$. Moreover, under this same measure, $Z_{t,T}^{{\mathbb{Q}}_1}\sim N\left(0,R_{t,T}\right)$. Thus, applying Proposition B.3 and Theorem B.1 in [9], we obtain ${\mathbb{Q}}_1\left\{S_T>K\left|{\mathcal{F}}_t\right.\right\}={\Xi }^{+}$.

For the second expectation, taking into account the properties of the T-forward measure, ${\mathbb{Q}}_T$, we have

$${\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}=P\left(t,T\right){\mathbb{Q}}_T\left\{S_T>K\left|{\mathcal{F}}_t\right.\right\}$$

In order to obtain this last probability, we need to know the conditional distribution of $S_T$ under ${\mathbb{Q}}_T$. From [19], we know that $d{\tilde{Z}}^{{\mathbb{Q}}_T}\left(t,y\right)=d\tilde{Z}\left(t,y\right)+dt{\int }_0^{T-t}c\left(t,u,y\right)\sigma \left(t,u\right)du$, where $d{\tilde{Z}}^{{\mathbb{Q}}_T}\left(t,y\right)$ is the stochastic string shock under the T-forward measure. Replacing this equality in Equation (3), and then in (17), we obtain

$$S_T={\overline{S}}_t^{DE}exp\left[{\int }_0^T{r}_{s}ds+Z_{t,T}^{{\mathbb{Q}}_T}-\frac{1}{2}{R}_{t,T}-{\int }_{s=t}^T{\int }_{u=0}^{\infty }{\int }_{y=0}^{T-t}{R}_t\left(u,y\right)dudyds+L_{t,T}\right]$$

(19)

where $Z_{t,T}^{{\mathbb{Q}}_T}\equiv {\int }_{s=t}^T{\int }_{u=0}^{\infty }\sigma \left(s,u\right)d{\tilde{Z}}^{{\mathbb{Q}}_T}\left(s,u\right)du$. Now, we need to study the behavior under ${\mathbb{Q}}_T$ of the terms in the exponential of (19), except for $-\frac{1}{2}{R}_{t,T}$, which is deterministic, and $L_{t,T}$, whose distribution is assumed to be the same, independent of the (equivalent) measure. Taking into account that $\sigma \left(s,u\right)$ is deterministic, we have that $Z_{t,T}^{{\mathbb{Q}}_T}\sim N\left(0,R_{t,T}\right)$ under ${\mathbb{Q}}_T$ on ${\mathcal{F}}_t$. In order to obtain the distribution of ${\int }_0^T{r}_{s}ds$ under that measure, we need to elaborate some more.

The dynamics of the discounted bond price in the stochastic string framework is given by ([19])

$$\begin{array}{ccc}\hfill \frac{d\overline{P}\left(t,T\right)}{\overline{P}\left(t,T\right)}& =& -{\int }_{y=0}^{T-t}\sigma \left(t,y\right)d\tilde{Z}\left(t,y\right)dy=\hfill \\ & =& \left[{\int }_{y=0}^{T-t}{\int }_{u=0}^{T-t}{R}_t\left(u,y\right)dudy\right]dt-{\int }_{y=0}^{T-t}\sigma \left(t,y\right)d{\tilde{Z}}^{{\mathbb{Q}}_T}\left(t,y\right)dy\hfill \end{array}$$

whose solution is

$$\begin{array}{cc}\frac{\overline{P}\left(t,T\right)}{\overline{P}\left(0,T\right)}\hfill & =exp\left\{-\frac{1}{2}{\int }_{s=0}^t{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds\right.\hfill \\ \hfill & \left.-{\int }_{s=0}^t{\int }_{y=0}^{T-s}\sigma \left(s,y\right)d\tilde{Z}\left(s,y\right)dy\right\}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & =exp\left\{\frac{1}{2}{\int }_{s=0}^t{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds\right.\hfill \\ \hfill & \left.-{\int }_{s=0}^t{\int }_{y=0}^{T-s}\sigma \left(s,y\right)d{\tilde{Z}}^{{\mathbb{Q}}_T}\left(s,y\right)dy\right\}\hfill \end{array}$$

(21)

With these preliminaries, we can now explicitly find that ${\mu }_{t,T}$, and the conditional expectation of ${\int }_0^T{r}_{s}ds$ under $\mathbb{Q}$, as

$$\begin{array}{ccc}\hfill {\mu }_{t,T}& =& {\mathbb{E}}^{\mathbb{Q}}\left\{{\int }_0^T{r}_{s}ds\left|{\mathcal{F}}_t\right.\right\}=-{\mathbb{E}}^{\mathbb{Q}}\left\{ln\overline{P}\left(T,T\right)\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& ln\overline{P}\left(t,T\right)+\frac{1}{2}{\int }_{s=t}^T{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds\hfill \end{array}$$

where we have used Equation (20). On the other hand, from [17], we know that, for Gaussian models, we have $P\left(t,T\right)=B_t{e}^{-{\mu }_{t,T}+\frac{1}{2}{\sigma }_{t,T}^2}$. Solving for ${\mu }_{t,T}$ and replacing in the previous equation, we can identify ${\sigma }_{t,T}^2={\int }_{s=t}^T{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds$, and write

$${\mu }_{t,T}=ln{\overline{P}}^{-1}\left(t,T\right)+\frac{1}{2}{\sigma }_{t,T}^2$$

Performing a similar calculation using Equation (21), we arrive at

$$\begin{array}{ccc}\hfill {\mathbb{E}}^{{\mathbb{Q}}_T}\left\{{\int }_0^T{r}_{s}ds\left|{\mathcal{F}}_t\right.\right\}& =& ln{\overline{P}}^{-1}\left(t,T\right)-\frac{1}{2}{\int }_{s=t}^T{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds\hfill \\ & =& {\mu }_{t,T}-{\sigma }_{t,T}^2\hfill \end{array}$$

Taking into account that, under our framework, conditional variances remain the same under a change of measure, we conclude that ${\int }_0^T{r}_{s}ds\sim N\left({\mu }_{t,T}-{\sigma }_{t,T}^2,{\sigma }_{t,T}^2\right)$ under ${\mathbb{Q}}_T$ on ${\mathcal{F}}_t$. The remaining term with the triple integral in (19) is equal to zero. To see this, considering (20) and the expression for $Y_{t,T}$, we can write

$$\begin{array}{ccc}\hfill {\int }_0^T{r}_{s}ds& =& -ln\overline{P}\left(0,T\right)+\frac{1}{2}{\int }_{s=t}^T{\int }_{y=0}^{T-s}{\int }_{u=0}^{T-s}{R}_s\left(u,y\right)dudyds\hfill \\ & & +{\int }_{s=0}^T{\int }_{y=0}^{T-s}\sigma \left(s,y\right)d\tilde{Z}\left(s,y\right)dy\hfill \end{array}$$

and

$$Y_{t,T}\equiv Z_{t,T}^{\mathbb{Q}}-\frac{1}{2}{R}_{t,T}={\int }_{s=t}^T{\int }_{y=0}^{\infty }\sigma \left(s,y\right)d\tilde{Z}\left(s,y\right)dy-\frac{1}{2}{R}_{t,T}$$

Taking conditional covariance and applying the conditional Itô isometry, after some calculations, we arrive at

$$Co{v}^{\mathbb{Q}}\left\{{\int }_0^T{r}_{s}ds,Y_{t,T}\left|{\mathcal{F}}_t\right.\right\}={\int }_{s=t}^T{\int }_{u=0}^{\infty }{\int }_{y=0}^{T-t}{R}_t\left(u,y\right)dudyds$$

which is zero by our independence assumption.

In summary, we have

$$S_T={\overline{S}}_t^{DE}exp\left[{\int }_0^T{r}_{s}ds+Z_{t,T}^{{\mathbb{Q}}_T}-\frac{1}{2}{R}_{t,T}+L_{t,T}\right]$$

(22)

and ${\int }_0^T{r}_{s}ds+Z_{t,T}^{{\mathbb{Q}}_T}-\frac{1}{2}{R}_{t,T}\sim N\left({\mu }_{t,T}-\frac{1}{2}{R}_{t,T}-{\sigma }_{t,T}^2,{\sigma }_{t,T}^2+R_{t,T}\right)$. Now, again applying Proposition B.3 and Theorem B.1 in [9], we arrive at $\mathbb{Q}\left\{S_T>K\left|{\mathcal{F}}_t\right.\right\}={\Xi }^{-}$, which concludes the proof. □

As in the case of lognormal jumps, we have the following consistency result.

Corollary 2.

In the cases of no jumps and Black–Scholes price dynamics with constant interest rates, the expression of Theorem 3 reduces to the call option price in the models of [9,17], respectively.

Proof. 

Working as in the proof of Corollary 1, in the no-jump case ($\lambda \to 0$), we obtain $\mathbb{Q}\left(N_{T-t}=n\right)=\frac{e^{-\lambda \left(T-t\right)}{\left[\lambda \left(T-t\right)\right]}^n}{n!}\to {\delta }_{n,0}$, ${\overline{S}}_t^{DE}\to {\overline{S}}_t$, and then $C_{GBS}^{DE}\left(t,S_t\right)\to C_{GBS}\left(t,T,S_t,R_{t,T}\right)$. For the Black–Scholes case with jumps, we had $R_{t,T}\to {\sigma }^2\left(T-t\right)$, ${\sigma }_{t,T}\to 0$, ${\mu }_{t,T}\to rT$ and $P\left(t,T\right)\to e^{-r\left(T-t\right)}$. Replacing them in the expression of Theorem 3, we arrive at the call option price in the model of [9]. □

5.3. Generalized Hyperbolic Lévy Motion

The generalized hyperbolic Lévy motion is a Lévy process, $L_t$, such that $\mathcal{L}\left(L_1\right)$, the distribution of $L_1$, has the density of a generalized hyperbolic distribution, ${\rho }_{GH}\left(x\right)$, given by

$$\begin{array}{ccc}\hfill {\rho }_{GH}\left(x;\lambda ,\alpha ,\beta ,\delta ,\mu \right)& =& a\left(\lambda ,\alpha ,\beta ,\delta \right){\left({\delta }^2+{\left(x-\mu \right)}^2\right)}^{\left(\lambda -\frac{1}{2}\right)/2}\hfill \\ & & \times K_{\lambda -\frac{1}{2}}\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)exp\left(\beta \left(x-\mu \right)\right)\hfill \end{array}$$

where

$$a\left(\lambda ,\alpha ,\beta ,\delta \right)=\frac{{\left({\alpha }^2-{\beta }^2\right)}^{\lambda /2}}{\sqrt{2\pi }{\alpha }^{\lambda -\frac{1}{2}}{\delta }^{\lambda }{K}_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\beta }^2}\right)}$$

and $K_{\nu }$ is the modified Bessel function of the third kind with index $\nu$:

$$K_{\nu }\left(z\right)=\frac{1}{2}{\int }_0^{\infty }{y}^{\nu -1}exp\left(-\frac{1}{2}z\left(y+y^{-1}\right)\right)dy$$

The parameter $\alpha >0$ determines the shape, $\beta$, with $0\le \left|\beta \right|<\alpha$, denotes the skewness, $\mu \in \mathbb{R}$ denotes the location, $\delta$ is a scaling parameter, and $\lambda \in \mathbb{R}$ characterizes certain subclasses. The value $\lambda =1$ corresponds to the subclass of hyperbolic distributions, which was the first used in finance. When $\lambda =-\frac{1}{2}$, we obtain the normal inverse Gaussian distribution, which has the gamma distribution as a special case ([11]).

The moment-generating function of the generalized hyperbolic distribution is given in ([11]) as

$$M_{GH}\left(u\right)=e^{\mu u}{\left(\frac{{\alpha }^2-{\beta }^2}{{\alpha }^2-{\left(\beta +u\right)}^2}\right)}^{\lambda /2}\frac{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +u\right)}^2}\right)}{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\beta }^2}\right)}$$

(23)

from which we obtain the characteristic function, ${\phi }_{GH}\left(u\right)$, as ${\phi }_{GH}\left(u\right)=M_{GH}\left(iu\right)$. From this expression, for ${\phi }_{GH}\left(u\right)$, we can easily obtain the characteristic exponent of the generalized hyperbolic Lévy motion as ${\Psi }_{L_1}\left(u\right)=ln\left({\phi }_{GH}\left(u\right)\right)$.

We know from the definition that ${\rho }_{GH}\left(x\right)$ is the density of $L_1$, but $L_t$ does not have the same density and it has to be obtained from its characteristic function, ${\phi }_{t,GH}\left(u\right)\equiv \mathbb{E}\left[e^{iu{L}_t}\right]={\left[{\phi }_{GH}\left(u\right)\right]}^t$ by Fourier inversion (The exception is the case $\lambda =-\frac{1}{2}$, for which we only have to make the substitutions $\delta \to t\delta$ and $\mu \to t\mu$ in the density ${\rho }_{GH}$).

Denoting by ${\rho }_{t,GH}\left(x\right)$ the density function associated with $L_t$, we can define, for each $\theta \in \mathbb{R}$, a new density, ${\rho }_{t,GH}\left(x;\theta \right)$, by

$${\rho }_{t,GH}\left(x;\theta \right)=\frac{e^{\theta x}}{{\int }_{-\infty }^{+\infty }{e}^{\theta y}{\rho }_{t,GH}\left(y\right)dy}$$

(24)

Under the corresponding probability, ${\mathbb{P}}_{\theta }$, the process is again a Lévy process, called the Esscher transform of the original process. We will choose the value ${\theta }^{*}$ of the parameter that makes ${\mathbb{P}}_{{\theta }^{*}}$ an equivalent martingale measure. Taking into account Expression (16), we have

$$\begin{array}{ccc}\hfill S_0& =& {\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left[{\overline{S}}_t\right]\hfill \\ & =& {\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left[S_{0}exp\left\{-\frac{1}{2}{R}_{0,t}+Z_{0,t}^{{\mathbb{P}}_{{\theta }^{*}}}-t{\Psi }_{L_1}(-i)+L_t\right\}\right]\hfill \\ & =& S_0{e}^{-t{\Psi }_{L_1}\left(-i\right)}{\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left[e^{L_t}\right]\hfill \\ & =& S_0{e}^{-tln\left({\phi }_{GH}\left(-i\right)\right)}{\int }_{-\infty }^{+\infty }{e}^x{\rho }_{t,GH}\left(x;{\theta }^{*}\right)dx\hfill \\ & =& S_0{\left({\phi }_{GH}\left(-i\right)\right)}^{-t}\frac{M_{L_t}\left({\theta }^{*}+1\right)}{M_{L_t}\left({\theta }^{*}\right)}\hfill \\ & =& S_0{\left({\phi }_{GH}\left(-i\right)\right)}^{-t}{\left[\frac{M_{GH}\left({\theta }^{*}+1\right)}{M_{GH}\left({\theta }^{*}\right)}\right]}^t\hfill \end{array}$$

where $M_{L_t}\left(u\right)$ is the moment generating function of $L_t$. Taking into account that ${\phi }_{GH}\left(-i\right)=M_{GH}\left(1\right)$, and using Equation (23), we can rewrite the previous expression as

$${\left[\frac{\left({\alpha }^2-{\beta }^2\right)\left({\alpha }^2-{\left(\beta +\left({\theta }^{*}+1\right)\right)}^2\right)}{\left({\alpha }^2-{\left(\beta +1\right)}^2\right)\left({\alpha }^2-{\left(\beta +{\theta }^{*}\right)}^2\right)}\right]}^{-\lambda t/2}{\left[\frac{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +\left({\theta }^{*}+1\right)\right)}^2}\right)}{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +{\theta }^{*}\right)}^2}\right)}\right]}^t=1$$

As this equality must hold for every $t>0$, we arrive at the following equation that determines the value of ${\theta }^{*}$:

$${\left[\frac{\left({\alpha }^2-{\beta }^2\right)\left({\alpha }^2-{\left(\beta +\left({\theta }^{*}+1\right)\right)}^2\right)}{\left({\alpha }^2-{\left(\beta +1\right)}^2\right)\left({\alpha }^2-{\left(\beta +{\theta }^{*}\right)}^2\right)}\right]}^{-\lambda /2}\frac{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +\left({\theta }^{*}+1\right)\right)}^2}\right)}{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +{\theta }^{*}\right)}^2}\right)}=1$$

(25)

The time-T stock price of Equation (16) can be written in this case as

$$S_T={\overline{S}}_t^{GH}exp\left\{{\int }_0^T{r}_{s}ds+Y_{t,T}+L_{t,T}\right\}$$

where we have defined

$$S_t^{GH}=S_{t}exp\left\{-\left(T-t\right)\left(\mu +\frac{\lambda }{2}ln\left(\frac{{\alpha }^2-{\beta }^2}{{\alpha }^2-{\left(\beta +1\right)}^2}\right)+ln\left(\frac{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\left(\beta +1\right)}^2}\right)}{K_{\lambda }\left(\delta \sqrt{{\alpha }^2-{\beta }^2}\right)}\right)\right)\right\}$$

Theorem 4.

The price, $C_{GBS}^{GH}$, of a European call option with pay-off ${\left[S_T-K\right]}_{+}$, is given by

$$\begin{array}{ccc}\hfill C_{GBS}^{GH}\left(t,S_t\right)& =& S_t^{GH}{\int }_{-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}+1\right)\Phi \left(d_1^{GH}\left(x\right)\right)dx\hfill \\ & & -KP\left(t,T\right){\int }_{-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right)\Phi \left(d_2^{GH}\left(x\right)\right)dx\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill d_1^{GH}\left(x\right)& =& \frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\hfill \\ \hfill d_2^{GH}\left(x\right)& =& d_1^{GH}\left(x\right)-\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\hfill \end{array}$$

Proof. 

The value of the call option, $C_{GBS}^{GH}$, is given by

$$\begin{array}{cc}{C}_{GBS}^{GH}\left(t,S_t\right)\hfill & ={\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\left[S_T-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ \hfill & ={\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}-K{\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \end{array}$$

(26)

For the first expectation, we obtain

$$\begin{array}{ccc}& & {\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{S}_T{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& {\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\overline{S}}_t^{GH}exp\left[{\int }_0^T{r}_{s}ds+Z_{t,T}^{{\mathbb{P}}_{{\theta }^{*}}}-\frac{1}{2}{R}_{t,T}+L_{t,T}\right]{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& S_t^{GH}{e}^{-\frac{1}{2}{R}_{t,T}}{\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{exp\left[Z_{t,T}^{{\mathbb{P}}_{{\theta }^{*}}}+L_{t,T}\right]{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& \frac{S_t^{GH}{e}^{-\frac{1}{2}{R}_{t,T}}}{2\pi {\sigma }_{t,T}\sqrt{R_{t,T}}}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{e}^x{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-\frac{{\left(r-{\mu }_{t,T}\right)}^2}{2{\sigma }_{t,T}^2}}{\int }_{z=ln\left(\frac{K}{{\overline{S}}_t^{GH}}\right)-r+\frac{1}{2}{R}_{t,T}-x}^{+\infty }{e}^z{e}^{-\frac{z^2}{2R_{t,T}}}dzdrdx\hfill \\ & =& \frac{S_t^{GH}}{{\sigma }_{t,T}\sqrt{2\pi }}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{e}^x{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-\frac{{\left(r-{\mu }_{t,T}\right)}^2}{2{\sigma }_{t,T}^2}}\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)+r+\frac{1}{2}{R}_{t,T}}{\sqrt{R_{t,T}}}\right)drdx\hfill \\ & =& \frac{S_t^{GH}}{\sqrt{2\pi }}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{e}^x{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-\frac{w^2}{2}}\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)+{\sigma }_{t,T}w+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{R_{t,T}}}\right)dwdx\hfill \\ & =& S_t^{GH}{\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}+1\right)\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\right)dx\hfill \end{array}$$

where, in the third step, we have used that $Z_{t,T}^{{\mathbb{P}}_{{\theta }^{*}}}\sim N\left(0,R_{t,T}\right)$ and ${\int }_0^T{r}_{s}ds\sim N\left({\mu }_{t,T},{\sigma }_{t,T}^2\right)$ under ${\mathbb{P}}_{{\theta }^{*}}$ on ${\mathcal{F}}_t$, and in the last step, we used the equality ${\int }_{-\infty }^{+\infty }\Phi \left(a+bx\right)\varphi \left(x\right)dx=\Phi \left(\frac{a}{\sqrt{1+b^2}}\right)$.

For the second expectation:

$$\begin{array}{ccc}& & {\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& B\left(t\right){\mathbb{E}}^{{\mathbb{P}}_{{\theta }^{*}}}\left\{e^{-{\int }_0^T{r}_{s}ds}{\mathbf{1}}_{S_T>K}\left|{\mathcal{F}}_t\right.\right\}\hfill \\ & =& \frac{B\left(t\right)}{2\pi {\sigma }_{t,T}\sqrt{R_{t,T}}}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-r}{e}^{-\frac{{\left(r-{\mu }_{t,T}\right)}^2}{2{\sigma }_{t,T}^2}}{\int }_{z=ln\left(\frac{K}{{\overline{S}}_t^{GH}}\right)-r+\frac{1}{2}{R}_{t,T}-x}^{+\infty }{e}^{-\frac{z^2}{2R_{t,T}}}dzdrdx\hfill \\ & =& \frac{B\left(t\right)}{2\pi {\sigma }_{t,T}}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-r}{e}^{-\frac{{\left(r-{\mu }_{t,T}\right)}^2}{2{\sigma }_{t,T}^2}}{\int }_{z=\frac{ln\left(\frac{K}{{\overline{S}}_t^{GH}}\right)-r+\frac{1}{2}{R}_{t,T}-x}{\sqrt{R_{t,T}}}}^{+\infty }{e}^{-\frac{w^2}{2}}dwdrdx\hfill \\ & =& \frac{B\left(t\right)}{{\sigma }_{t,T}\sqrt{2\pi }}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-r}{e}^{-\frac{{\left(r-{\mu }_{t,T}\right)}^2}{2{\sigma }_{t,T}^2}}\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)+r-\frac{1}{2}{R}_{t,T}}{\sqrt{R_{t,T}}}\right)drdx\hfill \\ & =& \frac{B\left(t\right)e^{\frac{1}{2}{\sigma }_{t,T}^2-{\mu }_{t,T}}}{\sqrt{2\pi }}\hfill \\ & & \times {\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right){\int }_{r=-\infty }^{+\infty }{e}^{-\frac{v^2}{2}}\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)-{\sigma }_{t,T}v-{\sigma }_{t,T}^2-\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{R_{t,T}}}\right)dv\hfill \\ & =& P\left(t,T\right){\int }_{x=-\infty }^{+\infty }{\rho }_{T-t,GH}\left(x;{\theta }^{*}\right)\Phi \left(\frac{x+ln\left(\frac{{\overline{S}}_t^{GH}}{K}\right)-{\sigma }_{t,T}^2-\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\right)dx\hfill \end{array}$$

Replacing both expectations in Equation (26), we obtain the desired result. □

Remark 1.

In order to apply the previous theorem, we first need to obtain the characteristic function ${\phi }_{t,GH}\left(u\right)={\left[{\phi }_{GH}\left(u\right)\right]}^t$, and then ${\rho }_{t,GH}\left(x\right)$ by Fourier inversion. Secondly, we solve the Equation (25) for ${\theta }^{*}$, and finally, we apply Expression (24) to obtain ${\rho }_{T-t,GH}\left(x;{\theta }^{*}\right)$.

The consistency result in this case is as follows.

Corollary 3.

In the cases of no jumps and only jumps, the expression of Theorem 3 reduces to the call option price in the models of [11,17], respectively.

Proof. 

If there are no jumps, then $S_t^{GH}\to S_t$, and ${\rho }_{T-t,GH}\left(x;{\theta }^{*}\right),{\rho }_{T-t,GH}\left(x;{\theta }^{*}+1\right)\to \delta \left(x\right)$, which lead to $C_{GBS}^{GH}\left(t,S_t\right)\to C_{GBS}\left(t,S_t\right)$. On the other hand, when we consider a pure jump model, as in [11], we have $S_T=S_{t}exp\left(L_{t,T}\right)$. Thus, we should make the replacements ${\sigma }_{t,T},{\mu }_{t,T},R_{t,T}\to 0$, but at the same time, maintaining $P\left(t,T\right)=e^{-r\left(T-t\right)}$. Then, $\Phi \left(d_1^{GH}\left(x\right)\right),\Phi \left(d_2^{GH}\left(x\right)\right)\to {\mathbf{1}}_{x>ln\left(\frac{K}{{\overline{S}}_t}\right)}$, and we recover the call option price of [11]. □

5.4. CGMY

In [10], an infinite activity pure-jump Lévy process called the CGMY process is introduced. This process has the Lévy measure given by

$${\nu }_{CGMY}\left(dx\right)=C\left(\frac{e^{Gx}}{{\left|x\right|}^{1+Y}}{1}_{x<0}+\frac{e^{-Mx}}{x^{1+Y}}{1}_{x>0}\right)dx$$

where $C>0$, $G>0$, $M>0$, and $Y<2$. Its characteristic function is

$${\phi }_{L_t}\left(u\right)=exp\left\{tC\Gamma (-Y)\left({(G+iu)}^Y-G^Y+{(M-iu)}^Y-M^Y\right)\right\}$$

from which we can obtain the characteristic exponent as

$${\Psi }_{L_1}\left(u\right)=C\Gamma (-Y)\left({(G+iu)}^Y-G^Y+{(M-iu)}^Y-M^Y\right)$$

(27)

where $\Gamma (·)$ is the gamma function. The CMGY process has Lévy–Itô decomposition given by

$$L_t=tE\left(X_1\right)+{\int }_0^t{\int }_{\mathbb{R}}x(N(dx,ds)-{\nu }_{CGMY}\left(dx\right)ds)$$

Using (27), we can rewrite Expression (16) for the CGMY process as

$$S_T={\overline{S}}_t^{CGMY}exp\left\{{\int }_0^T{r}_{s}ds+Y_{t,T}+L_{t,T}\right\}$$

with

$$S_t^{CGMY}=S_{t}exp\left\{-\left(T-t\right)C\Gamma (-Y)\left({(G+1)}^Y-G^Y+{(M-1)}^Y-M^Y\right)\right\}$$

The parameters C, G, M, and Y determine the properties of the process ([10]). Specifically, G and M control the rate of exponential decay on the right and the left of the Lévy density, respectively. The parameter $C$ is a measure of the overall level of activity, and in the case of symmetric Lévy measure ($G=M$), it determines the kurtosis of the distribution of the process. The parameter Y characterizes the fine structure of the stochastic process. For example, $Y>-1$ implies a completely monotone Lévy density, allowing large jumps to arrive less frequently than small jumps. When $Y=0$, we recover the variance gamma process, and if $1<Y<2$, the process has infinite activity and infinite variation.

In order to obtain a closed-form expression for the call option price under our model, we will follow [23], which, to the best of our knowledge, presents the unique closed-form solution of the CGMY model for European options. In the following theorem and proof, * denotes the convolution product.

Theorem 5.

The price, $C_{GBS}^{CGMY}$, of a European call option with payoff ${\left[S_T-K\right]}_{+}$, is given by

$$\begin{array}{cc}\hfill & C_{GBS}^{CGMY}(t,S_t)\hfill \\ \hfill & =\frac{k_0}{k_2^{2/Y}}\left\{B\left(t\right){\int }_{-\infty }^{+\infty }{e}^{x-s}\Phi \left(d_1\left(s,x\right)\right)\left(e^{-Gs}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)\ast \left(e^{Ms}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)ds\right.\hfill \\ \hfill & \left.-KP\left(t,T\right){\int }_{-\infty }^{+\infty }\Phi \left(d_2\left(s,x\right)\right)\left(e^{-Gs}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)\ast \left(e^{Ms}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)ds\right\}\hfill \end{array}$$

(28)

where $x=ln{\overline{S}}_t^{CGMY}$, $k_0=exp\left\{-\left(T-t\right)C\Gamma \left(-Y\right)\left(G^Y+M^Y\right)\right\}$, $k_2=\left(T-t\right)C\Gamma \left(-Y\right)$,

$$\begin{array}{ccc}\hfill d_1\left(s,x\right)& =& \frac{x-s-lnK+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\hfill \\ \hfill d_2\left(s,x\right)& =& d_1\left(s,x\right)-\sqrt{{\sigma }_{t,T}^2+R_{t,T}}\hfill \end{array}$$

and $f_{Y,0}\left(x\right)$ is the stable Lévy density given by

$$f_{Y,0}\left(x\right)=\frac{1}{Y}{H}_{2,2}^{1,1}\left[x\left|\begin{array}{cc}\left(1-\frac{1}{Y},\frac{1}{Y}\right)& \left(\frac{1}{2},\frac{1}{2}\right)\\ \left(0,1\right)& \left(\frac{1}{2},\frac{1}{2}\right)\end{array}\right.\right]$$

where $H_{p,q}^{m,n}\left(x\right)$ is the H-function introduced in [24], and given by

$$\begin{array}{ccc}\hfill H_{p,q}^{m,n}\left(x\right)& =& H_{p,q}^{m,n}\left[x\left|\begin{array}{c}\left(a_1,A_1\right),\dots ,\left(a_p,A_p\right)\\ \left(b_1,B_1\right),\dots ,\left(b_q,B_q\right)\end{array}\right.\right]\hfill \\ & =& \frac{1}{2\pi i}{\int }_L\frac{{\Pi }_{j=1}^m\Gamma \left(b_j+B_{j}s\right){\Pi }_{j=1}^n\Gamma \left(1-a_j-A_{j}s\right)}{{\Pi }_{j=m+1}^q\Gamma \left(1-b_j-B_{j}s\right){\Pi }_{j=n+1}^p\Gamma \left(a_j+A_{j}s\right)}{x}^{-s}ds\hfill \end{array}$$

with $A_i,B_j\in {\mathbb{R}}^{+}$, $a_i,b_j\in \mathbb{R}$ or $\mathbb{C}$, $i=1,\dots ,p$; $j=1,\dots ,q$; and L is a suitable contour separating the poles of the gamma functions $\Gamma \left(b_j+B_{j}s\right)$ and $\Gamma \left(1-a_j-A_{j}s\right)$ (a detailed study of the H-function can be found in the references [25,26,27]).

Proof. 

The call option value is given by

$$C_{GBS}^{CGMY}(t,S_t)={\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\left[S_T-K\right]}_{+}\left|{\mathcal{F}}_t\right.\right\}={\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}\Pi \left(V_T\right)\left|{\mathcal{F}}_t\right.\right\}$$

where $V_T\equiv ln{S}_T$ and $\Pi \left(x\right)\equiv {\left[e^x-K\right]}_{+}$.

Making use of the Fourier transform of the payoff

$$\tilde{\Pi }\left(\xi \right)={\int }_{-\infty }^{+\infty }{e}^{-i\xi u}\Pi \left(u\right)du$$

we can write

$$\begin{array}{cc}\hfill & C_{GBS}^{CGMY}(t,S_t)\hfill \\ \hfill & =\frac{1}{2\pi }{\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_t^T{r}_{s}ds}{\int }_{-\infty }^{+\infty }{e}^{i\xi V_T}\tilde{\Pi }\left(\xi \right)d\xi \left|{\mathcal{F}}_t\right.\right\}\hfill \\ \hfill & =\frac{B\left(t\right)}{2\pi }{\mathbb{E}}^{\mathbb{Q}}\left\{e^{-{\int }_0^T{r}_{s}ds}{\int }_{-\infty }^{+\infty }{e}^{i\xi V_T}\tilde{\Pi }\left(\xi \right)d\xi \left|{\mathcal{F}}_t\right.\right\}\hfill \\ \hfill & =\frac{B\left(t\right)}{2\pi }{\int }_{-\infty }^{+\infty }{\mathbb{E}}^{\mathbb{Q}}\left\{e^{i\xi \left(ln{\overline{S}}_t^{CGMY}+Y_{t,T}+L_{t,T}\right)+\left(i\xi -1\right){\int }_0^T{r}_{s}ds}\left|{\mathcal{F}}_t\right.\right\}\tilde{\Pi }\left(\xi \right)d\xi \hfill \end{array}$$

(29)

By the independence assumption, we can factorize the expectation and calculate the expectation of each factor separately. For the first expectation, we have

$${\mathbb{E}}^{\mathbb{Q}}\left\{e^{i\xi ln{\overline{S}}_t^{CGMY}}\left|{\mathcal{F}}_t\right.\right\}=e^{i\xi ln{\overline{S}}_t^{CGMY}}$$

For the second expectation, taking into account that $Y_{t,T}\sim N\left(-\frac{1}{2}{R}_{t,T},R_{t,T}\right)$ under $\mathbb{Q}$ on ${\mathcal{F}}_t$, we obtain

$${\mathbb{E}}^{\mathbb{Q}}\left\{e^{i\xi Y_{t,T}}\left|{\mathcal{F}}_t\right.\right\}=e^{-\frac{i}{2}\xi R_{t,T}-\frac{1}{2}{\xi }^2{R}_{t,T}}$$

The third expectation can be directly obtained as

$${\mathbb{E}}^{\mathbb{Q}}\left\{e^{i\xi L_{t,T}}\left|{\mathcal{F}}_t\right.\right\}=e^{\left(T-t\right){\Psi }_{L_1}\left(\xi \right)}=e^{\left(T-t\right)C\Gamma (-Y)\left({(G+i\xi )}^Y-G^Y+{(M-i\xi )}^Y-M^Y\right)}$$

For the fourth expectation, taking into account that ${\int }_0^T{r}_{s}ds\sim N\left({\mu }_{t,T},{\sigma }_{t,T}^2\right)$ under $\mathbb{Q}$ on ${\mathcal{F}}_t$, we obtain $\left(i\xi -1\right){\int }_0^T{r}_{s}ds\sim N\left(\left(i\xi -1\right){\mu }_{t,T},{\left(i\xi -1\right)}^2{\sigma }_{t,T}^2\right)$ under $\mathbb{Q}$ on ${\mathcal{F}}_t$, and then

$$\begin{array}{ccc}\hfill {\mathbb{E}}^{\mathbb{Q}}\left\{e^{\left(i\xi -1\right){\int }_0^T{r}_{s}ds}\left|{\mathcal{F}}_t\right.\right\}& =& e^{\left(i\xi -1\right){\mu }_{t,T}+\frac{1}{2}{\left(i\xi -1\right)}^2{\sigma }_{t,T}^2}\hfill \\ & =& \overline{P}\left(t,T\right)e^{i\xi \left({\mu }_{t,T}-{\sigma }_{t,T}^2\right)-\frac{1}{2}{\xi }^2{\sigma }_{t,T}^2}\hfill \end{array}$$

where we have used that, for Gaussian models, $\overline{P}\left(t,T\right)=e^{-{\mu }_{t,T}+\frac{1}{2}{\sigma }_{t,T}^2}$ ([17]).

Replacing these expectations in (29), we have

$$\begin{array}{cc}\hfill & C_{GBS}^{CGMY}(t,S_t)\hfill \\ \hfill & =\frac{k_{0}P\left(t,T\right)}{2\pi }\hfill \\ \hfill & \times {\int }_{-\infty }^{+\infty }{e}^{i\xi x}{e}^{{}^{i\xi \left({\mu }_{t,T}-{\sigma }_{t,T}^2-\frac{1}{2}{R}_{t,T}\right)-\frac{1}{2}{\xi }^2\left({\sigma }_{t,T}^2+R_{t,T}\right)}}{e}^{\left(T-t\right)C\Gamma (-Y)\left({(G+i\xi )}^Y+{(M-i\xi )}^Y\right)}\tilde{\Pi }\left(\xi \right)d\xi \hfill \\ \hfill & =k_{0}P\left(t,T\right)\hfill \\ \hfill & \times {\mathcal{F}}^{-1}\left[e^{i\xi k_1}\tilde{\Pi }\left(\xi \right)\right]\ast {\mathcal{F}}^{-1}\left[e^{-\frac{1}{2}{\xi }^2{k}_3}\right]\ast {\mathcal{F}}^{-1}\left[e^{k_2{\left(G+i\xi \right)}^Y}\right]\ast {\mathcal{F}}^{-1}\left[e^{k_2{\left(M-i\xi \right)}^Y}\right]\hfill \end{array}$$

(30)

where we have applied the convolution theorem for inverse Fourier transforms, and we defined $k_1={\mu }_{t,T}-{\sigma }_{t,T}^2-\frac{1}{2}{R}_{t,T}$, and $k_3={\sigma }_{t,T}^2+R_{t,T}$.

For the two first inverse transforms, we obtain

$${\mathcal{F}}^{-1}\left[e^{i\xi k_1}\tilde{\Pi }\left(\xi \right)\right]=\Pi \left(x+k_1\right)$$

and

$${\mathcal{F}}^{-1}\left[e^{-\frac{1}{2}{\xi }^2{k}_3}\right]=\frac{1}{\sqrt{2\pi k_3}}{e}^{-\frac{x^2}{2k_3}}$$

where, for the first one, we have applied the shift theorem for inverse Fourier transforms. For the two remaining transforms, we apply Lemma 3.1 of [23] to obtain

$${\mathcal{F}}^{-1}\left[e^{k_2{\left(G+i\xi \right)}^Y}\right]=e^{-Gx}\frac{1}{k_2^{1/Y}}{f}_{Y,0}\left(\frac{\left|x\right|}{k_2^{1/Y}}\right)$$

(31)

and

$${\mathcal{F}}^{-1}\left[e^{k_2{\left(M-i\xi \right)}^Y}\right]=e^{Mx}\frac{1}{k_2^{1/Y}}{f}_{Y,0}\left(\frac{\left|x\right|}{k_2^{1/Y}}\right)$$

(32)

With these inverse transforms, we can obtain explicitly the first convolution in (30) as

$$\begin{array}{ccc}& & {\mathcal{F}}^{-1}\left[e^{i\xi k_1}\tilde{\Pi }\left(\xi \right)\right]\ast {\mathcal{F}}^{-1}\left[e^{-\frac{1}{2}{\xi }^2{k}_3}\right]\hfill \\ & =& \Pi \left(x+k_1\right)\ast \frac{1}{\sqrt{2\pi k_3}}{e}^{-\frac{x^2}{2k_3}}\hfill \\ & =& \frac{1}{\sqrt{2\pi k_3}}{\int }_{-\infty }^{+\infty }{\left[e^{x-s+k_1}-K\right]}_{+}{e}^{-\frac{s^2}{2k_3}}ds\hfill \\ & =& e^{x+k_1+\frac{1}{2}{k}_3}\Phi \left(\frac{x+k_1-lnK+k_3}{\sqrt{k_3}}\right)-K\Phi \left(\frac{x+k_1-lnK}{\sqrt{k_3}}\right)\hfill \\ & =& e^x{\overline{P}}^{-1}\left(t,T\right)\Phi \left(\frac{x-lnK+\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\right)\hfill \\ & & -K\Phi \left(\frac{x-lnK-\frac{1}{2}{R}_{t,T}+{\mu }_{t,T}-{\sigma }_{t,T}^2}{\sqrt{{\sigma }_{t,T}^2+R_{t,T}}}\right)\hfill \end{array}$$

(33)

Substituting (31)–(33) in (30), we obtain the desired result. □

As in the previous examples of jump processes, we have a consistency result for stochastic strings plus CGMY jumps.

Corollary 4.

In the cases of no jumps and jumps without diffusion part, the expression of Theorem 5 reduces to the call option price in the models of [17,23], respectively.

Proof. 

The no jumps case is given by $C\to 0$, and then $k_0\to 1$ and

$$\begin{array}{ccc}& & \left(e^{-Gs}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)\ast \left(e^{Ms}{f}_{Y,0}\left(\frac{\left|s\right|}{k_2^{1/Y}}\right)\right)\hfill \\ & =& {\mathcal{F}}^{-1}\left[e^{k_2{\left(G+i\xi \right)}^Y}\right]\ast {\mathcal{F}}^{-1}\left[e^{k_2{\left(M-i\xi \right)}^Y}\right]\to \delta \left(s\right)\ast \delta \left(s\right)=\delta \left(s\right)\hfill \end{array}$$

Replacing these results in (28), and taking into account that $x=ln{\overline{S}}_t^{CGMY}\to ln{\overline{S}}_t$, we arrive at the call option price in [17]. On the other hand, in the case of CGMY jumps and no diffusion part, we have $R_{t,T}\to 0$, ${\mu }_{t,T}\to rT$, ${\sigma }_{t,T}\to 0$, $B_t\to e^{rt}$, $P\left(t,T\right)\to e^{-r\left(T-t\right)}$; and then, the two normal cumulative distribution functions in (28) tend to ${\mathbf{1}}_{s<ln\left(\frac{{\overline{S}}_t^{CGMY}}{K}\right)+rT}$. Introducing all these changes in (28), we recover the call option price expression in [23]. □

6. Alternative Models

In this section, we present, for completeness, a brief summary of each of the benchmark models that will be used for comparison with our GBS plus jumps model. These models are the Merton jump-diffusion model and the Heston stochastic volatility model with jumps.

6.1. Merton Jump-Diffusion Model

Starting from the price dynamics of [1]

$$\frac{d{S}_t}{S_t}=\mu dt+\sigma d{W}_t$$

where the parameters $\mu$ and $\sigma >0$ are constants, and $W_t$ is a standard Brownian motion under the physical probability measure $\mathbb{P}$, ref. [8] adds lognormally distributed jumps, so that

$$\frac{d{S}_t}{S_{t-}}=\mu dt+\sigma d{W}_t+d\left(\sum _{i=1}^{N_t}(e^{J_i}-1)\right),$$

(34)

where $N_t$ is a Poisson process with parameter $\lambda >0$, and $J\sim N(m,{\delta }^2)$. Under $\mathbb{Q}$, (34) becomes

$$\frac{d{S}_t}{S_{t-}}=rdt+\sigma d{W}_t^{*}+d\left(\sum _{i=1}^{N_t}(e^{J_i}-1)\right)-\lambda \left(e^{m+{\delta }^2/2}-1\right)dt,$$

(35)

where $r\in \mathbb{R}$ is the risk-free interest rate and $W_t^{*}$ is a Brownian motion under an equivalent martingale measure.

The solution to (35) is given by

$$S_t=S_{0}exp\left\{(r-{\sigma }^2/2-\lambda (e^{m+{\delta }^2/2}-1))t+\sigma W_t^{*}+\sum _{i=1}^{N_t}{J}_i\right\}$$

and the price of a European call option with exercise price K and time-to-maturity $\tau$ is given by ([22])

$$C_{MJD}(t,\tau ,S_t,K)=e^{-\lambda \tau }\sum _{n=0}^{\infty }\frac{{\left(\lambda \tau \right)}^n}{n!}{C}_{BS}(t,\tau ,S_n,{\sigma }_n),$$

where $C_{BS}(t,\tau ,S,\sigma )$ is the Black–Scholes price of a call option with the same strike, underlying price S, and volatility $\sigma$, $S_n=S_{t}exp\left[nm+\frac{n{\delta }^2}{2}-\lambda \tau (e^{m+{\delta }^2/2}-1)\right]$, and ${\sigma }_n=\sqrt{{\sigma }^2+n{\delta }^2/\tau }$.

The addition of jumps makes the model much more flexible and capable of capturing the skew and kurtosis observed in the empirical distributions of asset returns. The calibration of the Merton jump-diffusion model entails the estimation of four parameters ($\sigma$, $\lambda$, m, $\delta$).

6.2. Heston Model with Jumps

For the Heston stochastic volatility model with jumps, we have taken the stochastic volatility, stochastic interest rates and random jumps (SVSR-J) model of [14]. This model is an extension of the model in [7] with the addition of the CIR stochastic interest rates, and the compound Poisson jumps with the lognormal distribution of jump sizes. This model provides a semi-closed form expression for the price of a European call option, and generalizes most of the models with closed or semi-closed formulas for the option price.

The Heston model with jumps is defined by the following equations:

$$\begin{array}{ccc}\hfill \frac{d{S}_t}{S_t}& =& \left(r_t-\lambda {\mu }_J\right)dt+\sqrt{v_t}d{\tilde{W}}_{1t}+J_{t}d{q}_t\hfill \\ \hfill d{v}_t& =& {\varkappa }_v\left({\eta }_v-v_t\right)dt+{\sigma }_v\sqrt{v_t}d{\tilde{W}}_{2t}\hfill \\ \hfill d{r}_t& =& {\varkappa }_r\left({\eta }_r-r_t\right)dt+{\sigma }_r\sqrt{r_t}d{\tilde{W}}_{rt}\hfill \end{array}$$

where:

• $v_t$ is the instantaneous variance of returns (conditional of no jumps occurring);
• ${\tilde{W}}_{1t}$ and ${\tilde{W}}_{2t}$ are standard $\mathbb{Q}$-Brownian motions with $d{\left[{\tilde{W}}_{1·},{\tilde{W}}_{2·}\right]}_t=\rho dt$;
• $J_t$ is the percentage jump size (conditional on jump occurring), with $ln\left(1+J_t\right)\sim N\left(ln\left(1+{\mu }_J\right)-\frac{1}{2}{\sigma }_J^2,{\sigma }_J^2\right)$;
• $q_t$ is a Poisson jump counter with intensity $\lambda$;
• ${\varkappa }_v$, ${\eta }_v$ and ${\sigma }_v$ are, respectively, the mean-reversion speed, the long-run mean, and the volatility coefficient of the variance process;
• ${\varkappa }_r$, ${\eta }_r$ and ${\sigma }_r$ are, respectively, the mean-reversion speed, the long-run mean, and the volatility coefficient of the risk-free interest rate process;
• ${\tilde{W}}_{rt}$ is a $\mathbb{Q}$-Brownian motion uncorrelated with any other process.

The time-t price of a European call option written on S, with strike price K, and time-to-maturity $\tau$, $C_{HJ}(t,\tau ,S_t,K)$, is given by

$$C_{HJ}(t,\tau ,S_t,K)=S_t{\Pi }_1\left(t,\tau ;S_t,r_t,v_t\right)-KP\left(t,t+\tau \right){\Pi }_2\left(t,\tau ;S_t,r_t,v_t\right)$$

where the conditional probabilities ${\Pi }_1$ and ${\Pi }_2$ can be obtained from

$${\Pi }_j\left(t,\tau ;S_t,r_t,v_t\right)=\frac{1}{2}+\frac{1}{\pi }{\int }_0^{\infty }\mathcal{R}\left[\frac{e^{-i\varphi lnK}{f}_j\left(t,\tau ;S_t,r_t,v_t;\varphi \right)}{i\varphi }\right]d\varphi$$

and $f_j$ are the characteristic functions, whose explicit expression can be found in [14].

Setting $\lambda ={\varkappa }_r={\eta }_r={\sigma }_r=0$, we recover the [7] model and, therefore, the Heston model with jumps generalizes the original Heston model to stochastic interest rates and jumps. Moreover, when ${\varkappa }_v={\eta }_v={\sigma }_v=0$, the Heston model with jumps reduces to the Merton jump-diffusion model. Nevertheless, we will keep on using the Merton model in order to assess the effect of introducing stochastic volatility.

In order to calibrate the Heston model with jumps (excluding the CIR model), we need to estimate eight parameters ($\lambda ,{\varkappa }_v,{\eta }_v,{\sigma }_v,v_0,\rho ,{\mu }_J,{\sigma }_J$).

7. Empirical Tests

Given the characteristic functions for underlying returns processes, a variety of Fourier transform techniques, such as those developed in [28] or [29], can be used to price European options. Under the assumption that ${\int }_0^T{r}_{s}ds$, $Y_{t,T}$, and $L_t$ are mutually independent, we can calibrate a short-rate model for interest rates separately from the model for the returns of the underlying security, and then combine them when calculating option values. Specifically, we calibrate a CIR short-rate model using zero coupon bond (ZCB) yield curve data for each pricing date. Then, for each expiration, a constant yield-to-maturity (YTM) is implied from the calibrated ZCB price. Options are priced with the constant but expiration-dependent short-rate equal to the YTM. This is consistent with [14] in their SVJ case, but not with SVSR-J, which estimates the CIR parameters at the same time as the stock parameters. So our empirical study does not consider SR models, but rather models with constant but maturity-dependent interest rates estimated from an SR model. There are two main reasons for which we take this approach. First, estimating CIR and underlying stock parameters simultaneously adds three additional parameters to each optimization. Since our focus is on combining strings with jumps to model the dynamics of the underlying security, we prefer to keep the models as parsimonious as possible in order to focus on the marginal improvements achieved by changing the parameterizations of the underlying dynamics. Second, Ref. [14] finds little to no tangible improvement from including SR, either in their study of in-sample fit or in their study of hedging performance.

The characteristic function for the return process of the Gaussian string model under the $\mathbb{Q}$–measure, $Y_{t,T}$, is given by

$${\phi }_{\tau }^{GS}\left(u\right)=exp\left({\omega }_{1}iu\tau -\frac{{\sigma }_{TS}{\left(\tau \right)}^2}{2}\tau u^2\right)$$

(36)

with $\tau =T-t$, ${\sigma }_{TS}\left(\tau \right)$ as defined in (7), and ${\omega }_1=-{\sigma }_{TS}{\left(\tau \right)}^2/2$.

If ${\phi }_{\tau }^J$ is the characteristic function for $L_{\tau }$ under $\mathbb{Q}$, then the characteristic function for the string plus jumps process is given by the product,

$${\phi }_{\tau }\left(u\right)={\phi }_{\tau }^{GS}\left(u\right){\phi }_{\tau }^L\left(u\right).$$

(37)

The price of a European call option can then be expressed as a Black–Scholes style formula as derived in [29]

$$C(S_t,K,T)=S_t{\Pi }_1-KP(t,T){\Pi }_2$$

where

$${\Pi }_1=\frac{1}{2}+\frac{P(t,T)}{\pi }{\int }_0^{\infty }\mathcal{R}\left[\frac{e^{iulog(S_t/K)}{\phi }_T(u-i)}{iu}\right]du,$$

and

$${\Pi }_2=\frac{1}{2}+\frac{P(t,T)}{\pi }{\int }_0^{\infty }\mathcal{R}\left[\frac{e^{iulog(S_t/K)}{\phi }_T\left(u\right)}{iu}\right]du$$

7.1. The Carr and Madan Methods for Fourier Transform Pricing

In [28], Carr and Madan derive two Fourier transform pricing methods to be applied depending on the moneyness of the call option to be priced. For the in-the-money case, they begin with the following formula for the call option price

$$C\left(k\right)={\int }_k^{\infty }{e}^{-rT}(e^s-e^k)f\left(s\right)ds$$

where $k=ln\left(K\right)$ and $f_T\left(s\right)$ is the density function for $ln\left(S_T\right)$. However, $C\left(k\right)$ is not square-integrable, so they work with a dampened call price function

$$c\left(k\right)=e^{\alpha k}C\left(k\right),$$

where $\alpha >0$. Now, the Fourier transform $\widehat{c}\left(u\right)$ of $c\left(k\right)$ exists and it can be shown that

$$\widehat{c}\left(u\right)={\int }_{-\infty }^{\infty }{e}^{iuk}c\left(k\right)dk=\frac{e^{-rT}\phi (u-(\alpha +1)i)}{{\alpha }^2+\alpha -u^2+iu(2\alpha +1)}$$

For $\alpha =0$, $\widehat{c}\left(u\right)$ is undefined at $u=0$. Thus, $\alpha$ must be nonzero. The call price can then be found by Fourier inversion

$$C\left(k\right)=\frac{e^{-\alpha k}}{2\pi }{\int }_{-\infty }^{\infty }{e}^{iuk}\widehat{c}\left(u\right)du=\frac{e^{-\alpha k}}{\pi }{\int }_0^{\infty }\mathcal{R}\left[e^{-iuk}\widehat{c}\left(u\right)\right]du$$

For pricing out-of-the money options, [28] proceeded as follows. If $k=ln\left(K\right)$, define

$$z_T\left(k\right)=e^{-rT}{\int }_{-\infty }^{\infty }\left((e^k-e^s)1_{s<k,k<0}+(e^s-e^k)1_{s>k,k>0}\right)f_T\left(s\right)ds$$

which is the price of a call option if $k>ln\left(S_0\right)$ or the price of a put option if $k<ln\left(S_0\right)$. Now, let ${\zeta }_T$ be the Fourier transform of $z_T\left(k\right)$,

$$\begin{array}{cc}\hfill {\zeta }_T\left(v\right)& ={\int }_{-\infty }^{\infty }{e}^{ivk}{z}_T\left(k\right)dk\hfill \\ \hfill & ={\int }_{-\infty }^{\infty }{e}^{ivk}{e}^{-rT}{\int }_{-\infty }^{\infty }\left((e^k-e^s)1_{s<k,k<0}+(e^s-e^k)1_{s>k,k>0}\right)f_T\left(s\right)dkds\hfill \\ \hfill & ={\int }_{-\infty }^0{e}^{ivk}{e}^{-rT}{\int }_{-\infty }^k(e^k-e^s)q_T\left(s\right)dkds+{\int }_0^{\infty }{e}^{ivk}{e}^{-rT}{\int }_k^{\infty }(e^s-e^k)f_T\left(s\right)dkds\hfill \\ \hfill & ={\int }_{-\infty }^0{e}^{-rT}{q}_T\left(s\right){\int }_s^{\infty }(e^{(1+iv)k}-e^{s+ivk})dsdk+{\int }_0^{\infty }{e}^{-rT}{q}_T\left(s\right){\int }_0^s(e^{s+ivk}-e^{(1+iv)k})dsdk\hfill \\ \hfill & =e^{-rT}\left(\frac{1}{1+iv}-\frac{e^{rt}}{iv}-\frac{{\phi }_T(v-i)}{v^2-iv}\right).\hfill \end{array}$$

The pricing formula follows after applying the Fourier inversion to

$${\int }_{-\infty }^{\infty }{e}^{ivk}sinh\left(\alpha k\right)z_T\left(k\right)dk=\frac{{\zeta }_T(v-i\alpha )-{\zeta }_T(v+i\alpha )}{2},$$

to obtain

$$z_T\left(k\right)=\frac{1}{2\pi sinh\left(\alpha k\right)}{\int }_{-\infty }^{\infty }{e}^{-ivk}\left(\frac{{\zeta }_T(v-i\alpha )-{\zeta }_T(v+i\alpha )}{2}\right)dv.$$

In both moneyness cases, prices can be calculated easily using the fast Fourier transform and standard numerical integration procedures. Given the efficiency and accuracy of the Fourier methods developed in [28], we will use them to calculate option prices in the empirical part of this study.

7.2. Calibrating Models to Option Market Data

The calibration of option pricing models to market data constitutes a critical empirical test of their validity. We conduct a calibration exercise using actual option data to assess the value of including jumps in the stochastic string option pricing model. In this exercise, we utilize the Option Metrics volatility surface for the S&P 500 index to evaluate the merit of incorporating jumps into the stochastic string option pricing framework. We select three distinct dates— 3 February 2021, representing moderate market volatility; 30 August 2017, indicative of a low volatility climate; and 5 November 2008, a period of heightened volatility—to anchor our analysis. The corresponding 30-day 50 delta Black–Scholes-implied volatilities for these dates were 19.1%, 9.1%, and 45.5%, respectively. The closing levels of the S&P 500 index for these dates were, respectively, $3830.17$, $2457.59$, and $952.77$.

For each date, we extract option prices for eight expirations—spanning from one to twelve months (30, 60, 91, 122, 152, 182, 273, and 365 days)—and for each expiration, we collect call option prices (and strike prices) for calls with deltas ranging in moneyness from 10% to 90% in increments of 5%. Thus, for each trading date, we have 17 call option prices for each expiration, for a total of 136 call option prices per trading date. Our objective is, for each date, to estimate a single constant parameter set that optimally aligns model-generated option prices with prices observed in the market, minimizing the root mean square error (RMSE) across all expirations. The calibration is extended to include a Cox–Ingersoll–Ross (CIR) model for the short-rate, fitted to the zero-coupon U.S. Treasury yield curve on each selected date.

The RMSE, a goodness-of-fit measure for a calibrated model, is computed as:

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{(C_i-{\tilde{C}}_i)}^2}$$

(38)

where $C_i$ denotes the model-predicted call option prices, ${\tilde{C}}_i$ denotes the market-observed prices, and n denotes the aggregate count of observations. The minimization task is formulated as an optimization problem:

$$\underset{\theta \in \Theta }{min}RMSE({\left\{C_{\theta }(S_i,K_i,T_i)\right\}}_{i=1}^n,{\left\{\tilde{C}(S_i,K_i,T_i)\right\}}_{i=1}^n)$$

(39)

wherein:

• $\theta$ symbolizes the vector of model parameters, with $\Theta$ representing the permissible parameter space;
• $C_{\theta }(S_i,K_i,T_i)$, for $i=1\dots n$, is the price of a call option as determined by the model for given parameters $\theta$, underlying asset price $S_i$, strike price $K_i$, and time to expiration $T_i$;
• $\tilde{C}(S_i,K_i,T_i)$, for $i=1\dots n$, is the corresponding market price of the option.

For each market scenario, we calibrated a suite of option pricing models, considering various configurations of underlying price dynamics and jump components. The computational implementation of the Carr and Madan Fourier transform methods using the Python programming language is detailed in [30]. The optimization of the calibration function is carried out via the Nelder–Mead simplex algorithm, a gradient-free method well suited for application to nonlinear optimization problems. The algorithm employs a simplex or an N-dimensional polytope, to probe the objective function at its vertices and iteratively apply transformations—reflection, expansion, contraction, or shrinkage—based on the function’s evaluated results. This iterative process continues until some termination criterion is met, such as convergence to a specified tolerance level or the exhaustion of a maximum number allotted iterations.

The following list summarizes the various model configurations we will consider:

• Underlying models without jumps:

  –

  Black–Scholes (BS) ([1]) (one parameter);

  –

  Heston stochastic volatility (HSV) ([7]) (five parameters);

  –

  Generalized Black–Scholes (GBS) ([17]) (three parameters).

• Only jumps:

  –

  CGMY ([10]) (four parameters);

  –

  Generalized hyperbolic Lévy motion (GHLM) ([11]) (five parameters);

• Underlying models plus jumps:

  –

  Brownian motion plus CGMY (BM+CGMY) (five parameters);

  –

  Brownian motion plus GHLM (BM+GHLM) (six parameters);

  –

  BS plus lognormal compound Poisson (BS+LN) ([8]) (four parameters);

  –

  BS plus double-exponential compound Poisson (BS+DE) ([9] ) (five parameters);

  –

  HSV plus lognormal compound Poisson (HSV+LN) ([14]) (eight parameters);

  –

  GBS plus lognormal compound Poisson (GBS+LN) (five and six parameters);

  –

  GBS plus double-exponential compound Poisson (GBS+DE) (six and seven parameters);

  –

  GBS plus CGMY (GBS+CGMY) (six and seven parameters);

  –

  GBS plus GHLM (GBS+GHLM) (seven and eight parameters).

We report the RMSEs in

Table 1. Comparative analysis of in-sample model performance for different volatility conditions: This table reports the root mean square errors (RMSEs) of various option pricing models calibrated against market data from three distinct volatility environments. The RMSEs serve as a quantitative metric for the precision of each model’s fit to the observed option prices. The calibration exercise spans models without jump processes, models exclusively featuring jumps, and hybrid models combining the underlying price dynamics with jumps. Lower RMSE values denote a more accurate representation of market prices by the model. The volatility regimes are defined by the dates: low volatility on 3 February 2021, moderate volatility on 30 August 2017, and high volatility on 5 November 2008.

Underlying Models without Jumps	Parameters	Low Vol	Moderate Vol	High Vol
Black–Scholes (BS)	1	36.040	12.619	10.388
Heston stochastic volatility (HSV)	5	13.031	4.088	3.901
Generalized Black–Scholes (GBS(3))	3	35.794	11.879	8.117
Only Jumps	Parameters	Low Vol	Moderate Vol	High Vol
CGMY	4	9.563	1.284	8.339
Generalized hyperbolic Lévy motion (GHLM)	5	11.815	3.521	9.437
Underlying Models Plus Jumps	Parameters	Low Vol	Moderate Vol	High Vol
BS plus lognormal compound Poisson (BS+LN)	4	10.530	4.858	8.204
BS plus double-exponential compound Poisson (BS+DE)	5	8.491	1.138	8.606
Brownian motion plus CGMY (BM+CGMY)	5	9.410	1.165	8.280
Brownian motion plus GHLM (BM+GHLM)	6	10.268	2.907	8.867
GBS plus lognormal compound Poisson (GBS(3)+LN)	6	4.632	1.465	2.175
GBS plus double-exponential compound Poisson (GBS+DE)	7	3.700	1.083	3.577
GBS plus CGMY (GBS(3)+CGMY)	7	3.630	1.050	3.403
HSV plus lognormal compound Poisson (HSV+LN)	8	6.128	2.116	2.293
GBS plus GHLM (GBS(3)+GHLM)	8	6.512	2.892	4.828
GBS plus CGMY (GBS(4)+CGMY)	8	3.536	1.049	2.006

and optimal parameter values for the calibrated models are reported in

Table 2. Parameter estimation for calibrated models under high-volatility conditions: This table presents the estimated parameter values for a range of option pricing models calibrated to the S&P 500 index options data from 5 November 2008—a date characterized by high market volatility.

Underlying Models without Jumps	Parameters
Black–Scholes	1	$\sigma$
		0.356
Stochastic string model (GBS(3))	3	$\sigma$	$a_2$	$a_3$
		0.459	−0.222	0.117
Heston stochastic volatility (HSV)	5	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$
		1.216	0.073	0.420	−0.938	0.176
Only Jumps	Parameters
CGMY	4	C	G	M	Y
		71.234	12.546	5624.007	−0.574
GHLM	5	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		9.990	−2.819	1.428	−1.182	0.320
Underlying Models Plus Jumps	Parameters
Brownian motion plus CGMY (BM+CGMY)	5	$\sigma$	C	G	M	Y
		0.01	71.234	12.546	5624.007	−0.574
Brownian motion plus GHLM (BM+GHLM)	6	$\sigma$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		0.216	10.829	−3.785	0.169	4.258	0.272
BS plus lognormal compound Poisson (BS+LN)	4	$\sigma$	$\lambda$	$\mu$	$\delta$
		0.038	6.929	−0.130	0.062
BS plus double-exponential compound Poisson (BS+DE)	5	$\sigma$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.368	0.010	0.226	2.096	0.001
HSV plus lognormal compound Poisson (HSV+LN)	8	$\lambda$	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$	${\mu }_J$	${\sigma }_J$
		6.487	0.009	0.116	−1.000	0.189	0.522	−0.444	$1.591\times {10}^{-8}$
GBS plus lognormal compound Poisson (GBS(3)+LN)	6	$\sigma$	$a_2$	$a_3$	$\lambda$	$\mu$	$\delta$
		0.367	−0.264	0.142	0.947	-0.364	0.005
GBS plus double-exponential compound Poisson (GBS+DE)	7	$\sigma$	$a_2$	$a_3$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.486	0.001	−0.261	0.013	0.001	2.100	0.121
GBS plus CGMY (GBS(3)+CGMY)	7	$\sigma$	$a_2$	$a_3$	C	G	M	Y
		0.301	−0.112	0.023	477.300	11.023	78.750	−2.324
GBS plus GHLM (GBS(3)+GHLM)	8	$\sigma$	$a_2$	$a_3$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		−0.178	0.075	0.320	9.998	−4.818	0.028	3.352	0.357
GBS plus CGMY (GBS(4)+CGMY)	8	$\sigma$	$a_2$	$a_3$	$a_4$	C	G	M	Y
		0.362	−0.330	0.274	−0.076	1589.635	12.836	6343.843	−2.756

,

Table 3. Parameter estimation for calibrated models under low-volatility conditions: This table presents the estimated parameter values for a range of option pricing models calibrated to the S&P 500 index options data from 30 August 2017—a date characterized by low-market volatility.

Underlying Models without Jumps	Parameters
Black–Scholes	1	$\sigma$
		0.120
Stochastic string model (GBS(3))	3	$\sigma$	$a_2$	$a_3$
		0.086	0.020	−0.01
Heston stochastic volatility (HSV)	5	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$
		1.840	0.028	0.323	−1.000	0.013
Only Jumps	Parameters
CGMY	4	C	G	M	Y
		0.092	2.909	53.105	0.841
GHLM	5	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		8.571	4.610	0.064	0.072	0.069
Underlying Models Plus Jumps	Parameters
Brownian motion plus CGMY (BM+CGMY)	5	$\sigma$	C	G	M	Y
		0.046	0.470	5.282	1008.237	0.207
Brownian motion plus GHLM (BM+GHLM)	6	$\sigma$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		0.046	12.714	−7.301	0.002	0.639	0.076
BS plus lognormal compound Poisson (BS+LN)	4	$\sigma$	$\lambda$	$\mu$	$\delta$
		0.085	0.089	−0.498	0.001
BS plus double-exponential compound Poisson (BS+DE)	5	$\sigma$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.057	0.091	0.000	33.941	9.963
HSV plus lognormal compound Poisson (HSV+LN)	8	$\lambda$	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$	${\mu }_J$	${\sigma }_J$
		32.743	0.005	0.000	0.089	0.004	0.261	−0.252	0.001
GBS plus lognormal compound Poisson (GBS(3)+LN)	6	$\sigma$	$a_2$	$a_3$	$\lambda$	$\mu$	$\delta$
		0.060	0.001	−0.001	0.533	−0.153	0.106
GBS plus double-exponential compound Poisson (GBS+DE)	7	$\sigma$	$a_2$	$a_3$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.055	0.001	−0.001	0.091	0.000	26.920	9.938
GBS plus CGMY (GBS(3)+CGMY)	7	$\sigma$	$a_2$	$a_3$	C	G	M	Y
		0.048	$4.787\times {10}^{-4}$	−0.001	0.719	5.900	6199.174	0.044
GBS plus GHLM (GBS(3)+GHLM)	8	$\sigma$	$a_2$	$a_3$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		0.051	0.000	0.001	10.168	−5.410	−0.000	0.504	0.062
GBS plus CGMY (GBS(4)+CGMY)	8	$\sigma$	$a_2$	$a_3$	$a_4$	C	G	M	Y
		0.049	$4.953\times {10}^{-4}$	−0.001	$-2.000\times {10}^{-5}$	0.813	6.089	639.277	−0.004

and

Table 4. Parameter estimation for calibrated models under moderate-volatility conditions: This table presents the estimated parameter values for a range of option pricing models calibrated to the S&P 500 index options data from 3 February 2021—a date characterized by moderate market volatility.

Underlying Models without Jumps	Parameters
Black–Scholes	1	$\sigma$
		0.222
Stochastic string model (GBS(3))	3	$\sigma$	$a_2$	$a_3$
		0.086	0.020	−0.01
Heston stochastic volatility (HSV)	5	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$
		0.701	0.160	0.473	−1.000	0.042
Only Jumps	Parameters
CGMY	4	C	G	M	Y
		0.052	0.238	122.079	1.277
GHLM	5	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		8.571	4.610	0.064	0.072	0.069
Underlying Models Plus Jumps	Parameters
Brownian motion plus CGMY (BM+CGMY)	5	$\sigma$	C	G	M	Y
		0.050	0.043	0.083	201.355	1.291
Brownian motion plus GHLM (BM+GHLM)	6	$\sigma$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		0.086	4.653	−4.286	0.187	−0.987	0.151
BS plus lognormal compound Poisson (BS+LN)	4	$\sigma$	$\lambda$	$\mu$	$\delta$
		0.127	0.384	−0.380	0.202
BS plus double-exponential compound Poisson (BS+DE)	5	$\sigma$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.134	0.279	0.347	71.226	3.885
HSV plus lognormal compound Poisson (HSV+LN)	8	$\lambda$	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$	${\mu }_J$	${\sigma }_J$
		3.122	0.005	0.016	0.145	0.028	0.280	−0.494	0.006
GBS plus lognormal compound Poisson (GBS(3)+LN)	6	$\sigma$	$a_2$	$a_3$	$\lambda$	$\mu$	$\delta$
		0.150	−0.020	0.005	0.521	−0.302	0.240
GBS plus double-exponential compound Poisson (GBS+DE)	7	$\sigma$	$a_2$	$a_3$	$\lambda$	p	${\theta }_1$	${\theta }_2$
		0.147	−0.016	0.003	0.234	0.221	789.440	4.195
GBS plus CGMY (GBS(3)+CGMY)	7	$\sigma$	$a_2$	$a_3$	C	G	M	Y
		0.122	−0.015	0.000	0.227	1.537	74.000	0.550
GBS plus GHLM (GBS(3)+GHLM)	8	$\sigma$	$a_2$	$a_3$	$\lambda$	$\alpha$	$\beta$	$\delta$	$\mu$
		0.136	−0.016	0.004	7.491	−5.848	0.043	0.162	0.149
GBS plus CGMY (GBS(4)+CGMY)	8	$\sigma$	$a_2$	$a_3$	$a_4$	C	G	M	Y
		0.123	0.017	$1.128e^{-4}$	0.002	0.273	1.653	15.638	0.383

. In Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5, we plot option prices from calibrated models compared to market prices for a selection of the models.

The calibration exercise yields several key insights, not only about the comparative performance of the option pricing models under different market conditions, but also about the fine structure of the stochastic processes underlying equity price dynamics. Based on the RMSEs reported in Table 1, the addition of jump processes to models without jumps such as Black–Scholes (BS), Heston’s Stochastic Volatility (HSV), and generalized Black–Scholes (GBS) consistently improves the fit across all three volatility regimes. This fact supports the hypothesis of the presence of jumps in price dynamics.

Moreover, the GBS(3) model augmented with any jump process (lognormal, double-exponential, CGMY, GHLM) consistently outperforms its pure jump counterparts, as well as its continuous sample path alternatives. These results indicate that GBS plus jump models can indeed provide robust fits for the market prices of options across moneyness and expiry, simultaneously.

Specifically, the GBS(3) model with CGMY jumps provides superior fits of option skews compared to the alternative models considered in two out of the three volatility environments. The Heston model with jumps of [14] outperforms in the high volatility case, although it has one more parameter than GBS(3)+CGMY. Adding an additional stochastic string parameter to the GBS model with CGMY jumps, so that it has a total of eight parameters, results in dramatic reductions in RMSEs across all three volatility environments, with improvements exceeding 50% in the moderate- and high-volatility cases. In particular, the GBS(4)+CGMY model outperforms HSV+LN in each case. Close scrutiny of the plots in Figure 6 reveals that, in the high volatility case, the GBS(4)+CGMY is better able to fit the skews for out-of-the-money (OTM) calls with short expiries, while the HSV+LN appears to provide a better fit for deep OTM calls with 273 and 365 days until expiration. Overall, we can conclude that, on an equal parameter count basis, the GBS model with CGMY jumps provides the best in-sample fits, across all three volatility environments, compared to any of the alternative models reported in Table 2, Table 3, Table 4 and

Table 5. Model parameter estimates on 5 November 2008 for two 8-parameter models, HSV+LN and GBS(4)+CGMY.

Underlying Model	Parameters	$\mathit{RMSE}$
HSV plus lognormal compound Poisson (HSV+LN)	8		$\lambda$	${\varkappa }_v$	${\eta }_v$	${\sigma }_v$	$\rho$	$v_0$	${\mu }_J$	${\sigma }_J$
		2.293	6.487	0.009	0.116	−1.000	0.189	0.522	−0.444	$1.591\times {10}^{-8}$
GBS plus CGMY (GBS(4)+CGMY)	8		$\sigma$	$a_2$	$a_3$	$a_4$	C	G	M	Y
		2.006	0.366	−0.338	0.287	−0.082	1648.945	12.840	5503.566	−2.787

.

Parameter estimation from the calibrated models shows that the jump component in the GBS(4)+CGMY model, which has the best overall performance, exhibits a variance gamma process-like behavior $Y=0$ under moderate volatility conditions. The Lévy density is skewed in the three scenarios, with the left tail heavier than the right tail $G<M$, indicating a predominance of left-tail risk across all volatility scenarios. The parameter Y delineates distinct properties related to the fine structure of the stochastic process ([10]). In the high-volatility scenario ($Y<-1$), the process is not completely monotone and has finite activity, while for moderate volatility ($-1<Y<0$), the process is completely monotone with finite activity. In the low-volatility case ($0<Y<1$) the process has a completely monotone Lévy density with infinite activity and finite variation. The parameter C reflects the intensity of jumps within the CGMY process, with a markedly higher level of activity in the high volatility regime ($C=1589.635$) compared to the moderate and low volatility scenarios ($C=0.813$ and $C=0.273$, respectively), suggesting the existence of a volatility-dependent jump intensity in market dynamics.

8. Conclusions

In this paper, we present a new model for pricing European options that extends the Black–Scholes framework to include stochastic interest rates, stochastic string dynamics, and jump processes. Our model synthesizes the generalized Black–Scholes model proposed by [17] with an array of pure jump Lévy processes. Thus, our model is also a generalization of exponential Lévy models. We derive closed or semi-closed form pricing formulas for European call options in GBS plus jump models with four different types of jumps: compound Poisson processes with lognormal and double-exponential jumps ([8,9]), the CGMY process of [10], and the generalized hyperbolic Lévy motion of [11]. We also show that these expressions are reduced to those obtained in the literature for the cases of Black–Scholes dynamics or in the absence of jumps.

In order to assess the performance of the GBS model with jumps, we conduct in-sample empirical tests against several alternative models consisting of fitting the skews of options on S&P 500 under three different volatility scenarios. The results of the test can be summarized as follows: (a) including jumps improves the performance of all the models; (b) the GBS plus jumps model outperforms alternative models with the same jump distributions; (c) the GBS model with CGMY jumps (seven parameters) provides superior in-sample fits of option skews than the alternative models considered in two out of the three volatility regimes. The Heston model with jumps of [14] (eight parameters) is better in the high volatility regime, although it has one more parameter to estimate. Adding an additional stochastic string parameter to the GBS model with CGMY jumps, so that it has a total of eight parameters outperforms all alternative models in each of the three volatility regimes.

The findings of this paper not only firmly support the presence of random jump discontinuities in equity index dynamics, but also challenge the necessity of stochastic volatility in option pricing models. Our empirical results suggest that models combining continuous dynamics featuring a deterministic volatility function conditioned on option expirations, in addition to random jump processes, can yield equal or superior performance compared to alternatives that specify volatility as a stochastic process. The blend of simplicity and precision offered by the GBS model with jumps is well suited for practical application under diverse market conditions.

Building on the models developed in this paper, there are a number of interesting directions for future research. For instance, the utility of GBS plus jumps models for dynamic hedging applications could be investigated. Using the Black–Scholes hedge ratio combined with GARCH volatility estimates, Ref. [31] documented that delta-neutral portfolios consisting of long positions in index call options delta-hedged by short-index positions, on average, experience negative returns. This finding provides evidence of a negative volatility risk premium. Since a calibrated GBS model contains an estimated term-structure of volatility, conducting a dynamic hedging study using a GBS-derived hedge ratio could yield further insights on the pricing of volatility risk.

Additionally, the pricing of American and Asian options, as well as several types of exotic options, with a GBS plus jumps model could be investigated with simulation methods. Such a study would necessitate the development of Monte Carlo methods combining stochastic strings with pure jump Lévy processes. Ref. [18] simulated forward interest rates in a stochastic string model by discretizing stochastic string shocks as increments of correlated Brownian motion. They performed the simulation with three different parametric correlations. Efficient Monte Carlo methods for simulating the sample paths of the underlying price as exponential Lévy processes are developed in [32] and applied to price barrier options and fixed and floating strike lookback options. With the simulated paths of the underlying price, the least-squares Monte Carlo (LSM) technique of [33] can be applied to price American options for underlying dynamics given by a stochastic string plus jumps. In fact, the examples of derivative pricing using the LSM technique presented in [33] include both an American put when the underlying price follows a jump-diffusion, as well as the deferred American swaption in a 20-factor string model of the term structure. Moreover, as studied in [34], a modified LSM technique can be applied to price a variety of exotic path-dependent options.