A Flexible Generalized Hyperbolic Option Pricing Model and Its Special Cases

A Flexible Generalized Hyperbolic Option Pricing Model and Its Special Cases Abstract We study a generalized hyperbolic (GH) time-changed Lévy process for option pricing and examine six three-parameter special cases: the variance gamma (VG) model of Madan, Carr, and Chang (1998), t, hyperbolic (H), normal inverse Gaussian (NIG), reciprocal hyperbolic (RH), and normal reciprocal inverse Gaussian (NRIG) option pricing models. We study the GH model’s moment properties of the associated risk-neutral distribution of logarithmic spot returns, and obtain an explicit pricing formula for European options facilitated by the time-change Lévy process construction. Using S&P 500 Index European options during low and high volatility sample periods, we compare the GH model empirically with existing benchmark models such as the finite-moment log-stable model and the Black–Scholes model. The GH model offers the best in- and out-of-sample performance overall, and a proposed t model special case generally outperforms the existing VG special case. We also present a stochastic volatility extension of the GH model. Empirical option prices indicate that the likelihood of extreme logarithmic stock returns is higher than that implied by the Black–Scholes (BS) model. Option prices also reveal that market participants pay more to protect themselves from losses than to pursue gains of equivalent magnitude. The implication is that the risk-neutral distribution of log-returns exhibits excess kurtosis and negative skewness (Madan and Milne, 1991; Eberlein and Keller, 1995). These two digressions from the normality assumption (Black and Scholes, 1973; Merton, 1973) are in part responsible for the poor empirical pricing performance of the BS model. To combat this deficiency, the generalized hyperbolic (GH) distribution may be used to improve option pricing as it accommodates skewness and thicker, semi-heavy tails (Barndorff-Nielsen, 1977). In this paper, we focus on addressing the challenges faced by GH option pricing in the setting in which spot returns are independent. Our contribution to this field is to propose a new form of the GH option pricing model, the flexible GH model, which contains four free parameters and can be conveniently estimated. In addition, we present six three-parameter option pricing models as special cases of the flexible GH model, including the variance gamma (VG), t, hyperbolic (H), normal inverse Gaussian (NIG), reciprocal hyperbolic (RH), and normal reciprocal inverse Gaussian (NRIG) models. With the exception of the VG model proposed by Madan, Carr, and Chang (1998), the remaining five option pricing models are innovations of this paper. To construct the flexible GH option pricing model, we generalize the BS model by the time change (or subordination) method (Clark, 1973). More specifically, the risk-neutral dynamics of the log-returns is modeled by a drifted Brownian motion subordinated by a stochastic time deformation process, which is defined to be a generalized inverse Gaussian (GIG) process (Barndorff-Nielsen and Halgreen, 1977). This construction generates a class of pure-jump, infinite activity processes (Barndorff-Nielsen and Shephard, 2001).1 A typical sample path of the resulting GH process jumps infinitely many times within an infinitesimally small interval, thus enabling the model to capture both discrete and continuous asset price movements (Daal and Madan, 2005). The GH process constructed this way belongs to the time-changed Lévy process family (Barndorff-Nielsen and Shephard, 2001; Carr and Wu, 2004; Huang and Wu, 2004). Taking this time change approach is especially effective for developing more plausible and tractable models for asset pricing.2 In terms of derivative pricing, the methodology lends itself to the derivation of an explicit pricing formula for European options, as the characteristic function of the logarithmic spot price can be conveniently obtained. Furthermore, with a judicious choice of the Lévy process and the time deformation process, the modeler can easily incorporate stochastic volatility (SV) to relax the i.i.d. restriction on the spot return dynamics under the baseline models. Tractability in deriving the option pricing formula is preserved under the time-changed Lévy process framework, provided that the characteristic functions of the relevant processes are known. For example, with the use of a time-dependent time deformation process, the baseline GH option pricing model can be readily extended to allow for mean-reverting SV of spot returns. This extension is detailed in Section 4. As another example, short- and long-range dependence in volatility can be induced by choosing the time deformation process to be a superposition of many time-dependent processes (Finlay and Seneta, 2012), although its empirical performance is yet to be examined in the literature. Generalizing this idea further, it is possible to develop option pricing models with multi-factor volatility dynamics, such as the models in Andersen, Fusari, and Todorov (2015) and Gruber, Tebaldi, and Trojani (2015). The multi-factor structure provides the necessary flexibility in modeling the underlying risk-neutral spot dynamics. Empirical evidence suggests that they are capable of capturing option stylized facts (e.g., implied volatility surface) more accurately than the conventional models. Eberlein and Prause (2002) and Prause (1999) consider a GH model for option pricing. There are several important differences that separate our approach from theirs. First, using time-series data of stock returns, Eberlein and Prause (2002) and Prause (1999) estimate the GH model, which is then used to price the options and obtain the model-implied volatilities. This is made possible by means of Esscher transform, which allows them to recover the risk-neutral spot price dynamics from the statistical measure. We take a more direct route by specifying the spot price dynamics under the risk-neutral measure. The risk-neutral parameters of our flexible GH model are then estimated directly using option data. Our modeling approach is therefore in line with the risk-neutral option pricing literature (e.g., Heston, 1993; Bates, 1996; Madan, Carr, and Chang, 1998; Pan, 2002; Carr and Wu, 2003). Second, Eberlein and Prause (2002) and Prause (1999) express the European option price in terms of the integral of a τ-fold convolution of the GH densities (τ is the time-to-maturity) which cannot be further simplified as the GH distribution is not closed under convolution. To circumvent the numerical challenge, they employ Fourier inversion to obtain the cumulative probabilities from the characteristic function of the log-return process. By contrast, using the Carr and Madan (1999) approach, we obtain an explicit option pricing formula in terms of the characteristic function of the logarithmic spot price process. Third, in terms of numerical optimization, Eberlein and Prause (2002) and Prause (1999) estimate the risk-neural parameters by first keeping the index parameter fixed and then optimizing over the other parameters.3 We are able to overcome the numerical difficulty by estimating all four parameters of the flexible GH model jointly and efficiently. This is the combined result of a prudent choice of model parameterization, appropriate variable transformation, and a larger sample offered by the option panel compared with the time series of spot prices. In addition, we conduct an empirical study of the flexible GH model and its special cases by investigating its performance in modeling S&P 500 Index options. For a fair comparison, the GH model is evaluated against benchmark option pricing models that assume i.i.d. spot returns, including the finite moment log-stable (FMLS) and BS models.4 Among the special cases of GH, we will evaluate the NIG model, and additionally closely examine the VG and t models, which Eberlein and Prause (2002) and Prause (1999) precluded from being estimated under the risk-neutral setting. A direct comparison between the VG and t models is theoretically and practically motivated, since the models are complementary limiting cases of the flexible GH model (see Subsection 2.1), yet the VG model is much more widely used than the t model in empirical option pricing. The main empirical findings of our paper are summarized as follows. First, among all six models under consideration, the flexible GH model yields the best in-sample fitting and out-of-sample forecast across the entire panel of option prices; the added flexibility offered by the extra shape parameter is indispensable for the superior performance. Second, the t model often performs as good as, and sometimes even better than, the flexible GH model in fitting and predicting the prices of certain option categories. For example, the t model often provides an adequate fit to the left tail of the log-return distribution, as revealed by its superior performance in modeling the in-the-money (ITM) put options in the main estimation. Third, as shown in the robustness analysis, the superior overall out-of-sample performance of the flexible GH model is preserved over a turbulent sample period in 2011, although the category-specific results are mixed—the FMLS model surpasses the GH model in predicting the prices of at-the-money (ATM) and ITM options, thus pointing to the necessity of excess kurtosis in the spot price distribution for accurate option pricing. In summary, our contributions are three-fold. First, we study a flexible GH option pricing model in the time-changed Lévy process framework. Taking this dynamic perspective is fruitful in option pricing as it motivates the way the GH model is parameterized. The flexible GH model encompasses six option pricing models as special cases, including the popular VG model, as well as new models which employ certain familiar distributions, such as the t distribution. Second, we obtain an explicit, semi-closed-form option pricing formula under our GH model. The formula is derived using the characteristic function approach (Carr and Madan, 1999). Third, we investigate the empirical performance of the GH model and some of its special cases. By comparison with benchmark models such as the FMLS model (Carr and Wu, 2003) and the BS model, we illustrate the convenience of calibrating the GH model, and its consistent superiority in delivering in-sample fitting and out-of-sample prediction over the panel of S&P 500 Index options. To our knowledge, this paper provides the first systematic empirical analysis of the GH model and some of its special cases, including the VG, t, and NIG models. Although our empirical study focuses on Lévy-based option pricing models that assume i.i.d. spot returns, the empirical results would shed light on future work on evaluating the empirical performance of SV extensions of such models, some of which are developed in Section 4. It would be further elucidating, for example, to compare the GH–SV model with the conventional Heston (1993) model, whose underlying spot dynamics are characterized by a time-changed GH process and time-changed Brownian motion, respectively. The paper is structured as follows. Section 1 begins with a description of a subordinated process and presents a corresponding option pricing framework. In Section 2, we formulate the flexible GH option pricing model and its six special cases, including the two limiting cases. Section 3 comprises our empirical study. The data description, parameter estimates, in-sample fit, orthogonality test, and out-of-sample pricing results are provided therein. In Section 4, we discuss the extension of the GH model to incorporate SV. Section 5 concludes the paper. 1 Time-Changed Lévy Processes and Option Pricing We focus on a generic option pricing framework in which the log-return process of the underlying asset is driven by a time-changed Lévy process [Barndorff-Nielsen and Shephard, 2001; Carr et al., 2003 (hereafter CGMY, 2003); Carr and Wu, 2004]. Let Lt≡L(t) denote a Lévy process, and let Tt≡T(t) be an increasing, positively valued time deformation process. The resulting time-changed Lévy process is given by Zt≡Z(t)=L(T(t)). (1) The time deformation process Tt represents the total elapsed time in a transformed time scale which may reflect the varying rate at which economic and business activities occur. As will be explored in this paper, this modeling framework offers flexibility in the specification of the Lévy process and the time deformation process. It encompasses many option pricing models in the literature. In particular, the time deformation process Tt can be deterministic (e.g., an identity process Tt=t, as in BS model), a pure-jump Lévy process with stationary and independent increments (e.g., a gamma process, as in Madan, Carr, and Chang (1998), or a time-dependent process, possibly driven by the same random process that drives Lt (e.g., CGMY, 2003; Carr and Wu, 2004; Huang and Wu, 2004). To preserve the business clock interpretation, we require that Tt satisfies the following conditions: T0=0. Tt is non-negative and increasing with t. E[T1]=1. The third condition is a normalizing restriction on Tt. If Tt is a Lévy process, then, by infinite divisibility, this condition implies that E[Tt]=t for all positive time t. To derive the option pricing formula, it is instrumental that we obtain the characteristic function of the time-changed process Zt. Defining the characteristic functions of Lt and Tt by φLt(u)=E(eiuLt) and φTt(u)=E(eiuTt), respectively, we can compute the characteristic function of Zt as follows φZt(u)=φTt(−iφL1(u)). (2) We assume that Zt drives the log-return process of the underlying asset. More specifically, under the risk-neutral probability space5 (Ω,Q,{Ft}) equipped with the filtration {Ft}, the spot price of the underlying asset evolves according to the following dynamics: for all t  ≥  0, St=S0 exp{(r−q+ω)t+Zt}. (3) For simplicity, the risk-free rate r and the dividend yield q are assumed to be constant in our analysis. The drift adjustment ω is required to ensure that the discounted spot price process is a martingale under Q. It is given by ω=−1t log φZt(−i)=−log φZ1(−i) (see Appendix A). As a result, the characteristic function of the log-price is given by φ log St(u)=S0iu exp{iu(r−q+ω)t}φZt(u). (4) Using the characteristic function approach to option pricing (e.g., Scott, 1997; Bakshi and Madan, 2000; Carr and Wu, 2004), the price of an European option written on this underlying asset with the spot price process, defined in Equation (3), can be obtained. The time-t price of an European call with strike price K and time-to-maturity τ = T– t is Ct=EQ[e−rτ⁡max(ST−K,0)|Ft]=Ste−qτΠ1−Ke−rτΠ2, (5) where EQ[·] denotes the expectation taken under Q, and the quantities Π1 and Π2 are the prices of some Arrow–Debreu securities (Bakshi and Madan, 2000), given by Π1=12+1π∫0∞Re(e−iω log(K)φ log St(ω−i)iωφ log St(−i))dω,Π2=12+1π∫0∞Re(e−iω log(K)φ log St(ω)iω)dω. By the put-call parity, the price of the corresponding European put option with the same strike price and maturity date is Pt=Ke−rτ(1−Π2)−Ste−qτ(1−Π1). (6) 2 The Baseline Model In this section, we introduce the baseline GH option pricing model. Recall from Equation (3) that the spot price is driven by the time-changed Lévy process Zt. Under the flexible GH model specification, the time-changed Lévy process Zt is constructed by subordinating a Brownian motion with drift, B(t), to a GIG process g(t). More precisely, the time-changed Lévy process is defined by Zt:=Z(t)=B(g(t))=θg(t)+σW(g(t)), (7) where θ is the drift parameter, σ  >  0 is the volatility parameter, and W(t) is the standard Wiener process. The GIG process g(t) is a positively valued Lévy process such that g1:=g(1) follows the GIG distribution with parameters p, γ, and δ. The probability density function of the GIG distribution is provided in Appendix B.1. In particular, by defining the parameter ζ:=γδ, the mean is E[g1]=ζKp+1(ζ)γ2Kp(ζ), where Kj(·) is the modified Bessel function of the third kind with index j (Jørgensen, 1982). The mean of g1 is well defined provided that the parameters p and ζ satisfy the conditions: (i) ζ>0 if −1≤p≤0,(ii) ζ≥0  otherwise. By imposing the normalization E[g1]=1, we obtain the parameter restriction γ2=ζKp+1(ζ)Kp(ζ), (8) and hence the GIG distribution is characterized by two parameters: the index parameter p and the shape parameter ζ. From a static point of view, the process at time 1, Z1, follows the GH distribution, which can be viewed as a normal mean–variance mixture distribution with a GIG mixing density. For any fixed t, the distribution of Zt is skewed and leptokurtic, due to the mixing in the mean through the Brownian motion’s drift θ. From a dynamic point of view, the time-changed Lévy process Zt is constructed by subordinating a drifted Brownian motion to a GIG process. The distribution of g(t) is obtained as the t-fold convolution of the GIG distribution of g1. In particular, for any t  ≥  0, the increment g(t  +  1)–g(t) follows the same distribution as that of g1. Note that g(t) is not distributed as GIG for t  ≠  1, as the GIG distribution is not preserved under convolution. It follows from this dynamic construction that the marginal distribution of Zt is not GH in general except at t  =  1. Nevertheless, by infinite divisibility of the GIG distribution (Barndorff-Nielsen and Halgreen, 1977), the GIG process g(t) is infinitely divisible. As a result, the time-changed process Zt is a Lévy process, and so it is infinitely divisible and has stationary and independent increments. For a Lévy process, the characteristic functions of Zt and Z1 are related by φZt(u)=[φZ1(u)]t (e.g., Barndorff-Nielsen, Mikosch, and Resnick, 2001). The central moments of Zt grow over time just like any generic Lévy process. The mean, variance, and the third central moment of Zt increase linearly with t. In terms of standardized central moments, the skewnesses of Zt and Z1 are related by skew(Zt)=(1/t)skew(Z1) for t  >  0, and the kurtosis of Zt is a weighted average of the kurtosis of Z1 and that of a normal distribution, that is, kur(Zt)=(1/t)kur(Z1)+3(1−1/t) for t  >  0. Given that Z1 follows a GH distribution, we have skew(Z1)≠0 if the drift parameter θ is non-zero, and kur(Z1)  >  3. The same moment properties are therefore valid for Zt for all t > 0. A study of the moment properties is detailed in Appendix B. To price a European option under the flexible GH model, it is necessary to get the characteristic function of Zt. Using the characteristic function of the GH distribution (see Appendix B), we obtain φZt(u)=[1−2γ2(iuθ−12σ2u2)]−pt2[Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ)]t. Provided that θ<γ2−σ22, the drift adjustment is well-defined and is given by ω=−log φZ1(−i)=−p2log[1−2γ2(θ+σ22)]+log[Kp(ζ1−2γ2(θ+σ22))Kp(ζ)]. The European option under the flexible GH model can then be priced using Equations (5) and (6). 2.1 Special Cases of the Flexible GH Model In the family of risk-neutral spot price models with independent and stationary increments, the baseline GH option pricing model is quite flexible in that it encompasses a number of particular cases. In this section, we study six of its special cases, all of which are indexed by three parameters. These six models can be classified into two groups: four special cases obtained by restricting the value of the index parameter p, while the other two obtained as limiting cases as the shape parameter ζ:=δγ→0. Table 1 summarizes the parameter restrictions required to obtain the six special cases of the flexible GH model. Table 1 Special cases of the GH option pricing model GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) Notes: ζ=γδ. For the H, RH, NIG and NRIG models, γ is given by (8). For all models, E [g1]=1. Table 1 Special cases of the GH option pricing model GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) Notes: ζ=γδ. For the H, RH, NIG and NRIG models, γ is given by (8). For all models, E [g1]=1. The first group of three-parameter sub-models is obtained as follows. By setting p  =  1, the flexible GH model reduces to a H model, in which the mixing random variable g1 follows the positive H mixture distribution. Setting p  =  −1 leads to a RH model, in which the mixing distribution is reciprocal positive H. For p=−(1/2), we obtain the NIG distribution, generated by a normal mixture of the inverse Gaussian distribution. For p=1/2, we obtain the NRIG distribution, generated by a normal mixture of the reciprocal inverse Gaussian distribution. The second group of sub-models is obtained by restricting the shape parameter ζ to 0 and the sign of the index parameter p. By letting δ→0 and restricting p  >  0, we obtain the VG model (Madan, Carr, and Chang, 1998), in which g1 follows the gamma distribution with parameters p and γ2/2. By letting γ → 0 and restricting p < 0, we obtain the t model, in which g1 follows the reciprocal (or inverse) gamma distribution with parameters-p and δ2/2. With the normalization condition E[g1]=1 in place, we deduce that γ=2p in the VG case and δ=2(−p−1) in the t case. Figure 1 demonstrates the flexibility of the GH model in the modeling of option prices relative to its six special cases. The plotted surface represents the prices of an ITM call option under the GH model for varying p and ζ (the drift and volatility parameters are fixed). The option prices associated with the first group of sub-models (H, NRIG, NIG, and RH models) are obtained on the four dotted cross-sections of the surface at the respective fixed values of p, while those associated with the second group of sub-models (VG and t models) are found, respectively, along the positive and negative sections of the p axis at the boundary ζ = 0. When both the index and shape parameters are close to zero, the GH model implies a high kurtosis in the risk-neutral distribution of spot prices. This translates into more extreme option prices (lower prices under NIG, RH, and t with p close to 0, and higher prices under NRIG, H, and VG with p close to 0). The GH model offers the most versatility as manifested by the wide range of option prices it generates over all possible combinations of the index and shape parameters. Figure 1 View largeDownload slide Simulated option prices under the GH model. Notes: The price of an ITM call option with a strike price equal to 95% of the spot price, time-to-maturity equal to 3 weeks, θ=−0.04 and σ=0.17. The six special cases are RH at p = –1, NIG at p=−12, NRIG at p=12, hyperbolic (H) at p = 1, variance gamma at ζ = 0 and p > 0, and the t model at ζ = 0 and p<−1. At ζ = 0 and −1≤p≤0, the mean of the GH distribution is undefined. Figure 1 View largeDownload slide Simulated option prices under the GH model. Notes: The price of an ITM call option with a strike price equal to 95% of the spot price, time-to-maturity equal to 3 weeks, θ=−0.04 and σ=0.17. The six special cases are RH at p = –1, NIG at p=−12, NRIG at p=12, hyperbolic (H) at p = 1, variance gamma at ζ = 0 and p > 0, and the t model at ζ = 0 and p<−1. At ζ = 0 and −1≤p≤0, the mean of the GH distribution is undefined. 2.2 The VG Model We note that the characteristic function of the logarithmic spot price under the VG model can be simplified to φ log St(u)=S0iu exp{iu(r−q+ω)t}[1−ν(iuθ−12σ2u2)]−tν, (9) where ω=1νlog[1−ν(θ+12σ2)] for θ<(1ν−σ22). This formula can be obtained as the limit of the corresponding characteristic function under the GH model (see Appendix C). This is equivalent to the VG option pricing model of Madan, Carr, and Chang (1998), which contains three parameters: the drift θ, the volatility σ, and the mixing distribution variance ν (which is equal to 1/p using our index parameter). The distributional properties of the VG model are examined in Appendix B.2. 2.3 The t Model Let us turn to the t option pricing model.6 Under this model, the log-return of the underlying asset follows a skewed t distribution with –2p degrees of freedom. The characteristic function of the logarithmic spot price under the VG model is given by φ log St(u)=S0iu exp{iu(r−q+ω)t}×[2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν(−4(1ν−1)(iuθ−12σ2u2))]t, (10) where ω=−log[2(1ν−1)12νΓ(1ν)[−(θ+12σ2)]12νK1ν(−4(1ν−1)(θ+12σ2))]for θ<−σ22, and Kj(·) is the modified Bessel function of the third kind with index j. This formula can be obtained as the limit of the corresponding characteristic function under the flexible GH model (see Appendix C). The spot price dynamics under the t model is equivalent to the GH-skewed Student’s t model in Aas and Haff (2006) after re-parameterization, except that the latter specifies the spot price process under the statistical measure and is not used for option pricing. The distributional properties of the t model are examined in Appendix B.3. 3 Empirical Study 3.1 The Data The data used for our main estimation are the daily prices of S&P 500 Index European options. Put options are featured in our sample because, compared with their call counterparts, their prices carry more information about the left tail of the price distribution, corresponding to the side to which the price distribution is found to be skewed empirically (Madan and Milne, 1991). The sample period is chosen to be 2012 with 250 trading days. It contains some moderate upward and downward trends in the S&P 500 Index. The volatility index (VIX) fluctuates mildly between 13% and 27% over the period. Both series are displayed in the bottom two panels in Figure 2. In the robustness analysis in Section 3.8, we will choose the most turbulent period in 2011 as our alternative sample and repeat the analysis. Figure 2 View largeDownload slide Time-series plots of daily risk-neutral volatility (σ) and skewness (θ) parameter estimates under the GH option pricing model, and the plots of the S&P 500 Index and VIX. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The plots include the volatility parameter (top panel), the skewness parameter (second panel), the S&P 500 Index (third panel), and VIX (bottom panel), the VIX that measures the implied volatility of S&P 500 Index options. Figure 2 View largeDownload slide Time-series plots of daily risk-neutral volatility (σ) and skewness (θ) parameter estimates under the GH option pricing model, and the plots of the S&P 500 Index and VIX. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The plots include the volatility parameter (top panel), the skewness parameter (second panel), the S&P 500 Index (third panel), and VIX (bottom panel), the VIX that measures the implied volatility of S&P 500 Index options. The cross-sectional data structure for the put option sample used in the main estimation is summarized in Table 2. To retain a majority of put options that are more liquidly traded in the market, we focus on the puts whose moneyness (the strike-to-spot price ratio) lies between 0.94 and 1.06, and whose time-to-maturity is longer than a week. On each day, the sample retains only the three batches of put options belonging to the short-term, medium-term, and long-term categories, corresponding to the nearest-, second-nearest-, and third-nearest-to-maturity put options, respectively. This leaves us with 60 put options on average per day. In total, the option panel consists of 15,058 put prices. To price the European puts, we use the annualized yield of the 1-month U.S. Treasury bill as the risk-free rate. The time series of S&P 500 Index is adjusted by the mean annualized dividend yield of the stocks. Table 2 Options data characteristics of the main estimation sample Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Notes: A contingency table for the sample of S&P 500 Index European put option prices in 2012. The sample is selected according to the criteria in the Data section. n = 15,058. Average option prices are displayed in square brackets. K is the strike price, S is the spot price, and τ is the time-to-maturity. Table 2 Options data characteristics of the main estimation sample Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Notes: A contingency table for the sample of S&P 500 Index European put option prices in 2012. The sample is selected according to the criteria in the Data section. n = 15,058. Average option prices are displayed in square brackets. K is the strike price, S is the spot price, and τ is the time-to-maturity. Due to the way the sample is selected, the time-to-maturities exhibit some artificial cyclicity over time. This can be seen from the time-series plot of the average time-to-maturity at the bottom panel of Figure 3. Note that the S&P 500 Index options expire on the third Friday of each month. As the maturity date approaches, the time-to-maturities of the options in the sample decreases gradually over time until some of them (i.e., the time-to-maturities of those puts in the short-term category) falls below 1 week, at which point these nearest-to-maturity options will be excluded from the sample and replaced by a new batch of options that forms the long-term category. The original long-term (medium-term) category now becomes the medium-term (short-term) category. This roll-over of the sample, which occurs around the second Friday of each month, leads to an abrupt increase in the time-to-maturities by approximately 1 month across the three term categories. It is important to bear this in mind when we interpret the time-series variation of the parameter estimates. Figure 3 View largeDownload slide Time series plots of daily risk-neutral shape parameter estimates from different models, and the plot of daily averages of time-to-maturity in the main estimation sample in 2012. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The shape parameters include the tail index parameter α under the FMLS model (top panel), and the index parameter p under the GH model and its special cases (second panel), and the shape parameter ζ under the GH and NIG models (third panel), with ζ→0 under the VG and t models. The average time-to-maturity (τ, in days) is computed by averaging the time-to-maturity of the cross-section of puts each day (bottom panel). Figure 3 View largeDownload slide Time series plots of daily risk-neutral shape parameter estimates from different models, and the plot of daily averages of time-to-maturity in the main estimation sample in 2012. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The shape parameters include the tail index parameter α under the FMLS model (top panel), and the index parameter p under the GH model and its special cases (second panel), and the shape parameter ζ under the GH and NIG models (third panel), with ζ→0 under the VG and t models. The average time-to-maturity (τ, in days) is computed by averaging the time-to-maturity of the cross-section of puts each day (bottom panel). 3.2 Model Calibration and Parameter Estimation The option pricing models are calibrated daily by minimizing the sum of squared percentage pricing errors SSPEt(Θ) on each date t: SSPEt(Θ)=∑i=1nt(O^it(Θ)−OitOit)2, where Θ is the vector of model parameters, nt is the cross-sectional sample size on date t, and O^it(Θ) and Oit are the predicted and observed option prices (in dollars), respectively. The optimization search for Θ is carried out over the parameter space using the Nelder–Mead simplex algorithm (Gilli and Schumann, 2012). It is implemented through the MATLAB function fminsearch. To mitigate the sensitivity of the solution on the initial value of the parameter vector, a preliminary search is conducted as follows: starting with many different candidate vectors for Θ[0], we run a first-stage optimization with a fixed number of iterations (e.g., 10), then we select the optimal candidate vector that yields the smallest SSPEt(Θ). Starting at this optimal candidate vector, the numerical algorithm is then run until the local minimum is found. 3.3 In-Sample Estimation In this section, we present the main estimation results and examine the in-sample fit of various models. Table 3 summarizes the risk-neutral parameter estimates under the BS model, the FMLS model, as well as the GH model and three of its special cases (NIG, VG, and t). Table 3 Summary statistics of risk-neutral parameter estimates (in-sample) Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Notes: Parameter estimates and root mean-squared percentage error (RMSPE) are obtained daily for the year 2012, which includes 15,058 put prices over 250 days. The average daily sample size is 60. Table 3 Summary statistics of risk-neutral parameter estimates (in-sample) Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Notes: Parameter estimates and root mean-squared percentage error (RMSPE) are obtained daily for the year 2012, which includes 15,058 put prices over 250 days. The average daily sample size is 60. Starting with the volatility parameter σ, which appears in all models being considered, the FMLS model has the smallest average equal to 9.9%, followed by the BS model (14.6%), the VG model (15.0%), the NIG model (15.5%), and then the t and GH models (both 15.8%). The t model yields the single largest estimate equal to 24.8%. The skewness parameter θ is contained in all models except for BS and FMLS.7 It is estimated to be negative on average with a similar value across models: −0.044 for the NIG and VG models, −0.047 for the t model, and −0.046 for the GH model. However, the most extreme estimate is obtained under the t model: −0.149, compared with −0.108 for the NIG model, −0.100 for the VG model, and −0.098 for the GH model. The time-series plots of the daily estimates of σ and θ under the GH model are displayed in the top two panels of Figure 2. They reveal a negative association between the volatility and skewness parameters of the log-return distribution, confirming the presence of leverage effect (Christie, 1982): the volatility tends to be higher when the log-returns are more negatively skewed. Let us turn to the shape parameters (α for FMLS, and p and ζ for GH and its special cases). Under the FMLS model, the estimates of the tail index α lie in a narrow range between 1.848 and 1.999, with an average estimate equal to 1.927. By the properties of the asymmetric stable distribution (Property 1.5 of Carr and Wu, 2003), the log-return distribution has infinite moments of order 2 and higher if α  <  2. Under the GH model and its special cases, the shape is controlled by p (the index parameter) and ζ. Recall that p is constrained to be positive under the VG model, negative under the t model, and equal to −0.5 under the NIG model, and that ζ → 0 under the VG and t models, but ζ is free under the NIG model. The smaller the magnitude of p and ζ, the fatter the tail of the log-return distribution. From the empirical estimation, the VG model yields an average estimate of p equal to 9.305, which is higher in magnitude than that under the t model (−5.550, which translates into 11.1 degrees of freedom). The p estimates under the GH model average to −5.151 and range over −30.180 to 8.014, with a minimum magnitude of 1.117. These estimates are more in line with the t model than with the VG model. The estimation of p under the flexible GH model indicates that the special cases of GH, including RH (p  =  –1), NIG ( p=−(1/2)), NRIG ( p=1/2), and H (p = 1), are not supported by the data. In terms of the shape parameter ζ, the average estimate is 4.910 under the NIG model, but the average is only 0.417 under the flexible GH model. The discrepancy indicates that, due to the restriction on p, ζ needs to be inflated substantially under the NIG model in order to capture the excess kurtosis of the risk-neutral distribution. The time series of daily estimates of p under the various GH sub-models is plotted in Figure 3. It reveals a cyclical behavior in the degree of tail-thickness in the log-return distribution that aligns with the option expiration schedule.8 The tail of log-return distribution is thick at the start of the trading month (low magnitude of p), and gradually becomes thinner (increasing magnitude of p) toward the second Friday of each month as there is diminishing information uncertainty. Then, when the option data set rolls over to exclude options with less than 1 week to expiry and include a new batch of options, the index parameter p falls sharply in magnitude and the log-return distribution becomes fat-tailed again. This pattern repeats itself every month. 3.4 In-Sample Pricing Performance We turn now to the first of the performance measures, in-sample pricing error, also reported in Table 3. As expected, the daily-averaged root mean-squared percentage error (RMSPE) for the BS model is the highest at 16.69%, since it is the most restrictive model with only one parameter controlling the dispersion of the log-return distribution. The FMLS model reduces the RMSPE to 12.00%, the reduction in error representing an improvement to in-sample fit offered by the kurtosis parameter and a fixed negative skewness parameter. The VG, NIG, and t models further reduce the RMSPE to 9.50%, 9.10%, and 8.38%, respectively. The further improved fit is due to the flexibility in the skewness parameter and the evidently more correct specification of the log-return distribution over the FMLS model. The superior performance of the t model over the VG model indicates that the t distribution is a more plausible model empirically for the risk-neutral log-return distribution.9 The superiority of the t model over other sub-models is corroborated by the GH model’s parameter estimates: the index parameter p estimates are predominantly negative, and the estimates of the shape parameter ζ are close to zero. The in-sample fit of the t model, according to the average RMSPE, is as good as that of the GH model (8.38%). From the main estimation results, the GH model family clearly outperforms the BS and FMLS models. Allowing p to be free is crucial to the observed improvement in the in-sample fit, although the benefit of having an additional shape parameter ζ is not fully exploited in this particular dataset. We however wish to emphasize that, while the results are valid for the baseline dataset with sample period spanning 2012, the conclusions may depend on the sample period used. For comparison, the results of the robustness analysis based on a different sample period are presented in Section 3.8. 3.5 In-Sample Performance by Option Categories While Table 3 offers insight to overall in-sample fit, in Table 4 we examine how the models perform in terms of in-sample fit for different types of put options in the cross section, classified according to different moneyness and time-to-maturity groups as defined in Table 2. Table 4 In-sample pricing RMSPEs and mean absolute percentage errors (MAPEs) In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 Notes: The smallest error measure within each group is italicized and underlined. The bold values show the overall error measure. Parameters are estimated using all options on a given day regardless of their time-to-maturity and moneyness (an average of 60 observations per day, and 250 days in year 2012). Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 250 testing days collectively. The sample size is n = 15,058. Table 4 In-sample pricing RMSPEs and mean absolute percentage errors (MAPEs) In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 Notes: The smallest error measure within each group is italicized and underlined. The bold values show the overall error measure. Parameters are estimated using all options on a given day regardless of their time-to-maturity and moneyness (an average of 60 observations per day, and 250 days in year 2012). Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 250 testing days collectively. The sample size is n = 15,058. Let us start by comparing the performance across the moneyness categories. For ATM options, the GH model achieves the best in-sample pricing performance compared with other models: its average RMSPE is 6.78%, compared with 17.15% for the BS model, 11.02% for the FMLS model, 8.07% for the VG model, 7.97% for the NIG model, and 6.79% for the t model. For ITM options, the t model offers the best in-sample fit, with an average RMSPE of 8.24%, outperforming the GH model (8.27%). As all options in our sample are European puts, this suggests that the t distribution provides an optimal fit of the left tail of the log-return distribution. As for out-of-the-money (OTM) options, the GH model ties with the t model, and both of the models clearly outperform the other competitors. Similarly, a comparison across the different time-to-maturity categories reveals that the GH and t models give the best in-sample fit among all models, with the t model comes slightly ahead in pricing the short-dated puts, and the GH model becomes slightly better in pricing the medium-dated puts. Given that the FMLS model is outside the GH model family, one may be curious whether its peculiar features (e.g., infinite kurtosis of the log-return distribution) would facilitate more accurate pricing of certain options. Table 4 reveals that the FMLS model offers the best fit to the long-dated, ITM puts, achieving the smallest average RMSPE of 9.31% among all other models. As the price of a long-dated, ITM put is sensitive to the extreme left tail of the distribution of log-returns from now until maturity, the empirical result points to the necessity of excess kurtosis to model the extreme negative returns over a long time horizon. In this sense, the negatively skewed α-stable distribution under the FMLS model provides a better fit than the GH model, in which the kurtosis of log-returns tends to that of a normal distribution as the time horizon increases (see Appendix B.1). Nevertheless, the FMLS model is not adequate in modeling other sections of the long-range log-return distribution. Its RMSPE for long-dated OTM puts (14.31%) is worse than all other models, including the BS model (12.73%). Even though the log-return distribution under the FMLS model has the fattest tails, the empirical fit in the parts other than the left tail is not the best compared with other models with finite kurtosis. For robustness analysis of the in-sample performance measure, the results for the mean absolute percentage pricing error (MAPE) are also computed. While the average MAPE measures are generally smaller than the average RMSPE measures, probably due to the effect of extreme pricing errors on RMSPE, changing the measure to MAPE retains the optimal model in any given option category. 3.6 Orthogonality Test A comprehensive cross-sectional analysis of the pricing errors leads us to the second yardstick for assessing model adequacy. A well-specified model should not only minimize some loss function on pricing errors (e.g., RMSPE and MAPE) but also yield individual pricing errors that are independent of moneyness, time-to-maturity, and interest rates (Rubinstein, 1985). Furthermore, the pricing error measure used for assessing model adequacy needs to accord with that used for parameter estimation (Christoffersen and Jacobs, 2004). To ensure consistency with the model estimation procedure in Section 3.2, we therefore evaluate the model fitting using the sequence of option pricing errors relative to the observed option prices: εit=Oit(Θ^t)−OitOit, where Oit is the observed option price in dollars, and Oit(Θ^t) is the option price estimated by the model. We investigate the orthogonality of the relative pricing errors by regressing εit on some option and market attributes, including moneyness (defined as the ratio of the strike to the spot price), squared moneyness, time-to-maturity, squared time-to-maturity, and the risk-free interest rate. More specifically, we consider the following regression model for the orthogonality test: εit=b0+b1(KiSt)+b2(KiSt)2+b3τit+b4τit2+b5rt+eit, where Ki and τit are the strike price and time-to-maturity of the ith option at time t, St is the S&P 500 Index level at time t, rt is the risk-free interest rate at time t, and eit is the random error term. If the model is well specified, we expect close-to-zero regression coefficients and a low regression R2, both signifying greater orthogonality of the relative pricing errors. We run the regression on the entire 250-day sample consisting of 15,058 observations of relative pricing error. The results are presented in Table 5. As revealed by the t tests on the regression coefficients and the F tests on the regression R2, all models are generally misspecified to various degrees. Comparing the extent of error orthogonality across models, the BS model yields the highest regression R2 (72.0%), followed by the FMLS model (46.5%), the GH model (34.6%), the t model (34.5%), the VG model (32.8%), and the NIG model (28.7%). While NIG model achieves the least predictable relative pricing errors, its total sum of squared (TSS) pricing errors is larger in magnitude compared with the t and GH models, which have the smallest TSS among all models being considered. There is a trade-off to be made—from a practitioner’s perspective the overall pricing accuracy is likely to take precedence over the orthogonality of pricing errors when choosing among the models. Table 5 Results of orthogonality test Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Notes: S&P 500 Index put options from 2012 are used, with sample size n = 15,058. Moneyness is the strike-to-spot price ratio, and time-to-maturity is in years. Heteroskedastic standard errors (White, 1980) are shown in parentheses, * indicates statistical significance at the 1% level of significance. The critical t-statistics are, respectively, t(0.005,15052)=±2.58 and t(0.025,15052)=±1.96. TSS is the total sum of squared errors. At the 1% level of significance, the critical F(0.01,5,15052)-statistic is 3.02. Table 5 Results of orthogonality test Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Notes: S&P 500 Index put options from 2012 are used, with sample size n = 15,058. Moneyness is the strike-to-spot price ratio, and time-to-maturity is in years. Heteroskedastic standard errors (White, 1980) are shown in parentheses, * indicates statistical significance at the 1% level of significance. The critical t-statistics are, respectively, t(0.005,15052)=±2.58 and t(0.025,15052)=±1.96. TSS is the total sum of squared errors. At the 1% level of significance, the critical F(0.01,5,15052)-statistic is 3.02. 3.7 Out-of-Sample Pricing Performance Finally, we examine the out-of-sample pricing performance of the option pricing models. This is motivated from both a practical (price prediction) as well as a statistical consideration. While in-sample pricing performance benefits from having additional parameters, this is usually not the case when it comes to out-of-sample prediction, as prediction accuracy can be penalized by overfitting a model (Bakshi, Cao, and Chen, 1997). In order to compute the out-of-sample pricing errors, we first estimate the model over a rolling 5-day training period. Using the calibrated model, we then compute the 1-day ahead forecast of the option prices. With the 2012 S&P 500 Index put options data, there are 245 testing samples and a total of 14,794 pricing errors. Each training sample has on average 301 data points. The out-of-sample performance results are presented in Table 6. Table 6 Out-of-sample pricing RMSPEs and MAPEs in the main estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample spans 5 days (on average containing 301 data points) and each testing sample is 1 day ahead of the training sample (on average containing 60 data points). There are 245 testing days. Parameters are estimated using all options in the entire cross-section of the sample. Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 245 testing days collectively. The sample size is n  =  14,794. Table 6 Out-of-sample pricing RMSPEs and MAPEs in the main estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample spans 5 days (on average containing 301 data points) and each testing sample is 1 day ahead of the training sample (on average containing 60 data points). There are 245 testing days. Parameters are estimated using all options in the entire cross-section of the sample. Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 245 testing days collectively. The sample size is n  =  14,794. First, we focus on the out-of-sample performance of models in predicting the entire cross-section. The GH model emerges as the sole superior model overall in terms of both RMSPE (15.38%) and MAPE (11.04%). Unlike in the in-sample context, the GH model achieves greater out-of-sample accuracy than the t model, which attains an RMSPE of 16.33% (almost a full percentage point higher than the GH model) and an MAPE of 11.16% (0.12 percentage point higher than the GH model). Among the special cases of the GH model being considered (NIG, VG, and t), the t model is superior. The VG model obtains a 16.73% RMSPE, and the NIG model obtains a 16.71% RMSPE. It is evident that the GH model family dominates the FMLS model (18.98% RMSPE) and the BS model (19.94% RMSPE) in terms of overall out-of-sample performance. Next, we compare the out-of-sample performance of the models in different parts of the cross-section. The GH model is the best model for OTM, ATM, short-term, and long-term option categories. In particular, it outperforms the next best model by a wide margin for short-term options (by 2.97 percentage points) and for OTM options (by 1.45 percentage points). As in the in-sample setting, the t model fits the ITM options better than the GH model.10 For the medium-term option prices, the VG, t, and GH models perform similarly well. In all option categories, the NIG model is strictly dominated by the GH model. This result highlights the importance of freeing up the index parameter p for optimal out-of-sample performance. The FMLS model does not exhibit any advantage over the GH model family in forecasting the option prices in the given sample. It performs especially poorly in predicting the prices of short-term, OTM puts (43.01% RMSPE). This suggests that the heavy tail assumption under the FMLS model may have introduced unnecessary noise that hinders out-of-sample performance, at least during the sample period being considered. Interestingly, the BS model provides the best out-of-sample fit for long-term, OTM puts (14.34% RMSPE), implying that the additional complexity beyond normality as offered by the FMLS and GH models does not facilitate price prediction for this class of put options. 3.8 Robustness Analysis The above estimation results are based on the dataset from 2012 in which the U.S. stock market was relatively stable. Indeed, the VIX remains below 30 during the sample period. In this section, we want to investigate how the results will change under a more turbulent market condition. The robustness analysis will be based on a sample consisting of daily put option prices from August 4, 2011 to October 13, 2011. The period, which spans 50 trading days, was selected according to the criterion that the VIX Index exceeded and remained above the level 30 over the entire duration. The peak of VIX occurred at 48.00 on August 8, 2011, dubbed the “Black Monday” in the U.S. stock market as it reacted sharply to the Standard and Poor credit agency’s decision to downgrade the credit rating of the U.S. sovereign debt. The sample period is shown in Figure 4. Figure 4 View largeDownload slide The sample periods for main estimation and robustness estimation. Notes: The period used for main estimation spans over 2012 and is shaded in yellow. The period used for robustness estimation spans from August 4 to October 13 (50 trading days) and is shaded in red. The period used for robustness estimation was chosen such that the daily VIX (plotted) stays above 30 continuously over the period. Figure 4 View largeDownload slide The sample periods for main estimation and robustness estimation. Notes: The period used for main estimation spans over 2012 and is shaded in yellow. The period used for robustness estimation spans from August 4 to October 13 (50 trading days) and is shaded in red. The period used for robustness estimation was chosen such that the daily VIX (plotted) stays above 30 continuously over the period. In the robustness analysis, we revisit the out-of-sample model evaluation based on the new dataset. The sampling and estimation strategy are identical as before (using a 5-day rolling window and the cross-section as training samples to estimate the model and predict the next-day put option prices). The robustness estimation results are presented in Table 7 and summarized in Table 8. Given the superior out-of-sample performance of the flexible GH model in Section 3.7 (see Table 8) and to simplify the comparison, we do not repeat generating out-of-sample results for the special cases of the GH model in our robustness analysis. Table 7 Out-of-sample pricing RMSPEs and MAPEs in the robustness estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample is 5 days (on average containing 289 data points) and each testing sample is 1 day ahead of the training sample (on average containing 58 data points). There are 45 testing days. The sample size is n  =  2598. Table 7 Out-of-sample pricing RMSPEs and MAPEs in the robustness estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample is 5 days (on average containing 289 data points) and each testing sample is 1 day ahead of the training sample (on average containing 58 data points). There are 45 testing days. The sample size is n  =  2598. Table 8 A summary of the superior model(s) based on RMSPE by option category In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH Notes: A summary of the superior model(s) based on RMSPE in each of the option categories for the in-sample estimation and two out-of-sample predictions as reported in Tables 4, 6, and 7. For the out-of-sample (robustness) analysis, the comparison is done between only the BS, FMLS, and GH models, while for the in-sample and out-of-sample (main) analyses, the GH sub-models are also compared. See Table 2 for the definition of the option categories for moneyness and time-to-maturity. Table 8 A summary of the superior model(s) based on RMSPE by option category In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH Notes: A summary of the superior model(s) based on RMSPE in each of the option categories for the in-sample estimation and two out-of-sample predictions as reported in Tables 4, 6, and 7. For the out-of-sample (robustness) analysis, the comparison is done between only the BS, FMLS, and GH models, while for the in-sample and out-of-sample (main) analyses, the GH sub-models are also compared. See Table 2 for the definition of the option categories for moneyness and time-to-maturity. A number of observations are in order. First, the GH model remains the best model in terms of overall out-of-sample performance (15.21% RMSPE and 11.97% MAPE). Second, the GH model outperforms the other models in four individual option categories (OTM, short-term, medium-term, and long-term). In the ATM and ITM option groups, the FMLS model, which allows for infinite kurtosis in the log-return distribution, surpasses the GH model to be the best predictor. Third, focusing on the nine moneyness and time-to-maturity combinations, the GH model is the best model for OTM put options uniformly over all maturity groups, but the FMLS model performs uniformly the best for ITM puts; the model ranking is somewhat mixed for ATM puts across different maturities. From these observations, we learn that, among the models under consideration, not a single one of them eclipses the others in predicting every part of the cross-section of put option prices. The flexibility of the GH model contributes to the best overall out-of-sample performance and in particular to the prediction of the OTM put prices. Nevertheless, under turbulent market conditions, the advantage of the heavy-tailed log-return distribution under the FMLS model is manifested when predicting the ITM put prices, which are sensitive to the extreme left tail of the spot price distribution. We note that all models perform noticeably more poorly during the turbulent period in terms of out-of-sample RMSPE. Furthermore, the differential in the predictive ability of the more complex models relative to BS becomes narrower (e.g., during the normal sample period in the main analysis, the GH and FMLS models outperform the BS model overall by 4.56 and 0.96 percentage points, respectively, in out-of-sample RMSPE; the gaps shrink to 0.93 and 0.56 percentage points, respectively, during the turbulent period). The weaker pricing performance may be attributable to the stringent i.i.d. assumption, which is likely violated during turbulent market conditions. The heightened persistence of high volatilities as perceived by investors in turbulent times is associated with a higher option price; however, this is not directly modeled by the models we considered thus far. This motivates us to consider some of the possible SV extensions to the GH option pricing model, to be discussed in the following section. 4 SV Extensions So far, we have formulated and empirically tested the GH model by assuming stationary and independent increments in the spot price dynamics. In the literature, however, it was documented that the volatility of log-returns displays strong persistence and stochastic behaviors to the extent that the i.i.d. assumption is rejected empirically. This suggests that a GH model that could account for serial dependence in volatility is likely to reduce the pricing errors as observed in the previous empirical study. Building on the time-changed Lévy process modeling framework in Section 1, it is possible to introduce SV in a similar spirit of Carr et al. (2003), Huang and Wu (2004), and Carr and Wu (2004). We start with the baseline GH model, in which the source of randomness is the subordinated Brownian motion B(g(t)), where g(t) is the GIG process. Next, noting that B(g(t)) is a (GH) Lévy process with stationary and independent increments, we may construct a time-changed Lévy process as in Equation (1) by setting L(t)=B(g(t)) and define a time deformation process T(t) that satisfies the three conditions in Section 1. The resultant time-changed Lévy process Z(t)=L(T(t))—a subordinated GH process—is then treated as the random source that drives the log-return process in Equation (3). To induce stochastic volatilities, the process T(t) should exhibit time-dependent dynamics. One way to achieve this is to set T(t) to be the integral of a mean-reverting process. More precisely, it is given by T(t)=∫0th(u)du,dh(t)=κ[1−h(t)]dt+λh(t)dW˜(t). (11) The process h(t) represents the activity rate process which determines the speed at which the business clock runs. It is modeled by the square-root mean-reverting process that solves the stochastic differential equation (11), where κ controls the rate of mean reversion, λ represents the volatility of the activity rate process, and W˜(t) is a standard Wiener process possibly correlated with the standard Wiener process W(t) in Equation (7).11 The long-run activity rate is set to one so that the normalizing restriction on T(t) is satisfied. The characteristic function of T(t) is given by CGMY (2003) φTt(u)=exp(κ2tλ2)+2iuh(0)κ+ψ coth(ψt/2)[cosh⁡(ψt2)+κψsinh⁡(ψt2)]−2κλ2, where ψ=κ2−2λ2iu. This allows us to obtain the characteristic function of the log-price using Equations (2) and (4). The European option prices under the SV flexible GH model (GH–SV model) can then be calculated using Equations (5) and (6). There is an alternative modeling strategy that allows for long-range dependence in SV within the subordinated Brownian motion framework introduced in Section 1.12 Under the GH option pricing model of Finlay and Seneta (2012), serial dependence of stock returns is modeled by a background-driving Lévy process (BDLP) (Barndorff-Nielsen and Shephard, 2001) construction on the subordinator. More specifically, the subordinator g(t) is defined as a superposition of OU processes, each of which is driven by an independent BDLP and has a GIG marginal distribution. The important feature is that the increments of g(t) form an autocorrelated sequence whose marginal distribution remains to be GIG by self-decomposability (Halgreen, 1979). The autocorrelation in the increments of g(t) may induce short- or long-range dependence in squared stock returns. Comparing Finlay and Seneta’s model to our extended GH model, their European option pricing formula involves two-dimensional integrations. We also note that even in the i.i.d. case, the two models have different martingale adjustments and different parameterizations under the risk-neutral measure.13 One future direction of research is to reconcile the models and investigate their empirical performance in the presence of SV. 5 Conclusion Insight to the risk-neutral distribution of logarithmic stock returns is vital to the fitting and prediction of option prices. In this paper, we develop a flexible, four-parameter GH option pricing model in the time-changed Lévy process framework, and present six three-parameter special cases: the VG, t, H, RH, NIG, and NRIG option pricing models. In respect of the seven models’ properties, the flexible GH model and its special cases generalize the BS model by allowing the passage of economic time to be stochastic. As such, the GH model family forms a class of time-subordinated models that can cope with yet another facet of the unpredictable financial market. In addition, the subordination to a drifted Brownian motion entails that the class of GH processes are able to capture excess kurtosis and skewness. Using the characteristic function approach of Carr and Madan (1999), we obtain an explicit pricing formula for European options under the baseline GH model. The option pricing formula takes a convenient semi-closed form. With judicious choices of parameterization and variable transformation, all the four model parameters can be estimated freely. These features greatly facilitate the empirical study on the GH option pricing model. Using S&P 500 Index options, we conduct an extensive empirical study on the GH option pricing model and some of its special cases, and compare the empirical performance against some benchmark option pricing models that assume i.i.d. spot returns, including the BS model and the FMLS model (Carr and Wu, 2003). Our empirical analysis suggests that the GH model generally delivers the best in-sample fit and out-of-sample prediction for all put option prices in the cross-section. The same conclusion holds for the robustness analysis in which the sample spans a more turbulent period, although the heavy-tail FMLS option pricing model yields more accurate out-of-sample forecast for some specific option categories. Among the NIG, VG, and t special cases, the t model seems to be the most tenable model for pricing options. Remarkably, the t model’s average in-sample fit is better than that of the VG model (Madan, Carr, and Chang, 1998) for all option types—a result which is corroborated by the GH model’s parameter estimates. In this paper, we concentrate on the empirical study of the baseline GH model, which belongs to a class of the Lévy-based option pricing models that assumes independent and stationary increments in the underlying spot price process. However, the i.i.d. assumption on the spot returns can be crippling in view of the bulk of well-known stylized facts regarding financial asset returns. Leveraging the proliferous time-changed Lévy process framework (CGMY, 2003; Carr and Wu, 2004; Huang and Wu, 2004), we present an SV extension to the baseline model that employs a mean-reverting process for time deformation. Empirical study of this extension is intended for future work. Appendix A: The Characteristic Function of the Logarithmic Stock Price First, we derive the value of the martingale correction factor, denoted by ωt, as −log(φZt(−i)). Assuming no arbitrage and by risk-neutral pricing, the discounted stock price process {Ste−(r−q)t} follows a martingale under the risk-neutral probability measure Q. This means that, for any positive t, S0=EQ[Ste−(r−q)t|F0], (A.1) where F0 is the information set at time 0, and E Q[·|F0] represents the conditional expectation taken with respect to Q. Given the log stock price model in Equation (3), the martingale condition is satisfied by applying a martingale correction factor ωt to the drift of the log-return process, as follows: S0=EQ[S0 exp{(r−q)t+ωt+Zt−(r−q)t}|F0]. (A.2) Solving for ωt, we obtain: 1=EQ[exp{ωt+Zt}|F0]1=eωtEQ[eZt|F0]e−ωt=EQ[eZt|F0]ωt=−log(EQ[eZt|F0])ωt=−log(φZt(−i)) (A.3) This martingale correction factor is the same as the one in CGMY (2003, equation (4.7), p. 358) using ordinary exponential transform. We then obtain Equation (4), by expressing the characteristic function of the log spot price φ log St(u) in terms of φZt(u), the characteristic function of Z(t), as follows: φ log St(u)=EQ[exp{iu log(St)}|F0]=EQ[Stiu|F0]=EQ[S0iu exp{iu[(r−q)t−log(φZt(−i))+Zt]}]=S0iu exp{iu(r−q)t−log(φZt(−i))}φZt(u). (A.4) Appendix B: Derivations for the Flexible GH, VG, and t Option Pricing Models B.1 The Flexible GH Model B.1.1 Characteristic function used for option pricing The characteristic function of the GIG distribution is given by φg1(u)=(γ2γ2−2iu)p2Kp(δ2(γ2−2iu))Kp(γδ)=(1−2γ2iu)−p2Kp(ζ1−2γ2iu)Kp(ζ), (A.5) where ζ=γδ, γ2 is given by Equation (8) and Kj(·) is the modified Bessel function of the third kind with index j. Using Equation (2), we obtain φZ1(u)=φg1(uθ+iσ2u22)=(1−2γ2(iuθ−12σ2u2))−p2Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ). (A.6) It follows from Equation (4) that φ log St(u)=S0iu exp{iu(r−q+ω)t}φX1(u)t=S0iu exp{iu(r−q+ω)t}[1−2γ2(iuθ−12σ2u2)]−pt2×[Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ)]t, (A.7) where the unit time drift adjustment can be computed as ω=−log φX1(−i)=−p2 log[(1−2γ2(θ+12σ2))]+log[Kp(ζ1−2γ2(θ+12σ2))Kp(ζ)], (A.8) for θ<(γ22−σ22). B.1.2 Distributional properties of Z1 under the GH model Let us start by studying the statistical properties of Z1. First, we derive the probability density function of the GH distribution, which is the distribution of the random variable Z1 in the log-price specification Equation (3) under the flexible GH option pricing model. Let N (x;μ,σ2) denote the density of a normal distribution with mean μ and variance σ2. Suppose, conditional on g1, that X1 follows a normal distribution with mean θg1 and variance σ2g1, and that g1 follows a GIG with probability density function fg1(g), given by GIG(p,δ,γ):                 fg1(g)=(γδ)p2Kp(γδ)gp−1 exp{−12(γ2g+δ2g)}, (A.9) where Kj(·) is the modified Bessel function of the third kind with index j. Then, Z1 will follow a GH distribution with its density function fZ1(x) given by, for any real value x, fZ1(x)=∫0∞N(x;θg,σ2g)fg1(g)dg=∫0∞12πσ2g exp{−(x−θg)22σ2g}·(γδ)p2Kp(γδ)gp−1 exp{−12(γ2g+δ2g)}dg=(γδ)p22πσ2Kp(γδ)∫0∞gp−32 exp{−12[(x−θg)2σ2g+γ2g+δ2g]}dg (A.10) =(γδ)p22πσ2Kp(γδ)eθσ2x∫0∞gp−32 exp{−12[(θ2σ2+γ2)g+(x2σ2+δ2)1g]}dg=(γδ)p2πσ2Kp(γδ)eθσ2x(x2σ2+δ2θ2σ2+γ2)p2−14Kp−12((θ2σ2+γ2)(x2σ2+δ2)). The last equality uses the fact that, for all η,υ>0, ∫0∞xj−1e−12(ηx+υx)dx=2(υη)j2Kj(ηυ). The parameter restriction Equation (8) that γ2=ζKp+1(ζ)Kp(ζ) (where ζ=γδ) coming from the normalization E [g1]=1 still holds under the time-changed Lévy process framework. The first four standardized central moments of Z1 under the GH model are given by E(Z1)=θδγKp+1(δγ)Kp(δγ),Var(Z1)=σ2δγKp+1(δγ)Kp(δγ)+θ2δ2γ2[Kp+2(δγ)Kp(δγ)−(Kp+1(δγ)Kp(δγ))2],skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2{θ3δ3γ3[Kp+3(δγ)Kp(δγ)−3Kp+2(δγ)Kp+1(δγ)Kp2(δγ)+2(Kp+1(δγ)Kp(δγ))3]+3σ2θδ2γ2[Kp+2(δγ)Kp(δγ)−(Kp+1(δγ)Kp(δγ))2]},kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=Var(Z1)−2{θ4δ4γ4[Kp+4(δγ)Kp(δγ)−4Kp+3(δγ)Kp+1(δγ)Kp2(δγ)+6Kp+2(δγ)Kp+12(δγ)Kp3(δγ)−3(Kp+1(δγ)Kp(δγ))4]+σ2θ2δ3γ3[6Kp+3(δγ)Kp(δγ)−12Kp+2(δγ)Kp+1(δγ)Kp2(δγ)+6(Kp+1(δγ)Kp(δγ))3]+3σ4δ2γ2Kp+2(δγ)Kp(δγ)}. (A.11) Applying the normalization E[g1]=1, and denoting ζ=δγ, we obtain the parameter restriction γ2=ζKp+1(ζ)Kp(ζ). The four standardized central moments of Z1 then become E(Z1)=θ,Var(Z1)=σ2+θ2[Kp+2(ζ)Kp(ζ)Kp+12(ζ)−1],skew(Z1)=Var(Z1)−3/2{θ3[Kp+3(ζ)Kp2(ζ)Kp+13(ζ)−3Kp+2(ζ)Kp(ζ)Kp+12(ζ)+2]+3σ2θ[Kp+2(ζ)Kp(ζ)Kp+12(ζ)−1]},kur(Z1)=Var(Z1)−2{θ4[Kp+4(ζ)Kp3(ζ)Kp+14(ζ)−4Kp+3(ζ)Kp2(ζ)Kp+13(ζ)+6Kp+2(ζ)Kp(ζ)Kp+12(ζ)−3]+σ2θ2[6Kp+3(ζ)Kp2(ζ)Kp+13(ζ)−12Kp+2(ζ)Kp(ζ)Kp+12(ζ)+6]+3σ4Kp+2(ζ)Kp(ζ)Kp+12(ζ)}. (A.12) Note that in the symmetric case of θ = 0, we have kur(Z1)>3, by the Turán-type inequality Kp+2(y)Kp(y)>Kp+12(y) (Ismail and Muldoon, 1978). B.1.3 Distributional properties of Zt under the GH model Next, let us study the statistical properties of Zt. From the infinitely divisibility of a general Lévy process Zt, the characteristic functions of Zt and Z1 are related by φZt(u)=[φZ1(u)]t. It follows that their cumulant generating functions (i.e., the log-characteristic functions) are related by ΨZt(u)=tΨZ1(u). By repeated differentiation, the following relations regarding the cumulants of Zt and Z1 hold for all positive integers n: cn(Zt):=1in∂n∂unΨZt(0)=t1in∂n∂unΨZ1(0)=:tcn(Z1). Using the links between the cumulants and the central moments, the first four central moments of Zt and Z1 can be related as follows: E(Zt)=c1(Zt)=tc1(Z1)=tE(Z1),Var(Zt)=c2(Zt)=tc2(Z1)=tVar(Z1),E[(Zt−E(Zt))3]=c3(Zt)=tc3(Z1)=tE[(Z1−E(Z1))3],E[(Zt−E(Zt))4]=c4(Zt)+3c2(Zt)2=tc4(Z1)+3[tc2(Z1)]2=tE[(Z1−E(Z1))4]+3t(t−1)Var(Z1)2. (A.13) The skewness and kurtosis of Zt are then given by: skew(Zt)=E[(Zt−E(Zt))3]Var(Zt)3/2=tE[(Z1−E(Z1))3][tVar(Z1)]3/2=1tskew(Z1),kur(Zt)=E[(Zt−E(Zt))4]Var(Zt)2=tE[(Z1−E(Z1))4]+3t(t−1)Var(Z1)2[tVar(Z1)]2=1tkur(Z1)+3(1−1t). (A.14) We thus see that, for a general Lévy process, its mean and variance increase with t, its skewness diminishes with t at the square-root rate, and its kurtosis approaches that of a normal distribution in the limit as t→∞. In the symmetric case of θ = 0, we see that kur(Z1)>3, which implies that kur(Zt)>3 for all finite t. B.2 The VG Model and its Statistical Properties B.2.1 Characteristic function used for option pricing The characteristic function of the gamma distribution is given by φg1(u)=(1−iuν)−1ν such that φX1(u)=[1−ν(iuθ−12σ2u2)]−1ν. (A.15) It follows then that the VG option price can be computed using φ log St(u)=S0iu exp{iu(r−q+ω)t}[1−ν(iuθ−12σ2u2)]−tν, (A.16) with unit time drift adjustment, ω=1ν log[1−ν(θ+12σ2)], for θ<(1ν−σ22). B.2.2 Distributional properties of Z1 under the VG model Let us study the statistical properties of Z1. The density function of Z1 under the VG model can be obtained from Equation (A.10) by letting δ→0 and setting p=1ν>0, so that γ=2p=2ν. fZ1(x)=(γδ)1ν2πσ2Γ(1ν)21ν−1(γδ)−1νeθσ2x(x2σ2θ2σ2+2ν)12ν−14K1ν−12((θ2σ2+2ν)x2σ2)=γ2ν2πσ2Γ(1ν)21ν−1eθσ2x(x2θ2+2σ2ν)12ν−14K1ν−12(xσ2θ2+2σ2ν). (A.17) This recovers the density function in Madan, Carr, and Chang (1998, equation (23)). Imposing the same limits on the parameters, we obtain the first four standardized central moments under the VG model for ν>0: E(Z1)=θ,Var(Z1)=σ2+θ2ν,skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2(2θ3ν2+3σ2θν)kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=3Var(Z1)−2[θ4ν2(1+2ν)+2σ2θ2ν(1+2ν)+σ4(1+ν)]. (A.18) In the symmetric case of θ = 0, the kurtosis is always greater than 3 as ν>0. B.3 The t Model B.3.1 Characteristic function used for option pricing The characteristic function of the reciprocal gamma distribution, R Γ(a,b) is given by φg1(u)=2(−biu)a/2Γ(a)Ka[−4biu]=2(1ν−1)12νΓ(1ν)(−iu)12νK1ν[−4(1ν−1)iu], (A.19) where a=1ν and b=1ν−1 for 0<ν<1. Kj(·) is the modified Bessel function of the third kind with index j. Using the same procedure as for the GH and VG derivations, the characteristic function of X1 is given by φX1(u)=2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)]. (A.20) The characteristic function of log St is thus given by φ log St(u)=S0iu exp{iu(r−q+ω)t}×[2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν(−4(1ν−1)(iuθ−12σ2u2))]t, (A.21) and the unit time drift adjustment, ω=−log[2(1ν−1)12νΓ(1ν)[−(θ+12σ2)]12νK1ν(−4(1ν−1)(θ+12σ2))], (A.22) for θ<−σ22. B.3.2 Distributional properties of Z1 under the t model The density function of Z1 under the t model can be obtained from Equation (A.10) by letting γ→0 and setting p=−1ν<0, so that δ=2(−p−1)=2(1ν−1). fZ1(Z)=(1ν−1)1νσ2πΓ(1ν)eθσ2Z(Z2+2(1ν−1)σ2θ)−1ν−12K1ν+12(θσ2Z2+2(1ν−1)σ2). (A.23) By the variable and parameter transformation, y=Z/σ+μ, β=θ/σ, α=β2+γ2, λ=−1/ν, and δ=2(1ν−1), we recover the density function of the GH–skew–t distribution in Aas and Haff (2006, equation (3)). Imposing the same limits on the parameters, we obtain the first four standardized central moments under the t model and the associated domains of ν in which the moments exist: E(Z1)=θ,Var(Z1)=σ2+θ2ν1−2ν for ν∈(0,12),skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2[θ34ν2(1−2ν)(1−3ν)+3σ2θν1−2ν] for ν∈(0,13),kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=3Var(Z1)−2[θ4ν2(1+5ν)(1−2ν)(1−3ν)(1−4ν)+2σ2θ2ν(1+ν)(1−2ν)(1−3ν)+σ41−ν1−2ν] for ν∈(0,14). (A.24) In the symmetric case of θ = 0, the kurtosis is always greater than 3 as ν>0. Appendix C: The VG and t Limiting Cases of the GH Model The following lemma shows how we may obtain the characteristic functions of Z1 under the VG and t models as appropriate limits of that of Z1 under the GH model. Lemma 1 Suppose p=1ν>0, δ→0 and γ=2ν. Then, the GH model reduces to the VG model. Suppose p=−1ν<−1, γ→0 and δ=2(1ν−1). Then, the GH model reduces to the t model. Proof Either δ→0 or γ→0 implies that ζ→0. Making use of the properties of the modified Bessel function of the third kind with index j, Kj(·), that Kj(w)∼Γ(|j|)2|j|−1w−|j|    for w↓0, and that Kj(w)=K−j(w), the characteristic function of X1 under the GH model (given by Equation (A.6)) can be further simplified as follows: Under the assumptions in (i), we have φX1(u)=(1−ν(iuθ−12σ2u2))−12νlim⁡δ→0Γ(1ν)21ν−1(γδ1−ν(iuθ−12σ2u2))−1νΓ(1ν)21ν−1(γδ)−1ν=(1−ν(iuθ−12σ2u2))−12ν(1−ν(iuθ−12σ2u2))−12ν=(1−ν(iuθ−12σ2u2))−1ν, which is the corresponding characteristic function under the VG model, given in Equation (A.15). Under the assumptions in (ii), we have φX1(u)=lim⁡γ→0(1−2γ2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))K1ν(δγ)=lim⁡γ→0(1−2γ2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))Γ(1ν)21ν−1(γδ)−1ν=lim⁡γ→0δ1νΓ(1ν)21ν−1(γ2−2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))=δ1νΓ(1ν)21ν−1(−2(iuθ−12σ2u2))12νK1ν(δ−2(iuθ−12σ2u2)) =212ν(1ν−1)12νΓ(1ν)21ν−1212ν[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)]=2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)], which is the corresponding characteristic function under the t model, given in Equation (A.20). Footnotes * We deeply thank Fabio Trojani, the co-editor, and two anonymous referees for their constructive suggestions which led to substantial improvement on the paper. We also thank those who gave valuable comments to our work at the EcoSta2017 Conference and at a University of Sydney Business School seminar. All remaining errors are ours. 1 The pure-jump and infinite activity property of the GH process and special cases is a distinguishing property from Heston’s (1993) model, the BS model, and from other seminal models, such as Merton’s (1976) jump-diffusion model, Bates, (1996) model, and Pan’s (2002) model, as noted by Carr and Wu (2004). As a purely discontinuous process, the GH process is also different from the GH diffusion process of Bibby and Sørensen (1997) and Rydberg (1999). 2 The use of time deformation processes that induce persistent variations in drift and volatility is supported by a general equilibrium argument. For example, see Shaliastovich and Tauchen (2011). We thank one referee who brings this literature to our attention. 3 Admitting the computational difficulty in estimating the four parameters jointly, Prause (1999) only presents the risk-neutral estimation results for the NIG model, a special case of the GH option pricing model obtained by fixing the index parameter [or shape parameter according to Prause (1999)] to p=−1/2 [see Table 2.28 and the discussion on p. 61 of Prause (1999) for details]. Nevertheless, Eberlein and Prause (2002) and Prause (1999) estimate the GH model parameters under the physical measure using spot data. 4 The FMLS model of Carr and Wu (2003) features a heavy-tailed, negatively skewed log-return distribution. While not a time-changed Lévy process model, the FMLS model is chosen as a benchmark model rather than the log-stable model of Hurst, Platen, and Rachev (1997) because the former guarantees finite moments of the spot price and hence finite option prices. 5 In general, when there are a finite number of assets in the market, the risk-neutral probability measure Q may not be unique. When the price of the underlying asset exhibits jumps, it typically requires infinite many assets to complete the market (an exception is when the asset price contains Poisson jumps with a finite number of different jump sizes). On the technical level, sufficient and necessary conditions for market completeness are known for a general Lévy process (Sato, 1999). More recently, Li and Mendoza-Arriaga (2015) provide a set of sufficient conditions for equivalent measure changes for a subordinated diffusion model. 6 A similar version of the t option pricing model was proposed in Yeap (2014), except that a different parameterization ν=Var(g1) was adopted. In this paper, we enforce the normalization E[g1]=1 so that Var(g1)=ν1−2ν. Our model allows for skewness in the log-return distribution, whereas a symmetric t option pricing model has been proposed (Cassidy, Hamp, and Ouyed, 2010). 7 The log-returns (net of the risk-free rate and dividend yield) follow a symmetric normal distribution under BS and a negatively skewed stable distribution under the FMLS model (with skewness parameter −1). 8 Alternatively, we may smooth out the cyclicity by using a dynamically weighted sample in which different weights are assigned to options with different time-to-maturities (Pan, 2002). The weights are determined such that the average time-to-maturity of the weighted sample is maintained at a constant level. The advantage of our current method (by retaining the raw data) is that the term structure of the risk-neutral distribution, and hence the relationship between the shape parameters and the time-to-maturity becomes more apparent. 9 In particular, the empirical skewness and kurtosis of log-returns as implied from the mean parameter estimates under the t model are equal to −0.25 and 4.04, respectively, which are both bigger in magnitude than those under the VG model (−0.09 and 3.33). The empirical mean and variance are similar under the two models (−0.04 and 0.02 under the VG model; −0.05 and 0.03 under the t model). 10 We note that some pricing improvements are observed in the out-of-sample context compared with the in-sample context for ITM options (e.g., ITM option prices as fit by the VG model). This can be explained by the relatively small representation of ITM options (16%) in the data compared with ATM (41%) and OTM (43%) options (see Table 2), such that the ITM options have a relatively weak influence over the in-sample fitting. 11 We may set W˜(t)=ρW(t)+1−ρ2W⌣(t), where W(t) and W⌣(t) are independent standard Wiener processes. 12 The option pricing models for assets with long-range dependence were pioneered by Heyde (1999), and Heyde and Liu (2001). Special cases of the GH distribution have since been employed in short- and long-range dependence models, such as the t distribution by Heyde and Leonenko (2005); Finlay and Seneta (2006); and Leonenko, Petherick, and Sikorskii (2011); the VG distribution by Finlay and Seneta (2006) and Leonenko, Petherick, and Sikorskii (2012b); and the NIG distribution by Leonenko, Petherick, and Sikorskii (2012a,b). 13 Unlike our model, Finlay and Seneta (2012) do not impose the normalization restriction E[g1]=1 on the subordinator. Instead, their martingale condition implies a restriction on the skewness parameter θ=−12σ2, while θ is a free parameter under our model. References Aas K. , Haff I. H. . 2006 . The Generalised Hyperbolic Skew–Student’s t Distribution . Journal of Financial Econometrics 4 : 275 – 309 . Google Scholar CrossRef Search ADS Andersen T. G. , Fusari N. , Todorov V. . 2015 . The Risk Premia Embedded in Index Options . Journal of Financial Economics 117 : 558 – 584 . Google Scholar CrossRef Search ADS Bakshi G. , Cao C. , Chen Z. . 1997 . Empirical Performance of Alternative Option Pricing Models . Journal of Finance 52 : 2003 – 2049 . Google Scholar CrossRef Search ADS Bakshi G. , Madan D. . 2000 . Spanning and Derivative-Security Valuation . Journal of Financial Economics 55 : 205 – 238 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. 1977 . Exponentially Decreasing Distributions for the Logarithm of Particle Size . Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 353 : 401 – 419 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Halgreen C. . 1977 . Infinite Divisibility of the Hyperbolic and Generalised Inverse Gaussian Distributions . Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 38 : 309 – 311 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Shephard N. . 2001 . Non-Gaussian OU Based Models and some of their Uses in Financial Economics . Journal of Royal Statistical Society: Series B 63 : 167 – 241 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Mikosch T. , Resnick S. . 2001 . Lévy Processes—Theory and Applications . Boston : Birkhauser . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Shephard N. . 2012 . “Basics of Lévy Processes.” Economics Series Working papers, University of Oxford (No. 610) . Bates D. S. 1996 . Jumps and Stochastic Volatility: Exchange Rate Processes Implicit in Deutsche Mark Options . Review of Financial Studies 9 : 69 – 107 . Google Scholar CrossRef Search ADS Bibby B. M. , Sørensen M. . 1997 . A Hyperbolic Diffusion Model for Stock Prices . Finance & Stochastics 1 : 25 – 41 . Google Scholar CrossRef Search ADS Black F. , Scholes M. . 1973 . The Pricing of Options and Corporate Liabilities . Journal of Political Economy 81 : 637 – 654 . Google Scholar CrossRef Search ADS Carr P. , Madan D. . 1999 . Option Valuation Using Fast Fourier Transform . Journal of Computational Finance 2 : 61 – 73 . Google Scholar CrossRef Search ADS Carr P. , Wu L. . 2003 . The Finite Moment Log Stable Process and Option Pricing . Journal of Finance 58 : 753 – 777 . Google Scholar CrossRef Search ADS Carr P. , Wu L. . 2004 . Time-Changed Lévy Processes and Option Pricing . Journal of Financial Economics 71 : 113 – 141 . Google Scholar CrossRef Search ADS Carr P. , Geman H. , Madan D. B. , Yor M. . 2003 . Stochastic Volatility for Lévy Processes . Mathematical Finance 13 : 345 – 382 . Google Scholar CrossRef Search ADS Cassidy D. T. , Hamp M. J. , Ouyed R. . 2010 . Pricing European Options with a Log Student’s t-distribution: A Gosset Formula . Physica A 389 : 5736 – 5748 . Google Scholar CrossRef Search ADS Christie A. A. 1982 . The Stochastic Behavior of Common Stock Variances: Value, Leverage and Interest Rate Effects . Journal of Financial Economics 10 : 407 – 432 . Google Scholar CrossRef Search ADS Christoffersen P. , Jacobs K. . 2004 . The Importance of the Loss Function in Option Valuation . Journal of Financial Economics 72 : 291 – 318 . Google Scholar CrossRef Search ADS Clark P. K. 1973 . A Subordinated Stochastic Process Model with Finite Variance for Speculative Prices . Econometrica 41 : 135 – 155 . Google Scholar CrossRef Search ADS Daal E. A. , Madan D. B. . 2005 . An Empirical Examination of the Variance–Gamma model for Foreign Currency Options . The Journal of Business 78 : 2121 – 2152 . Google Scholar CrossRef Search ADS Eberlein E. , Keller U. . 1995 . Hyperbolic Distributions in Finance . Bernoulli 1 : 281 – 299 . Google Scholar CrossRef Search ADS Eberlein E. , Prause K. . 2002 . “The Generalized Hyperbolic Model: Financial Derivatives and Risk Measures.” In Mathematical Finance—Bachelier Congress 2000 , pp. 245 – 267 . Heidelberg : Springer-Verlag . Google Scholar CrossRef Search ADS Finlay R. , Seneta E. . 2006 . Stationary-Increment Student and Variance–Gamma Processes . The Journal of Applied Probability 43 : 441 – 453 . Google Scholar CrossRef Search ADS Finlay R. , Seneta E. . 2012 . A Generalized Hyperbolic Model for a Risky Asset with Dependence . Statistics and Probability Letters 82 : 2164 – 2169 . Google Scholar CrossRef Search ADS Gilli M. , Schumann E. . 2012 . Calibrating Option Pricing Models with Heuristics . Natural Computing in Computational Finance 380 : 9 – 37 . Google Scholar CrossRef Search ADS Gruber P. H. , Tebaldi C. , Trojani F. . 2015 . “The Price of the Smile and Variance Risk Premia.” Swiss Finance Institute Research Paper No. 15-36 . Halgreen C. 1979 . Self-decomposability of the Generalized Inverse Gaussian and Hyperbolic Distributions . Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 47 : 13 – 17 . Google Scholar CrossRef Search ADS Heston S. L. 1993 . A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options . The Review of Financial Studies 6 : 327 – 343 . Google Scholar CrossRef Search ADS Heyde C. C. 1999 . A Risky Asset Model with Strong Dependence through Fractal Activity Time . The Journal of Applied Probability 36 : 1234 – 1239 . Google Scholar CrossRef Search ADS Heyde C. C. , Leonenko N. N. . 2005 . Student Processes . Advances in Applied Probability 37 : 342 – 365 . Google Scholar CrossRef Search ADS Heyde C. C. , Liu S. . 2001 . Empirical Realities for a Minimal Description Risky Asset Model. The Need for Fractal Features . Journal of Korean Mathematical Society 38 : 1047 – 1059 . Huang J. , Wu L. . 2004 . Specification Analysis of Option Pricing Models Based on Time-Changed Levy Processes . The Journal of Finance 59 : 1405 – 1439 . Google Scholar CrossRef Search ADS Hurst S. R. , Platen E. , Rachev S. . 1997 . Subordinated Market Index Models: A Comparison . Financial Engineering and the Japanese Markets 4 : 97 – 124 . Google Scholar CrossRef Search ADS Ismail M. , Muldoon M. . 1978 . Monotonicity of the Zeros of a Cross-Product of Bessel Functions . SIAM Journal on Mathematical Analysis 9 : 759 – 767 . Google Scholar CrossRef Search ADS Jørgensen B. 1982 . Statistical Properties of the Generalised Inverse Gaussian Distribution . New York : Springer-Verlag . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2011 . The Student Subordinator Model with Dependence for Risky Asset Returns . Communications in Statistics—Theory and Methods 40 : 3509 – 3523 . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2012a . A Normal Inverse Gaussian Model for a Risky Asset with Dependence . Statistics and Probability Letters 82 : 109 – 115 . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2012b . Fractal Activity Time Models for Risky Asset with Dependence and Generalised Hyperbolic Distributions . Stochastic Analysis and Applications 30 : 476 – 492 . Google Scholar CrossRef Search ADS Li L. , Mendoza-Arriaga R. . 2015 . “Equivalent Measure Changes for Subordinate Diffusions.” Working paper . Madan D. B. , Carr P. P. , Chang E. C. . 1998 . The Variance Gamma Process and Option Pricing . European Finance Review 2 : 79 – 105 . Google Scholar CrossRef Search ADS Madan D. B. , Milne F. . 1991 . Option Pricing with VG Martingale Components . Mathematical Finance 1 : 39 – 55 . Google Scholar CrossRef Search ADS Merton R. C. 1973 . The Theory of Rational Option Pricing . Bell Journal of Economics and Management Science 4 : 141 – 183 . Google Scholar CrossRef Search ADS Merton R. C. 1976 . Option Pricing When Underlying Stock Returns Are Discontinuous . Journal of Financial Economics 3 : 125 – 144 . Google Scholar CrossRef Search ADS Pan J. 2002 . The Jump-Risk Premium Implicit in Options: Evidence from an Integrated Time-Series Study . Journal of Financial Economics 63 : 3 – 50 . Google Scholar CrossRef Search ADS Prause K. 1999 . The Generalized Hyperbolic Model: Estimation, Financial Derivatives, and Risk Measures . Thesis, Albert–Ludwigs–Universitẗ Freiburg . Rubinstein M. 1985 . Non-Parametric Tests of Alternative Option Pricing Models using All Reported Trades and Quotes on the 30 Most Active CBOE Option Classes from August 23, 1976 through August 31, 1978 . Journal of Finance 40 : 455 – 480 . Google Scholar CrossRef Search ADS Rydberg T. 1999 . Generalized Hyperbolic Diffusion Processes with Applications in Finance . Mathematical Finance 9 : 183 – 201 . Google Scholar CrossRef Search ADS Scott L. 1997 . Pricing Stock Options in a Jump-Diffusion Model with Stochastic Volatility and Interest Rates: Applications of Fourier Inversion Methods . Mathematical Finance 7 : 413 – 426 . Google Scholar CrossRef Search ADS Sato K. 1999 . Lévy Processes and Infinitely Divisible Distributions . Cambridge : Cambridge University Press . Shaliastovich I. , Tauchen G. . 2011 . Pricing of the Time-Change Risks . Journal of Economic Dynamics and Control 35 : 843 – 858 . Google Scholar CrossRef Search ADS White H. 1980 . A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity . Econometrica 48 : 817 – 838 . Google Scholar CrossRef Search ADS Yeap C. 2014 . The Skew-t Option Pricing Model . Thesis , The University of Sydney Business School . © The Author, 2017. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Financial Econometrics Oxford University Press

A Flexible Generalized Hyperbolic Option Pricing Model and Its Special Cases

Loading next page...
 
/lp/ou_press/a-flexible-generalized-hyperbolic-option-pricing-model-and-its-special-uz00hbnN1J
Publisher
Oxford University Press
Copyright
© The Author, 2017. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
1479-8409
eISSN
1479-8417
D.O.I.
10.1093/jjfinec/nbx030
Publisher site
See Article on Publisher Site

Abstract

Abstract We study a generalized hyperbolic (GH) time-changed Lévy process for option pricing and examine six three-parameter special cases: the variance gamma (VG) model of Madan, Carr, and Chang (1998), t, hyperbolic (H), normal inverse Gaussian (NIG), reciprocal hyperbolic (RH), and normal reciprocal inverse Gaussian (NRIG) option pricing models. We study the GH model’s moment properties of the associated risk-neutral distribution of logarithmic spot returns, and obtain an explicit pricing formula for European options facilitated by the time-change Lévy process construction. Using S&P 500 Index European options during low and high volatility sample periods, we compare the GH model empirically with existing benchmark models such as the finite-moment log-stable model and the Black–Scholes model. The GH model offers the best in- and out-of-sample performance overall, and a proposed t model special case generally outperforms the existing VG special case. We also present a stochastic volatility extension of the GH model. Empirical option prices indicate that the likelihood of extreme logarithmic stock returns is higher than that implied by the Black–Scholes (BS) model. Option prices also reveal that market participants pay more to protect themselves from losses than to pursue gains of equivalent magnitude. The implication is that the risk-neutral distribution of log-returns exhibits excess kurtosis and negative skewness (Madan and Milne, 1991; Eberlein and Keller, 1995). These two digressions from the normality assumption (Black and Scholes, 1973; Merton, 1973) are in part responsible for the poor empirical pricing performance of the BS model. To combat this deficiency, the generalized hyperbolic (GH) distribution may be used to improve option pricing as it accommodates skewness and thicker, semi-heavy tails (Barndorff-Nielsen, 1977). In this paper, we focus on addressing the challenges faced by GH option pricing in the setting in which spot returns are independent. Our contribution to this field is to propose a new form of the GH option pricing model, the flexible GH model, which contains four free parameters and can be conveniently estimated. In addition, we present six three-parameter option pricing models as special cases of the flexible GH model, including the variance gamma (VG), t, hyperbolic (H), normal inverse Gaussian (NIG), reciprocal hyperbolic (RH), and normal reciprocal inverse Gaussian (NRIG) models. With the exception of the VG model proposed by Madan, Carr, and Chang (1998), the remaining five option pricing models are innovations of this paper. To construct the flexible GH option pricing model, we generalize the BS model by the time change (or subordination) method (Clark, 1973). More specifically, the risk-neutral dynamics of the log-returns is modeled by a drifted Brownian motion subordinated by a stochastic time deformation process, which is defined to be a generalized inverse Gaussian (GIG) process (Barndorff-Nielsen and Halgreen, 1977). This construction generates a class of pure-jump, infinite activity processes (Barndorff-Nielsen and Shephard, 2001).1 A typical sample path of the resulting GH process jumps infinitely many times within an infinitesimally small interval, thus enabling the model to capture both discrete and continuous asset price movements (Daal and Madan, 2005). The GH process constructed this way belongs to the time-changed Lévy process family (Barndorff-Nielsen and Shephard, 2001; Carr and Wu, 2004; Huang and Wu, 2004). Taking this time change approach is especially effective for developing more plausible and tractable models for asset pricing.2 In terms of derivative pricing, the methodology lends itself to the derivation of an explicit pricing formula for European options, as the characteristic function of the logarithmic spot price can be conveniently obtained. Furthermore, with a judicious choice of the Lévy process and the time deformation process, the modeler can easily incorporate stochastic volatility (SV) to relax the i.i.d. restriction on the spot return dynamics under the baseline models. Tractability in deriving the option pricing formula is preserved under the time-changed Lévy process framework, provided that the characteristic functions of the relevant processes are known. For example, with the use of a time-dependent time deformation process, the baseline GH option pricing model can be readily extended to allow for mean-reverting SV of spot returns. This extension is detailed in Section 4. As another example, short- and long-range dependence in volatility can be induced by choosing the time deformation process to be a superposition of many time-dependent processes (Finlay and Seneta, 2012), although its empirical performance is yet to be examined in the literature. Generalizing this idea further, it is possible to develop option pricing models with multi-factor volatility dynamics, such as the models in Andersen, Fusari, and Todorov (2015) and Gruber, Tebaldi, and Trojani (2015). The multi-factor structure provides the necessary flexibility in modeling the underlying risk-neutral spot dynamics. Empirical evidence suggests that they are capable of capturing option stylized facts (e.g., implied volatility surface) more accurately than the conventional models. Eberlein and Prause (2002) and Prause (1999) consider a GH model for option pricing. There are several important differences that separate our approach from theirs. First, using time-series data of stock returns, Eberlein and Prause (2002) and Prause (1999) estimate the GH model, which is then used to price the options and obtain the model-implied volatilities. This is made possible by means of Esscher transform, which allows them to recover the risk-neutral spot price dynamics from the statistical measure. We take a more direct route by specifying the spot price dynamics under the risk-neutral measure. The risk-neutral parameters of our flexible GH model are then estimated directly using option data. Our modeling approach is therefore in line with the risk-neutral option pricing literature (e.g., Heston, 1993; Bates, 1996; Madan, Carr, and Chang, 1998; Pan, 2002; Carr and Wu, 2003). Second, Eberlein and Prause (2002) and Prause (1999) express the European option price in terms of the integral of a τ-fold convolution of the GH densities (τ is the time-to-maturity) which cannot be further simplified as the GH distribution is not closed under convolution. To circumvent the numerical challenge, they employ Fourier inversion to obtain the cumulative probabilities from the characteristic function of the log-return process. By contrast, using the Carr and Madan (1999) approach, we obtain an explicit option pricing formula in terms of the characteristic function of the logarithmic spot price process. Third, in terms of numerical optimization, Eberlein and Prause (2002) and Prause (1999) estimate the risk-neural parameters by first keeping the index parameter fixed and then optimizing over the other parameters.3 We are able to overcome the numerical difficulty by estimating all four parameters of the flexible GH model jointly and efficiently. This is the combined result of a prudent choice of model parameterization, appropriate variable transformation, and a larger sample offered by the option panel compared with the time series of spot prices. In addition, we conduct an empirical study of the flexible GH model and its special cases by investigating its performance in modeling S&P 500 Index options. For a fair comparison, the GH model is evaluated against benchmark option pricing models that assume i.i.d. spot returns, including the finite moment log-stable (FMLS) and BS models.4 Among the special cases of GH, we will evaluate the NIG model, and additionally closely examine the VG and t models, which Eberlein and Prause (2002) and Prause (1999) precluded from being estimated under the risk-neutral setting. A direct comparison between the VG and t models is theoretically and practically motivated, since the models are complementary limiting cases of the flexible GH model (see Subsection 2.1), yet the VG model is much more widely used than the t model in empirical option pricing. The main empirical findings of our paper are summarized as follows. First, among all six models under consideration, the flexible GH model yields the best in-sample fitting and out-of-sample forecast across the entire panel of option prices; the added flexibility offered by the extra shape parameter is indispensable for the superior performance. Second, the t model often performs as good as, and sometimes even better than, the flexible GH model in fitting and predicting the prices of certain option categories. For example, the t model often provides an adequate fit to the left tail of the log-return distribution, as revealed by its superior performance in modeling the in-the-money (ITM) put options in the main estimation. Third, as shown in the robustness analysis, the superior overall out-of-sample performance of the flexible GH model is preserved over a turbulent sample period in 2011, although the category-specific results are mixed—the FMLS model surpasses the GH model in predicting the prices of at-the-money (ATM) and ITM options, thus pointing to the necessity of excess kurtosis in the spot price distribution for accurate option pricing. In summary, our contributions are three-fold. First, we study a flexible GH option pricing model in the time-changed Lévy process framework. Taking this dynamic perspective is fruitful in option pricing as it motivates the way the GH model is parameterized. The flexible GH model encompasses six option pricing models as special cases, including the popular VG model, as well as new models which employ certain familiar distributions, such as the t distribution. Second, we obtain an explicit, semi-closed-form option pricing formula under our GH model. The formula is derived using the characteristic function approach (Carr and Madan, 1999). Third, we investigate the empirical performance of the GH model and some of its special cases. By comparison with benchmark models such as the FMLS model (Carr and Wu, 2003) and the BS model, we illustrate the convenience of calibrating the GH model, and its consistent superiority in delivering in-sample fitting and out-of-sample prediction over the panel of S&P 500 Index options. To our knowledge, this paper provides the first systematic empirical analysis of the GH model and some of its special cases, including the VG, t, and NIG models. Although our empirical study focuses on Lévy-based option pricing models that assume i.i.d. spot returns, the empirical results would shed light on future work on evaluating the empirical performance of SV extensions of such models, some of which are developed in Section 4. It would be further elucidating, for example, to compare the GH–SV model with the conventional Heston (1993) model, whose underlying spot dynamics are characterized by a time-changed GH process and time-changed Brownian motion, respectively. The paper is structured as follows. Section 1 begins with a description of a subordinated process and presents a corresponding option pricing framework. In Section 2, we formulate the flexible GH option pricing model and its six special cases, including the two limiting cases. Section 3 comprises our empirical study. The data description, parameter estimates, in-sample fit, orthogonality test, and out-of-sample pricing results are provided therein. In Section 4, we discuss the extension of the GH model to incorporate SV. Section 5 concludes the paper. 1 Time-Changed Lévy Processes and Option Pricing We focus on a generic option pricing framework in which the log-return process of the underlying asset is driven by a time-changed Lévy process [Barndorff-Nielsen and Shephard, 2001; Carr et al., 2003 (hereafter CGMY, 2003); Carr and Wu, 2004]. Let Lt≡L(t) denote a Lévy process, and let Tt≡T(t) be an increasing, positively valued time deformation process. The resulting time-changed Lévy process is given by Zt≡Z(t)=L(T(t)). (1) The time deformation process Tt represents the total elapsed time in a transformed time scale which may reflect the varying rate at which economic and business activities occur. As will be explored in this paper, this modeling framework offers flexibility in the specification of the Lévy process and the time deformation process. It encompasses many option pricing models in the literature. In particular, the time deformation process Tt can be deterministic (e.g., an identity process Tt=t, as in BS model), a pure-jump Lévy process with stationary and independent increments (e.g., a gamma process, as in Madan, Carr, and Chang (1998), or a time-dependent process, possibly driven by the same random process that drives Lt (e.g., CGMY, 2003; Carr and Wu, 2004; Huang and Wu, 2004). To preserve the business clock interpretation, we require that Tt satisfies the following conditions: T0=0. Tt is non-negative and increasing with t. E[T1]=1. The third condition is a normalizing restriction on Tt. If Tt is a Lévy process, then, by infinite divisibility, this condition implies that E[Tt]=t for all positive time t. To derive the option pricing formula, it is instrumental that we obtain the characteristic function of the time-changed process Zt. Defining the characteristic functions of Lt and Tt by φLt(u)=E(eiuLt) and φTt(u)=E(eiuTt), respectively, we can compute the characteristic function of Zt as follows φZt(u)=φTt(−iφL1(u)). (2) We assume that Zt drives the log-return process of the underlying asset. More specifically, under the risk-neutral probability space5 (Ω,Q,{Ft}) equipped with the filtration {Ft}, the spot price of the underlying asset evolves according to the following dynamics: for all t  ≥  0, St=S0 exp{(r−q+ω)t+Zt}. (3) For simplicity, the risk-free rate r and the dividend yield q are assumed to be constant in our analysis. The drift adjustment ω is required to ensure that the discounted spot price process is a martingale under Q. It is given by ω=−1t log φZt(−i)=−log φZ1(−i) (see Appendix A). As a result, the characteristic function of the log-price is given by φ log St(u)=S0iu exp{iu(r−q+ω)t}φZt(u). (4) Using the characteristic function approach to option pricing (e.g., Scott, 1997; Bakshi and Madan, 2000; Carr and Wu, 2004), the price of an European option written on this underlying asset with the spot price process, defined in Equation (3), can be obtained. The time-t price of an European call with strike price K and time-to-maturity τ = T– t is Ct=EQ[e−rτ⁡max(ST−K,0)|Ft]=Ste−qτΠ1−Ke−rτΠ2, (5) where EQ[·] denotes the expectation taken under Q, and the quantities Π1 and Π2 are the prices of some Arrow–Debreu securities (Bakshi and Madan, 2000), given by Π1=12+1π∫0∞Re(e−iω log(K)φ log St(ω−i)iωφ log St(−i))dω,Π2=12+1π∫0∞Re(e−iω log(K)φ log St(ω)iω)dω. By the put-call parity, the price of the corresponding European put option with the same strike price and maturity date is Pt=Ke−rτ(1−Π2)−Ste−qτ(1−Π1). (6) 2 The Baseline Model In this section, we introduce the baseline GH option pricing model. Recall from Equation (3) that the spot price is driven by the time-changed Lévy process Zt. Under the flexible GH model specification, the time-changed Lévy process Zt is constructed by subordinating a Brownian motion with drift, B(t), to a GIG process g(t). More precisely, the time-changed Lévy process is defined by Zt:=Z(t)=B(g(t))=θg(t)+σW(g(t)), (7) where θ is the drift parameter, σ  >  0 is the volatility parameter, and W(t) is the standard Wiener process. The GIG process g(t) is a positively valued Lévy process such that g1:=g(1) follows the GIG distribution with parameters p, γ, and δ. The probability density function of the GIG distribution is provided in Appendix B.1. In particular, by defining the parameter ζ:=γδ, the mean is E[g1]=ζKp+1(ζ)γ2Kp(ζ), where Kj(·) is the modified Bessel function of the third kind with index j (Jørgensen, 1982). The mean of g1 is well defined provided that the parameters p and ζ satisfy the conditions: (i) ζ>0 if −1≤p≤0,(ii) ζ≥0  otherwise. By imposing the normalization E[g1]=1, we obtain the parameter restriction γ2=ζKp+1(ζ)Kp(ζ), (8) and hence the GIG distribution is characterized by two parameters: the index parameter p and the shape parameter ζ. From a static point of view, the process at time 1, Z1, follows the GH distribution, which can be viewed as a normal mean–variance mixture distribution with a GIG mixing density. For any fixed t, the distribution of Zt is skewed and leptokurtic, due to the mixing in the mean through the Brownian motion’s drift θ. From a dynamic point of view, the time-changed Lévy process Zt is constructed by subordinating a drifted Brownian motion to a GIG process. The distribution of g(t) is obtained as the t-fold convolution of the GIG distribution of g1. In particular, for any t  ≥  0, the increment g(t  +  1)–g(t) follows the same distribution as that of g1. Note that g(t) is not distributed as GIG for t  ≠  1, as the GIG distribution is not preserved under convolution. It follows from this dynamic construction that the marginal distribution of Zt is not GH in general except at t  =  1. Nevertheless, by infinite divisibility of the GIG distribution (Barndorff-Nielsen and Halgreen, 1977), the GIG process g(t) is infinitely divisible. As a result, the time-changed process Zt is a Lévy process, and so it is infinitely divisible and has stationary and independent increments. For a Lévy process, the characteristic functions of Zt and Z1 are related by φZt(u)=[φZ1(u)]t (e.g., Barndorff-Nielsen, Mikosch, and Resnick, 2001). The central moments of Zt grow over time just like any generic Lévy process. The mean, variance, and the third central moment of Zt increase linearly with t. In terms of standardized central moments, the skewnesses of Zt and Z1 are related by skew(Zt)=(1/t)skew(Z1) for t  >  0, and the kurtosis of Zt is a weighted average of the kurtosis of Z1 and that of a normal distribution, that is, kur(Zt)=(1/t)kur(Z1)+3(1−1/t) for t  >  0. Given that Z1 follows a GH distribution, we have skew(Z1)≠0 if the drift parameter θ is non-zero, and kur(Z1)  >  3. The same moment properties are therefore valid for Zt for all t > 0. A study of the moment properties is detailed in Appendix B. To price a European option under the flexible GH model, it is necessary to get the characteristic function of Zt. Using the characteristic function of the GH distribution (see Appendix B), we obtain φZt(u)=[1−2γ2(iuθ−12σ2u2)]−pt2[Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ)]t. Provided that θ<γ2−σ22, the drift adjustment is well-defined and is given by ω=−log φZ1(−i)=−p2log[1−2γ2(θ+σ22)]+log[Kp(ζ1−2γ2(θ+σ22))Kp(ζ)]. The European option under the flexible GH model can then be priced using Equations (5) and (6). 2.1 Special Cases of the Flexible GH Model In the family of risk-neutral spot price models with independent and stationary increments, the baseline GH option pricing model is quite flexible in that it encompasses a number of particular cases. In this section, we study six of its special cases, all of which are indexed by three parameters. These six models can be classified into two groups: four special cases obtained by restricting the value of the index parameter p, while the other two obtained as limiting cases as the shape parameter ζ:=δγ→0. Table 1 summarizes the parameter restrictions required to obtain the six special cases of the flexible GH model. Table 1 Special cases of the GH option pricing model GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) Notes: ζ=γδ. For the H, RH, NIG and NRIG models, γ is given by (8). For all models, E [g1]=1. Table 1 Special cases of the GH option pricing model GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) GH (p,ζ,θ,σ) parameter conditions Density of X1 Density of g1 Hyperbolic model Positive hyperbolic GH (1,ζ≥0,θ,σ), δ≥0, γ>0 H (ζ,θ,σ) PH (δ,γ) Reciprocal hyperbolic model Reciprocal positive hyperbolic GH (−1,ζ>0,θ,σ), δ>0, γ>0 RH (ζ,θ,σ) RPH (δ,γ) Normal inverse Gaussian model Inverse Gaussian GH (−12,ζ>0,θ,σ), δ>0, γ>0 NIG (ζ,θ,σ) IG (δ,γ) Normal reciprocal inverse Gaussian model Reciprocal inverse Gaussian GH (12,ζ≥0,θ,σ), δ≥0, γ>0 NRIG (ζ,θ,σ) RIG (δ,γ) Variance gamma model Gamma GH (p>0,0,θ,σ), δ = 0, γ=2p VG (p,θ,σ) Γ(p,p) t model Reciprocal gamma GH (p<−1,0,θ,σ), δ=2(−p−1), γ = 0 t−2p(−p,θ,σ) R Γ(−p,−p−1) Notes: ζ=γδ. For the H, RH, NIG and NRIG models, γ is given by (8). For all models, E [g1]=1. The first group of three-parameter sub-models is obtained as follows. By setting p  =  1, the flexible GH model reduces to a H model, in which the mixing random variable g1 follows the positive H mixture distribution. Setting p  =  −1 leads to a RH model, in which the mixing distribution is reciprocal positive H. For p=−(1/2), we obtain the NIG distribution, generated by a normal mixture of the inverse Gaussian distribution. For p=1/2, we obtain the NRIG distribution, generated by a normal mixture of the reciprocal inverse Gaussian distribution. The second group of sub-models is obtained by restricting the shape parameter ζ to 0 and the sign of the index parameter p. By letting δ→0 and restricting p  >  0, we obtain the VG model (Madan, Carr, and Chang, 1998), in which g1 follows the gamma distribution with parameters p and γ2/2. By letting γ → 0 and restricting p < 0, we obtain the t model, in which g1 follows the reciprocal (or inverse) gamma distribution with parameters-p and δ2/2. With the normalization condition E[g1]=1 in place, we deduce that γ=2p in the VG case and δ=2(−p−1) in the t case. Figure 1 demonstrates the flexibility of the GH model in the modeling of option prices relative to its six special cases. The plotted surface represents the prices of an ITM call option under the GH model for varying p and ζ (the drift and volatility parameters are fixed). The option prices associated with the first group of sub-models (H, NRIG, NIG, and RH models) are obtained on the four dotted cross-sections of the surface at the respective fixed values of p, while those associated with the second group of sub-models (VG and t models) are found, respectively, along the positive and negative sections of the p axis at the boundary ζ = 0. When both the index and shape parameters are close to zero, the GH model implies a high kurtosis in the risk-neutral distribution of spot prices. This translates into more extreme option prices (lower prices under NIG, RH, and t with p close to 0, and higher prices under NRIG, H, and VG with p close to 0). The GH model offers the most versatility as manifested by the wide range of option prices it generates over all possible combinations of the index and shape parameters. Figure 1 View largeDownload slide Simulated option prices under the GH model. Notes: The price of an ITM call option with a strike price equal to 95% of the spot price, time-to-maturity equal to 3 weeks, θ=−0.04 and σ=0.17. The six special cases are RH at p = –1, NIG at p=−12, NRIG at p=12, hyperbolic (H) at p = 1, variance gamma at ζ = 0 and p > 0, and the t model at ζ = 0 and p<−1. At ζ = 0 and −1≤p≤0, the mean of the GH distribution is undefined. Figure 1 View largeDownload slide Simulated option prices under the GH model. Notes: The price of an ITM call option with a strike price equal to 95% of the spot price, time-to-maturity equal to 3 weeks, θ=−0.04 and σ=0.17. The six special cases are RH at p = –1, NIG at p=−12, NRIG at p=12, hyperbolic (H) at p = 1, variance gamma at ζ = 0 and p > 0, and the t model at ζ = 0 and p<−1. At ζ = 0 and −1≤p≤0, the mean of the GH distribution is undefined. 2.2 The VG Model We note that the characteristic function of the logarithmic spot price under the VG model can be simplified to φ log St(u)=S0iu exp{iu(r−q+ω)t}[1−ν(iuθ−12σ2u2)]−tν, (9) where ω=1νlog[1−ν(θ+12σ2)] for θ<(1ν−σ22). This formula can be obtained as the limit of the corresponding characteristic function under the GH model (see Appendix C). This is equivalent to the VG option pricing model of Madan, Carr, and Chang (1998), which contains three parameters: the drift θ, the volatility σ, and the mixing distribution variance ν (which is equal to 1/p using our index parameter). The distributional properties of the VG model are examined in Appendix B.2. 2.3 The t Model Let us turn to the t option pricing model.6 Under this model, the log-return of the underlying asset follows a skewed t distribution with –2p degrees of freedom. The characteristic function of the logarithmic spot price under the VG model is given by φ log St(u)=S0iu exp{iu(r−q+ω)t}×[2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν(−4(1ν−1)(iuθ−12σ2u2))]t, (10) where ω=−log[2(1ν−1)12νΓ(1ν)[−(θ+12σ2)]12νK1ν(−4(1ν−1)(θ+12σ2))]for θ<−σ22, and Kj(·) is the modified Bessel function of the third kind with index j. This formula can be obtained as the limit of the corresponding characteristic function under the flexible GH model (see Appendix C). The spot price dynamics under the t model is equivalent to the GH-skewed Student’s t model in Aas and Haff (2006) after re-parameterization, except that the latter specifies the spot price process under the statistical measure and is not used for option pricing. The distributional properties of the t model are examined in Appendix B.3. 3 Empirical Study 3.1 The Data The data used for our main estimation are the daily prices of S&P 500 Index European options. Put options are featured in our sample because, compared with their call counterparts, their prices carry more information about the left tail of the price distribution, corresponding to the side to which the price distribution is found to be skewed empirically (Madan and Milne, 1991). The sample period is chosen to be 2012 with 250 trading days. It contains some moderate upward and downward trends in the S&P 500 Index. The volatility index (VIX) fluctuates mildly between 13% and 27% over the period. Both series are displayed in the bottom two panels in Figure 2. In the robustness analysis in Section 3.8, we will choose the most turbulent period in 2011 as our alternative sample and repeat the analysis. Figure 2 View largeDownload slide Time-series plots of daily risk-neutral volatility (σ) and skewness (θ) parameter estimates under the GH option pricing model, and the plots of the S&P 500 Index and VIX. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The plots include the volatility parameter (top panel), the skewness parameter (second panel), the S&P 500 Index (third panel), and VIX (bottom panel), the VIX that measures the implied volatility of S&P 500 Index options. Figure 2 View largeDownload slide Time-series plots of daily risk-neutral volatility (σ) and skewness (θ) parameter estimates under the GH option pricing model, and the plots of the S&P 500 Index and VIX. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The plots include the volatility parameter (top panel), the skewness parameter (second panel), the S&P 500 Index (third panel), and VIX (bottom panel), the VIX that measures the implied volatility of S&P 500 Index options. The cross-sectional data structure for the put option sample used in the main estimation is summarized in Table 2. To retain a majority of put options that are more liquidly traded in the market, we focus on the puts whose moneyness (the strike-to-spot price ratio) lies between 0.94 and 1.06, and whose time-to-maturity is longer than a week. On each day, the sample retains only the three batches of put options belonging to the short-term, medium-term, and long-term categories, corresponding to the nearest-, second-nearest-, and third-nearest-to-maturity put options, respectively. This leaves us with 60 put options on average per day. In total, the option panel consists of 15,058 put prices. To price the European puts, we use the annualized yield of the 1-month U.S. Treasury bill as the risk-free rate. The time series of S&P 500 Index is adjusted by the mean annualized dividend yield of the stocks. Table 2 Options data characteristics of the main estimation sample Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Notes: A contingency table for the sample of S&P 500 Index European put option prices in 2012. The sample is selected according to the criteria in the Data section. n = 15,058. Average option prices are displayed in square brackets. K is the strike price, S is the spot price, and τ is the time-to-maturity. Table 2 Options data characteristics of the main estimation sample Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Short-term Medium-term Long-term All 1 week <τ< 1 month 1 month ≤τ≤ 3 months τ> 3 months Out-of-the-money 8% 24% 11% 43% 0.94<KS<0.98 [$3.49] [$13.40] [$26.86] [$14.87] At-the-money 8% 22% 11% 41% 0.98≤KS≤1.02 [$16.53] [$29.20] [$43.80] [$30.55] In-the-money 4% 9% 3% 16% 1.02<KS<1.06 [$51.81] [$58.79] [$71.32] [$59.01] All 20% 55% 25% 100% [$18.29] [$27.48] [$39.17] [$28.42] Notes: A contingency table for the sample of S&P 500 Index European put option prices in 2012. The sample is selected according to the criteria in the Data section. n = 15,058. Average option prices are displayed in square brackets. K is the strike price, S is the spot price, and τ is the time-to-maturity. Due to the way the sample is selected, the time-to-maturities exhibit some artificial cyclicity over time. This can be seen from the time-series plot of the average time-to-maturity at the bottom panel of Figure 3. Note that the S&P 500 Index options expire on the third Friday of each month. As the maturity date approaches, the time-to-maturities of the options in the sample decreases gradually over time until some of them (i.e., the time-to-maturities of those puts in the short-term category) falls below 1 week, at which point these nearest-to-maturity options will be excluded from the sample and replaced by a new batch of options that forms the long-term category. The original long-term (medium-term) category now becomes the medium-term (short-term) category. This roll-over of the sample, which occurs around the second Friday of each month, leads to an abrupt increase in the time-to-maturities by approximately 1 month across the three term categories. It is important to bear this in mind when we interpret the time-series variation of the parameter estimates. Figure 3 View largeDownload slide Time series plots of daily risk-neutral shape parameter estimates from different models, and the plot of daily averages of time-to-maturity in the main estimation sample in 2012. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The shape parameters include the tail index parameter α under the FMLS model (top panel), and the index parameter p under the GH model and its special cases (second panel), and the shape parameter ζ under the GH and NIG models (third panel), with ζ→0 under the VG and t models. The average time-to-maturity (τ, in days) is computed by averaging the time-to-maturity of the cross-section of puts each day (bottom panel). Figure 3 View largeDownload slide Time series plots of daily risk-neutral shape parameter estimates from different models, and the plot of daily averages of time-to-maturity in the main estimation sample in 2012. Notes: The sample data are the S&P 500 Index European put prices in 2012, selected according to the criteria detailed in the Data section. The shape parameters include the tail index parameter α under the FMLS model (top panel), and the index parameter p under the GH model and its special cases (second panel), and the shape parameter ζ under the GH and NIG models (third panel), with ζ→0 under the VG and t models. The average time-to-maturity (τ, in days) is computed by averaging the time-to-maturity of the cross-section of puts each day (bottom panel). 3.2 Model Calibration and Parameter Estimation The option pricing models are calibrated daily by minimizing the sum of squared percentage pricing errors SSPEt(Θ) on each date t: SSPEt(Θ)=∑i=1nt(O^it(Θ)−OitOit)2, where Θ is the vector of model parameters, nt is the cross-sectional sample size on date t, and O^it(Θ) and Oit are the predicted and observed option prices (in dollars), respectively. The optimization search for Θ is carried out over the parameter space using the Nelder–Mead simplex algorithm (Gilli and Schumann, 2012). It is implemented through the MATLAB function fminsearch. To mitigate the sensitivity of the solution on the initial value of the parameter vector, a preliminary search is conducted as follows: starting with many different candidate vectors for Θ[0], we run a first-stage optimization with a fixed number of iterations (e.g., 10), then we select the optimal candidate vector that yields the smallest SSPEt(Θ). Starting at this optimal candidate vector, the numerical algorithm is then run until the local minimum is found. 3.3 In-Sample Estimation In this section, we present the main estimation results and examine the in-sample fit of various models. Table 3 summarizes the risk-neutral parameter estimates under the BS model, the FMLS model, as well as the GH model and three of its special cases (NIG, VG, and t). Table 3 Summary statistics of risk-neutral parameter estimates (in-sample) Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Notes: Parameter estimates and root mean-squared percentage error (RMSPE) are obtained daily for the year 2012, which includes 15,058 put prices over 250 days. The average daily sample size is 60. Table 3 Summary statistics of risk-neutral parameter estimates (in-sample) Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Parameter Mean Standard deviation Minimum Maximum BS  σ 0.146 0.017 0.116 0.206  RMSPE 16.69% 4.49% 8.85% 34.65% FMLS model  σ 0.099 0.013 0.076 0.142  α 1.927 0.040 1.848 1.999  RMSPE 12.00% 5.18% 4.81% 30.63% NIG model  σ 0.155 0.022 0.121 0.231  θ −0.044 0.019 −0.108 −0.007  ζ 4.910 5.791 0.006 33.744  RMSPE 9.10% 4.45% 3.78% 24.24% VG model  σ 0.150 0.022 0.114 0.224  θ −0.044 0.016 −0.100 −0.009  p 9.305 9.752 1.777 62.350  RMSPE 9.50% 3.81% 4.42% 26.52% t model  σ 0.158 0.025 0.113 0.248  θ −0.047 0.021 −0.149 −0.013  p −5.550 5.569 −42.824 −1.448  RMSPE 8.38% 3.28% 3.51% 23.98% GH model  σ 0.158 0.025 0.113 0.225  θ −0.046 0.020 −0.098 −0.012  p −5.151 5.465 −30.180 8.014  ζ 0.417 1.483 0.000 18.778  RMSPE 8.38% 3.27% 3.51% 23.98% Notes: Parameter estimates and root mean-squared percentage error (RMSPE) are obtained daily for the year 2012, which includes 15,058 put prices over 250 days. The average daily sample size is 60. Starting with the volatility parameter σ, which appears in all models being considered, the FMLS model has the smallest average equal to 9.9%, followed by the BS model (14.6%), the VG model (15.0%), the NIG model (15.5%), and then the t and GH models (both 15.8%). The t model yields the single largest estimate equal to 24.8%. The skewness parameter θ is contained in all models except for BS and FMLS.7 It is estimated to be negative on average with a similar value across models: −0.044 for the NIG and VG models, −0.047 for the t model, and −0.046 for the GH model. However, the most extreme estimate is obtained under the t model: −0.149, compared with −0.108 for the NIG model, −0.100 for the VG model, and −0.098 for the GH model. The time-series plots of the daily estimates of σ and θ under the GH model are displayed in the top two panels of Figure 2. They reveal a negative association between the volatility and skewness parameters of the log-return distribution, confirming the presence of leverage effect (Christie, 1982): the volatility tends to be higher when the log-returns are more negatively skewed. Let us turn to the shape parameters (α for FMLS, and p and ζ for GH and its special cases). Under the FMLS model, the estimates of the tail index α lie in a narrow range between 1.848 and 1.999, with an average estimate equal to 1.927. By the properties of the asymmetric stable distribution (Property 1.5 of Carr and Wu, 2003), the log-return distribution has infinite moments of order 2 and higher if α  <  2. Under the GH model and its special cases, the shape is controlled by p (the index parameter) and ζ. Recall that p is constrained to be positive under the VG model, negative under the t model, and equal to −0.5 under the NIG model, and that ζ → 0 under the VG and t models, but ζ is free under the NIG model. The smaller the magnitude of p and ζ, the fatter the tail of the log-return distribution. From the empirical estimation, the VG model yields an average estimate of p equal to 9.305, which is higher in magnitude than that under the t model (−5.550, which translates into 11.1 degrees of freedom). The p estimates under the GH model average to −5.151 and range over −30.180 to 8.014, with a minimum magnitude of 1.117. These estimates are more in line with the t model than with the VG model. The estimation of p under the flexible GH model indicates that the special cases of GH, including RH (p  =  –1), NIG ( p=−(1/2)), NRIG ( p=1/2), and H (p = 1), are not supported by the data. In terms of the shape parameter ζ, the average estimate is 4.910 under the NIG model, but the average is only 0.417 under the flexible GH model. The discrepancy indicates that, due to the restriction on p, ζ needs to be inflated substantially under the NIG model in order to capture the excess kurtosis of the risk-neutral distribution. The time series of daily estimates of p under the various GH sub-models is plotted in Figure 3. It reveals a cyclical behavior in the degree of tail-thickness in the log-return distribution that aligns with the option expiration schedule.8 The tail of log-return distribution is thick at the start of the trading month (low magnitude of p), and gradually becomes thinner (increasing magnitude of p) toward the second Friday of each month as there is diminishing information uncertainty. Then, when the option data set rolls over to exclude options with less than 1 week to expiry and include a new batch of options, the index parameter p falls sharply in magnitude and the log-return distribution becomes fat-tailed again. This pattern repeats itself every month. 3.4 In-Sample Pricing Performance We turn now to the first of the performance measures, in-sample pricing error, also reported in Table 3. As expected, the daily-averaged root mean-squared percentage error (RMSPE) for the BS model is the highest at 16.69%, since it is the most restrictive model with only one parameter controlling the dispersion of the log-return distribution. The FMLS model reduces the RMSPE to 12.00%, the reduction in error representing an improvement to in-sample fit offered by the kurtosis parameter and a fixed negative skewness parameter. The VG, NIG, and t models further reduce the RMSPE to 9.50%, 9.10%, and 8.38%, respectively. The further improved fit is due to the flexibility in the skewness parameter and the evidently more correct specification of the log-return distribution over the FMLS model. The superior performance of the t model over the VG model indicates that the t distribution is a more plausible model empirically for the risk-neutral log-return distribution.9 The superiority of the t model over other sub-models is corroborated by the GH model’s parameter estimates: the index parameter p estimates are predominantly negative, and the estimates of the shape parameter ζ are close to zero. The in-sample fit of the t model, according to the average RMSPE, is as good as that of the GH model (8.38%). From the main estimation results, the GH model family clearly outperforms the BS and FMLS models. Allowing p to be free is crucial to the observed improvement in the in-sample fit, although the benefit of having an additional shape parameter ζ is not fully exploited in this particular dataset. We however wish to emphasize that, while the results are valid for the baseline dataset with sample period spanning 2012, the conclusions may depend on the sample period used. For comparison, the results of the robustness analysis based on a different sample period are presented in Section 3.8. 3.5 In-Sample Performance by Option Categories While Table 3 offers insight to overall in-sample fit, in Table 4 we examine how the models perform in terms of in-sample fit for different types of put options in the cross section, classified according to different moneyness and time-to-maturity groups as defined in Table 2. Table 4 In-sample pricing RMSPEs and mean absolute percentage errors (MAPEs) In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 Notes: The smallest error measure within each group is italicized and underlined. The bold values show the overall error measure. Parameters are estimated using all options on a given day regardless of their time-to-maturity and moneyness (an average of 60 observations per day, and 250 days in year 2012). Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 250 testing days collectively. The sample size is n = 15,058. Table 4 In-sample pricing RMSPEs and mean absolute percentage errors (MAPEs) In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 In-sample RMSPE (%) In-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 33.03 15.26 12.73 19.57 25.78 12.02 10.12 14.23  FMLS model 26.94 12.51 14.31 16.73 20.09 9.80 11.89 12.34  NIG model 16.33 11.14 12.13 12.57 10.35 8.45 9.47 9.08  VG model 15.72 11.31 12.71 12.64 10.98 8.55 10.01 9.39  t model 10.45 10.80 12.19 11.10 8.06 8.04 9.46 8.40  GH model 10.47 10.78 12.19 11.10 8.07 8.04 9.46 8.40 At-the-money  BS 25.31 16.34 9.55 17.15 23.37 15.29 8.38 15.06  FMLS model 19.13 8.53 6.24 11.02 15.77 7.31 5.13 8.40  NIG model 12.69 6.55 5.76 7.97 9.26 5.30 4.70 5.92  VG model 12.83 6.66 5.80 8.07 9.99 5.42 4.75 6.14  t model 9.63 6.03 5.56 6.79 7.93 5.03 4.51 5.46  GH model 9.61 6.03 5.55 6.78 7.91 5.03 4.50 5.46 In-the-money  BS 7.78 12.29 12.47 11.34 6.57 11.77 12.20 10.51  FMLS model 6.85 9.03 9.31 8.57 5.90 8.66 8.96 8.00  NIG model 6.25 9.09 9.91 8.60 5.46 8.72 9.63 8.03  VG model 6.66 9.37 10.09 8.88 5.83 9.07 9.77 8.35  t model 6.03 8.69 9.49 8.24 5.30 8.36 9.20 7.71  GH model 6.13 8.70 9.54 8.27 5.35 8.37 9.25 7.73 All  BS 26.60 15.26 11.38 17.48 20.99 13.30 9.56 13.98  FMLS model 21.17 10.49 10.86 13.47 15.55 8.60 8.55 10.02  NIG model 13.40 9.18 9.54 10.27 8.94 7.22 7.35 7.61  VG model 13.19 9.35 9.89 10.38 9.56 7.37 7.63 7.89  t model 9.39 8.79 9.47 9.08 7.46 6.88 7.21 7.08  GH model 9.40 8.78 9.47 9.08 7.46 6.88 7.21 7.08 Notes: The smallest error measure within each group is italicized and underlined. The bold values show the overall error measure. Parameters are estimated using all options on a given day regardless of their time-to-maturity and moneyness (an average of 60 observations per day, and 250 days in year 2012). Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 250 testing days collectively. The sample size is n = 15,058. Let us start by comparing the performance across the moneyness categories. For ATM options, the GH model achieves the best in-sample pricing performance compared with other models: its average RMSPE is 6.78%, compared with 17.15% for the BS model, 11.02% for the FMLS model, 8.07% for the VG model, 7.97% for the NIG model, and 6.79% for the t model. For ITM options, the t model offers the best in-sample fit, with an average RMSPE of 8.24%, outperforming the GH model (8.27%). As all options in our sample are European puts, this suggests that the t distribution provides an optimal fit of the left tail of the log-return distribution. As for out-of-the-money (OTM) options, the GH model ties with the t model, and both of the models clearly outperform the other competitors. Similarly, a comparison across the different time-to-maturity categories reveals that the GH and t models give the best in-sample fit among all models, with the t model comes slightly ahead in pricing the short-dated puts, and the GH model becomes slightly better in pricing the medium-dated puts. Given that the FMLS model is outside the GH model family, one may be curious whether its peculiar features (e.g., infinite kurtosis of the log-return distribution) would facilitate more accurate pricing of certain options. Table 4 reveals that the FMLS model offers the best fit to the long-dated, ITM puts, achieving the smallest average RMSPE of 9.31% among all other models. As the price of a long-dated, ITM put is sensitive to the extreme left tail of the distribution of log-returns from now until maturity, the empirical result points to the necessity of excess kurtosis to model the extreme negative returns over a long time horizon. In this sense, the negatively skewed α-stable distribution under the FMLS model provides a better fit than the GH model, in which the kurtosis of log-returns tends to that of a normal distribution as the time horizon increases (see Appendix B.1). Nevertheless, the FMLS model is not adequate in modeling other sections of the long-range log-return distribution. Its RMSPE for long-dated OTM puts (14.31%) is worse than all other models, including the BS model (12.73%). Even though the log-return distribution under the FMLS model has the fattest tails, the empirical fit in the parts other than the left tail is not the best compared with other models with finite kurtosis. For robustness analysis of the in-sample performance measure, the results for the mean absolute percentage pricing error (MAPE) are also computed. While the average MAPE measures are generally smaller than the average RMSPE measures, probably due to the effect of extreme pricing errors on RMSPE, changing the measure to MAPE retains the optimal model in any given option category. 3.6 Orthogonality Test A comprehensive cross-sectional analysis of the pricing errors leads us to the second yardstick for assessing model adequacy. A well-specified model should not only minimize some loss function on pricing errors (e.g., RMSPE and MAPE) but also yield individual pricing errors that are independent of moneyness, time-to-maturity, and interest rates (Rubinstein, 1985). Furthermore, the pricing error measure used for assessing model adequacy needs to accord with that used for parameter estimation (Christoffersen and Jacobs, 2004). To ensure consistency with the model estimation procedure in Section 3.2, we therefore evaluate the model fitting using the sequence of option pricing errors relative to the observed option prices: εit=Oit(Θ^t)−OitOit, where Oit is the observed option price in dollars, and Oit(Θ^t) is the option price estimated by the model. We investigate the orthogonality of the relative pricing errors by regressing εit on some option and market attributes, including moneyness (defined as the ratio of the strike to the spot price), squared moneyness, time-to-maturity, squared time-to-maturity, and the risk-free interest rate. More specifically, we consider the following regression model for the orthogonality test: εit=b0+b1(KiSt)+b2(KiSt)2+b3τit+b4τit2+b5rt+eit, where Ki and τit are the strike price and time-to-maturity of the ith option at time t, St is the S&P 500 Index level at time t, rt is the risk-free interest rate at time t, and eit is the random error term. If the model is well specified, we expect close-to-zero regression coefficients and a low regression R2, both signifying greater orthogonality of the relative pricing errors. We run the regression on the entire 250-day sample consisting of 15,058 observations of relative pricing error. The results are presented in Table 5. As revealed by the t tests on the regression coefficients and the F tests on the regression R2, all models are generally misspecified to various degrees. Comparing the extent of error orthogonality across models, the BS model yields the highest regression R2 (72.0%), followed by the FMLS model (46.5%), the GH model (34.6%), the t model (34.5%), the VG model (32.8%), and the NIG model (28.7%). While NIG model achieves the least predictable relative pricing errors, its total sum of squared (TSS) pricing errors is larger in magnitude compared with the t and GH models, which have the smallest TSS among all models being considered. There is a trade-off to be made—from a practitioner’s perspective the overall pricing accuracy is likely to take precedence over the orthogonality of pricing errors when choosing among the models. Table 5 Results of orthogonality test Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Notes: S&P 500 Index put options from 2012 are used, with sample size n = 15,058. Moneyness is the strike-to-spot price ratio, and time-to-maturity is in years. Heteroskedastic standard errors (White, 1980) are shown in parentheses, * indicates statistical significance at the 1% level of significance. The critical t-statistics are, respectively, t(0.005,15052)=±2.58 and t(0.025,15052)=±1.96. TSS is the total sum of squared errors. At the 1% level of significance, the critical F(0.01,5,15052)-statistic is 3.02. Table 5 Results of orthogonality test Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Explanatory variable BS FMLS NIG VG t GH Intercept −103.14 −54.58 −13.11 −18.96 −10.78 −10.92 (0.943)* (0.943)* (0.856)* (0.827)* (0.742)* (0.742)* Moneyness 204.36 107.88 24.95 36.72 20.37 20.63 (1.887)* (1.880)* (1.709)* (1.654)* (1.483)* (1.483)* Moneyness2 −101.07 −53.16 −11.75 −17.65 −9.48 −9.61 (0.944)* (0.937)* (0.853)* (0.827)* (0.741)* (0.740)* Time-to-maturity 0.322 −0.329 −0.550 −0.609 −0.632 −0.630 (0.047)* (0.049)* (0.038)* (0.036)* (0.027)* (0.027)* Time-to-maturity2 −1.204 −0.020 0.893 0.940 1.070 1.062 (0.103)* (0.107) (0.085)* (0.081)* (0.064)* (0.064)* Interest rate 5.851 6.953 5.493 5.653 5.354 4.802 (2.249)* (2.415)* (2.157) (2.136)* (1.831)* (1.826)* TSS 438.39 268.73 157.32 160.16 122.63 122.62 R2 72.0% 46.5% 28.7% 32.8% 34.5% 34.6% F-statistic (7735.1)* (2615.1)* (1211.9)* (1466.7)* (1586.2)* (1595.3)* Notes: S&P 500 Index put options from 2012 are used, with sample size n = 15,058. Moneyness is the strike-to-spot price ratio, and time-to-maturity is in years. Heteroskedastic standard errors (White, 1980) are shown in parentheses, * indicates statistical significance at the 1% level of significance. The critical t-statistics are, respectively, t(0.005,15052)=±2.58 and t(0.025,15052)=±1.96. TSS is the total sum of squared errors. At the 1% level of significance, the critical F(0.01,5,15052)-statistic is 3.02. 3.7 Out-of-Sample Pricing Performance Finally, we examine the out-of-sample pricing performance of the option pricing models. This is motivated from both a practical (price prediction) as well as a statistical consideration. While in-sample pricing performance benefits from having additional parameters, this is usually not the case when it comes to out-of-sample prediction, as prediction accuracy can be penalized by overfitting a model (Bakshi, Cao, and Chen, 1997). In order to compute the out-of-sample pricing errors, we first estimate the model over a rolling 5-day training period. Using the calibrated model, we then compute the 1-day ahead forecast of the option prices. With the 2012 S&P 500 Index put options data, there are 245 testing samples and a total of 14,794 pricing errors. Each training sample has on average 301 data points. The out-of-sample performance results are presented in Table 6. Table 6 Out-of-sample pricing RMSPEs and MAPEs in the main estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample spans 5 days (on average containing 301 data points) and each testing sample is 1 day ahead of the training sample (on average containing 60 data points). There are 245 testing days. Parameters are estimated using all options in the entire cross-section of the sample. Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 245 testing days collectively. The sample size is n  =  14,794. Table 6 Out-of-sample pricing RMSPEs and MAPEs in the main estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 36.35 19.36 14.34 23.35 26.95 15.38 11.58 17.04  FMLS model 43.01 19.20 15.81 24.93 33.59 15.53 13.03 18.35  NIG model 37.13 17.60 15.85 22.35 26.14 14.20 13.05 16.19  VG model 33.77 17.48 15.74 21.87 22.12 13.99 13.13 15.60  t model 33.63 17.27 15.35 21.66 20.68 13.86 12.66 15.09  GH model 29.31 17.28 15.32 20.21 19.52 13.89 12.65 14.85 At-the-money  BS 27.98 17.61 10.06 18.61 23.11 15.38 8.49 15.06  FMLS model 21.66 12.81 8.49 14.10 16.04 10.46 6.90 10.60  NIG model 16.29 10.64 8.41 11.46 12.05 8.50 6.53 8.67  VG model 20.44 9.76 7.28 12.10 14.52 7.87 5.73 8.59  t model 16.51 10.23 7.64 11.17 11.77 8.27 6.11 8.37  GH model 15.18 10.21 7.50 10.77 11.56 8.26 6.03 8.31 In-the-money  BS 7.14 12.10 11.59 10.98 5.92 11.26 10.93 9.87  FMLS model 6.36 9.64 9.63 8.94 5.49 8.79 8.77 7.96  NIG model 5.93 9.39 9.43 8.67 5.07 8.57 8.66 7.71  VG model 6.45 9.35 9.00 8.66 5.50 8.61 8.29 7.78  t model 5.64 9.08 9.29 8.39 4.79 8.27 8.58 7.45  GH model 5.70 9.13 9.28 8.43 4.80 8.32 8.57 7.48 All  BS 29.88 17.60 12.30 19.94 21.70 14.68 10.13 15.10  FMLS model 30.71 15.47 12.39 18.98 21.05 12.34 9.83 13.48  NIG model 25.90 13.89 12.38 16.71 16.37 10.94 9.67 11.72  VG model 26.01 13.55 11.96 16.73 16.31 10.61 9.30 11.52  t model 24.86 13.54 11.85 16.33 14.53 10.65 9.29 11.16  GH model 21.89 13.54 11.79 15.38 13.93 10.67 9.25 11.04 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample spans 5 days (on average containing 301 data points) and each testing sample is 1 day ahead of the training sample (on average containing 60 data points). There are 245 testing days. Parameters are estimated using all options in the entire cross-section of the sample. Pricing errors are classified by time-to-maturity and moneyness in accordance with the categories in Table 2. The errors are then averaged across the 245 testing days collectively. The sample size is n  =  14,794. First, we focus on the out-of-sample performance of models in predicting the entire cross-section. The GH model emerges as the sole superior model overall in terms of both RMSPE (15.38%) and MAPE (11.04%). Unlike in the in-sample context, the GH model achieves greater out-of-sample accuracy than the t model, which attains an RMSPE of 16.33% (almost a full percentage point higher than the GH model) and an MAPE of 11.16% (0.12 percentage point higher than the GH model). Among the special cases of the GH model being considered (NIG, VG, and t), the t model is superior. The VG model obtains a 16.73% RMSPE, and the NIG model obtains a 16.71% RMSPE. It is evident that the GH model family dominates the FMLS model (18.98% RMSPE) and the BS model (19.94% RMSPE) in terms of overall out-of-sample performance. Next, we compare the out-of-sample performance of the models in different parts of the cross-section. The GH model is the best model for OTM, ATM, short-term, and long-term option categories. In particular, it outperforms the next best model by a wide margin for short-term options (by 2.97 percentage points) and for OTM options (by 1.45 percentage points). As in the in-sample setting, the t model fits the ITM options better than the GH model.10 For the medium-term option prices, the VG, t, and GH models perform similarly well. In all option categories, the NIG model is strictly dominated by the GH model. This result highlights the importance of freeing up the index parameter p for optimal out-of-sample performance. The FMLS model does not exhibit any advantage over the GH model family in forecasting the option prices in the given sample. It performs especially poorly in predicting the prices of short-term, OTM puts (43.01% RMSPE). This suggests that the heavy tail assumption under the FMLS model may have introduced unnecessary noise that hinders out-of-sample performance, at least during the sample period being considered. Interestingly, the BS model provides the best out-of-sample fit for long-term, OTM puts (14.34% RMSPE), implying that the additional complexity beyond normality as offered by the FMLS and GH models does not facilitate price prediction for this class of put options. 3.8 Robustness Analysis The above estimation results are based on the dataset from 2012 in which the U.S. stock market was relatively stable. Indeed, the VIX remains below 30 during the sample period. In this section, we want to investigate how the results will change under a more turbulent market condition. The robustness analysis will be based on a sample consisting of daily put option prices from August 4, 2011 to October 13, 2011. The period, which spans 50 trading days, was selected according to the criterion that the VIX Index exceeded and remained above the level 30 over the entire duration. The peak of VIX occurred at 48.00 on August 8, 2011, dubbed the “Black Monday” in the U.S. stock market as it reacted sharply to the Standard and Poor credit agency’s decision to downgrade the credit rating of the U.S. sovereign debt. The sample period is shown in Figure 4. Figure 4 View largeDownload slide The sample periods for main estimation and robustness estimation. Notes: The period used for main estimation spans over 2012 and is shaded in yellow. The period used for robustness estimation spans from August 4 to October 13 (50 trading days) and is shaded in red. The period used for robustness estimation was chosen such that the daily VIX (plotted) stays above 30 continuously over the period. Figure 4 View largeDownload slide The sample periods for main estimation and robustness estimation. Notes: The period used for main estimation spans over 2012 and is shaded in yellow. The period used for robustness estimation spans from August 4 to October 13 (50 trading days) and is shaded in red. The period used for robustness estimation was chosen such that the daily VIX (plotted) stays above 30 continuously over the period. In the robustness analysis, we revisit the out-of-sample model evaluation based on the new dataset. The sampling and estimation strategy are identical as before (using a 5-day rolling window and the cross-section as training samples to estimate the model and predict the next-day put option prices). The robustness estimation results are presented in Table 7 and summarized in Table 8. Given the superior out-of-sample performance of the flexible GH model in Section 3.7 (see Table 8) and to simplify the comparison, we do not repeat generating out-of-sample results for the special cases of the GH model in our robustness analysis. Table 7 Out-of-sample pricing RMSPEs and MAPEs in the robustness estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample is 5 days (on average containing 289 data points) and each testing sample is 1 day ahead of the training sample (on average containing 58 data points). There are 45 testing days. The sample size is n  =  2598. Table 7 Out-of-sample pricing RMSPEs and MAPEs in the robustness estimation Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Out-of-sample RMSPE (%) Out-of-sample MAPE (%) Time-to-maturity Time-to-maturity Short Medium Long All Short Medium Long All Out-of-the-money  BS 27.04 18.52 15.94 19.94 21.63 14.80 12.31 15.51  FMLS model 26.08 17.98 15.84 19.39 19.97 14.49 12.27 15.00  GH model 24.46 17.13 15.26 18.41 19.10 13.88 12.12 14.46 At-the-money  BS 18.00 13.11 14.61 14.59 14.67 10.72 11.76 11.77  FMLS model 15.79 13.42 14.29 14.16 12.60 10.98 11.47 11.43  GH model 14.68 13.96 14.43 14.24 11.67 11.55 11.71 11.62 In-the-money  BS 9.83 11.18 15.35 11.92 8.28 9.24 12.47 9.71  FMLS model 8.12 10.41 14.79 11.05 7.08 8.42 11.84 8.84  GH model 8.46 10.52 14.94 11.21 7.36 8.57 12.16 9.06 All  BS 20.03 14.84 15.25 16.14 15.24 11.78 12.10 12.57  FMLS model 18.65 14.55 14.97 15.58 13.57 11.53 11.83 12.03  GH model 17.58 14.39 14.84 15.21 13.03 11.57 11.95 11.97 Notes: S&P 500 Index put options during 2012 are used. The smallest error measure within each group is italicized and underlined. The bold values show the overall error measures. Each training sample is 5 days (on average containing 289 data points) and each testing sample is 1 day ahead of the training sample (on average containing 58 data points). There are 45 testing days. The sample size is n  =  2598. Table 8 A summary of the superior model(s) based on RMSPE by option category In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH Notes: A summary of the superior model(s) based on RMSPE in each of the option categories for the in-sample estimation and two out-of-sample predictions as reported in Tables 4, 6, and 7. For the out-of-sample (robustness) analysis, the comparison is done between only the BS, FMLS, and GH models, while for the in-sample and out-of-sample (main) analyses, the GH sub-models are also compared. See Table 2 for the definition of the option categories for moneyness and time-to-maturity. Table 8 A summary of the superior model(s) based on RMSPE by option category In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH In-sample Out-of-sample (Main) Out-of-sample (Robustness) Out-of-the-money GH, t GH GH At-the-money GH GH FMLS In-the-money t GH, t FMLS Short-term t GH GH Medium-term GH GH, t GH Long-term GH, t GH GH All GH, t GH GH Notes: A summary of the superior model(s) based on RMSPE in each of the option categories for the in-sample estimation and two out-of-sample predictions as reported in Tables 4, 6, and 7. For the out-of-sample (robustness) analysis, the comparison is done between only the BS, FMLS, and GH models, while for the in-sample and out-of-sample (main) analyses, the GH sub-models are also compared. See Table 2 for the definition of the option categories for moneyness and time-to-maturity. A number of observations are in order. First, the GH model remains the best model in terms of overall out-of-sample performance (15.21% RMSPE and 11.97% MAPE). Second, the GH model outperforms the other models in four individual option categories (OTM, short-term, medium-term, and long-term). In the ATM and ITM option groups, the FMLS model, which allows for infinite kurtosis in the log-return distribution, surpasses the GH model to be the best predictor. Third, focusing on the nine moneyness and time-to-maturity combinations, the GH model is the best model for OTM put options uniformly over all maturity groups, but the FMLS model performs uniformly the best for ITM puts; the model ranking is somewhat mixed for ATM puts across different maturities. From these observations, we learn that, among the models under consideration, not a single one of them eclipses the others in predicting every part of the cross-section of put option prices. The flexibility of the GH model contributes to the best overall out-of-sample performance and in particular to the prediction of the OTM put prices. Nevertheless, under turbulent market conditions, the advantage of the heavy-tailed log-return distribution under the FMLS model is manifested when predicting the ITM put prices, which are sensitive to the extreme left tail of the spot price distribution. We note that all models perform noticeably more poorly during the turbulent period in terms of out-of-sample RMSPE. Furthermore, the differential in the predictive ability of the more complex models relative to BS becomes narrower (e.g., during the normal sample period in the main analysis, the GH and FMLS models outperform the BS model overall by 4.56 and 0.96 percentage points, respectively, in out-of-sample RMSPE; the gaps shrink to 0.93 and 0.56 percentage points, respectively, during the turbulent period). The weaker pricing performance may be attributable to the stringent i.i.d. assumption, which is likely violated during turbulent market conditions. The heightened persistence of high volatilities as perceived by investors in turbulent times is associated with a higher option price; however, this is not directly modeled by the models we considered thus far. This motivates us to consider some of the possible SV extensions to the GH option pricing model, to be discussed in the following section. 4 SV Extensions So far, we have formulated and empirically tested the GH model by assuming stationary and independent increments in the spot price dynamics. In the literature, however, it was documented that the volatility of log-returns displays strong persistence and stochastic behaviors to the extent that the i.i.d. assumption is rejected empirically. This suggests that a GH model that could account for serial dependence in volatility is likely to reduce the pricing errors as observed in the previous empirical study. Building on the time-changed Lévy process modeling framework in Section 1, it is possible to introduce SV in a similar spirit of Carr et al. (2003), Huang and Wu (2004), and Carr and Wu (2004). We start with the baseline GH model, in which the source of randomness is the subordinated Brownian motion B(g(t)), where g(t) is the GIG process. Next, noting that B(g(t)) is a (GH) Lévy process with stationary and independent increments, we may construct a time-changed Lévy process as in Equation (1) by setting L(t)=B(g(t)) and define a time deformation process T(t) that satisfies the three conditions in Section 1. The resultant time-changed Lévy process Z(t)=L(T(t))—a subordinated GH process—is then treated as the random source that drives the log-return process in Equation (3). To induce stochastic volatilities, the process T(t) should exhibit time-dependent dynamics. One way to achieve this is to set T(t) to be the integral of a mean-reverting process. More precisely, it is given by T(t)=∫0th(u)du,dh(t)=κ[1−h(t)]dt+λh(t)dW˜(t). (11) The process h(t) represents the activity rate process which determines the speed at which the business clock runs. It is modeled by the square-root mean-reverting process that solves the stochastic differential equation (11), where κ controls the rate of mean reversion, λ represents the volatility of the activity rate process, and W˜(t) is a standard Wiener process possibly correlated with the standard Wiener process W(t) in Equation (7).11 The long-run activity rate is set to one so that the normalizing restriction on T(t) is satisfied. The characteristic function of T(t) is given by CGMY (2003) φTt(u)=exp(κ2tλ2)+2iuh(0)κ+ψ coth(ψt/2)[cosh⁡(ψt2)+κψsinh⁡(ψt2)]−2κλ2, where ψ=κ2−2λ2iu. This allows us to obtain the characteristic function of the log-price using Equations (2) and (4). The European option prices under the SV flexible GH model (GH–SV model) can then be calculated using Equations (5) and (6). There is an alternative modeling strategy that allows for long-range dependence in SV within the subordinated Brownian motion framework introduced in Section 1.12 Under the GH option pricing model of Finlay and Seneta (2012), serial dependence of stock returns is modeled by a background-driving Lévy process (BDLP) (Barndorff-Nielsen and Shephard, 2001) construction on the subordinator. More specifically, the subordinator g(t) is defined as a superposition of OU processes, each of which is driven by an independent BDLP and has a GIG marginal distribution. The important feature is that the increments of g(t) form an autocorrelated sequence whose marginal distribution remains to be GIG by self-decomposability (Halgreen, 1979). The autocorrelation in the increments of g(t) may induce short- or long-range dependence in squared stock returns. Comparing Finlay and Seneta’s model to our extended GH model, their European option pricing formula involves two-dimensional integrations. We also note that even in the i.i.d. case, the two models have different martingale adjustments and different parameterizations under the risk-neutral measure.13 One future direction of research is to reconcile the models and investigate their empirical performance in the presence of SV. 5 Conclusion Insight to the risk-neutral distribution of logarithmic stock returns is vital to the fitting and prediction of option prices. In this paper, we develop a flexible, four-parameter GH option pricing model in the time-changed Lévy process framework, and present six three-parameter special cases: the VG, t, H, RH, NIG, and NRIG option pricing models. In respect of the seven models’ properties, the flexible GH model and its special cases generalize the BS model by allowing the passage of economic time to be stochastic. As such, the GH model family forms a class of time-subordinated models that can cope with yet another facet of the unpredictable financial market. In addition, the subordination to a drifted Brownian motion entails that the class of GH processes are able to capture excess kurtosis and skewness. Using the characteristic function approach of Carr and Madan (1999), we obtain an explicit pricing formula for European options under the baseline GH model. The option pricing formula takes a convenient semi-closed form. With judicious choices of parameterization and variable transformation, all the four model parameters can be estimated freely. These features greatly facilitate the empirical study on the GH option pricing model. Using S&P 500 Index options, we conduct an extensive empirical study on the GH option pricing model and some of its special cases, and compare the empirical performance against some benchmark option pricing models that assume i.i.d. spot returns, including the BS model and the FMLS model (Carr and Wu, 2003). Our empirical analysis suggests that the GH model generally delivers the best in-sample fit and out-of-sample prediction for all put option prices in the cross-section. The same conclusion holds for the robustness analysis in which the sample spans a more turbulent period, although the heavy-tail FMLS option pricing model yields more accurate out-of-sample forecast for some specific option categories. Among the NIG, VG, and t special cases, the t model seems to be the most tenable model for pricing options. Remarkably, the t model’s average in-sample fit is better than that of the VG model (Madan, Carr, and Chang, 1998) for all option types—a result which is corroborated by the GH model’s parameter estimates. In this paper, we concentrate on the empirical study of the baseline GH model, which belongs to a class of the Lévy-based option pricing models that assumes independent and stationary increments in the underlying spot price process. However, the i.i.d. assumption on the spot returns can be crippling in view of the bulk of well-known stylized facts regarding financial asset returns. Leveraging the proliferous time-changed Lévy process framework (CGMY, 2003; Carr and Wu, 2004; Huang and Wu, 2004), we present an SV extension to the baseline model that employs a mean-reverting process for time deformation. Empirical study of this extension is intended for future work. Appendix A: The Characteristic Function of the Logarithmic Stock Price First, we derive the value of the martingale correction factor, denoted by ωt, as −log(φZt(−i)). Assuming no arbitrage and by risk-neutral pricing, the discounted stock price process {Ste−(r−q)t} follows a martingale under the risk-neutral probability measure Q. This means that, for any positive t, S0=EQ[Ste−(r−q)t|F0], (A.1) where F0 is the information set at time 0, and E Q[·|F0] represents the conditional expectation taken with respect to Q. Given the log stock price model in Equation (3), the martingale condition is satisfied by applying a martingale correction factor ωt to the drift of the log-return process, as follows: S0=EQ[S0 exp{(r−q)t+ωt+Zt−(r−q)t}|F0]. (A.2) Solving for ωt, we obtain: 1=EQ[exp{ωt+Zt}|F0]1=eωtEQ[eZt|F0]e−ωt=EQ[eZt|F0]ωt=−log(EQ[eZt|F0])ωt=−log(φZt(−i)) (A.3) This martingale correction factor is the same as the one in CGMY (2003, equation (4.7), p. 358) using ordinary exponential transform. We then obtain Equation (4), by expressing the characteristic function of the log spot price φ log St(u) in terms of φZt(u), the characteristic function of Z(t), as follows: φ log St(u)=EQ[exp{iu log(St)}|F0]=EQ[Stiu|F0]=EQ[S0iu exp{iu[(r−q)t−log(φZt(−i))+Zt]}]=S0iu exp{iu(r−q)t−log(φZt(−i))}φZt(u). (A.4) Appendix B: Derivations for the Flexible GH, VG, and t Option Pricing Models B.1 The Flexible GH Model B.1.1 Characteristic function used for option pricing The characteristic function of the GIG distribution is given by φg1(u)=(γ2γ2−2iu)p2Kp(δ2(γ2−2iu))Kp(γδ)=(1−2γ2iu)−p2Kp(ζ1−2γ2iu)Kp(ζ), (A.5) where ζ=γδ, γ2 is given by Equation (8) and Kj(·) is the modified Bessel function of the third kind with index j. Using Equation (2), we obtain φZ1(u)=φg1(uθ+iσ2u22)=(1−2γ2(iuθ−12σ2u2))−p2Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ). (A.6) It follows from Equation (4) that φ log St(u)=S0iu exp{iu(r−q+ω)t}φX1(u)t=S0iu exp{iu(r−q+ω)t}[1−2γ2(iuθ−12σ2u2)]−pt2×[Kp(ζ1−2γ2(iuθ−12σ2u2))Kp(ζ)]t, (A.7) where the unit time drift adjustment can be computed as ω=−log φX1(−i)=−p2 log[(1−2γ2(θ+12σ2))]+log[Kp(ζ1−2γ2(θ+12σ2))Kp(ζ)], (A.8) for θ<(γ22−σ22). B.1.2 Distributional properties of Z1 under the GH model Let us start by studying the statistical properties of Z1. First, we derive the probability density function of the GH distribution, which is the distribution of the random variable Z1 in the log-price specification Equation (3) under the flexible GH option pricing model. Let N (x;μ,σ2) denote the density of a normal distribution with mean μ and variance σ2. Suppose, conditional on g1, that X1 follows a normal distribution with mean θg1 and variance σ2g1, and that g1 follows a GIG with probability density function fg1(g), given by GIG(p,δ,γ):                 fg1(g)=(γδ)p2Kp(γδ)gp−1 exp{−12(γ2g+δ2g)}, (A.9) where Kj(·) is the modified Bessel function of the third kind with index j. Then, Z1 will follow a GH distribution with its density function fZ1(x) given by, for any real value x, fZ1(x)=∫0∞N(x;θg,σ2g)fg1(g)dg=∫0∞12πσ2g exp{−(x−θg)22σ2g}·(γδ)p2Kp(γδ)gp−1 exp{−12(γ2g+δ2g)}dg=(γδ)p22πσ2Kp(γδ)∫0∞gp−32 exp{−12[(x−θg)2σ2g+γ2g+δ2g]}dg (A.10) =(γδ)p22πσ2Kp(γδ)eθσ2x∫0∞gp−32 exp{−12[(θ2σ2+γ2)g+(x2σ2+δ2)1g]}dg=(γδ)p2πσ2Kp(γδ)eθσ2x(x2σ2+δ2θ2σ2+γ2)p2−14Kp−12((θ2σ2+γ2)(x2σ2+δ2)). The last equality uses the fact that, for all η,υ>0, ∫0∞xj−1e−12(ηx+υx)dx=2(υη)j2Kj(ηυ). The parameter restriction Equation (8) that γ2=ζKp+1(ζ)Kp(ζ) (where ζ=γδ) coming from the normalization E [g1]=1 still holds under the time-changed Lévy process framework. The first four standardized central moments of Z1 under the GH model are given by E(Z1)=θδγKp+1(δγ)Kp(δγ),Var(Z1)=σ2δγKp+1(δγ)Kp(δγ)+θ2δ2γ2[Kp+2(δγ)Kp(δγ)−(Kp+1(δγ)Kp(δγ))2],skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2{θ3δ3γ3[Kp+3(δγ)Kp(δγ)−3Kp+2(δγ)Kp+1(δγ)Kp2(δγ)+2(Kp+1(δγ)Kp(δγ))3]+3σ2θδ2γ2[Kp+2(δγ)Kp(δγ)−(Kp+1(δγ)Kp(δγ))2]},kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=Var(Z1)−2{θ4δ4γ4[Kp+4(δγ)Kp(δγ)−4Kp+3(δγ)Kp+1(δγ)Kp2(δγ)+6Kp+2(δγ)Kp+12(δγ)Kp3(δγ)−3(Kp+1(δγ)Kp(δγ))4]+σ2θ2δ3γ3[6Kp+3(δγ)Kp(δγ)−12Kp+2(δγ)Kp+1(δγ)Kp2(δγ)+6(Kp+1(δγ)Kp(δγ))3]+3σ4δ2γ2Kp+2(δγ)Kp(δγ)}. (A.11) Applying the normalization E[g1]=1, and denoting ζ=δγ, we obtain the parameter restriction γ2=ζKp+1(ζ)Kp(ζ). The four standardized central moments of Z1 then become E(Z1)=θ,Var(Z1)=σ2+θ2[Kp+2(ζ)Kp(ζ)Kp+12(ζ)−1],skew(Z1)=Var(Z1)−3/2{θ3[Kp+3(ζ)Kp2(ζ)Kp+13(ζ)−3Kp+2(ζ)Kp(ζ)Kp+12(ζ)+2]+3σ2θ[Kp+2(ζ)Kp(ζ)Kp+12(ζ)−1]},kur(Z1)=Var(Z1)−2{θ4[Kp+4(ζ)Kp3(ζ)Kp+14(ζ)−4Kp+3(ζ)Kp2(ζ)Kp+13(ζ)+6Kp+2(ζ)Kp(ζ)Kp+12(ζ)−3]+σ2θ2[6Kp+3(ζ)Kp2(ζ)Kp+13(ζ)−12Kp+2(ζ)Kp(ζ)Kp+12(ζ)+6]+3σ4Kp+2(ζ)Kp(ζ)Kp+12(ζ)}. (A.12) Note that in the symmetric case of θ = 0, we have kur(Z1)>3, by the Turán-type inequality Kp+2(y)Kp(y)>Kp+12(y) (Ismail and Muldoon, 1978). B.1.3 Distributional properties of Zt under the GH model Next, let us study the statistical properties of Zt. From the infinitely divisibility of a general Lévy process Zt, the characteristic functions of Zt and Z1 are related by φZt(u)=[φZ1(u)]t. It follows that their cumulant generating functions (i.e., the log-characteristic functions) are related by ΨZt(u)=tΨZ1(u). By repeated differentiation, the following relations regarding the cumulants of Zt and Z1 hold for all positive integers n: cn(Zt):=1in∂n∂unΨZt(0)=t1in∂n∂unΨZ1(0)=:tcn(Z1). Using the links between the cumulants and the central moments, the first four central moments of Zt and Z1 can be related as follows: E(Zt)=c1(Zt)=tc1(Z1)=tE(Z1),Var(Zt)=c2(Zt)=tc2(Z1)=tVar(Z1),E[(Zt−E(Zt))3]=c3(Zt)=tc3(Z1)=tE[(Z1−E(Z1))3],E[(Zt−E(Zt))4]=c4(Zt)+3c2(Zt)2=tc4(Z1)+3[tc2(Z1)]2=tE[(Z1−E(Z1))4]+3t(t−1)Var(Z1)2. (A.13) The skewness and kurtosis of Zt are then given by: skew(Zt)=E[(Zt−E(Zt))3]Var(Zt)3/2=tE[(Z1−E(Z1))3][tVar(Z1)]3/2=1tskew(Z1),kur(Zt)=E[(Zt−E(Zt))4]Var(Zt)2=tE[(Z1−E(Z1))4]+3t(t−1)Var(Z1)2[tVar(Z1)]2=1tkur(Z1)+3(1−1t). (A.14) We thus see that, for a general Lévy process, its mean and variance increase with t, its skewness diminishes with t at the square-root rate, and its kurtosis approaches that of a normal distribution in the limit as t→∞. In the symmetric case of θ = 0, we see that kur(Z1)>3, which implies that kur(Zt)>3 for all finite t. B.2 The VG Model and its Statistical Properties B.2.1 Characteristic function used for option pricing The characteristic function of the gamma distribution is given by φg1(u)=(1−iuν)−1ν such that φX1(u)=[1−ν(iuθ−12σ2u2)]−1ν. (A.15) It follows then that the VG option price can be computed using φ log St(u)=S0iu exp{iu(r−q+ω)t}[1−ν(iuθ−12σ2u2)]−tν, (A.16) with unit time drift adjustment, ω=1ν log[1−ν(θ+12σ2)], for θ<(1ν−σ22). B.2.2 Distributional properties of Z1 under the VG model Let us study the statistical properties of Z1. The density function of Z1 under the VG model can be obtained from Equation (A.10) by letting δ→0 and setting p=1ν>0, so that γ=2p=2ν. fZ1(x)=(γδ)1ν2πσ2Γ(1ν)21ν−1(γδ)−1νeθσ2x(x2σ2θ2σ2+2ν)12ν−14K1ν−12((θ2σ2+2ν)x2σ2)=γ2ν2πσ2Γ(1ν)21ν−1eθσ2x(x2θ2+2σ2ν)12ν−14K1ν−12(xσ2θ2+2σ2ν). (A.17) This recovers the density function in Madan, Carr, and Chang (1998, equation (23)). Imposing the same limits on the parameters, we obtain the first four standardized central moments under the VG model for ν>0: E(Z1)=θ,Var(Z1)=σ2+θ2ν,skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2(2θ3ν2+3σ2θν)kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=3Var(Z1)−2[θ4ν2(1+2ν)+2σ2θ2ν(1+2ν)+σ4(1+ν)]. (A.18) In the symmetric case of θ = 0, the kurtosis is always greater than 3 as ν>0. B.3 The t Model B.3.1 Characteristic function used for option pricing The characteristic function of the reciprocal gamma distribution, R Γ(a,b) is given by φg1(u)=2(−biu)a/2Γ(a)Ka[−4biu]=2(1ν−1)12νΓ(1ν)(−iu)12νK1ν[−4(1ν−1)iu], (A.19) where a=1ν and b=1ν−1 for 0<ν<1. Kj(·) is the modified Bessel function of the third kind with index j. Using the same procedure as for the GH and VG derivations, the characteristic function of X1 is given by φX1(u)=2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)]. (A.20) The characteristic function of log St is thus given by φ log St(u)=S0iu exp{iu(r−q+ω)t}×[2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν(−4(1ν−1)(iuθ−12σ2u2))]t, (A.21) and the unit time drift adjustment, ω=−log[2(1ν−1)12νΓ(1ν)[−(θ+12σ2)]12νK1ν(−4(1ν−1)(θ+12σ2))], (A.22) for θ<−σ22. B.3.2 Distributional properties of Z1 under the t model The density function of Z1 under the t model can be obtained from Equation (A.10) by letting γ→0 and setting p=−1ν<0, so that δ=2(−p−1)=2(1ν−1). fZ1(Z)=(1ν−1)1νσ2πΓ(1ν)eθσ2Z(Z2+2(1ν−1)σ2θ)−1ν−12K1ν+12(θσ2Z2+2(1ν−1)σ2). (A.23) By the variable and parameter transformation, y=Z/σ+μ, β=θ/σ, α=β2+γ2, λ=−1/ν, and δ=2(1ν−1), we recover the density function of the GH–skew–t distribution in Aas and Haff (2006, equation (3)). Imposing the same limits on the parameters, we obtain the first four standardized central moments under the t model and the associated domains of ν in which the moments exist: E(Z1)=θ,Var(Z1)=σ2+θ2ν1−2ν for ν∈(0,12),skew(Z1)=E[(Z1−E(Z1))3]Var(Z1)3/2=Var(Z1)−3/2[θ34ν2(1−2ν)(1−3ν)+3σ2θν1−2ν] for ν∈(0,13),kur(Z1)=E[(Z1−E(Z1))4]Var(Z1)2=3Var(Z1)−2[θ4ν2(1+5ν)(1−2ν)(1−3ν)(1−4ν)+2σ2θ2ν(1+ν)(1−2ν)(1−3ν)+σ41−ν1−2ν] for ν∈(0,14). (A.24) In the symmetric case of θ = 0, the kurtosis is always greater than 3 as ν>0. Appendix C: The VG and t Limiting Cases of the GH Model The following lemma shows how we may obtain the characteristic functions of Z1 under the VG and t models as appropriate limits of that of Z1 under the GH model. Lemma 1 Suppose p=1ν>0, δ→0 and γ=2ν. Then, the GH model reduces to the VG model. Suppose p=−1ν<−1, γ→0 and δ=2(1ν−1). Then, the GH model reduces to the t model. Proof Either δ→0 or γ→0 implies that ζ→0. Making use of the properties of the modified Bessel function of the third kind with index j, Kj(·), that Kj(w)∼Γ(|j|)2|j|−1w−|j|    for w↓0, and that Kj(w)=K−j(w), the characteristic function of X1 under the GH model (given by Equation (A.6)) can be further simplified as follows: Under the assumptions in (i), we have φX1(u)=(1−ν(iuθ−12σ2u2))−12νlim⁡δ→0Γ(1ν)21ν−1(γδ1−ν(iuθ−12σ2u2))−1νΓ(1ν)21ν−1(γδ)−1ν=(1−ν(iuθ−12σ2u2))−12ν(1−ν(iuθ−12σ2u2))−12ν=(1−ν(iuθ−12σ2u2))−1ν, which is the corresponding characteristic function under the VG model, given in Equation (A.15). Under the assumptions in (ii), we have φX1(u)=lim⁡γ→0(1−2γ2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))K1ν(δγ)=lim⁡γ→0(1−2γ2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))Γ(1ν)21ν−1(γδ)−1ν=lim⁡γ→0δ1νΓ(1ν)21ν−1(γ2−2(iuθ−12σ2u2))12νK1ν(δγ2−2(iuθ−12σ2u2))=δ1νΓ(1ν)21ν−1(−2(iuθ−12σ2u2))12νK1ν(δ−2(iuθ−12σ2u2)) =212ν(1ν−1)12νΓ(1ν)21ν−1212ν[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)]=2(1ν−1)12νΓ(1ν)[−(iuθ−12σ2u2)]12νK1ν[−4(1ν−1)(iuθ−12σ2u2)], which is the corresponding characteristic function under the t model, given in Equation (A.20). Footnotes * We deeply thank Fabio Trojani, the co-editor, and two anonymous referees for their constructive suggestions which led to substantial improvement on the paper. We also thank those who gave valuable comments to our work at the EcoSta2017 Conference and at a University of Sydney Business School seminar. All remaining errors are ours. 1 The pure-jump and infinite activity property of the GH process and special cases is a distinguishing property from Heston’s (1993) model, the BS model, and from other seminal models, such as Merton’s (1976) jump-diffusion model, Bates, (1996) model, and Pan’s (2002) model, as noted by Carr and Wu (2004). As a purely discontinuous process, the GH process is also different from the GH diffusion process of Bibby and Sørensen (1997) and Rydberg (1999). 2 The use of time deformation processes that induce persistent variations in drift and volatility is supported by a general equilibrium argument. For example, see Shaliastovich and Tauchen (2011). We thank one referee who brings this literature to our attention. 3 Admitting the computational difficulty in estimating the four parameters jointly, Prause (1999) only presents the risk-neutral estimation results for the NIG model, a special case of the GH option pricing model obtained by fixing the index parameter [or shape parameter according to Prause (1999)] to p=−1/2 [see Table 2.28 and the discussion on p. 61 of Prause (1999) for details]. Nevertheless, Eberlein and Prause (2002) and Prause (1999) estimate the GH model parameters under the physical measure using spot data. 4 The FMLS model of Carr and Wu (2003) features a heavy-tailed, negatively skewed log-return distribution. While not a time-changed Lévy process model, the FMLS model is chosen as a benchmark model rather than the log-stable model of Hurst, Platen, and Rachev (1997) because the former guarantees finite moments of the spot price and hence finite option prices. 5 In general, when there are a finite number of assets in the market, the risk-neutral probability measure Q may not be unique. When the price of the underlying asset exhibits jumps, it typically requires infinite many assets to complete the market (an exception is when the asset price contains Poisson jumps with a finite number of different jump sizes). On the technical level, sufficient and necessary conditions for market completeness are known for a general Lévy process (Sato, 1999). More recently, Li and Mendoza-Arriaga (2015) provide a set of sufficient conditions for equivalent measure changes for a subordinated diffusion model. 6 A similar version of the t option pricing model was proposed in Yeap (2014), except that a different parameterization ν=Var(g1) was adopted. In this paper, we enforce the normalization E[g1]=1 so that Var(g1)=ν1−2ν. Our model allows for skewness in the log-return distribution, whereas a symmetric t option pricing model has been proposed (Cassidy, Hamp, and Ouyed, 2010). 7 The log-returns (net of the risk-free rate and dividend yield) follow a symmetric normal distribution under BS and a negatively skewed stable distribution under the FMLS model (with skewness parameter −1). 8 Alternatively, we may smooth out the cyclicity by using a dynamically weighted sample in which different weights are assigned to options with different time-to-maturities (Pan, 2002). The weights are determined such that the average time-to-maturity of the weighted sample is maintained at a constant level. The advantage of our current method (by retaining the raw data) is that the term structure of the risk-neutral distribution, and hence the relationship between the shape parameters and the time-to-maturity becomes more apparent. 9 In particular, the empirical skewness and kurtosis of log-returns as implied from the mean parameter estimates under the t model are equal to −0.25 and 4.04, respectively, which are both bigger in magnitude than those under the VG model (−0.09 and 3.33). The empirical mean and variance are similar under the two models (−0.04 and 0.02 under the VG model; −0.05 and 0.03 under the t model). 10 We note that some pricing improvements are observed in the out-of-sample context compared with the in-sample context for ITM options (e.g., ITM option prices as fit by the VG model). This can be explained by the relatively small representation of ITM options (16%) in the data compared with ATM (41%) and OTM (43%) options (see Table 2), such that the ITM options have a relatively weak influence over the in-sample fitting. 11 We may set W˜(t)=ρW(t)+1−ρ2W⌣(t), where W(t) and W⌣(t) are independent standard Wiener processes. 12 The option pricing models for assets with long-range dependence were pioneered by Heyde (1999), and Heyde and Liu (2001). Special cases of the GH distribution have since been employed in short- and long-range dependence models, such as the t distribution by Heyde and Leonenko (2005); Finlay and Seneta (2006); and Leonenko, Petherick, and Sikorskii (2011); the VG distribution by Finlay and Seneta (2006) and Leonenko, Petherick, and Sikorskii (2012b); and the NIG distribution by Leonenko, Petherick, and Sikorskii (2012a,b). 13 Unlike our model, Finlay and Seneta (2012) do not impose the normalization restriction E[g1]=1 on the subordinator. Instead, their martingale condition implies a restriction on the skewness parameter θ=−12σ2, while θ is a free parameter under our model. References Aas K. , Haff I. H. . 2006 . The Generalised Hyperbolic Skew–Student’s t Distribution . Journal of Financial Econometrics 4 : 275 – 309 . Google Scholar CrossRef Search ADS Andersen T. G. , Fusari N. , Todorov V. . 2015 . The Risk Premia Embedded in Index Options . Journal of Financial Economics 117 : 558 – 584 . Google Scholar CrossRef Search ADS Bakshi G. , Cao C. , Chen Z. . 1997 . Empirical Performance of Alternative Option Pricing Models . Journal of Finance 52 : 2003 – 2049 . Google Scholar CrossRef Search ADS Bakshi G. , Madan D. . 2000 . Spanning and Derivative-Security Valuation . Journal of Financial Economics 55 : 205 – 238 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. 1977 . Exponentially Decreasing Distributions for the Logarithm of Particle Size . Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 353 : 401 – 419 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Halgreen C. . 1977 . Infinite Divisibility of the Hyperbolic and Generalised Inverse Gaussian Distributions . Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 38 : 309 – 311 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Shephard N. . 2001 . Non-Gaussian OU Based Models and some of their Uses in Financial Economics . Journal of Royal Statistical Society: Series B 63 : 167 – 241 . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Mikosch T. , Resnick S. . 2001 . Lévy Processes—Theory and Applications . Boston : Birkhauser . Google Scholar CrossRef Search ADS Barndorff-Nielsen O. E. , Shephard N. . 2012 . “Basics of Lévy Processes.” Economics Series Working papers, University of Oxford (No. 610) . Bates D. S. 1996 . Jumps and Stochastic Volatility: Exchange Rate Processes Implicit in Deutsche Mark Options . Review of Financial Studies 9 : 69 – 107 . Google Scholar CrossRef Search ADS Bibby B. M. , Sørensen M. . 1997 . A Hyperbolic Diffusion Model for Stock Prices . Finance & Stochastics 1 : 25 – 41 . Google Scholar CrossRef Search ADS Black F. , Scholes M. . 1973 . The Pricing of Options and Corporate Liabilities . Journal of Political Economy 81 : 637 – 654 . Google Scholar CrossRef Search ADS Carr P. , Madan D. . 1999 . Option Valuation Using Fast Fourier Transform . Journal of Computational Finance 2 : 61 – 73 . Google Scholar CrossRef Search ADS Carr P. , Wu L. . 2003 . The Finite Moment Log Stable Process and Option Pricing . Journal of Finance 58 : 753 – 777 . Google Scholar CrossRef Search ADS Carr P. , Wu L. . 2004 . Time-Changed Lévy Processes and Option Pricing . Journal of Financial Economics 71 : 113 – 141 . Google Scholar CrossRef Search ADS Carr P. , Geman H. , Madan D. B. , Yor M. . 2003 . Stochastic Volatility for Lévy Processes . Mathematical Finance 13 : 345 – 382 . Google Scholar CrossRef Search ADS Cassidy D. T. , Hamp M. J. , Ouyed R. . 2010 . Pricing European Options with a Log Student’s t-distribution: A Gosset Formula . Physica A 389 : 5736 – 5748 . Google Scholar CrossRef Search ADS Christie A. A. 1982 . The Stochastic Behavior of Common Stock Variances: Value, Leverage and Interest Rate Effects . Journal of Financial Economics 10 : 407 – 432 . Google Scholar CrossRef Search ADS Christoffersen P. , Jacobs K. . 2004 . The Importance of the Loss Function in Option Valuation . Journal of Financial Economics 72 : 291 – 318 . Google Scholar CrossRef Search ADS Clark P. K. 1973 . A Subordinated Stochastic Process Model with Finite Variance for Speculative Prices . Econometrica 41 : 135 – 155 . Google Scholar CrossRef Search ADS Daal E. A. , Madan D. B. . 2005 . An Empirical Examination of the Variance–Gamma model for Foreign Currency Options . The Journal of Business 78 : 2121 – 2152 . Google Scholar CrossRef Search ADS Eberlein E. , Keller U. . 1995 . Hyperbolic Distributions in Finance . Bernoulli 1 : 281 – 299 . Google Scholar CrossRef Search ADS Eberlein E. , Prause K. . 2002 . “The Generalized Hyperbolic Model: Financial Derivatives and Risk Measures.” In Mathematical Finance—Bachelier Congress 2000 , pp. 245 – 267 . Heidelberg : Springer-Verlag . Google Scholar CrossRef Search ADS Finlay R. , Seneta E. . 2006 . Stationary-Increment Student and Variance–Gamma Processes . The Journal of Applied Probability 43 : 441 – 453 . Google Scholar CrossRef Search ADS Finlay R. , Seneta E. . 2012 . A Generalized Hyperbolic Model for a Risky Asset with Dependence . Statistics and Probability Letters 82 : 2164 – 2169 . Google Scholar CrossRef Search ADS Gilli M. , Schumann E. . 2012 . Calibrating Option Pricing Models with Heuristics . Natural Computing in Computational Finance 380 : 9 – 37 . Google Scholar CrossRef Search ADS Gruber P. H. , Tebaldi C. , Trojani F. . 2015 . “The Price of the Smile and Variance Risk Premia.” Swiss Finance Institute Research Paper No. 15-36 . Halgreen C. 1979 . Self-decomposability of the Generalized Inverse Gaussian and Hyperbolic Distributions . Zeitschrift fur Wahrscheinlichkeitstheorie und Verwandte Gebiete 47 : 13 – 17 . Google Scholar CrossRef Search ADS Heston S. L. 1993 . A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options . The Review of Financial Studies 6 : 327 – 343 . Google Scholar CrossRef Search ADS Heyde C. C. 1999 . A Risky Asset Model with Strong Dependence through Fractal Activity Time . The Journal of Applied Probability 36 : 1234 – 1239 . Google Scholar CrossRef Search ADS Heyde C. C. , Leonenko N. N. . 2005 . Student Processes . Advances in Applied Probability 37 : 342 – 365 . Google Scholar CrossRef Search ADS Heyde C. C. , Liu S. . 2001 . Empirical Realities for a Minimal Description Risky Asset Model. The Need for Fractal Features . Journal of Korean Mathematical Society 38 : 1047 – 1059 . Huang J. , Wu L. . 2004 . Specification Analysis of Option Pricing Models Based on Time-Changed Levy Processes . The Journal of Finance 59 : 1405 – 1439 . Google Scholar CrossRef Search ADS Hurst S. R. , Platen E. , Rachev S. . 1997 . Subordinated Market Index Models: A Comparison . Financial Engineering and the Japanese Markets 4 : 97 – 124 . Google Scholar CrossRef Search ADS Ismail M. , Muldoon M. . 1978 . Monotonicity of the Zeros of a Cross-Product of Bessel Functions . SIAM Journal on Mathematical Analysis 9 : 759 – 767 . Google Scholar CrossRef Search ADS Jørgensen B. 1982 . Statistical Properties of the Generalised Inverse Gaussian Distribution . New York : Springer-Verlag . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2011 . The Student Subordinator Model with Dependence for Risky Asset Returns . Communications in Statistics—Theory and Methods 40 : 3509 – 3523 . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2012a . A Normal Inverse Gaussian Model for a Risky Asset with Dependence . Statistics and Probability Letters 82 : 109 – 115 . Google Scholar CrossRef Search ADS Leonenko N. N. , Petherick S. , Sikorskii A. . 2012b . Fractal Activity Time Models for Risky Asset with Dependence and Generalised Hyperbolic Distributions . Stochastic Analysis and Applications 30 : 476 – 492 . Google Scholar CrossRef Search ADS Li L. , Mendoza-Arriaga R. . 2015 . “Equivalent Measure Changes for Subordinate Diffusions.” Working paper . Madan D. B. , Carr P. P. , Chang E. C. . 1998 . The Variance Gamma Process and Option Pricing . European Finance Review 2 : 79 – 105 . Google Scholar CrossRef Search ADS Madan D. B. , Milne F. . 1991 . Option Pricing with VG Martingale Components . Mathematical Finance 1 : 39 – 55 . Google Scholar CrossRef Search ADS Merton R. C. 1973 . The Theory of Rational Option Pricing . Bell Journal of Economics and Management Science 4 : 141 – 183 . Google Scholar CrossRef Search ADS Merton R. C. 1976 . Option Pricing When Underlying Stock Returns Are Discontinuous . Journal of Financial Economics 3 : 125 – 144 . Google Scholar CrossRef Search ADS Pan J. 2002 . The Jump-Risk Premium Implicit in Options: Evidence from an Integrated Time-Series Study . Journal of Financial Economics 63 : 3 – 50 . Google Scholar CrossRef Search ADS Prause K. 1999 . The Generalized Hyperbolic Model: Estimation, Financial Derivatives, and Risk Measures . Thesis, Albert–Ludwigs–Universitẗ Freiburg . Rubinstein M. 1985 . Non-Parametric Tests of Alternative Option Pricing Models using All Reported Trades and Quotes on the 30 Most Active CBOE Option Classes from August 23, 1976 through August 31, 1978 . Journal of Finance 40 : 455 – 480 . Google Scholar CrossRef Search ADS Rydberg T. 1999 . Generalized Hyperbolic Diffusion Processes with Applications in Finance . Mathematical Finance 9 : 183 – 201 . Google Scholar CrossRef Search ADS Scott L. 1997 . Pricing Stock Options in a Jump-Diffusion Model with Stochastic Volatility and Interest Rates: Applications of Fourier Inversion Methods . Mathematical Finance 7 : 413 – 426 . Google Scholar CrossRef Search ADS Sato K. 1999 . Lévy Processes and Infinitely Divisible Distributions . Cambridge : Cambridge University Press . Shaliastovich I. , Tauchen G. . 2011 . Pricing of the Time-Change Risks . Journal of Economic Dynamics and Control 35 : 843 – 858 . Google Scholar CrossRef Search ADS White H. 1980 . A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity . Econometrica 48 : 817 – 838 . Google Scholar CrossRef Search ADS Yeap C. 2014 . The Skew-t Option Pricing Model . Thesis , The University of Sydney Business School . © The Author, 2017. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Journal of Financial EconometricsOxford University Press

Published: Oct 17, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off