1. Introduction
Stochastic volatility (SV) models have been proposed to overcome the failure of the BlackScholesMerton (BSM) model in explaining nonconstant 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 heavytailed distribution generated by a continuoustime 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 fixedincome market traders, the normal model has long been a popular alternative to the BSM model for quoting and riskmanaging 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 alphabetarho (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.1.2. Skewed and heavytailed distribution
The second motivation of our study is that the normal SV models can serve as a means to generate distributions with skewness and heavytail, generalizing the normal distribution. Heavytailed 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 wellknown 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 continuoustime processes whose transition probabilities correspond to the following wellknown probability distributions: the arithmetic BM to a normal distribution (by definition), geometric BM to a lognormal distribution, and CEV and CIR processes to noncentral
distributions.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 ^{1}^{1}1The 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 closedform MonteCarlo (MC) simulation scheme that no longer needs a timediscretized 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 heavytailed 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 recognizedas an options pricing model—various quantities of interest, such as the vanilla option price, density functions, skewness, exkurtosis, valueatrisk, and expected shortfall, have closedform 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.
2. Models and Preliminaries
2.1. NSVh Model
The NSVh model is introduced as
(1) 
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 RadonNikodym 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 threedimensional hyperbolic geometry ().
To provide a background for the main result in section 3, we simplify the SDEs into the canonical forms,
(2) 
where the following changes of variables are used:
Here, the new variable is the integrated variance of the log volatility, and are the nondimensionalized 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 timetoexpiry 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
(3) 
It must be noted that the integral of , with respect to , is further simplified to the BM timechanged with an exponential functional of the BM defined by
(4) 
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 timeintegral 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 (
3).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:
(5) 
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 onefitsall 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 smalltime 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:
(6) 
where is the strike price, and denotes the timetoexpiry. 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 predetermined
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) pathdependent 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 outofthemoney 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 zerocorrelation case, the price process can be transformed to the CEV process timechanged with in a manner similar to that of Equation (3). Thus, the option price is expressed by a multidimensional 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 twodimensional 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é halfplane. Since the NSVh model also benefits from the same tools, we briefly introduce these tools. The dimensional Poincaré halfplane 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 LaplaceBeltrami 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.
Dimension  

Metric  
Volume element  
Geodesic distance  
to  
LaplaceBeltrami  
operator  
Standard BM  
Heat kernel  
for or  
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. HenryLabordère (2005, 2008) express the vanilla options price under the normal SABR model with a twodimensional integral, although it is later corrected by Korn and Tang (2013). Antonov et al. (2015) further simplify the price to a onedimensional 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 timeintegrated variance , conditional on the terminal volatility
. The cumulative distribution function (CDF) of the quantity is obtained from the Laplace transform of
, which has a closedform 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 rootsolving 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 hyperbolicgeometry extension of the BoxMuller 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 :
(7) 
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 , semibounded 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 heavytailed data set, and has been adopted in various fields; refer to Jones (2014) and the references therein. Examples in finance includes heavytailed innovation in the GARCH model (Choi and Nam, 2008), prediction of valueatrisk (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 heavytailed 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 GramCharlier or CornishFisher 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 quantilebased 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 firstclasscitizen status to the distribution among other heavytailed distributions by showing that it is a solution of a continuoustime 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:
(8) 
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 nonzero 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. MonteCarlo 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
(9) 
then the following is equal in distribution, conditional on :
(10) 
where is a twodimensional Bessel process, that is, the radius of a BM in twodimensional 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
(11)  
Furthermore, the three independent random variables can be simulated as
(12) 
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 pathdependent 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:
(13) 
Proof.
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 leftmost expression is shown by the timereversal . 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 reparametrized distribution:
(14) 
where the original parameters are mapped by
It also admits a simpler form:
(15) 
Proof.
The results are easily proved from the following hyperbolic function identities,
∎
Although Proposition 1 is a wellknown 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 pathdependent claims.
The key consequence of Corollary 3 is that the NSVh model bridges the SDEbased SV model and heavytailed 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 wellknown PDF and CDF of the distribution are respectively expressed by
(16) 
Conversely, from the distribution, the NSVh model obtains analytic tractability for option price and risk measures. Below are the closedform expressions for the quantities of interest:
Corollary 4.
For an asset price following the NSVh process with , option price, valueatrisk, and expected shortfall have the following closedform solutions.

The undiscounted price of a vanilla option with strike price :
(17) where is defined in Equation (16) and indicates call/put options, respectively.

Valueatrisk for quantile :
(18) 
Expected shortfall for quantile :
(19)
Proof.
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 closedform 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
(20)  
The central moments of the original form can be scaled as , and the skewness and exkurtosis are given as and , respectively. For the normal SABR (), further simplified expressions are obtained:
(21) 
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 exkurtosis. The result for large , at least for and 4 implies that the leading order of skewness and exkurtosis are independent of as and , as indicated in Equation (21). To illustrate that, in Figure 1, the contours of skewness and exkurtosis for and 1, as functions of (exkurtosis) and (skewness) for , are plotted. Although the parameters for is slightly higher than those for to obtain the same skewness and exkurtosis levels, the contours for and are very similar implying the similarity between the normal SABR and the distributions.
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 :
(22) 
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 rootfinding. 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 rootfinding, 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.
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) 
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 heavytailed 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 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 ()  

Parameters  
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 
(bps)  

200  2.275E2  4.1E5 1.8E5  2.274E2  6.4E7 1.8E5 
100  1.506E2  2.1E5 1.6E5  1.506E2  6.0E7 1.6E5 
0  9.083E3  1.2E5 1.3E5  9.083E3  5.3E7 1.3E5 
100  5.108E3  2.6E5 1.1E5  5.108E3  3.3E7 1.1E5 
200  2.807E3  4.8E5 8.8E6  2.804E3  3.9E7 9.1E6 
300  1.567E3  6.0E5 7.0E6  1.559E3  3.7E7 7.3E6 
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 12year 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 momentmatching methods of Tuenter (2001) and section 3.3 are used. Based on these values, the valueatrisk 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 valueatrisk 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 
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 
Finally, the goodnessoffit 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 ().
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 wellknown option pricing model in financial engineering, and the Johnson’s distribution is a popular skewed and heavytailed distribution defined by a transformation from the normal random variable. From the generalizations of Bougerol’s identity, the NSVh model is equipped with closedform MC simulation for general cases and closedform 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.
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 AsiaPacific 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. 2017001779, A machine learning and statistical inference framework for explainable artificial intelligence).
References
 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
Comments
There are no comments yet.