Discussion Paper No. 117
Computation of Greeks using Binomial Tree

Yoshifumi Muroi and Shintaro Suda
2014/10/8

TOHOKU MANAGEMENT & ACCOUNTING RESEARCH GROUPGRADUATE SCHOOL OF ECONOMICS AND
MANAGEMENT TOHOKU UNIVERSITY

KAWAUCHI, AOBA-KU, SENDAI
980-8576 JAPANComputation of Greeks using Binomial Tree
Yoshifumi Muroi ∗ and Shintaro Suda †Abstract: This paper proposes a new eﬃcient algorithm for the computation of Greeks for options
using the binomial tree. We also show that Greeks for European options introduced in this article are
asymptotically equivalent to the discrete version of Malliavin Greeks. This fact enables us to show that
our Greeks converge to Malliavin Greeks in the continuous time model. The computation algorithmsof Greeks for American options using the binomial tree is also given in this article. There are threeadvantageous points to use binomial tree approach for the computation of Greeks. First, mathematics
is much simpler than using the continous time Malliavin calculus approach. Second, we can constructa simple algorithm to obtain the Greeks for American options. Third, this algorithm is very eﬃcientbecause one can compute the price and Greeks (delta, gamma, vega, and rho) at once. In spite of its
importance, only a few previous studies on the computation of Greeks for American options exist, because
performing sensitivity analysis for the optimal stopping problem is diﬃcult. We believe that our method
will become one of the popular ways to compute Greeks for options.

Keywords: Options, Greeks, Binomial Tree

1

Introduction

Greeks are quantities that represent the sensitivity of the price of derivative securities with respect tochanges in the price of underlying assets or parameters. They are deﬁned by derivatives of the option∗ Graduate School of Economics and Management, Tohoku University, 27-1 Kawauchi, Aoba-Ku, Sendai City 980-8576

e-mail:ymuroi@econ.tohoku.ac.jp tel,FAX:+81-22-795-6308† Mitsubishi UFJ Trust Technology Institute Co., Ltd.

(MTEC), 4-2-6 Akasaka, Minato-Ku, Tokyo 107-0052, Japane-mail:suda@mtec-institute.co.jp

1price function with respect to parameters such as the price of underlying assets, volatility level, andspot interest rate. The computation of these sensitivities is very important for risk management. For an
example on this and similar topics, refer to Hull (2008). The recent development of the Malliavin calculus
approach in ﬁnancial mathematics enables us to compute Greeks for a various kinds of contingent claims.
For a previous research on this topic, see Fournié et al. (1999) and Kohatsu-Higa and Montero (2003)
among others. See also Bernis et al. (2003) for the computational methods of Greeks for exotic options
such as knock-out options and lookback options. In these studies, Greeks were usually represented by
expectation formulas derived from the Malliavin calculus; these expectations are computed using Monte
Carlo simulations. Many previous studies have focused on the computational methods of Greeks for
European options and exotic options such as Asian options. In recent times, several studies on thecomputation of Greeks for American options have been reported. See Bally et al. (2005) for example.
Despite recent advances of the Monte Carlo simulations approach to the pricing of American options, suchas Longstaﬀ and Schwartz (2001), it appears that the backward induction approach is still the natural
approach for computing the price of American options. Because of these reasons, we exploit the binomial
tree approach for the computation of Greeks for American options as well as European options. Many
studies have been carried out on the computation of Greeks using Malliavin calculus and Monte Carlosimulations for European options; however, there are only a few studies on the computation of Greeks
for American options. This is because the computation of the Greeks for American options with Monte
Carlo simulation are diﬃcult. If one uses the Malliavin calculus and Monte Carlo simulations to obtain
Greeks for American options, one has to manage these diﬃculties.
Muroi and Suda (2013) proposed new methods for the computation of Greeks for European options
using the binomial tree methods of Cox et al. (1979) and the discrete Malliavin calculus introduced by
Leitz-Martini (2000) and Privault (2008, 2009). Although the computation of Greeks (delta and gamma)
using the binomial tree has been previously discussed by Pelsser and Vorst (1994), no research havebeen reported on the computational methods of vega and rho. Hence, we propose numerical methods forcomputing rho and vega for options using the binomial tree. Recently, Chung et al. (2010) proved thebinomial delta for plain vanilla European options is converging to delta in the continuous time model. In
this article, we also show Greeks for European options are asymptotically equivalent to discrete version2

of the Malliavin Greeks and the binomial Greeks are converging to Greeks in the continuous time model.
The remainder of this paper is organized as follows. In section 2, we give brief explanations on thebinomial tree methods of Cox et al. (1979). The computational methods of Greeks are introduced in
section 3. Computation algorithms of Greeks for American options are discussed in section 4. The numerical results are given in section 5 followed by the concluding remarks in section 6. In the appendix, wederived the closed form formulas for Greeks for European options and we found these are asymptoticallyequivalent to the discrete version of the Malliavin Greeks.

2Binomial TreeIn this section, we brieﬂy explain the binomial tree model of Cox et al. (1979), which has been explainedin many textbooks such as Hull (2008). The basic idea for computation of Greeks for European options
is also presented. On the basis of these explanations, the computational methods of Greeks using thebinomial tree are explained in detail in the next section.
Introduce independently and identically distributed random variables {ϵi }i=1,...,N on the probability
space (Ω, F, Q), where ϵi is a random variable with probability Q[ϵi =√∆t] = p,√Q[ϵi = − ∆t] = 1−p.

The time step ∆t is ﬁxed as T = N ∆t. We introduce a random walk process, {Wi∆t }i=1,...,N , given byWi∆t = ϵ1 + · · · + ϵi .
One can regard the stochastic process {Wi∆t }i=1,...,N with p = 1/2 as an approximation of the standard
Brownian motion. On the other hand, the stochastic process {Wi∆t }i=1,...,N is generally regarded as an√in the general case (p ̸= 12 ).
approximation of the Brownian motion with drift 2p−1∆t

The binomial tree is a computational method for pricing options on securities whose price process isgoverned by the geometric Brownian motion

dPt = Pt (rdt + σdZt ), P0 = s ,(2.1)

where {Zt }t≥0 is a standard Brownian motion under the risk-neutral measure Q. A binomial tree isconstructed in the following manner. We consider a model with N periods and assume that the maturity
date of the options is ﬁxed as T = N ∆t. If the price of underlying assets at time i∆t (i = 0, 1, . . . , N − 1)3

is given by Si∆t , the price moves to uSi∆t or dSi∆t (d < 1 < u) in the next time period (i + 1)∆t.
The probability of the price moving upward to (uSi∆t ) is p and downward to (dSi∆t ) is 1 − p. In order
to construct the binomial tree so that the expectation and variance are consistent with the geometric
Brownian motion (2.1), we ﬁx the parameters u, d, p as√√u = eσ ∆t , d = e−σ ∆t , p =er∆t − du−dwhere we assume an additional relationship, u = 1/d. Consider European options with a pay-oﬀ function
Φ(·) and a maturity date T . The price of options at time i∆t is denoted by E(x, i∆t) if the price ofunderlying assets is x at time i∆t. The price is given by the backward induction algorithm
E(x, i∆t) = e−r∆t (pE(xu, (i + 1)∆t) + (1 − p)E(xd, (i + 1)∆t)),

E(x, N ∆t) = Φ(x) .

(2.2)

One can derive the closed-form formulas for Greeks for European options using the discrete Malliavin
calculus (See Muroi and Suda (2013)).

3

Computation of Greeks for European Options

New computational methods of Greeks for European options are presented in this section. Althoughcomputational methods of delta and gamma using the binomial tree have already been proposed by Pelsser
and Vorst (1994), computational methods of other Greeks such as vega and rho using the binomial tree
have not been deeply studied, except in the ﬁnite diﬀerence approach. We make a following assumption.
Assumption 3.1 We assume that the pay-oﬀ function Φ(·) is a smooth function.
In fact, this assumption is too strong for use in ﬁnancial mathematics. Therefore, we have to relax thisassumption. This is examined in section 7.2.

3.1

Computation of Delta

In this subsection, we calculate delta, which is used to measure the sensitivity of the option price withrespect to changes in the price of underlying assets. Delta is given by the ﬁrst derivative of the option
value function with respect to the price of underlying assets. Delta for European options is computed as∆E (s, 0)= e−r∆t {p√√∂∂E(seσ ∆t , ∆t) + (1 − p) E(se−σ ∆t , ∆t)}
∂s
∂s

4

= e−r∆t {p

√
√
√
√
d
d
E(seσ ∆t , ∆t)eσ ∆t + (1 − p) E(se−σ ∆t , ∆t)e−σ ∆t } .
dx
dx

(3.3)

∂
Now, we have two approaches to compute ∂z
E(seσz , ∆t)|z=±√∆t . The ﬁrst one is given by
√

√

√
∂
E(seσ ∆t , ∆t) − E(se−σ ∆t , ∆t)√E(seσz , ∆t)|z=±√∆t =+ O( ∆t) ,
∂z
2 ∆t

(3.4)

and the second one is given by∂dE(seσz , ∆t)|z=±√∆t =E(x, ∆t)|x=se±σ√∆t × σseσz |z=±√∆t .
∂zdxdE(x, ∆t)|x=se±σ√∆t :These two formulas yield the approximation formula for the derivative dx√√√
E(seσ ∆t , ∆t) − E(se−σ ∆t , ∆t)d1√√
E(x, ∆t)|x=se±σ√∆t =+ O( ∆t) .
dx
2 ∆t
σse±σ ∆t

(3.5)

Substituting this formula in (3.3) leads to the following relation
√
√
√
√
√
e−r∆t E(seσ ∆t , ∆t) ∆t + E(se−σ ∆t , ∆t)(− ∆t)∆E (s, 0) =+ O( ∆t) .
sσ∆t
2

(3.6)

Taylor expansion yields the approximation formula for p, and it is given by

p=

1 1 √
+ µ ∆t + O(∆t3/2 ) ,
2 2

(3.7)

2

where µ is a constant given by µ = σ1 (r − σ2 ) . Let us compute the expectation E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)].
The above approximation formula leads to√√√√E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] = pE(seσ ∆t , ∆t)( ∆t − µ∆t) + (1 − p)E(se−σ ∆t , ∆t)(− ∆t − µ∆t)
√
√
√
√
∆t
− ∆t
σ ∆t
=[
E(se
, ∆t) +E(se−σ ∆t , ∆t)] + O(∆t2 ) .
(3.8)
2
2

Note that we have used the identities√√√
1(E(seσ ∆t , ∆t) + E(se−σ ∆t , ∆t)) = E(s, ∆t) + O( ∆t)
2
√
√
√
√
∂
E(seσ ∆t , ∆t) − E(se−σ ∆t , ∆t) = 2 ∆t E(seσz , ∆t)|z=0 + O(∆t) = O( ∆t)∂zto derive the last equality. This formula allows further computation of (3.6). It is given by∆E (s, 0) =
≈

√
e−r∆t
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] + O( ∆t)
sσ∆t
e−r∆t
S
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] ≡ ∆M
E (s, 0)
sσ∆t
5

(3.9)

We call this MS delta, which is the one-step version of discrete Malliavin delta given in Muroi and Suda
(2013). Actually, we can show stronger result,S∆M
E (s, 0) = ∆E (s, 0) + O(∆t) ,(3.10)

if we assume smoothness to the pay-oﬀ function. This is shown in subsection 7.1 in Appendix. If the
pay-oﬀ function is not smooth, we can show,
√
S
BS
∆M
E (s, 0) = ∆E (s, 0) + O( ∆t) .
See subsection 7.2 in Appendix. Discrete Malliavin Greeks for European options using an N -steps binomial tree were obtained by Muroi and Suda (2013); they used the discrete Malliavin derivatives introduced
by Leitz-Martini (2000). See also Privault (2008, 2009). Discrete Malliavin delta ∆D
E is given by∆D
E =e−rN ∆tE[Φ(seσWN ∆t )(WN ∆t − µN ∆t)] .
sσN ∆tAs discussed in Appendix 7.2, MS delta and discrete Mallaivin delta is actually equivalent. Because
Muroi and Suda (2013) used the discrete version of the Malliavin calculus approach, one cannot use their
method to derive Greeks for American options. In our study, we exploit the stepwise approach to derive
Greeks for American options. See section 4 for a detailed discussion on the computation of Greeks for
American options.
dE(x, ∆t)|x=se±σ√∆t , we can obtain
If we exploit another approximation formula for the derivative dx

another approximation formula for delta. If we use a more direct formula for the derivative,
√

√

√

√

√
d
E(se∓σ ∆t , ∆t) − E(se±σ ∆t , ∆t)E(seσ ∆t , ∆t) − E(se−σ ∆t , ∆t)√√√,
E(se±σ ∆t , ∆t) ≈≈dx
2sσ ∆tse∓σ ∆t − se±σ ∆t(3.11)

we get another approximation formula for delta for European options in the one-period time model.
Substituting (3.11) in (3.3) yields ∆Hull, and it is given as
E
E(se
∆Hull
(s, 0) ≈
E

√
σ ∆t

√

, ∆t) − E(se−σ ∆t , ∆t)√.
2sσ ∆tThis is the delta introduced by Hull (2008). This fact reveals that the two diﬀerent approximation√dformulas for dx
E(se±σ ∆t , ∆t) given by (3.5) and (3.11) yield two diﬀerent approximation formulasfor delta 1 . Our computational results indicate that MS delta converges a little bit faster than delta
1 Actually, MS delta is deeply related to Hull’s delta.

√S−r∆t ∆Hull (s, 0) + O( ∆t),
We have a relation, ∆ME (s, 0) = e
Ebecause we have formulas (3.8) and (3.9).

6

S
D
introduced by Hull (2008). Moreover, as shown in the appendix, the formula, ∆ME (s, 0) = ∆E (s, 0), is

actually satisﬁed. This means that MS delta is more natural representation of delta than Hull’s delta.

3.2

Computation of Gamma

In this subsection, we calculate gamma, which is used to measure the sensitivity of delta with respect tochanges in the price of underlying assets. Gamma is given by the second derivative of the option valuefunction with respect to the price of underlying assets. The pricing algorithm for European options (2.2)
yields delta for European options:∆E (s, 0)= e−r∆t {p

√
√
√
√
∂
∂
E(seσ ∆t , ∆t)eσ ∆t + (1 − p) E(se−σ ∆t , ∆t)e−σ ∆t }
∂x
∂x
√

√

√

√

= e−r∆t {p∆E (seσ ∆t , ∆t)eσ ∆t + (1 − p)∆E (se−σ ∆t , ∆t)e−σ ∆t } .
Applying chain rule to E(seσz , ∆t) leads to√∂E∆E (se±σ ∆t , ∆t) = ∂z

(seσz , ∆t)|z=±√∆t .
σseσzDelta for European options is given by∆E (s, 0) =

e−r∆t
∂E
∂E
{p[ (seσz , ∆t)]z=√∆t + (1 − p)[ (seσz , ∆t)]z=−√∆t } .
σs
∂z
∂z

Taking derivative with respect to the price of underlying assets yields gamma for European options as1
ΓE (s, 0) = − ∆E (s, 0)
s
e−r∆t
∂
∂
+
{p (∆E (seσz , ∆t)eσz )|z=√∆t + (1 − p) (∆E (seσz , ∆t)eσz )|z=−√∆t } .
σs∂z∂z
Using the approximation formula
√

√

√

√

√
∆E (seσ ∆t , ∆t)eσ ∆t − ∆E (se−σ ∆t , ∆t)e−σ ∆t∂√∆E (seσz , ∆t)eσz |z=±√∆t =+ O( ∆t)∂z
2 ∆t
allows further calculation of gamma as
√
e−r∆t
e−r∆t
ΓE (s, 0) = − 2
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] +E[∆E (seσϵ1 , ∆t)eσϵ1 (ϵ1 − µ∆t)] + O( ∆t)
s σ∆t
σs∆t
−r∆t
√
e−r∆t
e
S
σϵ1
= − 2
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] +E[∆M, ∆t)eσϵ1 (ϵ1 − µ∆t)] + O( ∆t)
E (se
s σ∆t
σs∆t
−r∆t
e
e−r∆t
S
σϵ1
≈ − 2
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] +E[∆M, ∆t)eσϵ1 (ϵ1 − µ∆t)]
E (se
s σ∆t
σs∆t
S
≡ ΓM
E (s, 0)(3.12)
7

We call this formula MS gamma. This formula is valid only if the pay-oﬀ function Φ(·) is a smooth√√Sfunction. In this case, the order of the error term is O( ∆t), i.e. ΓE (s, 0) = ΓME (s, 0) + O( ∆t).
Second equality in (3.12) is shown by using the formula (3.10). We will also prove that MS gamma is
asymptotically equivalent to the discrete Malliavin Gamma in the appendix.

3.3

Computation of Vega

The computational method of vega for European options is presented in this subsection. Vega is thesensitivity of the option pricing formula with respect to changes in volatility level, σ. It appears that thecomputational methods of vega and rho using the binomial tree have not yet been considered seriously,
except for the ﬁnite diﬀerence approach. See Hull (2008) for computation of vega using the ﬁnite diﬀerence
approach with the binomial tree. Let us assume that the price for underlying assets at time i∆t isgiven by si∆t . The price and vega for European options at time i∆t are denoted by E(si∆t , i∆t; σ) andVE (si∆t , i∆t; σ), respectively. Vega for European options is given by taking derivatives to the pricing
formula (2.2):
VE (si∆t , i∆t; σ) =

e−r∆t

√
√
∂
{pE(si∆t eσ ∆t , (i + 1)∆t; σ) + (1 − p)E(si∆t e−σ ∆t , (i + 1)∆t; σ)}∂σ= V1i + V2i + V3i ,where V1i , V2i , V3i are given by
V1i
V2i

√
√
∂p
∂(1 − p)E(si∆t eσ ∆t , (i + 1)∆t; σ) +E(si∆t e−σ ∆t , (i + 1)∆t; σ)}
∂σ
∂σ
√
√ √
∂
= e−r∆t {p E(si∆t eσ ∆t , (i + 1)∆t; σ)si∆t eσ ∆t ∆t
∂x
√
√
√
∂
+(1 − p) E(si∆t e−σ ∆t , (i + 1)∆t; σ)si∆t e−σ ∆t (− ∆t)}∂x

= e−r∆t {= e−r∆t E[∆E (si∆t eσϵi+1 , (i + 1)∆t; σ)si∆t eσϵi+1 ϵi+1 ]V3i= e−r∆t E[VE (si∆t eσϵi+1 , (i + 1)∆t; σ)] .

The derivative of p with respect to σ is given by√√∂
2 − er∆t (eσ ∆t + e−σ ∆t ) √
2r √
1
√
√
p=
∆t = − (1 + 2 ) ∆t + O(∆t3/2 ) .
∂σ
4σ(eσ ∆t − e−σ ∆t )2This result and the approximation formula (3.8) lead to the approximation formula for V1i :
√
√
√
√
σ ∆t
−σ ∆t
1
2r
E(s
e
,
(i
+
1)∆t;
σ)
∆t
+
E(s
e
, (i + 1)∆t; σ)(− ∆t)
i∆t
i∆t
i
−r∆t
V1 = − (1 + 2 )e2σ28+O(∆t3/2 )=E[e−r∆t {−(1 +2r ϵi+1 − µ∆t)}E(si∆t eσϵi+1 , (i + 1)∆t; σ)] + O(∆t3/2 ) .
σ2
2

For further computation of V2i , one needs delta at time (i + 1)∆t. Moreover, it is necessary to computevega for European options at time (i + 1)∆t to evaluate V3i . This implies that one has to use thebackward induction algorithm to compute vega for European options. Vega at the maturity date is given
by VE (sN ∆t , N ∆t; σ) = 0. These results are combined to form the computational formulas for vega:
VE (si∆t , i∆t; σ)E[e−r∆t {−(1 +=2r ϵi+1 − µ∆tE(si∆t eσϵi+1 , (i + 1)∆t)])}σ22+E[e−r∆t ∆E (si∆t eσϵi+1 , (i + 1)∆t; σ)si∆t eσϵi+1 ϵi+1 ]+E[e−r∆t VE (si∆t eσϵi+1 , (i + 1)∆t; σ)] + O(∆t3/2 ) .

(3.13)

Delta and vega at time (i + 1)∆t are substituted by MS delta and vega to get MS vega at time i∆t,respectively. In other words, MS vega at time i∆t is deﬁned recursively by the backward algorithm,MS
VE(si∆t , i∆t; σ) ≡ E[e−r∆t {−(1 +2r ϵi+1 − µ∆t)}E(si∆t eσϵi+1 , (i + 1)∆t)]
σ2
2

S
σϵi+1
+E[e−r∆t ∆M, (i + 1)∆t; σ)si∆t eσϵi+1 ϵi+1 ]E (si∆t eMS+E[e−r∆t VE(si∆t eσϵi+1 , (i + 1)∆t; σ)] .

(3.14)

Note that we formally assume V2N −1 = 0 if the pay-oﬀ function Φ(·) is not smooth one. Although theorder for the error term of equation (3.13) is equal to O(∆t3/2 ), the relation√MS(s, 0; σ) + O( ∆t)VE (s, 0; σ) = VE

(3.15)

has to be satisﬁed. This is shown in (7.21) in appendix. As an alternative approach, the ﬁnite diﬀerence
approach is used to obtain vega for practical purposes. In many cases, it works well: however, in some itdoes not, for example one cannot obtain a stable estimator of vega for digital options.

3.4

Computation of Rho

The computational method of rho for European options is discussed in this subsection. Rho measures
the sensitivity of the option price with respect to changes in the spot interest rate level, r. It is deﬁned
by the ﬁrst derivative of the option value function with respect to the spot interest rate, r. The price and9rho for European options are denoted by E(si∆t , i∆t; r) and ρE (si∆t , i∆t; r), respectively, if the price ofunderlying assets at time i∆t is given by si∆t . Simple calculation yields
ρE (si∆t , i∆t; r) = ρi1 + ρi2 + ρi3 ,where ρi1 , ρi2 , ρi3 are given byρi1= −∆te−r∆t E[E(si∆t eσϵi+1 , (i + 1)∆t; r)]ρi2= e−r∆t {ρi3= e−r∆t E[ρE (si∆t eσϵi+1 , (i + 1)∆t; r)] .

√√∂(1 − p)∂p
E(si∆t eσ ∆t , (i + 1)∆t; r) +E(si∆t e−σ ∆t , (i + 1)∆t; r)}∂r∂ri
The derivative ∂p∂r in ρ2 is given by

√
er∆t
∆t
∂
√
p= √
+ O(∆t3/2 ) .
∆t =
−σ
σ
∆t
∆t
∂r
2σ
−e
e
This relation and the formula (3.8) leads to the further calculation,
ρi2

=
=

√
√
√
√
e−r∆t
∆t
− ∆t
σ ∆t
{
E(si∆t e, (i + 1)∆t; r) +E(si∆t e−σ ∆t , (i + 1)∆t; r)} + O(∆t3/2 )
σ
2
2
e−r∆t
E[E(si∆t eσϵi+1 , (i + 1)∆t; r)(ϵi+1 − µ∆t)] + O(∆t3/2 ) .
σ

In order to compute ρi3 , one needs to compute rho at time (i + 1)∆t, i.e., ρE (si∆t eσϵi+1 , (i + 1)∆t; r).
This implies that one has to use the backward induction approach to evaluate rho. These results are
combined to form the computational formulas of rho for European options:ρE (si∆t , i∆t; r) =e−r∆t E[(

ϵi+1
µ∆t
−
− ∆t)E(si∆t eσϵi+1 , (i + 1)∆t; r)σσ+ρE (si∆t eσϵi+1 , (i + 1)∆t; r)] + O(∆t3/2 ) .

(3.16)In order to compute MS rho at time i∆t, rho at time (i + 1)∆t is substituted by MS rho at time (i + 1)∆t.
In other words, MS rho is deﬁned recursively by the backward algorithm,
S
−r∆t
ρM
E[(
E (si∆t , i∆t; r) ≡ e

ϵi+1 µ∆tSσϵi+1−− ∆t)E(si∆t eσϵi+1 , (i + 1)∆t; r) + ρM
, (i + 1)∆t; r)] .
E (si∆t eσ
σ

Although the order for the error term of equation (3.16) is equal to O(∆t3/2 ), the relation√SρE (s, 0; r) = ρME (s, 0; r) + O( ∆t)
has to be satisﬁed. This is also shown in appendix.
104

Computation of Greeks for American Options

We suggest a new algorithm for computation of Greeks for American options. An American option is acontingent claim that its holder can exercise at any time before its maturity date. Consider an American
option with a pay-oﬀ function Φ(·) and a maturity date T = N ∆t. If the price of underlying assets attime i∆t is given by x, the price of these options at time i∆t is denoted by A(x, i∆t). The price ofAmerican options is given by the backward induction algorithm:
A(x, i∆t) = max{e−r∆t (pA(xu, (i + 1)∆t) + (1 − p)A(xd, (i + 1)∆t)), Φ(x)} ,where the price of options at the maturity date is given by A(x, N ∆t) = Φ(x). Introduce new sets,
Si = {x|x ∈ (0, ∞), A(x, i∆t) = Φ(x)} and Ci = (0, ∞) \ Si(i = 0, 1, . . . , N ). The continuous regionand stopping region for American options are deﬁned by C = ∪N
i=0 CiandS = ∪N
i=0 Si , respectively.

We will present an intuitive way for the computation of Greeks for American options. If the price ofunderlying assets at the initial time satisﬁes s ∈ S0 , Greeks are simply the derivative of the pay-oﬀfunction. Sensitivity for American options at the initial time respect to a parameter θ, is computed as∂θ A(s, 0) = ∂θ Φ(s)if the function Φ(·) is diﬀerentiable at s ∈ S0 . (Second derivatives such as gamma is given by ∂ss A(s, 0) =∂ss Φ(s).) If the price of underlying assets at the initial time satisﬁes s ∈ C0 , Greeks are also given by thederivative of the price function of American options, and it is given by2∂θ A(s, 0) =∂ −r∆t[e{pA(su, ∆t) + (1 − p)A(sd, ∆t)}] .
∂θ(4.17)

Notice that the formula (4.17) has a same functional form to the European options Greeks. Delta, forexample, is given byS∆M
A (s, 0) = ∂s Φ(s)if s ∈ S0


−r∆t

σϵ1
 e
, ∆t)(ϵ1 − µ∆t)] if s ∈ C0 .
sσ∆t E[A(se2 Strictly speaking, one cannot use our computational formula for delta if the initial price is on the early exercise boundary,and one cannot use our gamma, vega, and rho formulas if the nodes on the binomial tree are on the early exercise boundary.
However, numerical results show that our formula works very well when we compute Greeks for American put options withthe pay-oﬀ Φ(x) = (K − x)+ if the delta is ﬁxed at ∆A (s, 0) = −1 on the stopping region and the boundary. This is because
we have a smooth-ﬁt condition in the continuous model, and our model is an approximation of the Black and Scholes model.

11One can compute other Greeks using same method. Numerical demonstrations are shown in Section 5.
The extended binomial tree of Pelsser and Vorst (1994) is one of the most suitable alternative methods toderive delta and gamma using binomial trees. One can eﬃciently and accurately derive Greeks eﬃciently
and accurately as discussed in Section 5. However, one cannot apply this method to derive vega andrho.

5

Numerical Results

In this section, we demonstrate the numerical results for the new computational methods of Greeks thatwere introduced in previous sections. In order to check the eﬀectiveness of our approach, Greeks for
American put options are computed by the newly proposed approach and the ﬁnite diﬀerence approach.
We also demonstrate the extended binomial tree approach of Pelsser and Vorst (1994) to compute delta
and gamma. It is well known that the extended binomial tree approach of Pelsser and Vorst (1994) yields
very accurate and fast algorithms to compute delta and gamma for options. On the other hand, the ﬁnite
diﬀerence approach is a very popular approach for computing vega and rho, as discussed in Hull (2008).
We also compare to the existent other tree methods for computations of Greeks.
Figures 2, 3, 4 and 6 plot the values of Greeks (delta, gamma, vega, and rho, respectively) for
American put options computed using our approach and the ﬁnite diﬀerence approach. MS Greeks,
extended binomial Greeks (EB Greeks) calculated by the extended binomial tree of Pelsser and Vorst(1994), Greeks introduced by Hull (Hull’s Greeks) are also plotted in Figures 2 and 3. The extended(N-step) binomial tree is a N + 2 binomial tree starting from −2∆t, as shown in Figure 1. The EB deltaand EB gamma are given by
∆EB
A

=

ΓEB
A

=

A(s(2, 0), 0) − A(s(−2, 0), 0)s(2, 0) − s(−2, 0)A(s(2,0),0)−A(s(0,0),0)− A(s(0,0),0)−A(s(−2,0),0)s(2,0)−s(0,0)
s(0,0)−s(−2,0)s(2, 0) − s(−2, 0)

where s(i, j) = s ∗ ui . These results are computed using binomial trees with N = 1, 2, · · · , 50 steps for
one year. The parameter values assumed for these numerical studies weres = 100,

σ = 0.3,

r = 0.05,

12

K = 100,

T =1(year) .

The MS Greeks (delta, gamma, vega, and rho) are represented by the real lines, and two kinds of dottedlines represent the ﬁnite diﬀerence Greeks (FD Greeks) and EB Greeks. The horizontal lines in Figures
2 and 3 are the EB delta and EB gamma computed by the extended binomial tree with 100, 000 steps for
one year. The horizontal lines in Figures 4 and 6 are FD vega and FD rho computed using the binomial
trees with 100, 000 steps for one year. Because we use very ﬁne meshes for the computation of thesehorizontal lines, these numerical results are expected to be very accurate. If the initial underlying asset
price s is in the continuous region, i.e. s ∈ C0 , the FD delta and gamma are given by
D
∆F
A

=

D
ΓF
A

=

A(s + ∆s, 0) − A(s − ∆s, 0)2∆sA(s + ∆s, 0) − 2A(s, 0) + A(s − ∆s, 0)(∆s)2

and the FD vega and rho are given byFD(s, 0; σ) =VADρF
A (s, 0; r) =A(s, 0; σ + ∆σ) − A(s, 0; σ − ∆σ)2∆σA(s, 0; r + ∆r) − A(s, 0; r − ∆r)
2∆r

where ∆s, ∆σ, and ∆r are small parameters. We take ∆s = s×∆h, ∆σ = σ∆h, and ∆r = r×∆h, (∆h =10−3 ) for our computations. Figure 2 shows that MS delta converges much faster than the FD delta.
As shown, the oscillation phenomenon is observed for FD delta, whereas MS delta converges very fast.
Figure 3 shows that FD gamma is not stable, and we do not recommend the use of the ﬁnite diﬀerence
approach to compute gamma. This phenomenon has already been reported by Pelsser and Vorst (1994).
Moreover, MS delta and MS gamma converge slightly faster than EB delta, Hull’s delta, EB gamma,
and Hull’s gamma. Figure 4 reveals that MS vega converges slower than FD vega. However, this is nota universal result. MS vega and FD vega for American put options are plotted in Figures 5 with strikeprices of K = 105 (in-the-money case). The oscillation phenomenon for FD vega is observable for the
options with the strike price K = 105. The behaviors of MS rho and FD rho demonstrated in Figure 6
are almost the same, and we found this to be a universal relationship in our numerical experience. It
is important to note the backward induction algorithm needs to be used only once to obtain MS rho,
whereas it has to be used twice to obtain FD rho. Hence, the computational time for MS rho is expected
to be shorter than that for FD rho. Table 1 lists the computational time and results for MS rho and13s(4,2)s(3,1)s(2,2)s(2,0)s(1,1)

s(1,‑1)s(0,‑2)s(0,2)s(0,0)s(‑1,1)s(‑1,‑1)s(‑2,2)s(‑2,0)s(‑3,1)s(i,j)=s*u^is(‑4,2)

Figure 1: Extended Binomial Tree
Table 1: Computational time for Rho
MS rhoFD rhoValue-34.8461-34.8460

Time

14.77 (s)18.78 (s)FD rho computed by a binomial tree with 10, 000 steps 3 . It was found that the computational time forMS rho was about 20 % shorter even though the computational results obtained were almost the same.
Hence, computing MS rho rather than the FD rho is more advantageous.
Figures 7 to 10 present Greeks (delta, gamma, vega, and rho, respectively) for American put options as
a function of the price of underlying assets. The price range of underlying assets is from 50 to 200. Other
parameters used for these numerical studies are same as those used in the previous numerical studies.
Figures 7 and 8 plot the MS, FD, and EB delta and gamma computed using the 100 step binomial trees,
respectively. The curves of all the MS Greeks and EB Greeks are very smooth, whereas those of FD
delta and FD gamma are unstable. Figures 9 and 10 plot the MS and FD vega and rho using a 100 step
binomial tree, respectively. As shown in Figure 9, the shape of MS vega is very smooth, whereas theoscillation phenomenon is observed for FD vega. The oscillation phenomenon for FD vega is especiallystrong when the strike price is higher than the initial price of underlying assets. Figure 10 reveals thatthe numerical results of MS rho and FD rho are almost same.

3 It is enough to use binomial trees with 100 steps to obtain Greeks. Then, one can compute in an instant.

14

Finally, we compared our new methods with other existing tree methods. We compare MS delta withdelta for European options computed by other kinds of binomial tree, namely tree methods introducedby Chung and Shackelton (2002), Tian (1993), and Leisen and Reimer (1996). These are summarized
in Figure 11. In order to obtain delta using the binomial trees of Chung and Shackelton (2002) andTian (1993), we employed the extended binomial tree approach. On the other hand, we used the ﬁnite
diﬀerence approach to the binomial tree for Leisen and Reimer, because we wanted to implement simple
calculations. Leisen and Reimer (1996) introduced a new kind of binomial tree, which computes the priceof options eﬃciently. They construct two kinds of trees using two diﬀerent transform formulas. Notethat because no signiﬁcant diﬀerence is observed in two methods of Leisen and Reimer (1996), we used
”Method-1” described in their article. As shown in Figure 11, Greeks calculated by trees introducedby Chung and Shackelton (2002) and Tian (1993) converges to the real value smoothly, however, MSdelta converges faster than these methods. Delta computed by the tree introduced by Leisen and Reimer(1996) converges considerably fast, if one uses trees with odd steps.. It should be pointed out that it is
not easy to compute vega and rho by Leisen and Reimer’s binomial tree.

6

Conclusion

This paper presented new computational methods of Greeks using the binomial tree. There are two
important results in this paper. First, we obtain a very eﬃcient algorithm to evaluate Greeks. It isespecially eﬃcient to compute Greeks for American options. Although many studies have been conducted
for the computation of Greeks for European options, few papers have examined the computation of Greeks
for American options. We introduce the binomial tree approach to overcome these problems and conﬁrm
its eﬀectiveness for computing Greeks for American options very quickly and accurately. Numerical results
indicate that Greeks converge faster when computed using our method than when computed using theextended binomial tree approach of Pelsser and Vorst (1994). Second, we show that Greeks computed by
our algorithm converge to the Greeks in the continuous time model. We also showed the relation between
MS Greeks and discrete Malliavin Greeks. We are now preparing an article on computations of Greeks inthe jump diﬀusion models using the binomial tree approach (Muroi and Suda (2013) and Suda and Muroi

15(2013)). One can use our approach to compute Greeks only if one can construct a recombining binomial
tree. In contrast, it is not clear how to compute Greeks using our approach if one cannot construct a
recombining binomial tree (e.g. stochastic volatility model). We also mention that sensitivity analysis
using the Malliavin calculus is important in many other research areas; see Privault and Wei (2004) and
Loisel and Privault (2009) for applications in ruin theory, and Corcuera and Kohatsu=Higa (2011) for
applications in statistics. Applications for the density estimation method have also been examined by
Privault and Wei (2007). We hope that our approach is applied to other research areas. We also believethat computation of Greeks for exotic options using the binomial tree approach is remained as important
topics as future researches.

7

Appendix: Closed-Form Formulas for Option Greeks

There are two aims of this appendix; Firstly, we show that the order of error term for MS delta andvega are given by (3.10) and (3.15). This is presented in section 7.1. Secondly, we prove that MS greeks
converge to greeks for continous time Black and Scholes model. This is shown in section 7.2. Closed-form
expectation formulas for MS Greeks (delta, gamma, vega, rho) for European options are investigated inAppendix 7.2. We found that MS Greeks are approximations of the discrete version of the Malliavin
Greeks in the continuous time model and these results indicate that MS Greeks converge to Greeks for acontinuous time model (Black and Scholes model) when we take a limit, ∆t → 0.

7.1

Error Terms for MS Greeks

Let us assume that the pay-oﬀ function Φ(·) is a smooth function in this subsection.
[Delta]In this subsection, error terms of the MS Greeks are computed precisely. We still assume that the pay-oﬀfunction Φ(·) is a smooth function in this subsection. First we compute a higher order term for the error
term for MS delta. Taylor expansion,√E(se∓σ ∆t , ∆t) =√∂E(seσz , ∆t)(∓2 ∆t)|z=±√∆t
∂z
√
1 ∂2
+
E(seσz , ∆t)(∓2 ∆t)2 |z=±√∆t + O(∆t3/2 ) ,
2
2 ∂z
√

E(se±σ ∆t , ∆t) +16

yields the higher order term for equality (3.4). It is given by
√

√

√
∂2
∂
E(seσ ∆t , ∆t) − E(se−σ ∆t , ∆t)√E(seσz , ∆t)|z=±√∆t =± 2 E(seσz , ∆t) ∆t|z=±√∆t + O(∆t) .
∂z
∂z
2 ∆t
This leads to the higher order expansion formula for the delta given by (3.6):∆E (s, 0) =
+

√
√
√
√
e−r∆t E(seσ ∆t , ∆t) ∆t + E(se−σ ∆t , ∆t)(− ∆t)
sσ∆t
2
2
−r∆t
√
∂
e
∂2
(p 2 E(seσz , ∆t)|z=√∆t − (1 − p) 2 E(seσz , ∆t)|z=−√∆t ) ∆t + O(∆t) .(7.18)
sσ
∂z
∂z

The second term in (7.18) is calculated as
p

∂2
∂2
σz
√
E(se
,
∆t)|
−
(1
−
p)
E(seσz , ∆t)|z=−√∆t
z= ∆t
∂z 2
∂z 2
√
√
∂2
∂3
= (µ 2 E(seσz , ∆t)|z=0 + 3 E(seσz , ∆t)|z=0 ) ∆t + O(∆t) = O( ∆t) ,∂z∂z

where we have used the approximation (3.7). This result enables us to compute delta given by (7.18). It
is computed as√√√√e−r∆t E(seσ ∆t , ∆t) ∆t + E(se−σ ∆t , ∆t)(− ∆t)∆E (s, 0) =+ O(∆t) .
sσ∆t2

Applying the formula (3.8), the higher order expansion for delta is ﬁnally given byS∆E (s, 0) = ∆ME (s, 0) + O(∆t) .

(7.19)√
Error term is not identical to O( ∆t), if we assumed smoothness to the pay-oﬀ function. This result isused to derive the closed-form expression for MS gamma and MS vega.
[Vega]
Insert (7.19) to the formula (3.13) leads to
VE (si∆t , i∆t; σ)=E[e−r∆t {−(1 +2r ϵi+1 − µ∆t)}E(si∆t eσϵi+1 , (i + 1)∆t)]
σ2
2

S
σϵi+1
+E[e−r∆t ∆M, (i + 1)∆t; σ)si∆t eσϵi+1 ϵi+1 ]E (si∆t e+E[e−r∆t VE (si∆t eσϵi+1 , (i + 1)∆t; σ)] + O(∆t3/2 ) .
This formula yields
VE (si∆t , i∆t; σ) = V1i (si∆t , i∆t; σ) + V2i (si∆t , i∆t; σ) + E[e−r∆t VE (si∆t eσϵi+1 , (i + 1)∆t; σ)] + O(∆t3/2 ) ,(7.20)
17where new variables V1i and V2i are deﬁned by
2r ϵi+1 − µ∆t)}E(si∆t eσϵi+1 , (i + 1)∆t)|Fi∆t ]σ22V1i (si∆t , i∆t; σ)=E[e−r∆t {−(1 +V2i (si∆t , i∆t; σ)

=

S
σϵi+1
E[e−r∆t ∆M, (i + 1)∆t; σ)si∆t eσϵi+1 ϵi+1 |Fi∆t ]E (si∆t e

V2N −1 (s(N −1)∆t , (N − 1)∆t; σ)=E[e−r∆t {(i ̸= N − 1)dΦ(x)|x=s(N −1)∆t eσϵN }s(N −1)∆t eσϵN ϵN |F(N −1)∆t ]dx

If the pay-oﬀ function Φ(x) is not smooth, we formally deﬁne V2N −1 (s(N −1)∆t , (N −1)∆t; σ) = 0 , although
we assume that the pay-oﬀ function Φ(x) is smooth in this subsection. If the pay-oﬀ function Φ(x) issmooth, vega given by (7.20) is further computed as
VE (si∆t , i∆t; σ) = V1i (si∆t , i∆t; σ)+V2i (si∆t , i∆t; σ)+E[e−r∆t VE (si∆t eσϵi+1 , (i+1)∆t; σ)|Fi∆t ]+O(∆t3/2 ) .
Vega for the binomial tree model at the initial time is given by
VE (s, 0; σ)=V10 (s, 0; σ) + V20 (s, 0; σ) + O(∆t3/2 ) + e−r∆t E[VE (seσW∆t , ∆t; σ)]= V10 (s, 0; σ) + V20 (s, 0; σ) + O(∆t3/2 )+ e−r∆t E[V11 (seσW∆t , ∆t; σ) + V21 (seσW∆t , ∆t; σ) + O(∆t3/2 ) + e−2∆t VE (seσW2∆t , 2∆t; σ)]
= ···
=

N
−1
∑

√
e−ri∆t E[V1i (seσWi∆t , 0; σ) + V2i (seσWi∆t , i∆t; σ)] + O( ∆t)

i=0

√
MS
(s, 0; σ) + O( ∆t) .
= VE(7.21)[Rho]
Because rho at the maturity date is equal to 0, rho given by the formula (3.16) is further calculated as

ρE (s, 0; r)=e−r∆t E[(

ϵ1 − µ∆t− ∆t)E(seσϵ1 , ∆t; r) + ρE (seσϵ1 , ∆t; r) + O((∆t)3/2 )]
σ

= ···
N
∑

√
ϵi − µ∆t− ∆t)E(seσWi∆t , ∆t; r)] + O( ∆t)
σ
i=1
√
S
= ρM
E (s, 0; r) + O( ∆t) .
=e−ri∆t E[(

Note that smoothness of the pay-oﬀ function is not used to derive the fromula (7.22).

18(7.22)

7.2

Convergence of MS Greeks to Black and Scholes ModelIn this subsection, convergence of the MS Greeks are shown. We assume that the pay-oﬀ function Φ(·)
is not a necessarily smooth function in this subsection. We remove the assumption of the smoothness
for the pay-oﬀ function given in Assumption 3.1. In section 7.2, instead of Assumption 3.1, we make thefollowing assumption.
Assumption 7.1 We assume that the pay-oﬀ function, Φ(·), is a function in the class K. K is theclass of real-valued functions on R that satisfy the following conditions: (i) ϕ(·) is piecewise C (2) , (ii)at each x, the function ϕ(·) satisﬁes ϕ(x) = 12 (ϕ(x+) + ϕ(x−)), and (iii) ϕ, ϕ′ , and ϕ′′ are polynomial
bounded. We assume that f (·) is a function in the class K. Then E[e−rT f (ST )] = E[e−rT f (seσWN ∆t )] →σ2

E[e−rT f (se(r− 2 )T +σZT )] = E[e−rT f (PT )] (N → ∞) is satisﬁed.
For example, the pay-oﬀ function for European call/put options is included in class K. Therefore, if
we demonstrate the convergence of MS Greeks to Greeks for the Black and Scholes model under thisassumption, we can say that MS Greeks for European call/put option is an approximation formula for
Greeks in the continuous time model.
[Delta]MS delta given by the formula (3.9) is further calculated asS∆M
E (s, 0) =

=

e−rT
e−r∆t
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] =E[Φ(seσWT )(ϵ1 − µ∆t)]
sσ∆t
sσ∆t
e−rT
E[Φ(seσWT )(ϵi − µ∆t)] ,sσ∆twhere we used the fact that ϵ1 , . . . , ϵN is an i.i.d. sequence to deduce the last equality. This yields
1 ∑ e−rTE[Φ(seσWT )(ϵi − µ∆t)]
N i=1 sσ∆t
NS∆M
E (s, 0) ==e−rTlog(ST /s) − µσT
e−rT
E[Φ(seσWT )(WT − µT )] =E[Φ(ST )
],
sσT
sσT
σ

(7.23)

where ST is given by ST = seσWN ∆t . Note that this formula indicates that the MS delta is identicalto the discrete Malliavin delta. See Kohatsu-Higa and Montero (2003) about the Malliavin delta, forexample. If the pay-oﬀ function Φ(·) is smooth, proposition 2.1 in Heston and Zhou (2000) leads toS∆M
E (s, 0) =e−rTlog(PT /s) − µσT
1
e−rT
1
1
E[Φ(PT )]+O( ) =E[Φ(PT )ZT ]+O( ) = ∆BS)E (s, 0)+O(
sσT
σ
N
sσT
N
N
19

where the stochastic process Pt is a geometric Brownian motion given in (2.1) and ∆BSE (s, 0) is delta in acontinuous time model (Black and Scholes model). Even though MS delta does not approximate delta forthe binomial tree model, if the pay-oﬀ function Φ(·) is not smooth, it still is an approximation formula forthe continuous time delta. Under Assumption 7.1, corollary 4.2 in Walsh (2003) shows that the option
price in the binomial tree model converges to the options price in the Black and Scholes model. This
result shows 4S∆M
E (s, 0) =e−rTlog(PT /s) − µσT
1
1E[Φ(PT )] + O( √ ) = ∆BSE (s, 0) + O( √ ) .
sσT
σ
N
N

As previously discussed, MS delta does not approximate delta for the binomial tree model, if the pay-oﬀfunction is not smooth. On the other hand, even if the pay-oﬀ function is not smooth, MS delta still is
an approximation for continuous delta.
[Gamma]MS gamma for European options is given by (3.12). Further computation of MS gamma yields,
e−r∆t
e−r∆t
σϵ1
S
S
, ∆t)eσϵ1 (ϵ1 − µ∆t)] .
ΓM
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)] +E[∆M
E (s, 0) = − 2
E (se
s σ∆t
σs∆t

(7.24)

S
The ﬁrst and second terms in (7.24) are denoted by G1 and G2 , i.e. ΓME (s, 0) = G1 + G2 :

G1
G2

e−r∆t
= − 2
E[E(seσϵ1 , ∆t)(ϵ1 − µ∆t)]
s σ∆t
e−r∆t
S
σϵ1
=
E[∆M
, ∆t)eσϵ1 (ϵ1 − µ∆t)] .
E (se
σs∆t

The ﬁrst term is given by
G1

N
e−rT ∑e−rTE[Φ(seσWT )(ϵ1 − µ∆t)] = − 2E[Φ(seσWT )(ϵi − µ∆t)]
= − 2
s σ∆t
s σ∆tN i=1

1
e−rT
S
E[Φ(seσWT )(WT − µT )] = − ∆M(s, 0) ,
= − 2
s σT
s E
and the second term is given by
G2

=
=

e−r∆t
e−r∆t
E[ σϵ1E[E(seσ(ϵ1 +ϵ2 ) , 2∆t)(ϵ2 − µ∆t)|F∆t ]eσϵ1 (ϵ1 − µ∆t)]
σs∆t se σ∆t1
E[e−rN ∆t Φ(seσWT )(ϵ1 − µ∆t)(ϵ2 − µ∆t)] .
s2 σ 2 (∆t)2

4 The order of error term is O( √1 ), if the discontinuity for the pay-oﬀ function Φ(·) is not on a lattice point. However,N1
if all discontinuities are on lattice points, the order of the error term is O( N). See Corollary 4.2 in Walsh (2003) for details.

Also note that Chung et al. (2011) show that the rate of error terms of binomial delta for European options is O(1/N ).

20

This formula is divided into three parts, G2 = G12 − 2G22 + G32 , where G12 , G22 and G32 areG12

1
1
E[e−rT Φ(seσWT )ϵ1 ϵ2 ] = 2 2 2 E[e−rT Φ(seσWT )ϵi ϵj(̸=i) ]s2 σ 2 ∆t2s σ ∆t∑1
NE[e−rT Φ(seσWT ){(ϵi )2 − N ∆t}]222
s σ T N −1
i

=
=

1
E[e−rT Φ(seσWT )(WT2 − T )] + O(1/N )s2 σ 2 T 2
1µE[e−rT Φ(seσWT )ϵ1 µ∆t] = 2 2 E[e−rT Φ(seσWT )WT ]s2 σ 2 ∆t2
s σ T
1
µ2
−rT
σWT
2
2
E[e
Φ(se
)µ
∆t
]
=
E[e−rN ∆t Φ(seσWT )] .
s2 σ 2 ∆t2
s2 σ 2

=
G22

=

G32

=

These results are combined into the closed-form formulas for MS Gamma:SΓM
E (s, 0) =(WT − µT )2
1
e−rT
σWT
E[Φ(se){− (WT − µT ) − }] + O(1/N ) .
s2 σT
σT
σ

As discussed in the delta case, Walsh (2003) yields,SΓM
E (s, 0) =

=
=

e−rT
(log(ST /s) − µσT )2log(ST /s) − µσT
1E[Φ(ST ){−− }] + O(1/N )
2
s σT
σ3 T
σ
σ
√
e−rT
(log(PT /s) − µσT )2log(PT /s) − µσT
1
E[Φ(P
){
−
− }] + O(1/ N )
T
2
3
s σT
σ T
σ
σ
√
√
e−rT
ZT2
1
E[Φ(PT ){− ZT − }] + O(1/ N ) = ΓBSE (s, 0) + O(1/ N ) ,
s2 σT
σT
σ

even if Φ(·) is not smooth.
[Vega]MS Vega for the binomial tree model is given by (3.14). The expectation E[V1i ] is given by
E[V1i ] = e−r(N −i)∆t {−(1 +2rϵi+1 − µ∆t)}E[Φ(seσWT )] .
σ22Under the condtion i ̸= N − 1, the expectation E[V2i ] is given by
E[V2i ] =E[e−r∆te−r∆tseσ(Wi∆t +ϵi+1 ) σ∆tE[E(seσWi∆t +ϵi+1 eσϵi+2 , (i + 2)∆t; σ)(ϵi+2 − µ∆t)|F(i+1)∆t ]×seσ(Wi∆t +ϵi+1 ) ϵi+1 ]=e−r(N −i)∆t{E[Φ(seσWT )ϵi+2 ϵi+1 ] − E[Φ(seσWT )ϵi+1 ]µ∆t} .
σ∆tFormulasE[Φ(seσWT )ϵi+2 ϵi+1 ]= E[Φ(seσWT )ϵi ϵj(̸=i) ] ==E[Φ(seσWT )ϵi+1 ]

=

N
∑
1
σWT
E[Φ(se){(ϵi )2 − N ∆t}]
N2 − Ni=11E[Φ(seσWT )(WT2 − T )]N2 − NN∑
1
1
E[Φ(seσWT )ϵi ] = E[Φ(seσWT )WT ]
N
N
i=1
21

lead toE[V2i ] =e−r(N −i)∆t1{E[Φ(seσWT )(WT2 − T )] − E[Φ(seσWT )WT ]µ∆t} (i ̸= N − 1)σT
N −1

If the pay-oﬀ function Φ(·) is a smooth one, the expectation e−r(N −1)∆t E[V2N −1 ] is calculated ase−r(N −1)∆t E[V2N −1 ]√= e−rT seσW(N −1)∆t ∆t
√

√

√

√

×[pΦ′ (seσ(W(N −1)∆t + ∆t) )eσ ∆t − (1 − p)Φ′ (seσ(W(N −1)∆t − ∆t) )e−σ ∆t ]= e−rT seσW(N −1)∆t ∆t ×

∂ ′ σ(W(N −1)∆t +z) σz
Φ (se)e |z=0 + O(∆t3/2 )∂z= e−rT sσeσW(N −1)∆t ∆t
×(Φ′′ (seσW(N −1)∆t )seσW(N −1)∆t + Φ′ (seσW(N −1)∆t )) + O(∆t3/2 ) = O(1/N ) ,If the pay-oﬀ function Φ(x) is not smooth, the relatione−r(N −1)∆t E[V2N −1 ] = 0 (= O(1/N )) ,still satisﬁed, because we formally assumed V2N −1 = 0. These results are combined intoMS(s, 0; σ)
VE

=

N
−1
∑

e−rT {−(1 +i=0+

ϵi+1 − µ∆t2r)}E[Φ(seσWT )]
2
σ
2

N
−2
∑

e−rT
1
{
E[Φ(seσWT )(WT2 − T )] − E[Φ(seσWT )WT ]µ∆t}
σT
N
−
1
i=0

= e−rT E[((WT − µT )2
1
1− (WT − µT ) − )Φ(seσWT )] + O( ) .
σTσ
NAs discussed in delta case, the approximation(log(ST /s) − µσT )2log(ST /s) − µσT11−− )Φ(ST )] + O( )
σ3 T
σ
σ
N
(log(PT /s) − µσT )2log(PT /s) − µσT11= e−rT E[(−− )Φ(PT )] + O( √ )
3
σ T
σ
σ
N
2
Z
1
1
= e−rT E[( T − ZT − )Φ(PT )] + O( √ )
σT
σ
N
1
BS
= VE
(s, 0; σ) + O( √ )NMS
VE(s, 0; σ) =e−rT E[{is valid, even if the pay-oﬀ function Φ(·) is not smooth.
[Rho]
In this subsection, we derive the closed-form formula of rho for European options. We also show that MS
rho converges to rho in the continuous time model. On the other hand, the formulaE[(ϵi − µ∆tϵi − µ∆t− ∆t)E(seσWi∆t , ∆t; r)] = e−r(N −i)∆t E[(− ∆t)Φ(seσWN ∆t )]σσ22

leads the further calculation of rho. This formula is plugged into the formula (7.22) and we have
S
−rT
ρM
E[(
E (s, 0; r) = e

WT − µTS− T )Φ(ST )] = T (s∆M
E (s, 0) − E(s, 0; r)) .
σ

The last equality comes from the closed-form formula for MS delta given by (3.9)(or (7.23)). As is thediscussions in the previous cases, the formula
√
√
S
BS
−rT
ρM
Φ(PT )]) + O(1/ N ) = ρBSE (s, 0; r) = T (s∆E (s, 0) − E[eE (s, 0) + O(1/ N )must be satisﬁed under certain conditions given by Walsh (2001).

Reference[1] Bally, V., L. Caramellino and A. Zanette (2005) Pricing and Hedging American Options by Monte
Carlo Methods Using a Malliavin Calculus Approach, Monte Carlo Methods and Applications 11,
97-133.
[2] Bernis, G., E. Gobet, and A. Kohatsu-Higa (2003) Monte Carlo Evaluation of Greeks for Multidimensional Barrier and Lookback Options, Mathematical Finance 13, 99-113.
[3] Chung, S. L., W. Hung, H. H. Lee, and P. T. Shih (2011) On the Rate of Convergence of Binomial
Greeks, Journal of Future Markets 31, 562-597.
[4] Corcuera, J. M. and A. Kohatsu-Higa (2011) Statistical Inference and Malliavin Calculus, Seminaron Stochastic Analysis, Random Fields and Applications VI.
[5] Chung, S. L. and M. Shackleton (2002) The Binomial Black-Scholes Model and the Greeks, Journalof Futures Markets 22, 143-153.
[6] Cox, J. C., S. A. Ross and M. Rubinstein (1979) Option Pricing: A Simpliﬁed Approach, Journal of
Financial Economics 7, 229-263.
[7] Fournié, E., J. Laszry, J. Lebuchoux, P. Lions, and N. Touzi (1999) Applications of Malliavin Calculus
to Monte Carlo Methods in Finance, Finance and Stochastics 3, 391-412.
[8] Heston, S. and G. Zhou (2000) On the Rate of Convergence of Discrete-time Contingent Claims,
Mathematical Finance 10, 53-75.
23

[9] Hull, J. (2008) Options, Futures, and Other Derivatives (7th ed.), Prentice Hall.
[10] Kohatsu-Higa, A. and M. Montero (2003) Malliavin Calculus Applied to Finance, Physica A 320,548-570.
[11] Loisel, S. and N., Privault (2009) Sensitivity Analysis and Densisty Estimation for Finite-Time Ruin
Probabilities, Journal of Computational and Applied Mathematics230, 107-120.
[12] Leisen, D. and M. Reimer (1996) Binomial Models for Option Valuation - Examining and Improving
Convergence, Applied Mathematical Finance 3, 319-346.
[13] Leitz-Martini, M. (2000) A Discrete Clark-Ocone Formula, Maphysto Research Report 29.
[14] Longstaﬀ, F. and E. S. Schwartz (2001) Valuing American Options by Simulation: A Simple Least
Squares Approach, The Review of Financial Studies 14, 113-147.
[15] Muroi, Y. and S. Suda (2013) Discrete Malliavin Calculus and Computations of Greeks in the
Binomial Tree, European Journal of Operational Research 231, 349-361.
[16] Pelsser, A. and T. Vorst, (1994) The Binomial Model and the Greeks, Journal of Derivatives 1,
45-49.
[17] Privault, N. (2008) Stochastic Analysis of Bernoulli Processes, Probability Surveys 5, 435-483.
[18] Privault, N. (2009) Stochastic Analysis in Discrete and Continuous Settings, Springer.
[19] Privault, N. and X. Wei (2004) A Malliavin Calculus Approach to Sensitivity Analysis in Insurance,
Insurance: Mathematics and Economics35, 679-690.
[20] Privault, N. and X. Wei (2007) Integration by Parts for Point Processes and Monte Carlo Estimastion,
Journal od Applied Probability44, 806-823.
[21] Suda, S. and Y. Muroi (2013) Computation of Greeks using Binomial Trees in a Jump-Diﬀusion
Model, submitted.
[22] Tian, Y. (1993) A Modiﬁed Lattice Approach to Option Pricing, The Journal of Futures Markets
13, 563-577.
24

[23] Walsh, J. (2003) The Rate of Convergence of the Binomial Tree Scheme, Finance and Stochastics 7,
337-361.

25

Figure 2: Delta for American put options (MS delta and Finite Diﬀerence delta, K=100)

Figure 3: Gamma for American put options (MS gamma and Finite Diﬀerence gamma, K=100)

26

Figure 4: Vega for American put options (MS vega and Finite Diﬀerence vega, K=100)

Figure 5: Vega for American put options (MS vega and Finite Diﬀerence vega, K=105)27

Figure 6: Rho for American put options (MS rho and Finite Diﬀerence rho, K=100)

Figure 7: Delta for American put options as function of price of underlying assets. The strike price ofoptions is ﬁxed at K=100.

28

Figure 8: Gamma for American put options as function of price of the underlying assets. The strike priceof options is ﬁxed at K=100.

Figure 9: Vega for American put options as function of price of underlying assets. The strike price ofoptions is ﬁxed at K=100.

29

Figure 10: Rho for American put options as function of price of underlying assets. The strike price ofoptions is ﬁxed at K=100.

Figure 11: Delta (MS delta and delta computed by various kinds of binomial trees.)
MS: MS delta. CS: delta computed by binomial tree introduced by Chung and Shackelton (2002). Tian:delta computed by binomial tree introduced by Tian (1993).
introduced by Leisen and Reimer (1996).

30

LR: delta computed by binomial tree