Next Article in Journal
A Fuzzy Imperfect Production Inventory Model Based on Fuzzy Differential and Fuzzy Integral Method
Next Article in Special Issue
The Pricing Model of Pension Benefit Guaranty Corporation Insurance with Regime-Switching Processes
Previous Article in Journal
The ESG Disclosure and the Financial Performance of Norwegian Listed Firms
Previous Article in Special Issue
Portfolio Optimization on Multivariate Regime-Switching GARCH Model with Normal Tempered Stable Innovation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Generalized Gamma Distribution as a Useful RND under Heston’s Stochastic Volatility Model

Department of Mathematical Sciences, Indiana University—Purdue University Indianapolis (IUPUI), Indianapolis, IN 46202, USA
J. Risk Financial Manag. 2022, 15(6), 238; https://doi.org/10.3390/jrfm15060238
Submission received: 1 May 2022 / Revised: 20 May 2022 / Accepted: 21 May 2022 / Published: 26 May 2022
(This article belongs to the Special Issue Mathematical and Empirical Finance)

Abstract

:
We present the Generalized Gamma (GG) distribution as a possible risk neutral distribution (RND) for modeling European options prices under Heston’s stochastic volatility (SV) model. We demonstrate that under a particular reparametrization, this distribution, which is a member of the scale-parameter family of distributions with the mean being the forward spot price, satisfies Heston’s solution and hence could be used for the direct risk-neutral valuation of the option price under Heston’s SV model. Indeed, this distribution is especially useful in situations in which the spot’s price follows a negatively skewed distribution for which Black–Scholes-based (i.e., the log-normal distribution) modeling is largely inapt. We illustrate the applicability of the GG distribution as an RND by modeling market option data on three large market-index exchange-traded funds (ETF), namely the SPY, IWM and QQQ as well as on the TLT (an ETF that tracks an index of long-term US Treasury bonds). As of the writing of this paper (August 2021), the option chain of each of the three market-index ETFs shows a pronounced skew of their volatility ‘smile’, which indicates a likely distortion in the Black–Scholes modeling of such option data. Reflective of entirely different market expectations, this distortion in the volatility ‘smile’ appears not to exist in the TLT option data. We provide a thorough modeling of the option data we have on each ETF (with the 15 October 2021 expiration) based on the GG distribution and compare it to the option pricing and RND modeling obtained directly from a well-calibrated Heston’s SV model (both theoretically and also empirically, using Monte Carlo simulations of the spot’s price). All three market-index ETFs exhibited negatively skewed distributions, which are well-matched with those derived under the GG distribution as RND. The inadequacy of the Black–Scholes modeling in such instances, which involves negatively skewed distribution, is further illustrated by its impact on the hedging factor, delta, and the immediate implications to the retail trader. Similarly, the closely related Inverse Generalized Gamma distribution (IGG) is also proposed as a possible RND for Heston’s SV model in situations involving positively skewed distribution. In all, utilizing the Generalized Gamma distributions as possible RNDs for direct option valuations under the Heston’s SV is seen as particularly useful to the retail traders who do not have the numerical tools or the know-how to fine-calibrate this SV model.

1. Introduction

One of the most widely celebrated option pricing models for equities (and beyond) is that of Black and Scholes (1973) and of Merton (1973) (abbreviated here as the BSM model). Their option pricing model was derived under some simple assumptions concerning the distribution of the asset’s returns, coupled with presumptive continuous hedging, self-financing, zero dividend, risk-free interest rate, r, and no cost of carry or transactions fees. In its standard form, the BSM model assumes that the spot’s price process S = { S t , t 0 } evolves with a constant volatility of the spot’s returns, σ , as a geometric Brownian motion (under a risk-neutral probability measure Q , say), leading to an exact solution for the price, C ( · ) , of an European call option. Specifically, given the current spot price S τ = S and the risk-free interest rate r, the price of the corresponding call option with price-strike K and duration T,
C S ( K ) = S Φ ( d 1 ) K e r t Φ ( d 2 ) ,
where t = T τ is the lremaining time to expiry. Here, we use the conventional notation to denote by Φ ( · ) and ϕ ( · ) the standard normal cumulative distribution function ( c d f ) and density function ( p d f ), respectively, and where
d 1 : = log ( K S ) + ( r + σ 2 2 ) t σ t and d 2 : = d 1 σ t .
Despite its wide acceptability in the retail trading world1, this model hinges on several incorrect assumptions and hence suffers from some notable deficiencies; see for example Black (1989) who pointed out ‘the holes in Black–Scholes’. Chief among these deficiencies is the fact that the volatility of a spot’s returns (i.e., σ ) appears not to be constant over the ‘life’ of the option but rather varying at random.
The efforts to incorporate a non-constant volatility term in the option valuation (e.g., Wiggins 1987 or Stein and Stein 1991) has culminated with stochastic volatility (SV) model introduced by Heston (1993)  (see (A1)). This SV model incorporates, aside from the dynamics of the spot’s price process S , also the dynamics of a corresponding, though unobservable (hence untradeable), volatility process V = { V t , t 0 } . Instructed by the form of the exact BSM solution in (1), Heston (1993) obtained that the solution to the system of PDE he obtained from the stochastic volatility model he constructed is given by
C S ( K ) = S P 1 K e r t P 2 ,
where P j j = 1 , 2 are two related (under Q ) conditional probabilities that the option will expire in-the-money, conditional on the given current stock price S τ = S and the current volatility, V τ = V 0 . However, unlike the explicit BSM solution in (1) which is given in terms of the normal (or log-normal) distribution, Heston (1993) provided (semi) closed-form solutions to these two probabilities, P 1 and P 2 , which were both given in terms of their characteristic functions (c.f.); see (A2). These characteristic functions depend on some parameters of the SV model, ϑ = ( κ , θ , η , ρ ) , and they may be evaluated numerically for any choice of the parameters ϑ , in addition to the given S , V 0 and r (for more details, see Appendix A). The components of ϑ have particular meaning in the context of Heston’s SV model: ρ is the correlation between the random components of the spot’s price and volatility processes, θ is the long-run average volatility, κ is the mean-reversion speed for the volatility dynamics and η 2 is the variance of the volatility V. It should be noted that different choices of ϑ will lead to different values C S ( K ) in (3) and hence, the value ϑ = ( κ , θ , η , ρ ) must be appropriately ‘calibrated’ first for C S ( K ) to actually match the option market data. However, this calibration process typically involves substantial numerical challenges (largely resulting from numerical issues involved in the required multi-dimensional optimization, see for example Bin 2007, or Section 2.1 in Romo and Ortiz-Gracia 2021). These challenges are an obvious hindrance to the retail option traders who do not have the numerical tools or the know-how to finely calibrate the Heston (1993) SV model, as needed in the evaluation of C S ( K ) in (3).
On the other hand, as was established by Cox and Ross (1976), the risk-neutral option valuation (under Q ) provides that for T > τ (with t = T τ ), the option price C S ( K ) must also satisfy
C S ( K ) = e r t K ( S T K ) q ( S T ) d S T , = e r t K S T q ( S T ) d S T K e r t K q ( S T ) d S T e r t K S T q ( S T ) d S T K e r t · ( 1 Q ( K ) ) ,
where q ( · ) is the density of some risk-neutral distribution (RND) Q ( · ) , under the probability Q , reflective of the conditional distribution of the spot price S T at time T, given the spot price, S τ at time τ < T , whose expected value is the future value of the spot’s price. Namely, the RND q ( · ) must also satisfy,
E ( S T | S τ = S ) = S T · q ( S T ) d S T = S · e r t .
This risk-neutral density (or distribution) links together for the option valuation (under Q ) the distribution of the spot’s price S T and the stochastic dynamics governing the underlying model. As was mentioned earlier, in the case of the BSM model in (2), the RND is unique and is given by the log-normal distribution. However, since Heston’s SV model involves the dynamics of two stochastic processes, one of which (the volatility process, V ) is untradeable and hence not directly observable, there are innumerable many possible choices of RNDs, q ( · ) , that would satisfy (4) and (5), and hence, the general solutions of Heston’s P 1 an P 2 in (3) are as given by means of their characteristic functions (A2), per each possible choice of the structural parameter ϑ = ( κ , θ , η , ρ ) .

1.1. Heston’s RND as a Class of Scale-Parameter Distributions

In the literature, one can find numerous papers dealing with the ‘extraction’, ‘recovery’, ‘estimation’ or ‘approximation’, in parametric or non-parametric frameworks, or even within an entropy-based pricing framework (e.g., Yu 2020) of the RND, q ( · ) from the available (market) option prices. Some comprehensive literature reviews of the subject can be found in Jackwerth (2004); Figlewski (2010); Girth and Krätschmer (2012) and Figlewski (2018). In particular, within the parametric approach, one attempts to estimate by various standard means (maximum likelihood, method of moments, least squares, etc.) the parameters of some assumed distribution so as to approximate available option data or implied volatility (cf. Jackwerth and Rubinstein 1996). This type of assumed multi-parameter distributions includes some mixtures of log-normal distribution Mizrach (2010); Girth and Krätschmer (2012), generalized gamma Girth and Krätschmer (2012), generalized extreme value Figlewski (2010), the gamma and the Weibull distributions (Savickas 2005), among others. While empirical considerations have often led to suggesting these particular parametric distributions as possible p d f in (4), the motivation for these considerations did not include any direct link to the governing pricing model and its dynamics especially in the case of Heston’s SV model in (A1). As was mentioned earlier, in the case in the BSM model, linking directly the log-normal distribution and the price dynamics reflected by the geometric Brownian motion led to the BSM formula in (1).
In the case of the Heston (1993) SV model, this direct link between the governing price and volatility dynamics of ( S , V ) in (A1) and the class of distributions that could serve as RNDs for it has been established in Boukai (2021).
Let Δ ( K ) be the so-called delta function (or hedging fraction) in the option valuation, as defined by
Δ ( K ) = C S ( K ) S .
It is well-known that for Heston’s solution for the call option price C S ( K ) , P 1 Δ ( K ) in (3); see for example Bakshi et al. (1997) or Boukai (2021). Hence, accounting also for (4), it follows that under the SV model (A1), Heston’s solution for the option price in (3) can be written in an equivalent form as
C S ( K ) S · Δ ( K ) K e r t · ( 1 Q ( K ) ) ,
It was shown in Boukai (2021) that the general class of scale-parameter distributions that satisfies Assumption 1 below, with the mean being the forward spot’s price, μ : = S · e r t would admit the presentation in (7) and hence would satisfy Heston’s solution for the option price in (3). In fact, it was also shown that the actual RNDs (see (A4)) that may be calculated directly from Heston’s characteristic functions (A2), of to P 1 and P 2 (see Appendix A) are members of this class of distributions as well. Accordingly, the main results of Boukai (2021) (as are summarized in Theorem 2 there) establish the direct link through the solution of Heston (1993) in (3) (or (7)) between this class of RNDs and the assumed stochastic volatility model governing the spot price and volatility dynamics.
To fix ideas, we set μ = S · e r t to denote the forward spot’s price, and correspondingly, we denote by Q μ ( · ) the RND with a corresponding p d f q μ ( · ) as in (4) and (5).
Assumption 1.
It is assumed that μ is a scale parameter of Q μ ( · ) so that for any x > 0 , Q μ ( x ) Q 1 ( x / μ ) and q μ ( x ) q 1 ( x / μ ) / μ for some c d f Q 1 ( · ) with a p d f q 1 ( · ) satisfying 0 x q 1 ( x ) d x = 1 and 0 x 2 q 1 ( x ) d x = 1 + ν 2 . Here, ν 2 is some exogenous parameter (to be specified later).
Note that if Q μ ( · ) satisfies Assumption 1, then, for any k > 0 ,
c μ ( k ) : = k ( x k ) q μ ( x ) d x = k ( 1 Q μ ( x ) ) d x = k ( 1 Q 1 ( x / μ ) ) d x μ c 1 ( k / μ ) .
Hence, c μ ( k ) is a homogeneous function2 of degree one in both μ and k, so that for k = a k and μ = a μ with a > 0 , we have c μ ( k ) a c μ ( k ) . This property and (8) are the key elements in the proof of Theorem 1 of Boukai (2021), which establishes that for this class of scale-parameter distributions satisfying Assumption 1, the delta function (6) may in fact be written as
Δ ( K ) Δ μ ( K ) : = 1 μ K x q μ ( x ) d x Δ 1 ( K / μ ) ,
where Δ 1 ( a ) : = a u q 1 ( u ) d u for any a > 0 . Accordingly, for any member of this scale parameter (in μ = S · e r t ) class of distributions defined by Q 1 ( · ) as in Assumption 1, the option price C S ( K ) in (7) may equivalently be written as
C S ( K ) S · Δ 1 ( K / μ ) K e r t · ( 1 Q 1 ( K / μ ) ) ,
and thus, it could be used for the direct risk-neutral valuation, as RND, of the option price under Heston’s SV model. The expression in (10) is our ‘working’ formula for the direct calculations of the option price C S ( K ) in the case of scale-parameter distribution defined by Q 1 ( · ) . This result was illustrated in great detail by Boukai (2021) for several well-known parametric (scale) distributions which under a particular parametrization satisfy Assumption 1. These distributions include one-parameter versions of the log-Normal (i.e., the BSM model), Inverse-Gaussian, Gamma, Inverse-Gamma, Weibull and the Inverse-Weibull distributions, all which provide explicit RNDs for Heston’s pricing model in various market circumstances (e.g., negatively skewed RND to match SPX option data (i.e., Bakshi et al. 1997), or ODAX option data (i.e., Mrá¡zek and Pospíšil 2017), or positively skewed RND to match AMD (say), option data).
Remark 1.
In the case in which the risk-neutral valuation of the option includes a dividend with a rate ℓ, then E ( S T | S ) = S e ( r ) t in (5), in which case, by applying μ = S e ( r ) t to (7), we obtain from (8) and (10),
C S ( K ) = e r t c μ ( K ) = S e t Δ 1 ( K / μ ) K e r t ( 1 Q 1 ( K / μ ) ) .

1.2. An Overview

In this paper, we focus attention on a two-parameter version of the Generalized Gamma (GG) distribution as is especially parametrized to satisfy Assumption 1 and hence, to serve as an RND under Heston’s SV option valuation model. The particular version of this distribution we consider here is characterized by two shape parameters α and ξ say, and it is general enough to admit either positively skewed distributions ( ξ < 0 ) or negatively skewed distributions ( ξ > 0 ) . Aside from this noted ‘elasticity’ to match well the varying characteristics of different spot’s RNDs under the SV model (A1) implied from different market scenarios (i.e., with different ϑ = ( κ , θ , η , ρ ) in (A1)), this distribution is especially useful in modeling option prices in situations that exhibit put-over-call skew and and hence admit negatively skewed distribution of the spot’s price, indeed with ξ > 0 .
In Section 2, we present the GG distribution and its required reparametrization as a RND for the SV model (3). Though not of immediate interest, we also present in Section 2.2 the case with ξ < 0 (the so-called Inverse Generalized Gamma (IGG) distribution) as a possible RND under Heston’s SV model that could be useful in modeling positively skewed (implied) distributions.
In Section 3, we apply the GG distribution (with ξ > 0 ) as RND to modeling current3 market option data on three large market-index ETFs, namely the SPY, IWM and QQQ as well as on the TLT (a large ETF that tracks an index of long-term US Treasury bonds). The current option chain of each of the three market ETFs exhibits a pronounced skew of their volatility ‘smile’, which indicates a likely distortion in the Black–Scholes modeling of such option data. Reflective of entirely different market expectations, this distortion appears not to exist in the TLT option data (see Figure 1 below). We provide a thorough modeling of the available option data we have on each ETF (with the 15 October 2021 with 63 days to expiration) based on the GG distribution (with ξ > 0 ) and compare it to the option pricing and RND modeling obtained directly from a well-calibrated Heston (1993) SV model (both theoretically and empirically, using Monte Carlo simulations of the spot’s price). All three market-index ETFs exhibit negatively skewed distributions, which are well-matched with those derived under the GG distribution as RND. The inadequacy of the classical Black–Scholes modeling in such instances which involve negatively skewed implied distribution is further illustrated by its impact on the hedging factor, delta, and the immediate implications to the retail trader. In contrast, for the TLT ETF, which exhibits no such distortion to the volatility ‘smile’, the three pricing models (i.e., Heston’s, Black–Scholes and Generalized Gamma) appear to yield very similar results. Technical notes are provided in Section 3.2, and some details on Heston’s SV model and related cf. that are used in the calculation of (3) are provided in Appendix A.

2. The Generalized Gamma Distribution as an RND for Heston’s SV Model

Introduced by Stacy (1962), the Generalized Gamma (GG) distribution is demonstrably highly versatile, with a vast number applications, from survival analysis to meteorology and beyond. It includes among many others the Weibull distribution ( α = 1 ), the Gamma distribution ( ξ = 1 ), and also the log-normal distribution as a limiting case, ( α ). In this section, we show that this GG distribution along with its counterpart, the so-called Inverse Generalized Gamma distribution (IGG), both satisfy under a particular reparametrization, Assumption 1, and hence could serve as RND (for direct option valuation using (10)) under the Heston (1993) SV model for option valuation. Though similar, we will present these two cases of the Generalized Gamma distribution separately as they do present different profiles of skewness and kurtosis. We will however focus our attention on the GG distribution (with ξ > 0 ), as we will use it for option pricing modeling in a situation that involved negatively skewed (implied) risk-neutral distributions.
We begin with some standard notations. We write Y G ( α , λ ) to indicate that the random variable Y has the gamma distribution with a scale parameter λ > 0 and a shape parameter α > 0 (so that its mean is E ( Y ) = α / λ ). We write g ( · ; α , λ ) and G ( · ; α , λ ) for the corresponding p d f and c d f of Y, respectively,
g ( y ; α , λ ) λ α y α 1 e λ y Γ ( α ) and G ( y ; α , λ ) Γ ( y λ ; α ) Γ ( α ) ,
where Γ ( α ) : = 0 x α 1 e x d x denotes the gamma function whose incomplete version is Γ ( s ; α ) : = 0 s x α 1 e x d x , is defined for any s > 0 .

2.1. The GG Distribution

The Generalized Gamma (GG) distribution is typically characterized by three parameters: a scale parameter, λ > 0 , and two shape parameters, α > 0 and ξ > 0 , and it is defined as follows. We say that W GG ( λ , ξ , α ) , if
Y W λ ξ G ( α , 1 ) .
In light of relation (12), the c d f and p d f of W GG ( λ , ξ , α ) are readily available in terms of the Gamma distribution in (11). More specifically, for any w > 0 ,
F W ( w ) : = P r ( W w ) = G ( w λ ) ξ ; α , 1 ,
and
f W ( w ) = ξ λ ( w λ ) ξ 1 · g ( w λ ) ξ ; α , 1 .
In addition, the j th , j = 0 , 1 , 2 , , moment of this distribution (see Stacy and Mihram 1965), whenever it exists, (i.e., whenever α + j / ξ > 0 ) is given by
E ( W j ) = λ j Γ ( α + j / ξ ) Γ ( α ) : = λ j · h j ( ξ ) , with α + j / ξ > 0 , and α > 0 .
Now, suppose that for a given α > 0 , a random variable U has the ’standardized’ version of the GG distribution, with mean E ( U ) = 1 and a variance V a r ( U ) = ν 2 , for some ν > 0 (in fact, we will later take ν = σ t for some σ > 0 ). Utilizing h j ( ξ ) as defined in (13), we let for a given α > 0 and ν > 0 , ξ * ξ ( ν ) be the (unique) solution of the equation
h 2 ( ξ ) h 1 2 ( ξ ) = 1 + ν 2 ,
in which case, h j * h j ( ξ * ) , j = 1 , 2 , λ * 1 / h 1 * and U GG ( λ * , ξ * , α ) . Accordingly, the c d f of U is given by
Q 1 ( u ) : = P r ( U u ) = G ( u λ * ) ξ * ; α , 1 ,
for any u > 0 , and its p d f is given by,
q 1 ( u ) : = ξ * λ * ( u λ * ) ξ * 1 · g ( u λ * ) ξ * ; α , 1 , u > 0 .
It follows immediately from (15) that if X μ · U for some μ > 0 , then the p d f , q μ ( · ) of X is the ’scaled’ version of q 1 ( · ) above. For this RND, the values of Δ 1 ( s ) in (9) can be calculated using (16) for any a > 0 as,
Δ 1 ( a ) = a u q 1 ( u ) d u = 1 G ( ( a / λ * ) ξ * ; α + 1 / ξ * , 1 ) ,
which, when combined in (10) with the expression of Q 1 ( · ) as given in (15) above, provides the values of
c μ ( k ) = μ × Δ 1 ( k / μ ) k μ × ( 1 Q 1 ( k / μ ) ) , = μ × 1 G ( ( k / μ λ * ) ξ * ; α + 1 / ξ * , 1 ) k × 1 G ( ( k / μ λ * ) ξ * ; α , 1 )
for any μ > 0 . Finally, to calculate under this Generalized Gamma RND the price of a call option at a strike K when the current price of the spot is S, we will utilize (18) with μ S e r t (being the forward price), k K and with λ * 1 / h 1 ( ξ * ) and ξ * ξ ( ν ) as is determined by Equation (14) above with ν σ t to obtain C S ( K ) = e r t c μ ( K ) as
C S ( K ) = S · 1 G ( d ; α + 1 / ξ * , 1 ) K e r t · 1 G ( d ; α , 1 ) ,
where
d = K e r t h 1 ( ξ * ) S ξ * , with ξ * ξ ( ν ) from ( 14 ) .
We point out that for given current spot’s price, S, a strike price K, risk-free interest rate, r, and the remaining option’s duration t, the option value C S ( K ) in (19) involves, through Equation (14) (with ν σ t ), with only two parameters, namely α and σ . Their values can easily be “calibrated” from the available market option data. Indeed, in the Generalized Gamma case, this calibration task is computationally much simpler than the direct calibration of four parameters, ϑ = ( κ , θ , η , ρ ) , of Heston’s pricing model, based on the characteristic functions (see (A2) in Appendix A), which also involves integration over the complex domain.

2.2. The IGG Distribution

For the sake of completeness, we also present the details of this variant to the Generalized Gamma distribution here as well. With some additional restrictions on ξ , one can similarly define the Inverse Generalized Gamma distribution (IGG). Namely, we say that W IGG ( λ , ξ , α ) , if
Y W λ ξ G ( α , 1 ) .
The option pricing model under the Inverse Generalized Gamma distribution as RND for the Heston’s SV for option valuation is constructed similarly to that of the GG in the previous section. By relation (20), if W IGG ( λ , ξ , α ) , then its c d f is given, for w > 0 ,
F W ( w ) : = P r ( W w ) = 1 G ( w λ ) ξ ; α , 1 .
In this case, too, the ’standardized’ IGG distribution of U is constrained to have mean 1 and variance ν 2 , which requires a restriction on the parameter ξ > 2 / α . It follows that with such a restriction, U IGG ( λ * , ξ * , α ) , but now, ξ * ξ ( ν ) is the (unique) solution of the equation
h ˜ 2 ( ξ ) h ˜ 1 2 ( ξ ) = 1 + ν 2 ,
where, h ˜ j ( ξ ) h j ( ξ ) = Γ ( α j / ξ ) / Γ ( α ) , j = 1 , 2 , provided that α j / ξ > 0 , in which case, h ˜ j * h ˜ j ( ξ * ) , j = 1 , 2 , λ * 1 / h ˜ 1 * . Accordingly, the c d f of U is given by
Q 1 ( u ) : = P r ( U u ) = 1 G ( u λ * ) ξ * ; α , 1 ,
for any u > 0 , and in similarity to (17), its corresponding delta function is given by
Δ 1 ( s ) = G ( ( s / λ * ) ξ * ; α 1 / ξ * , 1 ) ,
Again, by combining (22) and (23) in (10), we obtain that for any μ > 0 ,
c μ ( k ) = μ × G ( ( k / μ λ * ) ξ * ; α 1 / ξ * , 1 ) k × G ( ( k / μ λ * ) ξ * ; α , 1 ) .
Accordingly, in order to calculate under this Inverse Generalized Gamma RND the price of a call option at a strike K when the current price of the spot is S, we will utilize (24) with μ S e r t , k K and with λ * 1 / h ˜ 1 ( ξ * ) and ξ * ξ ( ν ) as is determined by Equation (21) above with ν σ t to obtain, C S ( K ) = e r t c μ ( K ) as,
C S ( K ) = S · G ( d ; α 1 / ξ * , 1 ) K e r t · G ( d ; α , 1 ) ,
where
d = K e r t h ˜ 1 ( ξ * ) S ξ * , with ξ * ξ ( ν ) from ( 21 ) .

2.3. Skew and Kurtosis

As can be seen from the above construction of the RNDs, both the GG and IGG distributions depend on two shape parameters ( α , ξ * ) , or equivalently ( α , ν ) , where ν σ t , that affect their features, such as kurtosis and skewness, and hence their suitability as RND for various particular scenarios of the SV model (A1), as is determined by the structural model parameter ϑ = ( κ , θ , η , ρ ) (more on this point in the next section). Unlike the standardized log-normal distribution which has a positive skew only, these two classes of distributions offer a range of RNDs with positive as well as negative skewness. This is a critical feature to have when modeling option prices for characteristically different spots such as an Index (SPX, say) as opposed to modeling option prices for a technology firm (such as AMD, say).
For a given ( α , ξ * ) , we denote these two measures as γ 1 ( ξ * ) for skew and γ 2 ( ξ * ) for the kurtosis. Then, with h j ( ξ ) : = Γ ( α + j / ξ ) , we have in the GG case that with ξ * = ξ ( ν ) which satisfies (14),
γ 1 ( ξ * ) = h 3 ( ξ * ) 3 ν 2 1 ν 3
and
γ 2 ( ξ * ) = h 4 ( ξ * ) 4 ν 3 γ 1 ( ξ * ) 6 ν 2 1 ν 4 .
For the IGG case, these two measure are similar and are given by γ 1 ( ξ * ) and γ 2 ( ξ * ) , provided that ξ * = ξ ( ν ) as is determined by (21) satisfies that ξ * > 4 / α .

3. Calibration, Validation and Examples

3.1. Observing the Skew

In this section, we demonstrate the usefulness of the Generalized Gamma distribution to serve as an RND under Heston’s Stochastic Volatility model in cases that exhibit a high put–call skew (i.e., OTM puts in the option series are far more expensive than equidistant OTM calls) and hence expressing a pronounced skew in the so-called “volatility smile” of the series. Cases in point are traded market indexes such as the S&P 500 (SPX), Russel 2000 (RUT) or Nasdaq 100 (NDX), which all are (along with their corresponding ETF surrogates, SPY, IWM and QQQ) currently at (or near) their all time high levels4. Market expectations of an eminent ‘correction’ are often seen as the culprits that affected the implied volatility surface associated with the corresponding option series of the index (see for example Bakshi et al. 1997).
Figure 1 below displays the calculated implied volatility smile of the 15 October 2021 option series for these three ETFs, SPY, IWM and QQQ, as quoted on 13 August 2021, (EOD), each with 63 days to expiration (DTE). Several days later, on 18 August 2021, we obtained the corresponding quote for the TLT, but now with 57 DTE. For each ETF, the EOD option’s market prices (for puts and calls) at the corresponding strikes were recorded along with the BSM-based calculated delta and implied volatility as provided by the brokerage firm.5 As a reference, we also marked on these plots (in red) the current spot’s (ETF) price S along with the ATM (BSM-based) calculated implied volatility (IV) for each ETF. As can be seen from these figures, the options of the three market index ETFs exhibit a highly pronounced skew in their volatility ‘smile’, whereas the option on the TLT ETF does not (likely only reflective of market’s expectations of actions by the Federal Reserve).
However, since typically in the retail world, the calculated option’s implied volatility (as well as other associated quantities, such as the option’s delta) is calculated based on the Black–Scholes formula in (1), the noted distortion in the volatility smile (or surface) is nonetheless also indicative that the assumed underlying log-normal distribution of the Black–Scholes model (with its distinctive positive skew measure) is a poor choice to serve as RND in such instances involving a stochastic volatility structure as that of Heston (1993) (see (A1) below), particularly in those instances that admit a negatively skewed RND. To illustrate the extent of the “inaptness” of using the log-normal distribution as RND (the BSM formula in (1)) for the option valuation in such skewed cases, we have calibrated for each of these four ETFs the appropriate Heston’s SV model to fit the observed market option data (i.e., on 15 October 2021 option series for each) and derived from it the implied RND of Heston’s model (HS). This RND, which is obtained both theoretically, using (A2) and (3), and also via Monte Carlo simulations of (A1), will serve as a benchmark for comparison.
For each option series, the available market data consist of the N strikes, K 1 , , K N with corresponding call option (market) prices C 1 , , C N 6. As a standard measure of the goodness-of-fit between the model-calculated option prices C M o d e l ( K i ) , i = 1 , , N and the given option market price C i , i = 1 , , N , we used the Mean Squared Error, MSE,
M S E ( M o d e l ) = 1 N i = 1 N ( C M o d e l ( K i ) C i ) 2 .
Clearly, it is expected that within the scope of the SV model described in (A1), the well-calibrated Heston model will result with a smaller MSE as compared to the Black–Scholes model, so that M S E ( H S ) M S E ( B S ) . However, as we will see below for the available ETF data, pricing the options by a well-calibrated Generalized Gamma (GG) model (19) also resulted with a smaller MSE. In fact, in all four cases, M S E ( G G ) M S E ( B S ) . To demonstrate this, we have taken for each ETF the following steps (conditional of course on the current spot’s price S and volatility V 0 ):
  • Model Calibration
    -
    For a given model’s parameter, ϑ = ( κ , θ , η , ρ ) in (A1), we use the callHestoncf function of the NMOF package (see Gilli et al. 2019 and Schumann 2011–2021) and the R software (R Core Team 2017) to calculate the Heston model’s option prices C i H S for each K i .
    -
    To calibrate the Heston SV model, we used the optim(·) function R to minimize M S E ( H S ) over the model’s parameter, ϑ = ( κ , θ , η , ρ ) .
    -
    For a given ( α , ν ) with ν = σ t , we use (19) to calculate the Generalized Gamma model option prices C i G G for each K i .
    -
    To calibrate the GG model, we used the optim(·) function of R to minimize M S E ( G G ) over the model’s parameters, ( α , ν ) .
    -
    For a given ν (where ν = σ t ), we use (1) and (2) to calculate the Black–Scholes model option prices C i B S for each K i .
    -
    To calibrate the BS model, we used the optimize(·) function of R to minimize M S E ( B S ) over the single model’s parameter ν (namely σ ).
  • Validation
    -
    Using the calibrated Heston parameters, ϑ ^ , we drew, utilizing a discretized version of Heston’s stochastic volatility process (A1), a large number (M = 30,000) of Monte Carlo simulations, observations on ( S T , V T ) to obtain the simulated rendition of the Heston’s RND of S t (conditional on S and V 0 , with t = T τ ).
    -
    Using the calibrated Heston’s parameters, ϑ ^ , in (A4), we obtain the calculated rendition of the Heston’s theoretical RND of S t (conditional on S and V 0 , with t = T τ ) directly from the characteristics function of P 2 (see Appendix A).
    -
    Finally, we compared all three calibrated risk-neutral distributions of the standardized spot’s price (the rescaled spot priced, S t * = S t / μ , where μ = S e r t ) as obtained under the Black–Scholes (BS), Generalized Gamma (GG) and Heston (1993) option pricing models (HS).

3.2. Calculating the Implied RND under the Volatility Skew

As we mentioned earlier, the data on the 15 October 2021 option series of the SPY, IWM and QQQ were retrieved as of the closing of trading on Friday 13 August 2021 with 63 days to expiration, so that t = 63 / 365 and the prevailing (risk-free) interest rate at that time is r = 0.0016 . There will be common values for these three highly liquid ETFs. The 15 October 2021 option series of the TLT was retrieved on 18 August 2021 with 57 days to expiration, so that t = 57 / 365 for that ETF. However, we begin our exposition with the details of the largest (volume-wise) of them, namely the SPY. The cases of the IWM, QQQ and TLT will be treated similarly below.
On that day, the closing price of the SPY was S = 445.92 , and the dividend it pays is at a rate of = 0.0123 . We incorporate the dividend in our calculations along the lines of Remark 1. The reported (BS-based) implied volatility was I V = 16.15 % , which we will use as our initial value for V 0 and for σ . This option series has N = 211 pairs of strike-price ( K i , C i ) , which were all used to calibrate Heston’s SV model over the model’s parameter, ϑ = ( κ , θ , η , ρ ) , with the initial values of ( 15 , ( 0.1 ) 2 , 0.1 , 0.65 ) and with V 0 = I V 2 = ( 0.1615 ) 2 . The results of the calibrated values are
ϑ ^ = ( 15.03132587 , 0.02793781 , 2 , 0.77469470 ) .
This calibrated parameter, ϑ ^ , was then used to calculate, using Heston’s characteristic function (i.e., (A2)), the option prices according to Heston’s SV model (3). This resulted with M S E ( H S ) = 0.2226429 . The calibrated (least squares estimate) value of σ that minimizes the MSE for the BS model is σ ^ = 0.137348 , so that M S E ( B S ) = 1.781981 . Accordingly, ν ^ B S = 0.137348 t = 0.0570619 is to be used for the calculation of the p d f of the N ( ν 2 / 2 , ν 2 ) distribution, which leads to the BS formula in (1) (see Example 3.1 in Boukai 2021 for more details). Next, we calibrated the General Gamma distribution according to the pricing model in (19), with initial values of α = 0.5 and σ = 0.1615 , which resulted with calibrated value of α ^ = 0.1554312 and σ ^ = 0.1483843 and an M S E ( G G ) = 0.339441 . Clearly, in this case of the SPY, the MSE of the GG pricing model is substantially smaller than the MSE of the Black–Scholes model and is similar to that of Heston’s SV pricing model. Indeed, the MSE of the BS model is over 500% as large as those of the GG and the HS models.
To compare the actual distributions, as were calculated under each of these three pricing models, we present in Figure 2 the three implied distributions (as RNDs), which were calculated based on their respective calibrated parameter values. As an added validation, we plotted these three density curves against the histogram of the Monte Carlo simulation of the standardized SPY prices using a discretized version of the pricing model in (A1) (using the calibrated Heston parameters with a seed = 452361). This figure clearly demonstrates the ‘inaptness’ of the standard BSM Formula (1) and hence the log-normal distribution for the direct (risk-neutral) modeling of option prices in cases which involve negatively skewed price distributions. In fact, the calculated values of the kurtosis and skewness measures of each of these distributions (see Table 1) are also indicative of the noted lack-of-fit of the BS model in these cases and the apparent close agreement of the GG distribution to the exact risk-neutral distribution of the Heston’s model and that of the simulated price data.
The impact of this model’s misspecification on the calculated delta values associated with the option series is also of interest. It is a standard practice of the retail brokerage houses to provide, along with the market prices for the option chain, also the BS-base calculated delta for each strike (using some ATM implied volatility value). For example, for the ATM strike of K = 445 , the quoted delta is Δ * = 0.497 with a quoted I V of 0.1489 , whereas under the BS model we calibrated here with σ ^ = 0.137348 , we obtained Δ B S = 0.506 . However, accounting for stochastic volatility in the pricing model, we calculate for this same strike, K = 445 , Δ H S = 0.663 by the (better fitting) Heston SV model, and Δ G G = 0.638 , by its close proxy, the GG model. Thus, in this case, the BS modeling at the ATM strikes will result with grossly understated delta values (of nearly 25.0%). Without doubt, the impact of this model’s misspecification would have profound hedging implications for the retail trader. To fully appreciate the extent of this impact, we present in Figure 3 the values of the delta function (9) as was calculated for the HS model (using P 1 and (A1)), for the GG model (using (17)) and for the BS model (using Φ ( d 1 ) from (1)), along with quoted delta values for the SPY chain.
Needless to say, the noted understatement of the quoted (BS -based) delta values as compared to those derived from the SV model also impact the trading strategies. For example, a trader that would sell a 25-delta strangle based on the quoted values will sell the k 1 = 424 put for $5.215 and the k 2 = 460 call for $2.685, collecting a total of $7.90 for it, which amounts to 21.9% of the spread between the strikes (for a discussion of this ratio, see Boukai 2020). On the other hand, if the trader would have priced the 25-delta strangle according to the GG model (which accounts for the skew), she will sell the k 1 = 435 put for $7.205 and the k 2 = 466 call for $1.435, collecting a total of $8.64 for it, which amounts to 27.9% of the spread between the strikes, clearly collecting a higher premium for the same 25-delta strangle.
The situation with the other two market index ETFs, IWM and QQQ, is very similar to the one describing the SPY—see the corresponding depiction of their volatility ‘smiles’ in Figure 1. Following a similar calibration and validation approach, we present (implied) the price distributions derived from the IWM option data shown in Figure 4a and the QQQ option data shown in Figure 5a. The calculated values of the corresponding delta functions are displayed in Figure 4b and Figure 5b. Furthermore, to serve as a contrasting illustration, we present in Figure 6 the three implied price distributions derived from the TLT ETF option series, along with the corresponding calculated delta functions for that ETF. The situation with the TLT ETF is clearly different, as compared to the three market index ETFs (SPY, IWM and QQQ) which exhibit a pronounced skew of their volatility ‘smile’. In the case of the TLT) ETF, with a relatively intact volatility ‘smile’ (see Figure 1), the implied RNDs are relatively symmetric, and the three option pricing models (HS, GG and BS) yield very similar results. In Table 2, we provide a summary of the goodness-of-fit of each of the pricing models as measured by the respective MSE for each of the four ETFs. A corresponding comparison of the ATM delta calculations under each of the option price models is presented in Table 3. In Table 3, we provide a summary of the goodness-of-fit as measured by the respective MSE for each of the ETFs. Some of the technical details are provided in Section 4.1 below.

4. Summary and Discussion

As was illustrated in all the above examples, the Heston (1993) option pricing model (as given in (A1) and (3)), which accounts for the presences of stochastic volatility, produces as expected the best results overall as compared to the Black–Scholes option pricing model (1) with its presumed constant volatility. Clearly, a well-calibrated Heston model will always result in a better fit to realistic market option data (indeed, resulting with M S E ( H S ) < M S E ( B S ) ) and would be the default modeling choice for the practitioner. Unfortunately, however, the numerical challenges involved in the calculations and calibration (or optimization of ϑ = ( κ , θ , η , ρ ) ) process of the Heston’s option pricing model (see for example Romo and Ortiz-Gracia 2021 or Lemaire et al. 2020) render it largely inaccessible to many of the retail option traders who do not possess the prerequisite skills or know-how to meet these numerical challenges. In comparison, the calculations and calibration process involved with the two-parameter Generalized Gamma distribution as RND for Heston’s SV option valuation are substantially simpler and more straightforward (and could potentially be accomplished within an Excel spreadsheet). As was demonstrated earlier, the GG model is significantly more accurate than the Black–Scholes model for the pricing of the options in a skewed stochastic volatility environments as those exhibited (at present times) by the three markets ETFs, SPY, IWM and QQQ. In fact, in situations that imply negatively skewed price-distributions as RND, the Black–Scholes pricing model, and hence the log-normal distribution as RND, will surely be inferior to the GG distribution as an RND and surely to Heston’s SV pricing model in fitting realistic option market data. In such situations, one would realize (as we did in these examples) M S E ( G G ) < M S E ( B S ) and would want to adopt the GG RND for the underlying pricing model. In contrast, in situations such as the one exhibited by the TLT ETF, one would realize M S E ( G G ) M S E ( B S ) , as all three option pricing models (including Heston’s) produce similar results. Although not expressly covered by the examples we included here, we have grounds to believe that the same conclusion could be arrived upon using the Inverse Generalized Gamma (see Section 2.2) in situations involving positively skewed (implied) RND in the option pricing model. Although the two-parameter versions of the GG and IGG distribution are similar, they differ in their implied parameter space (see restrictions in (13) and Section 2.3) and thus should be treated separately. In all, both of these versions of the Generalized Gamma distribution could serve as useful proxies to the exact Heston’s RND given in (A4) and hence produce superior results to those obtained by the Black–Scholes model in an environment involving stochastic volatility. Thus, given the market option data, one could simply calculate M S E ( B S ) , M S E ( G G ) , and if appropriate also M S E ( I G G ) , and adopt the RND for the option pricing model, which produces the smaller M S E and hence the better fit to the market data.

4.1. Some Technical Notes

  • The 15 October 2021 option series data files SPY_63.csv, IWM_63.csv, and QQQ_63.csv as were obtained on the EOD of 13 August 2021 and that of TLT_57.csv obtained at the EOD of 18 August 2021 are available from the author upon request. Their basic summary information is provided in Table 4 below.
    Table 4. Summary information of the four ETFs.
    Table 4. Summary information of the four ETFs.
    ETFSDTENQuoted IVDiv. Rate
    SPY445.926321116.15%1.23%
    IWM221.13639324.30%0.63%
    QQQ368.826316018.13%0.43%
    TLT149.35576615.71%1.46%
  • The standard R function dgamma and pgamma were used to calculate the p d f and c d f in (11) and hence used in the calculation of (19); see Appendix B.
  • The cfHeston and callHestoncf functions of the NMOF package of R were used in the calculation of (A2) and (3).
  • A modification of the callHestoncf function of the NMOF package of R was used to calculate (A4).
  • The optim and optimize functions of R were used in the calibration of the three models (HS, GG and BS) for the available option data.
  • The initial and the calibrated values of ϑ = ( κ , θ , η , ρ ) of the Heston’s model were:
    • SPY:   ( 15 , ( 0.1 ) 2 , 0.1 , 0.65 ) and ( 15.03132587 , 0.02793781 , 2 , 0.77469470 ) .
    • IWM:   ( 5 , ( 0.1 ) 2 , 0.6 , 0 ) and ( 4.97834286 , 0.04032166 , 1.09837930 , 0.59905916 ) .
    • QQQ:   ( 3.5 , ( 0.2 ) 2 , 0.5 , 0.5 ) and ( 3.47635183 , 0.06382197 , 1.13505528 , 0.69137767 ) .
    • TLT:   ( 3 , ( 0.1 ) 2 , 0.1 , 0.1 ) and ( 2.99997881 , 0.01459405 , 0.10011507 , 0.10007980 ) .
  • For the Monte Carlo simulation of (A1), we employed the (reflective version of) Mil’shtein (1975) discretization scheme (see also Gatheral 2006) with seeds = 4569 (QQQ), = 777999 (IWM), = 452361 (SPY) and = 121290 (TLT).

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data are available upon request from the author.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A. Heston’s 1993 Solution

Heston (1993) considered the stochastic volatility model describing the price-volatility dynamics (of S = { S t , t 0 } and V = { V t , t 0 } ) as described via a system of stochastic deferential equations (SDE) given by
d S t = r S t d t + V t S t d W 1 , t d V t = κ ( θ V t ) + η V t d W 2 , t ,
where r is the risk-free interest rate, κ , θ and η are some constants (as discussed in Section 1) and where W 1 = { W 1 , t , t 0 } and W 2 = { W 2 , t , t 0 } are two Brownian motion processes (under the risk neutral probability Q ) with d ( W 1 W 2 ) = ρ d t for some ρ 2 ( 0 , 1 ) ). Heston (1993) offered C S ( K ) in (3) as the solution to the option valuation under the above SDE and provided (semi) closed form expressions to the probabilities P 1 and P 2 that comprise it. These closed form expressions are given for j = 1 , 2 by
P j = 1 2 + 1 π 0 R e e i ω k ψ j ( ω , t , v , x ) i ω d ω ,
where with x : = log ( S ) , k : = log ( K ) , b 1 = κ ρ η , b 2 = κ and ψ j ( · ) is the characteristics function
ψ j ( ω , t , v , x ) : = e i ω s p j ( s ) d s e B j ( ω , t ) + D j ( ω , t ) v + i ω x + i ω r t .
Here, p j ( · ) is the p d f of s T = log ( S T ) corresponding to the probability P j , j = 1 , 2 and
B j ( ω , t ) = κ θ η 2 { ( b j + d j i ω ρ η ) t 2 log ( 1 g j e d j t 1 g j ) }
D j ( ω , t ) = b j + d j i ω ρ η η 2 ( 1 e d j t 1 g j e d j t )
g j = b j i ω ρ η + d j b j i ω ρ η d j
d j = ( i ω ρ η b j ) 2 η 2 ( 2 i ω u j ω 2 ) .
Now, by a standard application of the Fourier transform, we obtain (see for example Schmelzle 2010) that the p d f p j ( · ) of s T = log ( S T ) can be computed, for any s R , as
p j ( s ) = 1 π 0 R e e i ω s ψ j ( ω , t , v , x ) d ω .
Hence, it follows immediately that the p d f p ˜ j ( · ) of S T is given, for any u > 0 , by
p ˜ j ( u ) = 1 u × p j ( log ( u ) ) 1 π 0 R e e i ω log ( u ) ψ j ( ω , t , v , x ) u d ω .
Further, since the c . f . , ψ j , above are affine in x + r t = log ( S ) + r t log ( μ ) , we may rewrite p ˜ j ( u ) as
p ˜ j ( u ) = 1 μ π 0 R e e i ω log ( u / μ ) ψ ˜ j ( ω , t , v ) u / μ d ω ,
where log ( ψ ˜ j ( ω , t , v ) ) : = log ( ψ j ( ω , t , v , x ) i ω x i ω r t . Note that p ˜ 2 ( · ) in (A4) is the p . d . f corresponding to the RND, Q μ ( · ) , of S T (under Q ) for Heston’s (1993) model with P 2 K p ˜ 2 ( u ) d u = Q ( S T > K ) Q μ ( K ) , in (3) and (7). In addition, note that this Q μ ( · ) , constitutes a scale-family of distributions in μ = S e r t , so that it satisfies the terms of Assumption 1. As was mentioned in Section 4.1, the cfHeston and callHestoncf functions of the NMOF package of R are readily available to accurately compute the values of ψ j and hence of P j (as well as p ˜ j ( u ) ) as well as the call option values C S ( K ) in (3) for given t , s and v and any given choice of ϑ = ( κ , θ , η , ρ ) .

Appendix B. R Code for the GG Model

The following is the simple R code for calculating the option call price under the GG model as given in (19).
##
# s0   	#= current spot’s price
# k    	#= strike
# te   	#= days/365
# r    	#= interest rate
# q    	#= dividend rate
# sig  	#= volatility (sigma)
# alpha #= first shape parameter, alpha
##
GG.value<−function(s0, k, te, r, q, sig, alpha) {
  nu<−sig*sqrt(te)
  k<−k*exp(−r*te)
  s0<−s0*exp(−q*te)
  s1<−k/s0
  f0<−GG.call(s1, nu, alpha)
  return(f0)
}
###
GG.call<−function(s, v, alpha){
  xi<−seq(0.1, 100, length=10000)
  yy<− (gamma(alpha)*gamma(alpha+2/xi))/(gamma(alpha+1/xi))^2
  xi0<−min(xi[yy<1+v^2])       # second shape parameter
  lam0<−gamma(alpha)/gamma(a)  	# the scale parameter
    s1<−(s/lam0)^(xi0)
    delta<−1−pgamma(s1, alpha+1/xi0, 1)
    prob<−1−pgamma(s1, alpha,1)
    cc0<−delta−s*prob
    f0<−cbind(s, cc0, delta, prob, xi0, lam0)
  return(f0)
}
##

Notes

1
Nowadays, many of the retail brokerage houses operate entirely within the ‘Black–Scholes world’ and provide, aside from market option bid–ask prices, the ‘theoretical price’ and related ‘Greeks’, and implied volatility values as are derived from and calculated under the BSM Formula (1) and (2).
2
Note that c μ ( · ) is merely the undiscounted version of C S ( · ) in (4). For the linear homogeneity property of the European options, in general, see for example Theorems 6 & 9 of Merton (1973).
3
As of the original draft of this paper, 14 August 2021.
4
As of the writing of this paper, 14 August 2021.
5
Option chain quotes were retrieved from TD Ameritrade using the TOS platform.
6
These prices could be the actual market prices or the average between the bid and ask prices of the market.

References

  1. Bakshi, Gurdip, Charles Cao, and Zhiwu Chen. 1997. Empirical Performance of Alternative Option Pricing Models. The Journal of Finance LII: 2003–49. [Google Scholar] [CrossRef]
  2. Bin, Chen. 2007. Calibration of the Heston Model with Application in Derivative Pricing and Hedging. Master’s thesis, Department of Mathematics, Technical University of Delft, Delft, The Netherlands. Available online: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.377.2222&rep=rep1&type=pdf (accessed on 20 May 2022).
  3. Black, Fischer. 1989. How to use the holes in Black-Scholes. Journal of Applied Corporate Finance 1: 67–73. [Google Scholar] [CrossRef]
  4. Black, Fischer, and Myron Scholes. 1973. The pricing of options and corporate liabilities. The Journal of Political Economy 81: 637–54. [Google Scholar] [CrossRef] [Green Version]
  5. Boukai, Ben. 2020. How Much is your Strangle Worth? On the Relative Value of the delta-Symmetric Strangle under the Black-Scholes Model. Applied Economics and Finance 7. [Google Scholar] [CrossRef]
  6. Boukai, Ben. 2021. On the RND under Heston’s Stochastic Volatility Model. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3763494 (accessed on 20 May 2022).
  7. Cox, C. John, and Stephen A. Ross. 1976. The valuation of options for alternative stochastic processes. Journal of Financial Economics 3: 145–66. [Google Scholar] [CrossRef]
  8. Figlewski, Stephen. 2010. Estimating the Implied Risk Neutral Density for the U.S. Market Portfolio. In Volatility and Time Series Econometrics: Essays in Honor of Robert F. Engle. Edited by Tim Bollerslev, Jeffrey Russell and Mark Watson. Oxford: Oxford University Press. [Google Scholar]
  9. Figlewski, Stephen. 2018. Risk Neutral Densities: A Review. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3120028 (accessed on 20 May 2022).
  10. Gatheral, Jim. 2006. The Volatility Surface. Hoboken: John Wiley and Sons. [Google Scholar]
  11. Gilli, Manfred, Dietmar Maringer, and Enrico Schumann. 2019. Numerical Methods and Optimization in Finance, 2nd ed. Waltham: Elsevier/Academic Press, ISBN 978-0128150658. [Google Scholar]
  12. Grith, Maria, and Volker Krätschmer. 2012. Parametric Estimation of Risk Neutral Density Functions. In Handbook of Computational Finance. Edited by Duan Jin-Chuan, Härdle Wolfgang and Gentle James. Springer Handbooks of Computational Statistics. Berlin/Heidelberg: Springer. [Google Scholar]
  13. Heston, L. Steve. 1993. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Review of Financial Studies 6: 327–43. [Google Scholar] [CrossRef] [Green Version]
  14. Jackwerth, J. Carsten. 2004. Option-Implied Risk-Neutral Distributions and Risk Aversion. Charlotteville: Research Foundation of AIMR. [Google Scholar]
  15. Jackwerth, J. Carsten, and Mark Rubinstein. 1996. Recovering Probability Distributions from Option Prices. The Journal of Finance 51: 1611–31. [Google Scholar] [CrossRef]
  16. Lemaire, Vincent, Thibaut Montes, and Gilles Pagès. 2020. Stationary Heston Model: Calibration and Pricing of Exotics Using Product Recursive Quantization. Available online: https://arxiv.org/abs/2001.03101 (accessed on 20 May 2022).
  17. Merton, C. Robert. 1973. Theory of rational option pricing. The Bell Journal of Economics and Management Science 4: 141–83. [Google Scholar] [CrossRef] [Green Version]
  18. Mil’shtein, G. Noyhovich. 1975. Approximate Integration of Stochastic Differential Equations. Theory of Probability & Its Applications 19: 557–62. [Google Scholar]
  19. Mizrach, Bruce. 2010. Estimating Implied Probabilities from Option Prices and the Underlying. In Handbook of Quantitative Finance and Risk Management. Edited by Lee Cheng-Few, Lee Alice and Lee John. New York: Springer Science Business Media. [Google Scholar]
  20. Mrá¡zek, Milan, and Jan Pospíšil. 2017. Calibration and simulation of Heston model. Open Mathematics 15. [Google Scholar] [CrossRef]
  21. R Core Team. 2017. R: A Language and Environment for Statistical Computing. Vienna: R Core Team, Available online: https://www.R-project.org/ (accessed on 20 May 2022).
  22. Romo, Eudald, and Luis Ortiz-Gracia. 2021. SWIFT Calibration of the Heston Model. Mathematics 9: 529. [Google Scholar] [CrossRef]
  23. Schmelzle, Martin. 2010. Option Pricing Formulae Using Fourier Transform: Theory and Application. Technical Report. Available online: https://pfadintegral.com/articles/option-pricing-formulae-using-fourier-transform/ (accessed on 20 May 2022).
  24. Schumann, Enrico. 2011–2021. Numerical Methods and Optimization in Finance (NMOF), Manual. Package Version 2.5-0. Available online: http://enricoschumann.net/NMOF/ (accessed on 20 May 2022).
  25. Stein, M. Elias, and Stein C. Jeremy. 1991. Stock price distributions with stochastic volatility: An analytic approach. Review of Financial Studies 4: 727–52. [Google Scholar] [CrossRef]
  26. Savickas, Robert. 2005. Evidence on delta hedging and implied volatilities for the Black-Scholes, gamma, and Weibull option pricing models. The Journal of Financial Research 18: 299–317. [Google Scholar] [CrossRef]
  27. Stacy, Edney Webb. 1962. A generalization of the gamma distribution. The Annals of Mathematical Statistics 33: 1187–92. [Google Scholar] [CrossRef]
  28. Stacy, Edney Webb, and Mihram George Arthur. 1965. Parameter Estimation for a Generalized Gamma Distribution. Technometrics 7: 349–58. [Google Scholar] [CrossRef]
  29. Wiggins, B. James. 1987. Option values under stochastic volatility: Theory and empirical estimates. Journal of Financial Economics 19: 351–72. [Google Scholar] [CrossRef]
  30. Yu, Xisheng. 2020. Risk-Neutrality of RND and Option Pricing within an Entropy Framework. Entropy 22: 836. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The volatility ‘smiles’ of the the 15 October 2021 option series (calls) as observed and calculated on 13 August 2021 (EOD) for the three market index ETFs, SPY, IWM and QQQ and on 18 August 2021 (EOD) for the TLT ETF.
Figure 1. The volatility ‘smiles’ of the the 15 October 2021 option series (calls) as observed and calculated on 13 August 2021 (EOD) for the three market index ETFs, SPY, IWM and QQQ and on 18 August 2021 (EOD) for the TLT ETF.
Jrfm 15 00238 g001
Figure 2. The SPY case: the calculated HS, GG and BS implied RNDs along with the Monte Carlo distribution of the spot’s price S * and the corresponding values of the MSEs.
Figure 2. The SPY case: the calculated HS, GG and BS implied RNDs along with the Monte Carlo distribution of the spot’s price S * and the corresponding values of the MSEs.
Jrfm 15 00238 g002
Figure 3. The SPY case—the calculated delta functions under each of the pricing models, HS, GG and BS, along with the quoted delta per each strike K in the 15 October 2021 option series.
Figure 3. The SPY case—the calculated delta functions under each of the pricing models, HS, GG and BS, along with the quoted delta per each strike K in the 15 October 2021 option series.
Jrfm 15 00238 g003
Figure 4. The IWM case: (a) the HS, GG and BS implied RNDs along with the Monte Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Figure 4. The IWM case: (a) the HS, GG and BS implied RNDs along with the Monte Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Jrfm 15 00238 g004
Figure 5. The QQQ case: (a) the HS, GG and BS implied RNDs along with the Monte Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Figure 5. The QQQ case: (a) the HS, GG and BS implied RNDs along with the Monte Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Jrfm 15 00238 g005
Figure 6. The TLT case- (a) the HS, GG and BS implied RNDs along with the Monte-Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Figure 6. The TLT case- (a) the HS, GG and BS implied RNDs along with the Monte-Carlo distribution of the Spot’s price S * , and (b) the corresponding delta functions along with the quoted delta per each strike K in the option series.
Jrfm 15 00238 g006
Table 1. Calculated (excess) kurtosis and skewness measures for the three distributions depicted in Figure 1 for the SPY option data.
Table 1. Calculated (excess) kurtosis and skewness measures for the three distributions depicted in Figure 1 for the SPY option data.
MeasureHSGGBS
Kurtosis7.3026743.5364610.05234164
Skewness 2.050771 1.580122 0.1715114
Table 2. Model’s goodness-of-fit as measured the respective MSE for each of the four ETFs.
Table 2. Model’s goodness-of-fit as measured the respective MSE for each of the four ETFs.
ETFHSGGBS
SPY0.22264290.3394411.781981
IWM0.0019009680.014196280.3750478
QQQ0.024180130.065611341.193867
TLT0.034237480.046187250.04341321
Table 3. Comparison of the the quoted ATM delta Δ * of the four market ETF to those calculated under each of the three option pricing models.
Table 3. Comparison of the the quoted ATM delta Δ * of the four market ETF to those calculated under each of the three option pricing models.
ETFSATM K Δ * Δ BS Δ GG Δ HS
SPY445.924450.4970.5060.6380.663
IWM221.132210.5100.5160.5980.610
QQQ368.823690.5070.5030.6250.632
TLT149.351500.5110.4530.4770.467
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Boukai, B. The Generalized Gamma Distribution as a Useful RND under Heston’s Stochastic Volatility Model. J. Risk Financial Manag. 2022, 15, 238. https://doi.org/10.3390/jrfm15060238

AMA Style

Boukai B. The Generalized Gamma Distribution as a Useful RND under Heston’s Stochastic Volatility Model. Journal of Risk and Financial Management. 2022; 15(6):238. https://doi.org/10.3390/jrfm15060238

Chicago/Turabian Style

Boukai, Benzion. 2022. "The Generalized Gamma Distribution as a Useful RND under Heston’s Stochastic Volatility Model" Journal of Risk and Financial Management 15, no. 6: 238. https://doi.org/10.3390/jrfm15060238

APA Style

Boukai, B. (2022). The Generalized Gamma Distribution as a Useful RND under Heston’s Stochastic Volatility Model. Journal of Risk and Financial Management, 15(6), 238. https://doi.org/10.3390/jrfm15060238

Article Metrics

Back to TopTop