‚ 2 Ct (·)

= . (8.4)

‚K 2

Equating (8.3) and (8.4) across all states yields:

‚ 2 Ct (·)

= e’rt,„ „ ft— (ST )

‚K 2 K=ST

where rt,„ denotes the risk-free interest rate at time t with time to maturity „

and ft— (·) denotes the risk-neutral PDF or the SPD in t. Therefore, the SPD

is de¬ned as:

2

rt,„ „ ‚ Ct (·)

—

ft (ST ) = e . (8.5)

‚K 2 K=ST

This method constitutes a no-arbitrage approach to recover the SPD. No as-

sumption on the underlying asset dynamics are required. Preferences are not

restricted since the no-arbitrage method only assumes risk-neutrality with re-

spect to the underlying asset. The only requirements for this method are

that markets are perfect (i.e. no sales restrictions, transactions costs or taxes

and that agents are able to borrow at the risk-free interest rate) and that

C(·) is twice di¬erentiable. The same result can be obtained by di¬erentiat-

ing (8.1) twice with respect to K after setting for Z the call payo¬ function

Z(ST ) = (ST ’ K)+ .

8.2 Extracting the SPD using Call-Options 175

8.2.1 Black-Scholes SPD

The Black-Scholes call option pricing formula is due to Black and Scholes (1973)

and Merton (1973). In this model there are no assumptions regarding prefer-

ences, rather it relies on no-arbitrage conditions and assumes that the evolution

of the underlying asset price St follows a geometric Brownian motion de¬ned

through

dSt

= µdt + σdWt . (8.6)

St

Here µ denotes the drift and σ the volatility assumed to be constant.

The analytical formula for the price in t of a call option with a terminal date

T = t + „ , a strike price K, an underlying price St , a risk-free rate rt,„ , a

continuous dividend yield δt,„ , and a volatility σ, is:

∞

= e’rt,„ —

max(ST ’ K, 0)fBS,t (ST )dST

CBS (St , K, „, rt,„ , δt,„ ; σ)

0

’δt,„ „

¦(d1 ) ’ Ke’rt,„ „ ¦(d2 )

= St e

where ¦(·) is the standard normal cumulative distribution function and

1

log(St /K) + (rt,„ ’ δt,„ + 2 σ 2 )„

√

d1 = ,

σ„

√

= d1 ’ σ „ .

d2

As a consequence of the assumptions on the underlying asset price process the

Black-Scholes SPD is a log-normal density with mean (rt,„ ’ δt,„ ’ 1 σ 2 )„ and

2

variance σ 2 „ for log(ST /St ):

‚ 2 Ct

—

fBS,t (ST ) = ert,„ „

‚K 2 K=ST

[log(ST /St ) ’ (rt,„ ’ δt,„ ’ 1 σ 2 )„ ]2

1 2

√ exp ’

= .

2σ 2 „

2„

ST 2πσ

The risk measures Delta (∆) and Gamma (“) are de¬ned as:

‚CBS

def

∆BS = = ¦(d1 )

‚St

2

def ‚ CBS ¦(d1 )

√

“BS = =

2

‚St St σ „

176 8 Estimating State-Price Densities with Nonparametric Regression

The Black-Scholes SPD can be calculated in XploRe using the following quant-

let:

bsspd = spdbs(K,s,r,div,sigma,tau)

estimates the Black-Scholes SPD

The arguments are the strike prices (K), underlying price (s), risk-free interest

rate (r), dividend yields (div), implied volatility of the option (sigma), and

the time to maturity (tau). The output consist of the Black-Scholes SPD

(bsspd.fbs), ∆ (bsspd.delta), and the “ (bsspd.gamma) of the call options.

Please note that spdbs can be applied to put options by using the Put-Call

parity.

However, it is widely known that the Black-Scholes call option formula is not

valid empirically. For more details, please refer to Chapter 6. Since the Black-

Scholes model contains empirical irregularities, its SPD will not be consistent

with the data. Consequently, some other techniques for estimating the SPD

without any assumptions on the underlying di¬usion process have been devel-

oped in the last years.

8.3 Semiparametric estimation of the SPD

8.3.1 Estimating the call pricing function

The use of nonparametric regression to recover the SPD was ¬rst investigated

by A¨±t-Sahalia and Lo (1998). They propose to use the Nadaraya-Watson esti-

mator to estimate the historical call prices Ct (·) as a function of the following

state variables (St , K, „, rt,„ , δt,„ ) . Kernel regressions are advocated because

there is no need to specify a functional form and the only required assumption

is that the function is smooth and di¬erentiable, H¨rdle (1990). When the re-

a

gressor dimension is 5, the estimator is inaccurate in practice. Hence, there is

a need to reduce the dimension or equivalently the number of regressors. One

method is to appeal to no-arbitrage arguments and collapse St , rt,„ and δt,„

into the forward price Ft = St e(rt,„ ’δt,„ )„ in order to express the call pricing

function as:

C(St , K, „, rt,„ , δt,„ ) = C(Ft,„ , K, „, rt,„ ). (8.7)

8.3 Semiparametric estimation of the SPD 177

An alternative speci¬cation assumes that the call option function is homoge-

neous of degree one in St and K (as in the Black-Scholes formula) so that:

C(St , K, „, rt,„ , δt,„ ) = KC(St /K, „, rt,„ , δt,„ ). (8.8)

Combining the assumptions of (8.7) and (8.8) the call pricing function can be

further reduced to a function of three variables ( FK , „, rt,„ ).

t,„

Another approach is to use a semiparametric speci¬cation based on the Black-

Scholes implied volatility. Here, the implied volatility σ is modelled as a non-

parametric function, σ(Ft,„ , K, „ ):

C(St , K, „, rt,„ , δt,„ ) = CBS (Ft,„ , K, „, rt,„ ; σ(Ft,„ , K, „ )). (8.9)

Empirically the implied volatility function mostly depends on two parameters:

the time to maturity „ and the moneyness M = K/Ft,„ . Almost equivalently,

˜ ˜

one can set M = St /K where St = St ’ D and D is the present value of

the dividends to be paid before the expiration. Actually, in the case of a

dividend yield δt , we have D = St (1 ’ e’δt ). If the dividends are discrete, then

Dti e’rt,„i where ti is the dividend payment date of the ith dividend

D=

ti ¤t+„

and „i is its maturity.

Therefore, the dimension of the implied volatility function can be reduced to

σ(K/Ft,„ , „ ). In this case the call option function is:

C(St , K, „, rt,„ , δt,„ ) = CBS (Ft,„ , K, „, rt,„ ; σ(K/Ft,„ , „ )). (8.10)

ˆ

‚ Ct (·)

ˆ ˆ

Once a smooth estimate of σ (·) is obtained, estimates of Ct (·), ∆t =

ˆ ‚St ,

ˆ ˆ

‚ 2 Ct (·) ‚ 2 Ct (·)

ˆ

ˆ and ft— = ert,„ „

“t = , can be calculated.

2 ‚K 2

‚St

8.3.2 Further dimension reduction

The previous section proposed a semiparametric estimator of the call pricing

function and the necessary steps to recover the SPD. In this section the di-

mension is reduced further using the suggestion of Rookley (1997). Rookley

178 8 Estimating State-Price Densities with Nonparametric Regression

uses intraday data for one maturity and estimates an implied volatility surface

where the dimension are the intraday time and the moneyness of the options.

Here, a slightly di¬erent method is used which relies on all settlement prices

of options of one trading day for di¬erent maturities to estimate the implied

volatility surface σ(K/Ft,„ , „ ). In the second step, these estimates are used for

a given time to maturity which may not necessarily correspond to the maturity

of a series of options. This method allows one to compare the SPD at di¬erent

dates because of the ¬xed maturity provided by the ¬rst step. This is interesting

if one wants to study the dynamics and the stability of these densities.

Fixing the maturity also allows us to eliminate „ from the speci¬cation of the

implied volatility function. In the following part, for convenience, the de¬nition

˜

of the moneyness is M = St /K and we denote by σ the implied volatility. The

notation ‚f (x‚xi n ) denotes the partial derivative of f with respect to xi and

1 ,...,x

df (x)

the total derivative of f with respect to x.

dx

Moreover, we use the following rescaled call option function:

Cit

cit = ,

˜

St

˜

St

Mit = .

Ki

where Cit is the price of the ith option at time t and Ki is its strike price.

The rescaled call option function can be expressed as:

e’r„ ¦(d2 )

= c(Mit ; σ(Mit )) = ¦(d1 ) ’

cit ,

Mit

log(Mit ) + rt + 1 σ(Mit )2 „

2

√

d1 = ,

σ(Mit ) „

√

= d1 ’ σ(Mit ) „ .

d2

The standard risk measures are then the following partial derivatives (for no-

tational convenience subscripts are dropped):

‚C ‚C ˜ ‚c ,

∆= = = c(M, σ(M )) + S

˜ ˜

‚S ‚S ‚S

‚2C ‚2C 2

‚∆ ‚c ˜‚ c.

“= = = =2 +S

˜ ˜ ˜

‚S 2

‚S ‚S2 ‚S2

‚S

8.3 Semiparametric estimation of the SPD 179

where

‚c dc ‚M dc 1

= = ,

˜ ˜

dM ‚ S dM K

‚S

2

d2 c

‚2c 1

= .

˜ dM 2 K

‚S2

The SPD is then the second derivative of the call option function with respect

to the strike price:

2 2

r„ ‚ C r„ ˜ ‚ c

—

f (·) = e =e S . (8.13)

‚K 2 ‚K 2

The conversion is needed because c(·) is being estimated not C(·). The analyt-

ical expression of (8.13) depends on:

2

d2 c

‚2c M dc M

= +2

‚K 2 dM 2 dM K 2

K

dc

The functional form of is:

dM

dc dd1 ¦ (d2 ) dd2 ¦(d2 )

’ e’r„ + e’r„

= ¦ (d1 ) , (8.14)

M2

dM dM M dM

d2 c

while is:

dM 2

2

d2 c d2 d1 dd1

’ d1

= ¦ (d1 )

dM 2 dM 2 dM

2

e’r„ ¦ (d2 ) d2 d2 2 dd2 dd2

’ ’ ’ d2

2

M dM M dM dM

2e’r„ ¦(d2 )

’ (8.15)

M3

180 8 Estimating State-Price Densities with Nonparametric Regression

The quantities in (8.14) and (8.15) are a function of the following ¬rst deriva-

tives:

dd1 ‚d1 ‚d1 ‚σ

= + ,

dM ‚M ‚σ ‚M

dd2 ‚d2 ‚d2 ‚σ

= + ,

dM ‚M ‚σ ‚M

‚d1 ‚d2 1

√,

= =

‚M ‚M Mσ „

√

‚d1 log(M ) + r„ „

√

’

= + ,

σ2 „

‚σ 2

√

‚d2 log(M ) + r„ „

√

’ ’

= .

σ2 „

‚σ 2

For the remainder of this chapter, we de¬ne:

V = σ(M ),

‚σ(M )

V = ,

‚M

‚ 2 σ(M )

V = . (8.16)

‚M 2

The quantities in (8.14) and (8.15) also depend on the following second deriva-

tive functions:

√

d2 d1 1 1 V „ log(M ) + r„

√ √

=’ ’

+ +V

dM 2 σ2 „

Mσ „ M σ 2

log(M ) + r„ 1

√ √,

’

+ V 2V (8.17)

3„ M σ2 „

σ

√

d2 d2 1 1 V „ log(M ) + r„

√ √

=’ ’V

+ +

dM 2 σ2 „

Mσ „ M σ 2

log(M ) + r„ 1

√ √.

’

+ V 2V (8.18)

σ3 „ M σ2 „

Local polynomial estimation is used to estimate the implied volatility smile and

its ¬rst two derivatives in (8.16). A brief explanation will be described now.

8.3 Semiparametric estimation of the SPD 181

8.3.3 Local Polynomial Estimation

Consider the following data generating process for the implied volatilities:

σ = g(M, „ ) + σ — (M, „ )µ,

where E(µ) = 0, Var(µ) = 1. M, „ and µ are independent and σ — (m0 , „0 ) is

the conditional variance of σ given M = m0 , „ = „0 . Assuming that all third

derivatives of g exist, one may perform a Taylor expansion for the function g

in a neighborhood of (m0 , „0 ):

1 ‚2g

‚g

(m ’ m0 )2

g(m, „ ) ≈ g(m0 , „0 ) (m ’ m0 ) +

+

2 ‚M 2

‚M m0 ,„0 m0 ,„0

2

‚g 1‚ g

(„ ’ „0 )2

(„ ’ „0 ) +

+

2 ‚„ 2

‚„ m0 ,„0 m0 ,„0

2

1 ‚g

(m ’ m0 )(„ ’ „0 ).

+ (8.19)

2 ‚M ‚„ m0 ,„0

This expansion suggests an approximation by local polynomial ¬tting, Fan

and Gijbels (1996). Hence, to estimate the implied volatility at the target

point (m0 , „0 ) from observations σj (j = 1, . . . , n), we minimize the following

expression:

n

σj ’ β0 + β1 (Mj ’ m0 ) + β2 (Mj ’ m0 )2 + β3 („j ’ „0 )

j=1

(8.20)

2

2

+β4 („j ’ „0 ) + β5 (Mj ’ m0 )(„j ’ „0 ) KhM ,h„ (Mj ’ m0 , „j ’ „0 )

where n is the number of observations (options), hM and h„ are the bandwidth

controlling the neighborhood in each directions and KhM ,h„ is the resulting

kernel function weighting all observation points. This kernel function may be

a product of two univariate kernel functions.

For convenience use the following matrix de¬nitions:

(M1 ’ m0 )2 („1 ’ „0 )2

M1 ’ m0 „1 ’ „0 (M1 ’ m0 )(„1 ’ „0 )

«

1

(M2 ’ m0 )2 („2 ’ „0 )2

M2 ’ m0 „2 ’ „0 (M2 ’ m0 )(„2 ’ „0 ) ·

¬1

X = ¬. ·,

¬ ·

. . . . .

. . . . . .

. . . . . .

(Mn ’ m0 )2 („n ’ „0 )2

Mn ’ m0 „n ’ „0 (Mn ’ m0 )(„n ’ „0 )

1

182 8 Estimating State-Price Densities with Nonparametric Regression

« «

σ1 β0

σ = . , β = . .

W = diag{KhM ,h„ (Mj ’ m0 , „j ’ „0 )}

. .

and

¬· ¬·

. .

σn β5

Hence, the weighted least squares problem (8.20) can be written as

min (σ ’ Xβ) W (σ ’ Xβ) . (8.21)

β

and the solution is given by

’1

ˆ

β = X WX X W σ. (8.22)

A nice feature of the local polynomial method is that it provides the estimated

implied volatility and its ¬rst two derivatives in one step. Indeed, one has from

(8.19) and (8.20):

‚g ˆ

= β1 ,

‚M m0 ,„0

‚2g ˆ

= 2β2 .

‚M 2 m0 ,„0

One of the concerns regarding this estimation method is the dependence on the

bandwidth which governs how much weight the kernel function should place

on an observed point for the estimation at a target point. Moreover, as the

call options are not always symmetrically and equally distributed around the

ATM point, the choice of the bandwidth is a key issue, especially for estimation

at the border of the implied volatility surface. The bandwidth can be chosen

global or locally dependent on (M, „ ). There are methods providing ”optimal”

bandwidths which rely on plug-in rules or on data-based selectors.

In the case of the volatility surface, it is vital to determine one bandwidth for the

maturity and one for the moneyness directions. An algorithm called Empirical-

Bias Bandwidth Selector (EBBS) for ¬nding local bandwidths is suggested by

Ruppert (1997) and Ruppert, Wand, Holst and H¨ssler (1997). The basic idea

o

of this method is to minimize the estimate of the local mean square error at

each target point, without relying on asymptotic result. The variance and the

bias term are in this algorithm estimated empirically.

8.4 An Example: Application to DAX data 183

Using the local polynomial estimations, the empirical SPD can be calculated

with the following quantlet:

lpspd = spdbl(m,sigma,sigma1,sigma2,s,r,tau)

estimates the semi-parametric SPD.

The arguments for this quantlet are the moneyness (m), V (sigma), V (sigma1),

V (sigma2), underlying price (s) corrected for future dividends, risk-free in-

terest rate (r), and the time to maturity (tau). The output consist of the local

polynomial SPD (lpspd.fstar), ∆ (lpspd.delta), and the “ (lpspd.gamma)

of the call-options.

8.4 An Example: Application to DAX data

This section describes how to estimate the Black-Scholes and local polynomial

SPD using options data on the German DAX index.

8.4.1 Data

The dataset was taken from the ¬nancial database MD*BASE located at CASE

(Center for Applied Statistics and Economics) at Humboldt-Universit¨t zu

a

Berlin. Since MD*BASE is a proprietary database, only a limited dataset

is provided for demonstration purposes.

This database is ¬lled with options and futures data provided by Eurex. Daily

series of 1, 3, 6 and 12 months DM-LIBOR rates taken from the T homson

F inancial Datastream serve as riskless interest rates. The DAX 30 futures

and options settlement data of January 1997 (21 trading days) were used in this

study. Daily settlement prices for each option contract are extracted along with

contract type, maturity and strike. For the futures, the daily settlement prices,

maturities and volumes are the relevant information. To compute the interest

rates corresponding to the option maturities a linear interpolation between the

available rates was used.

The DAX is a performance index which means that dividends are reinvested.

However, assuming no dividend yields when inverting the Black-Scholes for-

mula results in di¬erent volatilities for pairs of puts and calls contrary to the

184 8 Estimating State-Price Densities with Nonparametric Regression

no-arbitrage assumption contained in the Put-Call parity. This remark can

be explained by the fact that until January 2002 domestic investors have an

advantage as they may receive a portion or all of the dividend taxes back de-

pending on their tax status. Dividend tax means here the corporate income

tax for distributed gains from the gross dividend.

Since the dividends are rebated to domestic investors the DAX should fall by

an amount contained between 0 and these dividend taxes. Indeed, the value of

this fall depends on the level of these taxes which may be equal to zero and on

the weights of domestic and foreign investors trading the DAX. These dividend

taxes have the same e¬ects as ordinary dividends and should therefore be used

for computing the implied volatilities and the future price implicit in the Black

Scholes formula.

Hafner and Wallmeier (2001) suggest a method in order to get around this

problem which consists in computing dividends implied by the Put-Call parity.

Indeed, combining the futures pricing formula

Ft,„F = St ert,„F „F ’ Dt,„f

and the Put-Call parity

Ct ’ Pt = St ’ Dt,„O ’ Ke’rt,„o „o

we obtain:

Ct ’ Pt = Ft,„F e’rt,„F + Dt,„F ,„O ’ Ke’rt,„O „O (8.23)

where „O is the maturity of the options, „F is the maturity of the nearest

forward whose volume is positive and Dt,„F ,„O = Dt,„F ’ Dt,„O is the di¬erence

between the present values of the dividends.

Using (8.23), implied dividends were computed for each pair of put and call

with the same strike. Theoretically, for a given time to maturity there must

be only one value for these implied dividends. For each maturity the average

of these implied dividends was used to compute the corrected price. Using this

method implied volatilities are more reliable as the systematic “gap” between

put and call volatilities disappears. The only uncertainty at this stage is due

to the interpolated rates for the maturity „O .

The dataset consists of one ¬le XFGData9701 with 11 columns.

8.4 An Example: Application to DAX data 185

1 Day

2 Month

3 Year

4 Type of option (1 for calls, 0 for puts)

5 Time to maturity (in calendar days)

6 Strike prices

7 Option prices

8 Corrected spot price (implied dividends taken into account)

9 Risk-free interest rate

10 Implied volatility

11 Non-corrected spot price

The data can be read into XploRe by loading the quantlib finance and then

issuing the following command:

data=read("XFGData9701.dat")

Next extract all call options on January 3, 1997 with the paf command:

data=paf(data,(data[,1]==3)&&(data[,4]==1))

8.4.2 SPD, delta and gamma

This section provides an example using XploRe to calculate the semiparametric

SPD using DAX index options data. It is assumed that the quantlib finance

has been loaded.

XFGSPDonematurity.xpl plots the SPD of the series of options closest to

maturity. This ¬rst example only uses smoothing method in one dimension.

XFGSPDoneday.xpl calculates and plots the local polynomial SPD for Jan-

uary 10, 1997 for di¬erent times to maturity („ = 0.125, 0.25, 0.375). After

loading the data, the implied volatility is estimated using the volsurf quantlet,

while the ¬rst and second derivatives are estimated using lpderxest quantlet.

In this example the grid size is 0.01. The bandwidth is chosen arbitrarily at 0.15

and 0.125 for the moneyness and maturity directions respectively. The criteria

used is a visual inspection of the ¬rst and second derivatives to ensure that

they are continuous and smooth. Next the quantlet spdbl is used to calculate

the SPD which is ¬nally displayed in Figure 8.1.

186 8 Estimating State-Price Densities with Nonparametric Regression

This ¬gure shows the expected e¬ect of time to maturity on the SPD, which

is a loss of kurtosis. The x-axis represents the terminal prices ST . The local

polynomial SPD displays a negative skew compared to a theoretical Black-

Scholes SPD. The major reason for the di¬erence is the measure of implied

volatility. Using the local polynomial estimators one captures the e¬ect of the

“volatility smile” and its e¬ects on the higher moments such as skewness and

kurtosis. This result is similar to what A¨

±t-Sahalia and Lo (1998) and Rookley

(1997) found in their study.

Semi-parametric SPD: 10-01-1997

5 10 15 20 25

Density*E-4

0

2500 3000 3500

Stock price at expiry

Figure 8.1. Local Polynomial SPD for „ = 0.125 (blue,¬lled), „ = 0.25

(black,dashed) and „ = 0.375 (red,dotted).

XFGSPDoneday.xpl

Figure 8.2 and Figure 8.3 show Delta and Gamma for the full range of strikes

and for three di¬erent maturities. This method allows the user to get in one

step both greeks in one estimation for all strikes and maturities.

A natural question that may arise is how do the SPDs evolve over time. In

this section an illustrative example is used to show the dynamics of the SPD

over the month of January 1997. XFGSPDonemonth.xpl estimates and plots

the SPD for each trading day in January 1997. The x-axis is the moneyness,

y-axis is the trading day, and the z-axis is the SPD. Figure 8.4 shows the local

polynomial SPD for the three ¬rst weeks of January, 1997.

8.4 An Example: Application to DAX data 187

Semi-parametric Delta: 10-01-1997

1

0.5

Delta

0

2500 3000 3500

Strike prices

Figure 8.2. Local Polynomial Delta for „ = 0.125 (blue,¬lled), „ = 0.25

(black,dashed) and „ = 0.375 (red,dotted).

Semi-parametric Gamma: 10-01-1997

5 10 15 20 25

Gamma*E-4

0

2500 3000 3500

Strike prices

Figure 8.3. Local Polynomial Gamma for „ = 0.125 (blue,¬lled), „ =

0.25 (black,dashed) and „ = 0.375 (red,dotted).

XFGSPDoneday.xpl

8.4.3 Bootstrap con¬dence bands

Rookley™s method serves to estimate the SPD, where V , V and V from (8.16)

are computed via local polynomials. The method is now applied to estimate

a SPD whose maturity is equal to the maturity of a series of options. In this

case, the nonparametric regression is a univariate one.

188 8 Estimating State-Price Densities with Nonparametric Regression

Local-Polynomial SPD: 01-1997, tau=0.250

5.97

24

23

4.48

22

21

20

3.00

17

16

15

14

1.51

13

10

9

8

7

6 1.20

1.10

1.00

0.90

Figure 8.4. Three weeks State-Price Densities on a moneyness scale.

XFGSPDonemonth.xpl

With a polynomial of order p = 2 and a bandwidth h = n’1/9 , it can be

shown that

ˆ—

E|fn ’ f — |2 = O n’4/9 ,

because

ˆ = O n’8/9 ,

E|Vn ’ V |2

ˆ = O n’4/9 ,

E|Vn ’ V |2

ˆ = O n’4/9 .

E|Vn ’ V |2

8.4 An Example: Application to DAX data 189

This result can be obtained using some theorems related to local polynomial

estimation, for example in Fan and Gijbels (1996), if some boundary conditions

are satis¬ed.

ˆ— ˆ—

An asymptotic approximation of fn is complicated by the fact that fn is a

non linear function of V , V and V . Analytical con¬dence intervals can be

obtained using delta methods proposed by A¨ ±t-Sahalia (1996). However, an

alternative method is to use the bootstrap to construct con¬dence bands. The

idea for estimating the bootstrap bands is to approximate the distribution of

ˆ

sup |f — (k) ’ f — (k)|.

k

The following procedure illustrates how to construct bootstrap con¬dence

bands for local polynomial SPD estimation.

1. Collect daily option prices from MD*BASE, only choose those options

with the same expiration date, for example, those with time to maturity

49 days on Jan 3, 1997.

2. Use the local polynomial estimation method to obtain the empirical SPD.

Notice that when „ is ¬xed the forward price F is also ¬xed. So that the

implied volatility function σ(K/F ) can be considered as a ¬xed design

situation, where K is the strike price.

3. Obtain the con¬dence band using the wild bootstrap method. The wild

bootstrap method entails:

• Suppose that the regression model for the implied volatility function

σ(K/F ) is:

Ki

+ µi , i = 1, · · · , n.

Yi = σ

F

• Choose a bandwidth g which is larger than the optimal h in or-

der to have oversmoothing. Estimate the implied volatility function

σ(K/F ) nonparametrically and then calculate the residual errors:

Ki

µi = Yi ’ σh

˜ ˆ .

F

• Replicate B times the series of the {˜i } with wild bootstrap ob-

µ

—,j

taining {µi } for j = 1, · · · , B, H¨rdle (1990), and build B new

a

190 8 Estimating State-Price Densities with Nonparametric Regression

bootstrapped samples:

Ki

Yi—,j = σg + µ—,j .

ˆ i

F

• Estimate the SPD f —,j using bootstrap samples, Rookley™s method

and the bandwidth h, and build the statistics

ˆ

Tf = sup |f —,j (z) ’ f — (z)|.

—

z

ˆ ˆ

• Form the (1 ’ ±) bands [f — (z) ’ tf — ,1’± , f — (z) + tf — ,1’± ],

—

where tf — ,1’± denotes the empirical (1 ’ ±)-quantile of Tf .

Two SPDs (Jan 3 and Jan 31, 1997) whose times to maturity are 49 days

were estimated and are plotted in Figure (8.5). The bootstrap con¬dence

band corresponding to the ¬rst SPD (Jan 3) is also visible on the chart. In

Figure (8.6), the SPDs are displayed on a moneyness metric. It seems that the

di¬erences between the SPDs can be eliminated by switching to the moneyness

metric. Indeed, as can be extracted from Figure 8.6, both SPDs lie within

the 95 percent con¬dence bands. The number of bootstrap samples is set to

B = 100. The local polynomial estimation was done on standardized data, h

is then set to 0.75 for both plots and g is equal to 1.1 times h. Notice that

greater values of g are tried and the conclusion is that the con¬dence bands

are stable to an increase of g.

8.4.4 Comparison to Implied Binomial Trees

In Chapter 7, the Implied Binomial Trees (IBT) are discussed. This method is a

close approach to estimate the SPD. It also recovers the SPD nonparametrically

from market option prices and uses the Black Scholes formula to establish the

relationship between the option prices and implied volatilities as in Rookely™s

method. In Chapter 7, the Black Scholes formula is only used for Barle and

Cakici IBT procedure, but the CRR binomial tree method used by Derman

and Kani (1994) has no large di¬erence with it in nature. However, IBT and

nonparametric regression methods have some di¬erences caused by di¬erent

modelling strategies.

The IBT method might be less data-intensive than the nonparametric regres-

sion method. By construction, it only requires one cross section of prices. In the

8.4 An Example: Application to DAX data 191

SPDs and bootstrap CB, tau= 49 days

10 15 20 25 30

density*E-4

5

0

2400 2600 2800 3000

Stock price at expiry S(T)

Figure 8.5. SPD estimation and bootstrap con¬dence band.

XFGSPDcb.xpl

SPDs and bootstrap CB, tau= 49 days

8

6

density

4

2

0

0.8 0.9 1 1.1

S(T)/F(t,T)

Figure 8.6. SPD estimation and bootstrap con¬dence band (moneyness

metric).

XFGSPDcb2.xpl

earlier application with DAX data, option prices are used with di¬erent times

to maturity for one day to estimate the implied volatility surface ¬rst in order

192 8 Estimating State-Price Densities with Nonparametric Regression

SPD estimations: 19970103, tau= 77 days

0.2

0.15

density*E-2

0.1

0.05

0

2000 2500 3000 3500

stock price

Figure 8.7. Comparison of di¬erent SPD estimations, by Rookley™s

method (blue) and IBT (black, thin).

XFGSPDcom.xpl

to construct the tree using the relation formula between option prices and risk-

neutral probabilities. The precision of the SPD estimation using IBT is heavily

a¬ected by the quality of the implied volatility surface and the choice of the

levels of the implied tree. Furthermore, from the IBT method only risk-neutral

probabilities are obtained. They can be considered as a discrete estimation of

the SPD. However, the IBT method is not only useful for estimating SPD, but

also for giving a discrete approximation of the underlying process.

The greatest di¬erence between IBTs and nonparametric regression is the re-

quirement of smoothness. The precision of Rookley™s SPD estimation is highly

dependent on the selected bandwidth. Even if very limited option prices are

given, a part of the SPD estimation still can be obtained using nonparametric

regression, while the IBT construction has to be given up if no further struc-

ture is invoked on the volatility surface. Rookley™s method has on ¬rst sight

no obvious di¬erence with A¨ ±t-Sahalia™s method theoretically, A¨±t-Sahalia and

Lo (1998). But investigating the convergence rate of the SPD estimation using

A¨±t-Sahalia™s method allows one to conduct statistical inference such as test of

the stability of the SPD and tests of risk neutrality.

8.4 An Example: Application to DAX data 193

The quantlet XFGSPDcom.xpl shows a comparison of the SPD estimates by

IBT and Rookley™s methods. The di¬erences between these two SPD estimates

may be due to the selection of the bandwidths in Rookley™s method, the choice

of steps in the construction of the IBT and the use of DAX implied dividends

in Rookley™s method. Figure 8.7 shows the implied binomial trees and the local

polynomial SPDs for January, 3 1997.

Both densities seems to be quiet di¬erent. Indeed, the IBTs SPD shows a fatter

left tail than the Rookley™s one and the Rookley™s SPD shows a larger kurtosis.

To test which of both densities is more reliable, a cross-validation procedure is

performed. The idea of this test is to compare the theoretical prices based on

(8.1) with those observed on the market. However, as the whole tails are not

available for the Rookley™s SPD, the test is done on butter¬‚y spreads de¬ned in

Section 8.2 since their prices should not be in¬‚uenced by the tails of the SPDs.

For cross-validation, we remove the three calls used to calculate the observed

butter¬‚y prices from the sample before estimating the SPD. Moreover, since

the largest di¬erence between both SPDs is observed at the ATM point (see

Figure 8.7), the test is applied only to the two butter¬‚y spreads whose centers

surround the ATM point. The width 2∆K of the butter¬‚y spread is set to 200.

This procedure is done for the 21 days of January 1997. Figure 8.8 displays

the results in term of relative pricing error E:

Pobserved ’ PSP D

E=

Pobserved

where Pobserved is the observed price of the butter¬‚y spread and PSP D is the

price computed using the SPD estimate and (8.1). It seems that both SPDs

have a too small kurtosis since the observed prices of butter¬‚y spreads are

larger than those of both SPDs in most of the cases. However, Rookley™s SPD

is in mean nearer to the observed price than the IBT™s one.

Bibliography

A¨

±t-Sahalia, Y. (1996). The Delta method for Nonparametric Kernel Function-

als, mimeo.

A¨

±t-Sahalia, Y. and Lo, A. W. (1998). Nonparametric estimation of state-price

densities implicit in ¬nancial asset prices, Journal of Finance 53: 499“547.

194 8 Estimating State-Price Densities with Nonparametric Regression

Pricing errors for butterfly spread

20

Pricing error in %

10

0

-20 -10

5 10 15 20 25 30

Date

Pricing errors for butterfly spread

20

Pricing error in %

10

0

-10

5 10 15 20 25 30

Date

Figure 8.8. The upper graph display the relative pricing errors for the

butter¬‚y spread centered on the nearest strike on the left side of the

ATM point. The second graph corresponds to the butter¬‚y spread

centered on the nearest strike on the right side of the ATM point. The

black lines represent the IBT™s pricing errors and the blue the Rookley™s

errors.

Arrow, K. (1964). The role of securities in the optimal allocation of risk bearing,

Review of Economic Studies 31: 91“96.

Bahra, B. (1997). Implied risk-neutral probability density functions from option

8.4 An Example: Application to DAX data 195

prices: theory and application. Bank of England Working Paper 66.

Black, F. and Scholes, M. (1973). The pricing of options and corporate liabili-

ties, Journal of Political Economy 81: 637“654.

Breeden, D. and Litzenberger, R. H. (1978). Prices of state-contingent claims

implicit in option prices, Journal of Business 51: 621“651.

Debreu, G. (1959). Theory of Value, John Wiley and Sons, New York.

Derman, E. and Kani, I. (1994). Riding on the smile, Risk 7: 32“39.

Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Apllications,

Chapman and Hall, New York. Vol. 66 of Monographs on Statistics and

Applied Probability.

Hafner, R. and Wallmeier, M. (2001). The dynamics of DAX implied volatili-

ties, Quarterly International Journal of Finance 1: 1“27.

H¨rdle, W. (1990). Applied Nonparametric Regression, Cambridge University

a

Press, New York.

H¨rdle, W., Hl´vka, Z. and Klinke, S. (2000). XploRe Application Guide,

a a

Springer-Verlag, Berlin.

Hutchinson, J., Lo, A. and Poggio, A. (1994). A nonparametric approach to the

pricing and hedging of derivative securities via learning networks, Journal

of Finance 49: 851“889.

Lucas, R. E. (1978). Asset prices in an exchange economy, Econometrica

46(1429-1446).

Melick, W. and Thomas, C. (1997). Recovering an Asset™s Implied PDF from

Option Prices: An application to Crude Oil During the Gulf Crisis, Jour-

nal of Financial and Quantitative Analysis 32: 91“115.

Merton, R. B. (1973). Rational theory of option pricing, Bell Journal of Eco-

nomics and Management Science 4: 141“183.

Rookley, C. (1997). Fully exploiting the information content of intra-day option

quotes: Applications in option pricing and risk management. mimeo.

Rubinstein, M. (1976). The valuation of uncertain income streams and the

pricing of options, Bell Journal of Economics 7(407-425).

196 8 Estimating State-Price Densities with Nonparametric Regression

Rubinstein, M. (1985). Nonparametric tests of alternative option pricing mod-

els using all reported trades and quotes on the 30 most active cboe op-

tion classes from august 23, 1976 to august 31, 1978, Journal of Finance

40: 455“480.

Rubinstein, M. (1994). Implied binomial trees, Journal of Finance 49: 771“818.

Ruppert, D. (1997). Empirical-bias bandwidths for local polynomial nonpara-

metric regression and density estimation, Journal of the American Statis-

tical Association 92: 1049“1062.

Ruppert, D., Wand, M. P., Holst, U. and H¨ssler, O. (1997). Local polynomial

o

variance-function estimation, Technometrics 39: 262“273.

Yatchew, A. and H¨rdle, W. (2002). Dynamic nonparametric state price density

a

estimation using constrained least squares and the bootstrap. Journal of

Econometrics, forthcoming.

9 Trading on Deviations of Implied

and Historical Densities

Oliver Jim Blaskowitz and Peter Schmidt

9.1 Introduction

In recent years a number of methods have been developed to infer implied

state price densities (SPD) from cross sectional option prices, Chapter 7 and

8. Instead of comparing this density to a historical density extracted from the

observed time series of the underlying asset prices, i.e. a risk neutral density to

an actual density, Ait“Sahalia, Wang and Yared (2000) propose to compare two

risk neutral densities, one obtained from cross sectional S&P 500 option data

and the other from the S&P 500 index time series. Furthermore, they propose

trading strategies designed to exploit di¬erences in skewness and kurtosis of

both densities. The goal of this article is to apply the procedure to the german

DAX index. While the option implied SPD is estimated by means of the Barle

and Cakici, Barle and Cakici (1998), implied binomial tree version, the time

series density is inferred from the time series of the DAX index by applying a

method used by Ait“Sahalia, Wang and Yared (2000). Based on the comparison

of both SPDs the performance of skewness and kurtosis trades is investigated.

We use options data included in MD*BASE. This is a database located at

CASE (Center for Applied Statistics and Economics) of Humboldt“Universit¨ta

zu Berlin. The time period is limited to data of the period between 01/01/97

and 12/31/99 for which MD*BASE contains daily closing prices of the DAX

index, EUREX DAX option settlement prices and annual interest rates which

are adjusted to the time to maturity of the above mentioned EUREX DAX

options.

While Section 9.2 applies the Barle and Cakici implied binomial tree algorithm

198 9 Trading on Deviations of Implied and Historical Densities

which estimates the option implied SPD using a two week cross section of DAX

index options, Section 9.3 explains and applies the method to estimate DAX

time series SPD from 3 months of historical index prices. Following, in Section

9.4 we compare the conditional skewness and kurtosis of both densities. Section

9.5 and 9.6 complete the chapter with the investigation of 4 trading strategies

and Section 9.7 completes with some critical remarks.

9.2 Estimation of the Option Implied SPD

Barle“Cakici™s modi¬cation of Derman“Kani™s Implied Binomial Tree (IBT)

yields a proxy for the option implied SPD, f — , see Chapter 7. XploRe provides

quantlets computing Derman“Kani™s and Barle“Cakici™s IBT™s. Since the latter

proved to be slightly more robust than the former, Jackwerth (1999), we decide

to use Barle“Cakici™s IBT to compute the option implied SPD. In the following

subsection, we follow closely the notation used in Chapter 7. That is, N denotes

the number of evenly spaced time steps of length ∆t in which the tree is divided

into (so we have N + 1 levels). Fn,i = er∆t sn,i is the forward price of the

underlying, sn,i , at node i at, level n. Each level n corresponds to time tn =

n∆t.

9.2.1 Application to DAX Data

Using the DAX index data from MD*BASE, we estimate the 3 month option

implied IBT SPD f — by means of the XploRe quantlets IBTbc and volsurf and

a two week cross section of DAX index option prices for 30 periods beginning

in April 1997 and ending in September 1999. We measure time to maturity

(TTM) in days and annualize it using the factor 360, giving the annualized

time to maturity „ = TTM/360. For each period, we assume a ¬‚at yield curve.

We extract from MD*BASE the maturity consistent interest rate.

We describe the procedure in more detail for the ¬rst period. First of all, we

estimate the implied volatility surface given the two week cross section of DAX

option data and utilizing the XploRe quantlet volsurf which computes the 3

dimensional implied volatility surface (implied volatility over time to maturity

and moneyness) using a kernel smoothing procedure. Friday, April 18, 1997

is the 3rd Friday of April 1997. On Monday, April 21, 1997, we estimate the

volatility surface, using two weeks of option data from Monday, April 7, 1997,

to Friday, April 18, 1997. Following, we start the IBT computation using the

9.2 Estimation of the Option Implied SPD 199

DAX price of this Monday, April 21, 1997. The volatility surface is estimated

for the moneyness interval [0.8, 1.2] and the time to maturity interval [0.0, 1.0].

Following, the XploRe quantlet IBTbc takes the volatility surface as input and

computes the IBT using Barle and Cakici™s method. Note that the observed

smile enters the IBT via the analytical Black“Scholes pricing formula for a call

C(Fn,i , tn+1 ) and for a put P (Fn,i , tn+1 ) which are functions of St1 = s1,1 ,

K = Fn,i , r, tn+1 and σimpl (Fn,i , tn+1 ). We note, it may happen that at the

edge of the tree option prices, with associated strike prices Fn,i and node prices

sn+1,i+1 , have to be computed for which the moneyness ratio sn+1,i+1 /Fn,i is

outside the intverall [0.8, 1.2] on which the volatility surface has been estimated.

In these cases, we use the volatility at the edge of the surface. Note, as well,

that the mean of the IBT SPD is equal to the futures price by construction of

the IBT.

Finally, we transform the SPD over sN +1,i into a SPD over log“returns uN +1,i =

ln(sN +1,i /s1,1 ) as follows:

sN +1,i x

P(sN +1,i = x) = P ln = ln = P uN +1,i = u

s1,1 s1,1

where u = ln(x/s1,1 ). That is, sN +1,i has the same probability as uN +1,i . See

Figure 9.1 for the SPD computed with parameters N = 10 time steps and

interest rate r = 3.23.

A crucial aspect using binomial trees is the choice of the number of time steps

N in which the time interval [t, T ] is divided. In general one can state, the more

time steps are used the better is the discrete approximation of the continuous

di¬usion process and of the SPD. Unfortunately, the bigger N , the more node

prices sn,i possibly have to be overridden in the IBT framework. Thereby we are

e¬ectively losing the information about the smile at the corresponding nodes.

Therefore, we computed IBT™s for di¬erent numbers of time steps. We found

no hint for convergence of the variables of interest, skewness and kurtosis. Since

both variables seemed to ¬‚uctuate around a mean, we compute IBT™s with time

steps 10, 20, . . . , 100 and consider the average of these ten values for skewness

and kurtosis as the option implied SPD skewness and kurtosis.

Applying this procedure for all 30 periods, beginning in April 1997 and ending

in September 1999, we calculate the time series of skewness and kurtosis of the

3 month implied SPD f — shown in Figures 9.3 and 9.4. We see that the implied

SPD is clearly negatively skewed for all periods but one. In September 1999 it

is slightly positively skewed. The pattern is similar for the kurtosis of f — which

is leptokurtic in all but one period. In October 1998 the density is platykurtic.

200 9 Trading on Deviations of Implied and Historical Densities

IBT SPD

0.2

0.15

SPD

0.1

0.05

0

-0.3 -0.2 -0.1 0 0.1 0.2

LogReturn

Figure 9.1. Option implied SPD estimated on April 21, 1997, by an

IBT with N = 10 time steps, S0 = 3328.41, r = 3.23 and „ = 88/360.

9.3 Estimation of the Historical SPD

While the previous section was dedicated to ¬nding a proxy for f — used by

investors to price options, this section approximates the historical underlyings™

density g — for date t = T using all the information available at date t = 0. Of

course, if the process governing the underlying asset dynamics were common

knowledge and if agents had perfect foresight, then by no arbitrage arguments

both SPDs should be equal. Following Ait“Sahalia, Wang and Yared (2000),

the density extracted from the observed underlyings™ data is not comparable to

the density implied by observed option data without assumptions on investor™s

preferences. As in H¨rdle and Tsybakov (1995), they apply an estimation

a

method which uses the observed asset prices to infer indirectly the time series

SPD. First, we will explain the estimation method for the underlyings™ SPD.

Second, we apply it to DAX data.

9.3 Estimation of the Historical SPD 201

9.3.1 The Estimation Method

Assuming the underlying S to follow an ˆ di¬usion process driven by a Brow-

Ito

nian motion W :

dSt = µ(St )dt + σ(St )dWt . (9.1)

Ait“Sahalia, Wang and Yared (2000) rely on Girsanov™s characterization of the

change of measure from the actual density to the SPD. It says the di¬usion

function of the asset™s dynamics is identical under both the risk neutral and

the actual measure and only the drift function needs to be adjusted, leading to

the risk neutral asset dynamics:

—

(rt,„ ’ δt,„ )St dt + σ(St )dWt— .

— —

dSt = (9.2)

Let gt (St , ST , „, rt,„ , δt,„ ) denote the conditional density of ST given St gen-

—

erated by the dynamics de¬ned in equation (9.1) and gt (St , ST , „, rt,„ , δt,„ )

denote the conditional density generated by equation (9.2) then f — can only be

compared to the risk neutral density g — and not to g.

A crucial feature of this method is that the di¬usion functions are identical

under both the actual and the risk neutral dynamics (which follows from Gir-

sanov™s theorem). Therefore, it is not necessary to observe the risk neutral path

of the DAX index {St }. The function σ(•) is estimated using N — observed

—

index values {St } and applying Florens“Zmirou™s (1993) (FZ) nonparametric

version of the minimum contrast estimators:

N — ’1 ’S

KF Z ( SiF Z )N — {S(i+1)/N — ’ Si/N — }2

i=1 h

σF Z (S)

ˆ = , (9.3)

N— ’S

KF Z ( SiF Z )

i=1 h

where KF Z (•) is a kernel function and hF Z is a bandwidth parameter such

that:

(N — hF Z )’1 ln(N — ) N — h4 Z

’ ’

0 and 0

F

as N — ’ ∞. Without imposing restrictions on the drift function σF Z is an

ˆ

unbiased estimator of σ in the model speci¬ed in equation (9.2). Since the DAX

index is a performance index (δt,„ = 0), the risk neutral drift rate of equation

(9.2) is equal to rt,„ .

Once σ(•) is estimated, the time series SPD g — can be computed by Monte

Carlo integration. Applying the Milstein scheme (Kloeden, Platen and Schurz

202 9 Trading on Deviations of Implied and Historical Densities

(1994)), we simulate M = 10, 000 paths of the di¬usion process:

—

= rt,„ St dt + σF Z (St )dWt—

— —

dSt ˆ (9.4)

for a time period of 3 months, starting value St=0 equal to the DAX index value

at the beginning of the period, collect the endpoints at T of these simulated

paths {ST,m : m = 1, . . . , M } and annualize the index log“returns. Then

g — is obtained by means of a nonparametric kernel density estimation of the

continuously compounded log“returns u:

M

um ’ u

1

—

pt (u)

ˆ = KM C (9.5)

M hM C hM C

m=1

where um is the log“return at the end of the mth path and KM C (•) is a kernel

function and hM C is a bandwidth parameter. The equation:

log(S/St )

p— (u)du

P(ST ¤ S) P(u ¤ log(S/St ))

= = t

’∞

with u = ln(ST /St ) relates this density estimator to the SPD g — :

p— (log(S/St ))

— ‚

¤ S)

gt (S) = ‚S P(ST = .

t

S

√

This method results in a nonparametric estimator g — which is N — “consistent

ˆ

as M ’ ∞ even though σF Z converges at a slower rate (Ait“Sahalia, Wang

ˆ

and Yared (2000)).

In the absence of arbitrage, the futures price is the expected future value of the

spot price under the risk neutral measure. Therefore the time series distribu-

tion is translated such that its mean matches the implied future price. Then the

bandwidth hM C is chosen to best match the variance of the IBT implied distri-

bution. In order to avoid over“ or undersmoothing of g — , hM C is constrained

to be within 0.5 to 5 times the optimal bandwidth implied by Silverman™s rule

of thumb. This procedure allows us to focus the density comparison on the

skewness and kurtosis of the two densities.

9.3.2 Application to DAX Data

Using the DAX index data from MD*BASE we estimate the di¬usion function

σ 2 (•) from equation (9.2) by means of past index prices and simulate (forward)

M = 10, 000 paths to obtain the time series density, g — .

9.3 Estimation of the Historical SPD 203

Time Series Density

3

2

Density

1

-0.3 -0.2 -0.1 0 0.1 0.2 0.3

LogReturns

Figure 9.2. Mean and variance adjusted estimated time series density

on Friday, April 18, 1997. Simulated with M = 10, 000 paths, S0 =

3328.41, r = 3.23 and „ = 88/360.

To be more precise, we explain the methodology for the ¬rst period in more

detail. First, note that Friday, April 18, 1997, is the 3rd Friday of April 1997.

Thus, on Monday, April 21, 1997, we use 3 months of DAX index prices from

Monday, January 20, 1997, to Friday, April 18, 1997, to estimate σ 2 . Following,

on the same Monday, we start the 3 months ˜forward™ Monte Carlo simulation.

The bandwidth hF Z is determined by Cross Validation applying the XploRe

quantlet regxbwcrit which determines the optimal bandwidth from a range

of bandwidths by using the resubstitution estimator with the penalty function

™Generalized Cross Validation™.

Knowing the di¬usion function it is now possible to Monte Carlo simulate the

index evolution. The Milstein scheme applied to equation (9.2) is given by:

Si/N —— = S(i’1)/N —— + rS(i’1)/N —— ∆t + σ(S(i’1)/N —— )∆Wi/N —— +

1 ‚σ