. 7
( 14)


‚ 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:
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
= µdt + σdWt . (8.6)

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,„ ; σ)
’δt,„ „
¦(d1 ) ’ Ke’rt,„ „ ¦(d2 )
= St e
where ¦(·) is the standard normal cumulative distribution function and
log(St /K) + (rt,„ ’ δt,„ + 2 σ 2 )„

d1 = ,

= d1 ’ σ „ .

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
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 „
ST 2πσ
The risk measures Delta (∆) and Gamma (“) are de¬ned as:
∆BS = = ¦(d1 )
def ‚ CBS ¦(d1 )

“BS = =
‚St St σ „
176 8 Estimating State-Price Densities with Nonparametric Regression

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

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
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-
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,„ ).

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
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

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.

Moreover, we use the following rescaled call option function:
cit = ,
Mit = .

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 ,
log(Mit ) + rt + 1 σ(Mit )2 „

d1 = ,
σ(Mit ) „

= d1 ’ σ(Mit ) „ .

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
8.3 Semiparametric estimation of the SPD 179

‚c dc ‚M dc 1
= = ,
˜ ˜
dM ‚ S dM K
d2 c
‚2c 1
= .
˜ dM 2 K
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:
d2 c
‚2c M dc M
= +2
‚K 2 dM 2 dM K 2

The functional form of is:

dc dd1 ¦ (d2 ) dd2 ¦(d2 )
’ e’r„ + e’r„
= ¦ (d1 ) , (8.14)
dM dM M dM
d2 c
while is:
dM 2

d2 c d2 d1 dd1
’ d1
= ¦ (d1 )
dM 2 dM 2 dM
e’r„ ¦ (d2 ) d2 d2 2 dd2 dd2
’ ’ ’ d2
M dM M dM dM
2e’r„ ¦(d2 )
’ (8.15)
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-
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 = ,
‚ 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
(m ’ m0 )2
g(m, „ ) ≈ g(m0 , „0 ) (m ’ m0 ) +
2 ‚M 2
‚M m0 ,„0 m0 ,„0
‚g 1‚ g
(„ ’ „0 )2
(„ ’ „0 ) +
2 ‚„ 2
‚„ m0 ,„0 m0 ,„0
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

σj ’ β0 + β1 (Mj ’ m0 ) + β2 (Mj ’ m0 )2 + β3 („j ’ „0 )
+β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 )
« 
(M2 ’ m0 )2 („2 ’ „0 )2
M2 ’ m0 „2 ’ „0 (M2 ’ m0 )(„2 ’ „0 ) ·
X = ¬. ·,
¬ ·
. . . . .
. . . . . .
. . . . . . 
(Mn ’ m0 )2 („n ’ „0 )2
Mn ’ m0 „n ’ „0 (Mn ’ m0 )(„n ’ „0 )
182 8 Estimating State-Price Densities with Nonparametric Regression

«  « 
σ1 β0
σ =  . , β =  . .
W = diag{KhM ,h„ (Mj ’ m0 , „j ’ „0 )}
. .
¬· ¬·
. .
σ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
β = 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
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
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:


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


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


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).

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


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


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).

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

6 1.20

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

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 ,

ˆ = 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)|.

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:
+ µi , i = 1, · · · , n.
Yi = σ
• 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:

µi = Yi ’ σh
˜ ˆ .

• Replicate B times the series of the {˜i } with wild bootstrap ob-
taining {µi } for j = 1, · · · , B, H¨rdle (1990), and build B new
190 8 Estimating State-Price Densities with Nonparametric Regression

bootstrapped samples:

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

• 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)|.


ˆ ˆ
• 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

2400 2600 2800 3000
Stock price at expiry S(T)

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

SPDs and bootstrap CB, tau= 49 days

0.8 0.9 1 1.1

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

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



2000 2500 3000 3500
stock price

Figure 8.7. Comparison of di¬erent SPD estimations, by Rookley™s
method (blue) and IBT (black, thin).

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
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.


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

±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

Pricing error in %

-20 -10

5 10 15 20 25 30

Pricing errors for butterfly spread
Pricing error in %


5 10 15 20 25 30

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

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
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
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
variance-function estimation, Technometrics 39: 262“273.
Yatchew, A. and H¨rdle, W. (2002). Dynamic nonparametric state price density
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
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 =

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




-0.3 -0.2 -0.1 0 0.1 0.2
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
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-
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

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

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:
um ’ u

pt (u)
ˆ = KM C (9.5)
M hM C hM C

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 = .

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


-0.3 -0.2 -0.1 0 0.1 0.2 0.3
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 ‚σ


. 7
( 14)