# Pricing Extendible Options Using the Fast Fourier Transform

Authors: Siti Nur Iqmal Ibrahim; John G. O'Hara; Nick Constantinou
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper applies the fast Fourier transform (FFT) approach, within the Black-Scholes framework, to the valuation of options whose time to maturity can be extended to a future date (extendible options). We determine the valuation of the extendible options as sums of expectations of indicator functions, leading to a semianalytic expression for the value of the options over a range of strikes. Compared to Monte Carlo simulation, numerical examples demonstrate that the FFT is both computationally more efficient and higher in accuracy.


---

## Body

## 1. Introduction

An extendible option allows its time to expiration to be extended to a future date for a fee with a different strike price at the extended expiration time. Longstaff [1] discusses a number of applications of extendible options, such as extendible bonds and warrants, and applications in real estate, demonstrating that extendible options are important in financial markets. For example, Dias and Rocha [2] price petroleum concessions using a mean-reverting framework including jumps for an extendible option. Extendible options can be classified as a holder- or writer-extendible call option or a holder- or writer-extendible put option. Options that are extended by their writer have no fee and they are more straightforward to price. To the best of our knowledge, Longstaff [1] was the first to study extendible options. Chung and Johnson [3] extend the work of Longstaff [1] to a general case where the holder or the writer can extend the option more than once and provide a general closed-form pricing solution for n-time extendible options.In order to price options, many numerical techniques are available, such as Monte Carlo simulation and the Fourier transform method. The fast Fourier transform (FFT) was first introduced by Cooley and Tukey [4] to obtain the discrete Fourier transform (DFT) and its inverse. Implementing the DFT involves N points where, for added efficiency, N is a power of 2 and it is essential to have an extensively large N data set, which reduces the number of computations from O
(
N
2
) to O
(
N
ln
⁡
⁡
N
). Monte Carlo simulation often requires substantially more time to compute option prices, while the FFT, which was first implemented to option pricing by Carr and Madan [5], is often more straightforward and effective. This has been demonstrated by a number of authors, such as Pillay and O’Hara [6] and Ibrahim et al. [7]. Hence, the FFT has gained popularity as an effective option pricing methodology (see [8, 9]), and a significant amount of research has extended this method to price various types of options, for instance Asian options [10], spread options [11], and power options [12]. The FFT technique is widely applicable and has been shown to perform well in terms of its speed and accuracy in some other affine models, such as stochastic volatility, stochastic interest rates, and jumps (see [13–15]).In this paper, we review the pricing framework for extendible call options under the Black and Scholes [16] model first developed by Longstaff [1]. By following the work of Carr and Madan [5], we apply the FFT method to numerically compute the prices of extendible call options. For an effective implementation of the FFT, we provide an analytical representation of the characteristic functions for the extendible options. The usefulness of the FFT is that this technique gives a range of option prices for a collection of strike prices for a given maturity. In addition, we estimate the price of extendible options via Monte Carlo simulation for comparison.The paper is organized as follows. Section2 reviews the closed-form formulas for the holder-extendible call options under the Black-Scholes model as provided in Longstaff [1]. Section 2.1 develops a pricing procedure for holder-extendible call options using the FFT approach where the characteristic functions are given, while Section 2.2 applies the Monte Carlo simulation using an Euler scheme and a Milstein scheme. Section 3 presents the numerical results, and Section 4 concludes the paper.


## 2. Extendible Option

In this section, we derive a closed-form solution for holder-extendible call options as given by Longstaff [1], which is the focus of this paper. An analogous approach can be used to value extendible put options.Let(
Ω
,
F
,
Q
) be a probability space on which a Brownian motion, W
t
,
0
≤
t
≤
T, is defined, where F
t, 0
≤
t
≤
T, is the natural filtration generated by the Brownian motion. Suppose that Q is a risk-neutral measure under which the asset price process S
t
,
0
≤
t
≤
T, is governed by the following dynamics:
(1)
d
S
t
=
r
S
t
d
t
+
σ
S
t
d
W
t
.
Within the Black-Scholes framework, the volatility σ and the risk-free interest rate r are assumed to be constant. Using Itô’s lemma, the process followed by the logarithmic asset price x
t
:
=
ln
⁡
S
t is given by
(2)
d
x
t
=
(
r
-
1
2
σ
2
)
d
t
+
σ
d
W
t
.
This demonstrates that x
t follows a generalized Wiener process with a constant drift (
r
-
(
1
/
2
)
σ
2
) and a constant volatility σ. Therefore, the change in x between time t and T is normally distributed with a mean of (
r
-
(
1
/
2
)
σ
2
)
(
T
-
t
) and a variance of σ
2
(
T
-
t
) or
(3)
ln
⁡
S
T
~
N
[
ln
⁡
S
t
+
(
r
-
1
2
σ
2
)
(
T
-
t
)
,
σ
2
(
T
-
t
)
]
.
Therefore, S
T has a log-normal distribution. The valuation formula for a vanilla call option with a given underlying asset price S
t, strike price K, expiration time T, risk-free interest rate r, and volatility σ is given as follows:
(4)
C
(
S
t
,
K
,
T
-
t
)
=
S
t
N
(
d
1
)
-
K
e
-
r
(
T
-
t
)
N
(
d
2
)
,
where
(5)
d
1
=
ln
⁡
(
S
t
/
K
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
-
t
)
σ
T
-
t
,
d
2
=
d
1
-
σ
T
-
t
.LetEC be the value of an extendible call option with time to expiration of T
1. At the time to expiration T
1, the holder of the extendible call can(1)
let the call expire worthless ifS
T
1
<
L, or
(2)
exercise the call and getS
T
1
-
K
1 if S
T
1
>
H, or
(3)
make a payment of an additional premiumA to extend the call to T
2 with a new strike of K
2 if L
≤
S
T
1
≤
H,where S
T
1 is the underlying asset price and strike price at time T
1
,
K
1 is the strike price at time T
1, and Longstaff [1] refers to L and H as critical values, where L
<
H.If at expiration timeT
1 the option is worth more than the extendible value with a new strike price of K
2 for a fee of A for extending the expiration time T
1 to T
2, then it is best to exercise; that is, S
T
1
-
K
1
≥
C
(
S
T
1
,
K
2
,
T
2
-
T
1
)
-
A. Otherwise, it is best to extend the expiration time of the option to T
2 and exercise when it is worth more than zero; that is, C
(
S
T
1
,
K
2
,
T
2
-
T
1
)
-
A
>
0. Moreover, the holder of the option should be impartial between extending and not exercising at value L and impartial between exercising and extending at value H. Therefore, the critical values L and H are unique solutions of H
-
K
1
=
C
(
H
,
K
2
,
T
2
-
T
1
)
-
A and C
(
L
,
K
2
,
T
2
-
T
1
)
-
A
=
0. See Longstaff [1] and Gukhal [17] for an analysis of the conditions.The price of the extendible call option is the discounted expectation of the payoff at timeT
1 and is given by the following:
(6)
EC
(
S
t
,
K
1
,
T
1
,
K
2
,
T
2
,
A
)
=
e
-
r
(
T
1
-
t
)
a
a
a
a
×
E
Q
{
max
⁡
⁡
[
S
T
1
-
K
1
,
C
(
S
T
1
,
K
2
,
T
2
-
T
1
)
-
A
,
0
]
}
(7)
=
e
-
r
(
T
1
-
t
)
E
Q
[
(
S
T
1
-
K
1
)
1
{
S
T
1
>
H
}
∣
F
t
]
+
e
-
r
(
T
1
-
t
)
E
Q
[
1
{
L
≤
S
T
1
≤
H
}
∣
F
t
(
C
(
S
T
1
,
K
2
,
T
2
-
T
1
)
-
A
)
a
a
a
a
a
a
a
a
a
a
a
a
a
×
1
{
L
≤
S
T
1
≤
H
}
∣
F
t
]
.
An example illustration of its payoff is provided in Figure 1.Figure 1
An example of the payoff function for an extendible call option forK
1
=
100
,
K
2
=
105
,
r
=
0.08
,
σ
=
0.25
,
T
1
=
0.5
,
T
2
=
0.75
,
A
=
1.Formulating the value of a call option,C at time T
1 with a time to expiration extended to T
2, as the discounted conditional expected payoff,
(8)
e
-
r
(
T
2
-
T
1
)
E
Q
[
max
⁡
⁡
(
S
T
2
-
K
2
,
0
)
∣
F
t
]
,
yields
(9)
EC
(
S
t
,
K
1
,
T
1
,
K
2
,
T
2
,
A
)
=
e
-
r
(
T
1
-
t
)
E
Q
[
(
S
T
1
-
K
1
)
1
{
S
T
1
>
H
}
∣
F
t
]
+
e
-
r
(
T
2
-
t
)
E
Q
×
[
(
S
T
2
-
K
2
)
1
{
S
T
2
>
K
2
}
(
1
{
S
T
1
≥
L
}
-
1
{
S
T
1
≥
H
}
)
∣
F
t
]
-
e
-
r
(
T
1
-
t
)
E
Q
[
A
(
1
{
S
T
1
≥
L
}
-
1
{
S
T
1
≥
H
}
)
∣
F
t
]
,or(10)
EC
(
S
t
,
K
1
,
T
1
,
K
2
,
T
2
,
A
)
=
S
t
E
Q
S
[
1
{
x
T
1
>
h
}
]
-
K
1
e
-
r
(
T
1
-
t
)
E
Q
[
1
{
x
T
1
>
h
}
]
+
S
t
E
Q
S
[
1
{
x
T
1
≥
l
,
x
T
2
≥
k
2
}
]
-
K
2
e
-
r
(
T
2
-
t
)
E
Q
[
1
{
x
T
1
≥
l
,
x
T
2
≥
k
2
}
]
-
S
t
E
Q
S
[
1
{
x
T
1
≥
h
,
x
T
2
≥
k
2
}
]
+
K
2
e
-
r
(
T
2
-
t
)
E
Q
[
1
{
x
T
1
≥
h
,
x
T
2
≥
k
2
}
]
-
A
e
-
r
(
T
1
-
t
)
{
E
Q
[
1
{
x
T
1
≥
l
}
]
-
E
Q
[
1
{
x
T
1
≥
h
}
]
}
,
where Q
S is the spot price measure, x
t
=
ln
⁡
S
t
,
k
1
=
ln
⁡
K
1
,
k
2
=
ln
⁡
K
2
,
l
=
ln
⁡
L, and h
=
ln
⁡
H.Corollary 1.
The price of an extendible call option with time to expirationT
1 and strike price K
1, whose expiration time may be extended to T
2 upon a payment of an additional premium A, with a new strike price K
2 is specified as
(11)
E
C
(
S
t
,
K
1
,
T
1
,
K
2
,
T
2
,
A
)
=
[
S
t
N
(
a
1
)
-
K
1
e
-
r
(
T
1
-
t
)
N
(
a
2
)
]
+
[
S
t
N
2
(
b
1
,
c
1
,
ρ
)
-
K
2
e
-
r
(
T
2
-
t
)
N
2
(
b
2
,
c
2
,
ρ
)
]
-
[
S
t
N
2
(
a
1
,
c
1
,
ρ
)
-
K
2
e
-
r
(
T
2
-
t
)
N
2
(
a
2
,
c
2
,
ρ
)
]
-
A
e
-
r
(
T
1
-
t
)
[
N
(
b
2
)
-
N
(
a
2
)
]
,
where N
(
x
) denotes the standard univariate cumulative normal distribution function, N
2
(
x
,
y
,
ρ
) denotes the standard bivariate cumulative normal distribution function with a correlation coefficient, ρ
=
(
T
1
-
t
)
/
(
T
2
-
t
), and
(12)
a
1
=
ln
⁡
⁡
(
S
t
/
H
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
1
-
t
)
σ
T
1
-
t
,
a
2
=
a
1
-
σ
T
1
-
t
,
b
1
=
ln
⁡
⁡
(
S
t
/
L
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
1
-
t
)
σ
T
1
-
t
,
b
2
=
b
1
-
σ
T
1
-
t
,
c
1
=
ln
⁡
⁡
(
S
t
/
K
2
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
2
-
t
)
σ
T
2
-
t
,
c
2
=
c
1
-
σ
T
2
-
t
.Longstaff [1] also discusses several special cases of an extendible option, such that if L
>
0 and H
=
∞, then the extendible option reduces to a compound option on C
(
S
t
,
K
2
,
T
2
-
T
1
) with strike price A, and if L
→
K
1, the value of the extension privilege approaches zero; hence the value of the extendible options is simply C
(
S
t
,
K
1
,
T
1
) (for more discussion on the special cases of the extendible options, see Longstaff [1]).


### 2.1. Fast Fourier Transform

The FFT approach evaluates an (output) array of the following form:(13)
Y
[
l
1
,
…
,
l
n
]
∶
=
∑
j
1
=
0
N
-
1
‍
⋯
∑
j
n
=
0
N
-
1
‍
e
-
i
(
2
π
/
N
)
j
1
l
1
-
⋯
-
i
(
2
π
/
N
)
j
n
l
n
X
[
j
1
,
…
,
j
n
]
,
for any complex (input) array {
X
[
j
1
,
…
,
j
n
]
∈
C
∣
j
1
,
…
,
j
n
=
0
,
…
,
N
-
1
}. In addition, to implement the FFT method to option pricing, a closed-form representation of the characteristic function is required. For the case of extendible options, we require univariate and bivariate characteristic functions, which are provided by the following lemmas.Lemma 2.
Under the risk-neutral measureQ, the univariate characteristic function is given by
(14)
φ
x
T
1
(
u
1
)
=
exp
⁡
⁡
{
i
u
1
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
1
-
t
)
]
-
1
2
σ
2
(
T
1
-
t
)
u
1
2
}
,
φ
x
T
2
(
u
2
)
=
exp
⁡
⁡
{
i
u
2
[
x
T
1
+
(
r
-
1
2
σ
2
)
(
T
2
-
T
1
)
]
aaaaaaaaa
-
1
2
σ
2
(
T
2
-
T
1
)
u
2
2
}
.Lemma 3.
Under the risk-neutral measureQ, the bivariate characteristic function is given by
(15)
φ
x
T
1
,
x
T
2
(
u
1
,
u
2
)
=
exp
⁡
⁡
{
i
u
1
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
1
-
t
)
]
-
1
2
u
1
2
σ
2
(
T
1
-
t
)
}
×
exp
⁡
⁡
{
i
u
2
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
2
-
t
)
]
-
1
2
u
2
2
σ
2
(
T
2
-
t
)
}
.The univariate characteristic function is used in Carr and Madan [5] to price a European option using FFT within the Black-Scholes dynamics. By employing the Radon-Nikodym derivatives, the characteristic function under measure Q and the spot price measure Q
S are related as follows:
(16)
φ
x
T
1
S
(
u
1
)
=
e
-
r
(
T
1
-
t
)
-
x
t
E
Q
[
e
i
(
u
1
-
i
)
x
T
1
∣
F
t
]
,
φ
x
T
1
,
x
T
2
S
(
u
1
,
u
2
)
=
e
-
r
(
T
2
-
t
)
-
x
t
E
Q
[
e
i
u
1
x
T
1
+
i
(
u
2
-
i
)
x
T
2
∣
F
t
]
.
Hence, from (10), we consider the following expectations which can be expressed in integral form as follows:
(17)
E
Q
[
1
{
x
T
1
>
h
}
]
=
∫
h
∞
p
(
x
T
1
)
d
x
T
1
,
(18)
E
Q
[
1
{
x
T
1
≥
l
,
x
T
2
≥
k
2
}
]
=
∫
l
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
(19)
E
Q
[
1
{
x
T
1
≥
h
,
x
T
2
≥
k
2
}
]
=
∫
h
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
(20)
E
Q
[
1
{
x
T
1
≥
l
}
]
=
∫
l
∞
p
(
x
T
1
)
d
x
T
1
,
(21)
E
Q
[
1
{
x
T
1
≥
h
}
]
=
∫
h
∞
p
(
x
T
1
)
d
x
T
1
,
where p
(
·
) is the probability density of X under the risk-neutral probability Q. It is clear that integrals (17), (20), and (21) are similar, so are integrals (18) and (19). Hence we first apply the FFT to integral (18). Following Carr and Madan [5], we multiply (18) by dampening factors α
1 and α
2 to ensure square-integrability as follows:
(22)
e
α
1
l
+
α
2
k
2
∫
l
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
for
α
1
,
α
2
>
0
.
Then, by applying the Fourier transform, we have the following:
(23)
ψ
(
u
1
,
u
2
)
=
∬
-
∞
∞
e
i
u
1
l
+
i
u
2
k
2
(
e
α
1
l
+
α
2
k
2
∫
l
∞
∫
k
2
∞
p
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
∫
l
∞
∫
k
2
∞
p
∫
l
∞
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
)
.
This can be written in closed-form in terms of the characteristic function as follows:
(24)
ψ
(
u
1
,
u
2
)
=
φ
x
T
1
,
x
T
2
[
u
1
-
i
α
1
,
u
2
-
i
α
2
]
(
i
u
1
+
α
1
)
(
i
u
2
+
α
2
)
.
Using the inverse Fourier transform, we can recover integral (18):
(25)
e
-
α
1
l
-
α
2
k
2
(
2
π
)
2
∬
-
∞
∞
e
-
i
u
1
l
-
i
u
2
k
2
ψ
(
u
1
,
u
2
)
d
u
2
d
u
1
.
Then, we numerically compute this by invoking the trapezoid rule (see [11]) which approximates the integral by the following sum:
(26)
e
-
α
1
l
-
α
2
k
2
(
2
π
)
2
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
-
i
u
2
,
m
k
2
ψ
(
u
1
,
j
,
u
2
,
m
)
Δ
2
Δ
1
,
where Δ
1,2 is the distance between the points of the integration grid, u
1
,
j
=
(
j
-
(
N
/
2
)
)
Δ
1
and
u
2
,
m
=
(
m
-
(
N
/
2
)
)
Δ
2
,
for
j
,
m
=
0
,
…
,
N
-
1. In order to progress, we define a grid of size (
N
×
N
) by L
2
=
{
(
l
u
,
k
2
,
p
)
∣
0
≤
u
,
p
≤
N
-
1
} where ω
1,2
>
0 is the distance between the logarithmic critical prices and the logarithmic strikes, respectively, as follows:
(27)
l
u
=
(
u
-
N
2
)
ω
1
,
k
2
,
p
=
(
p
-
N
2
)
ω
2
.
Hence, we can now evaluate the following sum: Γ
(
l
,
k
2
)
=
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
-
i
u
2
,
m
k
2
ψ
(
u
1
,
j
,
u
2
,
m
).By choosingω
1,2
Δ
1,2
=
2
π
/
N this yields the following values of sum Γ
(
l
,
k
2
) on L
2:
(28)
Γ
(
l
u
,
k
2
,
p
)
=
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
u
-
i
u
2
,
m
k
2
,
p
ψ
(
u
1
,
j
,
u
2
,
m
)
=
(
-
1
)
u
+
p
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
(
2
π
/
N
)
j
u
-
i
(
2
π
/
N
)
m
p
a
a
a
a
a
a
a
a
a
a
a
a
a
×
[
(
-
1
)
j
+
m
ψ
(
u
1
,
j
,
u
2
,
m
)
]
.
We take the (input) array as X
[
j
,
m
]
=
(
-
1
)
j
+
m
ψ
(
u
1
,
j
,
u
2
,
m
)
,
for
j
,
m
=
0
,
…
,
N
-
1. Hence, the approximation is at (
N
×
N
) different logarithmic critical prices and strike prices such that
(29)
e
-
α
1
l
u
-
α
2
k
2
,
p
(
2
π
)
2
Γ
(
l
u
,
k
2
,
p
)
Δ
2
Δ
1
,
for
0
≤
u
,
p
≤
N
-
1
.
In an analogous way, by replacing l with h, we can approximate integral (19). The two-dimensional FFT can also be reduced to a one-dimensional FFT which involves integrals (17), (20), and (21). Therefore for brevity, we exclude the development for the one-dimensional FFT in this paper.


### 2.2. Monte Carlo Simulation

In this section, we price the extendible call option using Monte Carlo simulation [18]. Given the pricing function of an extendible call option given in (6), the vanilla call option C can be priced at time T
1 with strike price K
2 and expiration time T
2. Thus we only need to compute the asset price until time T
1. Let S
T
1
,
j be the asset price at time T
1 on the jth path, and let w
T
1
,
j be the vanilla call price with strike price K
2 and expiration time T
2 at time T
1 for the jth path; that is, C
(
S
T
1
,
j
,
K
2
,
T
2
-
T
1
). Recall the closed-form solution for a vanilla call option, so we have
(30)
C
=
S
T
1
N
(
e
1
)
-
K
2
e
-
r
(
T
2
-
T
1
)
N
(
e
2
)
,
where
(31)
e
1
=
ln
⁡
⁡
(
S
T
1
/
K
2
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
2
-
T
1
)
σ
T
2
-
T
1
,
e
2
=
e
1
-
σ
T
2
-
T
1
.
Therefore, the Monte Carlo estimator for the price of an extendible call option is given by
(32)
EC
(
t
,
x
^
T
1
)
=
e
-
r
(
T
1
-
t
)
n
∑
j
=
1
n
‍
max
⁡
⁡
[
e
x
^
T
1
j
-
K
1
,
w
T
1
,
j
-
A
,
0
]
,
where n is the number of simulations. Suppose we have x
t
=
ln
⁡
S
t with the following process:
(33)
x
t
=
(
r
-
1
2
σ
2
)
d
t
+
σ
d
W
t
.
We discretize the asset path using the Euler scheme as follows:
(34)
x
j
+
1
=
x
j
+
(
r
-
1
2
σ
2
)
Δ
t
+
σ
Δ
W
j
,
t
j
=
j
Δ
t
,
Δ
W
j
=
W
t
j
+
1
-
W
t
j
=
Z
Δ
t
with
Z
~
N
(
0,1
)
.
We also apply the Milstein scheme [19] that is strongly convergent with order one such that
(35)
x
j
+
1
=
x
j
+
(
r
-
1
2
σ
2
)
Δ
t
+
σ
Δ
W
j
+
1
2
σ
2
(
(
Δ
W
)
2
-
Δ
t
)
,
t
j
=
j
Δ
t
,
Δ
W
j
=
W
t
j
+
1
-
W
t
j
=
Z
Δ
t
with
Z
~
N
(
0,1
)
.
We let [
t
=
t
0
<
t
1
<
⋯
<
t
m
=
T
] be a partition of the time interval [
t
,
T
] of m equal segments with length Δ
t
j, where Δ
t
j
=
t
j
+
1
-
t
j and t
j
=
j
T
/
m, for each j
=
1
,
…
,
M. We draw a random sample from a normal distribution with mean 0 and volatility 1 for each j
=
0,1
,
…
,
m to simulate (34) and (35), thus generating a sample path for x
T by simulating x
j for j
=
1
to
j
=
m. This step is repeated to generate many paths. We then estimate the price of the extendible call option with Monte Carlo simulation.


## 2.1. Fast Fourier Transform

The FFT approach evaluates an (output) array of the following form:(13)
Y
[
l
1
,
…
,
l
n
]
∶
=
∑
j
1
=
0
N
-
1
‍
⋯
∑
j
n
=
0
N
-
1
‍
e
-
i
(
2
π
/
N
)
j
1
l
1
-
⋯
-
i
(
2
π
/
N
)
j
n
l
n
X
[
j
1
,
…
,
j
n
]
,
for any complex (input) array {
X
[
j
1
,
…
,
j
n
]
∈
C
∣
j
1
,
…
,
j
n
=
0
,
…
,
N
-
1
}. In addition, to implement the FFT method to option pricing, a closed-form representation of the characteristic function is required. For the case of extendible options, we require univariate and bivariate characteristic functions, which are provided by the following lemmas.Lemma 2.
Under the risk-neutral measureQ, the univariate characteristic function is given by
(14)
φ
x
T
1
(
u
1
)
=
exp
⁡
⁡
{
i
u
1
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
1
-
t
)
]
-
1
2
σ
2
(
T
1
-
t
)
u
1
2
}
,
φ
x
T
2
(
u
2
)
=
exp
⁡
⁡
{
i
u
2
[
x
T
1
+
(
r
-
1
2
σ
2
)
(
T
2
-
T
1
)
]
aaaaaaaaa
-
1
2
σ
2
(
T
2
-
T
1
)
u
2
2
}
.Lemma 3.
Under the risk-neutral measureQ, the bivariate characteristic function is given by
(15)
φ
x
T
1
,
x
T
2
(
u
1
,
u
2
)
=
exp
⁡
⁡
{
i
u
1
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
1
-
t
)
]
-
1
2
u
1
2
σ
2
(
T
1
-
t
)
}
×
exp
⁡
⁡
{
i
u
2
[
x
t
+
(
r
-
1
2
σ
2
)
(
T
2
-
t
)
]
-
1
2
u
2
2
σ
2
(
T
2
-
t
)
}
.The univariate characteristic function is used in Carr and Madan [5] to price a European option using FFT within the Black-Scholes dynamics. By employing the Radon-Nikodym derivatives, the characteristic function under measure Q and the spot price measure Q
S are related as follows:
(16)
φ
x
T
1
S
(
u
1
)
=
e
-
r
(
T
1
-
t
)
-
x
t
E
Q
[
e
i
(
u
1
-
i
)
x
T
1
∣
F
t
]
,
φ
x
T
1
,
x
T
2
S
(
u
1
,
u
2
)
=
e
-
r
(
T
2
-
t
)
-
x
t
E
Q
[
e
i
u
1
x
T
1
+
i
(
u
2
-
i
)
x
T
2
∣
F
t
]
.
Hence, from (10), we consider the following expectations which can be expressed in integral form as follows:
(17)
E
Q
[
1
{
x
T
1
>
h
}
]
=
∫
h
∞
p
(
x
T
1
)
d
x
T
1
,
(18)
E
Q
[
1
{
x
T
1
≥
l
,
x
T
2
≥
k
2
}
]
=
∫
l
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
(19)
E
Q
[
1
{
x
T
1
≥
h
,
x
T
2
≥
k
2
}
]
=
∫
h
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
(20)
E
Q
[
1
{
x
T
1
≥
l
}
]
=
∫
l
∞
p
(
x
T
1
)
d
x
T
1
,
(21)
E
Q
[
1
{
x
T
1
≥
h
}
]
=
∫
h
∞
p
(
x
T
1
)
d
x
T
1
,
where p
(
·
) is the probability density of X under the risk-neutral probability Q. It is clear that integrals (17), (20), and (21) are similar, so are integrals (18) and (19). Hence we first apply the FFT to integral (18). Following Carr and Madan [5], we multiply (18) by dampening factors α
1 and α
2 to ensure square-integrability as follows:
(22)
e
α
1
l
+
α
2
k
2
∫
l
∞
∫
k
2
∞
p
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
,
for
α
1
,
α
2
>
0
.
Then, by applying the Fourier transform, we have the following:
(23)
ψ
(
u
1
,
u
2
)
=
∬
-
∞
∞
e
i
u
1
l
+
i
u
2
k
2
(
e
α
1
l
+
α
2
k
2
∫
l
∞
∫
k
2
∞
p
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
∫
l
∞
∫
k
2
∞
p
∫
l
∞
(
x
T
1
,
x
T
2
)
d
x
T
2
d
x
T
1
)
.
This can be written in closed-form in terms of the characteristic function as follows:
(24)
ψ
(
u
1
,
u
2
)
=
φ
x
T
1
,
x
T
2
[
u
1
-
i
α
1
,
u
2
-
i
α
2
]
(
i
u
1
+
α
1
)
(
i
u
2
+
α
2
)
.
Using the inverse Fourier transform, we can recover integral (18):
(25)
e
-
α
1
l
-
α
2
k
2
(
2
π
)
2
∬
-
∞
∞
e
-
i
u
1
l
-
i
u
2
k
2
ψ
(
u
1
,
u
2
)
d
u
2
d
u
1
.
Then, we numerically compute this by invoking the trapezoid rule (see [11]) which approximates the integral by the following sum:
(26)
e
-
α
1
l
-
α
2
k
2
(
2
π
)
2
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
-
i
u
2
,
m
k
2
ψ
(
u
1
,
j
,
u
2
,
m
)
Δ
2
Δ
1
,
where Δ
1,2 is the distance between the points of the integration grid, u
1
,
j
=
(
j
-
(
N
/
2
)
)
Δ
1
and
u
2
,
m
=
(
m
-
(
N
/
2
)
)
Δ
2
,
for
j
,
m
=
0
,
…
,
N
-
1. In order to progress, we define a grid of size (
N
×
N
) by L
2
=
{
(
l
u
,
k
2
,
p
)
∣
0
≤
u
,
p
≤
N
-
1
} where ω
1,2
>
0 is the distance between the logarithmic critical prices and the logarithmic strikes, respectively, as follows:
(27)
l
u
=
(
u
-
N
2
)
ω
1
,
k
2
,
p
=
(
p
-
N
2
)
ω
2
.
Hence, we can now evaluate the following sum: Γ
(
l
,
k
2
)
=
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
-
i
u
2
,
m
k
2
ψ
(
u
1
,
j
,
u
2
,
m
).By choosingω
1,2
Δ
1,2
=
2
π
/
N this yields the following values of sum Γ
(
l
,
k
2
) on L
2:
(28)
Γ
(
l
u
,
k
2
,
p
)
=
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
u
1
,
j
l
u
-
i
u
2
,
m
k
2
,
p
ψ
(
u
1
,
j
,
u
2
,
m
)
=
(
-
1
)
u
+
p
∑
j
=
0
N
-
1
‍
∑
m
=
0
N
-
1
‍
e
-
i
(
2
π
/
N
)
j
u
-
i
(
2
π
/
N
)
m
p
a
a
a
a
a
a
a
a
a
a
a
a
a
×
[
(
-
1
)
j
+
m
ψ
(
u
1
,
j
,
u
2
,
m
)
]
.
We take the (input) array as X
[
j
,
m
]
=
(
-
1
)
j
+
m
ψ
(
u
1
,
j
,
u
2
,
m
)
,
for
j
,
m
=
0
,
…
,
N
-
1. Hence, the approximation is at (
N
×
N
) different logarithmic critical prices and strike prices such that
(29)
e
-
α
1
l
u
-
α
2
k
2
,
p
(
2
π
)
2
Γ
(
l
u
,
k
2
,
p
)
Δ
2
Δ
1
,
for
0
≤
u
,
p
≤
N
-
1
.
In an analogous way, by replacing l with h, we can approximate integral (19). The two-dimensional FFT can also be reduced to a one-dimensional FFT which involves integrals (17), (20), and (21). Therefore for brevity, we exclude the development for the one-dimensional FFT in this paper.


## 2.2. Monte Carlo Simulation

In this section, we price the extendible call option using Monte Carlo simulation [18]. Given the pricing function of an extendible call option given in (6), the vanilla call option C can be priced at time T
1 with strike price K
2 and expiration time T
2. Thus we only need to compute the asset price until time T
1. Let S
T
1
,
j be the asset price at time T
1 on the jth path, and let w
T
1
,
j be the vanilla call price with strike price K
2 and expiration time T
2 at time T
1 for the jth path; that is, C
(
S
T
1
,
j
,
K
2
,
T
2
-
T
1
). Recall the closed-form solution for a vanilla call option, so we have
(30)
C
=
S
T
1
N
(
e
1
)
-
K
2
e
-
r
(
T
2
-
T
1
)
N
(
e
2
)
,
where
(31)
e
1
=
ln
⁡
⁡
(
S
T
1
/
K
2
)
+
(
r
+
(
1
/
2
)
σ
2
)
(
T
2
-
T
1
)
σ
T
2
-
T
1
,
e
2
=
e
1
-
σ
T
2
-
T
1
.
Therefore, the Monte Carlo estimator for the price of an extendible call option is given by
(32)
EC
(
t
,
x
^
T
1
)
=
e
-
r
(
T
1
-
t
)
n
∑
j
=
1
n
‍
max
⁡
⁡
[
e
x
^
T
1
j
-
K
1
,
w
T
1
,
j
-
A
,
0
]
,
where n is the number of simulations. Suppose we have x
t
=
ln
⁡
S
t with the following process:
(33)
x
t
=
(
r
-
1
2
σ
2
)
d
t
+
σ
d
W
t
.
We discretize the asset path using the Euler scheme as follows:
(34)
x
j
+
1
=
x
j
+
(
r
-
1
2
σ
2
)
Δ
t
+
σ
Δ
W
j
,
t
j
=
j
Δ
t
,
Δ
W
j
=
W
t
j
+
1
-
W
t
j
=
Z
Δ
t
with
Z
~
N
(
0,1
)
.
We also apply the Milstein scheme [19] that is strongly convergent with order one such that
(35)
x
j
+
1
=
x
j
+
(
r
-
1
2
σ
2
)
Δ
t
+
σ
Δ
W
j
+
1
2
σ
2
(
(
Δ
W
)
2
-
Δ
t
)
,
t
j
=
j
Δ
t
,
Δ
W
j
=
W
t
j
+
1
-
W
t
j
=
Z
Δ
t
with
Z
~
N
(
0,1
)
.
We let [
t
=
t
0
<
t
1
<
⋯
<
t
m
=
T
] be a partition of the time interval [
t
,
T
] of m equal segments with length Δ
t
j, where Δ
t
j
=
t
j
+
1
-
t
j and t
j
=
j
T
/
m, for each j
=
1
,
…
,
M. We draw a random sample from a normal distribution with mean 0 and volatility 1 for each j
=
0,1
,
…
,
m to simulate (34) and (35), thus generating a sample path for x
T by simulating x
j for j
=
1
to
j
=
m. This step is repeated to generate many paths. We then estimate the price of the extendible call option with Monte Carlo simulation.


## 3.  Numerical Results

In this section, we present a numerical comparison between the FFT and the Monte Carlo simulation techniques which have been described in previous sections. (The computations were implemented in MATLAB and conducted on an Intel Core 2 Quad processor Q8200 @ 2.33 GHz machine running under Windows 7 Professional with 4 GB RAM).We price the extendible call option for a range of underlying asset prices,S
t
=
80,90,100,110,120, and σ
=
0.25, with the inputs in Table 1. We found that the dampening factor α differs for a one-dimensional and a two-dimensional FFT. For the one-dimensional FFT, we choose α
=
1.75
and
Δ
=
0.3, while for the two-dimensional FFT, α
=
{
0.75,1.75
}
and
Δ
=
0.2. Table 2 provides the CPU time taken by the FFT, Euler scheme, and Milstein scheme to produce the prices for extendible call options. It can be seen that FFT takes less than half a second, on average, approximately 0.0869 seconds, to compute the N
=
2
7 prices in a single computation, while the exact solution computed each of the option prices instantaneously. Euler takes slightly longer than a minute, on average 75.4870 seconds, whereas Milstein scheme takes the longest time to produce results, 311.0502 seconds on average. In Table 3, we document the prices obtained using Monte Carlo simulations with Euler and Milstein schemes. The exact solution is used to compute the relative errors between these two schemes. The error is given by the following:
(36)
ϵ
≈
|
E
C
True
-
E
C
Approx
|
E
C
True
×
100
,
where E
C
True is the Black-Scholes pricing solution and E
C
Approx is the estimation from the simulation. From the table, it can be seen that the relative error for the Euler scheme is smaller than the Milstein scheme. We also plot the price convergence of the Monte Carlo simulation when the underlying asset S
=
90 in Figure 2. From Tables 2 and 3, it can be seen that the Euler scheme is more computationally efficient than the Milstein scheme. According to Platen and Bruti-Liberati [20], the Euler scheme has fewer terms and the same order of weak convergence; hence, apart from efficiency issue, the Milstein scheme also has numerical stability issue which can be worse than that of the Euler scheme. We then compare the extendible call prices computed by the FFT and the exact solution. Table 4 shows the prices and percentage differences taken relative to the FFT prices.Table 1
Inputs to price the extendible options.
Input
Value
Initial strike price,K
1
100
Initial expiration time,T
1
0.5
Extended strike price,K
2
105
Extended expiration timeT
2
0.75
Risk-free rate,r
0.08
Premium,A
1Table 2
CPU time (in seconds) for pricing an extendible call option: FFT, Euler scheme, and Milstein scheme.
Underlying asset,S
FFT
Euler scheme
Milstein scheme
80
0.0821
71.4400
289.2060
90
0.0927
82.0390
310.9170
100
0.1042
81.9460
331.771
110
0.0912
77.6070
359.385
120
0.0643
64.4030
263.972Table 3
Prices and relative error (in percent) for pricing an extendible call option: exact solution, Euler scheme, and Milstein scheme.
Underlying asset,S
Exact solution
Euler price
ϵ
Milstein Price
ϵ
80
1.4264
1.4294
0.2071
1.8025
26.3641
90
4.3154
4.3282
0.2982
4.6685
8.1827
100
9.4228
9.3505
0.7679
9.4874
0.6851
110
16.5358
16.5445
0.0524
16.3863
0.9040
120
25.0506
25.0283
0.0888
24.8268
0.8933Table 4
Prices and relative error (in percent) for pricing an extendible call option: exact solution and FFT.
Underlying asset,S
Exact solution
FFT price
ϵ
80
1.4264
1.4247
0.1193
90
4.3154
4.2908
0.5733
100
9.4228
9.5106
0.9232
110
16.5358
16.4352
0.6121
120
25.0506
25.1067
0.2234Figure 2
Price convergence of Monte Carlo simulation for extendible call options forS
=
90.


## 4. Conclusion

This paper studies the pricing of extendible options under the Black and Scholes [16] model by developing a pricing framework which utilizes the FFT technique. The availability of a closed-form representation of the characteristic function enables us to implement the FFT. The characteristic function for the extendible options involves a univariate and a bivariate characteristic function. Hence, the application of the FFT technique involves a one-dimensional and a two-dimensional FFT algorithm.Additionally, we computed the prices of the extendible options using Monte Carlo simulation and used the closed-form pricing solution as a benchmark to compare the numerical efficiency of the FFT and the Monte Carlo simulation. The FFT produces a range of option prices for a range of strike prices for a given maturity. Moreover, the numerical examples show that the FFT is more computationally efficient than the Monte Carlo simulation.Dias and Rocha [2] show that the framework of extendible options can be applied to model oil prices using mean reversion with jumps. Hence, further development such as pricing extendible options with other stochastic features such as stochastic interest rate can be applied. In the future, we plan to use the FFT to price extendible options under the Heston stochastic volatility model [21].The work presented here can also be extended to a more general case, that is, forn-extensions. The FFT has become well established as an option pricing tool since it was first introduced by Carr and Madan [5]. Fang and Oosterlee [22] have introduced a Fourier-Cosine also known as the COS method to price options, and this method has been shown to be faster under certain circumstances. Hence, another possible development of this work is to price extendible options using the COS method and compare the performance of the two Fourier methods in the case for extendible options.



---

*Source: 831470-2014-06-16.xml*