Hyperbolic normal stochastic volatility model

For option pricing models and heavy-tailed distributions, this study proposes a continuous-time stochastic volatility model based on an arithmetic Brownian motion: a one-parameter extension of the normal stochastic alpha-beta-rho (SABR) model. Using two generalized Bougerol's identities in the literature, the study shows that our model has a closed-form Monte-Carlo simulation scheme and that the transition probability for one special case follows Johnson's S_U distribution---a popular heavy-tailed distribution originally proposed without stochastic process. It is argued that the S_U distribution serves as an analytically superior alternative to the normal SABR model because the two distributions are empirically similar.



There are no comments yet.


page 1

page 2

page 3

page 4


An extension of Heston's SV model to Stochastic Interest Rates

In 'A Closed-Form Solution for Options with Stochastic Volatility with A...

Vector autoregression models with skewness and heavy tails

With uncertain changes of the economic environment, macroeconomic downtu...

A Novel Maneuvering Target Tracking Approach by Stochastic Volatility GARCH Model

In this paper, we introduce a new single model maneuvering target tracki...

Comparative Efficiency of Altruism and Egoism as Voting Strategies in Stochastic Environment

In this paper, we study the efficiency of egoistic and altruistic strate...

Simulation using random numbers

This article is devoted to methods of construction and study of stochast...

On the Distribution of Traffic Volumes in the Internet and its Implications

Getting good statistical models of traffic on network links is a well-kn...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. 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.

1.1. Options pricing model

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 


1.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 (Theodossiou, 1998), Gaussian mixture (Kon, 1984; Behr and Pötter, 2009), generalized hyperbolic (Eberlein and Keller, 1995; Behr and Pötter, 2009), Turkey’s - and - (Badrinath and Chatterjee, 1988; Mills, 1995), and Johnson’s (hereafter ) 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


1.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 111The class is named as an abbreviation in a manner similar to the way hyperbolic sine becomes 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 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 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

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

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 distribution. Section 3 describes the main results. Section 4 presents the numerical results with empirical data. Finally, Section 5 concludes the paper.

2. Models and Preliminaries

2.1. NSVh Model

The NSVh model is introduced as


where and are the processes for the price and volatility, respectively, is the volatility of the volatility parameter, denotes the instantaneous correlation between and , and . The BMs and are independent, and denotes BM with drift .

The role of the model parameters , and is discussed. Similar to the lognormal SV models, correlation 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 , although it is in the context of the normal volatility in the NSVh model. The parameter 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 regardless of . Therefore, affects both the skewness and heavy tail at the same time.

The parameter is present in the drift of for both and . With regards to the volatility process, controls the power of that becomes a martingale as a geometric BM:

For example, yields the volatility , following a driftless geometric BM, as in the SABR model, and yields the variance , 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 from being a martingale except for or, less importantly, , although the expectation is easily computed as . Therefore, the resulting process may not be desirable as a price process. The NSVh model for is understood as a probability distribution perturbed from the case, by applying the Radon-Nikodym derivative with respect to . Essentially, the introduction of does not significantly diversify the shape of the distribution, and, therefore is not meant for parameter estimation. As we shall see, however, plays an important role in model selection; it brings under one unified process the three subjects separately studied: the normal SABR model (), Johnson’s distribution (), and BM on three-dimensional hyperbolic geometry ().

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:

Here, the new variable is the integrated variance of the log volatility, and are the non-dimensionalized price and volatility processes, respectively, under . The scaling of and naturally follows from the time change of the BMs with the new variable . The time is any fixed time of interest such as the time-to-expiry of the vanilla option. The price is shifted by first to ensure that at the corresponding time . Therefore, the canonical NSVh distribution is effectively parametrized by the three parameters—. Additionally, the original distribution is recovered by and . While the original and canonical representations are explicitly distinguished, the variable is often used in the original form as well for the sake of concise notation.

The stochastic integrals of the canonical forms up to are, respectively, expressed as


It must be noted that the integral of , with respect to , is further simplified to the BM time-changed with an exponential functional of the BM defined by


This quantity has been the topic of extensive research; see Matsumoto and Yor (2005a, b); 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 can be defined with any standard BM, we implicitly assume that is tied to a particular BM, , throughout this study. Essentially, and

are closely intertwined, and the knowledge of their joint distribution of is the key to solve Equation (


2.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 and are independent BMs. As mentioned earlier, the normal SABR model with is equivalent to the NSVh model with .

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 () and lognormal () 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 without any boundary condition at zero. It should not be confused with the continuous limit of , 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 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 is presented for reference and later use:


where is the strike price, and denotes the time-to-expiry. The volatility approximation is an asymptotic expansion that is valid when is small; therefore, the accuracy of the approximation noticeably deteriorates with an increase in . Despite the shortcoming, the inaccurate approximation does not cause problems in pricing vanilla options because the model parameters and the pre-determined

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 (), (ii) lognormal SABR (), and (iii) normal SABR (). 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 in a manner similar to that of Equation (3). Thus, the option price is expressed by a multi-dimensional integral representation of 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 options pricing under the normal SABR model depends on the mathematical tools developed for the BMs on hyperbolic geometry, represented as Poincaré half-plane. Since the NSVh model also benefits from the same tools, we briefly introduce these tools. The -dimensional Poincaré half-plane is denoted by . Table 1 provides a quick reference for the properties of and . The standard BM in a geometry is defined to be the stochastic process whose infinitesimal generator is given by the Laplace-Beltrami operator of the geometry. The heat kernel is the fundamental solution of the diffusion equation , and hence the transition probability of the standard BM. Table 1 shows the heat kernels on  (McKean, 1970), often referred to as the McKean kernel, and (Debiard et al., 1976). The analytical expressions for the heat kernels are also known for , in general; see Grigor’yan and Noguchi (1998) for derivation.

Volume element
Geodesic distance
Standard BM
Heat kernel
for or
Table 1. Properties of the -Dimensional Hyperbolic Geometry Represented by -Dimensional Poincaré Half-Plane for and . Symbols and are the shortened notations for partial derivative operators and , respectively.

The standard BM on is equivalent to the normal SABR with in canonical form, where the -axis is for the price process and the -axis for the volatility process. Naturally, the heat kernel has been used for the analysis of the normal SABR. Henry-Labordère (2005, 2008) express the vanilla options price under the normal SABR model with a two-dimensional integral, although it is later corrected by Korn and Tang (2013). Antonov et al. (2015) further simplify the price to a one-dimensional integration with an approximation. However, in the absence of efficient numerical schemes to evaluate those integral representations, the normal volatility approximation in Equation (6) remains as a practical approach to price vanilla options under the normal SABR model.

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, a) have been proposed, an exact simulation method Cai et al. (2017) is available for the following three special cases: (i) , (ii) , and (iii) . The key element in this method is to simulate the time-integrated variance , conditional on the terminal volatility

. The cumulative distribution function (CDF) of the quantity is obtained from the Laplace transform of

, which has a closed-form expression (Matsumoto and Yor, 2005a). Given the exact random numbers of and , in the normal SABR model is easily simulated as for a standard normal variable . 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.

As opposed to the previous studies using , this study’s main result in section 3.1, based on Alili and Gruet (1997), exploits . The pairs and from the standard BM in are the two NSVh processes with and , and it should be noted that the introduction of to the NSVh model makes this connection possible. Despite a higher dimensionality, the heat kernel of is given without an integral, and Alili and Gruet (1997) shows that the squared radius and can be exactly simulated. Subsequently, (or ) is extracted via cosine (or sine) projection with random angle. Therefore, the study’s simulation scheme can be considered as a hyperbolic-geometry extension of the Box-Muller algorithm (Box and Muller, 1958) for generating normal random variable, in which the -axis of is additionally added.

2.3. Johnson’s Distribution Family

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


where and are location parameters and and are scaling parameters. Although not explicitly included, normal distribution can be considered as a special intersection of the three families in the limit of and proportionally going to infinity. Therefore, it is often included as (normal) family with . The range is unbounded for and , semi-bounded for , and bounded for . 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 , , and .

Particularly, the 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 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

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 distribution among other heavy-tailed distributions by showing that it is a solution of a continuous-time SV process, the NSVh model with . This partially explains why the distribution has been superior in modeling asset return distributions and risk metrics.

3. 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 1.

Corollary 1 (Bougerol’s identity).

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


where , , and are independent BMs, and is defined by Equation (4).

This identity is surprising in that the stochastic integral involving two independent BMs is equal in distribution to the 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 () at a fixed time , not as a process for . Moreover, it does not directly help to solve Equation (3). The identity must be generalized to non-zero drift, , and provide the joint distribution with , 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 , and the other to a special case .

3.1. Monte-Carlo Simulation Scheme

The first generalization makes use of the BMs in to obtain the distribution of , conditional on . We restate Proposition 3 in Alili and Gruet (1997) in a modified and stronger form: [Bougerol’s identity in hyperbolic geometry]propboughyp Let and be two independent BMs and the function defined by


then the following is equal in distribution, conditional on :


where is a two-dimensional Bessel process, that is, the radius of a BM in two-dimensional Euclidean geometry, and

is a uniformly distributed random angle. The three random variables

, , and are independent. Refer to Appendix A for the proof. Among the two original proofs in Alili and Gruet (1997), the proof exploiting the interpretation with is provided. The proposition in this study is stronger than the original because the distribution equality holds, conditional on , which is implied from the proof without difficulty. The study further simplifies the original proof by using the heat kernel. Proposition 3.1 can be directly applied to Equation (3) to obtain the transition equation and MC scheme:

Corollary 2.

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


Furthermore, the three independent random variables can be simulated as


where and are independent standard normals.

The simulation method using standard normal variables in Equation (12) is more efficient than drawing and independently because the costly evaluation is avoided. The idea is similar to the Marsaglia polar method (Marsaglia and Bray, 1964) for drawing normal random variables. It must be noted that three random numbers—, , and —generate two pairs of and . Therefore, one draw only requires one and a half (1.5) normal random variables, which is an unprecedented efficiency for any SV model simulation. Particularly, this method is more efficient than the exact SABR simulation of Cai et al. (2017), although it is limited to the normal case. This study’s method directly draws and , whereas Cai et al. (2017) first draws and , and subsequently . Although Equation (11) states the transition from to , it can handle any time interval from to (). Therefore, the scheme is ideal for pricing path-dependent claims.

3.2. Distributions for

This subsection shows that the NSVh distribution for is expressed by the 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.

Proposition 1 (Bougerol’s identity with an arbitrary starting point).

For a fixed time and independent BMs, , , and , the following is equal in distribution:


The two processes,

are equivalent because they start from the same starting point and follow the SDE:

Therefore, and have the same distribution for any time . The equality between and the left-most expression is shown by the time-reversal . For a fixed time , for is also a standard BM with the same ending point , and therefore it may be replaced with . ∎

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

Corollary 3.

The price of the NSVh model with at a fixed time follows a re-parametrized distribution:


where the original parameters are mapped by

It also admits a simpler form:


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

Although Proposition 1 is a well-known result, to the best of our knowledge, this is the first time that it is interpreted in the context of the distribution or SV model. Compared to Corollary 2, Corollary 3 is an even more efficient MC scheme requiring one normal random number for one draw of price, although for a special case . However, this is at the expense of the terminal volatility being lost. Unlike Corollary 2, Corollary 3 can only generate the final price , and therefore cannot be used for path-dependent claims.

The key consequence of Corollary 3 is that the NSVh model bridges the SDE-based SV model and heavy-tailed distribution. The encounter between the two subjects is mutually beneficial. First, as a solution of the NSVh process, the distribution obtains a more intuitive parametrization than the original in Equation (7). From Equation (15), it is clear that the symmetric heavy tail comes from the term, controlled by , and the asymmetric skewness from the term, controlled by . The new parametrization also helps to understand the relationship between Johnson family members. The lognormal family is recognized as a special case with () as . The normal family is obtained as in the limit of . Under the NSVh parameters, the well-known PDF and CDF of the distribution are respectively expressed by


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

Corollary 4.

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

  • The undiscounted price of a vanilla option with strike price :


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

  • Value-at-risk for quantile :

  • Expected shortfall for quantile :


The option prices are easily derived by integrating Equation (15) with the boundary obtained from Equation (14). The relationship between the put option value and the expected shortfall, , is useful, where vanishes in the final expression for . ∎

Later, it is argued that the NSVh distributions with different values of are close to each other, and thus the analytically tractable case can represent the rest including the normal SABR model (). 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.3. Moments Matching of the NSVh Distribution

The study derives the moments of the NSVh distribution for general to be used for parameter estimation. The study also proposes a moment matching in the reduced form for to complement that for by Tuenter (2001).

coronsvhmoments The central moments of the canonical NSVh distribution , for , are given as


The central moments of the original form can be scaled as , and the skewness and ex-kurtosis are given as and , respectively. For the normal SABR (), further simplified expressions are obtained:


Refer to Appendix B for detailed derivation. Corollary 3.3 generalizes the moments for () and () distributions.

The similarity of the NSVh distributions for different is inferred from skewness and ex-kurtosis. The result for large , at least for and 4 implies that the leading order of skewness and ex-kurtosis are independent of as and , as indicated in Equation (21). To illustrate that, in Figure 1, the contours of skewness and ex-kurtosis for and 1, as functions of (ex-kurtosis) and (skewness) for , are plotted. Although the parameters for is slightly higher than those for to obtain the same skewness and ex-kurtosis levels, the contours for and are very similar implying the similarity between the normal SABR and the distributions.

Figure 1. Contour Plot of Skewness (Red Dashed Line) and Excess Kurtosis (Blue Solid Line) for Varying versus . The upper left triangle is for () and the lower right triangle for (normal SABR). The values for skewness are 0, 1.5, 3, 4.5, 6, and 8, and the values for excess kurtosis are 2, 7, 16, 40, 100, and 200 from the lower left to the upper right corner.

Parallel to Tuenter (2001)’s reduced moment matching method for the 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 and in Equation (21) through , is expressed as a univariate function of :


for which the study numerically finds the root of . It can be shown that is monotonically increasing for . Therefore, the root would be unique if it exists. We can further bound by to expedite the numerical root-finding. Lower bound is the unique cubic root of (the case) for :

Upper bound is obtained by plugging into Equation (22), except the term:

The existence of is equivalent to . If exists and is found from the numerical root-finding, then the parameters can be solved as

3.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, , 0, and 1, are summarized for comparison.

Figure 2. The Overview of the NSVh Model in Relation to Other Previously Known Models and Distributions.
Drift parameter
Equivalent Standard BM in Normal SABR distribution
model/distribution Standard BM in
Drifted BM for
Terminal volatility
Integrated variance
Exact MC simulation Corollary 2 Corollary 2 & 3
Vanilla option price Equation (6) Equation (17)
(approximation) (exact)
Moments Corollary 3.3
Moment matching Equation (22) Tuenter (2001)
Table 2. Summary of the NSVh Model for the Three Key Special Cases: , and .

4. 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, (normal SABR) and (), yield very similar distributions, and thus can be used interchangeably, if calibrated to the same target, such as implied volatility or moments.

4.1. Swaption Volatility Smile

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. Swaption Volatility Smile in the Implied Normal Volatility (in annual basis points) Observed in the US Market on March 14, 2017: (a) 1 year into 1-year swap (1y1y) and (b) 10 years into 10-year swap (10y10y). The circles represent the volatilities implied from observed market prices, among which the black ones are ATM and ATM 1% points for calibration. The solid (blue) line represents the model-implied smile curve for () and the dashed (red) one for (normal SABR). The insets are the implied BSM volatilities (in annual %).

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, , , , and from the forward swap rates , the study uses only the three spreads—0 and for calibration—to ensure that the calibrated parameter set——reproduces option prices at those three strike prices. For calibration, the study uses Equations (6) for and (17) for . 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 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 samples is repeated 100 times. For , 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 , 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 distribution being used as a better alternative to the normal SABR model.

Calibrated 1y1y () 10y10y ()
33.503 32.244 1.697 1.580
61.962 62.181 22.372 22.196
0.533 0.477 0.691 0.609
2.0221 3.0673
Table 3. Parameters Calibrated to the US Swaption Volatility Smile on March 14, 2017. Refer to Figure 3 for the calibration points. Equations (6) and (17) are used for the NSVh models with and , respectively.
-200 2.275E-2 -4.1E-5 1.8E-5 2.274E-2 6.4E-7 1.8E-5
-100 1.506E-2 -2.1E-5 1.6E-5 1.506E-2 6.0E-7 1.6E-5
0 9.083E-3 -1.2E-5 1.3E-5 9.083E-3 5.3E-7 1.3E-5
100 5.108E-3 -2.6E-5 1.1E-5 5.108E-3 3.3E-7 1.1E-5
200 2.807E-3 -4.8E-5 8.8E-6 2.804E-3 3.9E-7 9.1E-6
300 1.567E-3 -6.0E-5 7.0E-6 1.559E-3 3.7E-7 7.3E-6
Table 4. Vanilla Options Pricing Tested Against the Parameters of Table 3 for 10y10y Swaption. Analytic prices () are computed from Equations (6) and (17) and MC price () is shown relative to

with the standard deviation. The MC simulation with

paths is repeated 100 times. Prices are in the unit of the annuity of the underlying swap.

4.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 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 , MC simulation is used to compute the risk measures for . 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 and 1 are very close to each other, confirming the similarity between the two distributions.

Summary statistics S&P 500 CSI 300
3020 2914
0.0282 0.0417
1.5154 3.4092
-0.0933 -0.5075
11.4454 3.3348
Fitted parameters
-2.042 -1.725 -20.454 -18.539
88.533 84.587 63.782 61.853
99.915 82.538 166.213 150.167
Table 5. Summary Statistics and Fitted Parameters of Daily Returns of S&P 500 and CSI 300 Indices from 2005 to 2016. The number of samples, mean, variance, skewness, and excess kurtosis are denoted by , , , , and , respectively. The mean and variance are computed from percentage returns. We assume for convenience.
Risk Measures S&P 500 CSI 300
Normal Sample Normal Sample
VaR () -1.997 -1.825 -1.824 -1.832 -2.995 -3.032 -3.036 -3.007
VaR () -2.836 -3.405 -3.432 -3.615 -4.254 -5.234 -5.246 -5.732
ES () -2.511 -2.857 -2.872 -3.042 -3.767 -4.433 -4.440 -4.745
ES () -3.253 -4.781 -4.820 -5.309 -4.879 -6.849 -6.857 -7.298
Table 6. Value-at-risk (VaR) and Expected Shortfall (ES) from the Normal Distribution (Normal), the Two NSVh Distributions ( and ), and the True Values from the Dataset (Sample).

Finally, the goodness-of-fit of the () is illustrated by using the probability plot. In Figure 4, the theoretical -score— for the -th ordered sample—is shown on the -axis and two different sample’s -scores on the -axis as follows: (i) from the estimated normal distribution and (ii) from the distribution computed as

Therefore, is understood as the probability plot, while is the normal probability plot in the usual definition. Figure 4 shows that the points under the probability plot are close to the line, indicating that the stock returns closely follow the distribution ().

Figure 4. Probability Plots of Daily returns of the S&P 500 (left) and CSI 300 Indexes (right). The -scores under the normal distribution (black dot) and the distribution (red circle) in -axis are plotted against the theoretical -score in the -axis. The line (dashed blue) is for reference.

5. 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 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 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 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.


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).


  • Alili and Gruet [1997] Larbi Alili and J. C. Gruet. An explanation of a generalized Bougerol’s identity in terms of hyperbolic Brownian motion. In Yor [1997], pages 15–33.
  • Alili et al. [1997] Larbi Alili, Daniel Dufresne, and Marc Yor. Sur l’identité de Bougerol pour les fonctionnelles exponentielles du mouvement brownien avec drift. In Yor [1997], pages 3–14.
  • Antonov et al. [2013] Alexander Antonov, Michael Konikov, and Michael Spector. SABR spreads its wings. Risk, 2013(8):58–63, 2013.
  • Antonov et al. [2015] Alexandre Antonov, Michael Konikov, and Michael Spector. Mixing SABR models for negative rates. Available at SSRN, 2015. URL https://ssrn.com/abstract=2653682.
  • Bachelier [1900] L. Bachelier. Théorie de la spéculation. Annales Scientifiques de L’Ecole Normale Supérieure, 17:21–88, 1900.
  • Badrinath and Chatterjee [1988] Swaminathan G Badrinath and Sangit Chatterjee. On measuring skewness and elongation in common stock return distributions: The case of the market index. Journal of Business, pages 451–472, 1988.
  • Balland and Tran [2013] Philippe Balland and Quan Tran. SABR goes normal. Risk, 2013(6):76–81, 2013.
  • Behr and Pötter [2009] Andreas Behr and Ulrich Pötter. Alternatives to the normal model of stock returns: Gaussian mixture, generalised logF and generalised hyperbolic models. Annals of Finance, 5(1):49–68, 2009.
  • Biller and Ghosh [2006] Bahar Biller and Soumyadip Ghosh. Multivariate input processes. In Shane G Henderson and Barry L Nelson, editors, Handbooks in operations research and management science: Simulation, volume 13, pages 123–153. Elsevier, 2006.
  • Bougerol [1983] Philippe Bougerol. Exemples de théorèmes locaux sur les groupes résolubles. Annales de l’Institut Henri Poincaré, B: Probability and Statistics, 19(4):369–391, 1983.
  • Box and Muller [1958] George EP Box and Mervin E Muller. A note on the generation of random normal deviates. The Annals of Mathematical Statistics, 29(2):610–611, 1