Hyperbolic normal stochastic volatility model
HHYPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL
JAEHYUK CHOI
Peking University HSBC Business School
CHENRU LIU
Department of Management Science and Engineering, Stanford University
BYOUNG KI SEO
Ulsan National Institute of Science and Technology
Abstract.
For option pricing models and heavy-tailed distributions, this study proposes acontinuous-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 gener-alized Bougerol’s identities in the literature, the study shows that our model has a closed-formMonte-Carlo simulation scheme and that the transition probability for one special case followsJohnson’s S U distribution—a popular heavy-tailed distribution originally proposed withoutstochastic process. It is argued that the S U distribution serves as an analytically superioralternative to the normal SABR model because the two distributions are empirically similar. 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 priceson 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 onthe geometric Brownian motion (BM) (hereafter, lognormal SV models). On the other hand, thestudies 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 SVmodels. 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-taileddistribution generated by a continuous-time stochastic process.
E-mail addresses : [email protected], [email protected], [email protected] . Date : September 7, 2018.
Key words and phrases. stochastic volatility, SABR model, Bougerol’s identity, Johnson’s S U distribution. a r X i v : . [ q -f i n . M F ] S e p CHOI, LIU, AND SEO
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 firsttime as an options pricing model by Bachelier [1900] (hereafter, normal model) and still providesmore relevant dynamics than the geometric BM for some financial asset classes. Refer to Brooksand Brooks [2017] and Schachermayer and Teichmann [2008] for recent surveys on the normalmodel. 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 theasset price as opposed to relative change. One example of the applications of the normal modelis its use for modeling the interest rate. The proportionality between the daily changes andlevel 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 popularalternative to the BSM model for quoting and risk-managing the options on interest rate swapand 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 theweighted average of the implied normal volatilities of the US Treasury bond options. It is alsoworth noting that the hedging ratio, delta, from the normal and BSM models can often besignificantly different, even after the volatilities of the corresponding models are calibrated tothe same option price observed on the market. Therefore, the normal model’s delta providesa more efficient hedge when the fluctuation of the underlying asset price is more consistent inabsolute term than in percentage term. The use of the normal model for the interest market isfurther justified by the negative policy rates observed in several developed economies after theglobal financial crisis of 2008. Other than the interest rate, the normal model is often used formodeling 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 normaldiffusion framework with the log price transformation; it implies that any existing results onthe 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 contextof the stochastic alpha-beta-rho (SABR) model [Hagan et al., 2002]—an SV model popularamong 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 modelprovides a range of backbone choices, including the normal and lognormal backbones. TheSABR model with normal backbone (hereafter, normal SABR) is an important motivation forthis study. A detailed review on the SABR model is provided in section 2.2.1.2.
Skewed and heavy-tailed distribution.
The second motivation of our study is thatthe 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 theirimportance cannot be emphasized enough. In this regard, the study of normal SV modelshas a much broader significance than that of the lognormal SV models. This is because the
YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 3 latter generalizes the lognormal distribution whose application is limited when compared to thenormal distribution.Several distribution families have been proposed in statistics to incorporate skewness andheavy tails into a normal distribution. Even if the focus is narrowed to the applications tofinance, it can be found that numerous distributions have been adopted to describe the statis-tics of asset return: generalized lambda [Corlu and Corlu, 2015], stable [Fama, 1965], skewed t [Theodossiou, 1998], Gaussian mixture [Kon, 1984, Behr and P¨otter, 2009], generalized hy-perbolic [Eberlein and Keller, 1995, Behr and P¨otter, 2009], Turkey’s g - and h - [Badrinathand Chatterjee, 1988, Mills, 1995], and Johnson’s S U (hereafter S U ) distribution [Shang andTadikamalla, 2004, Gurrola, 2007, Choi and Nam, 2008].However, the above distributions are neither defined from or associated with stochastic dif-ferential equations (SDEs), not to mention the SV models in particular. Those distributionsare 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 analyticallytractable solution. There are only a few examples of continuous-time processes whose transitionprobabilities correspond to the following well-known probability distributions: the arithmeticBM to a normal distribution (by definition), geometric BM to a lognormal distribution, andCEV and CIR processes to non-central χ distributions.1.3. Contribution of this study.
The study proposes and analyzes a class of normal SVmodels, which includes the normal SABR model as a special case. Since the mathematicsbehind our model involves the BMs in hyperbolic geometry and the results are expressed byhyperbolic functions, the class is named hyperbolic normal SV or NSVh model. The importantmathematical 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 thatno longer needs a time-discretized Euler scheme. The MC scheme requires merely one anda 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 theprevious exact MC scheme of Cai et al. [2017]. Additionally, the original proof of the firstgeneralization [Alili and Gruet, 1997] is simplified in this study. The second generalizationshows that a special case of the NSVh model— different from the normal SABR model—givesrise to the S U distribution [Johnson, 1949], one of the popular heavy-tailed distributions. Thisallows the study to add to the literature one rare example of analytically tractable SDEs. Thenormal 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 thepopular 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 shownthat the two distributions are very close to each other when parameters are estimated from the The class is named as an abbreviation in a manner similar to the way hyperbolic sine becomes sinh
CHOI, LIU, AND SEO same data set, and can thus be used interchangeably. Among the benefits of the interchangeableusage 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-formexpressions that are not available in other SV models. To facilitate the interchangeability, aquick method of moments to convert the equivalent parameter sets between the two distributionsis proposed.This remainder of this paper is organized as follows. Section 2 defines the NSVh model andreviews the SABR model and S U distribution. Section 3 describes the main results. Section 4presents the numerical results with empirical data. Finally, Section 5 concludes the paper.2. Models and Preliminaries
NSVh Model.
The NSVh model is introduced as dF t = σ t (cid:16) ρ dZ [ λα/ t + ρ ∗ dX t (cid:17) and dσ t σ t = α dZ [ λα/ t , (1)where F t and σ t are the processes for the price and volatility, respectively, α is the volatilityof the volatility parameter, ρ denotes the instantaneous correlation between F t and σ t , and ρ ∗ = (cid:112) − ρ . The BMs Z t and X t are independent, and Z [ µ ] t = Z t + µ t 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 theequity market—is achieved by a negative ρ , although it is in the context of the normal volatilityin the NSVh model. The parameter α accounts for the heavy tail, that is, excess kurtosis orvolatility smile. It can be easily seen that the process converges to an arithmetic BM in thelimit α → ρ . Therefore, α affects both the skewness and heavy tail at the sametime.The parameter λ is present in the drift of Z t for both F t and σ t . With regards to the volatilityprocess, λ controls the power of σ t that becomes a martingale as a geometric BM: d ( σ t ) − λ ( σ t ) − λ = (1 − λ ) α dZ t if λ (cid:54) = 1 ,d (log σ t ) = α dZ t if λ = 1 . For example, λ = 0 yields the volatility σ t , following a driftless geometric BM, as in the SABRmodel, and λ = − σ t , following a driftless geometric BM as in the SV modelof Hull and White [1987]. With regards to the price process, however, the drift prevents F t frombeing a martingale except for λ = 0 or, less importantly, ρ = 0, although the expectation iseasily computed as ¯ F t = F + ( σ ρ/α ) (cid:0) e λα t/ − (cid:1) . Therefore, the resulting process may notbe desirable as a price process. The NSVh model for λ (cid:54) = 0 is understood as a probabilitydistribution perturbed from the λ = 0 case, by applying the Radon-Nikodym derivative withrespect to Z t . Essentially, the introduction of λ does not significantly diversify the shape of the YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 5 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 subjectsseparately studied: the normal SABR model ( λ = 0), Johnson’s S U distribution ( λ = 1), andBM on three-dimensional hyperbolic geometry ( λ = − d ˜ F s = ˜ σ s ( ρ dZ [ λ/ s + ρ ∗ dX s ) and d ˜ σ s ˜ σ s = dZ [ λ/ s (˜ σ = 1) , (2)where the following changes of variables are used: s = α t, ˜ σ s = σ t σ , and ˜ F s = ασ (cid:0) F t − ¯ F T (cid:1) . Here, the new variable s is the integrated variance of the log volatility, ˜ F s and ˜ σ s are the non-dimensionalized price and volatility processes, respectively, under s . The scaling of ˜ F s and ˜ σ s naturally follows from the time change of the BMs with the new variable s . The time T is anyfixed time of interest such as the time-to-expiry of the vanilla option. The price F t is shifted by¯ F T first to ensure that E ( ˜ F S ) = 0 at the corresponding time S = α T . Therefore, the canonicalNSVh distribution is effectively parametrized by the three parameters—( S, ρ, λ ). Additionally,the original distribution is recovered by F T = ( σ /α ) ˜ F α T + ¯ F T and σ T = σ ˜ σ α T . While theoriginal and canonical representations are explicitly distinguished, the variable S is often usedin the original form as well for the sake of concise notation.The stochastic integrals of the canonical forms up to s = S are, respectively, expressed as˜ F S = ρ (cid:0) e Z [( λ − / S − e λS (cid:1) + ρ ∗ X A [( λ − / S and ˜ σ S = exp (cid:16) Z [( λ − / S (cid:17) . (3)It must be noted that the integral of ˜ σ s , with respect to X s , is further simplified to the BMtime-changed with an exponential functional of the BM defined by A [ µ ] T = (cid:90) Tt =0 e Z [ µ ] t dt ( A T = A [0] T ) . (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 continuouslyaveraged price under the BSM model in Asian options, it is used for the time-integral of thevariance in the context of this study. Although A [ µ ] T can be defined with any standard BM, weimplicitly assume that A [ µ ] T is tied to a particular BM, Z t , throughout this study. Essentially, A [ µ ] T and Z T are closely intertwined, and the knowledge of their joint distribution of is the keyto solve Equation (3).2.2. SABR Model and Hyperbolic Geometry.
The SABR model is reviewed with a focuson the normal backbone along with the BM on hyperbolic geometry, which serves as a mathe-matical tool for the NSVh and SABR models. The SABR model [Hagan et al., 2002] is an SVmodel with the backbone of the CEV process: dF t F βt = σ t ( ρdZ t + ρ ∗ dX t ) and dσ t σ t = α dZ t , (5) CHOI, LIU, AND SEO where X t and Z t are independent BMs. As mentioned earlier, the normal SABR model with β = 0 is equivalent to the NSVh model with λ = 0.The SABR model has been widely used in the financial industry, for covering fixed incomein particular, due to several merits: (i) arbitrary backbone choice, including normal ( β = 0)and lognormal ( β = 1) ones, (ii) availability of an approximate but fast vanilla options pricingmethod [Hagan et al., 2002], and (iii) parsimonious and intuitive parameters. The commentson those merits are presented in order. Regarding the CEV backbone, the popularity of theSABR model provides another evidence that the lognormal backbone of the BSM model isnot 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 limitof β → + , 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 canbe quickly computed through the BSM formula. However, it is worth noting that the normalvolatility is first obtained from the small-time perturbation of the normal diffusion even for β (cid:54) = 0 and subsequently it is converted to the BSM volatility by another approximation [Haganand Woodward, 1999]. Therefore, the option price computed from the normal volatility and thenormal model formula has been considered more accurate because the second approximation canbe avoided. Since the normal volatility is more appropriate for this study, the normal volatilityapproximation for β = 0 is presented for reference and later use: σ N ( σ , α, ρ, K ) = σ (cid:18) ζχ (cid:19) (cid:18) − ρ α T (cid:19) where ζ = ασ ( F − K ) and χ = log (cid:32) (cid:112) − ρζ + ζ − ρ + ζ − ρ (cid:33) , (6)where K is the strike price, and T denotes the time-to-expiry. The volatility approximationis an asymptotic expansion that is valid when α T is small; therefore, the accuracy of theapproximation noticeably deteriorates with an increase in α T . Despite the shortcoming, theinaccurate approximation does not cause problems in pricing vanilla options because the modelparameters σ , α, ρ and the pre-determined β are to be calibrated to the option prices observedfrom the market. In this regard, the implied volatility formula rather serves as an interpolationmethod for the volatility smile. Inaccurate approximation starts causing issues only when theusage of the model goes beyond vanilla options pricing. Two such cases are as follows:(i) claimswith 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 impliedfrom 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 isnot consistent with that from the MC simulation with the same parameters. Therefore, theparameter calibration for MC scheme should be performed with extra care. Hence, the researchon the accurate option analytics and efficient MC simulation methods come after the SABRmodel establishes its popularity among practitioners. YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 7
The study reviews prior works on the SABR model. Concerning vanilla options pricing, therehave been various improvements to Hagan et al. [2002]’s result. A few examples of such studiesare Ob(cid:32)l´oj [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 specialcases: (i) zero correlation ( ρ = 0), (ii) lognormal SABR ( β = 1), and (iii) normal SABR( β = 0). For the rest of the parameter ranges, no analytic solution is reported. Hence, the finitedifference method [Park, 2014, Le Floc’h and Kennedy, 2017] is considered the most practicalapproach. Concerning the zero-correlation case, the price process can be transformed to theCEV process time-changed with A [ − / T in a manner similar to that of Equation (3). Thus, theoption price is expressed by a multi-dimensional integral representation of A [ − / T over the CEVoption prices [Schroder, 1989]. Refer to Antonov et al. [2013] for the most simplified expressionbased on the heat kernel on the two-dimensional hyperbolic geometry [McKean, 1970], whichwe introduce below. The solution for the lognormal SABR is expressed in terms of the Gaussianhypergeometric series [Lewis, 2000].The options pricing under the normal SABR model depends on the mathematical tools devel-oped for the BMs on hyperbolic geometry, represented as Poincar´e half-plane. Since the NSVhmodel also benefits from the same tools, we briefly introduce these tools. The n -dimensionalPoincar´e half-plane is denoted by H n . Table 1 provides a quick reference for the properties of H and H . The standard BM in a geometry is defined to be the stochastic process whose infini-tesimal generator is given by the Laplace-Beltrami operator ∆ of the geometry. The heat kernel p ( t, D ) is the fundamental solution of the diffusion equation ( ∂ t − ∆) p ( t, D ) = 0, and hencethe transition probability of the standard BM. Table 1 shows the heat kernels on H [McKean,1970], often referred to as the McKean kernel, and H [Debiard et al., 1976]. The analyticalexpressions for the heat kernels are also known for H n , in general; see Grigor’yan and Noguchi[1998] for derivation.The standard BM on H is equivalent to the normal SABR with ρ = 0 in canonical form,where the x -axis is for the price process and the z -axis for the volatility process. Naturally, the H heat kernel has been used for the analysis of the normal SABR. Henry-Labord`ere [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 simplifythe price to a one-dimensional integration with an approximation. However, in the absenceof efficient numerical schemes to evaluate those integral representations, the normal volatilityapproximation in Equation (6) remains as a practical approach to price vanilla options underthe 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 beenproposed, an exact simulation method Cai et al. [2017] is available for the following threespecial cases: (i) ρ = 0, (ii) β = 0, and (iii) β = 1. The key element in this method is tosimulate the time-integrated variance A [ − / S , conditional on the terminal volatility Z S . Thecumulative distribution function (CDF) of the quantity is obtained from the Laplace transform CHOI, LIU, AND SEO
Table 1.
Properties of the n -Dimensional Hyperbolic Geometry Representedby n -Dimensional Poincar´e Half-Plane H n for n = 2 and 3. Symbols ∂ x and ∂ x are the shortened notations for partial derivative operators ∂∂x and ∂ ∂x , respec-tively.Dimension H = { ( x, z ) : z > } H = { ( x, y, z ) : z > } Metric ( ds ) ( dx + dz ) /z ( dx + dy + dz ) /z Volume element dV dx dz/z dx dy dz/z Geodesic distance D acosh (cid:16) ( x (cid:48) − x ) + z + z (cid:48) zz (cid:48) (cid:17) acosh (cid:16) ( x (cid:48) − x ) +( y (cid:48) − y ) + z (cid:48) + z zz (cid:48) (cid:17) ( x, · , z ) to ( x (cid:48) , · , z (cid:48) )Laplace-Beltrami z ( ∂ x + ∂ z ) z ( ∂ x + ∂ y + ∂ z ) − z ∂ z operator ∆ H n Standard BM dx t = z t dX t , dx t = z t dX t , dy t = z t dY t ,dz t /z t = dZ t dz t /z t = dZ t − dt/ p n ( t, D ) √ e − t/ (2 πt ) / (cid:90) ∞ D ds se − s / t √ cosh s − cosh D πt ) / D sinh D e − ( t + D ) / t for n = 2 or 3( ∂ t − ∆ H n ) p n = 0of (1 /A [ − / S ) | Z S , which has a closed-form expression [Matsumoto and Yor, 2005a]. Given theexact random numbers of Z S and A [ − / S , X A [ − / S in the normal SABR model is easily simulatedas X (cid:113) A [ − / S for a standard normal variable X . Although a heavy Euler scheme is avoided,the method of Cai et al. [2017] still incurs a moderate computation cost due to the numericalinversion of the Laplace transform and root-solving for the CDF inversion.As opposed to the previous studies using H , this study’s main result in section 3.1, based onAlili and Gruet [1997], exploits H . The pairs ( x t , z t ) and ( y t , z t ) from the standard BM in H are the two NSVh processes with λ = − ρ = 0, and it should be noted that the introductionof λ to the NSVh model makes this connection possible. Despite a higher dimensionality, theheat kernel of H is given without an integral, and Alili and Gruet [1997] shows that thesquared radius x t + y t and z t can be exactly simulated. Subsequently, x t (or y t ) is extracted viacosine (or sine) projection with random angle. Therefore, the study’s simulation scheme can beconsidered as a hyperbolic-geometry extension of the Box-Muller algorithm [Box and Muller,1958] for generating normal random variable, in which the z -axis of H is additionally added.2.3. Johnson’s Distribution Family.
Johnson [1949] proposes a system of distribution fami-lies in which a random variable X is represented by the transformations from a standard normalvariable Z : X − γ X δ X = f (cid:18) Z − γ Z δ Z (cid:19) for f ( x ) = / (1 + e − x ) for S B (bounded) family e x for S L (lognormal) familysinh x for S U (unbounded) family , (7) YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 9 where γ X and γ Z are location parameters and δ X and δ Z are scaling parameters. Althoughnot explicitly included, normal distribution can be considered as a special intersection of thethree families in the limit of δ X and δ Z proportionally going to infinity. Therefore, it is oftenincluded 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 uniquefamily is chosen for any mathematically feasible pair of skewness and kurtosis. For a fixed valueof 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 dataset, 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 assetreturn 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 accom-modate arbitrary high values of kurtosis, which is not feasible in the classical approaches togeneralize a normal distribution, such as the Gram-Charlier or Cornish-Fisher expansions. Sec-ond, 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 themoment 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 simula-tions, particularly in a multivariate setting [Biller and Ghosh, 2006]. In general, random numbersampling 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 tothe S U distribution among other heavy-tailed distributions by showing that it is a solution of acontinuous-time SV process, the NSVh model with λ = 1. This partially explains why the S U 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. Sincethe original identity is generalized later, we state it as a Corollary and defer the proof toProposition 2.
Corollary 1 (Bougerol’s identity) . For a fixed time T , the following is equal in distribution: (cid:90) T e Z t dX t d = X A T d = sinh( W T ) , (8) 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 BMsis 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 interpretedwith caution; the equality holds as distribution ( d = ) at a fixed time t = T , not as a process for ≤ t ≤ T . Moreover, it does not directly help to solve Equation (3). The identity must begeneralized to non-zero drift, A [ µ ] , and provide the joint distribution with Z [ µ ] , which is foundin Alili and Gruet [1997] and Alili et al. [1997]. In the following subsections, we apply twogeneralizations to the NSVh model; one to the general λ , and the other to a special case λ = 1.3.1. Monte-Carlo Simulation Scheme.
The first generalization makes use of the BMs in H to obtain the distribution of X A [ µ ] T , conditional on Z [ µ ] T . We restate Proposition 3 in Alili andGruet [1997] in a modified and stronger form: Proposition 1 ( Bougerol’s identity in hyperbolic geometry ) . Let X t and Z t be twoindependent BMs and the function φ defined by φ ( Z, D ) = e Z/ √ D − Z for Z ≤ D, (9) then the following is equal in distribution, conditional on Z [ µ ] T : (cid:90) T e Z [ µ ] t dX t d = X A [ µ ] T d = cos θ φ (cid:18) Z [ µ ] T , (cid:113) R T + ( Z [ µ ] T ) (cid:19) , (10) where R t is a two-dimensional Bessel process, that is, the radius of a BM in two-dimensionalEuclidean geometry, and θ is a uniformly distributed random angle. The three random variables R T , Z T , 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 H is provided. The proposition in this study isstronger than the original because the distribution equality holds, conditional on Z [ µ ] , whichis implied from the proof without difficulty. The study further simplifies the original proof byusing the H heat kernel. Proposition 1 can be directly applied to Equation (3) to obtain thetransition equation and MC scheme: Corollary 2.
The joint distribution of the NSVh model at a fixed time S is given as ˜ σ S = exp (cid:16) Z [( λ − / S (cid:17) and ˜ F S d = ρ (cid:16) e Z [( λ − / S − e λS/ (cid:17) + ρ ∗ cos θ φ (cid:18) Z [( λ − / S , (cid:113) R S + ( Z [( λ − / S ) (cid:19) . (11) Furthermore, the three independent random variables can be simulated as (cid:0) Z S , R S , cos θ (cid:1) d = (cid:32) Z √ S, ( X + Y ) S, X (or Y ) (cid:112) X + Y (cid:33) , (12) where X and Y are independent standard normals. The simulation method using standard normal variables in Equation (12) is more efficientthan drawing R S and θ independently because the costly cos θ evaluation is avoided. The ideais similar to the Marsaglia polar method [Marsaglia and Bray, 1964] for drawing normal randomvariables. It must be noted that three random numbers— X , Y , and Z —generate two pairsof ˜ F S and ˜ σ S . 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 YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 11 is more efficient than the exact SABR simulation of Cai et al. [2017], although it is limited tothe normal case. This study’s method directly draws X A [( λ − / S and Z S , whereas Cai et al.[2017] first draws A [1 / S and Z S , and subsequently X A [ − / S . Although Equation (11) states thetransition from s = 0 to s = S , it can handle any time interval from s = S to s = S ( S < S ).Therefore, the scheme is ideal for pricing path-dependent claims.3.2. S U Distributions for λ = 1 . This subsection shows that the NSVh distribution for λ = 1is expressed by the S U distribution and is related to Bougerol’s identity generalized to anarbitrary 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 inProposition 1 of Alili et al. [1997] (or Proposition 2.1 of Vakeroudis [2012]), and we follow theproof therein. Proposition 2 ( Bougerol’s identity with an arbitrary starting point ) . For a fixed time T and independent BMs, X t , Z t , and W t , the following is equal in distribution: sinh( a ) e Z T + (cid:90) T e Z t dX t d = sinh( a ) e Z T + X A T d = sinh( W T + a ) . (13) Proof.
The two processes, P t = sinh( W t + a ) and Q t = e Z t (cid:18) sinh( a ) + (cid:90) t e − Z s dX s (cid:19) , are equivalent because they start from the same starting point P = Q = sinh( a ) and followthe SDE: dP t = 12 P t dt + (cid:113) P t dW t and dQ t = 12 Q t dt + dX t + Q t dZ t = 12 Q t dt + (cid:113) Q t dW t . Therefore, P t and Q t have the same distribution for any time t . The equality between Q T andthe left-most expression is shown by the time-reversal s → T − s . For a fixed time T , Z T − Z T − s for 0 ≤ s ≤ T is also a standard BM with the same ending point Z T , and therefore it may bereplaced with Z s . (cid:3) The original Bougerol’s identity in Corollary 1 is a special case, with a = 0. Now, Proposition2 can be applied to further simplify the NSVh distribution for λ = 1. Corollary 3.
The price of the NSVh model with λ = 1 at a fixed time S follows a re-parametrized S U distribution: ˜ F S d = ρ ∗ sinh ( W S + atanh ρ ) − ρ e S/ , (14) where the original parameters are mapped by δ Z = 1 √ S , γ Z δ Z = − atanh ρ, δ X = σ ρ ∗ α , and γ X = ¯ F T − σ ρα e S/ . It also admits a simpler form: ˜ F S d = sinh( W S ) + ρ (cid:16) cosh( W S ) − e S/ (cid:17) . (15) Proof.
The results are easily proved from the following hyperbolic function identities,asinh (cid:18) ρρ ∗ (cid:19) = atanh ρ = 12 log (cid:18) ρ − ρ (cid:19) . (cid:3) Although Proposition 2 is a well-known result, to the best of our knowledge, this is the firsttime that it is interpreted in the context of the S U distribution or SV model. Compared toCorollary 2, Corollary 3 is an even more efficient MC scheme requiring one normal randomnumber for one draw of price, although for a special case λ = 1. However, this is at the expenseof the terminal volatility ˜ σ S being lost. Unlike Corollary 2, Corollary 3 can only generate thefinal price ˜ F S , and therefore cannot be used for path-dependent claims.The key consequence of Corollary 3 is that the NSVh model bridges the SDE-based SVmodel and heavy-tailed distribution. The encounter between the two subjects is mutuallybeneficial. First, as a solution of the NSVh process, the S U distribution obtains a more intuitiveparametrization than the original in Equation (7). From Equation (15), it is clear that thesymmetric heavy tail comes from the sinh term, controlled by S , and the asymmetric skewnessfrom the cosh term, controlled by ρ . The new parametrization also helps to understand therelationship between Johnson family members. The lognormal family S L is recognized as aspecial case with ρ = ± ρ ∗ = 0) as ˜ F S d = ± ( e Z S − e S/ ). The normal family S N is obtainedas ˜ F S / √ S in the limit of S →
0. Under the NSVh parameters, the well-known PDF and CDFof the S U distribution are respectively expressed by p λ =1 ( x ) = n ( d ) ρ ∗ σ √ T (cid:112) ξ and P λ =1 ( x ) = N ( − d )where d = 1 √ S (cid:18) asinh (cid:18) αρ ∗ σ ( ¯ F T − x ) − ρρ ∗ e S/ (cid:19) + atanh ρ (cid:19) , (16)Conversely, from the S U distribution, the NSVh model obtains analytic tractability for optionprice 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 λ = 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 : V ± ( x ) = σ α e S/ (cid:16) (1 + ρ ) N ( d + √ S ) − (1 − ρ ) N ( d − √ S ) − ρN ( d ) (cid:17) ± (cid:0) ¯ F T − K (cid:1) N ( ± d ) , (17) where d is defined in Equation (16) and ± indicates call/put options, respectively. • Value-at-risk for quantile p :VaR ( p ) = ¯ F T − σ α (cid:16) ρ ∗ sinh (cid:16) d √ S − atanh ρ (cid:17) + ρ e S/ (cid:17) for d = − N − ( p ) . (18) • Expected shortfall for quantile p :ES ( p ) = ¯ F T − σ e S/ αp (cid:16) (1 + ρ ) N ( d + √ S ) − (1 − ρ ) N ( d − √ S ) − ρ (1 − p ) (cid:17) . (19) YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 13
Proof.
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 ex-pected shortfall, ES( p ) = VaR( p ) − V − (VaR( p )) /p , is useful, where VaR( p ) vanishes in the finalexpression for ES( p ). (cid:3) Later, it is argued that the NSVh distributions with different values of λ are close to eachother, and thus the analytically tractable λ = 1 case can represent the rest including the normalSABR model ( λ = 0). The option price from the closed-form formula serves as a benchmarkagainst 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 ofthe NSVh distribution for general λ to be used for parameter estimation. The study alsoproposes a moment matching in the reduced form for λ = 0 to complement that for λ = 1 byTuenter [2001]. Corollary 5.
The central moments of the canonical NSVh distribution ˜ µ n = E ( ˜ F nS ) , for ≤ n ≤ , are given as ˜ µ = ρ w λ ( w −
1) + ρ ∗ w λ −
11 + λ for w = e S ≥ , ˜ µ = ρ w λ ( w − ( w + 2) + 3 ρρ ∗ w λ (cid:18) w λ −
13 + λ − w λ −
11 + λ (cid:19) , and ˜ µ = ρ w λ ( w − ( w + 2 w + 3 w −
3) + 6 ρ ρ ∗ w λ (cid:18) w w λ −
15 + λ − w λ −
13 + λ + w λ −
11 + λ (cid:19) + 32 ρ ∗ (cid:18) − w λ w λ −
15 + λ + ( w λ + 1) w λ −
13 + λ − w λ −
11 + λ (cid:19) . (20) The central moments of the original form can be scaled as µ n = E (( F T − ¯ F T ) n ) = ( σ /α ) n ˜ µ n ,and the skewness and ex-kurtosis are given as s = ˜ µ / ˜ µ / and κ = ˜ µ / ˜ µ − , respectively. Forthe normal SABR ( λ = 0 ), further simplified expressions are obtained: ˜ µ = w − , s = ρ ( w +2) √ w − , κ = ( w − (cid:18)(cid:18) ρ + 15 (cid:19) ( w + 3 w + 6 w + 5) + 1 (cid:19) . (21)Refer to Appendix B for detailed derivation. Corollary 5 generalizes the moments for S L ( ρ = ±
1) and S U ( λ = 1) distributions.The similarity of the NSVh distributions for different λ is inferred from skewness and ex-kurtosis. The result ˜ µ k = O ( w k ( λ + k − / ) for large w , at least for k = 2 , , and 4 impliesthat the leading order of skewness and ex-kurtosis are independent of λ as s = O ( w / ) and κ = O ( w ), as indicated in Equation (21). To illustrate that, in Figure 1, the contours ofskewness and ex-kurtosis for λ = 0 and 1, as functions of S (ex-kurtosis) and ρS (skewness) for ρ ≥
0, are plotted. Although the parameters for λ = 0 is slightly higher than those for λ = 1 toobtain the same skewness and ex-kurtosis levels, the contours for λ = 0 and 1 are very similarimplying the similarity between the normal SABR and the S U distributions.Parallel to Tuenter [2001]’s reduced moment matching method for the S U distribution, asimilar method is developed for the normal SABR model. Combined with Tuenter [2001], the Figure 1.
Contour Plot of Skewness (Red Dashed Line) and Excess Kurtosis(Blue Solid Line) for Varying S (= α T ) versus ρS . The upper left triangle( ρS, S ) is for λ = 1 ( S U ) and the lower right triangle ( S, ρS ) for λ = 0 (normalSABR). The values for skewness are 0, 1.5, 3, 4.5, 6, and 8, and the values forexcess kurtosis are 2, 7, 16, 40, 100, and 200 from the lower left to the upperright corner. . . . . . . . S ( l = S ( l = ) r S ( l = ) r S ( l = l = l = two methods can quickly find an equivalent parameter set of one distribution from the other.By joining the expressions for s and κ in Equation (21) through ρ , κ is expressed as a univariatefunction of w ≥ f ( w ) = 4 s ( w + 3 w + 6 w + 5)5( w + 2) + ( w − (cid:18) w + 3 w + 6 w + 5) (cid:19) , (22)for which the study numerically finds the root w ∗ of κ = f ( w ∗ ). It can be shown that f ( w ) ismonotonically increasing for w ≥
1. Therefore, the root w ∗ would be unique if it exists. Wecan further bound w ∗ by w m ≤ w ∗ ≤ w M to expedite the numerical root-finding. Lower bound w m is the unique cubic root of s = ( w − w + 2) (the ρ = ± w ≥ w m = 2 cosh (cid:18)
13 acosh (cid:18) s (cid:19)(cid:19) . Upper bound w M is obtained by plugging w = w m into Equation (22), except the ( w −
1) term: w M = 1 + κ − s ( w m + 3 w m + 6 w m + 5) / ( w m + 2) ( w m + 3 w m + 6 w m + 5) . YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 15
Figure 2.
The Overview of the NSVh Model in Relation to Other PreviouslyKnown Models and Distributions.
𝑁𝑆𝑉ℎ 𝑀𝑜𝑑𝑒𝑙
𝑆𝐴𝐵𝑅 𝑀𝑜𝑑𝑒𝑙 𝛽 = 0 𝜆 = 0𝑁𝑜𝑟𝑚𝑎𝑙 𝑆𝐴𝐵𝑅 𝛽 = 1
𝐿𝑜𝑔𝑛𝑜𝑟𝑚𝑎𝑙 𝑆𝐴𝐵𝑅𝐽𝑜ℎ𝑛𝑠𝑜𝑛 ′ 𝑠 𝐷𝑖𝑠𝑡𝑟𝑖𝑏𝑢𝑡𝑖𝑜𝑛 𝑋𝑋 − 𝛾 𝑋 𝛿 𝑋 = 𝑓 𝑍 − 𝛾 𝑍 𝛿 𝑍 𝑓 𝑥= sinh 𝑥𝜆 = 1𝑆 𝑈 𝑢𝑛𝑏𝑜𝑢𝑛𝑑𝑒𝑑 𝑓𝑎𝑚𝑖𝑙𝑦 Table 2.
Summary of the NSVh Model for the Three Key Special Cases: λ = − ,
0, and 1.Drift parameter λ = − λ = 0 λ = 1Equivalent Standard BM in H Normal SABR S U distributionmodel/distribution Standard BM in H Drifted BM for ˜ σ s Z [ − / s Z s Z [1 / s Terminal volatility ˜ σ S exp( Z [ − S ) exp( Z [ − / S ) exp( Z S )Integrated variance A [ − S A [ − / S A [0] S Exact MC simulation Corollary 2 Corollary 2 & 3Vanilla option price Equation (6) Equation (17)(approximation) (exact)Moments Corollary 5Moment matching Equation (22) Tuenter [2001]The existence of w ∗ is equivalent to f ( w m ) ≤ κ . If w ∗ exists and is found from the numericalroot-finding, then the parameters can be solved as S = log w ∗ , ρ = s ( w ∗ + 2) √ w ∗ − , and σ α = (cid:115) µ log w ∗ ( w ∗ − S .
Summary of results.
In Figure 2, the relationship of the NSVh model and other relatedmodels is shown. In Table 2, the results for the three important drift values, λ = −
1, 0, and 1,are summarized for comparison.4.
Parameter Estimation from Empirical Data
The NSVh distribution is calibrated to two empirical data sets— swaption volatility smileand daily stock index return. The purpose of this exercise is to demonstrate various numericalprocedures presented in this study rather than arguing that the NSVh model is superior to
Figure 3.
Swaption Volatility Smile in the Implied Normal Volatility (inannual basis points) Observed in the US Market on March 14, 2017: (a) 1 yearinto 1-year swap (1y1y) and (b) 10 years into 10-year swap (10y10y). The circlesrepresent the volatilities implied from observed market prices, among which theblack ones are ATM and ATM ±
1% points for calibration. The solid (blue) linerepresents the model-implied smile curve for λ = 1 ( S U ) and the dashed (red)one for λ = 0 (normal SABR). The insets are the implied BSM volatilities (inannual %). Strike (%) N o r m a l i m p li ed v o l ( bp s ) Strike (%) N o r m a l i m p li ed v o l ( bp s ) other SV models or heavy-tailed distributions in fitting these data. Additionally, the studyshows that the two NSVh models, that is, λ = 0 (normal SABR) and λ = 1 ( S U ), yield verysimilar 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 March14, 2017 from Reuters. The two heavily traded expiry–tenor pairs of the US swaption—1y1yand 10y10y—are chosen to illustrate different volatility smile shapes. To avoid the complicationof the annuity price of the underlying swap, the study computes the price in the unit of theannuity 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, ± . ± . ± . ± .
5% from the forward swap rates ¯ F T , the study uses only the three spreads—0 and ±
1% for calibration—to ensure that the calibrated parameter set—( σ , α, ρ )—reproduces optionprices at those three strike prices. For calibration, the study uses Equations (6) for λ = 0 and(17) for λ = 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 λ = 0 and 1 for reference.In Table 4, the option prices are compared from the MC simulation of Corollary 2 with theoption prices from the aforementioned analytic methods. In the experiment, the MC simulation YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 17
Table 3.
Parameters Calibrated to the US Swaption Volatility Smile on March14, 2017. Refer to Figure 3 for the calibration points. Equations (6) and (17)are used for the NSVh models with λ = 0 and λ = 1, respectively.Calibrated 1y1y ( T = 1) 10y10y ( T = 10)Parameters λ = 0 λ = 1 λ = 0 λ = 1 ρ (%) 33.503 32.244 1.697 1.580 α (%) 61.962 62.181 22.372 22.196 σ (%) 0.533 0.477 0.691 0.609¯ F T (%) 2.0221 3.0673 Table 4.
Vanilla Options Pricing Tested Against the Parameters of Table 3 for10y10y Swaption. Analytic prices ( P ana ) are computed from Equations (6) and(17) and MC price ( P mc ) is shown relative to P ana with the standard deviation.The MC simulation with 10 paths is repeated 100 times. Prices are in the unitof the annuity of the underlying swap. K − ¯ F T λ = 0 λ = 1(bps) P ana P mc − P ana P ana P mc − P ana -200 2.275E-2 -4.1E-5 ± ± ± ± ± ± ± ± ± ± ± ± samples is repeated 100 times. For λ = 1, both Equation (17) and the MC method areexact, and thus the option prices from the two methods have very little difference due to theMC noise. For λ = 0, however, Equation (6) is an approximation. The analytic prices showclear deviation from the accurate MC prices. Overall, this exercise reconfirms the possibility ofthe S U distribution being used as a better alternative to the normal SABR model.4.2. Daily Return of Stock Index.
We fit the NSVh distribution to the daily returns of twostock 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 outrightindex values, rather than as holding period return.The statistics of the daily returns are summarized in Table 5. S&P 500 shows heavier tailsbut less skewness than CSI 300. The table also shows the fitted parameters of the NSVhdistributions for λ = 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 Table 5.
Summary Statistics and Fitted Parameters of Daily Returns of S&P500 and CSI 300 Indices from 2005 to 2016. The number of samples, mean,variance, skewness, and excess kurtosis are denoted by n , ¯ F T , µ , s , and κ ,respectively. The mean and variance are computed from percentage returns. Weassume T = 1 for convenience.Summary statistics S&P 500 CSI 300 n F T µ s -0.0933 -0.5075 κ λ = 0 λ = 1 λ = 0 λ = 1 ρ (%) -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 6.
Value-at-risk (VaR) and Expected Shortfall (ES) from the NormalDistribution (Normal), the Two NSVh Distributions ( λ = 0 and λ = 1), and theTrue Values from the Dataset (Sample).Risk Measures S&P 500 CSI 300Normal λ = 0 λ = 1 Sample Normal λ = 0 λ = 1 SampleVaR ( p = 5%) -1.997 -1.825 -1.824 -1.832 -2.995 -3.032 -3.036 -3.007VaR ( p = 1%) -2.836 -3.405 -3.432 -3.615 -4.254 -5.234 -5.246 -5.732ES ( p = 5%) -2.511 -2.857 -2.872 -3.042 -3.767 -4.433 -4.440 -4.745ES ( p = 1%) -3.253 -4.781 -4.820 -5.309 -4.879 -6.849 -6.857 -7.298are computed and compared to those from the normal distribution assumption and the historicaldata. While Equations (18) and (19) are used for λ = 1, MC simulation is used to compute therisk measures for λ = 0. The results are shown in Table 6. As the return distribution has heavytail, the value-at-risk and expected shortfall from the NSVh distributions are much closer tothose from the historical data than those from the normal distribution. The risk measures from λ = 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 ( λ = 1) is illustrated by using the probability plot.In Figure 4, the theoretical Z -score— Z ( j )0 = N − (( j − / /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) fromthe estimated normal distribution Z ( j )1 = ( X j − ¯ F T ) / √ µ and (ii) from the S U distributioncomputed as Z ( j )2 = N − ( P λ =1 ( X j )) = 1 √ S (cid:18) asinh (cid:18) αρ ∗ σ ( X j − ¯ F T ) + ρρ ∗ e S/ (cid:19) − atanh ρ (cid:19) . YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 19
Figure 4.
Probability Plots of Daily returns of the S&P 500 (left) and CSI 300Indexes (right). The Z -scores under the normal distribution (black dot) and the S U distribution (red circle) in y -axis are plotted against the theoretical Z -scorein the x -axis. The y = x line (dashed blue) is for reference. −2 0 2 − − Theoretical Quantile S a m p l e Q uan t il e −2 0 2 − − Theoretical Quantile S a m p l e Q uan t il e Therefore, ( Z , Z ) is understood as the S U probability plot , while ( Z , Z ) is the normal proba-bility plot in the usual definition. Figure 4 shows that the points under the S U probability plotare close to the y = x line, indicating that the stock returns closely follow the S U distribution( λ = 1). 5. Conclusion
The study generalizes the arithmetic Brownian motion with stochastic volatility. The NSVhprocess 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 pric-ing model in financial engineering, and the Johnson’s S U distribution is a popular skewed andheavy-tailed distribution defined by a transformation from the normal random variable. Fromthe generalizations of Bougerol’s identity, the NSVh model is equipped with closed-form MCsimulation for general cases and closed-form option formula for the S U case. This study demon-strates the usage of the model with two empirical datasets—the US swaption and daily returnsdistribution of the S&P 500 and CSI 300 indexes. Acknowledgements
The authors are grateful to Larbi Alili for sharing manuscripts, Alili and Gruet [1997] andAlili et al. [1997], Robert Webb (editor), and Minsuk Kwak (discussant at the 2018 Asia-PacificAssociation of Derivatives conference in Busan, Korea). Jaehyuk Choi would like to expressgratitude to his previous employer, Goldman Sachs, for laying the foundation for the motivationof this study. Byoung Ki Seo was supported by the Institute for Information & communicationsTechnology 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 explainableartificial intelligence).
References
Larbi Alili and J. C. Gruet. An explanation of a generalized Bougerol’s identity in terms ofhyperbolic Brownian motion. In Yor [1997], pages 15–33.Larbi Alili, Daniel Dufresne, and Marc Yor. Sur l’identit´e de Bougerol pour les fonctionnellesexponentielles du mouvement brownien avec drift. In Yor [1997], pages 3–14.Alexander Antonov, Michael Konikov, and Michael Spector. SABR spreads its wings.
Risk ,2013(8):58–63, 2013.Alexandre Antonov, Michael Konikov, and Michael Spector. Mixing SABR models for negativerates.
Available at SSRN , 2015. URL https://ssrn.com/abstract=2653682 .L. Bachelier. Th´eorie de la sp´eculation.
Annales Scientifiques de L’Ecole Normale Sup´erieure ,17:21–88, 1900.Swaminathan G Badrinath and Sangit Chatterjee. On measuring skewness and elongation incommon stock return distributions: The case of the market index.
Journal of Business , pages451–472, 1988.Philippe Balland and Quan Tran. SABR goes normal.
Risk , 2013(6):76–81, 2013.Andreas Behr and Ulrich P¨otter. Alternatives to the normal model of stock returns: Gaussianmixture, generalised logF and generalised hyperbolic models.
Annals of Finance , 5(1):49–68,2009.Bahar Biller and Soumyadip Ghosh. Multivariate input processes. In Shane G Hendersonand Barry L Nelson, editors,
Handbooks in operations research and management science:Simulation , volume 13, pages 123–153. Elsevier, 2006.Philippe Bougerol. Exemples de th´eor`emes locaux sur les groupes r´esolubles.
Annales del’Institut Henri Poincar´e, B: Probability and Statistics , 19(4):369–391, 1983.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, 1958.Robert Brooks and Joshua A Brooks. An option valuation framework based on arithmeticBrownian motion: Justification and implementation issues.
Journal of Financial Research ,40(3):401–427, 2017. doi:10.1111/jfir.12129.Ning Cai, Yingda Song, and Nan Chen. Exact simulation of the SABR model.
OperationsResearch , 65(4):931–951, 2017.Bin Chen, Cornelis W Oosterlee, and Hans Van Der Weide. A low-bias simulation schemefor the SABR stochastic volatility model.
International Journal of Theoretical and AppliedFinance , 15(2):1250016, 2012. doi:10.1142/S0219024912500161.Pilsun Choi and Kiseok Nam. Asymmetric and leptokurtic distribution for heteroscedastic assetreturns: the S U -normal distribution. Journal of Empirical finance , 15(1):41–63, 2008.Canan G Corlu and Alper Corlu. Modelling exchange rate returns: which flexible distributionto use?
Quantitative Finance , 15(11):1851–1864, 2015.
YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 21
Am´ed´ee Debiard, Bernard Gaveau, and Edmond Mazet. Th´eoremes de comparaison eng´eom´etrie riemannienne.
Publications of the Research Institute for Mathematical Sciences ,12(2):391–425, 1976.Ernst Eberlein and Ulrich Keller. Hyperbolic distributions in finance.
Bernoulli , 1(3):281–299,1995.Eugene F Fama. The behavior of stock-market prices.
The Journal of Business , 38(1):34–105,1965.Alexander Grigor’yan and Masakazu Noguchi. The heat kernel on hyperbolic space.
Bulletinof the London Mathematical Society , 30(6):643–650, 1998.Pedro Gurrola. Capturing fat-tail risk in exchange rate returns using S U curves: A comparisonwith the normal mixture and skewed student distributions. The Journal of Risk , 10(2):73,2007.Patrick S Hagan and Diana E Woodward. Equivalent Black volatilities.
Applied MathematicalFinance , 6(3):147–157, 1999.Patrick S Hagan, Deep Kumar, Andrew S Lesniewski, and Diana E Woodward. Managing smilerisk.
Wilmott Magazine , 2002(9):84–108, 2002.Pierre Henry-Labord`ere. A general asymptotic implied volatility for stochastic volatility models.
Available at SSRN , 2005. URL https://ssrn.com/abstract=698601 .Pierre Henry-Labord`ere.
Analysis, geometry, and modeling in finance: Advanced methods inoption pricing . CRC Press, 2008.Steven L Heston. A closed-form solution for options with stochastic volatility with appli-cations to bond and currency options.
Review of Financial Studies , 6(2):327–343, 1993.doi:10.1093/rfs/6.2.327.John Hull and Alan White. The pricing of options on assets with stochastic volatilities.
TheFournal of Finance , 42(2):281–300, 1987.Norman L Johnson. Systems of frequency curves generated by methods of translation.
Biometrika , 36(1/2):149–176, 1949.David L Jones. Johnson curve toolbox for Matlab: Analysis of non-normal data using theJohnson family of distributions., 2014.Richard Jordan and Charles Tier. Asymptotic approximations to deterministic and stochasticvolatility models.
SIAM Journal on Financial Mathematics , 2(1):935–964, 2011.Joanne E Kennedy, Subhankar Mitra, and Duy Pham. On the approximation of the SABRmodel: A probabilistic approach.
Applied Mathematical Finance , 19(6):553–586, 2012.doi:10.1080/1350486X.2011.646523.Chris Kenyon. Inflation is normal.
Risk , 2008(9):54–60, 2008.Stanley J Kon. Models of stock returns: A comparison.
The Journal of Finance , 39(1):147–165,1984.Ralf Korn and Songyin Tang. Exact analytical solution for the normal SABR model.
WilmottMagazine , 2013(7):64–69, 2013. doi:10.1002/wilm.10235.
Fabien Le Floc’h and Gary Kennedy. Finite difference techniques for arbitrage free SABR.
Journal of Computational Finance , 20(3):51–79, 2017. doi:10.21314/JCF.2016.320.´Alvaro Leitao, Lech A Grzelak, and Cornelis W Oosterlee. On an efficient multiple time stepMonte Carlo simulation of the SABR model.
Quantitative Finance , 17(10):1549–1565, 2017a.doi:10.1080/14697688.2017.1301676.´Alvaro Leitao, Lech A Grzelak, and Cornelis W Oosterlee. On a one time-step Monte Carlo sim-ulation approach of the SABR model: Application to European options.
Applied Mathematicsand Computation , 293:461–479, 2017b. doi:10.1016/j.amc.2016.08.030.Alexander Levin. Interest rate model selection.
The Journal of Portfolio Management , 30(2):74–86, 2004.Alan L Lewis.
Option Valuation Under Stochastic Volatility . Finance Press, Newport Beach,CA, 2000.Matthew Lorig, Stefano Pagliarani, and Andrea Pascucci. Explicit implied volatilities for multi-factor local-stochastic volatility models.
Mathematical Finance , 2015. doi:10.1111/mafi.12105.George Marsaglia and Thomas A Bray. A convenient method for generating normal variables.
SIAM Review , 6(3):260–264, 1964.Hiroyuki Matsumoto and Marc Yor. Exponential functionals of Brownian motion, I: Probabilitylaws at fixed time.
Probability Surveys , 2:312–347, 2005a.Hiroyuki Matsumoto and Marc Yor. Exponential functionals of Brownian motion, II: Somerelated diffusion processes.
Probability Surveys , 2:348–384, 2005b.Henry P McKean. An upper bound to the spectrum of ∆ on a manifold of negative curvature.
Journal of Differential Geometry , 4(3):359–366, 1970.Terence C Mills. Modelling skewness and kurtosis in the London stock exchange FT-SE indexreturn distributions.
The Statistician , pages 323–332, 1995.Jan Ob(cid:32)l´oj. Fine-tune your smile: Correction to Hagan et al.
Available at arXiv , 2007. URL https://arxiv.org/pdf/0708.0998 .Hyukjae Park. SABR symmetry.
Risk , 2014(1):106–111, 2014.Geoffrey Poitras. Spread options, exchange options, and arithmetic Brownian motion.
Journalof Futures Markets , 18(5):487–517, 1998.Walter Schachermayer and Josef Teichmann. How close are the option pricing formulas ofBachelier and Black-Merton-Scholes?
Mathematical Finance , 18(1):155–170, 2008.Mark Schroder. Computing the constant elasticity of variance option pricing formula.
Journalof Finance , 44(1):211–219, 1989. doi:10.1111/j.1540-6261.1989.tb02414.x.Jen S Shang and Pandu R Tadikamalla. Modeling financial series distributions: A versatile datafitting approach.
International Journal of Theoretical and Applied Finance , 7(03):231–251,2004.Jean-Guy Simonato. The performance of Johnson distributions for computing value at risk andexpected shortfall.
The Journal of Derivatives , 19(1):7–24, 2011.Elias M Stein and Jeremy C Stein. Stock price distributions with stochastic volatility: ananalytic approach.
The Review of Financial Studies , 4(4):727–752, 1991.
YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 23
Panayiotis Theodossiou. Financial data and the skewed generalized t distribution. ManagementScience , 44(12-part-1):1650–1661, 1998.Hans JH Tuenter. An algorithm to determine the parameters of S U -curves in the Johnsonsystem of probabillity distributions by moment matching. Journal of Statistical Computationand Simulation , 70(4):325–347, 2001.Stavros Vakeroudis. Bougerol’s identity in law and extensions.
Probability Surveys , 9:411–437,2012.Sree Vinutha Venkataraman and SVD Nageswara Rao. Estimation of dynamic VaR using JSUand PIV distributions.
Risk Management , 18(2-3):111–134, 2016.Robert E Wheeler. Quantile estimators of Johnson curve parameters.
Biometrika , pages 725–728, 1980.Marc Yor, editor.
Exponential Functionals and Principal Values related to Brownian Motion.A collection of research papers . Biblioteca de la Revista Matematic´a Iberoamericana, 1997.Marc Yor.
Exponential functionals of Brownian motion and related processes . Springer Science& Business Media, 2012.
Appendix A. Proof of Proposition 1
Proposition 1 ( Bougerol’s identity in hyperbolic geometry ) . Let X t and Z t be twoindependent BMs and the function φ defined by φ ( Z, D ) = e Z/ √ D − Z for Z ≤ D, (9) then the following is equal in distribution, conditional on Z [ µ ] T : (cid:90) T e Z [ µ ] t dX t d = X A [ µ ] T d = cos θ φ (cid:18) Z [ µ ] T , (cid:113) R T + ( Z [ µ ] T ) (cid:19) , (10) where R t is a two-dimensional Bessel process, that is, the radius of a BM in two-dimensionalEuclidean geometry, and θ is a uniformly distributed random angle. The three random variables R T , Z T , and θ are independent.Proof. Let ( x t , y t , z t ) be a three-dimensional hyperbolic BM in H with a drift on z -axis, startingat ( x , y , z ) = (0 , , dx t = z t dX t , dy t = z t dY t , and dz t z t = dZ t + (cid:18)
12 + µ (cid:19) dt, where the drift µ will be replaced by ( λ − / H introduced in Table 1 corresponds to µ = λ = −
1. Evidently, x T d = y T d = X A [ µ ] T and z T = exp (cid:16) Z [ µ ] T (cid:17) . If we let D t be the hyperbolic distance between ( x t , y t , z t ) and the starting point (0 , , D t = acosh (cid:18) (cid:18) x t + y t z t + z t + 1 z t (cid:19)(cid:19) , the Euclidean radius of ( x t , y t ) is expressed by the function φ : r t = (cid:113) x t + y t = φ (cid:16) Z [ µ ] t , D t (cid:17) . (23)The underlying BM Z [ µ ] t can be also interpreted as the projection of ( x t , y t , z t ) on the z -axis, thatis, the signed hyperbolic distance from (0 , ,
1) to (0 , , z t ). Therefore, the restriction Z [ µ ] t ≤ D t is naturally satisfied.A critical step of the proof is to show, for a fixed time T , D T d = (cid:113) X T + Y T + ( Z [ µ ] T ) conditional on Z [ µ ] T , (24), which effectively means that the hyperbolic distance between ( x T , y T , z T ) and the startingpoint (0 , ,
1) has the same distribution as the Euclidean distance of the underlying BMs,( X T , Y T , Z [ µ ] T ), from (0 , , Z [ µ ] T . Based onthe identity, it follows that (cid:113) x T + y T d = φ (cid:18) Z [ µ ] T , (cid:113) X T + Y T + ( Z [ µ ] T ) (cid:19) , or x T d = cos θ φ (cid:18) Z [ µ ] T , (cid:113) X T + Y T + ( Z [ µ ] T ) (cid:19) , where θ is a uniformly distributed random angle.Proving Equation (24) for just one value of µ is enough because the rest follows from theGirsanov’s theorem. Deviating from the original proof, µ = − H heat kernel. From the derivative of Equation (23), r T dr T = z T sinh D T dD T , the jointPDF on r T and z T is obtained from p ( t, D ):Prob( r T ∈ dr T , z T ∈ dz T ) = p ( T, D T ) 2 πr T dr T dz T z T = 1 √ πT D T e − T ( T + D T ) dD T dz T z T . From dz T /z T = dZ [ − T , it is also seen thatProb( z T ∈ dz T ) = 1 √ πT e − T Z T · dz T z T . Therefore, the conditional probability is given asProb( r T ∈ dr T | Z T ) = Prob( r T ∈ dr T , z T ∈ dz T )Prob( z T ∈ dz T ) = D T T e − T ( D T + T − Z T ) dD T z T = D T T e − T (cid:16) D T − ( Z [ − T ) (cid:17) dD T . The probability can be interpreted as the conditional probability Prob( (cid:113) X T + Y T + ( Z [ − T ) ∈ dD T | Z T ). (cid:3) YPERBOLIC NORMAL STOCHASTIC VOLATILITY MODEL 25
Appendix B. Derivation of moments of NSVh distribution
Corollary 5.
The central moments of the canonical NSVh distribution ˜ µ n = E ( ˜ F nS ) , for ≤ n ≤ , are given as ˜ µ = ρ w λ ( w −
1) + ρ ∗ w λ −
11 + λ for w = e S ≥ , ˜ µ = ρ w λ ( w − ( w + 2) + 3 ρρ ∗ w λ (cid:18) w λ −
13 + λ − w λ −
11 + λ (cid:19) , and ˜ µ = ρ w λ ( w − ( w + 2 w + 3 w −
3) + 6 ρ ρ ∗ w λ (cid:18) w w λ −
15 + λ − w λ −
13 + λ + w λ −
11 + λ (cid:19) + 32 ρ ∗ (cid:18) − w λ w λ −
15 + λ + ( w λ + 1) w λ −
13 + λ − w λ −
11 + λ (cid:19) . (20) The central moments of the original form can be scaled as µ n = E (( F T − ¯ F T ) n ) = ( σ /α ) n ˜ µ n ,and the skewness and ex-kurtosis are given as s = ˜ µ / ˜ µ / and κ = ˜ µ / ˜ µ − , respectively. Forthe normal SABR ( λ = 0 ), further simplified expressions are obtained: ˜ µ = w − , s = ρ ( w +2) √ w − , κ = ( w − (cid:18)(cid:18) ρ + 15 (cid:19) ( w + 3 w + 6 w + 5) + 1 (cid:19) . (21) Proof.
We first compute the moments, conditional on Z [ µ ] T : E (cid:16) X nA [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) = E (cos n θ ) E (cid:18) φ n (cid:18) Z [ µ ] T , (cid:113) X T + Y T + ( Z [ µ ] T ) (cid:19)(cid:19) = (2 n − n ! T n e nu √ T (cid:90) ∞ u re − r (cid:16) cosh( r √ T ) − cosh( u √ T ) (cid:17) n dr, where u = | Z [ µ ] T | / √ T and (2 n − n − n − · · · ·
1. This formula is comparable tothe formula for E (cid:16)(cid:0) A [ µ ] T (cid:1) n (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) given in Proposition 5.3 in Matsumoto and Yor [2005a]. Thefirst two values are computed in closed form: E (cid:16) X A [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) = T e u √ T m ( u, √ T ) E (cid:16) X A [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) = 3 T e u √ T (cid:16) m ( u, √ T ) − cosh( u √ T ) m ( u, √ T ) (cid:17) where m ( u, (cid:15) ) = N ( u + (cid:15) ) − N ( u − (cid:15) )2 (cid:15) e − (cid:15) n ( u ) . From these, the first two conditional moments of A [ µ ] T are trivially obtained as E (cid:16) A [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) = E (cid:16) X A [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) and E (cid:16)(cid:0) A [ µ ] T (cid:1) (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) = 13 E (cid:16) X A [ µ ] T (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) . The same results for λ = 0 are derived in Kennedy et al. [2012] in the context of the normalSABR model. The unconditional moments of X A [ µ ] S are given as E (cid:16) X A [ µ ] S (cid:17) = w µ −
12 + 2 µE (cid:16) X A [ µ ] S (cid:17) = 32 (cid:18) − w µ w µ −
16 + 2 µ + ( w µ + 1) w µ −
14 + 2 µ − w µ −
12 + 2 µ (cid:19) , where w = e S . Subsequently, the expressions for the moments follow from E (cid:16) ˜ F S (cid:17) = E (cid:16) ρ ( e Z [ µ ] S − e λS/ ) + ρ ∗ E (cid:16) X A [ µ ] S (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) (cid:17) E (cid:16) ˜ F S (cid:17) = E (cid:16) ρ ( e Z [ µ ] S − e λS/ ) + 3 ρρ ∗ ( e Z [ µ ] S − e λS/ ) E (cid:16) X A [ µ ] S (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) (cid:17) E (cid:16) ˜ F S (cid:17) = E (cid:16) ρ ( e Z [ µ ] S − e λS/ ) + 6 ρ ρ ∗ ( e Z [ µ ] S − e λS/ ) E (cid:16) X A [ µ ] S (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) + ρ ∗ E (cid:16) X A [ µ ] S (cid:12)(cid:12)(cid:12) Z [ µ ] T (cid:17) (cid:17) and the substitution µ = ( λ − /2.