Subordination, Self-similarity, and Option Pricing

Recommended by Henry Schellhorn We use additive processes to price options on the Standard and Poor's 500 index SPX for the sake of comparison of pricing performance across both model class and family of time-one distribution. Each of the additive processes in this study is defined using one of the following: subordination, Sato's 2002 construction of self-similar additive processes from self-decomposable distributions, or both. We find that during the year 2005: 1 for a given family of time-one distributions, four-parameter self-similar additive models consistently yielded lower pricing errors than those of four-parameter subordinated, and time-inhomogeneous additive models, 2 for a given class of additive models, the time-one marginal given by the normal inverse Gaussian distribution consistently yielded lower pricing errors than those of the variance gamma distribution. Market and model benchmarks for the additive models under consideration are obtained via the bid-ask spreads of the options and Lévy stochastic volatility model prices, respectively.


Introduction
As noted by Carr and Wu 1 , Cont and Tankov 2, page 453 , and others, Lévy models are incapable of adequately fitting implied volatility surfaces of equity options across both strike and maturity.Although stochastic volatility models with jumps have success in pricing across both strike and maturity 3 , 2, page 470 , 4, page 98 , the use of such models typically incurs the cost of computing with a significantly larger number of parameters than in the Lévy case.By using certain additive models with nonstationary increments, one obtains flexibility in pricing across both strike and maturity, at the cost of one more parameter than a Lévy model of like time-one distribution.
One class of time-inhomogeneous additive processes which has recently received much attention is the class of the H-self-similar additive processes of Sato 5, Theorem 16.1 ii , page 99 .In 1998 Barndorff-Nielsen 6 used Sato's results to construct a self-similar, time-inhomogeneous additive process with marginal normal inverse Gaussian distributions.The financial application of self-similar additive processes determined from selfdecomposable distributions was later recommended by Nolder 7 .In 2007, Carr, Geman, Madan, and Yor CGMY 8 introduced the class of exponential H-self-similar additive models and showed that these could accurately price SPX and equity options across both strike and maturity with as few as four parameters.In their work, exponential self-similar additive models were constructed and calibrated using option data on several underlying assets.In other words, the self-similar additive model class was fixed, and the in-sample model pricing error was observed with respect to both the choice of the family of time-one distribution and the choice of underlying asset.In contrast, our study fixes the underlying asset SPX and observes pricing errors with respect to both the choice of the family of timeone distribution and choice of the additive model class.
The additive model classes used in this study correspond to the Lévy processes, Hself-similar additive processes, and those in which a Brownian motion is time-changed by an independent, increasing H-self-similar additive process with time-one self-decomposable distribution.For each of the three additive classes, two processes are constructed; one in which the time-one marginal is a variance gamma distribution, and the other, a normal inverse Gaussian distribution.The names for each of the additive processes to be compared are given in Table 1.
The time-inhomogeneous additive processes discussed here require only one more parameter than a Lévy model of like time-one distribution.This additional parameter is known as the Hurst exponent, H.In order that this research may be viewed in historical context, a brief history of option pricing with additive processes now follows.

Brief history of additive processes in option pricing
The use of additive processes to price options dates back to 1900 when Louis Bachelier, a student of Henri Poincaré, began the development of the theory of Brownian motion 9 .In his dissertation, Théorie de la Spéculation, changes in the price of a stock were modeled with what is now known as a zero mean, Brownian motion 10 , a 1/2-self-similar additive process.Mandelbrot's exponential stable-Paretian processes 10 and Samuelson's "economic" geometric Brownian motion 11 overcame the problem of negative stock prices associated with a Brownian motion model.The latter of the two processes was found extremely useful to the theory of option pricing developed by Sprenkle 12 , Samuelson 11 , and others.Finally, in 1973 the famous closed-form solution for the price of a European call option was introduced by Black and Scholes 13 .In the same year Merton 14 presented another derivation of the Black-Scholes formula, using weaker assumptions than in 13 , and extended the theory to include the case of stochastic interest rates.
Not long after the introduction of the Black-Scholes model, Merton 15 , and later, Rubinstein 16 , Bates 17 , and others documented inconsistencies of this model with market option prices.Merton 15 introduced the jump-diffusion model in which the logarithm of the return on the underlying asset follows a Lévy process whose one-dimensional marginals are given by the convolution of a Normal distribution with a compound Poisson distribution.Madan and Seneta 18 resp., Madan et al. 19 introduced the two resp., three parameter subordinated Lévy process, known as the variance gamma process, in which a Brownian motion without resp., with drift is time changed by an independent gamma process.Barndorff-Nielsen 6 introduced the three-parameter normal inverse Gaussian process in which a Brownian motion is time changed by an independent inverse Gaussian process.Other well-known Lévy processes include the following: Meixner 20 , KoBoL 21, 22 or equivalently, CGMY 23 , hyperbolic 24, 25 , generalized hyperbolic 25, 26 , and nonparametric Lévy of Cont and Tankov 27 .
Although Lévy processes with jumps could be used to accurately price options across strike price for a given maturity, they were inadequate to do so across both strike and maturity 1 and 2, page 453 .One answer to this problem came in 2003 when CGMY 3 showed that one may subordinate a Lévy process by the time integral of a mean-reverting, positive process in order to obtain a model which was quite flexible across both strike and maturity.Such models have been termed "Lévy stochastic volatility" models.Another answer to this problem came with the simpler Finite Moment Log Stable FMLS model of Carr and Wu 1 .In this model the logarithm of the return on the underlying followed a stable Paretian process in which the skewness parameter, β, was set to an extreme value of −1.With β fixed and the drift term set so that the martingale condition was satisfied, the FMLS model was able to fit relatively well, with two free parameters, the prices observed in SPX option data from April 1999 to May 2000.In an effort to find models which were more flexible than the FMLS model and more parsimonious than the stochastic volatility models, researchers pursued the development of time-inhomogeneous additive processes.
In 2003 Nolder 7 proposed the modeling of asset returns with the H-self-similar additive processes of Sato.Carr et al. 8 constructed exponential H-self-similar additive models and showed that they could successfully price equity options of various underlyings across both maturity and strike.Madan and Eberlein 28 demonstrated that H-self-similar additive models could be useful for pricing options on annualized realized variance and annualized realized volatility.Such suitability is attributed to the dependence of the variance of realized quadratic variation to time t per unit time on the Hurst exponent.

Purpose and organization
The purpose of this paper is to compare the option pricing performance of certain additive models with respect to the choice of time-one distribution VG or NIG and with respect to the choice of additive model class Lévy, self-similar, and subordinated additive .We use SPX option data from the year 2005 to calibrate the exponential additive models corresponding to the VG, VGSSD, VHG, NIG, NIGSSD, and NHIG processes, along with the benchmark exponential VG-OU-Γ and NIG-OU-Γ Lévy stochastic volatility models.The latter are constructed by using the time-integral of a process of O-U type with associated stationary gamma distribution to time change an independent VG or NIG process.
The following describe the organization of this paper.In Section 2 the six additive models used in this study are defined.The moments of the increments corresponding to the variance gamma and normal inverse Gaussian families of processes are given in Section 3. Also included is a discussion of the dependence of increment variance on the Hurst exponent for both time-inhomogeneous additive classes.In Section 4 the characteristic functions of the exponential additive models are given, along with those of the Lévy stochastic volatility models.Algorithmic and data specifications for the calibration procedure are discussed in Section 5.The calibration results are presented in Section 6.The conclusion is presented in Section 7.

Construction of additive processes
We first begin with a few definitions and facts associated with self-similar additive processes.
ii For any choice of n ∈ N and 0 iii {X t } is stochastically continuous.
iv There exists a Ω 0 ∈ F with P Ω 0 1 such that, for every ω ∈ Ω 0 , X t ω is rightcontinuous in t ≥ 0 and has left limits in t > 0.
Remark 2.2.A Lévy process is a time-homogeneous additive process-an additive process which satisfies where μ z is the characteristic function of μ.
Since a self-decomposable distribution μ is infinitely divisible, it generates a Lévy process {X t : t ≥ 0}.In this case, the characteristic function of the distribution of X t is μ t z for all t ≥ 0 5, Lemma 7.10 ii , page 35 .The following theorem shows that a self-decomposable distribution also generates a self-similar additive process.Theorem 2.5 Sato 5, Theorem 16.1 ii , page 99 .If μ is a nontrivial self-decomposable distribution on R d , then, for any H > 0, there exists, uniquely in law, a nontrivial H-self-similar additive process {X t : t ≥ 0} such that P X 1 μ.Remark 2.6 Sato 5, Proof of Theorem 16.1, page 100 .In the previous theorem, the characteristic function of X t is given by The following representation theorem is useful for the purpose of checking selfdecomposability of a distribution.
where and k x is increasing on −∞, 0 and decreasing on 0, ∞ .
The following terminology will be used in the application of this theorem.Definition 2.8.A k-function of the self-decomposable distribution μ 5, page 403 is an Rvalued function on R \ {0} which is decreasing and left-continuous on R and which is increasing and right-continuous on R − with the Lévy density of μ given by k x /|x|, x ∈ R \ {0}.
The notation below will be used in the theorems and derivations which follow.Notation 2.9.Given an R-valued process X ≡ {X t : t ≥ 0}, we define the following: If X is R -valued, then we also have Laplace exponent: Remark 2.10.Given a Lévy process X and an independent, increasing additive process Z, the characteristic function of the time-t marginal of X time changed by Z can easily be computed by conditioning on Z t , yielding This computation is done in 2, page 109 , where the directing process Z is assumed to be time homogeneous.However, the derivation of 2.7 does not use the property of time homogeneity of Z.
In the remainder of this section, three sets of additive processes will be defined.For each set, one process will have a time-one marginal which is variance gamma distributed.The other process will have a time-one marginal which is normal inverse Gaussian distributed.

L évy processes
The three-parameter variance gamma process of Madan et al. 19 can be defined as a Brownian motion with drift which is time changed by an independent gamma process.Let σ, ν > 0 and μ ∈ R. On the probability space Ω, F, P let X ≡ {X t : t ≥ 0} be an R-valued Brownian motion, where X 1 is normal μ, σ 2 distributed.On the same probability space let Z ≡ {Z t : t ≥ 0} be an R -valued, independent gamma process where, for any t ≥ 0, Z t is gamma t/ν, 1/ν distributed.The variance gamma process, {VG t }, can be constructed as which characterizes the distributions of a variance gamma process.Equation 2.7 of Remark 2.10, Remark 2.2, and the Laplace exponent of Z t , are used to obtain the characteristic function of the random variable VG t as follows:

2.10
The normal inverse Gaussian process of Barndorff-Nielsen 6 can be defined as a Brownian motion with drift which is time changed by an independent inverse Gaussian process.Fix δ, γ > 0, and β ∈ R. On the probability space Ω, F, P let X ≡ {X t : t ≥ 0} be an R-valued Brownian motion, where X 1 is normal βδ 2 , δ 2 distributed.On the same probability space let Z ≡ {Z t : t ≥ 0} be an R -valued, independent inverse Gaussian process where, for any t ≥ 0, Z t is IG t, δγ distributed.The normal inverse Gaussian process, {NIG t }, can be constructed as Since the Laplace exponent of Z t is given by we have, similar to the previous case, the characteristic function of the distribution of NIG t : 2.13

Self-decomposability of time-one distributions
H-self-similar additive processes are constructed using Remark 2.6.Consequently, selfdecomposability of the one-dimensional marginal distribution at time one must be verified for the variance gamma and normal inverse Gaussian processes.

Variance gamma case
As mentioned in the construction of the variance gamma process in Section 2.1, X 1 is Normal μ, σ 2 distributed while Z t is gamma t/ν, 1/ν distributed for all t ≥ 0. Let Z 1 denote the Lévy density of the distribution of Z 1 see Table 2 .By the subordination theorem for Lévy processes 5, Theorem 30.1, page 197 , the Lévy density of the distribution of the time-one variance gamma random variable, VG 1 , is given by where the representation of the modified Bessel function of the second kind of order ν in 2.15 2, page 500 is used in the third step with ν 1/2 and α, β > 0.
Identity 2.16 29, page 109 is used in the final step.
The product given by x → |x| VG 1 x is positive, decreasing on R and increasing on R − if ν > 0. By Theorem 2.7, the k-function given by verifies self-decomposability of the time-one distribution of the variance gamma process, provided ν > 0. For an alternative proof of self-decomposability using the C, G, M parameterization of the variance gamma process, see CGMY 8 .

Normal inverse Gaussian case
Using the same construction as in Section 2.1 for the development of the NIG process, X 1 is normal βδ 2 , δ 2 distributed while Z t is IG t, δγ distributed.Denoting the Lévy density of the distribution of Z 1 by Z 1 Table 2 , the Lévy density of the distribution of the time-one normal inverse Gaussian random variable, NIG 1 , is given by

2.19
Since x → |x| NIG 1 x is decreasing on R and increasing on R − if and only if x → e β/ √ β 2 γ 2 x K 1 |x| is decreasing on R and increasing on R − , we make use of identity 2.20 30, equation 9.6.24: 2.20 We have, in a manner similar to Laforgia 31 ,

2.21
It follows that the product of NIG 1 and |x| is positive, decreasing on R and increasing on R − for all δ > 0, β ∈ R. By Theorem 2.7, the k-function given by verifies self-decomposability of the time-one distribution of the normal inverse Gaussian process.For an alternative proof of self-decomposability with this parameterization see Barndorff-Nielsen 6 .CGMY 8 address self-decomposability using a different parameterization.Sato's construction of H-self-similar additive processes from self-decomposable distributions Remark 2.6 may now be applied to the variance gamma and normal inverse Gaussian cases.

Construction of H-self-similar additive processes
Using the parameterizations for this study, we reproduce the results of CGMY 8 and Barndorff-Nielsen 6 NIG case only for the characteristic functions of the time-t marginal distributions of the self-similar additive analogs of the variance gamma and normal inverse Gaussian processes denoted in 8 as "VGSSD" and "NIGSSD," resp. .Let σ, ν, H > 0, μ ∈ R, and m P VG 1 be the self-decomposable variance gamma σ, μ, ν distribution on R. By Theorem 2.5 and Remark 2.6, the characteristic function of VGSSD t is given by E e iξ VGSSD t m t H ξ; σ, μ, ν

2.23
For the characteristic function of the distribution of NIGSSD t , let δ, γ, H > 0, β ∈ R, and m P NIG 1 be the self-decomposable normal inverse Gaussian δ, β, γ distribution on R.

Journal of Applied Mathematics and Decision Sciences
Similar to the previous case, the characteristic function of the distribution of NIGSSD t is given by

2.24
By construction, the distributions given by 2.10 and 2.23 agree at time one.The same is also true for 2.13 and 2.24 .

Subordinated additive processes
The subordinated, time-inhomogeneous additive processes in this study are a special case of the class of subordinated processes discussed by Madan and Yor 32 .Such processes are defined by time-changing a Brownian motion by an independent, increasing Markov process with independent and time-inhomogeneous increments.The first process to be constructed is a Brownian motion time changed by an independent H-self-similar additive gamma process VHG , while the second is a Brownian motion time changed by an independent H-selfsimilar additive inverse Gaussian process NHIG .In order to construct the self-similar additive directing processes, self-decomposability of the gamma and inverse Gaussian distributions must first be verified.The Lévy densities for these distributions are given in Table 2 where a, b > 0 for both distributions 4 .By Theorem 2.7, the k-functions given by the following equations: verify self-decomposability of the time-one distributions of the gamma and inverse Gaussian processes, respectively.The characteristic function of the distribution of VHG t is determined as follows.Let σ, ν, H > 0, and μ ∈ R. On the probability space Ω, F, P let X ≡ {X t : t ≥ 0} be an Rvalued Brownian motion where X 1 is normal μ, σ 2 distributed.On the same probability space, using Theorem 2.5, define Z ≡ {Z t : t ≥ 0} to be an independent, increasing H-selfsimilar additive process such that Z 1 is gamma 1/ν, 1/ν distributed.The subordinated, time-inhomogeneous additive analog of the VG process, {VHG t }, can be constructed as Using 2.7 in Remark 2.10, Remark 2.6, and the Laplace exponent of Z 1 , the characteristic function of the distribution of the random variable VHG t is given by Consequently, the distributions given by 2.10 , 2.23 , and 2.28 all agree at time one.The characteristic function of the distribution of NHIG t is similarly defined.Fix δ, γ, H > 0, and β ∈ R. On the probability space Ω, F, P define X ≡ {X t : t ≥ 0} to be an Rvalued Brownian motion where X 1 is normal βδ 2 , δ 2 distributed.On the same probability space, using Theorem 2.5, define Z ≡ {Z t : t ≥ 0} to be an independent, increasing Hself-similar additive process such that Z 1 is IG 1, δγ distributed.The subordinated, timeinhomogeneous additive analog of the NIG process, {NHIG t }, can be constructed as The Laplace exponent of Z 1 is given by Similar to the previous case, the characteristic function of the distribution of NHIG t is obtained as follows:

2.31
Note that 2.31 agrees at time one with 2.13 and 2.24 .

Increment moments
Given an R-valued additive process Y ≡ {Y t } on Ω, F, P , let μ Y t −Y s denote the distribution of Y t − Y s where 0 ≤ s < t.By additivity of Y, the characteristic function of the increment is given by  Moment NIG Table 5: Self-similar variance gamma process: increment moments.
Moment VGSSD Moment NIGSSD 5, Theorem 9.7, page 51 .Consequently, the moments of the increment are easily computed by taking derivatives of the difference of characteristic exponents.The mean, variance, skewness, and excess kurtosis of the increment Y t −Y s are calculated for each of the previously defined additive processes in Tables 3-8.An interesting property among the three classes of additive processes is the dependence of the variance of the increment Y t δ − Y t on time t for any t ≥ 0 and any fixed δ > 0 .The property of time homogeneity of a Lévy process implies that the variance of a Lévy increment is constant with respect to time.See Tables 3 and 4 for VG and NIG cases.Using Remark 2.6 and 3.1 , the variance of the time-t increment of any H-self-similar additive process {X t } is given by Var  Moment NHIG provided the variance of X 1 exists.Consequently, these processes have an increment variance which is constant resp., decreasing, increasing with respect to time when H 1/2 resp., H < 1/2, H > 1/2 .See Tables 5 and 6 for self-similar VG and NIG cases.Using 2.7 in Remark 2.10, Remark 2.6, and 3.1 , the variance of the time-t increment of a Brownian motion {X t } time changed by an independent, increasing H-self-similar additive process {Z t } is given by provided the variance of Z 1 exists.Consequently, this class of subordinated additive processes has an increment variance which is decreasing resp., increasing in time for H ≤ 1/2 resp., H ≥ 1 .For H ∈ 1/2, 1 the signum of the rate of change of increment variance with respect to time depends on the parameter values of the distributions of X 1 and Z 1 , as well as time t.See Tables 7 and 8 for the increment moments associated with the time-inhomogeneous, subordinated analogs of the VG and NIG processes.

Exponential additive models
Let X ≡ {X t } t∈ 0,T be an additive process on the stochastic basis Ω, F, F {F t }, P , where X is adapted to the filtration F. The underlying of a European-exercised option is modeled by the price process {S t } t∈ 0,T , which is defined, path by path, in 4.1 for each t ∈ 0, T , 4.1 where r and q are the continuously compounded values of both interest rate and dividend yield, respectively.Note that the property of independent increments of {X t } implies that {e − r−q t S t } is an F − P martingale.

Journal of Applied Mathematics and Decision Sciences
The characteristic functions of the time-t distributions of the logarithm of the underlying for the additive processes defined in Sections 2.1-2.3 are given below: i exponential variance gamma ii exponential H-self-similar variance gamma iii exponential variance H-gamma

Exponential L évy stochastic volatility models
We now define, according to CGMY 3 and Schoutens 4 , two Lévy stochastic volatility processes to be used for the creation of model benchmarks for the previously defined models.The solution of 4.9 is given by y t y 0 e −λt t 0 e −λ t−s dL s .

4.10
We fix a, b > 0 and specify {L t } to be a compound Poisson process with exponentially distributed jump sizes where CGMY 3 and Cont and Tankov 2, page 488 show that such a choice for {L t } yields the following expression for the characteristic function of the distribution of Z t : We now conclude the definition of Lévy stochastic volatility processes via subordination.Let {X t } be a variance gamma process or normal inverse Gaussian process which is independent of the increasing process {Z t }, as determined by 4.8 , 4.10 , and 4.11 .The VG-OU-Γ and NIG-OU-Γ processes can be constructed such that, for each t ≥ 0 and ω in the sample space VG-OU-Γ t ω X Z t ω ω , X t is a VG process, 4.13 NIG-OU-Γ t ω X Z t ω ω , X t is a NIG process.4.14

Journal of Applied Mathematics and Decision Sciences
Let {X Z t } t∈ 0,T be an R-valued Lévy stochastic volatility process on Ω, F, F, P as in 4.13 or 4.14 .The underlying price process at time t is defined, as in 3, 4 , to satisfy for each t ∈ 0, T , 4.15 where {e −t r−q S t } is a F − P martingale on another filtered probability space Ω, F, F, P .
The exponential pricing models corresponding to 4.13 and 4.14 are equivalent to the "VGSG" and "NIGSG" models of CGMY 3 when the "VGSG" and "NIGSG" correlation parameters are set to zero.Unlike the exponential additive pricing models, these do not prohibit dynamic arbitrage.Rather, the intent of definition 4.15 is to prohibit arbitrages which result from a limited form of dynamic trading, known as "static arbitrages" 3, 33 .CGMY 3 demonstrated that this simple as opposed to stochastic exponentiation of Lévy stochastic volatility processes yielded models which priced out-of-the-money SPX and equity options quite well, and better than their stochastically exponentiated counterparts.We refer to Cont and Tankov 2, pages 492, 493 for an alternative formulation in which the discounted price process is a martingale with respect to the enlarged filtration determined from the time-changed Lévy process and the directing process.The characteristic functions of the time-t marginal distributions are given below for the exponential VG-OU-Γ and NIG-OU-Γ models.

4.16
The characteristic function of the time-one marginal of the VG process is given by 4.17 The number of parameters in 4.13 is reduced from seven to six via where see 4 and Φ Z t is defined in 4.12 .

(ii) Exponential normal inverse Gaussian Ornstein-Uhlenbeck-gamma
The characteristic function of the time-one marginal of the NIG process is The number of parameters in 4.14 is reduced from seven to six via Φ NIG-OU-Γ ξ; β, δ, γ, λ, a, b, y 0 Φ NIG-OU-Γ ξ; β, δy 0 , γ, λ, a, by 0 , 1 4.22 see 4 .It follows that the characteristic function of the time-t distribution of the logarithm of the underlying is given by where see 4 and Φ Z t is defined in 4.12 .

Model implementation and data specifications
In order to calculate European option prices, the "modified call" method of Carr and Madan 34 is used.This requires an analytical expression for the inverse Fourier transform of a damped call price function in terms of the characteristic function of the distribution of the logarithm of the risk-neutral underlying.Put prices are calculated via the call price using put-call parity.Let Θ denote the set of model parameter vectors such that for all t, both θ ∈ Θ and t determine the time-t distribution of the logarithm of the risk-neutral underlying, ln S t .Given a set of N market option prices {P i } i 1•••N , the set of corresponding model option prices is denoted by In order to estimate θ, the quadratic pricing error, E : Θ → R , is minimized with where w i , T i , K i are the weighting factor, maturity, and strike price of option i, respectively.Following Cont and Tankov 2 , the weighting factor of each option is given by the reciprocal of the square of the Black-Scholes vega evaluated at the associated Black-Scholes market implied volatility.Consequently, the sum of squared differences between market and model implied volatilities is approximated by the sum of weighted squared differences between model and market prices.Cont and Tankov 35 discuss some of the computational difficulties pertaining to the calibration of a given model, for example, multiple local minima.
For further discussion on sources of error in calibration, we refer to 36 .
The option data used to calibrate the risk-neutral models in this study were 4 p.m. EST time-stamped bid and ask quotes of out-of-the-money Standard and Poor's 500 index options.The selected quote dates coincided with the second Tuesday of the month for the time period spanning January 2005 through December 2005.Market option prices were calculated as the average of bid and ask, while the closing index value was used as the underlying price.The interest rate and dividend yield for each quote date were calculated using the daily updated rate on the three-month Treasury bill 37 and the quarterly updated SPX dividend yield 38 , respectively.
Several filters were applied to the data.First, the nearest-term SPX options were excluded from this study, leaving options of five maturities ranging from approximately 1 month to 1 year.Second, any option with an ask or bid quote equal to zero was discarded.Third, options were excluded which violated the requirement of monotonicity of price with respect to strike.Most of these deep out-of-the-money options were eliminated, as in 8 , by requiring the ratio of option price to underlying price to be bounded below by a sufficiently large positive number.Here, an empirically determined value of 9.5 × 10 −4 was used to filter out such options.After all of the filters were applied, the moneyness strike/underlying values for all options in this study ranged from 0.5893 to 1.2680, with an average of 117 options per quote date.

Performance measures
In order to measure the pricing performance of each of the eight previously defined models, the average percentage error APE is now defined.Denote by N the number of options used on a particular quote date, {P i } the set of observed bid-ask averages, and {π i } the set of corresponding calculated model prices.The average percentage error is defined in 4 as 6.1 Since the market option prices are taken as the set of averages of the bid and ask quotes on a given quote date, each model price may deviate from the market price by at most ask − bid /2 if it is to lie between the bid and ask quotes.Setting the price difference to be ask − bid /2 and the market price to be ask bid /2 in 6.1 yields the measure of uncertainty in market prices given below: We define the in-sample APE using 6.1 with bid-ask quotes and parameter estimates corresponding to the same quote date.The one-day-ahead resp., one-week-ahead out-ofsample APE is defined using 6.1 with in-sample parameter estimates and bid-ask quotes obtained one day resp., one week later than the in-sample quote date.The in-sample resp., one-day-ahead out-of-sample, one-week-ahead out-of-sample market APE is calculated using 6.2 with bid-ask quotes obtained zero days resp., one day, one week later than the in-sample quote date.

Individual model performance
In Figures 1, 2, 3 are the bar graphs of the in-sample, one-day-ahead, and one-week-ahead out-of-sample APE for each quote date and model.The time averages of APE, taken over all twelve-quote dates, are reported in Table 9.The worst performing model was the VG Lévy with 11.9% resp., 12.7%, 14.6% mean in-sample resp., one-day-ahead, one-week-ahead APE.The best performing model was the self-similar NIG NIGSSD with respective mean APEs of 4.06%, 5.13%, and 8.32%.Furthermore, this model was the only one with a mean in-sample pricing error below that of the market.
Table 10 contains values of the mean ratio of APE, taken over all twelve-quote dates, for a model of one class to that of another with a common associated distribution VG or NIG .As shown in the sixth row of Table 10, the NIGSSD model outperformed its model benchmark, with APE values which were, on average, 84% and 86% of those obtained by the NIG-OU-Γ model for in-sample and one-day-ahead out-of-sample cases, respectively.In contrast, the VGSSD model, on average, had pricing errors which were 115% and 107% of those obtained by the VG-OU-Γ model for in-sample and one-day-ahead out-of-sample cases, respectively.Below is a summary of the NIGSSD pricing performance statistics taken from Figures 1-3.
i The proportions of quote dates in which the NIGSSD model had a lower APE than that of the market for the in-sample, one-day-ahead, and one-week-ahead out-of-sample cases were 9/12, 6/12, and 2/12, respectively.
ii The proportions of quote dates in which the NIGSSD model had a lower APE than those of both stochastic volatility models for the in-sample, one-day-ahead, and one-weekahead out-of-sample cases were 11/12, 12/12, and 8/12, respectively.

Performance comparison of additive classes
Rows one and two of Table 10 indicate the benefits of time inhomogeneity, in that the insample pricing error of the subordinated additive class was approximately 60% of that corresponding to the Lévy class.The benefit of self-similarity in the log-return process over self-similarity of the directing process is observed via rows three and four of Table 10.For the in-sample case, the NIGSSD model mean ratio of 0.70 performed considerably better against its subordinated additive counterpart than did the VGSSD model mean ratio of 0.83 .Such differences in mean error ratio decreased with time lag, such that, regardless of time-one distribution, the one-week-ahead values were approximately 73% and 85% for the subordinated additive : Lévy and self-similar : subordinated additive comparisons, respectively.

Performance comparison of NIG and VG models
In each row of Table 11 a model class is fixed, and the mean ratio taken over all twelvequote dates of APE of the NIG model to that of the corresponding VG model is given.The mean ratio of in-sample APE of the NIG resp., NHIG, NIGSSD model to that of the VG resp., VHG, VGSSD model was approximately 86% resp., 80%, 67% .Furthermore, the variation in mean APE ratio across additive model class decreased with time lag.As shown in rows one through three of the third column of Table 11, the one-week-ahead out-of-sample mean ratios were approximately 87%, regardless of the additive model class under consideration.Finally, the consistency with which the NIG distributions outperformed those of the VG family may be observed in the following summary statistics taken from Figures 1-3.The proportions of quote dates in which each of the four models of the NIG family had a lower APE than that of the corresponding member of the VG family for insample, one-day-ahead, and one-week-ahead out-of-sample cases were 12/12, 11/12, and 8/12, respectively.

Model price behavior across strikes and maturities
In Figures 4-7 the market and model prices are plotted versus moneyness strike/underlying for the set of options available on April 12, 2005.Results for this month are displayed in order to provide a typical observation of option pricing errors for the self-similar NIG model.Figures 4 and 5 consist of the plots for the VG family of models, while Figures 6 and 7 consist of the plots for the NIG family.We define the relative model price difference for two models corresponding to a given option as the absolute value of the difference between model prices divided by the market option price.The relative model price differences between the Lévy and subordinated additive models were the greatest for deep out-of-the-money puts belonging to the two nearest term maturities see Figures 4 and 6 .The relative model price differences between the subordinated additive and self-similar additive models were the .This same observation is true for relative price differences between the self-similar and stochastic volatility models, as shown in Figures 5 and 7.

Term structure of moments
The first four moments of the risk-neutral distributions determined from the options on the April 12, 2005 quote date are plotted versus maturity in Figure 8.For the set of parameters chosen by the market in the subordinated additive models, the skewness resp., kurtosis  is decreasing resp., increasing with maturity.This is contrary to the case of the Lévy models in which the skewness resp., kurtosis is increasing resp., decreasing with maturity.Self-similar processes have one-dimensional marginal distributions with a rather striking constant skewness and kurtosis.Invariance with respect to maturity is due to the fact that the characteristic function at time t is obtained by composition of the characteristic function  at time one with the map ξ → t H ξ. When the chain rule is applied to the composition, a factor of t H is produced.Consequently, the moment number is matched by the same number of factors of t H , thereby allowing cancelation of maturity in the calculation of skewness and excess kurtosis.

Hurst exponent of additive models
Tables 12 and 13 contain the time-averaged parameter estimates, taken over the twelve-quote dates, for each additive model.The sample means of H were approximately 0.6 and 0.9 for the self-similar and subordinated additive classes, respectively.Furthermore, the relative difference between sample means of H in the two self-similar models was a mere 0.091%, compared with 8.9% in the subordinated additive case.In the next section, we will see that the modest variation of Hurst exponent in self-similar models during the year 2005 is not maintained over the course of the first half of the following year.

Conclusion
In this paper, we considered additive processes defined via subordination, Sato's construction of H-self-similar additive processes from self-decomposable distributions, and a combination of the two methods.For each of these three methods, or model classes, two processes were defined: one with a time-one marginal variance gamma VG distribution, and the other, a time-one marginal normal inverse Gaussian NIG distribution.The corresponding models for the price process were calibrated, along with those constructed from two benchmark Lévy stochastic volatility processes, using out-of-the-money Standard and Poor's 500 index options from the year 2005.The four-parameter subordinated additive models had in-sample pricing errors which were approximately 60% pooled across time and time-one distribution of those obtained by the corresponding three-parameter Lévy models.The four-parameter self-similar additive models had in-sample pricing errors which were approximately 77% pooled across time and time-one distribution of those obtained by the corresponding subordinated additive models.The results for the comparisons between self-similar and stochastic volatility models were mixed.The in-sample pricing errors of the self-similar VG and NIG models were 115% and 84% of those obtained by the VG and NIG Lévy stochastic volatility models, respectively.For a given model class, the normal inverse Gaussian model consistently outperformed its variance gamma counterpart for in-sample and one-day-ahead out-of-sample cases.The Hurst exponent estimates associated with the self-similar and subordinated classes were approximately 0.6 and 0.9, respectively.Although the estimates of Hurst exponent of self-similar models had a modest range during the year 2005, Hurst exponent estimates can vary greatly over the course of a year, as demonstrated in Section 6.8.

Figure 1 :
Figure 1: Average percentage error in-sample, out-of-the-money options .

Figure 5 :
Figure 5: Market and model option prices versus moneyness strike/underlying for VGSSD and VG-OUgamma models: April 2005 quote date.

Figure 6 :Figure 7 :
Figure 6: Market and model option prices versus moneyness strike/underlying for NIG and NHIG models: April 2005 quote date.

Figure 8 :
Figure 8: Risk-neutral mean, variance, skewness, and kurtosis term structures for April 12, 2005 quote date VG models-top row NIG models-bottom row .

For d ≥ 1 a probability measure μ on R d is self-decomposable if, for all b > 1, there exists a probability measure ρ b on R d such that
1
These processes are formed by time changing a Lévy process by the time integral of a process of Ornstein-Uhlenbeck type with associated stationary gamma distribution.Let {Z t } It follows that the characteristic function of the time-t distribution of the logarithm of the underlying is given by , a, b, y 0 Φ VG-OU-Γ ξ; Cy 0 , G, M, λ, a, by 0 , 1 4.18 see 4 .

Table 9 :
Mean of average percentage error: twelve quote date time average January-December 2005 .Standard deviation is given in parentheses.

Table 10 :
Model class comparison: twelve quote date time-averaged ratio of average percentage errors January-December 2005 .Standard deviation is given in parentheses.

Table 11 :
Variance gamma and normal inverse Gaussian distribution comparison: twelve quote date timeaveraged ratio of average percentage errors January-December 2005 .Standard deviation is given in parentheses.
Market and model option prices versus moneyness strike/underlying for VG and VHG models: April 2005 quote date.