LOG-NORMAL STOCHASTIC VOLATILITY MODEL WITH QUADRATIC DRIFT

1. Introduction

Empirical studies strongly support evidence that the volatility of price returns is itself a stochastic process (see, for example, Shephard (2005) for a comprehensive survey). It is accepted that the celebrated model by Black & Scholes (1973) and Merton (1973), which assumes constant volatility, cannot explain implied volatility surfaces observed in option markets, which are inhomogeneous in strike and maturity dimensions. In contrast, stochastic volatility (SV) models are able to fit to implied volatility surfaces and their dynamics.

1.1. Evidence of log-normality of implied and realized volatilities

Applications of log-normal SV models, which are based either on the log-normal dynamics for the instantaneous volatility or on the Gaussian dynamics for the logarithm of the volatility, are widespread in practice. Predominantly, SABR model by Hagan et al. (2002) is widely used among practitioners for fitting implied volatility surfaces. However, for modeling the dynamics of volatility surfaces, we need to account for the mean-reversion of the volatility process like in conventional SV models. The mean-reversion ensures that the distribution of the volatility has a stationary long-run distribution and volatility paths simulated in Monte Carlo (MC) do not explode. In another widely used specification, the so-called Exp-OU specification, the SV process is defined for the logarithm of the volatility.^b

In an excellent study Christoffersen et al. (2010) examine the empirical performance of Heston (see Heston (1993)), log-normal and 3/2 SV models using three sources of market data: the VIX index, the implied volatility for options of the S&P500 index, and the realized volatility of returns on the S&P500 index. They found that, for all three sources of data, the log-normal SV model outperforms its alternatives. Tegner & Poulsen (2018) report similar finding using realized volatilities.

1.2. Analytical tractability

Log-normal SV models cannot be handled by semi-analytic Fourier methods available for affine SV models. The analytical intractability may have lead to a slow adaptation of log-normal SV models in spite of their empirical support. On the one hand, the log-normal SV model with the linear drift does not belong to the family of affine diffusions because the so-called admissibility conditions are not satisfied for the variance term, see Duffie et al. (2003), Filipović & Mayerhofer (2009). On the other hand, the log-normal SV model with the quadratic drift does not belong to a general family of polynomial diffusions, see Filipović & Larsson (2016), Cuchiero et al. (2012). Carr & Willems (2019) introduce measure change to relate the model to a class of polynomial diffusions, which comes at cost of introducing an additional state variable. Authors then approximate the payoff function using an orthogonal polynomial expansion, see Ackerer & Filipović (2020), for option pricing applications. A related approach for computing moments under SV models using expansions is developed by Alòs et al. (2020).

We contribute to these studies by developing a closed-form approximate and very accurate solution for the log-normal SV model which enables analytical intractability of this model. We develop an affine expansion to derive closed-form solution for the moment generating function (MGF) of the log-price process, volatility and its quadratic variance (QV). We show that our solution to the MGF produces valid probability density functions of the state variables, which are matched with MC simulations.

1.3. Properties of log-normal SV model

In a related research, Lewis (2018) introduces a log-normal SV model with a quadratic drift without constant term. Lewis (2018) suggests to apply the model only when return–volatility correlation is negative because, otherwise, the model leads to the loss of martingality. Carr & Willems (2019) introduce the log-normal SV model with the quadratic drift, similar to our model. In our paper we augment the log-normal SV model with the linear drift term introduced in Karasinski & Sepp (2012) to include the quadratic drift term. We contribute to studies of Lewis (2018) and Carr & Willems (2019) by defining the martingale conditions under both MMA and inverse martingale measures. We also formulate the Feynman–Kac theorem linking the model dynamics with the quadratic drift to the partial differential equation (PDE) approach, because the classic Feynman–Kac does not apply to non-Lipschitz coefficients. We then introduce the MGF as the solution to the model PDE and develop an analytic transform-based approach for the valuation of vanilla and inverse options.

1.4. Assets with positive return–volatility correlation

Carr & Wu (2017) propose the following three economic factors behind the negative return–volatility correlation, which most typically observed in the dynamics of stock prices and stock indices: aggregate financial leverage, systematic risk, positive auto-correlation of negative shocks. However, in many markets we actually observe price dynamics with positive return–volatility correlation either on a permanent basis (for VIX index, short and leveraged short exchange-traded funds (ETFs)) or on a transitory regime-based basis (for some commodities, foreign currencies, and cryptocurrencies). In fact, exp-ou and log-normal SV models with the linear drift lead to the “loss of martingality” of the price process when the return–volatility correlation is positive (see, for example, Lewis (2000, 2018) and Lions & Musiela (2007)). We contribute by establishing conditions on model parameters of log-normal SV model which produce true martingale dynamics.

1.5. Inverse martingale measure

For inverse options, both option premium and final payoff are quoted in the units of the underlying asset. In particular, for options on cryptocurrencies, the advantage of inverse payoffs is that all option-related transactions can be handled using units of underlying cryptocurrencies, such as Bitcoin or Ether, without using fiat currencies (see Lucic (2021) and Lucic & Sepp (2023) for a connection between the valuation of the inverse cryptocurrency options and FX options). In case of cryptocurrencies, cash settled options would require using a stablecoin as a tradable equivalent of United States Dollar (USD) or any other fiat currency. Given inherent and frequent depegging risk of stablecoins, Alexander et al. (2023) argue for wider adoptation of inverse and quanto options. We also suggest that similar payoff can be advantageous for some commodity markets (such as gold) where option holders may opt for physical settlement of positive option premium through the delivery of the underlying commodity.

We contribute to the literature by incorporating the inverse measure in our valuation framework using log-normal SV model. We formulate and solve the joint valuation problem of vanilla and inverse options using Fourier transform methods.

1.6. Organization of the paper

In Sec. 2, we provide the framework for valuation under money-market account (MMA) and inverse measures. In Sec. 3, we introduce the log-normal beta SV model with quadratic drift and establish important properties of the volatility process. We further characterize martingale property for the drift-adjusted price process and its inverse under MMA and inverse measures, respectively. We also study moments of the volatility process by presenting an approximate solution of the finite-dimensional linear system of ordinary differential equations (ODEs). Finally, we present drift-implicit Euler scheme for MC simulation that preserves the positivity of the volatility process and study its convergence properties. In Sec. 4, we present first- and second-order affine expansion for the MGF of the log-price process, the QV and the volatility. In Sec. 5, we apply affine expansion to derive closed-form solutions for option valuation problem under MMA and inverse measures. In Sec. 6, we apply our method for model calibration to market data and demonstrate the accuracy of both the model specification and our proposed solution methods. We conclude in Sec. 7. Technical proofs are provided in Appendix A.

2. Framework and Definitions

2.1. Vanilla and inverse payoffs

The payoff of (cash-settled) vanilla options is settled in the units of cash. In contrast, the payoff of an inverse option is settled in the units of the underlying asset using the spot price observed at the option maturity. We denote the spot price at time T by $S_T$ and option strike price by K, with both quoted in United States Dollars (USD). We consider vanilla call and put payoffs settled in USD cash at maturity time T and denote them by

Payoffs of inverse call and put options are converted to units of the underlying asset at maturity T as follows :

2.2. Valuation under equivalent MMA and inverse measures

We consider a continuous-time market with a fixed horizon date $T^{\ast }>0$ and uncertainty modeled on probability space $(\mathrm{\Omega },𝔽,\mathbb{P})$ equipped with filtration $𝔽={\{{{ℱ}}_t\}}_{0\le t\le T^{\ast }}$. We assume that $𝔽$ is right-continuous and satisfies usual conditions.

2.2.1. Spot markets

Assumption 2.1. Risk-free rate $r(t)$ is deterministic with the value of one unit of the MMA $M(T)$ given by $M(T)=e^{\bar{r}(0,T)}$, $\bar{r}(t,T)={\int }_t^{T}r(s)ds$.

We assume that the market is complete and satisfies the so-called No Free Lunch with Vanishing Risk (NFLVR) and No Dominance (ND) conditions (see Jarrow & Larsson (2012)). By choosing M as a numéraire, we consider an equivalent martingale measure $\mathbb{Q}$, induced by M. Accordingly, the time-t value of an option, denoted by $U(t,S)$, with payoff function $u(S_T)$ at time T, equals

where expectation $𝔼$ is taken under the MMA martingale measure $\mathbb{Q}$.

Next, by choosing spot price S as a numéraire, we consider an inverse martingale measure $\tilde{\mathbb{Q}}$, induced by S. The option value function, denoted by $Ũ(t,S)$, with the payoff paid in units of S, equals

where expectation $\widetilde{𝔼}$ is taken under the inverse martingale measure $\tilde{\mathbb{Q}}$.

Theorem 2.1. We assume that market is complete and satisfies NFLVR and ND conditions. We further assume that both MMA measure $\mathbb{Q}$ and inverse measure $\tilde{\mathbb{Q}}$ are equivalent martingale measures. Then the values of options under MMA measure in Eq. (2.3) and under inverse measure in Eq. (2.4) are unique and satisfy

Proof. As shown in Theorem 3.2 of Jarrow & Larsson (2012) and Theorem 2.17 of Herdegen & Schweizer (2018), the market satisfies NFLVR and ND conditions if and only if $\mathbb{Q}$ is a true martingale measure. Hence, asset prices, discounted by numéraire M are true $\mathbb{Q}$-martingales. Consequently, due to the uniqueness of price, the equality in Eq. (2.5) holds. To justify the equality of both sides, we apply the existence and uniqueness of equivalent inverse (true) martingale measure $\tilde{\mathbb{Q}}$ for the numéraire S, see Theorem 4.4 and Corollary 3.11 in Herdegen & Schweizer (2018). □

2.2.2. Futures markets

Since most option markets for commodities and cryptocurrencies are based on futures, we also adopt our methodology for futures options. We consider a set of futures contracts settled at times ${\{T_k\}}_{k=1,\dots ,K}$ with prices ${\{F_t^{T_k}\}}_{k=1,\dots ,K}$.

Assumption 2.2. The convenience yield $q(t)$ of futures prices is deterministic with the integrated convenience yield given by

Definition 2.1. Measure ${\tilde{\mathbb{Q}}}^T$ induced by price $F_t^T$ with futures maturity at T, $T\in \{T_k\}$, chosen as a numéraire is called the inverse T-futures measure.

Corollary 2.1. We assume that market is complete and satisfies NFLVR and ND conditions. We further assume that both MMA measure $\mathbb{Q}$ and inverse futures measure ${\tilde{\mathbb{Q}}}^T$ are equivalent martingale measures. Then option values under MMA and inverse measures are unique and satisfy

For markets with both vanilla and inverse options being tradable, such as cryptocurrency markets, any SV model applied for joint valuation of options across several exchanges and contract specifications must satisfy the conditions (2.5) and (2.7). For the proposed log-normal SV dynamics in Eq. (3.12), we must impose certain conditions on model parameters, which we address in Theorem 3.7.

2.2.3. Put–call parity

We will discuss an important topic of the put–call parity for inverse options under a generic SV model with SV driver ${\sigma }_t$. We establish that the put–call parity is valid if volatility ${\sigma }_t$ is not explosive.

Lemma 2.1. We assume that volatility process ${\sigma }_t$ is nonexplosive under $\mathbb{Q}$. Then put–call parity for vanilla call option, $C(t,S,K),$ and put option, $P(t,S,K),$ defined by

holds irrespective whether the price process is a true martingale under $\mathbb{Q}$.

We assume that volatility process ${\sigma }_t$ is nonexplosive under $\tilde{\mathbb{Q}}$. Put–call parity for inverse call option, $\tilde{C}(t,S,K),$ and inverse put option, $\tilde{P}(t,S,K),$ defined by

holds irrespective whether (drift-adjusted) inverse price process is a true martingale under $\tilde{\mathbb{Q}}$.

Proof. We only provide a sketch of the proof and refer to Theorems 9.2 and 9.3 in Lewis (2000) and to Sin (1998) for details. As vanilla call option payoff is unbounded, the price of the call option is only local martingale under $\mathbb{Q}$. Taking into account “correction” term (so-called martingale defect) reflecting the difference between true martingale and strict local martingale properties, yields Eq. (2.8). Relationship in Eq. (2.9) follows from the same argument, but applied to $\tilde{\mathbb{Q}}$. □

In Theorem 3.2, we show that volatility process ${\sigma }_t$ does not explode under the proposed log-normal SV dynamics in Eq. (3.12).

3. Log-Normal Beta SV Model

3.1. Dynamics under $\mathbb{P}$ measure

We start with statistical $\mathbb{P}$ measure and we consider the log-normal beta SV model introduced by Karasinski & Sepp (2012). We augment the dynamics with the quadratic drift in the volatility process following Lewis (2018) and Carr & Willems (2019). We first introduce the dynamics of price process $S_t$ and volatility process ${\sigma }_t$ under statistical measure $\mathbb{P}$ as follows :

where ${\mu }_t$ is statistical drift, $W^{(0)}$, $W^{(1)}$ are uncorrelated Brownian motions under $\mathbb{P}$, ${\kappa }_1>0$ and ${\kappa }_2\ge 0$ are linear and quadratic mean-reversion rates, respectively, $𝜃>0$ is the mean of the volatility, $\beta \in ℝ$ is the volatility beta which measures the sensitivity of the volatility to changes in the spot price, and $𝜀>0$ is the volatility of residual volatility.

Empirically, Bakshi et al. (2006) find that the presence of nonlinear terms, including the quadratic term, in the volatility drift is significant when using econometric estimation of the dynamics of the VIX index.

We assume that statistical measure $\mathbb{P}$ and MMA measure $\mathbb{Q}$ are related by

where ${ℰ}(\cdot )$ denotes Itô exponential, and ${\lambda }_0(t)$ and ${\lambda }_1(t)$ are equity and volatility risk-premias, respectively.

Assuming that process $\zeta (t)$ is a true martingale, the following processes

define independent Brownian motions under $\mathbb{Q}$ by Girsanov theorem (we establish necessary and sufficient conditions for equivalence of measures in Theorem 3.1). Thus, the dynamics of price process in Eq. (3.1) under MMA measure $\mathbb{Q}$ becomes As discounted price process must be a martingale under $\mathbb{Q}$, we need to require that

As a result, the dynamics of volatility process under MMA measure $\mathbb{Q}$ become

where ${Ŵ}^{(0)}$, ${Ŵ}^{(1)}$ are independent Brownian motions under $\mathbb{Q}$.

As customary,^c we assume that the market price of risk is proportional to the volatility as follows :

where ${\bar{\lambda }}_0$ and ${\bar{\lambda }}_1$ are constants. As a result, we rewrite Eq. (3.5) as follows :

To retain the functional form of the volatility process, we rescale model parameters $\widehat{𝜃}$, ${\widehat{\kappa }}_1$, ${\widehat{\kappa }}_2$ under $\mathbb{Q}$ as follows :

Analytic inversion of Eq. (3.8) yields risk-transformed model parameters under $\mathbb{Q}$ :

Thus, under martingale measure $\mathbb{Q}$, the dynamics of price and volatility processes become

Theorem 3.1. We assume that ${\kappa }_1\ge 0,$ $𝜃>0$. Then measures $\mathbb{P}$ and $\mathbb{Q}$ are equivalent, iff

Proof. By analogy to the proof of Theorem 3.6. □

As a result, we obtain that the functional form of log-normal SV model with the quadratic drift is invariant under the change of measure from $\mathbb{P}$ to $\mathbb{Q}$. We note that the log-normal volatility model with linear drift does not share such invariance property due to a quadratic term arising under $\mathbb{Q}$. The model invariance under the change of measure is an important feature shared by the log-normal SV model with quadratic drift. In contrast, Sepp & Rakhmonov (2023) show that most conventional SV models do not satisfy the invariance of $\mathbb{P}$-model under the change of measure to $\mathbb{Q}$. In the following, we will consider model dynamics under $\mathbb{Q}$ and drop the hat notation in Eq. (3.10).

3.2. Dynamics under MMA measure

We consider price dynamics with the log-normal beta SV model in Eq. (3.1) under the MMA measure $\mathbb{Q}$. The joint dynamics are specified for the spot price $S_t$, the instantaneous volatility ${\sigma }_t$ and the QV $I_t$ under the MMA measure $\mathbb{Q}$ as follows :

where $r(t)$ is the deterministic risk-free rate; $W_t^{(0)}$ and $W_t^{(1)}$ are uncorrelated Brownian motions, and the model parameters are the same as defined in Eq. (3.1). Using volatility beta $\beta$ and volatility-of-volatility $𝜀$, we introduce the total instantaneous variance of the volatility process ${𝜗}^2$ as follows :

The mean-reversion of the volatility process is specified using the mean level of the volatility $𝜃$, $𝜃>0$, and the mean-reversion speed which is the linear function of the volatility ${\kappa }_1+{\kappa }_2{\sigma }_t$, ${\kappa }_1\ge 0,{\kappa }_2\ge 0$. As a result, the quadratic drift of the volatility process is given by

We introduce the zero-drift stochastic driver $Z_t$, $Z_t=e^{-{\int }_0^{t}r(t^{\prime })d{t}^{\prime }}{S}_t$ and its inverse $R_t$, $R_t=Z_t^{-1}$, with the following dynamics :

Next, we introduce the log of the zero-drift price process $X_t=ln{Z}_t$ :

Using $X_T$ under Assumption 2.1, we represent the spot price $S_T$ by

Similarly, under Assumption 2.2, using Eq. (3.17), we represent the futures settled at time T in $\mathbb{Q}$ as where the integrated convenience rate $\bar{q}(0,T)$ is defined in Eq. (2.6). Accordingly, for modeling purposes, we focus on the zero-drift price dynamics $Z_t$ and its log-price process $X_t$ defined in Eq. (3.16).

3.2.1. Properties of log-normal SV process

We establish important properties related to the regularity of the volatility process. We consider the volatility process in model dynamics (3.12) represented as follows :

where $W_t^{(\ast )}$ is a standard Brownian motion and $𝜗$ is the total volatility of volatility.

We notice that the process ${\sigma }_t$ is not a polynomial diffusion, see Filipović & Larsson (2016), Lemma 2.2, because the drift coefficient is a second-order polynomial in ${\sigma }_t$. Since the process in Eq. (3.19) has super-linearly growing drift, standard results for the uniqueness of strong solution are no longer applicable, see Sec. 5 of Karatzas & Shreve (1991) for details. We will establish the existence of the unique strong solution in Theorem 3.3. We first provide a boundary classification for the volatility process to show that log-normal volatility process ${\sigma }_t$ in dynamics (3.19) cannot reach 0 or explode to $+\infty$, which are important properties of the volatility dynamics. Our proof is based on existence of solutions for general stochastic differential equations (SDEs) with so-called monotonic coefficients, so that it is potentially applicable to wider class of diffusion processes.

Assumption 3.1. We assume that model parameters ${\kappa }_1$, ${\kappa }_2$, $𝜃$, $𝜗$ satisfy

Theorem 3.2 (Regularity). Under Assumption 3.1, the boundary points $\{0,+\infty \}$ are unattainable for process ${\sigma }_t$ in SDE (3.19).

Proof. Our proof is based on Feller boundary classification for one-dimensional diffusion processes with technical details given in Sec. A.1. □

Theorem 3.3 (Existence and uniqueness of solution). Under Assumption 3.1, SDE in Eq. (3.19) has the unique strong solution.

Proof. We will refer to following theorem that guarantees the existence and uniqueness of general SDEs with monotonic coefficients.

Theorem 3.4 (Gyöngy & Krylov (1980)). Let $a_t(x):[0,T]\times ℝ\to ℝ$, $b_t(x):[0,T]\times ℝ\to ℝ$ be functions, continuous in x, satisfying following conditions:

(1)	${\int }_0^T{sup}_{|x|\le R}(|a_t(x)|+|b_t(x){|}^2)dt<\infty$ for any $T\ge 0,$ $R\ge 0,$
(2)	for $x,y,z\in ℝ$ such that $|x|\le R,|y|\le R,|z|\le R$ following conditions hold: (a) (monotonicity condition) : $$2(x-y)(a_t(x)-a_t(y))+{(b_t(x)-b_t(y))}^2\le K_t(R){(x-y)}^2,$$(3.21) (b) (growth condition): $2z{a}_t(z)+{(b_t(z))}^2\le K_t(1){(1+z)}^2,$ where function $K_t(R):[0,T]\times {ℝ}_{+}\to {ℝ}_{+}$ satisfies ${\int }_0^t{K}_s(R)ds<\infty$ for any $T\ge 0,$ $R\ge 0$.

Then there exists unique solution of the SDE

We show that SDE (3.19), where ${\kappa }_1,{\kappa }_2,𝜃,𝜗$ satisfy (3.20), does indeed satisfy conditions 1–3. We emphasize that the positivity of the process demonstrated in Theorem 3.2 is important. Since in our case $a_t(x)=({\kappa }_1+{\kappa }_{2}x)(𝜃-x)$, $b_t(x)=𝜗x$, we get

when $|x|\le R,|y|\le R$. As for any $R>0$ we have ${\int }_0^t[{𝜗}^2-2({\kappa }_1-{\kappa }_2𝜃)+2{\kappa }_{2}R]ds<\infty$, we conclude that 2a is satisfied. For the growth condition, we have

Using that zero is unattainable boundary and simple inequalities $z\le 1+z^2$, $z^3<0$ valid for $z>0$, we can estimate the right hand side of (3.23) as follows :

We conclude that the growth condition 2b is also satisfied. □

Proposition 3.1 (Bounded moments). Assume that ${\kappa }_1𝜃>0,$ ${\kappa }_2>0$ and remaining parameters ${\kappa }_2,$ $𝜃,$ $𝜗$ satisfy conditions (3.20). Let $|p|\ge 1$ and $T>0$. Then we have

Proof. First we prove a slightly weaker result. □

Proposition 3.2. Assume that ${\kappa }_1𝜃>0,$ ${\kappa }_2>0$ and ${\kappa }_2,$ $𝜃,$ $𝜗$ satisfy conditions (3.20). Then, for every $|p|\ge 1:$ ${sup}_{t\in [0,T]}𝔼{\sigma }_t^p<\infty$ for any $T>0$.

Proof. Let us first consider the case $p\ge 1$. We use same approach as in Lions & Musiela (2007). Applying Itô’s lemma to ${\sigma }_t^p$, we have

The function $U(\tau )=𝔼({\sigma }_t^p)$ then solves the problem Let us find supersolution $Ū$ to (3.24) using ansatz $Ū(\sigma )={\sigma }^{p}exp(a\sigma +b\tau )$, $a>0$, $b>0$. Inserting it into (3.24), we see that its left-hand side equals Since the term in the brackets is a second-order polynomial in $\sigma$, $Ū$ is a supersolution as soon as ${\kappa }_2>0$.

The case $p\le -1$ is treated similarly: we denote $q:=-p\ge 1$ and introduce ${\tilde{\sigma }}_t:={\sigma }_t^{-1}$ which is well-defined due to Theorem 3.2. We find that $Ũ(\tau ):=𝔼({\tilde{\sigma }}_t^q)$ solves the problem

Inserting ansatz into (3.25), we see that its left-hand side equals where we omitted the terms not depending on b and $\tilde{\sigma }$ for brevity. Since the term in square brackets is a second-order polynomial in $\tilde{\sigma }$, $Ũ$ is a supersolution if ${\kappa }_1𝜃>0$. □

Similarly to the proof of Proposition 3.2, we apply Itô’s formula to ${\sigma }_t^p$

We obtain

Using Burkhölder–Davis–Gundy inequality and Proposition 3.2, we establish

3.3. Dynamics under inverse measure

First we define the model dynamics (3.12) under the inverse measure $\tilde{\mathbb{Q}}$. As a by-product, we show that in SDE (3.26) the quadratic term is important to preserve the functional specification of volatility drift under $\tilde{\mathbb{Q}}$.

Theorem 3.5. The dynamics of the volatility ${\sigma }_t$ in (3.12) under $\tilde{\mathbb{Q}}$ are given by

MMA measure $\mathbb{Q}$ and inverse measure $\tilde{\mathbb{Q}}$ are related by the density process

where zero-drift stochastic driver $Z_t$ is defined in Eq. (3.15).

Proof. We consider $R_t=Z_t^{-1}$. Applying Itô’s lemma, we see that $R_t$ satisfies

Using Girsanov’s theorem, we switch to the measure $\tilde{\mathbb{Q}}$ under which the processes ${\tilde{W}}_t^{(0)}$ and ${\tilde{W}}_t^{(1)}$ defined by ${\tilde{W}}_t^{(0)}=W_t^{(0)}-{\int }_0^t{\sigma }_{s}ds$, ${\tilde{W}}_t^{(1)}=W_t^{(1)}$ are uncorrelated Brownian motions. Then the dynamics of volatility process ${\sigma }_t$ under $\tilde{\mathbb{Q}}$ become

and Eq. (3.26) follows. We note that $\tilde{\mathbb{Q}}$ and $\mathbb{Q}$ are related by the density ${\mathrm{\Lambda }}_t$ as follows :  □

We highlight that the idea of using change of measure to study the martingale property is widely used in the literature, see Ruf (2015) and references therein. We now derive sufficient conditions for the equivalence of measures $\mathbb{Q}\sim \tilde{\mathbb{Q}}$ which is required for the application of Girsanov theorem in the proof of Theorem 3.5.

Theorem 3.6. Under Assumption 3.1, measures $\mathbb{Q}$ and $\tilde{\mathbb{Q}}$ are equivalent if and only if ${\kappa }_2\ge \beta$.

Proof. By Theorem 3.5, the density ${𝔼}_t[d\tilde{\mathbb{Q}}/d\mathbb{Q}]=Z_t/Z_0$. Therefore, the martingale property for $Z_t$ under $\mathbb{Q}$ guarantees equivalence $\mathbb{Q}\sim \tilde{\mathbb{Q}}$. Using Theorem 3.7(1) concludes the proof. □

Theorem 3.7 (Martingality of Log-normal SV model with quadratic drift). Under Assumption 3.1, the price processes $Z_t$ and $R_t$ in Eq. (3.15) for model dynamics (3.12) have following properties:

(1)	The process $Z_t$ is a martingale under the MMA measure $\mathbb{Q}$ iff ${\kappa }_2\ge \beta$.
(2)	The process $R_t$ is a martingale under the inverse measure $\tilde{\mathbb{Q}}$ iff ${\kappa }_2\ge 2\beta$.

Proof. (1) Theorem 9.2 in Lewis (2000) or Sin (1998) shows that $Z_t$ is $\mathbb{Q}$-martingale if and only if the process

has unique strong solution under $\tilde{\mathbb{Q}}$. Using Theorem 3.3, we immediately obtain that the uniqueness holds iff ${\kappa }_2\ge \beta$.

(2) According to Theorem 3.5, the dynamics of $\{{\sigma }_t,R_t\}$ under $\tilde{\mathbb{Q}}$ become

Using Theorems 2.4(i), 2.4(ii) of Lions & Musiela (2007), $R_t$ is a martingale if

which holds if ${\kappa }_2\ge 2\beta$. $R_t$ is not a martingale if for some smooth, positive, increasing function $\phi$ we have where ${\int }_{𝜀}^{+\infty }\phi {(\xi )}^{-1}d\xi <\infty$, $𝜀>0$. Choosing $\phi (\xi )={\xi }^2$, (3.31) holds if $\kappa <2\beta$. □

Corollary 3.1 (Martingality of Log-normal SV model with linear drift). For dynamics (3.12) with ${\kappa }_2=0$ we obtain the following:

(1)	The process $Z_t$ is martingale under $\mathbb{Q}$ iff $\beta \le 0$.
(2)	The process $R_t$ is martingale under $\tilde{\mathbb{Q}}$ iff $\beta \le 0$.

Proof. Follows by setting ${\kappa }_2=0$ in Theorem 3.7. □

Importantly we conclude that the log-normal SV model with linear drift leads to arbitrages when the return–volatility correlation is positive.

3.4. Joint dynamics of state variables under $\mathbb{Q}$ and $\tilde{\mathbb{Q}}$

We introduce the mean-adjusted volatility process $Y_t$

on the domain $(-𝜃,\infty )$. The marginal dynamics of $Y_t$ using SDE (3.19) become where $W_t^{(\ast )}$ is a standard Brownian motion.

Corollary 3.2 (Dynamics under the MMA measure ℚ). The joint dynamics of log-price $X_t,$ the mean-adjusted volatility process $Y_t,$ and the QV $I_t$ under the MMA measure $\mathbb{Q}$ are driven by

Itô’s lemma, applied to Eq. (3.33), yields marginal dynamics of Y under $\tilde{\mathbb{Q}}$ :

where ${\tilde{W}}_t^{(\ast )}$ is a standard Brownian motion under $\tilde{\mathbb{Q}}$.

Corollary 3.3 (Dynamics under the inverse measure ℚ̃). The joint dynamics of log-price $X_t$, the mean-adjusted volatility process $Y_t,$ and the QV $I_t$ under the inverse measure $\tilde{\mathbb{Q}}$ are driven by

3.5. Steady-state distribution of volatility process

It is important to establish that the steady-state density (SSD) of the volatility process ${\sigma }_t$ in SDE (3.19) exists. The SSD, denoted by $G(\sigma )$, solves the following ODE :

The solution is Generalized Inverse Gaussian distribution (Jørgensen (1982)) : where $b>0$, $q\ge 0$, $\eta \in ℝ$ and $K_{\eta }$ is a modified Bessel function of the second kind. For ${\kappa }_2=0$ we use $c(0)=q^{-\eta }\mathrm{\Gamma }(-\eta )$, where $\mathrm{\Gamma }(x)$ is the gamma function, and the steady state PDF is the inverse gamma function. The rth moment $m(r)$ associated to the density in Eq. (3.38) is given by We see that when ${\kappa }_2=0$, the volatility moments exist only to the order r less than $r<1+2{\kappa }_1/{𝜗}^2$.

In Fig. 1(a), we show the SSD for the three sets of model parameters. The case $({\kappa }_1=4,{\kappa }_2=0)$ corresponds to the linear mean-reversion; the case $({\kappa }_1=4,{\kappa }_2=4)$ (generally the case with ${\kappa }_2={\kappa }_1/𝜃$) corresponds to the pure quadratic drift; the case $({\kappa }_1=4,{\kappa }_2=8)$ corresponds to the quadratic drift. The linear model with ${\kappa }_2=0$ implies the heavy right tail for the distribution of the volatility. In Fig. 1(b), we show the skewness of the volatility as function of ${\kappa }_2$. As ${\kappa }_2$ increases, the skewness of the volatility declines.

3.5.1. Unconditional distribution of returns

We make an interesting connection between the skewness of the volatility process and the kurtosis of returns distribution. Following Chap. 9 in Barndorff-Nielsen & Shiryaev (2015), we assume that the distribution of returns over period $\delta t$ is Gaussian conditional on $\sigma$ :

The unconditional distribution of returns under the SDD of the volatility is

The above integral cannot be computed explicitly. Instead, we compute the central moments of returns under the unconditional distribution where we exchange the integration order. It is clear that for odd n, $m(n)=0$. For the second and the fourth moments, we have

We compute the excess kurtosis of returns distribution using $k=3m(4)/m^2(2)-3$ with Eq. (3.39) as follows :

It follows that the kurtosis of the unconditional distribution in Eq. (3.44) is maximal when $b=0$ (${\kappa }_2=0$). The kurtosis is finite only if liner mean-reversion parameter ${\kappa }_1$ is higher than the total volatility of variance. The kurtosis declines when rate b (${\kappa }_2$) increases. In Fig. 1(c), we show the kurtosis for the three choices of ${\kappa }_1=1,4,8$. We see that for the each choice, the parameter ${\kappa }_2$ can be thought as a dampening parameter to reduce both the skewness of the steady-state PDF of the volatility and the kurtosis of the unconditional PDF of returns.

3.6. Moments of volatility process

Computation of moments for the log-normal SV with the quadratic drift is an open problem with no exact solution. We show that there exists a closed-form solution based on a truncation method, which we further adopt for the affine expansion of the MGF. We consider the mean-adjusted volatility process defined in Eq. (3.32). We define its power function $m_t^{(n)}$, $m_t^{(n)}=Y_t^n$, and the mth moment, ${\overline{m}}_{\tau }^{(n)}$, $n=0,1,2,\dots ,$ as follows :

Applying Itô’s lemma to $m_t^{(n)}$ under the dynamics in Eq. (3.33), we obtain

with $m_0^{(n)}=Y_0^n$ and $c(n)=\frac{1}{2}{𝜗}^{2}n(n-1)$. We notice a pattern of the powers of n which allows for a recursive solution as follows.

Proposition 3.3 (Moments of volatility process). The solution to moments in Eq. (3.45) can be presented as a matrix equation for an infinite-dimensional vector

where ${\partial }_{\tau }$ is the derivative w.r.t. $\tau$ and

An approximate solution to ODEs (3.47) is obtained by truncating the number of terms to $k^{\ast }$ and defining a finite dimensional vector of mth moments $M^{(1,k^{\ast })},$ $m=1,\dots ,k^{\ast },$ as follows:

The analytic solution to ODE system (3.48) is given by

where $\mathrm{\text{expm}}()$ is the matrix exponent, ⋅ and ⁻¹are the matrix product and inverse, respectively, and $I^{(k^{\ast })}$ is $k^{\ast }\times k^{\ast }$ identity matrix.

Proof. We first present Eq. (3.46) as follows :

Thus, moment ${\overline{m}}_{\tau }^{(n)}$ solves recursive equation (omitting subscript t)

Using the fact that ${\overline{m}}^{(0)}=1$ and $c(1)=0$ we can present Eq. (3.50) for $n=0,1,\dots$ as a matrix equation for infinite-dimensional column vector in Eq. (3.47). Clearly, that under regularity conditions (real parts of eigenvalues of ${\mathrm{\Lambda }}^{(0,\infty )}$ are negative), we must have ${lim}_{k\to \infty }{\partial }_{\tau }{\overline{m}}^{(k)}=0$. Thus, to make a finite-dimensional approximation of Eq. (3.47), we fix $k^{\ast }$ and assume

Finally, by substituting fixed values for ${\overline{m}}^{(0)}$ and ${\overline{m}}^{(k^{\ast }+1)}$, we present (3.47) for finite-dimensional vector of mth moments, $m=1,\dots ,k^{\ast }$, in Eq. (3.48). Solution in Eq. (3.49) is obtained by integration. □

In Fig. 2, we plot first four moments of the volatility process computed using Eq. (3.49) with truncation order $k^{\ast }=4$ (Fig. 2(a)) and $k^{\ast }=8$ (Fig. 2(b)). For comparison, we show estimates and $95\%$ confidence intervals obtained by MC simulations. We see that the low-order truncation with $k^{\ast }=4$ is consistent with MC estimates for the first and second moments. The high-order truncation with $k^{\ast }=8$ is consistent with the first four moments compared to MC estimates.

Remark 3.1. We note that if ${\kappa }_2=0$, the system (3.48) becomes lower tridiagonal so that it can be solved exactly for all moments up to $k^{\ast }$ by the sequential integration from $n=1$ to $n=k^{\ast }$.

3.7. Properties of quadratic variance process

We consider the QV process $I_t$ defined by $d{I}_t=d{\sigma }_t^{2}dt$. The QV is an important quantity for option pricing. In particular, the expected QV equals model-based fair value of the continuous-time variance swap, which could be used for model calibration.

Proposition 3.4 (Bounded moments). Under Assumption 3.1 and ${\kappa }_1𝜃>0,$ the QV is well defined under the MMA measure and $𝔼[{(I_t)}^p]<\infty$ for any $p,$ $|p|\ge 1$. If ${\kappa }_2\ge \beta ,$ then the QV is well defined under the inverse measure and $Ẽ[{(I_t)}^p]<\infty$ for any $p,$ $|p|\ge 1$.

Proof. As $𝔼[{(I_t)}^p]=𝔼[{({\int }_0^t{\sigma }_s)}^p]\le t^p𝔼[{sup}_{s\in [0,t]}{\sigma }_s^p]<\infty$, the first statement follows from Proposition 3.1. The second statement follows using the volatility dynamics under $\tilde{\mathbb{Q}}$ in Eq. (3.28). □

We define the annualized expected QV as follows :

Corollary 3.4 (Expected value). The expected value of the QV using $k^{\ast }$th order truncation in Proposition 3.3 is given by

where ${\widehat{\overline{m}}}^{(1)}$ and ${\widehat{\overline{m}}}^{(2)}$ are the first and second elements in vector ${\widehat{M}}^{(1,k^{\ast })}(\tau )$ computed using Eq. (3.54), and $O(k^{\ast })$ is the truncation error. The integrated moments of the volatility with the $k^{\ast }$-order truncation are computed by

The approximation term $O(k^{\ast })$ includes only the truncation error $O(k^{\ast })$ because Eq. (3.54) is computed analytically using matrix exponentiation. In Fig. 3, we present comparison of the analytical solution for the expected QV using Eq. (3.53) to the estimate and $95\%$ confidence interval based on MC simulations. The truncation-based analytic solution is consistent with MC estimates.

3.8. Monte Carlo discretization

We consider the SDE (3.19) for the volatility process ${\sigma }_t$ and introduce the log-volatility process $L_t=ln{\sigma }_t$. Applying Itô’s formula, we obtain the following dynamics for $L_t$ :

We consider a time horizon $T>0$ and an equidistant time grid $t_k=k\mathrm{\Delta }$, $\mathrm{\Delta }=\frac{T}{n}$, $0\le k\le n$. Backward Euler–Maruyama (BEM) scheme for $L_t$ is defined by

BEM scheme (also known as drift implicit Euler–Maruyama scheme) is based on Lamperti transform where original SDE (3.19) is transformed into an SDE (3.55) with constant diffusion coefficient and nonlinearity in the diffusion coefficient shifted into the drift of the process. We refer to Neuenkirch & Szpruch (2014), Alfonsi (2013), and Chassagneux et al. (2016) for further details.

We prove that SDE (3.56) has a unique solution in $L_{t_{k+1}}$ in order to make BEM scheme well-defined. Although BEM scheme in Eq. (3.56) cannot be solved analytically in ${\widehat{L}}_{t_{k+1}}$, we still can solve Eq. (3.56) numerically due to the smoothness and monotonicity of function $G(L):=L-\zeta (L)\mathrm{\Delta }$, using just few iterations of Newton algorithm, see Proposition 3.5.

Proposition 3.5. We introduce $G(l):=l-\zeta (l)\mathrm{\Delta }$. Then, for every $c\in ℝ,$ the equation $G(l)=c$ has a unique solution in $l\in ℝ$.

Proof. As $G(l)$ is continuous, $G^{\prime }(l)>0$ and ${lim}_{l\to \pm \infty }G(l)=\pm \infty$, the result follows. □

Theorem 3.8. Backward Euler–Maruyama scheme for log-volatility process $L_t$ in SDE (3.55) has a strong convergence rate of 1.

Proof. The proof of the theorem is based on Proposition 3 of Alfonsi (2013). We present it in Proposition 3.6.

Proposition 3.6. We fix $p\ge 1$ and assume that

Then, there exists a constant $K_p>0$ such that ${(𝔼{sup}_{t\in [0,T]}|{\widehat{L}}_t-L_t{|}^p)}^{1/p}\le K_p\mathrm{\Delta }$.

Instead of Eq. (3.57), we prove slightly stronger estimate in Eq. (3.58). The latter follows from the combination of Hölder inequality and Proposition 3.1. We omit straightforward details to obtain

 □

Corollary 3.5. MC discretization scheme of the model dynamics (3.12) under the MMA measure using log-price $X_t$ in SDE (3.16) is given by

with $X_{t_0}=ln{S}_0,$ $L_{t_0}=ln{\sigma }_0,$ $I_{t_0}=I_0,$ and where $W^{(0)}$ and $W^{(1)}$ are independent Brownian motions.

We note that the first two dynamics are defined in $(-\infty ,+\infty )$, which does not require extra handling of boundary conditions, in contrast to some affine models for the variance process.

4. Affine Expansion of Moment Generating Function

We consider the valuation problem of a derivative security with payoff function $u(X,I)$ settled at maturity time T using log-price process $X_T$ defined in Eq. (3.16). We denote current time by t and introduce time-to-maturity variable $\tau =T-t$. We denote the undiscounted value function of this derivative security by $U(\tau ,X,I,Y)$.

To solve the general valuation problem under both $\mathbb{Q}$ and $\tilde{\mathbb{Q}}$, we parametrize the value function using binary parameter p, $p\in \{1,-1\}$, respectively, as follows :

where the expectation $𝔼$ under the MMA measure $\mathbb{Q}$ is computed using dynamics (3.34) and the expectation $\widetilde{𝔼}$ under the inverse measure $\tilde{\mathbb{Q}}$ is computed using dynamics (3.36). Here we apply the fact that the dynamics are Markovian under the both measures.

Theorem 4.1. The value function $U(\tau ,X,I,Y;p)$ solves the following PDE on the domain $[0,T]\times ℝ\times {ℝ}_{+}\times (-𝜃,\infty )$

where diffusive operators ${{ℒ}}^{(Y;p)},$ ${{ℒ}}^{(X;p)},$ and ${{ℒ}}^{(I;p)}$ are defined as follows:

Proof. We note that the classic Feynman–Kac formula (for an example, see Sec. A.17 in Musiela & Rutkowski (2009)) cannot be applied directly to Eq. (4.1) because the model dynamics (3.34) and (3.36) do not satisfy the linear growth condition. In Sec. A.2 we prove that the existence and the uniqueness is guaranteed by Theorem A.1. This theorem states that if the covariance matrix of diffusive operator ${ℒ}\equiv {{ℒ}}^{(Y;p)}+{{ℒ}}^{(X;p)}+{{ℒ}}^{(I;p)}$ is nonnegative definite and the state variables neither explode nor leave an open domain, then there exists a unique solution to the PDE representation of the conditional expectation under corresponding model dynamics.

Now using the diffusive operators in the valuation PDE (4.2), we obtain that the diffusion matrix a in Eq. (A.14) equals

with eigenvalues given by Since $\beta <𝜗$, we observe that $0={\lambda }_1<{\lambda }_2<{\lambda }_3$, hence matrix a is only nonnegative definite and its smallest eigenvalue is always zero.

We now set the domain D of state variables $\{X_t,Y_t,I_t\}$ in PDE (4.2) by $D=ℝ\times (-𝜃,+\infty )\times {ℝ}_{+}$. As we showed in Theorem (3.2), under the MMA measure (3.34) and the inverse measure (3.36), the boundary points $\{0,+\infty \}$ are unattainable for the volatility process ${\sigma }_t$, hence solution for mean-adjusted volatility process $Y_t$ and QV $I_t$ cannot explode or leave D before T. Same applies to $X_t$ as $\{0,+\infty \}$ are unattainable for ${\sigma }_t$. Thus, condition 1 of Theorem A.1 is satisfied. As boundedness condition A.1 is satisfied by construction, the statement of Theorem 4.1 follows from Theorem A.1 stated in Sec. A.2. □

4.1. Moment generating function

We introduce the MGF of the three model variables: the log-price $X_t$, the QV $I_t$, and the mean-adjusted volatility driver $Y_t={\sigma }_t-𝜃$. We extend our analysis using the transform variable in $Y_t$ to generalize our solution method for all the three state variables in the model. We denote the MGF by $G(\tau ,X,I,Y;\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta };p)$ using respective complex-valued transform variables $\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta }\in ℂ$

Using model dynamics (3.34) under the MMA measure and model dynamics (3.36) under the inverse measure, respectively, the MGF solves the following PDE on the domain $[0,T]\times ℝ\times {ℝ}_{+}\times (-𝜃,\infty )$ :

with operators defined in Eq. (4.3). First, we establish a sufficient condition for the existence of the MGF for the three state variables.

Theorem 4.2. Given the transform variable $\mathrm{\Phi }={\mathrm{\Phi }}_R+i{\mathrm{\Phi }}_I\in ℂ,$ the MGF of the log-price $X_{\tau }$

exists for ${\mathrm{\Phi }}_R\in (-1,0),$ if $Z_{\tau }$ is a martingale under the MMA measure. By Theorem 3.7(1), the necessary condition is ${\kappa }_2\ge \beta$.

Similarly, the MGF of the log-price $X_{\tau }$

exists for ${\mathrm{\Phi }}_R\in (0,1)$ if $R_{\tau }$ is a martingale under inverse measure. By Theorem 3.7(2), the necessary condition is ${\kappa }_2\ge 2\beta$.

Proof. As $Z_T$ is $\mathbb{Q}$-martingale, the result immediately follows from Jensen’s inequality, as the function $g(y):=y^{|{\mathrm{\Phi }}_R{|}^{-1}}$ is convex for $|{\mathrm{\Phi }}_R|\in (0,1)$. Result for the inverse measure follows analogously. □

Theorem 4.3. Given transform variable $\mathrm{\Psi }={\mathrm{\Psi }}_R+i{\mathrm{\Psi }}_I\in ℂ,$ the MGF of QV $I_{\tau }$

exists for ${\mathrm{\Psi }}_R<0$ if ${\kappa }_2>𝜗\sqrt{-2{\mathrm{\Psi }}_R}$.

Proof. We denote $K_t=e^{I_t}$ which satisfies SDE $d{K}_t=K_{t}d{I}_t=K_t{\sigma }_t^{2}dt$ and set $p=-{\mathrm{\Psi }}_R$. Feynman–Kac formula allows to recast the solution $U(\tau ,\sigma ,K)=𝔼(K_T^p|K_t=K,{\sigma }_t=\sigma )$ as the solution of the Cauchy problem

where operator ${{ℒ}}^{(\sigma ,K)}$ is defined by

We build a supersolution $Ū(\tau ,\sigma ,K)$ to (4.7) using the following ansatz :

which yields

We show there exists $a\ge 0$ such that

We notice that inequality (4.9) holds if there exist $a>0$ such that $\frac{1}{2}{𝜗}^2{a}^2-{\kappa }_{2}a+p<0$. As latter defines a parabola in a, we can always find required positive a, if ${\kappa }_2>𝜗\sqrt{2p}=𝜗\sqrt{-2{\mathrm{\Psi }}_R}.$ □

Theorem 4.4. Given the transform variable $\mathrm{\Theta }={\mathrm{\Theta }}_R+i{\mathrm{\Theta }}_I\in ℂ,$ the MGF of the mean-adjusted volatility process $Y_{\tau }$

exists for ${\mathrm{\Theta }}_R<0,$ if ${\kappa }_2>\frac{{𝜗}^2}{2}|{\mathrm{\Theta }}_R|,$ $p=1$ and ${\kappa }_2-\beta >\frac{{𝜗}^2}{2}|{\mathrm{\Theta }}_R|,$ $p=-1$.

Proof. We set $c=-{\mathrm{\Phi }}_R$ and consider

Using Feynman–Kac formula we can rewrite $U(\tau ,Y;c;p)$ as the solution of the Cauchy problem

where operator ${{ℒ}}^{(Y;p)}$ is defined by Eq. (4.3). We build a supersolution $Ū(\tau ,Y;c;p)$ to (4.11) using ansatz $Ū=e^{aY+b\tau },a\ge c,b\ge 0$. Inserting $Ū$ into Eq. (4.11), we have where coefficients are defined in (4.3). We show there exists $a\ge c$ such that We notice that last inequality is satisfied if there exist $a\ge c>0$ such that $-{\kappa }_2^{(p)}a+\frac{1}{2}{𝜗}^2{a}^2<0$. Dividing by $a>0$ we conclude that it is satisfied if ${\kappa }_2^{(p)}>\frac{{𝜗}^2}{2}c.$ Rewriting the last condition in terms of ${\kappa }_2$, $\beta$ using (4.3), completes the proof. □

4.2. Affine expansion

Given that coefficients in PDE (4.6) are nonaffine in state variable Y, the PDE cannot be solved by assuming an exponential affine-solution with a finite number of coefficients. Instead, we can show that the solution to the PDE (4.6) can be represented by an exponential-affine ansatz $E^{[\infty ]}$ using an infinite dimensional vector $A(\tau )=\{A^{(k)}(\tau ;\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta };p)\}$, $k=0,1,\dots ,\infty$, as

where vector function $A(\tau )$ solve an infinite-dimensional system of quadratic ODEs : where $M^{(k)}$, $L^{(k)}$, $H^{(k)}$ are infinite-dimensional sparse matrices and vectors.

We note the similarity between the current ansatz (4.13) for problem (4.6) and the problem of computing the moments of the volatility process addressed in Proposition 3.3. There we first represent the recursive solution to the moments as an infinite dimensional vector in Eq. (3.47). We then obtain an approximate solution by a truncation of infinite series and specifying a boundary condition to reduce the problem to a finite-dimensional vector (3.48), which can be solved using analytical methods. Here we follow this insight for solving PDE (4.6) and term our (truncated) solution as the affine expansion for a given order.

4.2.1. First-order expansion

Theorem 4.5 (First-order affine expansion). The MGF in Eq. (4.5) can be decomposed as follows:

where $E^{[1]}$ is the leading first-order term and $R^{[1]}$ is the remainder term, such that $E^{[1]}(\tau =0)=1$ and $R^{[1]}(\tau =0)=0$. The leading term $E^{[1]}$ is given by the exponential-affine formwhere vector function $A(\tau )=\{A^{(k)}(\tau ;\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta };p)\},$ $k=0,1,2,$ solve the quadratic differential system as a function of $\tau :$with the initial condition $A(0)={(0,-\mathrm{\Theta },0)}^{\top }$.

The remainder term $R^{[1]}$ solves the following PDE (omitting arguments) :

with operator ${{ℒ}}^{(Y)}$ defined in Eq. (4.3) and residual $F^{[1]}$ defined by

Proof. We take as an ansatz the leading term $E^{[1]}(\tau ,Y;\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta };p)$ in decomposition (4.16) and substitute it into the PDE (4.6). To obtain the ODEs (4.17) for $A^{(k)}$, we match terms with $Y^k$, $k=0,1,2$. To obtain the PDE for the remainder term $R^{[1]}$, we account for higher powers of $Y^k$, $k=0,1,2$. Finally, initial conditions for ODEs (4.17) and PDE (4.18) are set to reproduce conditions of Eq. (4.6). □

Proposition 4.1. We take $\mathrm{\Psi }=\mathrm{\Theta }=0$ and $\Re \mathrm{\Phi }=-p/2$, $p\in \{-1,1\}$. Assuming that hypotheses 1, 2, 3 listed in Theorem 4.7 hold, the continuous solution $A(\tau )$ in Eq. (4.17) exists on $[0,+\infty )$.

Proof. The proof is identical to that of Theorem 4.7 for the second-order expansion. We skip it for brevity. □

In Fig. 4, we show the real and imaginary parts of ODE solutions $A(\tau )$ and $E^{[1]}$. We see that functions $A^{(1)}$ and $A^{(2)}$ quickly reach an equilibrium point, while $A^{(0)}$ is the nonstationary part which contributes to the leading term $E^{[1]}$.

Corollary 4.1 (Estimate of the first-order remainder term). We obtain the following estimate for the remainder term $R^{[1]}$ in Eq. (4.18) :

where $M_{\sigma }^{(n)}(\tau )$ is the nth central moment of the volatility defined by

Proof. The proof is analogous to the second-order expansion in Proposition 4.3. □

Corollary 4.2. First-order affine expansion for the MGF (4.5) is obtained by using the leading term $E^{[1]}$ in Eq. (4.15) :

where the approximation error arises from the truncation of the infinite dimensional affine anzats in Eq. (4.13) with the error magnitude estimated using Eq. (4.20).

In Fig. 6, we illustrate that the first-order approximation (4.22) for the MGF produces valid and accurate PDFs for all three state variables.

4.2.2. Second-order expansion

Theorem 4.6 (Second-order affine expansion). The MGF in Eq. (4.5) can be decomposed as follows:

where $E^{[2]}$ is the leading term such that $E^{[2]}(\tau =0)=1$ and $R^{[2]}(\tau =0)=0$ is the remainder term.

The leading term $E^{[2]}(\tau ,Y;\mathrm{\Phi },\mathrm{\Psi };p)$ is given by the exponential-affine form

where vector function $A(\tau )=\{A^{(k)}(\tau ;\mathrm{\Phi },\mathrm{\Psi },\mathrm{\Theta };p)\},$ $k=0,\dots ,4,$ solves the system the quadratic differential system as a function of $\tau :$with initial condition $A(0)={(0,-\mathrm{\Theta },0,0,0)}^{\top }$.

The remainder term $R^{[2]}$ solves the following PDE (omitting arguments)

with operator ${{ℒ}}^{(Y)}$ defined in Eq. (4.3) and function $F^{[2]}$ defined by

Proof. The proof follows by analogy to the first-order expansion in Theorem 4.5. □

Remark 4.1. We emphasize that sparse matrices $M^{(k)}$ in the quadratic term do not depend on the measure parameter p, while the linear $L^{(k)}(p)$ and free terms $H^{(k)}(p)$ depend on p in ODEs (4.17), (4.25) for the first-and the second-order expansion.

We now present the result on the existence of continuous solution to Eq. (4.25) in Proposition 4.2 and discuss it further in Sec. 4.2.3.

Proposition 4.2. We take $\mathrm{\Psi }=\mathrm{\Theta }=0$ and $\Re \mathrm{\Phi }=-p/2$, $p\in \{-1,1\}$. Assuming that hypotheses 1,2,3 listed in Theorem 4.7 hold, the continuous solution $A(\tau )$ in Eq. (4.25) exists on $[0,+\infty )$.

Proof. See Theorem 4.7. □

In Fig. 5, we show the real and imaginary parts of $A(\tau )$ for the second-order expansion and the leading term $E^{[2]}$. We see that, similarly to the first order expansion illustrated in Fig. 4, functions $A^{(1)},A^{(2)},A^{(3)},A^{(4)}$ reach an equilibrium point, while $A^{(0)}$ is the nonstationary part which contributes to the leading term $E^{[2]}$. By comparing the low order terms $A^{(0)},A^{(1)},A^{(2)}$ between the first and the second expansions, we notice no visible difference, which suggest that the expansion is rather recursive and higher order terms make insignificant contribution to the low order terms. We also notice that terms $A^{(3)}$ and $A^{(4)}$ are very small in magnitude compared to the first three terms, which suggest that higher order term only provide a small marginal contribution.

Corollary 4.3 (Estimate of the second-order remainder term). We obtain the following estimate for the remainder term $R^{[2]}$ which solves Eq. (4.26) :

where $M_{\sigma }^{(n)}$ is the nth central moment of the volatility defined by (4.21).

Proof. The formal solution for the remainder term $R^{[2]}$ solving problem (4.26) is obtained by applying the Feynman–Kac formula and is given by

where $P(t,Y;Y^{\prime })$ is PDF of $Y^{\prime }=Y(t)$ conditional on Y. P solves PDE (4.2), subject to initial condition $P(0,Y;Y^{\prime })=\delta (Y^{\prime }-Y)$. Using Eq. (4.29), we obtain Here, we apply the bound on the MGF $|E^{[2]}(\tau -t,Y^{\prime };\mathrm{\Phi },\mathrm{\Psi })|\le 1$ and integrate out the log-return X and the QV I. Next we substitute the polynomial function $F^{[2]}$ in Eq. (4.27). We use the continuity of functions $A^{(k)}$, $k=1,2,3,4$, and apply the first mean value theorem for definite integrals to estimate the time integral. Finally, we approximate the expected moments of the mean-adjusted volatility by its steady-state nth order central moments $M_{\sigma }^{(n)}$ specified in Eq. (4.21). □

Corollary 4.4. Second-order affine expansion for the MGF (4.5) is obtained by using the leading term $E^{[2]}$ in Eq. (4.24) as

where the approximation error arises from the truncation of the infinite dimensional affine anzats in (4.13) with the error magnitude estimated by Eq. (4.28).

In Fig. 6, we illustrate that the second-order approximation (4.31) for the MGF produces valid and accurate PDFs for all three state variables.

4.2.3. Solution of quadratic differential systems

It is established that the global existence problem for the solution of quadratic differential systems is nontrivial, see Coppel (1966), Dickson & Perko (1970), Jacobson (1977), Baris et al. (2008) among others. We can show that the matrices of quadratic terms in Eqs. (4.17) and (4.25) are indefinite because their eigenvalues are of both signs. Thus, we can provide a result for global existence that is only conditional by nature. We focus on important case of option valuation on underlying assets with $\Re \mathrm{\Phi }=-p/2$, $p\in \{-1,1\}$, and $\mathrm{\Psi }=\mathrm{\Theta }=0$.

Dickson & Perko (1970) show that complex-valued solution of quadratic ODEs such as $A(\tau ;\mathrm{\Phi })$ in Eq. (4.17) or Eq. (4.25) has maximal interval of existence $[0,{\tau }_{+}(\mathrm{\Phi };A))$, where ${\tau }_{+}(\mathrm{\Phi };A)$ denotes its blow-up time.

Theorem 4.7. Assume following conditions are satisfied:

(1)	${\tau }_{+}(\mathrm{\Phi };\Re A)\ge {\tau }_{+}(\Re \mathrm{\Phi };A),$ $\mathrm{\Phi }\in ℂ,$ i.e. real part of $A(\tau ;\mathrm{\Phi })$ cannot blow up before $A(\tau ;\Re \mathrm{\Phi })$.
(2)	${\sum }_{k=0}^{2m}{A}_k(\tau ;\mathrm{\Phi };p)Y^k\to +\infty$ as $\tau \uparrow {\tau }_{+}(\mathrm{\Phi }),$ $\mathrm{\Phi }\in ℝ,$ i.e. leading term of the affine expansion does not vanish as we approach blow-up time when transform variable $\mathrm{\Phi }$ is real.
(3)	${lim}_{\tau \uparrow {\tau }_{+}(\mathrm{\Phi })}{R}^{[n]}(\tau ;\mathrm{\Phi })>-\infty ,$ $\mathrm{\Phi }\in ℝ,$ i.e. remainder term is uniformly bounded from below.

If MGF in Eq. (4.6) is finite with $G({\tau }_0;X,I,Y;\mathrm{\Psi }=\Re \mathrm{\Phi },\mathrm{\Psi }=0,\mathrm{\Theta }=0)<\infty ,$ then continuous solution $A(\tau ;\mathrm{\Phi },\mathrm{\Psi }=0,\mathrm{\Theta }=0)$ of Eq. (4.17) and of Eq. (4.25) exists on $[0,{\tau }_0)$.

Proof. See Sec. A.4. □

By Theorem 4.7, we obtain that continuous solution $A(\tau )$ for the first-and second-order affine expansions in Eqs. (4.16) and (4.24), respectively, exists on $[0,+\infty )$ if 1, 2, 3 hold.

4.3. Properties of affine expansions

We now consider the key properties of the affine expansion of the MGF. For brevity, we focus on the properties of the second-order expansion, which is our key development. For the first order expansion, martingale conditions stated in Proposition 4.3 hold, and the consistency with the expected values (but not with their variances) of three state variables is ascertained in the same way as for the second-order expansion.

Proposition 4.3 (Martingale conditions). Assuming that hypotheses 1, 2, 3 in Theorem4.7 hold, the second-order expansion of the MGF in (4.31) with $E^{[2]}$ defined in Eq. (4.24) satisfies the martingale conditions for log-return $X_{\tau }:$

Proof. For $\mathrm{\Phi }=\mathrm{\Psi }=\mathrm{\Theta }=0$, we obtain $H^{(p)}\equiv 0$ in (4.25). Given zero initial conditions, $A^{(k)}(\tau ;\mathrm{\Phi }=0,\mathrm{\Psi }=0,\mathrm{\Theta }=0;p)\equiv 0$ due to the uniqueness of the solution. Similarly, for $\mathrm{\Phi }=-1$ and $\mathrm{\Theta }=\mathrm{\Psi }=0$, we obtain that $H^{(1)}\equiv 0$, hence $A^{(k)}(\tau ;\mathrm{\Phi }=-1,\mathrm{\Psi }=0,\mathrm{\Theta }=0;p=1)\equiv 0$. Similarly, for $\mathrm{\Phi }=1$ and $\mathrm{\Psi }=0$ we obtain that $H^{(-1)}\equiv 0$, hence $A^{(k)}(\tau ;\mathrm{\Phi }=1,\mathrm{\Psi }=0,\mathrm{\Theta }=0;p=-1)\equiv 0$. □

Proposition 4.4 (Volatility moments). The MGF in Eq. (4.31) with the second-order affine term $E^{[2]}$ in (4.24) reproduces exactly the expected value and the variance of the volatility for ${\kappa }_2=0$.

Proposition 4.5 (Log-price moments). The MGF in Eq. (4.31) with the second-order affine term reproduces exactly the expected value and the variance of the log-price for ${\kappa }_2=0$.

Proposition 4.6 (QV moments). The MGF in Eq. (4.31) with the second-order affine term reproduces exactly the expected value and the variance of QV for ${\kappa }_2=0$.

Proof. Propositions 4.4 and 4.5 are proved in Secs. A.3.1 and A.3.2, respectively. Proof of Proposition 4.6 is similar to the proof of Proposition 4.5. □

Remark 4.2. In Remark 3.1 we noted that the system of volatility moments (3.48) is solved exactly for ${\kappa }_2=0$, otherwise the solution is approximate. Thus, we cannot generalize above statements for ${\kappa }_2>0$. However, it is our hypothesis that moments obtained using (4.24) correspond to moments computed in (3.48) using truncation order $k^{\ast }=4$ and that truncation error as function of ${\kappa }_2$ is very small.

In Fig. 6, we show PDFs computed using the inversion of first-order solution $E^{[1]}$ in Eq. (4.16) and second-order solution $E^{[2]}$ in Eq. (4.24) for the three state variables in model dynamics (3.34). The blue histogram is computed using MC simulations of the joint model dynamics. We see that both first-and second-order expansions produce valid PDFs and match consistently histograms produced by the MC simulations. The first-order expansion may not be very accurate for the QV and the volatility variables because it is not consistent with variances of these variables. The second-order expansion is consistent with variances and co-variances of all state variables and it accurately matches PDFs computed using MC simulations.

5. Option Valuation

We apply the zero-drift log-price process $X_t$ in Eq. (3.16) and Eqs. (3.17) and (3.18) for spot and futures prices, respectively, to consider option valuation on both spot and futures underlyings. We denote by $P_T$ the price of either spot or futures asset

where $\bar{\mu }(T)=\bar{r}(0,T)$ for spot price and $\bar{\mu }(T)=\bar{r}(0,T)-\bar{q}(0,T)$ for futures.

As we conclude in Corollaries 4.2 and 4.4, we apply affine expansion given by either the first order in Eq. (4.16) or the second order in Eq. (4.24) as analytic solutions to the MGF under the dynamics (3.16). Given that the MFG is available analytically, we then apply Lewis–Lipton approach for valuation of vanilla options using Fourier transform, see Lewis (2000), Lipton (2001, 2002).

5.1. Vanilla calls and puts

Using Eq. (5.1), we represent put and call payoff functions for strike price K using capped payoff by

Thus, we evaluate capped payoff $u(P_T)=min\{P_T,K\}$ under MMA measure $\mathbb{Q}$ by

Proposition 5.1. The valuation formula for capped payoff is given by

where $X^{\ast }=ln(S_0/K)+\overline{\mu }(T)$ is log-moneyness, $E^{[1]}$ and $E^{[2]}$ are given in Eqs. (4.16) and (4.24), respectively.

Proof. Using affine expansion for MGF G with either first-order in Eq. (4.22) or second-order in Eq. (4.31), the PDF of $X_{\tau }$ is computed by Fourier inversion

For valuation problem in Eq. (5.3), we obtain

Assuming that inner integrals are finite, we exchange the integration order as

where $û(\mathrm{\Phi })$ is the transformed payoff function where $X^{\ast }=ln(S/K)+\overline{\mu }(T)$ is log-moneyness, $k^{\ast }=lnK-\overline{\mu }(T)$ with the first integral being finite for $\Re [\mathrm{\Phi }]>-1$ and the second integral being finite for $\Re [\mathrm{\Phi }]<0$. Integral in Eq. (5.7) is finite for $-1<\Re [\mathrm{\Phi }]<0$. For $\mathrm{\Phi }=iy-1/2$, we get Eq. (5.4)  □

As a result, calls and puts on spot underlying are valued using Eqs. (5.2)

The value of call option on future underlying, in turn, equals where $f_0(T)$ is the futures and $X^{\ast }=ln(f_0(T)/K)$ denotes log-moneyness. Hence, given a set of futures prices quoted in the market, we can evaluate options on this set of futures without need to compute the convenience yield explicitly.

We emphasize that in valuation equations (5.9) we implicitly assume that the price process $Z_t=exp(X_t)$, $t\in (0,T]$, is $\mathbb{Q}$-martingale under the MMA measure, otherwise formulas (5.9) are not valid. For that we need a restriction ${\kappa }_2\ge \beta$, see Theorem 3.7.

5.2. Inverse calls and puts

Using generic price process of underlying spot and futures given in Eq. (5.1), we represent the payoffs of inverse calls and puts by

As a result, we need to evaluate option on the inverse capped payoff under $\tilde{\mathbb{Q}}$

Proposition 5.2. The valuation formula (5.6) for the inverse capped payoff becomes

where $X^{\ast }=ln(S_0/K)+\overline{\mu }(T)$ is log-moneyness, $E^{[1]}$ and $E^{[2]}$ are given in Eqs. (4.16) and (4.24), respectively.

Proof. As in Eq. (5.7), we compute the transform $û(\mathrm{\Phi })$ of the inverse capped option

where $X^{\ast }=ln(S/K)+\overline{\mu }(T)$ is the log-moneyness and $k^{\ast }=lnK-\overline{\mu }(T)$, with the first integral being finite for $\Re [\mathrm{\Phi }]>0$ and the second integral being finite for $\Re [\mathrm{\Phi }]<1$. Integral in Eq. (5.14) is finite for $0<\Re [\mathrm{\Phi }]<1$. We set $\mathrm{\Phi }=iy+1/2$, then  □

As a result, calls and puts are valued using capped payoff (5.12) :

Hereby, we need to compute

with the last equality following from Eq. (4.33). Given the term structure of futures prices $f_0(T)$ in Eq. (3.18), the put option on the future is with log-moneyness $X^{\ast }=ln(f_0(T)/K)$. We emphasize that Eq. (5.16) holds only if $\mathbb{Q}\sim \tilde{\mathbb{Q}}$, which requires restriction ${\kappa }_2\ge \beta$, see Theorem 3.6. In addition, we assume that the inverse process $R_t=exp(-X_t)$, $t\in (0,T]$, is a true $\tilde{\mathbb{Q}}$-martingale. For this we need ${\kappa }_2\ge 2\beta$ by Theorem 3.7.

5.3. Options on quadratic variance

We consider a call option on QV with the payoff function

where K is the strike in terms of the annualized variance.

Proposition 5.3. The value of the call option on the QV under $\mathbb{Q}$ is given by

whereprovided $\Re [\mathrm{\Psi }]<0,$ with $E^{[1]}$ and $E^{[2]}$ given in Eqs. (4.16) and (4.24), respectively.

Proof. Consider an option on QV $I_T$ with terminal payoff $u(I_T)$. Using Sepp (2008), we find that option value is given by Eq. (5.20), where $û(\mathrm{\Psi })$ is the transformed payoff function

For a call option on QV with strike K, provided $\Re [\mathrm{\Psi }]<0$, Eq. (5.22) becomes

 □

Sepp (2008) derives transforms $û(\mathrm{\Psi })$ for typical payoff on the QV including puts and calls on the square root of the QV. Assuming that the inverse measure $\tilde{\mathbb{Q}}$ is a martingale measure for generic spot and futures underlying $P_t$ defined in Eq. (5.1), the inverse option on the QV is given under $\tilde{\mathbb{Q}}$ by

Similarly to Eq. (2.5), values $Ũ(\tau ,X,Y,I)$ and $U(\tau ,X,Y,I)$ in Eqs. (5.23) and (5.20), respectively, are related by

Proposition 5.4. The value of the call option on the QV under the inverse measure $\tilde{\mathbb{Q}}$ is given by

provided $\Re [\mathrm{\Psi }]<0,$ with $E^{[m]}$ and $û(\mathrm{\Psi })$ defined as in Eq. (5.21).

Proof. Using affine expansion for MGF G with either first-order in Eq. (4.22) or second-order in Eq. (4.31), we compute the joint density of $\{X_{\tau },I_{\tau }\}$ using inverse Fourier transform by

As a result, we obtain

Exchanging the integration order, we have

where $û(\mathrm{\Phi },\mathrm{\Psi })$ is the transformed payoff function provided $\Re [\mathrm{\Psi }]<0$. We integrated out the complex exponential in Eq. (5.26) using results of Sepp (2007). Thus, Eq. (5.25) becomes  □

6. Model Implementation and Calibration to Options Data

6.1. Implementation

In practical applications, either for model calibration or for options valuation, we need to compute call and put option prices on a grid consisting of several maturities and strikes, also referred to as an option chain. We consider grid with M maturities and J strikes. For option valuation using the MGF, we apply the affine expansion given by either the first order in Eq. (4.16) or the second order in Eq. (4.24).

First, we fix the space grid with size P (typically, $P\approx 500$) of transform variable $\mathrm{\Phi }$ for the log-price. We solve the system of ODEs numerically for either the first-order expansion in Eq. (4.16) or for the second-order expansion in Eq. (4.24) in time up to the last maturity time. The computational cost of this step is

where $N_{max}$ is the number of time steps for ODE solver (typically we have one-two steps per one week) and ${𝒪}$ is the number of arithmetic (elementary) operations.

Second, we compute values of the MGF on the grid of $\mathrm{\Phi }$ for each maturity slice. We then compute values of capped payoffs for the MMA measure in Eq. (5.4) or for the inverse measure in Eq. (5.13) using Simpson rule for numerical integration. We then value all options in a given chain. The computational cost of this step is

The main difference with numerical implementation of affine models is that the first step can be computed with cost close to ${𝒪}(P\times M)$ because the solution to the MGF can be computed analytically without solving ODEs numerically. The cost ratio between the implementation of the log-normal SV model and an affine SV model compares favorably when the number of maturities becomes large with M close to $N_{max}$, and when the number of strikes J becomes large as well.

6.2. Model calibration

To illustrate that the log-normal SV model can handle different market regimes, we use time series of option prices on Bitcoin observed on Deribit exchange from April 2019 to October 2023. We focus on short-dated options with expiries including one and two weeks and one month, because these are the most liquid and actively traded expiries. We use time series of options chains traded on Deribit. We fix calibration every week with option prices samples each Friday at 10:00 UTC. We include puts and calls with absolute values of option deltas between $[-0.10,-0.50]$ and $[0.10,0.50]$, respectively.

For the stability of calibration, we estimate mean-reversion parameters ${\kappa }_1$ and ${\kappa }_2$ beforehand. The reason is that mean-reversion parameters and volatility-of-volatility with volatility beta are co-dependent so that it is possible to produce similar smiles and skews using different values of these parameters. We follow the procedure in Sepp & Rakhmonov (2023) for estimation of parameters ${\kappa }_1$ and ${\kappa }_2$ by fitting the model implied auto-correlation of the log-normal volatility process to the empirical auto-correlation. This is efficient because the model auto-correlation function is determined solely by mean-reversion parameters. The estimated values of ${\kappa }_1$ and ${\kappa }_2$ are the following ${\widehat{\kappa }}_1=2.21$ and ${\widehat{\kappa }}_2=2.18$, which are kept fixed through the entire period (in practice we would adjust these estimate regularly using most recent data).

Thus, we need to estimate four model parameters ${\sigma }_0,𝜃,\beta ,𝜀$. Parameters ${\sigma }_0$ and $𝜃$ control the level and the slope of at-the-money (ATM) implied volatilities, and $\beta$ and $𝜀$ control the skewness and the convexity of implied volatilities, respectively. We apply a stronger constraint ${\kappa }_2\ge 2\beta$ that guarantees martingale property for the inverse process $R_t$ under $\tilde{\mathbb{Q}}$ using Theorem 3.7. To compute model prices we apply the second-order affine expansion in Eq. (5.9). For each set of sampled option prices, we minimize the weighted mean squared error (WMSE) between model implied volatilities, ${\sigma }_n^{\mathrm{\text{model}}}(T,K)$, and market mid-point implied volatilities, ${\sigma }_n^{\mathrm{\text{implied}}}(T,K)$, where K and T are option’s strike and maturity, respectively, as follows^d :

with weight $w_n(T,K)={{𝒱}}_n(T,K)$ set to option Black–Scholes vega ${{𝒱}}_n(T,K)$.

In Fig. 7(a), we show the time series with the quality of the model fit which is computed using the square root of average squared differences between model and market implied volatilities. As a benchmark, we show the average spread between implied volatilities of ask and bid quotes for options in the calibration set. We observe that most of the times, the average model mis-pricing is within the average bid-ask spread. We see that year 2020 of COVID-19 pandemic is characterized by high volatility and high spreads, yet the model is able to fit these markets well. In recent years of 2022 and 2023, bid-ask spreads declined while the model error (the average difference between market and model implied volatility) became less than $1\%$ most of the times.

In Figs. 7(b)–7(d), we show the time series of calibrated values of initial volatility ${\sigma }_0$ and mean volatility $𝜃$, volatility beta $\beta$, and volatility-of-volatility $𝜀$, respectively. We observe downward trend in the level of implied volatility, falling from average levels of $80\%$ in years 2021 and 2022 to most recent averages of about $40\%$. At the same time, volatility beta and volatility-of-volatility remain range bound suggesting than relative value of calls and puts is intact despite falling level of the overall volatility. The volatility beta fluctuates between positive and negative values most of the times following the market sentiment. When the sentiment is bullish, calls are bid and volatility beta becomes positive and, vice versa, during bearish sentiment, puts are bid and the volatility beta becomes negative. Thus, as we discuss in Sec. 3.3, the conventional SV models may not be arbitrage-free for the valuation of options on cryptocurrencies when return–volatility correlation become positive.

In Fig. 8, we show the quality of model fit to Bitcoin options traded on Deribit exchange on 20 June 2023. On this day, fitted model parameters are the following :

In Fig. 8, we show the quality of model fit on this date. We observe positive skewness of implied volatilities, which results in positive volatility beta. Continuous black line shows model implied volatilities computed using the second-order affine expansion in Eq. (5.9). Bid and ask displays the market bid and ask quotes, the ATM is the ATM mid-point volatility. MSE denotes the mean squared error between the model implied volatilities and the mid of market bid-ask implied volatilities. We see that the log-normal SV model calibrated to Bitcoin options is able to capture the market implied skew very well across most liquid maturities with only four parameters. The average MSE is about $0.5\%$ for longer maturities, which is within the bid-ask spread. Calibration to the ATM region can be further improved using a term structure of the mean volatility $𝜃$. In a practical setting, the SV model dynamics can be augmented with a local volatility (Lipton (2002)) to accurately fit the implied volatility surface.

6.3. Comparison of valuation under MMA and inverse measures

Since Bitcoin options, which are traded on the Deribit exchange, have inverse payoffs, it is important that there is no-arbitrage between the MMA and the inverse measures. In Fig. 9, we compare the model implied volatilities computed from prices of vanilla and inverse options valued under the MMA and inverse measures, respectively. We use the second-order affine expansion with model parameters reported in Eq. (6.4) for computing model implied volatilities for the same maturities as in Fig. 8. As benchmark, we use the MC simulation of the model dynamics (3.12) using scheme in Eq. (3.59) with the number of paths equal $400,000$ and daily time steps. We show the 95% confidence interval of MC estimates as dashed lines like bid/ask lines

We observe that option values computed using the second order affine expansion under the MMA and inverse measures are very close to each other with the differences being within the numerical accuracy of an ODE solver and Fourier inversion (of order $10^{-4}$). The difference in terms of implied volatilities does not exceed $10^{-4}$.

In Fig. 10, we show model prices of call options on the QV computed using the second-order affine expansion under the MMA and inverse measures, and the comparison with MC estimates. We use the same expiries as for the model calibration and the same parameters reported in Eq. (6.4). We also observe that the solution under the both measures are consistent and agree with MC simulations. Interestingly, we observe that the implied volatility skew of options on the QV is upward sloping, which is observed in market data of options on the QV and listed options on VIX index. It is known that Heston model, even augmented with jumps, produces downward sloping skews on options on the QV (see for an example Drimus (2012) and Sepp (2012)). As a result, the log-normal SV model implies realistic patterns of implied volatility-of-volatility.

7. Conclusion

We have considered the log-normal SV model with the quadratic drift which plays an important role for the existence of martingale measures. We have shown that the log-normal volatility process has the strong solution in spite of having quadratic drift. We have derived admissible bounds of model parameters for the existence of both the MMA measure and the inverse measures. Based on that, we applied proposed SV model for modeling implied volatility surfaces of assets with positive return–volatility correlation.

Since the log-normal SV model is not affine, there is no analytical solution for the MGF in this model. To circumvent this, we have developed an affine expansion approach under which the joint MGF of three state variables (log-price, the QV, and the volatility) is decomposed into a leading term, which has exponential-affine form, and a residual term, whose estimate depends on the higher order moments of the volatility process. We have proved that the second-order leading term is consistent with the expected values and the variances of the state variables. We have applied Fourier inversion techniques for valuation of vanilla and inverse options on spot and futures underlying. By comparison of model values computed using our method with estimates obtained from MC simulations, we have shown that the second-order leading term is very accurate for option valuation. We have demonstrated that our model fits implied volatilities of assets with positive return–volatility correlation. We have shown that the log-normal SV model fits well different market regimes by calibrating model parameters to time series of options on Bitcoin for the past four years.

We have derived backward Euler–Maruyama scheme for discretization of the log-volatility process and show that it has a strong convergence rate of 1. Proposed scheme for the simulation of the log-normal SV model is robust, computationally efficient and does not require any special treatment in discretization.

As an extension of the framework^e we shall consider the applications to the rough SV models proposed by Gatheral et al. (2018) using log-normal volatility. We shall also consider a path-dependent model of the log-normal volatility following approaches in Guyon & Lekeufack (2023) and Lipton and Reghai (2023).^f

Acknowledgments

The authors are thankful to Vladimir Lucic for valuable insights and comments. We thank Alan Lewis, Johannes Muhle-Karbe, Blanka Horvath, Kalev Pärna and participants at Imperial College Finance and Stochastics Seminar, TU Munich Oberseminar, ETH Seminar, EAJ Conference 2022 in Tartu for useful discussion and comments. We are grateful for anonymous referees for helpful comments.

Notes

^a Github project https://github.com/ArturSepp/StochVolModels provides Python code with the implementation of option valuation using the affine expansion and examples of model calibration to implied volatility data applied in this paper.

^b Exp-OU SV models are studied in Fouque et al. (2000), Detemple & Osakwe (2000), Masoliver & Perelló (2006), Perelló et al. (2008), Muhle-Karbe et al. (2022), Drimus (2012). Exp-OU SV models are widespread in industry as a basis for local SV models (Lipton (2002)), see Bergomi (2015), Bühler (2006), Ren et al. (2007), Henry-Labordère (2009), Bergomi & Guyon (2012).

^c Such linear specification is the simplest functional form that allows to establish simple relationship between gains on a delta-hedged option portfolio and level of the market volatility. As shown in Bakshi & Kapadia (2003), such form allows to get precise information about the sign and strength of the risk premium.

^d We apply scipy implementation of Sequential Least Squares Programming (SLSQP) algorithm.

^e Sepp & Rakhmonov (2023a) apply the developed solution for modeling of SV in Cheyette-type interest rate model and obtain closed-form solution for valuation of swaptions under this model.

^f The volatility beta is an interesting variable to make path-dependent as well. We observe that in many markets the implied options skewness follows most recent price performance. For example, in options on Bitcoin, the skewness implied from option prices becomes positive following strong performance of Bitcoin.