Hyperbolic normal stochastic volatility model

Introduction

Stochastic volatility (SV) models have been proposed to overcome the failure of the Black-Scholes-Merton (BSM) model in explaining non-constant implied volatilities across strike prices on option markets, a phenomenon called volatility smile. Therefore, most previous studies (e.g., Hull and White 1987; Stein and Stein 1991; Heston 1993) discuss the SV models based on the geometric Brownian motion (BM) (hereafter, lognormal SV models). On the other hand, the studies on the SV models that are based on the arithmetic BM (hereafter, normal SV models) are scarce. This study aims to fill this gap by proposing and analyzing a class of normal SV models. Our motivation for choosing arithmetic BM as the backbone of the SV model is twofold: an options pricing model alternative to the lognormal SV model and a skewed and heavy-tailed distribution generated by a continuous-time stochastic process.

First, the study discusses the aspect of options pricing model. Although eclipsed by the success of the BSM model, the arithmetic BM is analyzed for the first time as an options pricing model by Bachelier 1900 (hereafter, normal model) and still provides more relevant dynamics than the geometric BM for some financial asset classes. Refer to Brooks and Brooks 2017 and Schachermayer and Teichmann 2008 for recent surveys on the normal model. An important difference between them is that the volatility under the normal model (hereafter, normal volatility) measures the uncertainty in terms of the absolute change in the asset price as opposed to relative change. One example of the applications of the normal model is its use for modeling the interest rate. The proportionality between the daily changes and level of interest rate—a key assumption of the BSM model—is empirically weak (Levin 2004). Therefore, among fixed-income market traders, the normal model has long been a popular alternative to the BSM model for quoting and risk-managing the options on interest rate swap and Treasury bonds (and futures). For example, the Merrill Lynch option volatility index (MOVE)—the bond market’s equivalent of the volatility index (VIX)—is calculated as the weighted average of the implied normal volatilities of the US Treasury bond options. It is also worth noting that the hedging ratio, delta, from the normal and BSM models can often be significantly different, even after the volatilities of the corresponding models are calibrated to the same option price observed on the market. Therefore, the normal model’s delta provides a more efficient hedge when the fluctuation of the underlying asset price is more consistent in absolute term than in percentage term. The use of the normal model for the interest market is further justified by the negative policy rates observed in several developed economies after the global financial crisis of 2008. Other than the interest rate, the normal model is often used for modeling the inflation rate (Kenyon 2008) and spread option (Poitras 1998).

Despite this background, it is difficult to find previous studies on the normal SV model. It is surprising, given that the lognormal SV models are often analyzed under the normal diffusion framework with the log price transformation; it implies that any existing results on the lognormal SV models can be effortlessly applied to the corresponding normal SV models. To the best of our knowledge, the only previous study on the normal SV model is in the context of the stochastic alpha-beta-rho (SABR) model (Hagan et al. 2002)—an SV model popular among practitioners. In the SABR model, the price follows a constant elasticity of variance (CEV) backbone, while the volatility follows a geometric BM. Therefore, the SABR model provides a range of backbone choices, including the normal and lognormal backbones. The SABR model with normal backbone (hereafter, normal SABR) is an important motivation for this study. A detailed review on the SABR model is provided in section 2.2.

2. Skewed and heavy-tailed distribution

The second motivation of our study is that the normal SV models can serve as a means to generate distributions with skewness and heavy-tail, generalizing the normal distribution. Heavy-tailed distributions are ubiquitous and their importance cannot be emphasized enough. In this regard, the study of normal SV models has a much broader significance than that of the lognormal SV models. This is because the latter generalizes the lognormal distribution whose application is limited when compared to the normal distribution.

Several distribution families have been proposed in statistics to incorporate skewness and heavy tails into a normal distribution. Even if the focus is narrowed to the applications to finance, it can be found that numerous distributions have been adopted to describe the statistics of asset return: generalized lambda (Corlu and Corlu 2015), stable (Fama 1965), skewed $t$ (Theodossiou 1998), Gaussian mixture (Kon 1984; Behr and Pötter 2009), generalized hyperbolic (Eberlein and Keller 1995; Behr and Pötter 2009), Turkey’s $g$- and $h$- (Badrinath and Chatterjee 1988; Mills 1995), and Johnson’s $S_{U}$ (hereafter $S_{U}$) distribution (Shang and Tadikamalla 2004; Gurrola 2007; Choi and Nam 2008).

However, the above distributions are neither defined from or associated with stochastic differential equations (SDEs), not to mention the SV models in particular. Those distributions are defined by the probability density function (PDF) or by the transformations of other well-known random variables. This is because it is usually difficult for an SDE to yield an analytically tractable solution. There are only a few examples of continuous-time processes whose transition probabilities correspond to the following well-known probability distributions: the arithmetic BM to a normal distribution (by definition), geometric BM to a lognormal distribution, and CEV and CIR processes to non-central $\chi^{2}$ distributions.

3. Contribution of this study

The study proposes and analyzes a class of normal SV models, which includes the normal SABR model as a special case. Since the mathematics behind our model involves the BMs in hyperbolic geometry and the results are expressed by hyperbolic functions, the class is named hyperbolic normal SV or NSVh The class is named as an abbreviation in a manner similar to the way hyperbolic sine becomes $\sinh$ model. The important mathematical tool to analyze the NSVh model comes from the two generalizations (Alili et al. 1997; Alili and Gruet 1997) of Bougerol’s identity (Bougerol 1983).

The first generalization leads us to a closed-form Monte-Carlo (MC) simulation scheme that no longer needs a time-discretized Euler scheme. The MC scheme requires merely one and a half (1.5) normal random numbers for a transition between time intervals of any length. Although limited to the normal SABR case, this study’s scheme is far more efficient than the previous exact MC scheme of Cai et al. 2017. Additionally, the original proof of the first generalization (Alili and Gruet 1997) is simplified in this study. The second generalization shows that a special case of the NSVh model— different from the normal SABR model—gives rise to the $S_{U}$ distribution (Johnson 1949), one of the popular heavy-tailed distributions. This allows the study to add to the literature one rare example of analytically tractable SDEs. The normal SV model provides a framework to understand the $S_{U}$ distribution in a better manner; the distribution can be parametrized more intuitively by using the NSVh parameters, and the popular use of the distribution is explained to some extent.

Importantly, under the NSVh model framework, two unrelated subjects are brought together, that is, the normal SABR model and the $S_{U}$ distribution. It is argued and empirically shown that the two distributions are very close to each other when parameters are estimated from the same data set, and can thus be used interchangeably. Among the benefits of the interchangeable usage of the two is the superior analytic tractability of the $S_{U}$ distribution when recognized as an options pricing model—various quantities of interest, such as the vanilla option price, density functions, skewness, ex-kurtosis, value-at-risk, and expected shortfall, have closed-form expressions that are not available in other SV models. To facilitate the interchangeability, a quick method of moments to convert the equivalent parameter sets between the two distributions is proposed.

This remainder of this paper is organized as follows. Section 2 defines the NSVh model and reviews the SABR model and $S_{U}$ distribution. Section 3 describes the main results. Section 4 presents the numerical results with empirical data. Finally, Section 5 concludes the paper.

Models and Preliminaries

where $F_{t}$ and $\sigma_{t}$ are the processes for the price and volatility, respectively, $\alpha$ is the volatility of the volatility parameter, $\rho$ denotes the instantaneous correlation between $F_{t}$ and $\sigma_{t}$, and $\rho_{\ast}=\sqrt{1-\rho^{2}}$. The BMs $Z_{t}$ and $X_{t}$ are independent, and $Z^{[\mu]}_{t}=Z_{t}+\mu\,t$ denotes BM with drift $\mu$.

The role of the model parameters $\rho,\alpha$, and $\lambda$ is discussed. Similar to the lognormal SV models, correlation $\rho$ accounts for the asymmetry in the distribution, that is, skewness or volatility skew. The leverage effect—the negative correlation between the spot price and volatility seen in the equity market—is achieved by a negative $\rho$, although it is in the context of the normal volatility in the NSVh model. The parameter $\alpha$ accounts for the heavy tail, that is, excess kurtosis or volatility smile. It can be easily seen that the process converges to an arithmetic BM in the limit $\alpha\rightarrow 0$ regardless of $\rho$. Therefore, $\alpha$ affects both the skewness and heavy tail at the same time.

The parameter $\lambda$ is present in the drift of $Z_{t}$ for both $F_{t}$ and $\sigma_{t}$. With regards to the volatility process, $\lambda$ controls the power of $\sigma_{t}$ that becomes a martingale as a geometric BM:

For example, $\lambda=0$ yields the volatility $\sigma_{t}$, following a driftless geometric BM, as in the SABR model, and $\lambda=-1$ yields the variance $\sigma_{t}^{2}$, following a driftless geometric BM as in the SV model of Hull and White 1987. With regards to the price process, however, the drift prevents $F_{t}$ from being a martingale except for $\lambda=0$ or, less importantly, $\rho=0$, although the expectation is easily computed as $\bar{F}_{t}=F_{0}+(\sigma_{0}\rho/\alpha)\big(e^{\lambda\alpha^{2}t/2}-1\big)$. Therefore, the resulting process may not be desirable as a price process. The NSVh model for $\lambda\neq 0$ is understood as a probability distribution perturbed from the $\lambda=0$ case, by applying the Radon-Nikodym derivative with respect to $Z_{t}$. Essentially, the introduction of $\lambda$ does not significantly diversify the shape of the distribution, and, therefore $\lambda$ is not meant for parameter estimation. As we shall see, however, $\lambda$ plays an important role in model selection; it brings under one unified process the three subjects separately studied: the normal SABR model ($\lambda=0$), Johnson’s $S_{U}$ distribution ($\lambda=1$), and BM on three-dimensional hyperbolic geometry ($\lambda=-1$).

To provide a background for the main result in section 3, we simplify the SDEs into the canonical forms,

where the following changes of variables are used:

The stochastic integrals of the canonical forms up to $s=S$ are, respectively, expressed as

This quantity has been the topic of extensive research; see Matsumoto and Yor 2005a; Matsumoto and Yor 2005b; Yor 2012 for a detailed review. While the functional is originally defined as the continuously averaged price under the BSM model in Asian options, it is used for the time-integral of the variance in the context of this study. Although $A^{[\mu]}_{T}$ can be defined with any standard BM, we implicitly assume that $A^{[\mu]}_{T}$ is tied to a particular BM, $Z_{t}$, throughout this study. Essentially, $A^{[\mu]}_{T}$ and $Z_{T}$ are closely intertwined, and the knowledge of their joint distribution of is the key to solve Equation (3).

2. SABR Model and Hyperbolic Geometry

The SABR model is reviewed with a focus on the normal backbone along with the BM on hyperbolic geometry, which serves as a mathematical tool for the NSVh and SABR models. The SABR model (Hagan et al. 2002) is an SV model with the backbone of the CEV process:

where $X_{t}$ and $Z_{t}$ are independent BMs. As mentioned earlier, the normal SABR model with $\beta=0$ is equivalent to the NSVh model with $\lambda=0$.

The SABR model has been widely used in the financial industry, for covering fixed income in particular, due to several merits: (i) arbitrary backbone choice, including normal ($\beta=0$) and lognormal ($\beta=1$) ones, (ii) availability of an approximate but fast vanilla options pricing method (Hagan et al. 2002), and (iii) parsimonious and intuitive parameters. The comments on those merits are presented in order. Regarding the CEV backbone, the popularity of the SABR model provides another evidence that the lognormal backbone of the BSM model is not a one-fits-all solution. The normal SABR in this study allows the negative value of $F_{t}$ without any boundary condition at zero. It should not be confused with the continuous limit of $\beta\rightarrow 0^{+}$, which does not allow a negative value. In the original article, Hagan et al. 2002 derives an approximate formula for the implied BSM volatility, from which the option price can be quickly computed through the BSM formula. However, it is worth noting that the normal volatility is first obtained from the small-time perturbation of the normal diffusion even for $\beta\neq 0$ and subsequently it is converted to the BSM volatility by another approximation (Hagan and Woodward 1999). Therefore, the option price computed from the normal volatility and the normal model formula has been considered more accurate because the second approximation can be avoided. Since the normal volatility is more appropriate for this study, the normal volatility approximation for $\beta=0$ is presented for reference and later use:

where $K$ is the strike price, and $T$ denotes the time-to-expiry. The volatility approximation is an asymptotic expansion that is valid when $\alpha^{2}T$ is small; therefore, the accuracy of the approximation noticeably deteriorates with an increase in $\alpha^{2}T$. Despite the shortcoming, the inaccurate approximation does not cause problems in pricing vanilla options because the model parameters $\sigma_{0},\alpha,\rho$ and the pre-determined $\beta$ are to be calibrated to the option prices observed from the market. In this regard, the implied volatility formula rather serves as an interpolation method for the volatility smile. Inaccurate approximation starts causing issues only when the usage of the model goes beyond vanilla options pricing. Two such cases are as follows:(i) claims with an exotic payout (e.g., quadratic) that require the knowledge of PDF, and (ii) path-dependent claims, which must resort to an MC simulation. In the first case, the PDF implied from Hagan et al. 2002’s formula often results in negative density at out-of-the-money strikes, thereby allowing arbitrage. In the second case, the vanilla option price from the formula is not consistent with that from the MC simulation with the same parameters. Therefore, the parameter calibration for MC scheme should be performed with extra care. Hence, the research on the accurate option analytics and efficient MC simulation methods come after the SABR model establishes its popularity among practitioners.

The study reviews prior works on the SABR model. Concerning vanilla options pricing, there have been various improvements to Hagan et al. 2002’s result. A few examples of such studies are Obłój 2007; Jordan and Tier 2011; Balland and Tran 2013; Lorig et al. 2015. However, they remain as approximations. The exact pricing is known only for the following three special cases: (i) zero correlation ($\rho=0$), (ii) lognormal SABR ($\beta=1$), and (iii) normal SABR ($\beta=0$). For the rest of the parameter ranges, no analytic solution is reported. Hence, the finite difference method (Park 2014; Le Floc’h and Kennedy 2017) is considered the most practical approach. Concerning the zero-correlation case, the price process can be transformed to the CEV process time-changed with $A^{[-1/2]}_{T}$ in a manner similar to that of Equation (3). Thus, the option price is expressed by a multi-dimensional integral representation of $A^{[-1/2]}_{T}$ over the CEV option prices (Schroder 1989). Refer to Antonov et al. 2013 for the most simplified expression based on the heat kernel on the two-dimensional hyperbolic geometry (McKean 1970), which we introduce below. The solution for the lognormal SABR is expressed in terms of the Gaussian hypergeometric series (Lewis 2000).

The development of MC simulation methods of the SABR dynamics is relatively recent. While several efficient approximations (Chen et al. 2012; Leitao et al. 2017b; Leitao et al. 2017a) have been proposed, an exact simulation method Cai et al. 2017 is available for the following three special cases: (i) $\rho=0$, (ii) $\beta=0$, and (iii) $\beta=1$. The key element in this method is to simulate the time-integrated variance $A^{[-1/2]}_{S}$, conditional on the terminal volatility $Z_{S}$. The cumulative distribution function (CDF) of the quantity is obtained from the Laplace transform of $(1/A^{[-1/2]}_{S})\,|\,Z_{S}$, which has a closed-form expression (Matsumoto and Yor 2005a). Given the exact random numbers of $Z_{S}$ and $A^{[-1/2]}_{S}$, $X_{A^{[-1/2]}_{S}}$ in the normal SABR model is easily simulated as $X_{1}\sqrt{A^{[-1/2]}_{S}}$ for a standard normal variable $X_{1}$. Although a heavy Euler scheme is avoided, the method of Cai et al. 2017 still incurs a moderate computation cost due to the numerical inversion of the Laplace transform and root-solving for the CDF inversion.

3. Johnson’s Distribution Family

Johnson 1949 proposes a system of distribution families in which a random variable $X$ is represented by the transformations from a standard normal variable $Z$:

where $\gamma_{X}$ and $\gamma_{Z}$ are location parameters and $\delta_{X}$ and $\delta_{Z}$ are scaling parameters. Although not explicitly included, normal distribution can be considered as a special intersection of the three families in the limit of $\delta_{X}$ and $\delta_{Z}$ proportionally going to infinity. Therefore, it is often included as $S_{N}$ (normal) family with $f(x)=x$. The $X$ range is unbounded for $S_{U}$ and $S_{N}$, semi-bounded for $S_{L}$, and bounded for $S_{B}$. The system is designed in such a way that a unique family is chosen for any mathematically feasible pair of skewness and kurtosis. For a fixed value of skewness, the kurtosis increases in the order of $S_{B}$, $S_{L}$, and $S_{U}$.

Particularly, the $S_{U}$ family has been an attractive choice for modeling a heavy-tailed data set, and has been adopted in various fields; refer to Jones 2014 and the references therein. Examples in finance includes heavy-tailed innovation in the GARCH model (Choi and Nam 2008), prediction of value-at-risk (Simonato 2011; Venkataraman and Rao 2016), and asset return distribution (Shang and Tadikamalla 2004; Corlu and Corlu 2015).

The $S_{U}$ distribution has several advantages over alternative heavy-tailed distributions. First, it explains a wide range of skewness and kurtosis. For a fixed value of skewness, it can accommodate arbitrary high values of kurtosis, which is not feasible in the classical approaches to generalize a normal distribution, such as the Gram-Charlier or Cornish-Fisher expansions. Second, many properties of the distributions are available in closed forms: PDF, CDF, skewness, and kurtosis. Third, the parameters are efficiently estimated—refer to Tuenter 2001 for the moment matching in the reduced form and Wheeler 1980 for the quantile-based estimation. Finally, drawing random numbers is easy, which makes $S_{U}$ distribution ideal for MC simulations, particularly in a multivariate setting (Biller and Ghosh 2006). In general, random number sampling is not trivial, even if the distribution functions are given in closed forms.

In addition to the existing merits, our result in section 3.2 gives a first-class-citizen status to the $S_{U}$ distribution among other heavy-tailed distributions by showing that it is a solution of a continuous-time SV process, the NSVh model with $\lambda=1$. This partially explains why the $S_{U}$ distribution has been superior in modeling asset return distributions and risk metrics.

Main Results

The study’s main results first present Bougerol 1983’s identity in the original form. Since the original identity is generalized later, we state it as a Corollary and defer the proof to Proposition 2.

For a fixed time $T$, the following is equal in distribution:

where $X_{t}$, $W_{t}$, and $Z_{t}$ are independent BMs, and $A_{T}$ is defined by Equation (4).

This identity is surprising in that the stochastic integral involving two independent BMs is equal in distribution to the $\sinh$ transformation of one BM. Refer to Matsumoto and Yor 2005a; Vakeroudis 2012 for a review and related topics. The identity should be interpreted with caution; the equality holds as distribution ($\,{\mathrel{\mathop{\kern 0.0pt=}\limits^{d}}}\,$) at a fixed time $t=T$, not as a process for $0\leq t\leq T$. Moreover, it does not directly help to solve Equation (3). The identity must be generalized to non-zero drift, $A^{[\mu]}$, and provide the joint distribution with $Z^{[\mu]}$, which is found in Alili and Gruet 1997 and Alili et al. 1997. In the following subsections, we apply two generalizations to the NSVh model; one to the general $\lambda$, and the other to a special case $\lambda=1$.

Let $X_{t}$ and $Z_{t}$ be two independent BMs and the function $\phi$ defined by

then the following is equal in distribution, conditional on $Z^{[\mu]}_{T}$:

where $R_{t}$ is a two-dimensional Bessel process, that is, the radius of a BM in two-dimensional Euclidean geometry, and $\theta$ is a uniformly distributed random angle. The three random variables $R_{T}$, $Z_{T}$, and $\theta$ are independent.

The joint distribution of the NSVh model at a fixed time $S$ is given as

Furthermore, the three independent random variables can be simulated as

where $X_{1}$ and $Y_{1}$ are independent standard normals.

2. SUS_{U} Distributions for λ=1\lambda=1

This subsection shows that the NSVh distribution for $\lambda=1$ is expressed by the $S_{U}$ distribution and is related to Bougerol’s identity generalized to an arbitrary starting point. In the following proposition, Proposition 4 of Alili and Gruet 1997 (or Theorem 3.1 of Matsumoto and Yor 2005a) is restated. More general results are found in Proposition 1 of Alili et al. 1997 (or Proposition 2.1 of Vakeroudis 2012), and we follow the proof therein.

For a fixed time $T$ and independent BMs, $X_{t}$, $Z_{t}$, and $W_{t}$, the following is equal in distribution:

are equivalent because they start from the same starting point $P_{0}=Q_{0}=\sinh(a)$ and follow the SDE:

Therefore, $P_{t}$ and $Q_{t}$ have the same distribution for any time $t$. The equality between $Q_{T}$ and the left-most expression is shown by the time-reversal $s\rightarrow T-s$. For a fixed time $T$, $Z_{T}-Z_{T-s}$ for $0\leq s\leq T$ is also a standard BM with the same ending point $Z_{T}$, and therefore it may be replaced with $Z_{s}$. ∎

The original Bougerol’s identity in Corollary 1 is a special case, with $a=0$. Now, Proposition 2 can be applied to further simplify the NSVh distribution for $\lambda=1$.

The price of the NSVh model with $\lambda=1$ at a fixed time $S$ follows a re-parametrized $S_{U}$ distribution:

where the original parameters are mapped by

The results are easily proved from the following hyperbolic function identities,

Conversely, from the $S_{U}$ distribution, the NSVh model obtains analytic tractability for option price and risk measures. Below are the closed-form expressions for the quantities of interest:

For an asset price following the NSVh process with $\lambda=1$, option price, value-at-risk, and expected shortfall have the following closed-form solutions.

The undiscounted price of a vanilla option with strike price $x$:

where $d$ is defined in Equation (16) and $\pm$ indicates call/put options, respectively.

The option prices are easily derived by integrating Equation (15) with the boundary $d$ obtained from Equation (14). The relationship between the put option value and the expected shortfall, $\text{ES}(p)=\text{VaR}(p)-V_{-}(\text{VaR}(p))/p$, is useful, where $\text{VaR}(p)$ vanishes in the final expression for $\text{ES}(p)$. ∎

Later, it is argued that the NSVh distributions with different values of $\lambda$ are close to each other, and thus the analytically tractable $\lambda=1$ case can represent the rest including the normal SABR model ($\lambda=0$). The option price from the closed-form formula serves as a benchmark against which the option price from the MC scheme of Corollary 2 is compared in Section 4.

3. Moments Matching of the NSVh Distribution

The study derives the moments of the NSVh distribution for general $\lambda$ to be used for parameter estimation. The study also proposes a moment matching in the reduced form for $\lambda=0$ to complement that for $\lambda=1$ by Tuenter 2001.

Refer to Appendix B for detailed derivation. Corollary 5 generalizes the moments for $S_{L}$ ($\rho=\pm 1$) and $S_{U}$ ($\lambda=1$) distributions.

Parallel to Tuenter 2001’s reduced moment matching method for the $S_{U}$ distribution, a similar method is developed for the normal SABR model. Combined with Tuenter 2001, the two methods can quickly find an equivalent parameter set of one distribution from the other. By joining the expressions for $s$ and $\kappa$ in Equation (21) through $\rho$, $\kappa$ is expressed as a univariate function of $w\geq 1$:

for which the study numerically finds the root $w_{*}$ of $\kappa=f(w_{*})$. It can be shown that $f(w)$ is monotonically increasing for $w\geq 1$. Therefore, the root $w_{*}$ would be unique if it exists. We can further bound $w_{*}$ by $w_{m}\leq w_{*}\leq w_{M}$ to expedite the numerical root-finding. Lower bound $w_{m}$ is the unique cubic root of $s^{2}=(w-1)(w+2)^{2}$ (the $\rho=\pm 1$ case) for $w\geq 1$:

Upper bound $w_{M}$ is obtained by plugging $w=w_{m}$ into Equation (22), except the $(w-1)$ term:

The existence of $w_{*}$ is equivalent to $f(w_{m})\leq\kappa$. If $w_{*}$ exists and is found from the numerical root-finding, then the parameters can be solved as

4. Summary of results

In Figure 2, the relationship of the NSVh model and other related models is shown. In Table 2, the results for the three important drift values, $\lambda=-1$, 0, and 1, are summarized for comparison.

Parameter Estimation from Empirical Data

The NSVh distribution is calibrated to two empirical data sets— swaption volatility smile and daily stock index return. The purpose of this exercise is to demonstrate various numerical procedures presented in this study rather than arguing that the NSVh model is superior to other SV models or heavy-tailed distributions in fitting these data. Additionally, the study shows that the two NSVh models, that is, $\lambda=0$ (normal SABR) and $\lambda=1$ ($S_{U}$), yield very similar distributions, and thus can be used interchangeably, if calibrated to the same target, such as implied volatility or moments.

The study obtains the US swaption market prices on March 14, 2017 from Reuters. The two heavily traded expiry–tenor pairs of the US swaption—1y1y and 10y10y—are chosen to illustrate different volatility smile shapes. To avoid the complication of the annuity price of the underlying swap, the study computes the price in the unit of the annuity from the BSM implied volatilities provided by Reuters, rather than raw dollar prices. Refer to the insets of Figure 3 for the BSM implied volatilities.

Figure 3 shows the normal volatility smile implied from the market and the calibrated models. While option prices are observable for strike prices with spreads of 0, $\pm 0.25\%$, $\pm 0.5\%$, $\pm 1.0\%$, and $\pm 1.5\%$ from the forward swap rates $\bar{F}_{T}$, the study uses only the three spreads—0 and $\pm 1\%$ for calibration—to ensure that the calibrated parameter set—$(\sigma_{0},\alpha,\rho)$—reproduces option prices at those three strike prices. For calibration, the study uses Equations (6) for $\lambda=0$ and (17) for $\lambda=1$. The volatility smile curves implied from the two models are indistinguishable, and thus close to each other in distributions. Table 3 shows the calibrated parameter values for $\lambda=0$ and 1 for reference.

In Table 4, the option prices are compared from the MC simulation of Corollary 2 with the option prices from the aforementioned analytic methods. In the experiment, the MC simulation of $10^{6}$ samples is repeated 100 times. For $\lambda=1$, both Equation (17) and the MC method are exact, and thus the option prices from the two methods have very little difference due to the MC noise. For $\lambda=0$, however, Equation (6) is an approximation. The analytic prices show clear deviation from the accurate MC prices. Overall, this exercise reconfirms the possibility of the $S_{U}$ distribution being used as a better alternative to the normal SABR model.

2. Daily Return of Stock Index

We fit the NSVh distribution to the daily returns of two stock indices—the US Standard & Poor’s 500 Index (S&P 500) and China Securities Index 300 (CSI 300). The data covers the 12-year period from the beginning of 2005 to the end of 2016. For the analysis to be easily reproducible, daily returns are computed as the return on outright index values, rather than as holding period return.

The statistics of the daily returns are summarized in Table 5. S&P 500 shows heavier tails but less skewness than CSI 300. The table also shows the fitted parameters of the NSVh distributions for $\lambda=0$ and 1. In fitting, the reduced moment-matching methods of Tuenter 2001 and section 3.3 are used. Based on these values, the value-at-risk and expected shortfall are computed and compared to those from the normal distribution assumption and the historical data. While Equations (18) and (19) are used for $\lambda=1$, MC simulation is used to compute the risk measures for $\lambda=0$. The results are shown in Table 6. As the return distribution has heavy tail, the value-at-risk and expected shortfall from the NSVh distributions are much closer to those from the historical data than those from the normal distribution. The risk measures from $\lambda=0$ and 1 are very close to each other, confirming the similarity between the two distributions.

Finally, the goodness-of-fit of the $S_{U}$ ($\lambda=1$) is illustrated by using the probability plot. In Figure 4, the theoretical $Z$-score—$Z_{0}^{(j)}=N^{-1}((j-1/2)/n)$ for the $j$-th ordered sample—is shown on the $x$-axis and two different sample’s $Z$-scores on the $y$-axis as follows: (i) from the estimated normal distribution $Z_{1}^{(j)}=(X_{j}-\bar{F}_{T})/\sqrt{\mu_{2}}$ and (ii) from the $S_{U}$ distribution computed as

Therefore, $(Z_{0},Z_{2})$ is understood as the $S_{U}$ probability plot, while $(Z_{0},Z_{1})$ is the normal probability plot in the usual definition. Figure 4 shows that the points under the $S_{U}$ probability plot are close to the $y=x$ line, indicating that the stock returns closely follow the $S_{U}$ distribution ($\lambda=1$).

Conclusion

The study generalizes the arithmetic Brownian motion with stochastic volatility. The NSVh process proposed in this study incorporates the SABR model and Johnson’s $S_{U}$ distribution, which have been studied in different contexts. The SABR model is a well-known option pricing model in financial engineering, and the Johnson’s $S_{U}$ distribution is a popular skewed and heavy-tailed distribution defined by a transformation from the normal random variable. From the generalizations of Bougerol’s identity, the NSVh model is equipped with closed-form MC simulation for general cases and closed-form option formula for the $S_{U}$ case. This study demonstrates the usage of the model with two empirical datasets—the US swaption and daily returns distribution of the S&P 500 and CSI 300 indexes.

Acknowledgements

The authors are grateful to Larbi Alili for sharing manuscripts, Alili and Gruet 1997 and Alili et al. 1997, Robert Webb (editor), and Minsuk Kwak (discussant at the 2018 Asia-Pacific Association of Derivatives conference in Busan, Korea). Jaehyuk Choi would like to express gratitude to his previous employer, Goldman Sachs, for laying the foundation for the motivation of this study. Byoung Ki Seo was supported by the Institute for Information & communications Technology Promotion (IITP) grant funded by the Ministry of Science and ICT (MSIT), Korea (No. 2017-0-01779, A machine learning and statistical inference framework for explainable artificial intelligence).

References

Appendix A Proof of Proposition 1

. If we let $D_{t}$ be the hyperbolic distance between $(x_{t},y_{t},z_{t})$ and the starting point $(0,0,1)$,then

the Euclidean radius of $(x_{t},y_{t})$ is expressed by the function $\phi$:

The underlying BM $Z^{[\mu]}_{t}$ can be also interpreted as the projection of $(x_{t},y_{t},z_{t})$ on the $z$-axis, that is, the signed hyperbolic distance from $(0,0,1)$ to $(0,0,z_{t})$. Therefore, the restriction $Z^{[\mu]}_{t}\leq D_{t}$ is naturally satisfied.

A critical step of the proof is to show, for a fixed time $T$,

, which effectively means that the hyperbolic distance between $(x_{T},y_{T},z_{T})$ and the starting point $(0,0,1)$ has the same distribution as the Euclidean distance of the underlying BMs, $(X_{T},Y_{T},Z^{[\mu]}_{T})$, from $(0,0,0)$. Furthermore, the identity holds, conditional on $Z^{[\mu]}_{T}$. Based on the identity, it follows that

where $\theta$ is a uniformly distributed random angle.

From $dz_{T}/z_{T}=dZ^{}_{T}$, it is also seen that

Therefore, the conditional probability is given as

The probability can be interpreted as the conditional probability $\text{Prob}(\sqrt{X_{T}^{2}+Y_{T}^{2}+(Z^{}_{T})^{2}}\in dD_{T}\;|\;Z_{T})$. ∎

Appendix B Derivation of moments of NSVh distribution

We first compute the moments, conditional on $Z^{[\mu]}_{T}$:

where $u=|Z^{[\mu]}_{T}|/\sqrt{T}$ and $(2n-1)!!=(2n-1)(2n-3)\cdots 3\cdot 1$. This formula is comparable to the formula for $E\Big(\big({A^{[\mu]}_{T}}\big)^{n}\Big|Z^{[\mu]}_{T}\Big)$ given in Proposition 5.3 in Matsumoto and Yor 2005a. The first two values are computed in closed form:

From these, the first two conditional moments of $A^{[\mu]}_{T}$ are trivially obtained as

The same results for $\lambda=0$ are derived in Kennedy et al. 2012 in the context of the normal SABR model.

The unconditional moments of $X_{A^{[\mu]}_{S}}$ are given as

where $w=e^{S}$. Subsequently, the expressions for the moments follow from

and the substitution $\mu=(\lambda-1)/2$. ∎