Mathematical Finance, Vol. 10, No.

2 (April 2000), 157–177

PRICING AMERICAN OPTIONS FITTING THE SMILE

M. A. H. Dempster and D. G. Richards

University of Cambridge

This paper is a compendium of results—theoretical and computational—from a series of recent

papers developing a new American option valuation technique based on linear programming (LP).
Some further computational results are included for completeness. A proof of the basic analytical
theorem is given, as is the analysis needed to solve the inverse problem of determining local
(one-factor) volatility from market data. The ideas behind a fast accurate revised simplex method,
whose performance is linear in time and space discretizations, are described and the practicalities
of fitting the volatility smile are discussed. Numerical results are presented which show the LP
valuation technique to be extremely fast—lattice speed with PDE accuracy. American options valued
in the paper range from vanilla, through exotic with constant volatility, to exotic options fitting the
volatility smile.
Key Words: American options, exotic options, linear programming, volatility smile, inverse
problem

1. INTRODUCTION
This article surveys the linear programming (LP) approach to fast valuation of American
options through the use of a special version of the revised simplex method. This special
algorithm makes use of the tridiagonal structure of the finite-difference discretization of
the Black–Scholes partial differential equation, a novel basis factorization, and the nature
of the optimal exercise boundary, to create a pricing algorithm that is essentially linear in
the number of discretization steps in space or time with the other held fixed. The method
is applicable to a variable coefficient version of the Black–Scholes equation which allows
the treatment and rapid solution of path-dependent exotic options taking account of local
volatility. One emphasis of the paper is a complete overview of the theory underlying
the LP approach to these difficult problems. Another is a representative survey of the
numerical results obtained to date.
The paper is structured as follows. In Section 2 we review the formulation of the
American put option valuation problem presented in Dempster and Hutton (1997a, 1999)
and Dempster, Hutton, and Richards (1998). The problem is a classical optimal stopping
problem that may be formulated as a free-boundary problem by considering the domain
partition of optimal stopping. Removing any explicit reference to the free boundary,
the option value may be seen to be the unique solution of an order complementarity
problem (Borwein and Dempster 1989) by considering its equivalent formulation as a
variational inequality and utilizing standard results for coercive operators. Finally, the
value is the solution of an abstract linear program (for which a simpler proof than previ-
ously is given) which can be solved with standard LP techniques upon suitable domain
Address correspondence to the author at the Centre for Financial Research, Judge Institute of Manage-
ment Studies, University of Cambridge, Cambridge CB2 1AG, United Kingdom (01223) 339700; e-mail:
m.dempster@jims.cam.ac.uk & d.richards@jims.cam.ac.uk; www-cfr.jims.cam.ac.uk.

© 2000 Blackwell Publishers, 350 Main St., Malden, MA 02148, USA, and 108 Cowley Road, Oxford,
OX4 1JF, UK.

157
158 m. a. h. dempster and d. g. richards

truncation and discretization. In Section 2 basic numerical methods and our variable
coefficient tridiagonal simplex algorithm, together with the degenerate PDE approach
(Wilmott, Howison, and Dewynne 1993; Barraquand and Pudet 1996) to valuing market-
traded, discretely sampled exotic options, are also reviewed.
Section 3 discusses both theory and basic numerical methods for pricing exotic options
with the local volatility surface implied by market values of European options. This is
an area pioneered in Rubinstein (1994) and Dupire (1997) and, although some new
theoretical proofs are given in Dempster and Richards (1999) and reliable numerical
methods are developed in this paper, its treatment here is mainly seen as a vehicle to
demonstrate the generality and efficiency of the LP valuation algorithm in Section 2.
In Section 4 results for FTSE 100 exotic American index options—fixed-strike Asian
puts—are presented to substantiate these claims. First, representative numerical results
for vanilla and constant volatility path-dependent American options are presented. Then
vanilla European and American options fitting the smile are studied in order to evaluate
potential pricing errors in fitting the local volatility surface used to finally price American
exotics. Conclusions are drawn in Section 5 and some directions of current and future
work are indicated.

2. PRICING AMERICAN OPTIONS USING LINEAR PROGRAMMING

Consider the standard (Black and Scholes 1973) economy, in which there are two finan-
cial instruments: a “risky” asset with price S modeled by a geometric Brownian motion
(GBM) and a savings account whose balance is continuously compounded at a constant
risk-free rate r ≥ 0. Define an equivalent martingale (or risk neutral) probability measure
(EMM) Q (see Harrison and Kreps (1979) and Harrison and Pliska (1981)) under which
the discounted stock price process e−rt S(t) is a martingale and the stochastic differential
equation (SDE) for the stock price process becomes the GBM dS = rS dt + σ S d W̃,
where t ∈ [0, T ], S(0) > 0, σ > 0 is the constant volatility of the stock price, and W̃ is
a Wiener process under Q.
A European (vanilla) call or put option confers the right (but not the obligation) to the
holder to buy or sell respectively one unit of the asset for a fixed strike price K exactly
at a maturity date T . The American equivalent on the other hand may be exercised at any
exercise time τ ∈ [0, T ]. Since under these assumptions an American call option will
optimally be held to maturity, we shall formulate a version of the American put problem
which is suitable for numerical solution. Define the value function v : R+ × [0, T ] → R
giving an option’s fair value v(x, t) to the holder at stock price x > 0 and time t ∈ [0, T ].
This value is partially determined by the payoff function ψ : R+ → R, which for the
American put is defined as ψ(S(τ )) := (K −S(τ ))+ received by the holder upon exercise
at a general stopping time τ ∈ [0, T ].
The value function of an American put option can be considered to be the solution
to a classical optimal stopping problem—choose the stopping time ρ(t) that maximizes
the conditional expectation under Q of the discounted payoff—and may be shown to
be the first time the value falls to the payoff at exercise, that is, ρ(t) := inf{s ∈
[t, T ] : v(S(s), s) = ψ(S(s))}. The domain of the value function can thus be parti-
tioned into a continuation region C, on which the option has value greater than the
payoff for early exercise, and a stopping region S, where the value equals the payoff.
Hence C := {(x, t) ∈ R+ × [0, T ) : v(x, t) > ψ(x)} and S := {(x, t) ∈ R+ × [0, T ) :
v(x, t) = ψ(x)}.
 pricing american options fitting the smile 159

On the continuation region, the value function satisfies the Black–Scholes parabolic
partial differential equation (PDE) LBS v + (∂v/∂t) = 0 for (x, t) ∈ R+ × [0, T ], where
the elliptic operator LBS := 21 σ 2 x 2 (∂ 2 /∂x 2 )+rx(∂/∂x)−r, since the discounted stopped
price process of the option is a martingale, while as soon as the process crosses into S,
v = ψ, and to preclude arbitrage LBS v + (∂v/∂t) ≤ 0. Hence we have
 
∂v
(2.1) −LBS v − ∧ (v − ψ) = 0
∂t

on the whole domain R+ × [0, T ), where ∧ denotes the pointwise minimum of two
functions. We now have a free-boundary formulation where v(x, t) = ψ(x, t) for (x, t)
on the optimal stopping or exercise boundary. We can remove any reference to the
optimal stopping boundary by formulating the problem in terms of (2.1) as a linear
order complementarity problem (OCP) using the log-transformed stock price variable
ξ := ln x, with respect to which the Black–Scholes operator is given
by −Lv − (∂v/∂t)
with constant coefficient elliptic part L := 21 σ 2 (∂ 2 /∂ξ 2 ) + r − 21 σ 2 (∂/∂ξ ) − r and v
is now the option value as a function of ξ . The various inequalities carry through
the
domain transformation and the new payoff function is given by ψ(ξ ) := K − eξ . It
will also be convenient to reverse time—to remaining time to maturity, for which we
will again use the symbol t for simplicity—so that the payoff function ψ becomes an
initial condition for the Black–Scholes PDE. The American put value function is then
the unique solution to


v(·, 0) = ψ


v ≥ ψ
(2.2) (OCP)
−Lv + ∂v


 ∂t ≥ 0


−Lv + ∂v ∂t ∧ (v − ψ) = 0 a.e. in R × [0, T ]

posed in a suitable vector lattice Hilbert space, which is a Hilbert space H with inner
product h·, ·i and partial order defined by a positive cone P such that for any points x
and y the maximum x ∨ y and the minimum x ∧ y exist in the given order (Cryer and
Dempster 1980; Borwein and Dempster 1989). Dempster and Hutton (1999) (see also
Jaillet, Lamberton, and Lapeyre (1990)) use another equivalent formulation of the value
function problem as a variational inequality (VI) to show the uniqueness of the solution
to (OCP) if the differential operator −L is coercive; that is, ∃α ∈ R+ s.t. hu, −Lui ≥
αkuk2 ∀u ∈ H . They show that the value function, as the unique solution to (OCP),
can be expressed as the unique solution of an abstract linear program given by

(2.3) (LP) inf hu, ci s.t. u ∈ F for any c > 0 a.e. on R × [0, T ],
u

where
 
∂u
(2.4) F := u ∈ H : u(·, 0) = ψ, u ≥ ψ, −Lu + ≥0 ⊂P
∂t

since the linear operator −L on the Hilbert space H is of type-Z; that is, v ∧ y = 0 (⇒
hv, yi = 0) ⇒ hv, −Lyi ≤ 0 ∀v, y ∈ H . This basic result giving equivalence between
160 m. a. h. dempster and d. g. richards

(OCP) and (LP) is an extension to the parabolic case of a result of Cryer and Dempster
(1980) for elliptic partial differential operators. To see this we shall need some notation.
Define the Sobolev space W m, p, µ (R2 ) as the space of functions u ∈ Lp (R2 , e−µ|x| dx)
whose weak derivatives of order not exceeding m ∈ N exist and are also in
Lp (R2 , e−µ|x| dx), for p ∈ [0, ∞] and µ ∈ (0, ∞). (Here | · | denotes the L1 norm
on R2 and dx denotes the Lebesgue measure on R2 , and it should be noted that the
extension of the results in the sequel to Rn+1 , for arbitrary n ∈ N, is completely straight-
forward.) Consider the Hilbert space H 1 (R2 ) := W 1, 2, µ (R2 ), for some fixed µ > 0, of
square integrable functions with square integrable derivatives defined on R2 . The Hilbert
space H 1 (R2 ) has as Banach dual the Sobolev space H −1 (R2 ) := W −1, 2, µ (R2 ), also a
Hilbert space of Radon measures, with which it may be identified. Consider the pairing
h·, ·i : H 1 × H −1 → R between dual spaces given by
Z
(2.5) hu, vi := u(ξ, t)v(ξ, t)e−µ(|ξ |+|t|) dξ dt,
R2

where we may interpret v ∈ H −1 as the density function of the Radon measure ele-
ment of the dual space H −1 of H 1 with respect to e−µ(|ξ |+|t|) dξ dt. Alternatively, we
may consider h·, ·i given by (2.5) as an inner product on the Hilbert space H 0 (R2 ) :=
L2 (R2 , e−µ|x| dx) by virtue of the canonical injections H 1 ,→ H 0 ,→ H −1 , see Baioc-
chi and Capelo (1984, p. 79). In this setting the partial differential operator L may be
interpreted either as a map H 1 → H −1 or as an operator on H 1 . Consider also the
bilinear form a(·, ·) : H 1 × H 1 → R given by
Z     
σ2 σ2 σ2 ξ
(2.6) a(u, v) := uξ vξ − r− +µ uξ v + ruv
R2 2 2 2 |ξ |
× e−µ(|ξ |+|t|) dξ dt.

which can be chosen to satisfy a(u, v) = hv, −Lui, u, v ∈ H 1 .

Finally, note that H 1 (and hence H −1 ) is a vector lattice Hilbert space (but not a
Hilbert lattice) with positive cone defined in terms of (Lebesgue) almost everywhere
nonnegativity. See Cryer and Dempster (1980); Baiocchi and Capelo (1984) and Borwein
and Dempster (1989, p. 553–554) for more details on these ideas, which have been
adapted here to match the more general setting of Jaillet et al. (1990). In particular,
we shall assume all functions in H 1 (∼ = H −1 ) considered to be defined as u(·, |t|) on
R × (−∞, 0) and as u(·, T ) on R × [T , ∞) (see Cryer and Dempster (1980, pp. 89 ff)).
The following results (Cryer and Dempster 1980; Dempster and Hutton 1999) relate
the bilinear form a(·, ·) to the elliptic part of the partial differential operator L and show
variational inequality (VI) and (OCP) equivalence.

Theorem 2.1. The variational inequality (VI) given by



v(·, T ) = ψ

(2.7) (VI) v≥ψ


u ≥ ψ a.e. ⇒ a(v, u − v) + u − v, ∂v
∂t ≥ 0 a.e. in [0, T ]

is equivalent to the order complementarity problem (OCP) and is uniquely solvable.

 pricing american options fitting the smile 161

It can be shown that a given by (2.6), and hence −L, is coercive (see Jaillet et al.
(1990, p. 267), whose spaces L2 ([0, T ], Vµ ) and L2 ([0, T ], Hµ ) may be considered
restrictions respectively of our spaces H 1 and H 0 ). Then the Lions–Stampacchia theorem
(Baiocchi and Capelo 1984) implies that the solution to (VI) is unique. The formulation
(VI) is a type of classical physical problem, termed the (Stefan) obstacle problem, where
the payoff function ψ is the obstacle below which the solution cannot fall.
Next define, for a closed subset F ⊆ P ⊆ H of a vector lattice Hilbert space H with
positive cone P := {u ∈ H : u ≥ 0}, the least element problem

(LE) v∈F s.t. v ≤ u, u ∈ F.

The least element v is denoted by LE(F ). Note that if it exists, the least element is
always unique since, if v1 and v2 are least elements of F , then v1 ≤ v2 and v2 ≤ v1 , so
from the vector lattice property v1 = v2 .

Theorem 2.2. In the setting described above, if T is a coercive type Z temporally

homogeneous elliptic differential operator and F is given by (2.4), then there exists a
unique solution v to the following equivalent problems:


v(·, 0) = ψ



v ≥ ψ
(OCP) ∂v

T v + ∂t ≥ 0



(T v + ∂v ) ∧ (v − ψ) = 0 a.e. R × [0, T ],
∂t

(LE) find v = LE(F ),

(LP) inf hv, ci s.t. v ∈ F , for any c > 0 a.e. on R × [0, T ]
v

Proof. We first prove the equivalence between (OCP) and (LE), after making the
trivial domain extensions of the problem functions given above to set them in H 1 . Let
L denote the Laplace transform operator with respect to the measure e−µ|t| , so that, for
(ξ, λ) ∈ R2 , the Laplace transform û ∈ H 1 of a function v ∈ H 1 is defined by
Z ∞
(2.8) û(ξ, λ) := Lu(ξ, ·)(λ) := e−|λ|t u(ξ, t)e−µt dt.
0

As noted above, we have extended the temporal domain of our value functions v to [0, ∞)
as constant on (T , ∞), so that this generalized Laplace transform is well defined. L is
a linear operator and T is temporally homogeneous (i.e., it has time-independent coef-
ficients), and therefore commutes with the Laplace operator, so that taking the Laplace
transform of the operator T + ∂t∂ gives T L+L ∂t∂ . The Laplace transform of the first-order
time derivative is given by
  Z ∞
∂v ∂v
(2.9) L (ξ, λ) := e−|λ|t (ξ, t)e−µt dt
∂t 0 ∂t

= −v(ξ, 0) + (|λ| + µ)v̂(ξ, λ)

and v(ξ, 0) is given by the initial condition v(·, 0) = ψ (in backwards time).
162 m. a. h. dempster and d. g. richards

Now, note that the Laplace transform is positivity-preserving in the sense that v ≥ 0 ⇒
v̂ ≥ 0 a.e. on R2 . Then, writing the initial condition, constant in λ, as q̂(·, λ) :≡ −v(·, 0)
to agree with the notation of Borwein and Dempster (1989), (OCP) is equivalent to the
transformed order complementarity problem (OCP), \ also posed in H 1 , given by


v̂ ≥ ψ̂
\
(OCP) T v̂ + q̂ ≥ 0


(T v̂ + q̂) ∧ (v̂ − ψ̂) = 0 a.e. on R2 ,

where T̂ := T + |λ| + µ and ψ̂ is the Laplace transform of the log-transformed payoff

function ψ, given by ψ̂(ξ, λ) = ψ(ξ )/(|λ| + µ). T̂ remains coercive and it is easy
to check from the definition that it remains type Z. Consider solutions to the pro-
jected (OCP) obtained from (OCP) \ by fixing λ ∈ R. We can now apply the order
complementarity–least element equivalence result of Borwein and Dempster for coercive
type Z elliptic operators, so that for each λ ∈ R one such, v̂(·, λ), is the solution (neces-
sarily unique) to the least element problem defined by LE(F̂λ ), where Fˆλ is defined by

(2.10) Fˆλ := v̂(·, λ) : v̂(·, λ) ≥ ψ(·, λ), T v̂(·, λ) + q̂(·, λ) ≥ 0 .

It follows that v̂ is the unique solution to the least element problem (LE)d defined by
LE(F̂ ), where F̂ is defined by F̂ := {v̂ : v̂ ≥ ψ̂, T v̂ + q̂ ≥ 0}. Applying the inverse
Laplace transform L−1 to v̂ shows that v = L−1 v̂ solves both (LE), given by LE(F ),
and (OCP), as required. Indeed suppose the contrary: that there exists u ∈ F such that
u ≤ v, u 6= v. Then it follows since L is positivity preserving that û ∈ F̂ and û ≤ v̂,
û 6= v̂, a contradiction to v̂ = LE(F̂ ).
With this least element result, the LP equivalence is immediate: v is the least element
of F ⇐⇒ v ≤ u for all u ∈ F , and so hc, vi ≤ hc, ui for all u ∈ F and any vector
c > 0. Therefore u minimizes hc, vi over all v in F and is thus the solution to the abstract
linear program (LP). Restricting to the original problem domain yields the result. 2

It should be noted that the above proof depends on time running “backwards” since
otherwise we cannot substitute ψ for v(·, 0) in (2.9). The least element result tells us
that the linear constraint set lies within the positive cone translated so that its apex lies
at v. We can pick out the least element of the constraint set by minimizing hc, ui over
the set u ∈ F , where c > 0; specifically in R2 by minimizing the intercept of negatively
sloped lines defined by c0 u with normal c > 0 intersecting F .
Theorem 2.2 gives equivalence between (VI), (OCP), (LE), and (LP) for the American
put, since −L is coercive type Z (see Jaillet et al. 1990). It should be stressed that the
result is very general and applies to virtually any parabolic partial differential operator
with a temporally homogeneous coercive type Z elliptic part, and virtually any payoff
function. For example, it may be applied to the Black–Scholes operator −LBS directly
without prior logarithmic transformation of the space variable. A more delicate argument
involving step function coefficient approximation and a suitable passage to the limit can
be used to establish the results of Theorem 2.2 for time-dependent coefficient operators
(the details will appear elsewhere) which are required in this paper to handle the time-
varying nature of local volatility coefficients. Theorem 2.2 also suggests a simple way
to solve the equivalent problems numerically: by a suitable discretization the infinite-
dimensional abstract linear program (LP) reduces to an ordinary linear program with
solutions in Rn . Thus we next discretize the problem and consider efficient LP numerical
solution methods.
 pricing american options fitting the smile 163

2.1. Computation
We shall approximate the value function by a function that is piecewise constant
on rectangular intervals between points in a regular lattice. Approximating the partial
derivatives by standard Crank–Nicolson finite differences (see, e.g., Wilmott et al. 1993)
we obtain a discrete form of (OCP) that can be rewritten in matrix form upon collapsing
the space index. A matrix is type-Z if it has nonnegative off-diagonal elements (see,
e.g., Borwein and Dempster 1989), which in the case of the matrix derived from the
discretized negative Black–Scholes operator occurs when |r − σ 2 /2| ≤ σ 2 /1ξ and can
be satisfied by adjusting the number of space steps I in the discretization (Dempster
and Hutton 1999). From this condition it can also be shown that this matrix is coercive
(Jaillet et al. 1990; Hutton 1995).
The LP formulation can be solved in backwards time either directly or iteratively
and the interested reader can find comparisons of solution methods in Hutton (1995) and
Dempster and Hutton (1999). Here we solve the discretized (LP) using time-stepping and
a simplified revised simplex algorithm that takes advantage of the tridiagonal structure
of the constraint matrix formed from standard Crank–Nicolson finite difference approx-
imations to produce a fast accurate direct solution method detailed in Dempster et al.
(1998) and Richards (1999). This procedure is suitable for any standard constant param-
eter Black–Scholes type formulation, but also yields significant computational savings
for valuation problems with volatility and drift parameters which are functions of time. It
incorporates a technique for the solution of problems with nonconstant constraint matrix
coefficients such as those involving the untransformed Black–Scholes PDE, which has
coefficients given by functions of the underlying asset price, or for exotic option pric-
ing problems, where the coefficients vary with the third variable representing the path-
dependency. In Dempster et al. (1998) results are presented for this updating procedure
which show that even for a general constraint matrix the procedure outperforms standard
commercial LP solvers by orders of magnitude.
To understand the exact computational savings of these simplex methods, first consider
the complexity of the vanilla American put option valuation problem after transformation
to the constant-coefficient Black–Scholes operator. At each time-step the maximum num-
ber of real variables which can enter the simplex basis is O(I ) and hence we have O(I )
iterations at each time step, where I is the number of points in the spatial discretization.
In fact, after the first few time steps—where the exercise boundary has greatest curvature
away from ln K (see Figure 4.3 later)—at most one new basic variable enters at each
time step. Far from maturity, calculations for several time steps may even utilize the
same basis. Each iteration requires O(n) operations to solve, where n ≤ I , giving O(I )
operations at each time step. Hence the space complexity of the algorithm is linear and
the total operation count is O(MI ), where M is the number of time steps.
For the basis factorization updating technique required by each simplex iteration (space
step) the calculations result in a similar complexity, but can be performed in three floating
point operations, although extra computation time is needed for the dynamic allocation
of the upper-lower (UL) factorization. Results for the constant coefficient method and
for the nonconstant coefficient updating technique are reported in Dempster et al. (1998),
along with results for a complete calculation of the full LU factorization at each iteration
to highlight the overheads of using general commercial solvers.

2.2. Extension to Exotic Options

An exotic option is any derivative security that has a path-dependent component in its
payoff at exercise. We may formulate discretely sampled exotic option valuation problems
164 m. a. h. dempster and d. g. richards

as linear programs through state augmentation, particularly for American discretely sam-
pled lookback and Asian options. Exotic option values are dependent on the underlying
stock price, (forward) time, and an additional “independent” variable that encapsulates
the required path information.
We now outline the formulation of a generic American exotic option in a discretely
sampled setting using the unifying framework of Wilmott et al. (1993). Denote by
V (S, M, t) the value function of the option with V : R+ × R+ × [0, T ] → R, where S
denotes the asset price and M denotes the value of a path-dependent variable, such as the
average in the case of Asian options, or the maximum/minimum for lookback options.
Assuming that the asset price is sampled on N occasions during the life of the option
with maturity T , denote by Mn the observed value of the augmented variable at the
sampling date tn , n = 0, . . . , N − 1, so that sampling begins at time 0 and t0 = 0. For
completeness, define tN := T . The variable Mn is a constant value throughout the period
[tn , tn+1 ), since no sampling takes place until time tn+1 . Effectively Mn is a parameter in
the formulation during this period and any randomness in the model is due to the asset
price process. The Black–Scholes PDE will thus be satisfied within the period with jump
conditions applied at sampling dates (see Wilmott et al. 1993 for more details). Across
a general sampling date tn the augmented variable is updated from a value Mn−1 just
prior to the date to a value Mn
at the sample date. No-arbitrage arguments lead to the
jump condition V S, Mn−1 , tn− = V (S, Mn , tn ) n = 1, . . . , N − 1, where tn− is a time
immediately before the sampling date tn .
In the time interval [tn , tn+1 ) the European value V satisfies the augmented Black–
Scholes PDE defined by LBS + f (S, t)∂V /∂M + 2V /2t, where f (S, t) is a function
specified for the option of interest. We consider the final period [tN −1 , T ] and use
a dynamic programming algorithm to determine values for earlier periods. As in the
vanilla American put case treated above, the American exotic valuation domain in R3
can be partitioned into a continuation region CN and a stopping region SN and we can
establish the existence of an optimal exercise boundary (Dempster and Richards 1999).
To complete the formulation of the discretely sampled exotic value in the final period
we require a terminal condition V (S, MN−1 , T ) = ψ(S, MN −1 ) for all S, MN −1 ∈ R+
and boundary conditions in S at S = 0 and as S → ∞ which are option dependent
and are discussed in more detail in the cited paper. If again we log-transform the prim-
itive variables (ξ := ln S, ζN−1 := ln MN−1 ) and formulate the valuation problem
with fixed ζN−1 as an OCP with respect to the transformed operator L, we may define
a new partition with regions CN and SN . Thus the American exotic valuation prob-
lem in the final period
may be formulated in terms of the transformed value function
V := V eξ , eζN −1 , t as the unique solution of the order complementarity problem of


the form (2.2) over R2 × (tN−1 , T ] involving payoff ψ(ξ, ζN −1 ) := max eζN −1 − eξ , 0
with V also denoting the option value as a function of ξ and ζN −1 . This puts us in a
framework equivalent to that of the vanilla American put of Section 2 but with the addi-
tional parameter ζN−1 , and hence we have equivalence to an abstract LP for each value of
the augmented variable ζN−1 ∈ (−∞, ∞). This problem must be solved for all possible
values of the parameter ζN−1 . Applying the jump conditions at tN −1 to obtain the termi-
−
nal value V (S, MN−2 , tN−1 ), the argument may be repeated for the period [tN −2 , tN −1 ]
and, by backwards recursion, eventually for the period [0, t1 ].

3. THE INVERSE PROBLEM OF OPTION PRICING

The Black–Scholes economy has one unobservable quantity—the volatility parameter
σ —which must be inferred. Particularly since the stock market crash of 1987 the volatility
 pricing american options fitting the smile 165

of equity options has exhibited variation both with the strike price—the volatility smile
(or skew)—and the option’s maturity—the volatility term structure. These effects high-
light the market’s deviation from the Black–Scholes assumption that the future asset
price has a constant variance lognormal conditional probability density.
Several approaches have been suggested in the literature to model this behavior. One
approach is to treat the volatility as a second stochastic factor with the aim of spec-
ifying the time-varying volatility from the model (Hull and White 1987; Clarke and
Parrott 1998). However this approach is difficult to fit to market data, is not arbitrage
free, and introduces an additional dimension to the pricing problem. An alternative one-
factor approach allows the volatility σ := σ (S, t) to be a variable that is both state and
time dependent. By starting from the market data and finding the local volatilities that
are consistent with the market—commonly termed an inverse problem (Andersen and
Brotherton-Ratcliffe 1997; Lagnado and Osher 1997; Bouchouev and Isakov 1999)—
this model can be made to price the market nearly exactly. The most popular structures
on which this local volatility is determined are binomial or trinomial trees, which allow
specification of nodal transitional probabilities to fit the smile (Rubinstein 1994; Der-
man and Kani 1998). All methods used to fit market data are prone, however, to the
same instabilities. Generally the data can imply unreasonable (e.g., negative or large)
values of the local volatility, which may create negative transitional probabilities and
allow arbitrage possibilities in the model. The general inverse problem is ill-posed since
the number of volatility parameters to be found far outnumbers the limited number of
available option prices in the market. Therefore, it is often assumed that a continuum of
European call option prices C(K, T ) are available for all strikes and maturities. Alterna-
tively, recent papers have implemented regularization methods (Jackwerth and Rubinstein
1996; Bodurtha 1997; Lagnado and Osher 1997; Coleman, Li, and Verma 1999) to make
the inverse problem stable.
All approaches to fitting the smile suffer from the consequences of inconsistent data
and will not price correctly all options—particularly those far out of the money—in the
face of these data problems. The modeling approach used in this paper is not claimed
to be the most accurate or efficient, but it is quick and highlights the versatility of
the LP pricing method in the face of a degenerate ill-posed problem with nonconstant
coefficients.

3.1. Continuous-Time Volatility Theory

An arbitrage-free local volatility surface in continuous time can be inferred from mar-
ket data, in particular the prices of European call options. This theory was first derived
by Dupire (1997) and has since been given a more formal treatment in Derman and
Kani (1998). The main idea is that there exists an adjoint or dual PDE to the Black–
Scholes PDE with strike price K and maturity T as the independent variables, which
can be derived through consideration of the conditional probability distribution of the
underlying stochastic process and the forward PDE satisfied by this probability density.
We assume that a continuum of European call option prices C(K, T ) for all strikes
and maturities K, T ∈ R+ are available from the market and any gaps in the data can
be filled by interpolation or extrapolation techniques. We will deal with any arbitrage
violations in these approximated values later. The underlying asset price is assumed to
follow a diffusion process under the risk-neutral measure Q, but now with nonconstant
volatility; that is, dS = Sr(t) dt + Sσ (S, t) dW. The price of a European call option can
be written in terms of an expectation under Q with respect to the conditional probability
166 m. a. h. dempster and d. g. richards

density function (pdf) p(s, T |S, t) of the underlying asset S having value s at time T
given that the asset price at time t is S as
Z ∞
(3.11) C(S, t; K, T ) = P (t, T ) p(s, T |S, t)(s − K)+ ds,
0
RT

where t ≤ T and the discount factor P (t, T ) := exp − t r(u) du .
Breeden and Litzenberger (1978) showed that the European call option price and this
conditional probability density were related as p(K, T |S, t) = P (t, T )−1 ∂ 2 C(S, t; K, T )/
∂K 2 by differentiation of (3.11). The function p(K, T |S, t) is also the Green’s function
(or fundamental solution) of the Black–Scholes PDE. Thus it satisfies this PDE with
terminal condition (t = T ) p(K, T |S, t) = δ(S − K), where δ(·) is the Dirac delta func-
tion. Since C is assumed known, this density function can be found from the idealized
market data. The function p(K, T |S, t) can also be shown to satisfy the Fokker–Planck
(or forward Kolmogorov) PDE utilizing the following theorem (see, e.g., Jackwerth and
Rubinstein 1996).

Theorem 3.1. The conditional pdf p(y, τ |x, t) of a general stochastic process X(t)
where t ≥ 0 given by dX(t) = µ(X, t) dt + σ (X, t) dW(t) satisfies the Fokker–Planck
or forward Kolmogorov equation


∂p(y, τ |x, t) ∂ (µ(y, τ )p(y, τ |x, t)) 1 ∂ 2 σ 2 (y, τ )p(y, τ |x, t)
(3.12) + − = 0,
∂τ ∂y 2 ∂y 2

for fixed (x, t) ∈ R × R+ with initial condition p(y, t|x, t) = δ(x − y).

Corollary 3.2. The transitional probability density function p(S 0 , T |S, t) of the
stock price process dS/S = r(t) dt + σ (S, t) dW satisfies the PDE

∂p(S 0 , T |S, t) 1 ∂ 2 (σ 2 (S 0 , T )S 0 2 p(S 0 , T |S, t) ∂

(3.13) = 2
− 0 rS 0 p(S 0 , T |S, t)
∂T 2 ∂S 0 ∂S
with initial condition (T = t) p(S 0 , T |S, t) = δ(S 0 − S).

The following corollary is due to Dupire (1997).

Corollary 3.3. Given that the underlying asset price process in Corollary 3.2 the
price of a European call option C(S, t; K, T ) solves the partial differential equation

∂C σ (K, T )2 K 2 ∂ 2 C ∂C
(3.14) = − r(T )K
∂T 2 ∂K 2 ∂K
with boundary condition C(S, t; S, t) = 0 (see also Derman and Kani 1998).

3.2. Computation
The local volatility function σ (S, t) can be fully determined from (3.14) since all other
terms in the equation can be found from market data. However, since the inverse problem
is ill-posed there may not be a unique local volatility functional that fits the market data.
To obtain an approximately arbitrage-free pricing algorithm the implied volatility (or,
alternatively, call price) data available from the market must be fitted exactly. To this
 pricing american options fitting the smile 167

end we apply cubic splines to approximate the implied volatility surface. The ease with
which approximations to first and second derivatives can be obtained from a spline fit
to market data is a major advantage in reducing the computational complexity of our
method since the construction of a length N cubic spline is the solution of a tridiagonal
system and hence an O(N) calculation. The complexity of the problem can be further
reduced by precomputing derivative information, but at the expense of additional memory
requirements.
We outline now the approach used to obtain a consistent local volatility surface that
fits the market-implied volatility data for European call options. This will be used in the
linear programming approach of Section 2 to price American exotic options consistent
with the volatility smile. Since we consider the implied volatility nonconstant in the
underlying diffusion, assume that the quoted European call prices are Black–Scholes
prices with implied volatilities as an additional parameter and define the call prices in
terms of the implied volatility ν(K, T ) for the call option of strike K and maturity T as
C(K, T ) := CBS (S, t; K, T , ν(K, T )), where CBS is the Black–Scholes European call
price. Using this formulation we can write the local volatility in terms of derivatives of
the Black–Scholes implied volatility as follows.1

Theorem 3.4. Let the asset price S follow the diffusion process in Corollary 3.2.
Then the local volatility function σ (S, t) consistent with the arbitrage-free European
call prices is given uniquely, in the absence of dividends, by
  −1
∂C ∂C ∂ 2C
(3.15) σ 2 (K, T ) = 2 + r(T )K K2
∂T ∂K ∂K 2

with S = K. In terms of the implied volatility function ν(K, T ) this can be written as

(3.16)  
∂ν 1ν ∂ν
2 γ ∂T + 2γ + r(T )K ∂K
σ 2 (K, T ) =   
,
1 ∂ν 1 ∂ν 1 ∂2ν γ ∂ν
K2 ν ∂K ((T − t) − d1 ) − Kγ −d1 ∂K γ− K + ∂K 2
γ + K ∂K

where


ln(S/K) + r + 21 ν(K, T )2 γ 2 √
d1 := and γ = T − t.
ν(K, T )γ

Proof. Equation (3.15) follows immediately from (3.14). By portfolio dominance

arguments (see Andersen and Brotherton-Ratcliffe 1997) it can be shown that σ 2 (K, T )
is nonnegative in the absence of arbitrage if the numerator of (3.15) is nonnegative. The
denominator is the Breeden–Litzenberger conditional (transitional) probability density
which must be nonnegative. Equation (3.16) follows from (3.15) after much calculation
(see Richards 1999). 2

We will use (3.16) to obtain the local volatilities. For our method the implied volatili-
ties are relatively stable and the expression for the local volatility involves relatively few
1
This result was derived following discussions with S. H. Babbs (of Bank One) and has since been inde-
pendently presented in Andersen and Brotherton-Ratcliffe (1997).
168 m. a. h. dempster and d. g. richards

numerical calculations and only one computationally expensive logarithmic calculation

for each evaluation point. Splines are fitted to the market implied volatility data, with the
calculated second derivatives with respect to the maturity stored in an array. By fitting
the splines across maturities for each strike first we obtain approximations for the first
and second derivatives with respect to the strike. However the first-order derivative of the
volatility with respect to maturity is not a natural by-product of the interpolation and is
specified here by a simple first-order approximation. Given strike and maturity values at
a mesh point the cubic spline interpolation is sufficient to supply all the values required
for the calculation of the local volatility in (3.16).
Since the option valuation takes place on a log-transformed grid for consistency with
the earlier developed methods, the strike prices need to be transformed so that the same
grid can be utilized for calibration and pricing. At each node (i, m) of the transformed
grid the local volatility is calculated from the spline approximation and the array of
calculated volatilities is stored to be used later in the pricing algorithm.
The pricing procedure follows from that developed for the Black–Scholes model in
Section 2 except that the volatility in the formulation is no longer assumed constant.
This changes the constraint matrix in the order complementarity problem, and thus in
the linear programming formulation; the matrix now has nonconstant diagonals but is still
tridiagonal in nature. However, for American and exotic options this is the ideal problem
to be solved by the nonconstant tridiagonal simplex method described in Section 2.
Several studies have aimed at detailing the instabilities in market data and why they
occur. If a model produces arbitrage opportunities we must use some regularization or
filtering procedure applied to the original market data to remove them. Some practition-
ers suggest setting any negative probability densities to zero (implying an infinite local
volatility) or restricting any local volatilities to lie within a range (σmin , σmax ), where
the bounds are supplied somewhat arbitrarily. We filter the data using methods derived
from the underlying theory. When upon discretization a negative denominator occurs
in (3.16) we necessarily have an arbitrage opportunity appearing since this implies a
negative value of the transitional probability density. We correct this value using put-call
parity to consider the prices of European put options implied by the call option market
data. Indeed, if the market data imply a negative value of p(Kj ) for some j = 0, . . . , N,
where p(·) := p(·, T |S, t) now denotes the discrete form of the conditional density, we
consider a European put option of strike Kj +1 with value P (S, t; Kj +1 , T ). Then

j
X
(3.17) P (S, t; Kj +1 , T ) ≈ P (t, T ) (Kj +1 − Ki )p(Ki )1K
i=0
j −1
X
= P (t, T ) (Kj +1 − Ki )p(Ki )1K + P (t, T )(Kj +1 − Kj )p(Kj )1K,
i=0

which given that the P (S, t; Ki , T ) are known for all i ≤ j (from put-call parity) allows
us to find a value of the discrete probability p(Kj ) consistent with the market data. If
this probability is also negative, its value is set to zero. The other possible inconsistency
in the data occurs when the numerator of (3.16) is negative, when we set σ 2 (K, T ) to
zero.
 pricing american options fitting the smile 169

4. NUMERICAL RESULTS
In this section we present empirical results for the procedure outlined in Section 3 for
pricing American options consistent with the observed market volatility smile, together
with benchmark results for their constant volatility equivalents. To facilitate comparisons
we use as our underlying implied volatility surfaces market data appearing previously
in the literature and we price European, American, and exotic options with respect to it.
All solution times quoted are for calculations on an IBM RS6000/590 workstation with
1 GB RAM running under AIX 4.3, although only a small proportion of this memory is
utilized. Results are quoted in Dempster et al. (1998) for solution on a Pentium II 400
Mhz PC which gives significant speed-ups over those presented here for most levels of
domain discretization.

4.1. Vanilla American Put

Table 4.1 illustrates the savings that the new tridiagonal simplex solver makes over
the PSOR algorithm2 (Cryer 1971). The timings in Table 4.1 are CPU times, including
all data initialization for the value at time 0 of an at-the-money American put with
parameters K := 1.0, T := 1.0, σ := 20%, and r := 10%. The log stock price was
bounded above by U := 2 and below by L := −1, giving the range in untransformed
variables as [0.37, 7.39]. The number of time steps M was set at 1000 and the number
of space steps I varied.
We see from the table that all of our tridiagonal revised simplex methods are linear in
space and give impressive speed-ups over PSOR, with the constant-coefficient method
(column 2) ranging from 4–500 times faster. For comparison, results are included for
the UL update technique (column 3) and solution times for the tridiagonal solver with
recalculation of the whole decomposition at each iteration (column 4). The slowest of
our solution methods is faster than PSOR at all but the lower levels of discretization and
all our methods are able to price 4 or 5 options per second accurately.

Table 4.1
Comparison of Tridiagonal Simplex Solvers with PSOR

Tridiagonal simplex
Space Constant Recalculation
steps coefficients UL update simplex PSOR
75 0.02 0.05 0.10 0.07
150 0.05 0.08 0.17 0.13
300 0.10 0.14 0.33 0.27
600 0.19 0.24 0.65 1.25
1200 0.38 0.47 1.26 6.37
2400 0.77 1.00 2.47 37.55
4800 1.61 2.24 5.12 255.06
9600 3.71 5.09 10.77 1856.91
Solution CPU times in seconds for M := 1000.

2
See Dempster and Hutton (1999) for comparison of the commercial IBM Optimization Subroutine Library
(OSL) with PSOR, interior-point, and simplex methods.
170 m. a. h. dempster and d. g. richards

4.2. Discretely Sampled Lookback Strike Options

Next we compare our results for discretely sampled American lookback strike options
against some of those in the literature. A lookback strike option has a payoff similar to
the corresponding (put or call) vanilla option, but with the strike price replaced by the
maximum or minimum realized asset price. For example, an American lookback strike
put might have payoff ψ(S, M) := max(M −S, 0) at exercise, where M is the maximum
asset price over the life of the option until exercise, and S is the asset price at exercise as
usual. We concentrate on strike options here rather than on rate options (Wilmott et al.
1993) which are treated below. The valuation of American discretely sampled lookbacks
is usually achieved using tree methods, though some closed-form solutions are available
for the European continuously sampled case.
Our main source of comparison for our numerical results is the PSOR method utilized
by Wilmott et al. (1993) after a similarity transform of the augmented Black–Scholes
PDE. The sampling schemes investigated, A, B, and C, correspond approximately to
samples every 1, 2, and 3 months respectively. The lookback option is of one-year
maturity with σ := 20% and r := 10% and the option is valued at-the-money with initial
stock price S := 100. Further details and additional results can be found in Dempster
et al. (1998) and Richards (1999). Table 4.2 shows the CPU timings and solutions for
the sampling schemes. These results are a good comparison because of the explicit
description of the sampling schemes employed, but are quoted in terms of the similarity
reduced variable to only one decimal place. We have agreement in the results to the
accuracy quoted, with solution times of much less than 6 seconds using the constant-
coefficients simplex solver, which far outperforms the PSOR method. The higher-order
discretizations for the former are given to show convergence. Further comparison results
for discretely sampled Asian strike options may be found in Dempster et al. (1998).

4.3. Fitting the Volatility Smile

To test the numerical procedure described in Section 3 we use real FTSE 100 index
volatility values implied from the European call option data described in Duan (1995).
The implied volatility data are shown graphically in Figure 4.1 and correspond to FTSE
100 volatility values for 31 March 1995. The initial index level was 3129.5 and data were
quoted for eight strike prices and five maturities, although for the last two maturities
prices were quoted for different strikes. The data were interpolated to fill the gaps as

Table 4.2
Comparison of Lookback Valuation Results for the Tridiagonal Method against the
Similarity Transformed Method

Tridiagonal
Similarity
200 × 300 × 200 400 × 600 × 400 800 × 1200 × 800
Transform
Scheme Solution Time Solution Time Solution Time Solution
A 10.532 6.46 10.550 53.88 10.555 477.32 10.5
B 9.445 6.69 9.454 59.05 9.457 508.95 9.5
C 8.114 6.59 8.116 59.32 8.116 489.60 8.1
 pricing american options fitting the smile 171

Figure 4.1. Implied volatility surface for European call options on the FTSE 100 index.

described in Section 3. We also assume that a constant rate of interest, r = 10%, applies
throughout for all maturities.
As described in Section 3 the pricing procedure occurs in two stages. The first—
calibrating the local volatility—uses the estimates of the implied volatility derivatives to
calculate the surface σ (S, t) for use in pricing. In the sequel this calibration takes place
on the same mesh that is used in the pricing procedure. There are significant compu-
tational savings to be made by calculating the local volatilities on a coarser mesh and
performing some type of interpolation between the calibrated values. The initial cali-
bration of the surface is computationally expensive. To calculate a volatility surface for
a discretization of 800 time-steps and 200 space-steps takes approximately 20 seconds,
although the time required for subsequent valuations on the same grid is a fraction of
this if the surface is stored for future use. We quote solution times for the option val-
uation only, assuming that the precomputed local volatility surface is in memory. After
calibrating the volatility surface, the options are valued using the UL update algorithm
of Section 2 for options with American exercise and a simple linear equation solver
for European options. At each time-step for the former, the basis decomposition was
calculated and the UL update applied only when new variables entered the basis.
Empirical tests showed that the underlying discretization error due to the Crank–
Nicolson numerical approximation of the derivatives was less than 8 basis points on a
grid of 3200 time steps and 200 space steps; this must be accounted for in any accuracy
comparisons.
172 m. a. h. dempster and d. g. richards

4.4. European Option Results

We now use our pricing algorithm to recover market option values in order to assess
the fit of the calibrated surface. Table 4.3 contains European call option prices for the
strikes and maturities given by the market. The reference (actual) values are calculated
through the use of the Black–Scholes pricing equation using the FTSE 100 implied
volatility values. This requires an approximation for the cumulative normal distribution
which is accurate to 10 decimal places. The errors can be seen to be comparable to—and
in some cases more accurate than—the baseline numerical accuracy described earlier; in
fact all pricing errors are less than 5 basis points at the discretization level used. We
conclude that fitting the volatility smile does not induce significant errors above the
baseline accuracy into the option values for European call options. Since the original
European call options can be seen to be accurately priced, the volatility surface fit to the
data is consistent with the local volatilities implied by the market through quoted prices.
Solution time for a European option at the discretization employed is approximately 0.6
seconds.
The calculated local volatility surface is displayed in Figure 4.2 and a comparison
with the implied volatility surface of Figure 4.1 shows significant differences. The local
volatility surface is represented in Figure 4.2 on a truncated strike domain. At short
maturities a spike of local volatility occurs at a strike value of approximately 3250, which
distorts any graphical representation of the local volatility surface but does not cause
instabilities in the calculated option values. For ease of representation, the local volatility
is shown for all strikes less than the level at which the spike occurs and illustrates that
for short maturities the local volatility behaves like the reciprocal of the conditional
transitional probability density, smoothing out somewhat for higher maturities.

4.5. Pricing American Options

Next we introduce an added dimension to the problem by pricing American options,
since there is a very real possibility that the vanilla optimal exercise boundary (see
Figure 4.3) will be moved by fluctuations in the local volatility surface. The solution
is a modification of the valuation of European options in the previous section, with the
updating tridiagonal LP solver introduced in place of the tridag linear equation solver.
The boundary conditions used in this section are as described in Dempster et al. (1998).

Table 4.3
Pricing errors ×100 of FTSE 100 European Call Options Fitting the Smile

April May June

Strike
Price Value Actual Error Value Actual Error Value Actual Error

2975 158.84 158.87 3.50 189.08 189.09 1.00 221.23 221.25 2.30
3025 114.92 114.88 4.40 148.82 148.80 1.80 181.84 181.82 1.70
3075 74.54 74.52 1.63 112.51 112.52 0.70 146.10 146.11 1.20
3125 44.43 44.46 2.89 80.69 80.72 3.32 114.15 114.19 4.00
3175 23.50 23.51 1.09 54.92 54.91 0.54 86.11 86.13 1.93
3225 9.77 9.79 2.34 35.13 35.08 4.88 62.02 61.98 3.99
3275 4.15 4.17 2.10 21.24 21.21 2.62 42.85 42.84 0.94
3325 1.09 1.06 2.79 11.85 11.82 3.48 27.83 27.80 3.36

M = 3200, I = 200. Actual value given by the Black–Scholes call pricing equation.
 pricing american options fitting the smile 173

Figure 4.2. Local volatilities σ (S, t) for FTSE 100 index options on a truncated strike
domain. Parameters M = 200, I = 200 for the full strike range.

The results for the FTSE 100 index in Table 4.4 show American put valuations for the
LP tridiagonal solver. As a benchmark we use the original LP valuation on a very fine
solution mesh using the at-the-money BS implied volatility. American option values are
higher than the corresponding constant volatility values to within the numerical tolerance
previously evaluated, with the LP solution time being 0.6 seconds. Since both solution
algorithms converge to within the same numerical accuracy, the discrepancies between
the smile-fitting and at-the-money implied values are due to the volatility surface. It was
noted in Section 2 that the convex shape of the optimal exercise boundary for the vanilla
American put problem was useful in increasing the efficiency of the tridiagonal solver.
Figure 4.3 highlights the reason why the volatility smile fitting option value is different
from its vanilla counterpart by illustrating the shifted nonconvex shape of its optimal
exercise boundary. When we take account of local volatility, the exercise boundary is
no longer a convex function of the asset price, but is shifted horizontally by changes in
local volatility. Although this poses no problem for the accuracy of the pricing algorithm,
it does radically affect the realized option price. Similar results for the S&P 500 index
options are given in Dempster and Richards (1999).

4.6. Pricing American Discretely Sampled Asian Fixed-Strike Options

Finally we price Asian fixed-strike options fitting the smile. An Asian rate (or fixed-
strike) option has the vanilla payoff with the asset price at exercise replaced by the
average. For example, the Asian rate put has payoff ψ(S, A) := max(K − A, 0) for the
fixed strike K.
Table 4.5 contains option values for the American Asian put option with fixed-strike
equal to the initial asset price. All results correspond to options of 0.211 year maturity
with the risk-free rate assumed constant at 10 percent. In the table the “Implied value”
columns refer to the value found using the LP approach with constant volatility set to the
174 m. a. h. dempster and d. g. richards

Figure 4.3. American put optimal exercise boundary for the smile fit compared with the
1-factor LP method using ATM BS implied volatility. Parameters: M = 1000, I = 3000,
S0 = 3129.5, K = 3125, σBS = 15.89%. Dotted lines show market European call
maturities.

Black–Scholes at-the-money implied volatility for the options in question. The columns
labeled “Smile value” illustrate the results obtained by fitting the volatility smile and
term structure.
As can be seen, the option price fitting the volatility smile is significantly different
from the constant volatility price. Unlike the vanilla American put the American fixed-
strike Asian option results on the FTSE 100 index (Table 4.5) show a mixed effect. For
a low arithmetic average sampling rate the smile-fitting value is less than the Black–
Scholes implied value, but the converse relationship holds for higher sampling rates.
This effect is likely due to the time-to-maturity variability of implied volatility at strike
2975 depicted in Figure 4.1.

5. CONCLUSION
After surveying earlier results, we applied a fast accurate linear programming valua-
tion algorithm to pricing American exotic options fitting the volatility smile implied by
the market prices of vanilla European call options. We have demonstrated first that the
basic Crank–Nicolson finite difference methods have low discretization error and that
the quoted vanilla options are accurately priced by the fitted local volatility surface.
Subsequently we have seen that, due to local volatility effects on the computed optimal
exercise boundary, prices of American options fitted to the smile differ significantly from
 pricing american options fitting the smile 175

Table 4.4
American Put Option Valuation Results Fitting the FTSE 100 Volatility Smile

April May June

Strike
price Smile fit LP(σatm ) Smile fit LP(σatm ) Smile fit LP(σatm )
2975 9.31 7.12 18.54 15.01 29.86 24.62
3025 15.72 12.77 28.66 24.05 40.69 34.13
3075 26.03 21.86 42.92 37.11 55.40 47.43
3125 46.40 41.77 62.28 54.65 74.50 64.45
3175 76.53 70.38 88.59 78.01 98.70 85.33
3225 116.51 106.80 122.79 107.71 129.19 110.15
3275 166.50 150.90 165.84 143.13 167.70 139.95
3325 215.50 197.59 215.50 183.13 215.50 173.86
M = 3200, I = 200. Solution times approximately 0.6 seconds. LP(σatm ) is the
value calculated from the tridiagonal solver with discretization M = 10000, I =
10000.

those with constant volatilities. Finally, we have seen similar effects for American exotic
options, as represented by discretely sampled fixed-strike Asian options.
Current research extends the testing of these methods to lookbacks and barriers, includ-
ing both digitals and knock-in and knock-out features for Asians and lookbacks. An inter-
esting area of related research involves the Kalman filtering of local volatility surfaces—
as for example computed in this paper—from one market epoch (day) to the next in order
to achieve better long-run hedging. Another line of our current research with PDE-based
valuation methods concerns wavelet basis techniques for discontinuous option payoffs,
including barriers (Dempster et al. 1999), and high-dimensional Bermudan and American
fixed income derivatives (see Dempster and Hutton (1997b)).

Table 4.5
Discretely Sampled American Fixed-Strike Asian Put Option Results Fitting the
FTSE 100 Smile

2 Samples 12 Samples
Implied Smile Implied Smile
M I J value value M I J value value
400 200 200 0.360 0.295 270 200 200 3.134 3.710
800 200 200 0.360 0.296 540 200 200 3.135 3.759
800 400 400 0.353 0.278 540 400 400 3.127 3.737
Parameters: K = 2975, S0 = 3129.5 and T = 0.211 years. Solution time for
M = 200, I = J = 200 is approximately 18 seconds.
176 m. a. h. dempster and d. g. richards

references

Andersen, L. B. G., and R. Brotherton-Ratcliffe (1997): The Equity Option Volatility

Smile: An Implicit Finite-Difference Approach, J. Computat. Finance 1, 5–37.
Baiocchi, C., and A. Capelo (1984): Variational and Quasivariational Inequalities. New
York: Wiley.
Barraquand, J., and T. Pudet (1996): Pricing of American Path-Dependent Contingent
Claims, Math. Finance 6, 17–51.
Black, F., and M. Scholes (1973): The Pricing of Options and Corporate Liabilities, J. Polit-
ical Econ. 81, 637–659.
Bodurtha, J., Jr. (1997): A Linearization-Based Solution to the Ill-Posed Local Volatility
Estimation Problem, Working paper, Georgetown University.
Borwein, J. M., and M. A. H. Dempster (1989): The Linear Order Complementarity Prob-
lem, Math. Operations Res. 14, 534–558.
Bouchouev, I., and V. Isakov (1999): Uniqueness, Stability and Numerical Methods for the
Inverse Problem that Arises in Financial Markets, Inverse Prob. 15, R95–R116.
Breeden, D., and R. Litzenberger (1978): Prices of State-Contingent Claims Implicit in
Option Prices, J. Business 51, 621–651.
Clarke, N., and K. Parrott (1998): Multigrid for American Option Pricing with Stochastic
Volatility, Working paper, Mathematical Sciences, University of Greenwich.
Coleman, T. F., Y. Li, and A. Verma (1999): Reconstructing the Unknown Local Volatility
Function, J. Computat. Finance 2, 77–102.
Cryer, C. W. (1971): The Solution of a Quadratic Programme Using Systematic Overrelax-
ation, SIAM J. Control Optim. 9, 385–392.
Cryer, C. W., and M. A. H. Dempster (1980): Equivalence of Linear Complementarity
Problems and Linear Programs in Vector Lattice Hilbert Spaces, SIAM J. Control Optim.
18, 76–90.
Dempster, M. A. H., A. Eswaran, and D. G. Richards (1999): Wavelets and Their Applica-
tion in Valuing Derivatives by PDE Methods, Working paper, Centre for Financial Research,
University of Cambridge.
Dempster, M. A. H., and J. P. Hutton (1997a): Fast Numerical Valuation of American,
Exotic and Complex Options, Appl. Math. Finance 4, 1–20.
Dempster, M. A. H., and J. P. Hutton (1997b): Numerical Valuation of Cross-Currency
Swaps and Swaptions; in Mathematics of Derivative Securities, M. A. H. Dempster and S.
R. Pliska, eds. Cambridge: Cambridge University Press, 473–503.
Dempster, M. A. H., and J. P. Hutton (1999): Pricing American Stock Options by Linear
Programming, Math. Finance 9, 229–254.
Dempster, M. A. H., J. P. Hutton, and D. G. Richards (1998): LP Valuation of Exotic
American Options Exploiting Structure, J. Computat. Finance 2, 61–84.
Dempster, M. A. H., and D. G. Richards (1999): Pricing Exotic American Options Fit-
ting the Volatility Smile, Working paper 17/99, Judge Institute of Management Studies,
University of Cambridge, submitted.
Derman, E., and I. Kani (1998): Stochastic Implied Trees: Arbitrage Pricing with Stochastic
Term and Strike Structure of Volatility, Inter. J. Theoret. Appl. Finance 1, 61–110.
Duan, J.-C. (1995): The Garch Option Pricing Model, J. Finance 5, 13–32.
 pricing american options fitting the smile 177

Dupire, B. (1997): Pricing and Hedging with Smiles; in Mathematics of Derivative Securities,
M. A. H. Dempster and S. R. Pliska, eds. Cambridge: Cambridge University Press, 103–
111.
Harrison, J. M., and D. Kreps (1979): Martingales and Arbitrage in Multi-Period Security
Markets, J. Econ. Theory 20, 381–408.
Harrison, J. M., and S. R. Pliska (1981): Martingales and Stochastic Integrals in the Theory
of Continuous Trading. Stoch. Process. Appl. 11, 215–260.
Hull, J., and A. White (1987): The Pricing of Options on Assets with Stochastic Volatilities,
J. Financial Quant. Anal. 3, 281–300.
Hutton, J. P. (1995): Fast Pricing of Derivative Securities, Ph.D. thesis, Department of
Mathematics, University of Essex.
Jackwerth, J. C., and M. Rubinstein (1996): Recovering Probability Distributions from
Option Prices, J. Finance 51, 62–74.
Jaillet, P., D. Lamberton, and B. Lapeyre (1990): Variational Inequalities and the Pricing
of American Options, Acta Appl. Math. 21, 263–289.
Lagnado, R., and S. Osher (1997): A Technique for Calibrating Derivative Security Pricing
Models: Numerical Solution of an Inverse Problem, J. Computat. Finance 1, 13–25.
Richards, D. G. (1999): Pricing American Exotic Options, Ph.D. thesis, Judge Institute of
Management Studies, University of Cambridge.
Rubinstein, M. (1994): Implied Binomial Trees, J. Finance 49, 771–818.
Wilmott, P., S. Howison, and J. Dewynne (1993): Option Pricing: Mathematical Models
and Computation. Oxford: Oxford Financial Press.