Option Pricing by Transform Methods:

Extensions, Unification, and Error Control

Roger W. Lee ∗

Stanford University, Department of Mathematics

NYU, Courant Institute of Mathematical Sciences

Journal of Computational Finance (2004), 7(3):51–86

This version: February 2, 2005

Abstract

We extend and unify Fourier-analytic methods for pricing a wide class of options on any underlying
state variable whose characteristic function is known. In this general setting, we bound the numerical
pricing error of discretized transform computations, such as DFT/FFT. These bounds enable algorithms
to select efficient quadrature parameters and to price with guaranteed numerical accuracy.

∗ This work was partially supported by an NSF Mathematical Sciences Postdoctoral Fellowship. I thank Marco Avellaneda,
Peter Carr, Amir Dembo, Darrell Duffie, Esben Hedegaard, George Papanicolaou, Liuren Wu, and seminar participants at Banc of
America Securities, the Courant Institute, and Stanford University. Email: rogerlee@[Link].

1
1 Introduction

In a large and growing family of financial models, explicit formulas exist for the characteristic functions of
the state variables. Given any such characteristic function, our project is to compute, efficiently and accu-
rately, the prices of a wide class of options on those underlying state variables. Fourier-analytic solutions to
various forms of this problem have appeared in the finance literature. They express option prices in terms of
Fourier-inversion integrals, which are in practice evaluated numerically.
This paper extends and unifies those ideas. In a general setting, we bound the error in the numeri-
cal evaluation of these integrals as N-point sums, of the kind that may be computed as a discrete Fourier
transform (DFT) by schemes including the fast Fourier transform (FFT). Then we show how these bounds
lead to algorithms that make efficient choices of quadrature parameters and compute prices with guaranteed
numerical accuracy.

1.1 Outline

This paper generalizes Carr-Madan (1999); unifies it with extensions of the relevant elements of Duffie-
Pan-Singleton (2000), Lewis (2001), and Bakshi-Madan (2000); and develops error bounds and error mini-
mization strategies. Carr-Madan’s underlying random variable X is the logarithm of a terminal stock price,
and their objective is to compute the call price, as a function of log-strike. In terms of the characteristic
function of X, they calculate analytically the Fourier transform of the call price function, damped to enforce
integrability. Inverting this Fourier transform by FFT and then undamping, they recover simultaneously the
call prices at many strikes.
We begin by setting forth the option pricing problem and defining the options to be priced. Our scope
includes not only vanilla calls on variables exponential in a single state variable, but also three other classes
of payoffs. These extended payoff classes contain all of the derivative structures treated in Duffie-Pan-
Singleton (2000), and in particular they allow payoffs dependent on multidimensional state variables.
Next, we derive upper bounds on option prices, intended for use at extreme strikes. These bounds will
become relevant to discretization errors in transform-pricing of options at all strikes.
In Section 4 we extend, to all four payoff classes, Carr-Madan’s analytic calculation of Fourier trans-
forms, as well as their inversion formula recovering the option price. Also, by taking as given the Bakshi-
Madan (2000) discounted characteristic function, we extend Carr-Madan to allow stochastic interest rates.
As Lewis (2001) observes, transform representations of option prices may be interpreted as contour
integrals in the complex plane; shifts of the contours generate alternative pricing formulas. Applying this

2
idea, we prove a unified pricing formula encompassing not just our original four but also ten complementary
formulas, including as special cases some well-known transform formulas.
These formulas involve integrals over (a translate of) the real line, so approximation by an N-point sum
is subject to two forms of error: sampling error because the integrand is evaluated numerically only at the
grid points, and truncation error because the upper limit of numeric summation is finite. We then establish
bounds for both kinds of error, in all four payoff classes.
Section 7 addresses strategic issues in error bound minimization. From an error-management perspec-
tive, we apply our bounds analysis to argue in favor of the Carr-Madan one-integral approach to call pricing,
and against the traditional two-integral approach. Then we make recommendations for choosing among our
five one-integral call formulas. For choosing quadrature parameters, we offer a simple algorithm as a robust
alternative to the specific constant parameters suggested in Carr-Madan.
The first appendix facilitates truncation error calculations by providing bounds on the decay of charac-
teristic functions in two prominent models. The second appendix gives sampling error bounds, for subcases
deferred from the main text. The third appendix deals with specific DFT/FFT implementation issues.

1.2 Guiding Principles

Wherever possible, we observe the following principles.

First, we take as primitive the discounted characteristic function. From there, our analysis proceeds to
the computation of option prices. We do not derive any characteristic functions; other papers have already
taken the responsibility of finding characteristic functions given, for example, SDE or generating triplet
specifications of the underlying financial dynamics; and indeed others take the characteristic function as
the specification of the underlying dynamics. Duplication of research effort will be reduced, one hopes,
by the emerging division of labor between, on one hand, those projects that specify or derive characteristic
functions; and, on the other hand, projects such as this one, which derive option pricing formulas, given
arbitrary characteristic functions.
Second, we strive to maintain generality. We do not assume that the underlying state variable is, say,
a jump-diffusion or Lévy process. We do not assume that its probability distribution has a density. Time
and the state space may be continuous or discrete. The state variables may be one-dimensional or multi-
dimensional. Interest rates and dividends may be deterministic or stochastic. As long as the discounted
characteristic function for such dynamics is known, option prices are computable. Technical restrictions do
apply, which brings us to the next point.
Third, we formulate our technical conditions with the view that they should facilitate the design of

3
provably robust pricing algorithms. So we place a premium on expressing assumptions in a complete,
concise, rigorous, and readily testable way.

2 The Option Pricing Problem

Working in a filtered probability space (Ω, P∗ , {Ft }), we intend to calculate numerically the time-0 price C0
of an option paying at time T the FT -measurable random variable CT .
Let rt be the interest rate process, possibly stochastic.
Let Mt := exp( 0t rs ds) be the time-t value of a money market account.
R

Let Bt be the time-t value of a discount bond maturing at T .

2.1 Numeraires and Martingale Measures

Assuming that the prices (of C, M, B, and any other assets under consideration) admit no arbitrage, there
must exist a risk-neutral probability measure P under which asset prices, discounted by M, are martingales.
See Harrison and Kreps (1979) or Delbaen and Schachermayer (1994) for technical definitions of “admit no
arbitrage” that make this statement true. Let E denote expectation with respect to P. Then the option price
and bond price satisfy

C0 = E[MT−1CT ]

B0 = E[MT−1 ].

The positive price process M is an example of a numeraire. For any numeraire N there exists a probability
measure PN , said to be risk neutral with respect to N, meaning that the Nt -discounted price of any asset is a
PN -martingale; see El Karoui, Geman, and Rochet (1995). The change of measure from P to PN is given by

dPN NT /N0
= .
dP FT MT /M0

When the numeraire is chosen to be the price Bt of a T -maturity discount bond, the risk-neutral measure
PB is known as the T -forward measure. Let us write E for expectation with respect to PB . The option price
satisfies, therefore,
C0 = B0 ECT .

In the case of deterministic interest rates, the forward measure is identical to the usual risk-neutral measure.
In our setting, however, interest rates may be stochastic, and the measures are not necessarily identical; the
forward measure has the advantage of discounting outside the expectation.

4
2.2 Options

Let the state variable X be an FT -measurable random variable with values in Rn . For a payoff function
G : Rn × R → R, define
CG (k) := B0 E(G(X, k)),

which is the time-0 price of an option on X, paying G(X, k) at time T . The trigger k is some contract
variable, such as a strike, or the logarithm of a strike.
Our goal is accurate numerical computation of CG (k) for these cases of G:

G1 (x, k) := (exp(x) − exp(k))+ b0 := 1, b1 := 1, x ∈ R

G2 (x, k) := (x − k)+ b0 := 1, b1 := 0, x ∈ R

G3 (x, k) := exp(b1 · x)I(b0 · x > k) x ∈ Rn

G4 (x, k) := (b2 · x) exp(b1 · x)I(b0 · x > k) x ∈ Rn

where I is the indicator function, so I(b0 · x > k) equals 1 if b0 · x > k, but 0 otherwise. In payoffs G3 and
G4 , the b0 , b1 , b2 ∈ Rn are arbitrary constants. When it is clear what payoff(s) is/are under discussion, we
may suppress the subscript of G or C.
We choose these four functional forms because they include a wide family of payoffs of practical interest.
For example, with payoff G1 , if one chooses X to be bond yield, or the logarithm of a stock price or FX rate,
then one obtains a call on a stock, bond, or currency. With payoff G2 , if one chooses X to be an interest rate,
or a time-averaged interest rate, then one obtains respectively a European or an Asian option on an interest
rate.
Our G3 and G4 are the payoff classes treated in Duffie-Pan-Singleton (2000). With payoff G3 , if one
chooses b0 and b1 appropriately, then one can obtain asset-or-nothing, binary, equity-linked FX, and two-
asset exchange/maximum options, all on the exponentials of components of X, which could be stock price
logarithms or bond yields or FX-rate logarithms. With payoff G4 , if one chooses b1 = 0 and b0 and b2
appropriately, then one can obtain basket or spread options on the components of X, which could be interest
rates or their time-averages, for example.

3 Upper Bounds on Option Prices at Extreme Strikes

For practical use in bounding numerical transform-inversion errors, it is important that CG be dominated by
an expression that is easily evaluated in terms of the characteristic function of X.

5
 For each G = G1 , . . . , G4 , we give two bounds; both bounds are valid for all k, but the first is intended for
use with large positive k, whereas the second is intended for use with large negative k. The usual conventions
about ∞ are in force, so each of Theorems 3.1-3.4 holds automatically if the expectation on the right-hand
side is infinite.
The first of these four results is nearly identical to a bound obtained in Broadie-Cvitanic-Soner (1998).
The differences, though minor, make it appropriate to present briefly a full proof.

Theorem 3.1. For any p > 0,

 p
B0 E exp((p + 1)X) p
CG1 (k) 6 and CG1 (k) 6 B0 E exp(X).
(p + 1) exp(pk) p+1

Proof. For all s > 0 we have p

s p+1

k p
s−e 6
(p + 1) exp(pk) p+1
because the left-hand and right-hand sides, as functions of s, have equal values and first derivatives at
s = (p + 1) exp(k)/p, but the right-hand side has everywhere a positive second derivative. Moreover, since
the right-hand side is positive, the left-hand side can be improved to (s − exp(k))+ .
Now substitute s = exp(X), take expectations, and multiply by B0 to obtain the first bound. The second
bound is obvious.

Remark 3.1. Therefore, if ST is a nonnegative random variable with ESTp+1 < ∞ for some p > 0, then calls
on ST must have prices that decay as O(K −p ) for strikes K → ∞.
A corresponding fact for puts follows from Theorem 6.4: if EST−q < ∞ for some q > 0, then puts on on
ST must have prices that decay as O(K q+1 ) for strikes K → 0.
Lee (2003) uses these bounds to derive an explicit “moment formula” for the growth of implied volatility
at extreme strikes.

Theorem 3.2. For any p > 0,

B0 E exp(pX)
CG2 (k) 6 .
p exp(pk + 1)
For any q > 0,
B0 E exp(−qX)
CG2 (k) 6 B0 (EX − k) + .
q exp(1 − qk)
Proof. For all x ∈ R we have
exp(px)
x−k 6
p exp(pk + 1)

6
because the left-hand and right-hand sides, as functions of x, have equal values and first derivatives at
x = k + 1/p, but the right-hand side has everywhere a positive second derivative. Substitute X for x, take
expectations, and multiply by B0 to obtain the first bound.
A similar argument shows that for all x,

exp(−qx)
(x − k)+ = x − k + (k − x)+ 6 x − k + ,
q exp(1 − qk)

which implies the second bound.

Theorem 3.3. For any p > 0,

B0 E exp((pb0 + b1 ) · X)
CG3 (k) 6 and CG3 (k) 6 B0 E exp(b1 · X).
exp(pk)

Proof. For all x ∈ Rn we have

exp(pb0 · x)
I(b0 · x > k) 6 ,
exp(pk)
which implies the first bound. The second bound is obvious.

Theorem 3.4. For any p0 > 0 and p2 > 0,

B0 E exp((p0 b0 + p2 b2 + b1 ) · X)
CG4 (k) 6 .
p2 exp(p0 k + 1)

For any q0 > 0 and q2 > 0,

B0 E exp((−q0 b0 − q2 b2 + b1 ) · X)
CG4 (k) 6 B0 E((b2 · X) exp(b1 · X)) + .
q2 exp(−q0 k + 1)

Proof. For all x ∈ Rn we have

exp(p2 b2 · x) exp(p0 b0 · x)
(b2 · x)I(b0 · x > k) 6
p2 e exp(p0 k)

and
exp(−q2 b2 · x) exp(−q0 b0 · x)
(b2 · x)I(b0 · x > k) = (b2 · x)(1 − I(b0 · x 6 k)) 6 b2 · x + ,
q2 e exp(−q0 k)
implying the two bounds.

Remark 3.2. To bound |CG4 |, apply Theorem 3.4 to (b0 , b1 , b2 ) and (b0 , b1 , −b2 ), and take the larger of the
two bounds. To bound C|G4 | , take the sum of those two bounds, because |b2 · X| is the sum of (b2 · X)+ and
(−b2 · X)+ .

7
4 From Characteristic Functions to Option Prices

Our starting point is the discounted characteristic function f of the state variable X. Unlike the usual charac-
teristic functions of probability theory, the definition of f includes a discount factor inside the expectation,
which is essential for pricing under stochastic interest rates.
We produce formulas for prices of each of the four option classes, by expressing option price transforms
in terms of f , and then inverting the transforms.

4.1 The Discounted Characteristic Function

Let X be an Rn -valued random variable. Let AX denote the interior of the set

v ∈ Rn : Eev·X < ∞ .


The complex vectors whose negated imaginary parts are in AX form a “strip” or “tube”

ΛX := {ζ ∈ Cn : −Im(ζ ) ∈ AX }.

Adopting the terminology suggested in Bakshi-Madan (2000), define the discounted characteristic function
rt dt), to be the function f : ΛX → C where
RT
of X, with respect to a discount factor exp(− 0

RT
f (ζ ) := E e− rt dt iζ ·X


0 e .

Note that the expectation is with respect to P, but f is also related to the forward measure PB , because

f (ζ )/ f (0) = Eeiζ ·X ,

which is (for ζ restricted to Rn ) the usual characteristic function of X with respect to PB .

Theorem 4.1. The discounted characteristic function f is well-defined and analytic in ΛX , which is a convex
set. Partial derivatives of f may be taken through the expectation.

Proof. This follows from Zemanian (1966), Theorems 4 and 5.

In certain models, one can derive the discounted characteristic function from an SDE specification of
state variable dynamics. For example, affine jump-diffusion specifications give rise to tractable characteristic
functions, as shown in Heston (1993), Bates (1996, 2000), Bakshi-Cao-Chen (1997), Bakshi-Madan (2000),
Duffie-Pan-Singleton (2000), and Chacko-Das (2002). Outside of that family, Lewis (2000), Schöbel-Zhu
(1999), and Zhu (2000) obtain characteristic functions also for non-affine volatility and interest rate models.

8
 In other models, the state variables follow Lévy processes, and one can derive the characteristic func-
tion from a specification of the generating triplet, or directly take the characteristic function to define the
dynamics. Examples include the Finite Moment Log Stable model in Carr-Wu (2003), the Normal In-
verse Gaussian model in Barndorff-Nielsen (1998), the Generalized Hyperbolic model in Eberlein-Prause
(2002), the Variance Gamma model in Madan-Carr-Chang (1998), and the CGMY and KoBoL models in
Carr-Geman-Madan-Yor (2002) and Boyarchenko-Levendorskiǐ (2002). Extensions of Lévy process mod-
els which introduce stochastic time changes also have, in certain cases, explicit solutions for characteristic
functions; see Barndorff-Nielsen/Nicolato/Shephard (2002), Carr-Wu (2002), and Carr-Geman-Madan-Yor
(2003).
Appendix A gives details of the characteristic functions in two models – one in the affine class, and one
in the Lévy class.
Note that if discounted characteristic functions are available not just for state variables but also for path
functionals of the state variables, then our pricing and error control results will apply not just to European
options, but also to path-dependent options. For example, in affine models, the availability of characteristic
functions for time-averages enables us to price Asian options (on the state variables, not on their expo-
nentials). Such availability is, however, the exception rather than the rule. Transform-based pricing of
exotic options is feasible even without a readily computable characteristic function for the path-dependent
quantity, provided that the dynamics are simple enough (under geometric Brownian motion, for example,
see Fu-Madan-Wang (1999) or Carr-Schröder (2003) for Asian options, and Geman-Yor (1996) or Pelsser
(2000) for barrier options); but this falls outside the scope of our pricing and error control results, which
assume the availability of the characteristic function.

4.2 Fourier Transform of the Damped Option Price

The usual Fourier transform of CG itself does not exist, because CG (k) does not decay as k → −∞.
Following Carr-Madan, then, for each damping constant α > 0, we define the damped option price
function cα,G : R → R by
cα,G (k) := exp(αk)CG (k).

We will show that the damped option price cα,G does have a Fourier transform ĉα,G : R → C, well-defined
by
Z ∞
ĉα,G (u) := eiuk cα,G (k)dk,
−∞

provided that α is chosen appropriately.

9
Theorem 4.2. Assume that G satisfies b1 ∈ AX . Then there exists α > 0 with αb0 + b1 ∈ AX . For any such
α the Fourier transform ĉα,G of cα,G exists and

f (u − (α + 1)i) f (u − αi)
ĉα,G1 (u) = ĉα,G2 (u) =
α2 + α − u2 + i(2α + 1)u (α + iu)2

f (ub0 − (αb0 + b1 )i) −ib2 · ∇ f (ub0 − (αb0 + b1 )i)

ĉα,G3 (u) = ĉα,G4 (u) = .
α + iu α + iu

Proof. There exists p > α such that pb0 + b1 ∈ AX . So Theorem 3.1, 3.2, 3.3, or 3.4 implies that c(k)
decays exponentially for |k| → ∞. Also c(k) is bounded. Therefore c(k) is L1 and has a Fourier transform;
moreover, the use of Fubini in the following computation of ĉ is justified:
Z ∞ Z ∞ Z ∞
ĉα,G (u) := eiuk cα,G (k)dk = eiuk eαk B0 EG(X, k)dk = f (0)E G(X, k)e(α+iu)k dk.
−∞ −∞ −∞

Evaluating the integral,

f (0)Ee(α+1+iu)X f (0)Ee(α+iu)X
ĉα,G1 (u) = ĉα,G2 (u) =
α 2 + α − u2 + i(2α + 1)u (α + iu)2

f (0)E(eb1 ·X e(α+iu)b0 ·X ) f (0)E((b2 · X)eb1 ·X e(α+iu)b0 ·X )

ĉα,G3 (u) = ĉα,G4 (u) = .
α + iu α + iu

The result follows because ub0 − (αb0 + b1 )i ∈ ΛX .

4.3 Fourier Inversion

Option prices may be recovered via Fourier inversion.

Theorem 4.3. Suppose G and α satisfy the hypotheses of Theorem 4.2.

In cases G = G1 , G2 , the option price is given by
Z ∞ Z ∞
e−αk −iuk e−αk
Re e−iuk ĉα,G (u) du.
 
CG (k) = e ĉα,G (u)du = (4.1)
2π −∞ π 0

In cases G = G3 , G4 , define the average of left and right limits C̄(k) := [C(k + 0) −C(k − 0)]/2. Then
Z ∞
e−αk e−αk
Z R
−iuk
Re e−iuk ĉα,G (u) du,
 
C̄G (k) = lim e ĉα,G (u)du = (4.2)
2π R→∞ −R π 0

which can be strengthened to (4.1) if ĉ is L1 .

10
Proof. In all cases, the damped option price c(k) is L1 , as argued in Theorem 4.2.
In cases G = G1 and G = G2 , the damped price c(k) is continuous, by the dominated convergence
theorem; and the transform is L1 , because |ĉ(u)| 6 | f (−(α + b1 )i)|/(u2 + α 2 ). Therefore the usual Fourier
inversion recovers c; see, for example, Champeney (1987) Theorem 8.2. Undamping with a factor of e−αk
yields (4.1).
In cases G = G3 and G = G4 , the damped price c(k) is locally of bounded variation, because C(k) is the
difference of two monotonic functions, exp(αk) is monotonic, and both C(k) and exp(αk) are bounded on
any finite interval. By, say, Champeney (1987) Theorem 8.12, we have (4.2).

5 The Pricing Formula for General α

Transform representations of option prices can be viewed as contour integrals in the complex plane. Shifting
the contour across a pole of the integrand changes the value of the integral, a technique which Lewis (2001)
exploits, as will we.
Lewis differs from our approach in that he derives formulas for the transforms of option prices with
respect to the spot variable X0 ; whereas we, like Carr-Madan and Duffie-Pan-Singleton, transform with
respect to the trigger variable k. His assumptions require that the option be written on the exponential of a
variable XT where the distribution of XT − X0 is not permitted to depend on X0 . Our formulas are not subject
to this restriction and apply to a wider class of underlying state variables X, including those exhibiting
mean-reversion.
One can modify the formulas of Lewis for non-independent-increments. However, the resulting formulas
in that case do not allow the direct application of FFT to calibrate parameters to the prices of options at
multiple strikes. For that purpose one needs transform-in-strike formulas, which we now derive.
Specifically, let Γ := ΓX,G := {z ∈ C : −Im(z)b0 + b1 ∈ AX }, and define CˆG : ΓX,G → C by
f (z − i) − f (z)
CˆG1 (z) := CˆG2 (z) :=
iz − z2 z2

f (b0 z − b1 i) −b2 · ∇ f (b0 z − b1 i)

CˆG3 (z) := CˆG4 (z) := . (5.1)
iz z
Theorem 4.2 proves that for positive α with αb0 + b1 ∈ AX , we have

ĉα,G (u) = CˆG (u − αi),

and hence, for z ∈ Γ such that −Im(z) > 0,

Z ∞
CˆG (z) = eizkCG (k)dk. (5.2)
−∞

11
Thus, in this region, CˆG (z) is the complex Fourier transform of the unmodified option price CG (k). Equiv-
alently (modulo rotation by a factor of i), CˆG (z) is the bilateral Laplace transform of CG (k). Rewriting the
conclusion of Theorem 4.3 shows that CG may be inverted by integrating along the contour Im(z) = −α in
the complex plane:

1 ∞−αi ˆ 1 ∞−αi
Z Z
CG (k) = CG (z)e−ikz dz = Re[CˆG (z)e−ikz ]dz (G = G1 , G2 )
2π −∞−αi π 0−αi
(5.3)
1 ∞−αi
Z R−αi
1
Z
C̄G (k) = lim CˆG (z)e−ikz dz = Re[CˆG (z)e−ikz ]dz (G = G3 , G4 ),
2π R→∞ −R−αi π 0−αi
again assuming positive α with αb0 + b1 ∈ AX .
For negative α, the transform ĉα,G does not exist (for G = G1 , . . . , G4 ); likewise, the integral in (5.2)
does not exist for −Im(z) < 0. Nonetheless, the definitions (5.1) do make sense, and the integrals in (5.3)
do exist for α < 0, but they do not recover C̄G , because the integration path has shifted across the pole z = 0;
instead they recover C̄G less the contribution of the residue of CˆG at z = 0. In each case G = G1 , . . . G4 , this
generates one additional pricing formula. In case G = G1 , it generates a second additional formula, because
CˆG1 has a second pole at z = i.
For zero α (and for α = −1 in case G = G1 ), the final integrals in (5.3) are again well-defined, but now
the integration contour passes through a pole, and the contribution from the residue is cut in half. (The only
exception is in the case G = G2 which has a double pole at z = 0; this case calls for introducing into the
integrand a term that tames the singularity, without affecting the value of the integral.) This generates two
additional pricing formulas for payoff G1 , and one additional formula for the other payoffs.
Theorem 5.1 makes this discussion precise. Note that by taking α = 0 in cases G = G3 and G = G4 ,
we recover both of Duffie-Pan-Singleton’s (2000, Prop 2 and Eqn 3.8) pricing formulas. Taking α > 0 in
case G = G1 recovers Carr-Madan’s damped-call pricing formula. Taking α = 0 in two instances of case
G = G3 recovers the traditional two-integral call-pricing formulas, which we discuss further in Section 7.1.
Our central pricing result is as follows.

Theorem 5.1. Assume that b1 ∈ AX . Let α be any real number such that αb0 + b1 ∈ AX . Then in all cases
except (G = G2 ; α = 0),
Z ∞−αi
1
C̄G (k) = Rα,G + Re[CˆG (z)e−ikz ]dz (5.4)
π 0−αi

12
where



 f (−i) − ek f (0) α < −1

 
f (−i) − ek f (0)/2 α = −1 −i f 0 (0) − k f (0) α <0

 

 

 
Rα,G1 := f (−i) −1 < α < 0 Rα,G2 := (−i f 0 (0) − k f (0))/2 α = 0

 

 



 f (−i)/2 α =0  0 α >0


 0 α >0

 


 f (−b1 i) α <0 

 −ib2 · ∇ f (−b1 i) α <0
 
Rα,G3 := f (−b1 i)/2 α = 0 Rα,G4 := −ib2 · ∇ f (−b1 i)/2 α = 0

 

 
 0 α >0  0 α >0

and CˆG is given in (5.1). In cases G = G1 and G = G2 , the C̄G can be replaced by CG .

We will prove simultaneously Theorem 5.1 and the following (G = G2 ; α = 0) theorem.

Theorem 5.2. Assume that b1 ∈ AX . Then we have the (G = G2 ; α = 0) formula

Z ∞
1 1
CG2 (k) = R0,G2 + Re[CˆG2 (z)e−ikz ] + 2 dz. (5.5)
π 0 z

Proofs. For α > 0, see Theorem 4.3.

For α < 0 (except for α = −1 in case G = G1 ): Note that each CˆG is analytic in the strip ΓX,G , except
for a pole at z = 0 (and also z = i in case G = G1 ). The residue theorem applies to any rectangular path
with horizontal segments on Im(z) = −α1 and Im(z) = −α2 , and vertical segments on Re(z) = ±R. Since
the integrals over the vertical segments approach 0 as R → ∞, it follows that shifting a horizontal contour
across the pole changes the value of the integral by 2πi times the residue at that pole. Residue calculation is
straightforward.
For α = 0 (including α = −1 in case G = G1 ): Our proof will be for α = 0; a similar argument proves
the (G = G1 ; α = −1) formula. Define the functions SG1 (z) := SG3 (z) := SG4 (z) := 0 and SG2 (z) := 1/z2 .
On ΓX,G ∩ −ΓX,G the function
 
1 ˆ −izk
h(z) := SG (z) + CG (z)e + CˆG (−z)e izk
2

(modulo a removable singularity in case G2 ) is analytic. Choose ε > 0 such that b1 ± εb0 ∈ AX . Applying

13
Cauchy’s Theorem to the appropriate rectangle, and then using the relevant α 6= 0 results, we have
Z R R−εi R−εi
1 1 1
Z Z
lim h(z)dz = lim h(z)dz = lim h(z) − SG (z)dz
2π R→∞ −R 2π R→∞ −R−εi 2π R→∞ −R−εi
 Z R−εi Z R+εi 
1 −izk −izk
= lim CˆG (z)e dz + CˆG (z)e dz
4π R→∞ −R−εi −R+εi
 
1 R−ε,G
= C̄G (k) + (C̄G (k) − R−ε,G ) = C̄G (k) − ,
2 2

as claimed.
Theorem 5.1’s final assertion is by continuity of CG1 and CG2 .

A single piece of numerical integration code (coupled with the appropriate Rα,G adjustment) can eval-
uate, for example, all five formulas for payoff G1 ; the only difference is the value of α passed into the
procedure. Thus, without writing additional code, one gains the flexibility to choose, say, a negative or zero
α if the integrand should happen to behave better there than it does along positive α. The extent to which
an integrand is “well-behaved” can be quantified by the error bounds that arise from that particular choice
of α. This is the subject of the next section.
Also in the next section we give alternative proofs for many of the formulas in Theorem 5.1. The
contour-shift proof, given above, has the purpose of unifying the various Fourier pricing formulas; but for
the purpose of deriving error bounds, it will be useful to reinterpret the results. For example, our α < 0
bounds will exploit the equivalence between contour shifts and parity relations, such as put/call.

6 Bounds for Sampling and Truncation Errors

The Fourier inversion (5.4) can be approximated discretely via an N-point sum with a grid spacing of ∆ in
the Fourier domain. This quadrature introduces two forms of error (aside from roundoff error): truncation
error because the upper limit of the numeric integration is finite, and sampling error because the integrand
is evaluated numerically only at the grid points. Our bounds will account for both sources of error.
The total error is defined as the absolute difference between the true value
Z ∞
e−αk
Re e−iuk ĉα,G (u) du
 
CG (k) = Rα,G +
π 0

and the discrete approximation given by the N-point sum

−αk ∆
 N−1 
Σ (k) :=
N
ΣN,∆
α,G (k) := Rα,G + e
π
Re ∑ ĉα,G ((n + 1/2)∆)e −i(n+1/2)k∆
. (6.1)
n=0

14
The total error is bounded by the sum of the sampling error and the truncation error

|C − ΣN | 6 |C − Σ∞ | + |Σ∞ − ΣN |,

where Σ∞ is defined as ΣN is, except with an infinite upper limit of summation.

Truncation errors can be bounded by a formula that applies regardless of the sign of α.
Sampling errors, however, will require treatment that depends on the sign of α. Our strategy is based on
Davies (1973), but he restricts attention to the inversion of characteristic functions to recover probabilities,
which is not always appropriate for us; we extend his approach to the inversion of option price transforms.

6.1 Truncation Error

Carr-Madan and Pan each suggest bounds on the tails of certain Fourier inversion integrals, but our specific
need is to bound the tails of the infinite discrete sums that approximate our Fourier integrals.

Theorem 6.1. Assume the hypotheses of Theorem 5.1.

If f is such that ĉα,G decays as a power |ĉα,G (u)| 6 Φ(u)/u1+γ for all u > u0 , where γ ≡ γα,G > 0 and
Φ(u) ≡ Φα,G (u) is decreasing in u, then the truncation error
Φ(N∆)
|Σ∞ (k) − ΣN (k)| 6 , (6.2)
πeαk γ(N∆)γ
provided that N∆ > u0 .
If f is such that ĉα,G decays exponentially |ĉα,G (u)| 6 Φ(u)e−γu for all u > u0 , where γ > 0 and Φ(u) is
decreasing in u, then the truncation error
∆Φ((N + 1/2)∆)
|Σ∞ (k) − ΣN (k)| 6 , (6.3)
2πeαk+γN∆ sinh(γ∆/2)
provided that N∆ > u0 .

Proof. In any case,

∆ ∞
|Σ∞ (k) − ΣN (k)| 6 e−αk ∑ ĉ((n + 1/2)∆) .
π n=N
In the power decay case,
∞ ∞ Z ∞
Φ(n∆) Φ(N∆) dx Φ(N∆)
∆ ∑ ĉ((n + 1/2)∆) 6 ∆ ∑
[(n + 1/2)∆]γ+1
6
∆γ N xγ+1
=
γ(N∆)γ
,
n=N n=N

where the middle step uses the convexity of 1/x2 .

In the exponential decay case,
∞ ∞
Φ((N + 1/2)∆)
∑ ĉ((n + 1/2)∆) 6 Φ((N + 1/2)∆) ∑ e−γ(n+1/2)∆ 6 eγ(N+1/2)∆ − eγ(N−1/2)∆ ,
n=N n=N

as claimed.

15
Remark 6.1. As observed by Carr and Madan in case G = G1 , and as one can verify also in case G = G2 ,
the power decay hypothesis always holds; specifically, let γ = 1 and Φ = | f (−(α + b1 )i)|. However, the
resulting bound is typically poor. For practical purposes it is desirable to improve the power γ or to establish
exponential decay, by factoring in the contribution from the large-u decay of | f (u − (α + b1 )i)|. The nature
of this decay presents itself in the explicit expression for f ; examples appear in Appendix A.

Remark 6.2. The requirement that N∆ > u0 can be dropped, by modifying the right-hand sides of (6.2)
and (6.3) as follows: first replace each N∆ by u0 (thus bounding the n > u0 /∆ terms of truncation error).
Then add a second term, to bound the N 6 n < u0 /∆ terms of the truncation error, by integrating, over the
appropriate finite interval, a bound on ĉ(u) valid for u < u0 , such as the quadratically decaying bound of
Remark 6.1.

6.2 Sampling Error: Positive α

A form of the “aliasing” effect is at work here; by sampling ĉ only at regular discrete intervals, one recovers
not c but rather a periodic function equal to a combination of c and infinitely many shifted copies of c. The
unwanted copies are shifted farther away as ∆ → 0, so the extreme-strike bounds of Section 3 come into
play.
In the main text, our sampling error analysis will focus on the payoff classes of greatest practical interest,
G1 and G3 . For sampling error in cases G2 and G4 , see Appendix B.

Theorem 6.2. Assume that b1 ∈ AX and αb0 + b1 ∈ AX with α > 0.

In case G = G1 we have
p
e−2πα/∆ f (−i) e2π(α−p)/∆ f (−i(p + 1))
 
p
|CG1 − Σ∞
α,G1 | 6 inf + .
p>α: p+1∈AX 1 − e−4πα/∆ (p + 1)e pk (1 − e4π(α−p)/∆ ) p+1

In case G = G3 , assume also that ĉ(u) = O(u−1−γ ) as u → ∞, where γ > 0. Then

 −2πα/∆
f (−ib1 ) e2π(α−p)/∆ f (−i(pb0 + b1 ))

∞ e
|CG3 − Σα,G3 | 6 inf + .
p>α: pb0 +b1 ∈AX 1 − e−4πα/∆ e pk (1 − e4π(α−p)/∆ )

Proof. For any ∆ > 0 and any positive integer j,

Z ∞
1
c(k − 2π j/∆) + c(k + 2π j/∆) = e−iuk ĉ(u) cos(2π ju/∆)du
π −∞
Z ∆
=2 F(u) cos(2π ju/∆)du,
0

where
∞
1
F(u) :=
2π ∑ ĉ(u + n∆)e−i(u+n∆)k .
n=−∞

16
Since F is Lipschitz, the Fourier cosine series may be summed:
∞ h
c(k) + ∑ c(k − 2π j/∆) + c(k + 2π j/∆) cos(2π ju/∆) = F(u)∆.
i

j=1

In particular, taking u = ∆/2, we have

∞
∑
h i
j
c(k) − F(∆/2)∆ = (−1) c(k − 2π j/∆) + c(k + 2π j/∆) .
j=1

Multiplying by exp(−αk) to undamp the call prices,

∞
C(k) − Σ∞ (k) = ∑ (−1) j e−2π jα/∆C(k − 2π j/∆) + e2π jα/∆C(k + 2π j/∆) .
h i

j=1

Therefore, Theorems 3.1-3.4 imply that

∞ p
e2π j(α−p)/∆ Ee(p+1)X
 
p
|CG1 (k) − Σ (k)| 6 B0 ∑ e
∞ −2π jα/∆ X
Ee + (6.4)
j=1 (p + 1)ekp p+1
j odd

and
∞
e2π j(α−p)/∆ Ee(pb0 +b1 )·X
 
|CG3 (k) − Σ (k)| 6 B0 ∑ e
∞ −2π jα/∆
Eeb1 ·X
+ . (6.5)
j=1 ekp
j odd
The results follow from computing the sums.

Note that application of these bounds does not require actual computation of infimums. For example,
in case G = G1 , any choice of p > α with p + 1 ∈ AX produces a valid upper bound, which is subject to
improvement by taking more trial values of p.

6.3 Sampling Error: Negative α

The Theorem 6.2 error bounds assumed that α > 0, and must be modified for α < 0.
We have seen that shifting a Fourier inversion contour across a pole of the integrand changes the value
of the integral. Sampling error bounds will now follow from the fact that the new integral value is the price
of an contract related to the original option via a parity identity, such as put/call.
Specifically, for each G = G1 , . . . , G4 , define one “complementary” payoff G∗ , and for case G1 define a
second complementary payoff G∗∗ by

G∗1 (x, k) := min(exp(x), exp(k)) b0 := 1, b1 := 1, x ∈ R

G∗∗
1 (x, k) := (exp(k) − exp(x))
+
b0 := 1, b1 := 1, x ∈ R

G∗2 (x, k) := (k − x)+ b0 := 1, b1 := 0, x ∈ R

G∗3 (x, k) := exp(b1 · x)I(b0 · x 6 k) x ∈ Rn

G∗4 (x, k) := (b2 · x) exp(b1 · x)I(b0 · x 6 k) x ∈ Rn .

17
These payoffs have the following time-0 values.

Theorem 6.3. Assume that b1 ∈ AX and αb0 + b1 ∈ AX .

In cases G = G2 , G3 , G4 with α < 0; or in case G = G1 with −1 < α < 0, we have
Z ∞−αi
1
C̄G∗ (k) = Re[CˆG (k)e−izk ]dz.
π 0−αi

In case G = G1 with α < −1, this holds after replacing the G∗ with G∗∗ .

Proof. Subtract from each original payoff function G its complementary payoff G∗ ; then take expectations
to verify the parity relation
B0 EG(X, k) − B0 EG∗ (X, k) = R0− ,G (6.6)

and similarly for G∗∗ . Theorem 5.1 now implies the result.
Alternatively, without using Theorem 5.1, one may adapt Theorem 4.2 and compute directly the complex
Fourier transforms of each CG∗ . Inverting as in Theorem 4.3 finishes the proof. Moreover, the negative-α
formulas in Theorem 5.1 would then follow from (6.6).

This equivalence between contour shifts and parity relations allows us to control the negative-α sampling
error by bounding the extreme-strike values of the complementary payoffs. In particular, we state explicitly
the complementary bounds for cases G = G1 and G = G3 .

Theorem 6.4. In case G = G∗1 we have

CG∗1 (k) 6 B0 ek and CG∗1 (k) 6 B0 EeX .

In case G = G∗∗
1 we have, for any q > 0,
q
B0 Ee−qX

q
CG∗∗ (k) 6 e(1+q)k and CG∗∗ (k) 6 B0 ek .
1
1+q 1+q 1

In case G = G∗3 we have, for any q > 0,

B0 Ee(−qb0 +b1 )·X

CG∗3 (k) 6 and CG∗3 (k) 6 B0 Eeb1 ·X .
e−qk

Proof. Adapt the reasoning in Theorems 3.1 and 3.3. We omit the details.

The sampling error bounds now follow.

18
Theorem 6.5. Assume that b1 ∈ AX and αb0 + b1 ∈ AX .
In case G = G1 , for α ∈ (−1, 0):

ek−2π(α+1)/∆ f (0) e2πα/∆ f (−i)

|CG1 − Σ∞
α,G1 | 6 + .
1 − e−4π(α+1)/∆ 1 − e4πα/∆

and for α < −1:

q 
ek+2π(1+α)/∆ f (0) e(1+q)k e−2π(1+q+α)/∆ f (iq)
 
q
|CG1 − Σ∞
α,G1 | 6 inf + .
q>−(α+1): 1 − e4π(1+α)/∆ (1 + q)(1 − e−4π(1+q+α)/∆ ) 1 + q
−q∈AX

In case G = G3 , assume also ĉ(u) = O(u−1−γ ) as u → ∞, where γ > 0. Then for α < 0,
 −2π(α+q)/∆
f (−i(−qb0 + b1 )) e2πα/∆ f (−ib1 )

∞ e
|CG3 − Σα,G3 | 6 inf + .
q>−α: e−qk (1 − e−4π(α+q)/∆ ) 1 − e4πα/∆
−qb0 +b1 ∈AX

Proof. Adapt the reasoning in Theorem 6.2. We omit the details.

6.4 Sampling Error: Zero α

Here we bound the sampling error along contours that pass through a pole. This means α = 0 and, in case
G = G1 , also α = −1.
We present results for cases G = G1 and G = G3 . In each case the option price function can be inter-
preted, after normalization, as a cumulative distribution function, so bounds from the probability literature
apply directly, and we avoid reinvention of the wheel.
Our proof strategy yields, as a by-product, complete alternative proofs of 3 of the 13 formulas in The-
orem 5.1, including the (G = G3 ; α = 0) case, which was Duffie-Pan-Singleton’s (2000, Prop 2) pricing
formula; their proof influenced ours but lacks the highly convenient normalization step.

Theorem 6.6. Assume that b1 ∈ AX and αb0 + b1 ∈ AX .

Then the (α, G) ∈ {(0, G1 ), (−1, G1 ), (0, G3 )} subcases of Theorem 5.1 all hold.
In case G = G1 , the α = −1 and α = 0 sampling errors are bounded by
  q (1+q)k 
∞ f (iq) q e
CG1 − Σ−1,G1 6 max inf , f (−i)e−2π/∆
q>0:−q∈AX 1 + q 1 + q e2πq/∆
  p
∞ f (−i(p + 1)) p
CG1 − Σ0,G1 6 max f (0)e k−2π/∆
, inf .
p>0: p+1∈AX (p + 1)e p(k+2π/∆) p + 1

In case G = G3 assume also that ĉ(u) = O(u−1−γ ) as u → ∞, where γ > 0. Then

 
∞ f (i(qb0 − b1 )) f (−i(pb0 + b1 ))
CG3 − Σ0,G3 6 inf max , .
p>0: pb0 +b1 ∈AX e−q(k−2π/∆) e p(k+2π/∆)
q>0:−qb0 +b1 ∈AX

19
Proof. Case G = G1 , α = −1:
On some probability space (Ω1 , P1 , F ) there exists a real-valued random variable Y with density ϕ(y) :=
e−y E[eX I(X > y)]. It is easy to verify that

CG∗∗
1
(k) = f (0)ek P1 (Y < k),

and that Y has P1 -characteristic function f (u)/[ f (0)(1 − iu)]. By the Gil-Pelaez (1951) formula,

1 1 ∞ f (0)ek 1 ∞+i
     
f (u)/ f (0) −iuk f (z − i) −izk
Z Z
k
CG∗∗ (k) = f (0)e − Re e du = + Re e dz.
1
2 π 0 iu(1 − iu) 2 π 0+i iz(iz + 1)

Davies (1973) now directly implies the sampling error bound

|CG1 − Σ∞
h i
k
−1,G1 | 6 f (0)e max P1 (Y < k − 2π/∆), P1 (Y > k + 2π/∆) .

So for any q > 0,

q
B0 Ee−qX e(1+q)k −2π/∆
  
q
|CG1 − Σ∞
−1,G1 | 6 max ,e X
B0 e ,
1+q 1+q e2πq/∆

as claimed.
Case G = G1 , α = 0:
On some probability space (Ω1 , P1 , F ) there exists a real-valued random variable Y with density ϕ(y) :=
ey PB (X > y) f (0)/ f (−i). It is easy to verify that

CG1 (k) = f (−i)P1 (Y > k),

and that Y has P1 -characteristic function f (u − i)/[ f (−i)(iu + 1)]. By the Gil-Pelaez formula,

1 1 ∞
   
f (z − i)/ f (−i) −izk
Z
CG1 (k) = f (−i) + Re e dz .
2 π 0 iz(iz + 1)

By Davies (1973),

|CG1 − Σ∞
h i
0,G1 | 6 f (−i) max P1 (Y < k − 2π/∆), P1 (Y > k + 2π/∆) .

So for any p > 0,

p
B0 Ee(p+1)X
 
p
|CG1 − Σ∞
0,G1 | 6 max f (0)e k−2π/∆
, .
(p + 1)e p(k+2π/∆) p+1

as claimed.
Case G = G3 , α = 0:

20
 Define the probability measure P3 by dP3 /dP := MT−1 exp(b1 · X)/E(MT−1 exp(b1 · X)). Then
 
−1 b1 ·X dP3
CG3 = E[MT e I(b0 · X > k)] = f (−b1 i)E I(b0 · X > k) = f (−b1 i)P3 (b0 · X > k),
dP

and f (b0 z − b1 i)/ f (−b1 i) is the characteristic function of b0 · X with respect to P3 . So

1 1 ∞
   
f (b0 z − b1 i)/ f (−b1 i) −izk
Z
C̄G3 = f (−b1 i) + Re e dz ,
2 π 0 iz

according to Gil-Pelaez. By Davies (1973),

|CG3 − Σ∞
h i
0,G3 | 6 f (−b 1 i) max P (b
3 0 · X < k − 2π/∆), P (b
3 0 · X > k + 2π/∆) . (6.7)

Therefore, writing E3 for expectation with respect to P3 ,

E3 e−qb0 ·X E3 e pb0 ·X
 
∞
|CG3 − Σ0,G3 | 6 f (−b1 i) max −q(k−2π/∆) , p(k+2π/∆)
e e

for any positive p and q, as claimed.

Remark 6.3. In the case of G = G3 with b0 = b1 = 1 and α = 0, Pan (2002), following Davies, gives
sampling error bounds.
Our G = G3 proof extends Pan, because an alternative way to proceed from (6.7) is (writing H f for the
Hessian matrix of f ):

0 H f (−b1 i)b0 + 2ikb0 · ∇ f (−b1 i) + k f (−b1 i)

(b0 · X − k)2 −b> 2
|CG3 − Σ∞
0,G3 | 6 f (−b1 i)E3 = ,
(2π/∆)2 (2π/∆)2

which improves her bound. Our other incremental contributions here include the complete explicit formu-
lation of technical assumptions and the generality of vectors b0 , b1 ∈ Rn .
This “quadratic” strategy also applies in case G = G1 with α = 0 or −1. However, we prefer sampling
error bounds that go to zero exponentially in −1/∆, rather than quadratically in ∆, so Theorem 6.6 reports
only the exponential results.

6.5 Overall Error Bound: an Example

Consider a call on a stock, under Variance Gamma dynamics, as described in Section A.1. Of the five
formulas in case G = G1 , we choose α > 0. The domain condition αb0 + b1 ∈ AX entails the restriction
α + 1 < a+ , where a+ is defined in (A.1). By (A.2) and Theorem 4.2,

| f (u − (α + 1)i)| exp(−rT + (α + 1)(log S0 + µT ))

|ĉα,G1 (u)| 6 6 .
u2 (u2 νσ 2 /2)T /ν u2

21
By Theorem 6.1, truncation error is bounded by

exp(−rT + (α + 1)(log S0 + µT ))
.
πe (νσ 2 /2)T /ν (1 + 2T /ν)(N∆)1+2T /ν
αk

By Theorem 6.2, sampling error is bounded by

p
e−2πα/∆ f (−i) e2π(α−p)/∆ f (−i(p + 1))

p
+
1 − e−4πα/∆ (p + 1)ekp (1 − e4π(α−p)/∆ ) p+1

for any p > α such that p + 1 < a+ .

Summing the sampling and truncation bounds gives an overall error bound.

7 How To Minimize Error Bounds?

We propose some strategies for choosing α and the quadrature parameters N and ∆ to obtain small error
bounds, given limited computational resources.
Throughout this section, our illustrative problem is to price a vanilla call on a non-dividend-paying
stock whose terminal price is ST = exp(X). According to Theorem 5.1, we may price using any α such that
α + 1 ∈ AX .
Assuming that 0 ∈ AX and 1 ∈ AX (mild assumptions since Ee0·X < ∞ and by no-arbitrage Ee1·X < ∞),
we may write
AX = (−q̄, 1 + p̄),

where p̄ and q̄ are positive. One determines p̄ and q̄ from the explicit expression for the characteristic
function; see Appendix A for examples.
So α may be chosen anywhere in (−q̄ − 1, p̄). This interval comprises five subintervals, corresponding
to the five G1 formulas in Theorem 5.1.
A central question is how to choose from among these five G1 formulas. While the pricing algorithm is
invariant across all five α regimes, the fundamental nature of the bounds differs across the α regimes.
Before addressing this question, let us reject a sixth alternative.

7.1 How Not To Minimize Error Bounds

Instead of pricing a call as a G1 payoff, one can price it as the difference of two G3 payoffs. Indeed, the
latter approach has dominated the literature (exceptions include Carr-Madan and Lewis).

22
 Specifically, most authors have priced the call by decomposing it as long an asset-or-nothing call and
short a binary call. Writing K = exp(k) for the strike,

C(k) = E[MT−1 (ST − K)+ ] = E[MT−1 ST I(ST > K)] − KE[MT−1 I(ST > K)]
   
ST /S0 1/B0
= S0 E I(ST > K) − KB0 E I(ST > K)
MT MT

= S0 PS (ST > K) − KB0 PB (ST > K).

They calculate both pseudo-probabilities via Gil-Pelaez inversions of the PS -characteristic function and the
PB -characteristic function of XT .
In other words,
CG1 (k) = CGb03=b1 =1 (k) − ekCGb03=1,b1 =0 (k),

and each CG3 is evaluated according to the α = 0 formula proved in Theorem 5.1 and again in Theorem 6.6
(the popular proof corresponds to our second proof). This G3 approach to call-pricing has some merits, but
from a computational point of view, it has significant disadvantages.
The generalized Carr-Madan approach of directly pricing G1 has the computational advantage that we
need invert only one Fourier transform, instead of two distinct characteristic functions.
Moreover, our direct G1 error bounds have several advantages over combining two G3 error bounds.
The first is in truncation error control: ĉα,G1 , unlike ĉα,G3 , decays like f divided by the square of u, instead
of u. The second is in sampling error control: the strategy of using exponential functions to dominate
payoff functions produces tighter bounds when the payoff is a call than when the payoff is a binary. Third,
note that summing two G3 error bounds does not take advantage of possible error cancellation between the
two components; in contrast, with only one integral to bound, the direct G1 approach is not subject to this
inefficiency.

7.2 Choice of Quadrature Parameters, Given an α Regime

Deferring once again the discussion of how to choose from among the five α regimes, we address here the
question of how to choose quadrature parameters α, ∆, and N, given an α regime.
In two of the five regimes, the α interval consists of a single point, so the only question is how to choose
N and ∆. For illustrative purposes, suppose we seek quadrature parameters in the α > 0 regime, so all three
quadrature parameters are in question.
The computational burden of the numeric Fourier inversion is determined by the grid point count N.
First suppose that N is fixed and the goal is to find α and ∆ that minimize the total error bound for a given

23
contract k.
Theorem 6.1 gives truncation error bounds; to be concrete, let us suppose the discounted characteristic
function has power decay, as defined there, for u > 0. Theorem 6.2 gives a sampling error bound. Combining
the two, we have the optimization problem

Φα,G1 (N∆)
p
e−2πα/∆ f (−i) e2π(α−p)/∆ f (−i(p + 1))
 
p
min + + . (7.1)
(α,∆,p): ∆>0, πeαk γ(N∆)γ 1 − e−4πα/∆ (p + 1)ekp (1 − e4π(α−p)/∆ ) p + 1
0<α<p< p̄

The choice of (α, ∆, p) can be automated by commonly available simplex optimization algorithms.
To modify (7.1) for models where the desired decay in ĉ(u) is guaranteed only for u > u0 > 0 (see
Theorem 6.1), various options exist. The simplest is to change the ∆ constraint to ∆ > u0 /N, but a more
flexible solution is to allow also ∆ 6 u0 /N, but use the Remark 6.2 bound instead.
Instead of minimizing error bounds for a given computational budget N, an alternative goal would be
to specify a desired error tolerance η, and to find the smallest N for which we can guarantee total error
bounded by η. One strategy here is to choose as trial values for N successively increasing powers of 2. For
each trial N, optimize, over (α, ∆, p), the total error bound in (7.1). Terminate the N loop when the error
bound is smaller than the target η. Pricing can then proceed using the optimal α, ∆, and N.

7.3 Choice of Regime for α

The remaining question is how to choose from among the five regimes.
With unlimited resources, the answer is simple: compute error bounds in all five α regimes, and choose
the one with the smallest error bound. Indeed, even with limited resources, this may prove to be a workable
solution.
However, when the potential benefits of testing all five regimes do not justify the computing or program-
ming effort, it is useful to have rules of thumb regarding which of the five formulas to implement. In any
event, these rules also embody a comparative summary of our various bounds for pricing G1 payoffs.
First consider sampling error. In each of the five regimes, our sampling error bound has ∆ → 0 decay of
order exp(−2πρ/∆), where the “decay rate” ρ is apparent from the relevant Theorem.
Since a greater decay rate yields a better sampling error bound for ∆ sufficiently small, Table 1 suggests
the following sampling error guideline. The “call” regime and “put” regime are intervals with widths p̄
and q̄ respectively; choose α from inside the wider of these two intervals – unless both widths p̄ and q̄ are
smaller than 2. In that case, choose α = −1 or α = 0, according to whether q̄ or p̄ is the larger – unless both
p̄ and q̄ are smaller than 1/2. In that case choose α ∈ (−1, 0).

24
 Table 1: Decay rates of sampling error bounds

α range Value of integral Decay rate ρ of bound Reference

1 R ∞−αi
π 0−αi Re[CG1 (k)e
ˆ −izk ]dz on sampling error

(−q̄ − 1, −1) put min(−α − 1, 1 + q + α) 6 q̄/2 Thm 6.5

−1 half-cash-secured put min(q, 1) 6 min(q̄, 1) Thm 6.6
(−1, 0) cash-secured put = covered call min(α + 1, −α) 6 1/2 Thm 6.5
0 half-covered call min(p, 1) 6 min( p̄, 1) Thm 6.6
(0, p̄) call min(α, p − α) 6 p̄/2 Thm 6.2

In other words, a high p̄ = sup{p : E exp((p + 1)X) < ∞} indicates an X distribution with thin right-hand
tail. Similarly, a high q̄ indicates an X distribution with thin left-hand tail. To control sampling error, the
guideline is to price the call if the right-hand tail is thinner in the sense that p̄ > q̄, but otherwise price the
corresponding put. However, if both tails are sufficiently thick, then instead price either a covered call or
one of the “hybrid” payoffs induced by α = −1 or 0.
Now consider truncation error. To develop intuition, we treat only the trivial case where X has zero
variance; it should be understood that the resulting rule of thumb will lose accuracy for X with high variance.
Writing F0 = S0 /B0 for the T -forward price, we have

| f (u − (α + 1)i)| |B0 Eei(u−(α+1)i)X | B0 F0α+1 S0 F0α

|ĉα,G1 (u)| 6 = = = 2 .
u2 u2 u2 u

So Theorem 6.1 gives the truncation error bound

S0 (F0 /K)α
, (7.2)
πN∆

which is increasing in α if K < F0 , but decreasing in α if K > F0 . The rule of thumb, therefore, is that
to minimize truncation error at low strikes, price the put; at high strikes, price the call. Specifically, in this
zero-variance case, the rule for a given strike is to price whichever contract (put or call) is out-of-the-money-
forward at that strike.
By combining these sampling error and truncation error heuristics, we will generate overall recommen-
dations.

25
7.4 Recommendations

In this subsection, assume that we are to price options on equities, whose (risk-adjusted) return distributions
typically exhibit significant negative skew, consistent with p̄ > q̄.
For the common task of pricing a set of contracts with strikes nearly at-the-money-forward, the sampling-
error heuristics (Table 1) tend to outweigh the truncation-error heuristics (7.2), which are in this case rela-
tively insensitive to α. Therefore we recommend that the default procedure be to take α > 0 and price the
call.
The primary exception to this rule occurs at strikes away from the money, where truncation error tends
to become more α-sensitive. At large strikes, the effect of (7.2) still favors the α > 0 strategy of pricing
the call. However, for small strikes, it favors the opposite strategy; indeed for strikes sufficiently small, this
effect can swamp the sampling error effect, resulting in the opposite recommendation: take α < −1 and
price the put.
A second exception arises when the X distribution’s tails are thick, in the sense that p̄ < 2 (implying
that stock prices have infinite third moment). In this case the default procedure should be to price the half-
covered call, by taking α = 0, which outperforms α = −1 on sampling error bounds. However, if the X
distribution’s tails are very thick, in the sense that p̄ < 1/2, then the default procedure should be to price
price the covered call, by taking α = −1/2, which uniquely in (−1, 0) attains the sampling bound decay
rate of 1/2.
A third exception could arise when one wishes to avoid optimizing α, possibly because of the com-
puting, programming, or mathematical effort involved (where “mathematical” effort refers to the analytic
determination of p̄ and q̄, given an unfamiliar characteristic function). Suppose one needs only a simple
choice for α, that still guarantees the error bounds will go to zero for large N. Since AX contains the interval
[0, 1], the three choices α ∈ {0, −1/2, −1} are all acceptable.
Some caveats apply to the “simple” choices α ∈ {0, −1/2, −1}. If one declines to optimize over α > 0
or α < −1, then one relinquishes the possibility of obtaining better error bounds – possibly much better,
especially for N small and p̄ large. Moreover, we put “simple” in quotation marks, because such α still
require choices for ∆ and N; and making those choices in an efficient way still requires optimization in
some sense. Note that these caveats apply also to the traditional approach of computing a difference of two
integrals, which has furthermore the error-management disadvantages of Section 7.1.

26
7.5 Numerical Examples and Discussion

For numerical examples we take the Variance Gamma model in Table 2 and the Heston model in Table 3.
The parameters come from empirical studies of S&P 500 futures options: VG parameters from Madan-Carr-
Chang (1998) which uses data from 1992–1994, and Heston parameters from Bates (2000) which uses data
from 1988–1993.
For each model we generate two sub-tables: one for options at T = 1 month, and one for options at
T = 4 months to expiry. Each one shows calls with strikes ranging from 80 to 120, on an underlying with
value 100.
For each model and expiry, we choose the number of quadrature points N large enough to guarantee
error smaller than one penny (0.01) at all listed strikes. For each strike and each of the five α regimes, we
choose α and ∆ to minimize our error bounds. The tables report the a priori error bounds and the realized
errors.

Remark 7.1. These tables demonstrate that in certain examples with plausible parameters, we can guarantee
accuracy of within one penny (which is 0.0001 times the underlying S0 = 100), by sampling at a number of
points N not in the thousands, but instead well under one hundred, and indeed in some cases under ten.

Remark 7.2. For each contract, our recommended quadrature parameters delivered a realized accuracy of
within one-tenth of a penny (which is 0.00001 times the underlying).

Remark 7.3. The effect of increasing the time to expiry T depends on the model. To maintain the same
accuracy in the VG model the computational burden decreased (from N = 32 to N = 8), whereas in the
Heston model the burden increased (from N = 8 to N = 16).
Consider the following two pieces of intuition. One effect of increasing T is that return densities become
smoother, which thins the tails of the characteristic function, hence decreasing truncation error; this effect
is more significant in VG than Heston, because the former characteristic function decays polynomially, but
the latter decays exponentially. Another effect of increasing T , however, is that return densities have fatter
tails, which tends to make the characteristic function less smooth, hence increasing sampling error; this
effect is more significant in Heston than VG, because the measure of tail-thinness relevant to our bounds
is the number of finite moments. Under VG, the returns follow a Lévy process so the number of moments
is invariant to time horizon, unlike Heston, where volatility is persistent, and hence works to decrease the
number of moments and increase sampling error as T increases. To see this numerically in Tables 2 and 3,
refer to AX , which shows the number of moments to be T -dependent under Heston, but not under VG.

Remark 7.4. The optimal choice of α regime in our examples agrees with the rules of thumb proposed

27
in Sections 7.3 and 7.4. Near-the-money and out-of-the-money, the best error bounds arise from choosing
α > 0 and pricing the call. However, for strikes sufficiently deep-in-the-money, the best error bounds arise
from choosing α < −1 and pricing the out-of-the-money put. The other choices α ∈ [−1, 0] underperformed,
as we would anticipate, given the sufficiently thin tails of the return distributions; the thickest (Heston at 4
months) had p̄ = 24.32 and q̄ = 9.97, well above 2.

Remark 7.5. The numerics reflect a general viewpoint of this paper, which holds that the freedom to choose
integration path, via the α parameter, plays an essential role in the accuracy, efficiency, and robustness of
the transform approach.

28
 Table 2: Realized errors and a priori bounds in five α regimes, under VG dynamics.

2(a). T = 1 month. N = 32 points.

Strike: 80 90 100 110 120

Call price: 20.0057 10.0877 1.2678 0.0138 0.0004

α Error Bound Error Bound Error Bound Error Bound Error Bound

< −1 0.0000 0.0006 -0.0001 0.0032 -0.0061 0.0128 -0.0008 0.0370 -0.0115 0.0829
= −1 0.1045 0.5271 0.0880 0.5771 0.4670 0.6258 -0.0109 0.6731 0.2862 0.7193
∈ (−1, 0) 0.5403 1.8857 0.0419 2.0033 1.4950 2.1127 0.1028 2.2151 0.2985 2.3115
=0 0.0872 0.5987 0.1248 0.6157 0.4632 0.6312 -0.0271 0.6453 0.2713 0.6584
>0 -0.0147 0.1056 -0.0019 0.0342 0.0005 0.0058 -0.0000 0.0006 -0.0000 0.0001

Optimal α -14.98 -13.63 21.25 25.56 28.56

2(b). T = 4 months. N = 8 points.

Strike: 80 90 100 110 120

Call price: 20.0565 10.4903 2.8992 0.2310 0.0129

α Error Bound Error Bound Error Bound Error Bound Error Bound

< −1 -0.0000 0.0013 0.0004 0.0057 -0.0001 0.0191 -0.0138 0.0505 -0.0338 0.1109
= −1 4.3683 7.4092 4.3888 7.7501 5.0984 8.0651 4.8697 8.3584 4.9340 8.6331
∈ (−1, 0) 21.7868 36.2630 23.1801 38.5159 25.3678 40.6206 25.9590 42.5988 26.4961 44.4674
=0 4.3509 7.1282 4.3498 7.6780 5.1691 8.2023 5.1130 8.7045 5.1763 9.1873
>0 -0.0144 0.0923 0.0017 0.0259 -0.0001 0.0055 -0.0003 0.0009 -0.0000 0.0001

Optimal α -12.85 -11.73 17.55 20.66 23.31

Underlying: S0 = 100.

Variance Gamma model parameters: σ = 0.1213, ν = 0.1686, θ = −0.1436.

The interval of permissible α + 1 values is AX = (−20.26, 39.78) at both time horizons.

CPU time for the N-point quadrature evaluations totaled 0.25 sec. for all of 2(a), and 0.08 sec. for all of 2(b).

29
 Table 3: Realized errors and a priori bounds in five α regimes, under Heston dynamics.

3(a). T = 1 month. N = 8 points.

Strike: 80 90 100 110 120

Call price: 20.0043 10.1213 1.8314 0.0150 0.0001

α Error Bound Error Bound Error Bound Error Bound Error Bound

< −1 -0.0000 0.0003 -0.0005 0.0034 -0.0012 0.0225 -0.0066 0.0903 -0.0365 0.2529
= −1 10.1124 15.2531 10.1303 15.7479 10.8452 16.2042 10.5570 16.6283 11.1423 17.0249
∈ (−1, 0) 39.6391 66.4649 41.6391 70.5012 44.6160 74.3163 46.0391 77.9426 48.3202 81.4054
=0 9.0364 13.8209 9.6049 15.0579 10.8497 16.2574 11.1304 17.4243 12.1809 18.5626
>0 -0.0547 0.2081 -0.0021 0.0383 0.0001 0.0031 -0.0000 0.0001 -0.0000 0.0000

Optimal α -22.20 -19.34 33.12 44.03 52.80

3(b). T = 4 months. N = 16 points.

Strike: 80 90 100 110 120

Call price: 20.3808 11.2277 3.7412 0.5343 0.0770

α Error Bound Error Bound Error Bound Error Bound Error Bound

< −1 -0.0004 0.0078 -0.0009 0.0157 -0.0057 0.0284 -0.0186 0.0473 -0.0170 0.0735
= −1 0.4526 0.9876 0.2692 1.0609 0.6069 1.1306 0.4485 1.1972 0.3552 1.2611
∈ (−1, 0) 2.4046 5.6995 2.5307 6.0546 3.5693 6.3855 3.3662 6.6960 2.7545 6.9886
=0 0.5299 1.0613 0.2947 1.1131 0.5986 1.1610 0.4855 1.2058 0.3176 1.2479
>0 -0.0017 0.0107 -0.0010 0.0040 -0.0001 0.0015 -0.0001 0.0005 -0.0000 0.0002

Optimal α -6.11 9.84 10.96 11.98 12.91

Underlying: S0 = 100.

Heston model parameters: κ = 1.49, θ = 0.0671, σ = 0.742, ρ = −0.571. State variable: V0 = 0.0262.

Interval of permissible α + 1 values is AX = (−38.41, 89.59) at T = 1, and AX = (−9.97, 25.32) at T = 4 months.

CPU time for the N-point quadrature evaluations totaled 0.16 sec. for all of 3(a), and 0.29 sec. for all of 3(b).

30
A Appendix: Examples of Known Characteristic Functions

A.1 Variance Gamma

Reference: Madan-Carr-Chang (1998). The VG model has parameters σ , θ , ν.

The log price X = log ST has discounted characteristic function

exp(−rT + iζ [log S0 + µT ])
f (ζ ) = ,
(1 − iνθ ζ + νσ 2 ζ 2 /2)T /ν
where µ := r + (1/ν) log(1 − θ ν − σ 2 ν/2); its domain is the strip ΛX induced by AX = (a− , a+ ), where
r
θ 2 θ2
a± = − 2 ± + . (A.1)
σ νσ 2 σ 4
We have the following bound on the large-u decay of f . For u > 0,

| f (u + wi)| 6 φ (w)u−2T /ν , (A.2)

where
φ (w) := exp(−rT − w[log S0 + µT ])(νσ 2 /2)−T /ν . (A.3)

Hence the VG model’s discounted characteristic function satisfies the power decay condition of Theorem
6.1. So in, for example, the case G = G1 , one can take Φα,G1 (u) = φ (−(α + 1)) and 1 + γ = 2 + 2T /ν.

A.2 Square-Root Stochastic Volatility

Reference: Heston (1993). The model has parameters κ, θ , σ , ρ, and a state variable V0 .
The log price X = log ST has discounted characteristic function

f (ζ ) = exp[−rT + iζ (log S0 + rT ) +C(ζ ) + D(ζ )V0 ],

where

1 − gedT
  
κθ
C(ζ ) := (κ − ρσ ζ i + d)T − 2 log
σ2 1−g
dT
 
κ − ρσ ζ i + d 1 − e
D(ζ ) :=
σ2 1 − gedT
κ − ρσ ζ i + d
g := g(ζ ) :=
κ − ρσ ζ i − d
q
d := d(ζ ) := (ρσ ζ i − κ)2 + σ 2 (ζ i + ζ 2 ).

The square root and the complex logarithm are multi-valued functions. For the square root here, either of
the two values may be chosen, because f is even in d. For the logarithm, however, choosing the wrong value

31
can lead to wildly incorrect answers. To define f (wi) for real w, the correct choice of log(z) is the principal
branch log |z| + arg(z), where −π < arg(z) < π. However, as pointed out by Schöbel and Zhu (1999), to
define f (ζ ) for general ζ , the correct choice of log is not necessarily the principal branch. Instead, the value
of log when ζ = u + wi is determined by the analyticity of f , which implies that log must vary continuously
as ζ varies from 0 + wi to u + wi.
This issue presents a challenge to the traditional approach of taking the Fourier integrals in Heston
and simply passing the integrands into a numerical integration routine from a standard software library.
Enforcing the required continuity of the log is tricky if the integration routine samples the integrand at an
unpredictable sequence of points. On the other hand, for a method, such as ours, that samples the integrand
at an increasing sequence of points with spacing ∆, enforcing continuity typically does not present any
difficulty.
The domain of f is the strip ΛX induced by AX = (a− , a+ ), where a− < 0 and a+ > 1 solve

g(−ia) exp(d(−ia)T ) = 1.

Specifically, if we assume κ − ρσ > 0, then a− is the largest (closest to 0) solution in (−∞, y− ), and a+ is
the smallest solution in (y+ , ∞), where
p
σ − 2κρ ± σ 2 − 4κρσ + 4κ 2
y± := .
2σ (1 − ρ 2 )
For ζ = u + wi we bound the large-u decay of f , as follows. Define

HR1 (u) := u2 σ 2 (1 − ρ 2 )

HR2 (w) := w2 σ 2 (1 − ρ 2 ) − w(2κρσ − σ 2 ) − κ 2

HR (u, w) := Re(d 2 ) = HR1 (u) − HR2 (w)

HI (u, w) := Im(d 2 ) = σ u(2wσ (1 − ρ 2 ) + σ − 2κρ)

p
h(u, w) := HR (u, w),

and define
√
∗ κ |σ − 2κρ| + κ 2 /(σ u2 + w2 )
g (u, w) := √ + p
σ u2 + w2 h(u, w) + (u2 − w2 )σ 2 (1 − ρ 2 )

g(u, w) := (1 − g∗ (u, w)) (1 + g∗ (u, w))



  
1 1
J(u, w) := 1 + 1+ .
g(u) g(u) exp(T h(u, w)) − 1

32
Let u0 > |w| satisfy 1 > g∗ (u0 , w) and h(u0 , w) > (1/T ) max(log(1/g(u0 , w)), 1) and HR1 (u0 ) > |HR2 (w)|.
Then for all u > u0 , we have
 
p V0 + κθ T
| f (u + wi)| 6 φ (u, w) exp − 1 − ρ2 u ,
σ

where (suppressing the arguments (u, w) for convenience) we let

 
2κθ /σ 2 V0 + κθ T p
φ (u, w) := J exp − rT − (log S0 + rT )w + (κ + ρσ w + max(0, HR2 ))
σ2
 
V0 J  p p
× exp 2 κ + |ρσ u| max(1, HR /HR1 ) + |ρσ w| + HR + |HI | .
σ exp(T h)

Hence the square-root stochastic volatility model’s discounted characteristic function satisfies the expo-
nential decay condition of Theorem 6.1. So in, for example, the case G = G1 , one can take Φα,G1 (u) =
φ (u, −(α + 1))/u2 .

B Appendix: Sampling Error Bounds for Payoffs G2 and G4

Sections 6.2-6.4 gave sampling error bounds for G = G1 and G = G3 , which are the payoff classes of
greatest practical interest. The other cases G = G2 and G = G4 can be treated by similar techniques, albeit
with messier results. Specifically, the G2 /G4 version of Theorem 6.2 is as follows.

Theorem B.1. Assume that b1 ∈ AX and αb0 + b1 ∈ AX with α > 0.

In case G = G2 we have

e−2πα/∆ (−i f 0 (0) − k) 4πe−2πα/∆ f (0)


|CG2 − Σ∞
α,G2 | 6 inf +
p>α: p∈AX 1 − e−4πα/∆ ∆(1 − e−4πα/∆ )2
q>0:−q∈AX

e−2π(q+α)/∆ f (iq) e2π(α−p)/∆ f (−ip)


+ + .
qe1−qk (1 − e−4π(q+α)/∆ ) pe pk+1 (1 − e4π(α−p)/∆ )

In case G = G4 , assume also that ĉ(u) = O(u−1−γ ) as u → ∞, where γ > 0. Then

χb2 · ∇ f (−ib1 )
 −2πα/∆
∞ e
|CG4 − Σα,G4 | 6 max inf
χ=±1 (p0 ,p2 ,q0 ,q2 ) 1 + e−2πα/∆
f (−i(−q0 b0 − q2 χb2 + b1 )) + e−2π(α+q0 )/∆ f (−i(−q0 b0 + q2 χb2 + b1 ))
+
q2 e−q0 k+1 (e2π(α+q0 )/∆ − e−2π(α+q0 )/∆ )
f (−i(p0 b0 + p2 χb2 + b1 )) + e2π(α−p0 )/∆ f (−i(p0 b0 − p2 χb2 + b1 ))

+ ,
p2 e p0 k+1 (e−2π(α−p0 )/∆ − e2π(α−p0 )/∆ )

where the inf is over all positive p0 , p2 , q0 , q2 such that p0 > α and p0 b0 ± p2 b2 + b1 ∈ AX and −q0 b0 ∓
q2 b2 + b1 ∈ AX .

33
Proof. In the proof of Theorem 6.2, replace equations (6.4) and (6.5) with
∞ 
Ee−qX e2π j(α−p)/∆ Ee pX
  
|CG2 (k) − Σ (k)| 6 B0 ∑ e
∞ −2π jα/∆
EX − k + 2π( j + 1)/∆ + 1−q(k−2π j/∆) +
j=1 qe pekp+1
j odd

and
∞ j
Ee(−q0 b0 −χ(−1) q2 b2 +b1 )·X

|CG4 (k) − Σ (k)| 6 B0 max ∑ e
∞ −2π jα/∆ j
E((χ(−1) b2 · X)e ) + b1 ·X
χ=±1 j=1 q2 e−q0 k+1 e2π j(α+q0 )/∆
j
e2π j(α−p0 )/∆ Ee(p0 b0 +χ(−1) p2 b2 +b1 )·X

+ .
p2 e p0 k+1

The rest of the proof holds.

C Appendix: Application of DFT/FFT

We treat here two issues: a recipe for DFT evaluation of the quadrature scheme, and modifications to the
Section 7.2 optimization problem so that the DFT output has the desired contract spacing.
First, our Davies-style discretization samples ĉ at the midpoints (n + 1/2)∆ of intervals of length ∆,
whereas Carr-Madan’s sampling scheme applies Simpson’s rule to the endpoints of those intervals. For
completeness we describe how to adapt their formulas to midpoint sampling.
Define the discrete Fourier transform (DFT) of an N-vector x to be the vector X where
N
Xm = ∑ e−i(2π/N)(n−1)(m−1) xn , m = 1, . . . , N.
n=1

Other definitions exist; this is the one in Carr-Madan, and in a number of standard software packages,
including Matlab.
Using a spacing of λ = 2π/(N∆) between consecutive triggers, we want to use DFT to compute prices
ΣN (km ) at triggers
km := k1 + λ (m − 1), m = 1, . . . , N

for arbitrary k1 . Typically one would choose k1 such that the interval [k1 , k1 + (N − 1)λ ] contains all of the
contracts to be priced. By (6.1),
∆
 N 
Σ (km ) = αk Re ∑ ĉ((n − 1/2)∆)e
N −i(n−1/2)(k1 +λ (m−1))∆
πe m n=1
∆ N
 
= αk Re e
πe m
−i(m−1)λ ∆/2
∑ ĉ((n − 1/2)∆)e −i(n−1)(m−1)λ ∆ −i(n−1/2)k1 ∆
e
n=1
∆ N
 
= αk Re e−iπ(m−1)/N ∑ e−i(2π/N)(n−1)(m−1) ĉ((n − 1/2)∆)e−ik1 (n−1/2)∆ ,
πe m n=1

34
where the sum is computable as the m-th component of the DFT of the vector whose n-th component is
ĉ((n − 1/2)∆) exp(−ik1 (n − 1/2)∆).
The second issue is the reciprocity relation λ ∆ = 2π/N. For N fixed, a decrease in Fourier-domain grid
spacing ∆ would cause the contract spacing λ to increase.
If one wishes to impose an upper bound λ̄ on spacing between contracts, then the minimum in (7.1)
should be taken over ∆ > 2π/(N λ̄ ) instead of ∆ > 0. Moreover, in certain instances it is desirable to
constrain ∆ to be an integer times 2π/(N λ̄ ), because this forces λ to divide λ̄ , so that a set of contracts with
trigger spacings of λ̄ can be priced in a single DFT, without interpolation.
We deferred this material to an Appendix to emphasize that the analysis in the body of this paper does
not make any assumption on how the sum in (6.1) is computed (aside from absence of roundoff error). One
can use the DFT; or its efficient implementation the fast Fourier transform (FFT); or, perhaps even more
efficiently (if few enough strikes need to be simultaneously priced), simple direct summation.

References

Bakshi, G., C. Cao, and Z. Chen (1997). Empirical performance of alternative option pricing models.
Journal of Finance 52, 2003–2049.

Bakshi, G. and D. Madan (2000). Spanning and derivative-security valuation. Journal of Financial Eco-
nomics 55, 205–238.

Barndorff-Nielsen, O. (1998). Processes of normal inverse Gaussian type. Finance and Stochastics 2,
41–68.

Barndorff-Nielsen, O., E. Nicolato, and N. Shephard (2002). Some recent developments in stochastic
volatility modelling. Quantitative Finance 2, 11–23.

Bates, D. (1996). Jumps and stochastic volatility: Exchange rate processes implicit in deutsche mark
options. Review of Financial Studies 9(1), 69–107.

Bates, D. (2000). Post-’87 crash fears in the S&P 500 futures option market. Journal of Econometrics 94,
181–238.

Boyarchenko, S. I. and S. Z. Levendorskiǐ (2002). Non-Gaussian Merton-Black-Scholes Theory. World

Scientific.

Broadie, M., J. Cvitanic, and H. M. Soner (1998). Optimal replication of contingent claims under portfo-
lio constraints. Review of Financial Studies 11(1), 59–79.

35
Carr, P., H. Geman, D. Madan, and M. Yor (2002, April). The fine structure of asset returns: An empirical
investigation. Journal of Business 75(2), 305–332.

Carr, P., H. Geman, D. Madan, and M. Yor (2003). Stochastic volatility for Lévy processes. Mathematical
Finance 13(3), 345–382.

Carr, P. and D. Madan (1999). Option valuation using the fast Fourier transform. Journal of Computa-
tional Finance 3, 463–520.

Carr, P. and M. Schröder (2003). Bessel processes, the integral of geometric Brownian motion, and Asian
options. Theory of Probability and Its Applications. Forthcoming.

Carr, P. and L. Wu (2002). Time-changed Lévy processes and option pricing. Journal of Financial Eco-
nomics. Forthcoming.

Carr, P. and L. Wu (2003). The finite moment log stable process and option pricing. Journal of Fi-
nance 58(2), 753–777.

Chacko, G. and S. Das (2002). Pricing interest rate derivatives: A general approach. Review of Financial
Studies 15(1), 195–241.

Champeney, D. (1987). A Handbook of Fourier Theorems. Cambridge University Press.

Davies, R. (1973). Numerical inversion of a characteristic function. Biometrika 60(2), 415–417.

Delbaen, F. and W. Schachermayer (1994). A general version of the fundamental theorem of asset pricing.
Mathematische Annalen 2(4), 61–73.

Duffie, D., J. Pan, and K. Singleton (2000). Transform analysis and option pricing for affine jump-
diffusions. Econometrica 68(6), 1343–1376.

Eberlein, E. and K. Prause (2002). The generalized hyperbolic model: Financial derivatives and risk
measures. In H. Geman, D. Madan, S. Pliska, and T. Vorst (Eds.), Mathematical Finance – Bachelier
Congress 2000, pp. 245–267. Springer.

El Karoui, N., H. Geman, and J.-C. Rochet (1995). Changes of numeraire, changes of probability mea-
sure, and option pricing. Journal of Applied Probability 32, 443–458.

Fu, M., D. Madan, and T. Wang (1999). Pricing continuous Asian options: A comparison of Monte Carlo
and Laplace transform inversion methods. Journal of Computational Finance 2(2), 49–74.

Geman, H. and M. Yor (1996). Pricing and hedging double-barrier options: A probabilistic approach.
Mathematical Finance 6(4), 365–378.

36
Gil-Pelaez, J. (1951). Note on the inversion theorem. Biometrika 38(3/4), 481–482.

Harrison, J. M. and D. Kreps (1979). Martingales and arbitrage in multiperiod securities markets. Journal
of Economic Theory 20, 381–408.

Heston, S. L. (1993). A closed-form solution for options with stochastic volatility and applications to
bond and currency options. Review of Financial Studies 6(2), 327–343.

Lee, R. (2003). The moment formula for implied volatility at extreme strikes. Mathematical Finance.
Forthcoming.

Lewis, A. (2000). Option Valuation under Stochastic Volatility. Finance Press.

Lewis, A. (2001). A simple option formula for general jump-diffusion and other exponential Lévy pro-
cesses. Envision Financial Systems and [Link].

Madan, D., P. Carr, and E. Chang (1998). The variance gamma process and option pricing. European
Finance Review 2, 79–105.

Pan, J. (2002). The jump-risk premia implicit in options: Evidence from an integrated time-series study.
Journal of Financial Economics 63, 3–50.

Pelsser, A. (2000). Pricing double barrier options using Laplace transforms. Finance and Stochastics 4,
95–104.

Schöbel, R. and J. Zhu (1999). Stochastic volatility with an Ornstein-Uhlenbeck process: An extension.
European Finance Review 3, 23–46.

Zemanian, A. (1966). The distributional Laplace and Mellin transformations. SIAM Journal on Applied
Mathematics 14(1), 41–59.

Zhu, J. (2000). Modular Pricing of Options. Springer.

37