Computation of Greeks Using Binomial Tree

Yoshifumi Muroi¹, Shintaro Suda²

¹Graduate School of Economics and Management, Tohoku University, Sendai, Japan.
²Mitsubishi UFJ Trust Investment Technology Institute Co., Ltd. (MTEC), Tokyo, Japan.

DOI: 10.4236/jmf.2017.73031                

Abstract

This paper proposes a new efficient 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 algorithm of Greeks for American options using the binomial tree is also given in this article. There are three advantageous points to use binomial tree approach for the computation of Greeks. First, mathematics is much simpler than using the continuous time Malliavin calculus approach. Second, we can construct a simple algorithm to obtain the Greeks for American options. Third, this algorithm is very efficient because 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 difficult. 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 to changes in the price of underlying assets or parameters. They are defined by derivatives of the option price function with respect to parameters such as the price of underlying assets, volatility level, and spot interest rate. The computation of these sensitivities is very important for risk management, for an example on this and similar topics, refer to Hull [1] . The recent development of the Malliavin calculus approach in financial mathematics enables us to compute Greeks for a various kinds of contingent claims, for a previous research on this topic, see Fournié et al. [2] and Kohatsu-Higa and Montero [3] among others. See also Bernis et al. [4] 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 the computation of Greeks for American options have been reported. See Gobet [5] and Bally et al. [6] , for example. One can obtain Greeks by simply taking the differential of the option pricing function if the closed-form pricing formula is known. However, one generally cannot compute Greeks for American options using this direct method because the explicit formula for the price of American options is not known. The most direct approximation method to derive Greeks for American options is the finite difference method. However, it is well known that the finite difference method sometimes leads to unstable results. In order to overcome this short- coming, Malliavin calculus approaches combined with Monte Carlo methods are widely used. Despite recent progress of the Monte Carlo simulation approach to the pricing of American options, such as Longstaff and Schwartz [7] , it is not easy to derive the price and Greeks for American options using this approach. Although Gobet [5] , Bally et al. [6] considered the Monte Carlo simulation an approach to compute Greeks for American options, the approach is mathema- tically difficult to understand. Therefore, the finite difference approach is still widely used in practical purposes to compute vega and rho. Therefore, we introduce a new approach in this paper.

Muroi and Suda [8] [9] proposed new methods for the computation of Greeks for European options using the binomial tree methods of Cox et al. [10] and the discrete Malliavin calculus introduced by Leitz-Martini [11] and Privault [12] [13] . See also Privault and Schoutens [14] . In that article, we took derivatives of the expectation form-formula for European options directly and computed it further using relationships between the discrete Malliavin calculus and the discrete Skorohod integrals. Although this is an elegant way to derive Greeks for European options, this method still requires a closed-form formula for the option price to derive option Greeks. Muroi and Suda [8] [9] took derivatives of the pricing formula for European options, however, in this article we take derivative at each node on the binomial tree to derive Greeks for American options. In other words, we employed a step-by-step approach. In this article, we also show Greeks for European options are asymptotically equivalent to discrete version of the Malliavin Greeks and the binomial Greeks are converging to Greeks in the continuous time model. This model is also extended to the jump- diffusion model in Suda and Muroi [15] . Chung et al. [16] also proved the binomial delta for plain vanilla European options is converging to delta in the continuous time model.

The remainder of this paper is organized as follows. In Section 2, we give brief explanations on the binomial tree methods of Cox et al. [10] . 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, we show the derived Greeks converge to the Malliavin Greeks in Black and Scholes model.

2. Binomial Tree

In this section, we briefly explain the binomial tree model of Cox et al. [10] , which has been explained in many textbooks such as Hull [1] . 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 the binomial tree are explained in detail in the next section.

Introduce independently and identically distributed random variables ${\left\{{ϵ}_i\right\}}_{i=\mathrm{1,}\cdots \mathrm{,}N}$ on the probability space $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}Q\right)$ , where ${ϵ}_i$ is a random variable with probability $Q\left[{ϵ}_i=\sqrt{\Delta t}\right]=p$ , $Q\left[{ϵ}_i=-\sqrt{\Delta t}\right]=1-p$ . The time step $\Delta t$ is fixed as $T=N\Delta t$ . We introduce a random walk process, ${\left\{W_{i\Delta t}\right\}}_{i=1,\cdots ,N}$ , given by

$W_{i\Delta t}={ϵ}_1+\cdots +{ϵ}_i.$

One can regard the stochastic process ${\left\{W_{i\Delta t}\right\}}_{i=1,\cdots ,N}$ with $p=1/2$ as an ap- proximation of the standard Brownian motion. On the other hand, the stocha- stic process ${\left\{W_{i\Delta t}\right\}}_{i=1,\cdots ,N}$ is generally regarded as an approximation of the

Brownian motion with drift $\frac{2p-1}{\sqrt{\Delta t}}$ in the general case ( $p\ne \frac{1}{2}$ ).

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

$\text{d}{P}_t=P_t\left(r\text{d}t+\sigma \text{d}{Z}_t\right),P_0=s,$ (1)

where ${\left\{Z_t\right\}}_{t\ge 0}$ is a standard Brownian motion under the risk-neutral measure Q. A binomial tree is constructed in the following manner. We consider a model with N periods and assume that the maturity date of the options is fixed as $T=N\Delta t$ . If the price of underlying assets at time $i\Delta t\left(i=0,1,\cdots ,N-1\right)$ is given by $S_{i\Delta t}$ , the price moves to $u{S}_{i\Delta t}$ or $d{S}_{i\Delta t}\left(d<1<u\right)$ in the next time period $\left(i+1\right)\Delta t$ . The probability of the price moving upward to ( $u{S}_{i\Delta t}$ ) is p and downward to ( $d{S}_{i\Delta t}$ ) is $1-p$ . In order to construct the binomial tree so that the expectation and variance are consistent with the geometric Brownian motion (1), we fix the parameters u, d, p as

$u=e^{\sigma \sqrt{\Delta t}},d=e^{-\sigma \sqrt{\Delta t}},p=\frac{e^{r\Delta t}-d}{u-d}$

where we assume an additional relationship, $u=1/d$ . Consider European options with a pay-off function $\Phi (\cdot )$ and a maturity date T. The price of options at time $i\Delta t$ is denoted by $\mathcal{E}\left(x\mathrm{,}i\Delta t\right)$ if the price of underlying assets is x at time $i\Delta t$ . The price is given by the backward induction algorithm

$\mathcal{E}\left(x,i\Delta t\right)=e^{-r\Delta t}\left(p\mathcal{E}\left(xu,\left(i+1\right)\Delta t\right)+\left(1-p\right)\mathcal{E}\left(xd,\left(i+1\right)\Delta t\right)\right)$

$\mathcal{E}\left(x,N\Delta t\right)=\Phi \left(x\right).$ (2)

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

3. Computation of Greeks for European Options

New computational methods of Greeks for European options are presented in this section. Although computational methods of delta and gamma using the binomial tree have already been proposed by Pelsser and Vorst [17] , computa- tional methods of other Greeks such as vega and rho using the binomial tree have not been deeply studied, except in the finite difference approach.

Because the option price is given by the weighted sum of the pay-off function, we can compute Greeks if we assume that the pay-off function is smooth. However, the pay-off function for, say, European call options is not smooth. Therefore, the smoothness of the pay-off function is too strong to be assumed for practical purposes. In this section, however, we compute Greeks under the following assumption. We will show that the obtained formulas in this section converge to Greeks in the continuous time model under the milder conditions. It is discussed in Section 7.2.

Assumption 1. We assume that the pay-off function $\Phi (\cdot )$ is a smooth function.

3.1. Computation of Delta

In this subsection, we calculate delta, which is used to measure the sensitivity of the option price with respect to changes in the price of underlying assets. Delta is given by the first derivative of the option value function with respect to the price of underlying assets. Delta for European options is computed as

$\begin{array}{c}{\Delta }_E\left(s,0\right)=e^{-r\Delta t}\left\{p\frac{\partial }{\partial s}\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)+\left(1-p\right)\frac{\partial }{\partial s}\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\right\}\\ =e^{-r\Delta t}\left\{p\frac{\text{d}}{\text{d}x}\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)e^{\sigma \sqrt{\Delta t}}+\left(1-p\right)\frac{\text{d}}{\text{d}x}\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)e^{-\sigma \sqrt{\Delta t}}\right\}.\end{array}$ (3)

Now, we have two approaches to compute ${\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z}\mathrm{,}\Delta t\right)|}_{z=\pm \sqrt{\Delta t}}$ . The first one is given by

${\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=\pm \sqrt{\Delta t}}=\frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)}{2\sqrt{\Delta t}}+O\left(\sqrt{\Delta t}\right),$ (4)

and the second one is given by

${\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=\pm \sqrt{\Delta t}}={\frac{\text{d}}{\text{d}x}\mathcal{E}\left(x,\Delta t\right)|}_{x=s{e}^{\pm \sigma \sqrt{\Delta t}}}\times {\sigma s{e}^{\sigma z}|}_{z=\pm \sqrt{\Delta t}}.$

These two formulas yield the approximation formula for the derivative ${\frac{\text{d}}{\text{d}x}\mathcal{E}\left(x\mathrm{,}\Delta t\right)|}_{x\mathrm{<mo>\; =\; </mo>}s{e}^{\pm \sigma \sqrt{\Delta t}}}$ :

${\frac{\text{d}}{\text{d}x}\mathcal{E}\left(x,\Delta t\right)|}_{x=s{e}^{\pm \sigma \sqrt{\Delta t}}}=\frac{1}{\sigma s{e}^{\pm \sigma \sqrt{\Delta t}}}\frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)}{2\sqrt{\Delta t}}+O\left(\sqrt{\Delta t}\right).$ (5)

Substituting this formula in (3) leads to the following relation

#Math_42# (6)

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

$p=\frac{1}{2}+\frac{1}{2}\mu \sqrt{\Delta t}+O\left(\Delta t^{3/2}\right),$ (7)

where $\mu$ is a constant given by $\mu =\frac{1}{\sigma }\left(r-\frac{{\sigma }^2}{2}\right).$ Let us compute the expec- tation $E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1}\mathrm{,}\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]$ . The above approximation formula leads to

$\begin{array}{l}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ =p\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)\left(\sqrt{\Delta t}-\mu \Delta t\right)+\left(1-p\right)\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\left(-\sqrt{\Delta t}-\mu \Delta t\right)\\ =\left[\frac{\sqrt{\Delta t}}{2}\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)+\frac{-\sqrt{\Delta t}}{2}\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\right]+O\left(\Delta t^2\right).\end{array}$ (8)

Note that we have used the identities

$\frac{1}{2}\left(\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)+\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\right)=\mathcal{E}\left(s,\Delta t\right)+O\left(\sqrt{\Delta t}\right)$

$\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)={2\sqrt{\Delta t}\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=0}+O\left(\Delta t\right)=O\left(\sqrt{\Delta t}\right)$

to derive the last equality. This formula allows further computation of (6). It is given by

$\begin{array}{c}{\Delta }_E\left(s,0\right)=\frac{e^{-r\Delta t}}{s\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]+O\left(\sqrt{\Delta t}\right)\\ \approx \frac{e^{-r\Delta t}}{s\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\equiv {\Delta }_E^{MS}\left(s,0\right)\end{array}$ (9)

We call this MS delta, which is the one-step version of discrete Malliavin delta given in Muroi and Suda [8] . Actually, we can show the stronger result,

${\Delta }_E^{MS}\left(s,0\right)={\Delta }_E\left(s,0\right)+O\left(\Delta t\right),$ (10)

if we assume smoothness to the pay-off function. This is organized as the following theorem, because we use the Formula (10) to show the accuracy of Vega.

Theorem 1. We assume that the pay-off function $\Phi (\cdot )$ as a smooth function. We can estimate the accuracy of MS delta as

${\Delta }_E^{MS}\left(s\mathrm{,0}\right)={\Delta }_E\left(s\mathrm{,0}\right)+O\left(\Delta t\right)\mathrm{.}$

[Proof] First we compute a higher order term for the error term for MS delta. Taylor expansion,

$\begin{array}{c}\mathcal{E}\left(s{e}^{\mp \sigma \sqrt{\Delta t}},\Delta t\right)=\mathcal{E}\left(s{e}^{\pm \sigma \sqrt{\Delta t}},\Delta t\right)+{\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)\left(\mp 2\sqrt{\Delta t}\right)|}_{z=\pm \sqrt{\Delta t}}\\ +{\frac{1}{2}\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right){\left(\mp 2\sqrt{\Delta t}\right)}^2|}_{z=\pm \sqrt{\Delta t}}+O\left(\Delta t^{3/2}\right),\end{array}$

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

$\begin{array}{c}{\frac{\partial }{\partial z}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=\pm \sqrt{\Delta t}}=\frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)}{2\sqrt{\Delta t}}\\ \pm {\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)\sqrt{\Delta t}|}_{z=\pm \sqrt{\Delta t}}+O\left(\Delta t\right).\end{array}$

This leads to the higher order expansion formula for the delta given by (6):

$\begin{array}{l}{\Delta }_E\left(s,0\right)=\frac{e^{-r\Delta t}}{s\sigma \Delta t}\frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)\sqrt{\Delta t}+\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\left(-\sqrt{\Delta t}\right)}{2}\\ +\frac{e^{-r\Delta t}}{s\sigma }\left({p\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=\sqrt{\Delta t}}-\left(1-p\right){\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=-\sqrt{\Delta t}}\right)\sqrt{\Delta t}+O\left(\Delta t\right).\end{array}$ (11)

The second term in (11) is calculated as

$\begin{array}{l}p{\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=\sqrt{\Delta t}}-{\left(1-p\right)\frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=-\sqrt{\Delta t}}\\ =\left({\mu \frac{{\partial }^2}{\partial z^2}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=0}+{\frac{{\partial }^3}{\partial z^3}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)|}_{z=0}\right)\sqrt{\Delta t}+O\left(\Delta t\right)\\ =O\left(\sqrt{\Delta t}\right),\end{array}$

where we have used the approximation (7). This result enables us to compute delta given by (11). It is computed as

${\Delta }_E\left(s,0\right)=\frac{e^{-r\Delta t}}{s\sigma \Delta t}\frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)\sqrt{\Delta t}+\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)\left(-\sqrt{\Delta t}\right)}{2}+O\left(\Delta t\right).$

Applying the Formula (8), the higher order expansion for delta is finally given by

${\Delta }_E\left(s\mathrm{,0}\right)={\Delta }_E^{MS}\left(s\mathrm{,0}\right)+O\left(\Delta t\right)\mathrm{.}$ (12)

This shows the theorem.

Discrete Malliavin Greeks for European options using an N-steps binomial tree were obtained by Muroi and Suda [8] ; they used the discrete Malliavin derivatives introduced by Leitz-Martini [11] . See also Privault [12] [13] . Discrete Malliavin delta ${\Delta }_E^D$ is given by

${\Delta }_E^D=\frac{e^{-rN\Delta t}}{s\sigma N\Delta t}E\left[\Phi \left(s{e}^{\sigma W_{N\Delta t}}\right)\left(W_{N\Delta t}-\mu N\Delta t\right)\right].$

As discussed in Appendix, MS delta and discrete Malliavin delta is actually equivalent. Because Muroi and Suda [8] used the discrete version of the Mal- liavin calculus approach, one cannot use their method to derive Greeks for Ame- rican 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.

If we exploit another approximation formula for the derivative ${\frac{\text{d}}{\text{d}x}\mathcal{E}\left(x\mathrm{,}\Delta t\right)|}_{x=s{e}^{\pm \sigma \sqrt{\Delta t}}}$ , we can obtain another approximation formula for delta. If we use a more direct formula for the derivative,

$\begin{array}{c}\frac{\text{d}}{\text{d}x}\mathcal{E}\left(s{e}^{\pm \sigma \sqrt{\Delta t}},\Delta t\right)\approx \frac{\mathcal{E}\left(s{e}^{\mp \sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{\pm \sigma \sqrt{\Delta t}},\Delta t\right)}{s{e}^{\mp \sigma \sqrt{\Delta t}}-s{e}^{\pm \sigma \sqrt{\Delta t}}}\\ \approx \frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)}{2s\sigma \sqrt{\Delta t}},\end{array}$ (13)

we get another approximation formula for delta for European options in the one-period time model. Substituting (13) in (3) yields ${\Delta }_E^{Hull}$ , and it is given as

${\Delta }_E^{Hull}\left(s\mathrm{,0}\right)\approx \frac{\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)-\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)}{2s\sigma \sqrt{\Delta t}}\mathrm{.}$

This is the delta introduced by Hull [1] . This fact reveals that the two different approximation formulas for $\frac{\text{d}}{\text{d}x}\mathcal{E}\left(s{e}^{\pm \sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)$ given by (5) and (13) yield two

different approximation formulas for delta¹. Our computational results indicate that MS delta converges a little bit faster than delta introduced by Hull [1] . Moreover, as shown in the Appendix, the formula, ${\Delta }_E^{MS}\left(s\mathrm{,0}\right)={\Delta }_E^D\left(s\mathrm{,0}\right)$ , is actually satisfied. This means that MS delta is more natural representation of delta than Hull’s delta.

Lastly, we can show the convergence of MS delta to delta in Black and Scholes model, even if we do not assume the smoothness to the pay-off function. This is shown in Appendix.

3.2. Computation of Gamma

¹Actually, MS delta is deeply related to Hull’s delta. We have a relation, ${\Delta }_E^{MS}\left(s\mathrm{,0}\right)=e^{-r\Delta t}{\Delta }_E^{Hull}\left(s\mathrm{,0}\right)+O\left(\sqrt{\Delta t}\right)$ , because we have Formulas (8) and (9).

In this subsection, we calculate gamma, which is used to measure the sensitivity of delta with respect to changes in the price of underlying assets. Gamma is given by the second derivative of the option value function with respect to the price of underlying assets. The pricing algorithm for European options (2) yields delta for European options:

$\begin{array}{c}{\Delta }_E\left(s\mathrm{,0}\right)=e^{-r\Delta t}\left\{p\frac{\partial }{\partial x}\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)e^{\sigma \sqrt{\Delta t}}+\left(1-p\right)\frac{\partial }{\partial x}\mathcal{E}\left(s{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)e^{-\sigma \sqrt{\Delta t}}\right\}\\ =e^{-r\Delta t}\left\{p{\Delta }_E\left(s{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)e^{\sigma \sqrt{\Delta t}}+\left(1-p\right){\Delta }_E\left(s{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)e^{-\sigma \sqrt{\Delta t}}\right\}\mathrm{.}\end{array}$

Applying chain rule to $\mathcal{E}\left(s{e}^{\sigma z}\mathrm{,}\Delta t\right)$ leads to

${\Delta }_E\left(s{e}^{\pm \sigma \sqrt{\Delta t}}\mathrm{,}\Delta t\right)={\frac{\frac{\partial \mathcal{E}}{\partial z}\left(s{e}^{\sigma z}\mathrm{,}\Delta t\right)}{\sigma s{e}^{\sigma z}}|}_{z=\pm \sqrt{\Delta t}}\mathrm{.}$

Delta for European options is given by

${\Delta }_E\left(s,0\right)=\frac{e^{-r\Delta t}}{\sigma s}\left\{p{\left[\frac{\partial \mathcal{E}}{\partial z}\left(s{e}^{\sigma z},\Delta t\right)\right]}_{z=\sqrt{\Delta t}}+\left(1-p\right){\left[\frac{\partial \mathcal{E}}{\partial z}\left(s{e}^{\sigma z},\Delta t\right)\right]}_{z=-\sqrt{\Delta t}}\right\}.$

Taking derivative with respect to the price of underlying assets yields gamma for European options as

$\begin{array}{c}{\Gamma }_E\left(s,0\right)=-\frac{1}{s}{\Delta }_E\left(s,0\right)+\frac{e^{-r\Delta t}}{\sigma s}\{{p\frac{\partial }{\partial z}\left({\Delta }_E\left(s{e}^{\sigma z},\Delta t\right)e^{\sigma z}\right)|}_{z=\sqrt{\Delta t}}\\ +{\left(1-p\right)\frac{\partial }{\partial z}\left({\Delta }_E\left(s{e}^{\sigma z},\Delta t\right)e^{\sigma z}\right)|}_{z=-\sqrt{\Delta t}}\}.\end{array}$

Using the approximation formula

$\begin{array}{l}{\frac{\partial }{\partial z}{\Delta }_E\left(s{e}^{\sigma z},\Delta t\right)e^{\sigma z}|}_{z=\pm \sqrt{\Delta t}}\\ =\frac{{\Delta }_E\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right)e^{\sigma \sqrt{\Delta t}}-{\Delta }_E\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right)e^{-\sigma \sqrt{\Delta t}}}{2\sqrt{\Delta t}}+O\left(\sqrt{\Delta t}\right)\end{array}$

allows further calculation of gamma as

$\begin{array}{c}{\Gamma }_E\left(s,0\right)=-\frac{e^{-r\Delta t}}{s^2\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({\epsilon }_1-\mu \Delta t\right)\right]\\ +\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[{\Delta }_E\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right]+O\left(\sqrt{\Delta t}\right)\\ =-\frac{e^{-r\Delta t}}{s^2\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ +\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[{\Delta }_E^{MS}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right]+O\left(\sqrt{\Delta t}\right)\\ \approx -\frac{e^{-r\Delta t}}{s^2\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ +\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[{\Delta }_E^{MS}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right]\\ \equiv {\Gamma }_E^{MS}\left(s,0\right)\end{array}$ (14)

We call this formula MS gamma. This formula is valid only if the pay-off function $\Phi (\cdot )$ is a smooth function. In this case, the order of the error term is $O\left(\sqrt{\Delta t}\right)$ , i.e. ${\Gamma }_E\left(s\mathrm{,0}\right)={\Gamma }_E^{MS}\left(s\mathrm{,0}\right)+O\left(\sqrt{\Delta t}\right)$ . Second equality in (14) is shown by using the Formula (10). We will also prove that MS gamma is asymptotically equivalent to the Gamma in the Black and Scholes model in the Appendix.

3.3. Computation of Vega

The computational method of vega for European options is presented in this subsection. Vega is the sensitivity of the option pricing formula with respect to changes in volatility level, $\sigma$ . It appears that the computational methods of vega and rho using the binomial tree have not yet been considered seriously, except for the finite difference approach. See Hull [1] for computation of vega using the finite difference approach with the binomial tree. Let us assume that the price for underlying assets at time $i\Delta t$ is given by $s_{i\Delta t}$ . The price and vega for European options at time $i\Delta t$ are denoted by $\mathcal{E}\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)$ and ${\mathcal{V}}_E\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)$ , respectively. Vega for European options is given by taking de- rivatives to the pricing Formula (2):

$\begin{array}{c}{\mathcal{V}}_E\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)=e^{-r\Delta t}\frac{\partial }{\partial \sigma }\{p\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\\ +\left(1-p\right)\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\}\\ ={\mathcal{V}}_1^i+{\mathcal{V}}_2^i+{\mathcal{V}}_3^i\mathrm{,}\end{array}$

where ${\mathcal{V}}_1^i\mathrm{,}{\mathcal{V}}_2^i\mathrm{,}{\mathcal{V}}_3^i$ are given by

$\begin{array}{l}{\mathcal{V}}_1^i=e^{-r\Delta t}\{\frac{\partial p}{\partial \sigma }\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\\ +\frac{\partial \left(1-p\right)}{\partial \sigma }\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\}\end{array}$

$\begin{array}{c}{\mathcal{V}}_2^i=e^{-r\Delta t}\{p\frac{\partial }{\partial x}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\sqrt{\Delta t}\\ +\left(1-p\right)\frac{\partial }{\partial x}\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\left(-\sqrt{\Delta t}\right)\}\\ =e^{-r\Delta t}E\left[{\Delta }_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}{ϵ}_{i+1}\right]\end{array}$

${\mathcal{V}}_3^i=e^{-r\Delta t}E\left[{\mathcal{V}}_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\right]\mathrm{.}$

The derivative of p with respect to $\sigma$ is given by

$\frac{\partial }{\partial \sigma }p=\frac{2-e^{r\Delta t}\left(e^{\sigma \sqrt{\Delta t}}+e^{-\sigma \sqrt{\Delta t}}\right)}{{\left(e^{\sigma \sqrt{\Delta t}}-e^{-\sigma \sqrt{\Delta t}}\right)}^2}\sqrt{\Delta t}=-\frac{1}{4}\left(1+\frac{2r}{{\sigma }^2}\right)\sqrt{\Delta t}+O\left(\Delta t^{3/2}\right).$

This result and the approximation Formula (8) lead to the approximation formula for ${\mathcal{V}}_1^i$ :

$\begin{array}{c}{\mathcal{V}}_1^i=-\frac{1}{2}\left(1+\frac{2r}{{\sigma }^2}\right)e^{-r\Delta t}\\ \times \frac{\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\sqrt{\Delta t}+\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\left(-\sqrt{\Delta t}\right)}{2}+O\left(\Delta t^{3/2}\right)\\ =E\left[e^{-r\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\right]+O\left(\Delta t^{3/2}\right)\mathrm{.}\end{array}$

For further computation of ${\mathcal{V}}_2^i$ , one needs delta at time $\left(i+1\right)\Delta t$ . Moreover, it is necessary to compute vega for European options at time $\left(i+1\right)\Delta t$ to evaluate ${\mathcal{V}}_3^i$ . This implies that one has to use the backward induction algorithm to compute vega for European options. Vega at the maturity date is given by ${\mathcal{V}}_E\left(s_{N\Delta t}\mathrm{,}N\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)=0$ . These results are combined to form the computational formulas for vega:

$\begin{array}{c}{\mathcal{V}}_E\left(s_{i\Delta t},i\Delta t;\sigma \right)=E\left[e^{-r\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t\right)\right]\\ +E\left[e^{-r\Delta t}{\Delta }_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}{ϵ}_{i+1}\right]\\ +E\left[e^{-r\Delta t}{\mathcal{V}}_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)\right]+O\left(\Delta t^{3/2}\right).\end{array}$ (15)

Delta and vega at time $\left(i+1\right)\Delta t$ are substituted by MS delta and vega to get MS vega at time $i\Delta t$ , respectively. In other words, MS vega at time $i\Delta t$ is defined recursively by the backward algorithm,

$\begin{array}{c}{\mathcal{V}}_E^{MS}\left(s_{i\Delta t},i\Delta t;\sigma \right)\equiv E\left[e^{-r\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t\right)\right]\\ +E\left[e^{-r\Delta t}{\Delta }_E^{MS}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}{ϵ}_{i+1}\right]\\ +E\left[e^{-r\Delta t}{\mathcal{V}}_E^{MS}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)\right].\end{array}$ (16)

We have the following theorem on MS Vega.

Theorem 2 Let us assume that the pay-off function $\Phi (\cdot )$ is a smooth function. The accuracy of the MS vega is given by

${\mathcal{V}}_E\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma \right)={\mathcal{V}}_E^{MS}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma \right)+O\left(\sqrt{\Delta t}\right)\mathrm{.}$

We need a little bit more effort show this fact and it is shown in Appendix. We also show the convergence of MS vega to the vega in the Black and Scholes model even if we do not assume the smoothness of the pay-off function in the Appendix. It should be noted that in order to compute Vega, we formally assume ${\mathcal{V}}_2^{N-1}=0$ if the pay-off function $\Phi (\cdot )$ is not smooth one. As an alter- native approach, the finite difference approach is used to obtain vega for prac- tical purposes. In many cases, it works well: however, in some it does 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 defined by the first derivative of the option value function with respect to the spot interest rate, r. The price and rho for European options are denoted by $\mathcal{E}\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}r\right)$ and ${\rho }_E\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}r\right)$ , respec- tively, if the price of underlying assets at time $i\Delta t$ is given by $s_{i\Delta t}$ . Simple cal- culation yields

${\rho }_E\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}r\right)={\rho }_1^i+{\rho }_2^i+{\rho }_3^i\mathrm{,}$

where ${\rho }_1^i\mathrm{,}{\rho }_2^i\mathrm{,}{\rho }_3^i$ are given by

${\rho }_1^i=-\Delta t{e}^{-r\Delta t}E\left[\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)\right]$

${\rho }_2^i=e^{-r\Delta t}\left\{\frac{\partial p}{\partial r}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;r\right)+\frac{\partial \left(1-p\right)}{\partial r}\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;r\right)\right\}$

${\rho }_3^i=e^{-r\Delta t}E\left[{\rho }_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)\right].$

The derivative $\frac{\partial p}{\partial r}$ in ${\rho }_2^i$ is given by

$\frac{\partial }{\partial r}p=\frac{e^{r\Delta t}}{e^{\sigma \sqrt{\Delta t}}-e^{-\sigma \sqrt{\Delta t}}}\Delta t=\frac{\sqrt{\Delta t}}{2\sigma }+O\left(\Delta t^{3/2}\right).$

This relation and the Formula (8) leads to the further calculation,

$\begin{array}{c}{\rho }_2^i=\frac{e^{-r\Delta t}}{\sigma }\{\frac{\sqrt{\Delta t}}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;r\right)\\ +\frac{-\sqrt{\Delta t}}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;r\right)\}+O\left(\Delta t^{3/2}\right)\\ =\frac{e^{-r\Delta t}}{\sigma }E\left[\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)\left({ϵ}_{i+1}-\mu \Delta t\right)\right]+O\left(\Delta t^{3/2}\right).\end{array}$

In order to compute ${\rho }_3^i$ , one needs to compute rho at time $\left(i+1\right)\Delta t$ , i.e., ${\rho }_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}r\right)$ . This implies that one has to use the backward induc- tion approach to evaluate rho. These results are combined to form the computa- tional formulas of rho for European options:

$\begin{array}{c}{\rho }_E\left(s_{i\Delta t},i\Delta t;r\right)=e^{-r\Delta t}E[\left(\frac{{ϵ}_{i+1}}{\sigma }-\frac{\mu \Delta t}{\sigma }-\Delta t\right)\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)\\ \begin{array}{c}\\ \end{array}+{\rho }_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)]+O\left(\Delta t^{3/2}\right).\end{array}$ (17)

In order to compute MS rho at time $i\Delta t$ , rho at time $\left(i+1\right)\Delta t$ is substituted by MS rho at time $\left(i+1\right)\Delta t$ . In other words, MS rho is defined recursively by the backward algorithm,

$\begin{array}{c}{\rho }_E^{MS}\left(s_{i\Delta t},i\Delta t;r\right)\equiv e^{-r\Delta t}E[\left(\frac{{ϵ}_{i+1}}{\sigma }-\frac{\mu \Delta t}{\sigma }-\Delta t\right)\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)\\ \begin{array}{c}\\ \end{array}+{\rho }_E^{MS}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;r\right)].\end{array}$

Because rho at the maturity date is equal to 0, rho given by the Formula (17) is further calculated as

$\begin{array}{c}{\rho }_E\left(s,0;r\right)=e^{-r\Delta t}E\left[\left(\frac{{ϵ}_1-\mu \Delta t}{\sigma }-\Delta t\right)\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t;r\right)+{\rho }_E\left(s{e}^{\sigma {ϵ}_1},\Delta t;r\right)+O\left({\left(\Delta t\right)}^{3/2}\right)\right]\\ =\cdots \\ ={\displaystyle \underset{i=1}{\overset{N}{\sum }}}{e}^{-ri\Delta t}E\left[\left(\frac{{ϵ}_i-\mu \Delta t}{\sigma }-\Delta t\right)\mathcal{E}\left(s{e}^{\sigma W_{i\Delta t}},\Delta t;r\right)\right]+O\left(\sqrt{\Delta t}\right)\\ ={\rho }_E^{MS}\left(s,0;r\right)+O\left(\sqrt{\Delta t}\right).\end{array}$ (18)

This shows the final result,

${\rho }_E\left(s\mathrm{,0\; <mo>\; ;\; </mo>}r\right)={\rho }_E^{MS}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}r\right)+O\left(\sqrt{\Delta t}\right)\mathrm{.}$

MS rho is also asymptotically equivalent to the rho in the Black and Scholes model.

4. Computation of Greeks for American Options

We suggest a new algorithm for computation of Greeks for American options. An American option is a contingent claim that its holder can exercise at any time before its maturity date. Consider an American option with a pay-off function $\Phi (\cdot )$ and a maturity date $T=N\Delta t$ . If the price of underlying assets at time $i\Delta t$ is given by x, the price of these options at time $i\Delta t$ is denoted by $\mathcal{A}\left(x\mathrm{,}i\Delta t\right)$ . The price of American options is given by the backward induction algorithm:

$\mathcal{A}\left(x\mathrm{,}i\Delta t\right)=max\left\{e^{-r\Delta t}\left(p\mathcal{A}\left(xu\mathrm{,}\left(i+1\right)\Delta t\right)+\left(1-p\right)\mathcal{A}\left(xd\mathrm{,}\left(i+1\right)\Delta t\right)\right)\mathrm{,}\Phi \left(x\right)\right\}\mathrm{,}$

where the price of options at the maturity date is given by $\mathcal{A}\left(x\mathrm{,}N\Delta t\right)=\Phi \left(x\right)$ . Introduce new sets,

${\mathcal{S}}_i=\left\{x\mathrm{<mo>\; |\; </mo>}x\in \left(\mathrm{0,}\infty \right)\mathrm{,}\mathcal{A}\left(x\mathrm{,}i\Delta t\right)=\Phi \left(x\right)\right\}$ and ${\mathcal{C}}_i=\left(0,\infty \right)\backslash {\mathcal{S}}_i\left(i=0,1,\cdots ,N\right)$

The continuous region and stopping region for American options are defined by $\mathcal{C}={\displaystyle {\cup }_{i=0}^N{\mathcal{C}}_i}$ and $\mathcal{S}={\displaystyle {\cup }_{i=0}^N{\mathcal{S}}_i}$ , respectively. We will present an intuitive way for the computation of Greeks for American options. If the price of underlying assets at the initial time satisfies $s\in {\mathcal{S}}_0$ , Greeks are simply the derivative of the pay-off function. Sensitivity for American options at the initial time respect to a parameter $\theta$ , is computed as

${\partial }_{\theta }\mathcal{A}\left(s\mathrm{,0}\right)={\partial }_{\theta }\Phi \left(s\right)$

if the function $\Phi (\cdot )$ is differentiable at $s\in {\mathcal{S}}_0$ . (Second derivatives such as gamma is given by ${\partial }_{ss}\mathcal{A}\left(s\mathrm{,0}\right)={\partial }_{ss}\Phi \left(s\right)$ .) If the price of underlying assets at the initial time satisfies $s\in {\mathcal{C}}_0$ , Greeks are also given by the derivative of the price function of American options, and it is given by²

${\partial }_{\theta }\mathcal{A}\left(s\mathrm{,0}\right)=\frac{\partial }{\partial \theta }\left[e^{-r\Delta t}\left\{p\mathcal{A}\left(su\mathrm{,}\Delta t\right)+\left(1-p\right)\mathcal{A}\left(sd\mathrm{,}\Delta t\right)\right\}\right]\mathrm{.}$ (19)

Notice that the Formula (19) has a same functional form to the European options Greeks. Delta, for example, is given by

${\Delta }_A^{MS}\left(s\mathrm{,0}\right)=\{\begin{array}{ll}{\partial }_s\Phi \left(s\right)\hfill & \text{if}s\in {\mathcal{S}}_0\hfill \\ \frac{e^{-r\Delta t}}{s\sigma \Delta t}E\left[\mathcal{A}\left(s{e}^{\sigma {ϵ}_1}\mathrm{,}\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\hfill & \text{if}s\in {\mathcal{C}}_0\mathrm{.}\hfill \end{array}$

²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 with the pay-off $\Phi \left(x\right)={\left(K-x\right)}^{+}$ if the delta is fixed at ${\Delta }_A\left(s,0\right)=-1$ on the stopping region and the boundary. This is because we have a smooth-fit condition in the continuous model, and our model is an approximation of the Black and Scholes model.

One can compute other Greeks using same method. Numerical demonstra- tions are shown in Section 5. The extended binomial tree of Pelsser and Vorst [17] is one of the most suitable alternative methods to derive delta and gamma using binomial trees. One can efficiently and accurately derive Greeks efficiently and accurately as discussed in Section 5. However, one cannot apply this method to derive vega and rho.

5. Numerical Results

In this section, we demonstrate the numerical results for the new computational methods of Greeks that were introduced in previous sections. In order to check the effectiveness of our approach, Greeks for American put options are com- puted by the newly proposed approach and the finite difference approach. We also demonstrate the extended binomial tree approach of Pelsser and Vorst [17] to compute delta and gamma. It is well known that the extended binomial tree approach of Pelsser and Vorst [17] yields very accurate and fast algorithms to compute delta and gamma for options. On the other hand, the finite difference approach is a very popular approach for computing vega and rho, as discussed in Hull [1] . We also compare to the existent other tree methods for computations of Greeks.

Figures 2-4 and Figure 6 plot the values of Greeks (delta, gamma, vega, and rho, respectively) for American put options computed using our approach and the finite difference approach. MS Greeks, extended binomial Greeks (EB Greeks) calculated by the extended binomial tree of Pelsser and Vorst [17] , Greeks introduced by Hull (Hull’s Greeks) are also plotted in Figure 2 and Figure 3. The extended (N-step) binomial tree is a $N+2$ binomial tree starting from $-2\Delta t$ , as shown in Figure 1. The EB delta and EB gamma are given by

${\Delta }_A^{EB}=\frac{\mathcal{A}\left(s\left(2,0\right),0\right)-\mathcal{A}\left(s\left(-2,0\right),0\right)}{s\left(2,0\right)-s\left(-2,0\right)}$

${\Gamma }_A^{EB}=\frac{\frac{\mathcal{A}\left(s\left(2,0\right),0\right)-\mathcal{A}\left(s\left(0,0\right),0\right)}{s\left(2,0\right)-s\left(0,0\right)}-\frac{\mathcal{A}\left(s\left(0,0\right),0\right)-\mathcal{A}\left(s\left(-2,0\right),0\right)}{s\left(0,0\right)-s\left(-2,0\right)}}{s\left(2,0\right)-s\left(-2,0\right)}$

where $s\left(i\mathrm{,}j\right)=s\ast u^i$ . These results are computed using binomial trees with $N=4,8,\cdots ,100$ steps for one year. The parameter values assumed for these numerical studies were

$s=100,\sigma =0.3,r=0.05,K=100,T=1\left(\text{year}\right).$

Figure 1. Extended binomial tree.

The MS Greeks (delta, gamma, vega, and rho) are represented by the real lines and the dotted line (without any mark) represents the finite difference Greeks (FD Greeks). Two other kinds of dotted lines, dotted lines with a circle and square, represent EB Greeks and Hull’s Greeks (delta and gamma), respectively. The horizontal lines in Figure 2 and Figure 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 Figure 4 and Figure 6 are FD vega and FD rho computed using the binomial trees with 100,000 steps for one year. Because we use very fine meshes for the computation of these horizontal lines, these numerical

Figure 2. Delta for American put options (MS delta, EB delta, FD delta, and Hull delta, K = 100) computed by binomial trees with $N=4,8,\cdots ,100$ steps.

Figure 3. Gamma for American put options (MS gamma, EB gamma, FD gamma, and Hull gamma, K = 100) computed by binomial trees with $N=4,8,\cdots ,100$ steps.

results are expected to be very accurate. If the initial underlying asset price $s$ is in the continuous region, i.e. $s\in {\mathcal{C}}_0$ , the FD delta and gamma are given by

${\Delta }_A^{FD}=\frac{\mathcal{A}\left(s+\Delta s,0\right)-\mathcal{A}\left(s-\Delta s,0\right)}{2\Delta s}$

${\Gamma }_A^{FD}=\frac{\mathcal{A}\left(s+\Delta s,0\right)-2\mathcal{A}\left(s,0\right)+\mathcal{A}\left(s-\Delta s,0\right)}{{\left(\Delta s\right)}^2}$

and the FD vega and rho are given by

${\mathcal{V}}_A^{FD}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma \right)=\frac{\mathcal{A}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma +\Delta \sigma \right)-\mathcal{A}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma -\Delta \sigma \right)}{2\Delta \sigma }$

${\rho }_A^{FD}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}r\right)=\frac{\mathcal{A}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}r+\Delta r\right)-\mathcal{A}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}r-\Delta r\right)}{2\Delta r}$

where $\Delta s$ , $\Delta \sigma$ , and $\Delta r$ are small parameters. We take $\Delta s=s\times \Delta h,\Delta \sigma =\sigma \Delta h$ , and $\Delta r=r\times \Delta h,\left(\Delta h={10}^{-3}\right)$ for our computations. Figure 2 shows that MS delta converges much faster than the FD delta. Figure 3 shows that FD gamma does not appear in the picture, and we do not recommend the use of the finite difference approach to compute gamma. Figure 4 reveals that MS vega conver- ges slower than FD vega. However, this is not a universal result. MS vega and FD vega for American put options are plotted in Figures 5 with strike prices 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

Figure 4. Vega for American put options (MS vega and FD gamma, K = 100) computed by binomial trees with $N=4,8,\cdots ,500$ steps.

Figure 5. Vega for American put options (MS vega and FD gamma, K = 105) computed by binomial trees with $N=4,8,\cdots ,500$ steps.

Figure 6. Rho for American put options (MS rho and FD rho, K = 100) computed by binomial trees with $N=4,8,\cdots ,500$ steps.

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

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 and FD rho computed by a binomial tree with 10,000 steps³. It was found that the computational time for MS 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 advanta- geous.

Figures 7-10 present Greeks (delta, gamma, vega, and rho, respectively) for

Figure 7. Delta for American put options as a function of the undelying asset price (MS delta, EB delta, and FD delta, K = 100).

Figure 8. Gamma for American put options as a function of the undelying asset price (MS gamma, EB gamma, and FD gamma, K = 100).

Table 1. Computational time for Rho.

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

Figure 9. Vega for American put options as a function of the undelying asset price (MS vega and FD vega, K = 100).

Figure 10. Rho for American put options as a function of the undelying asset price (MS rho and FD rho, K = 100).

and EB Greeks are very smooth, whereas those of FD delta and FD gamma are unstable. Figure 9 and Figure 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 the oscillation phenomenon is observed for FD vega. The oscillation phenomenon for FD vega is especially strong when the strike price is higher than the initial price of underlying assets. Figure 10 reveals that the numerical results of MS rho and FD rho are almost same.

Finally, we compared our new methods with other existing tree methods. We compare MS delta with delta for European options computed by other kinds of binomial tree, namely tree methods introduced by Chung and Shackelton [18] , Tian [19] , and Leisen and Reimer [20] . These are summarized in Figure 11. In order to obtain delta using the binomial trees of Chung and Shackelton [18] and Tian [19] , we employed the extended binomial tree approach. On the other hand, we used the finite difference approach to the binomial tree for Leisen and Reimer, because we wanted to implement simple calculations. Leisen and Reimer [20] introduced a new kind of binomial tree, which computes the price of options efficiently. They construct two kinds of trees using two different transform formulas. Note that because no significant difference is observed in two methods of Leisen and Reimer [20] , we used “Method-1” described in their article. As shown in Figure 11, Greeks calculated by trees introduced by Chung and Shackelton [18] and Tian [19] converges to the real value smoothly, however, MS delta converges faster than these methods. Delta computed by the tree introduced by Leisen and Reimer [20] 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 efficient algorithm to evaluate Greeks. It is especially efficient 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

Figure 11. Delta (MS delta and delta computed by various kinds of binomial trees.) computed by binomial trees with $N=4,8,\cdots ,100$ steps. MS: MS delta. CS: delta com- puted by binomial tree introduced by Chung and Shackelton [18] . Tian: delta computed by binomial tree introduced by Tian [19] . LR: delta computed by binomial tree intro- duced by Leisen and Reimer [20] .

computation of Greeks for American options. We introduce the binomial tree approach to overcome these problems and confirm its effectiveness for comput- ing Greeks for American options very quickly and accurately. Numerical results indicate that Greeks converge faster when computed using our method than when computed using the extended binomial tree approach of Pelsser and Vorst [17] . 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 in the jump diffusion models using the binomial tree approach (Muroi and Suda [9] and Suda and Muroi [15] ).

Acknowledgements

This work was supported by Japan Society for the Promotion of Science. Grant- in-Aid for Scientific Research (C) 16K03731.This paper was written while Shintaro Suda was a graduate student at Tohoku University. The views ex- pressed in this paper are those of the authors and do not necessarily reflect the official views of MTEC.

Appendix: Closed-Form Formulas for Option Greeks

We first show the error estimate of MS vega given by (16). This is presented in Section A. We also prove that MS Greeks converge to Greeks for continuous time Black and Scholes model. This is shown in Section B. Closed-form expectation formulas for MS Greeks (delta, gamma, vega, rho) for European options are investigated in Appendix B. We found that MS Greeks are approxi- mations of the discrete version of the Malliavin Greeks in the continuous time model and these results indicate that MS Greeks converge to Greeks for a continuous time model (Black and Scholes model) when we take a limit, $\Delta t\to 0$ .

A. Error terms for MS vega

We present theorem 2 again as theorem 3.

Theorem 3. Let us assume that the pay-off function $\Phi (\cdot )$ is a smooth function. The accuracy of the MS vega is given by

${\mathcal{V}}_E\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma \right)={\mathcal{V}}_E^{MS}\left(s\mathrm{,0\; <mo>\; ;\; </mo>}\sigma \right)+O\left(\sqrt{\Delta t}\right)\mathrm{.}$

[Proof]

Insert (12) to the Formula (15) leads to

$\begin{array}{c}{\mathcal{V}}_E\left(s_{i\Delta t},i\Delta t;\sigma \right)=E\left[e^{-r\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}\frac{{\epsilon }_{i+1}-\mu \Delta t}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t\right)\right]\\ +E\left[e^{-r\Delta t}{\Delta }_E^{MS}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}{ϵ}_{i+1}\right]\\ +E\left[e^{-r\Delta t}{\mathcal{V}}_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)\right]+O\left(\Delta t^{3/2}\right).\end{array}$

This formula yields

$\begin{array}{c}{\mathcal{V}}_E\left(s_{i\Delta t},i\Delta t;\sigma \right)={\mathcal{V}}_1^i\left(s_{i\Delta t},i\Delta t;\sigma \right)+{\mathcal{V}}_2^i\left(s_{i\Delta t},i\Delta t;\sigma \right)\\ +E\left[e^{-r\Delta t}{\mathcal{V}}_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)\right]+O\left(\Delta t^{3/2}\right),\end{array}$ (20)

where new variables ${\mathcal{V}}_1^i$ and ${\mathcal{V}}_2^i$ are defined by

${\mathcal{V}}_1^i\left(s_{i\Delta t},i\Delta t;\sigma \right)=E\left[e^{-r\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\mathcal{E}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t\right)|{\mathcal{F}}_{i\Delta t}\right]$

$\begin{array}{l}{\mathcal{V}}_2^i\left(s_{i\Delta t},i\Delta t;\sigma \right)\\ =E\left[e^{-r\Delta t}{\Delta }_E^{MS}\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}{ϵ}_{i+1}|{\mathcal{F}}_{i\Delta t}\right]\left(i\ne N-1\right)\end{array}$

$\begin{array}{l}{\mathcal{V}}_2^{N-1}\left(s_{\left(N-1\right)\Delta t},\left(N-1\right)\Delta t;\sigma \right)\\ =E\left[e^{-r\Delta t}\left\{{\frac{\text{d}}{\text{d}x}\Phi \left(x\right)|}_{x=s_{\left(N-1\right)\Delta t}{e}^{\sigma {ϵ}_N}}\right\}{s}_{\left(N-1\right)\Delta t}{e}^{\sigma {ϵ}_N}{ϵ}_N|{\mathcal{F}}_{\left(N-1\right)\Delta t}\right]\end{array}$

If the pay-off function $\Phi \left(x\right)$ is not smooth, we formally define ${\mathcal{V}}_2^{N-1}\left(s_{\left(N-1\right)\Delta t},\left(N-1\right)\Delta t;\sigma \right)=\mathrm{0,}$ although we assume that the pay-off function $\Phi \left(x\right)$ is smooth in this subsection. If the pay-off function $\Phi \left(x\right)$ is smooth, vega given by (20) is further computed as

$\begin{array}{c}{\mathcal{V}}_E\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)={\mathcal{V}}_1^i\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)+{\mathcal{V}}_2^i\left(s_{i\Delta t}\mathrm{,}i\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\\ +E\left[e^{-r\Delta t}{\mathcal{V}}_E\left(s_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}}\mathrm{,}\left(i+1\right)\Delta t\mathrm{<mo>\; ;\; </mo>}\sigma \right)\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_{i\Delta t}\right]+O\left(\Delta t^{3/2}\right)\mathrm{.}\end{array}$

Vega for the binomial tree model at the initial time is given by

$\begin{array}{c}{\mathcal{V}}_E\left(s,0;\sigma \right)={\mathcal{V}}_1^0\left(s,0;\sigma \right)+{\mathcal{V}}_2^0\left(s,0;\sigma \right)+O\left(\Delta t^{3/2}\right)+e^{-r\Delta t}E\left[{\mathcal{V}}_E\left(s{e}^{\sigma W_{\Delta t}},\Delta t;\sigma \right)\right]\\ ={\mathcal{V}}_1^0\left(s,0;\sigma \right)+{\mathcal{V}}_2^0\left(s,0;\sigma \right)+O\left(\Delta t^{3/2}\right)+e^{-r\Delta t}E[{\mathcal{V}}_1^1\left(s{e}^{\sigma W_{\Delta t}},\Delta t;\sigma \right)\\ +{\mathcal{V}}_2^1\left(s{e}^{\sigma W_{\Delta t}},\Delta t;\sigma \right)+O\left(\Delta t^{3/2}\right)+e^{-2\Delta t}{\mathcal{V}}_E\left(s{e}^{\sigma W_{2\Delta t}},2\Delta t;\sigma \right)]\\ =\cdots \\ ={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}}{e}^{-ri\Delta t}E\left[{\mathcal{V}}_1^i\left(s{e}^{\sigma W_{i\Delta t}},0;\sigma \right)+{\mathcal{V}}_2^i\left(s{e}^{\sigma W_{i\Delta t}},i\Delta t;\sigma \right)\right]+O\left(\sqrt{\Delta t}\right)\\ ={\mathcal{V}}_E^{MS}\left(s,0;\sigma \right)+O\left(\sqrt{\Delta t}\right).\end{array}$ (21)

This shows the result.

B. Convergence of MS greeks to black and scholes model

In Section 3, we derive Greeks under the Assumption 1. As discussed in Section 3, the smoothness of the pay-off function is too strong to be assumed. Therefore, in this section, we assume that the pay-off function $\Phi (\cdot )$ is not a necessarily smooth function. We show that MS Greeks obtained in Section 3 converge to Greeks in the continuous time model under the following assum- ption.

Assumption 2. We assume that the pay-off function, $\Phi (\cdot )$ , is a function in the class $\mathcal{K}$ . $\mathcal{K}$ is the class of real-valued functions on R that satisfy the following conditions: (i) $\varphi (\cdot )$ is piecewise $C^{\left(2\right)}$ , (ii) at each x, the function

$\varphi (\cdot )$ satisfies $\varphi \left(x\right)=\frac{1}{2}\left(\varphi \left(x+\right)+\varphi \left(x-\right)\right)$ , and (iii) $\varphi$ , ${\varphi }^{\prime }$ , and ${\varphi }^{″}$ are poly-

nomial bounded. We assume that $f(\cdot )$ is a function in the class $\mathcal{K}$ .

For example, the pay-off function for European call/put options is included in class $\mathcal{K}$ .

Theorem 4. We assume that the pay-off function is in class $\mathcal{K}$ . We also assume the number of steps of the binomial tree give by N to be even. Then, MS delta, gamma vega and rho converge to delta, gamma, vega, and rho in Black and Scholes model.

[Proof] We show that MS Greeks (delta, gamma, vega, and rho) are asympto- tically equivalent to the Malliavin Greeks. Malliavina Greeks are Greeks calcu- lated using the Malliavin calculus. See Kohatsu-Higa and Montero [3] for detail.

1. Delta

MS delta given by the Formula (9) is further calculated as

$\begin{array}{c}{\Delta }_E^{MS}\left(s,0\right)=\frac{e^{-r\Delta t}}{s\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ =\frac{e^{-rT}}{s\sigma \Delta t}E[\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_1-\mu \Delta t\right))]\\ =\frac{e^{-rT}}{s\sigma \Delta t}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_i-\mu \Delta t\right)\right],\end{array}$

where we used the fact that ${ϵ}_1\mathrm{,}\cdots \mathrm{,}{ϵ}_N$ is an i.i.d. sequence to deduce the last equality. This yields

$\begin{array}{c}{\Delta }_E^{MS}\left(s,0\right)=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}\frac{e^{-rT}}{s\sigma \Delta t}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_i-\mu \Delta t\right)\right]\\ =\frac{e^{-rT}}{s\sigma T}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T-\mu T\right)\right]\\ =\frac{e^{-rT}}{s\sigma T}E\left[\Phi \left(S_T\right)\frac{\log \left(S_T/s\right)-\mu \sigma T}{\sigma }\right],\end{array}$ (22)

where $S_T$ is given by $S_T=s{e}^{\sigma W_{N\Delta t}}$ . Note that this formula indicates that the MS delta is identical to the discrete Malliavin delta. See Kohatsu-Higa and Montero [3] about the Malliavin delta, for example. If the pay-off function $\Phi (\cdot )$ is smooth, proposition 2.1 in Heston and Zhou [21] leads to

$\begin{array}{c}{\Delta }_E^{MS}(s,0)=\frac{e^{-rT}}{s\sigma T}E\left[\Phi \left(P_T\right)\frac{\log \left(P_T/s\right)-\mu \sigma T}{\sigma }\right]+O\left(\frac{1}{N}\right)\\ =\frac{e^{-rT}}{s\sigma T}E\left[\Phi \left(P_T\right)Z_T\right]+O\left(\frac{1}{N}\right)\\ ={\Delta }_E^{BS}\left(s,0\right)+O\left(\frac{1}{N}\right)\end{array}$

⁴The order of error term is $O\left(\frac{1}{\sqrt{N}}\right)$ , if the discontinuity for the pay-off function $\Phi (\cdot )$ is not on a lattice point. However, if all discontinuities are on lattice points, the order of the error term is $O\left(\frac{1}{N}\right)$ . See Corollary 4.2 in Walsh [22] for details. Also note that Chung et al. [16] show that the rate of error terms of binomial delta for European options is $O\left(1/N\right)$ .

where the stochastic process $P_t$ is a geometric Brownian motion given in (1) and ${\Delta }_E^{BS}\left(s\mathrm{,0}\right)$ is delta in a continuous time model (Black and Scholes model). Even though MS delta does not approximate delta for the binomial tree model, if the pay-off function $\Phi (\cdot )$ is not smooth, it still is an approximation formula for the continuous time delta. Under Assumption 2, corollary 4.2 in Walsh [22] shows that the option price in the binomial tree model converges to the options price in the Black and Scholes model. Note that Walsh show the convergence only on the binomial tree with even numbers. This result shows⁴

$\begin{array}{c}{\Delta }_E^{MS}\left(s,0\right)=\frac{e^{-rT}}{s\sigma T}E\left[\Phi \left(P_T\right)\frac{\log \left(P_T/s\right)-\mu \sigma T}{\sigma }\right]+O\left(\frac{1}{\sqrt{N}}\right)\\ ={\Delta }_E^{BS}\left(s,0\right)+O\left(\frac{1}{\sqrt{N}}\right).\end{array}$

As previously discussed, MS delta does not approximate delta for the binomial tree model, if the pay-off function is not smooth. On the other hand, even if the pay-off function is not smooth, MS delta still is an approximation for continuous delta.

2. Gamma

MS gamma for European options is given by (14). Further computation of MS gamma yields,

$\begin{array}{c}{\Gamma }_E^{MS}\left(s,0\right)=-\frac{e^{-r\Delta t}}{s^2\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ +\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[{\Delta }_E^{MS}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right].\end{array}$ (23)

The first and second terms in (23) are denoted by $G_1$ and $G_2$ , i.e. ${\Gamma }_E^{MS}\left(s,0\right)=G_1+G_2$ :

$G_1=-\frac{e^{-r\Delta t}}{s^2\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)\left({ϵ}_1-\mu \Delta t\right)\right]$

$G_2=\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[{\Delta }_E^{MS}\left(s{e}^{\sigma {ϵ}_1},\Delta t\right)e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right].$

The first term is given by

$\begin{array}{c}{G}_1=-\frac{e^{-rT}}{s^2\sigma \Delta t}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_1-\mu \Delta t\right)\right]\\ =-\frac{e^{-rT}}{s^2\sigma \Delta tN}{\displaystyle \underset{i=1}{\overset{N}{\sum }}}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_i-\mu \Delta t\right)\right]\\ =-\frac{e^{-rT}}{s^2\sigma T}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T-\mu T\right)\right]\\ =-\frac{1}{s}{\Delta }_E^{MS}\left(s,0\right),\end{array}$

and the second term is given by

$\begin{array}{c}{G}_2=\frac{e^{-r\Delta t}}{\sigma s\Delta t}E\left[\frac{e^{-r\Delta t}}{s{e}^{\sigma {ϵ}_1}\sigma \Delta t}E\left[\mathcal{E}\left(s{e}^{\sigma \left({ϵ}_1+{ϵ}_2\right)},2\Delta t\right)\left({ϵ}_2-\mu \Delta t\right)|{\mathcal{F}}_{\Delta t}\right]e^{\sigma {ϵ}_1}\left({ϵ}_1-\mu \Delta t\right)\right]\\ =\frac{1}{s^2{\sigma }^2{\left(\Delta t\right)}^2}E\left[e^{-rN\Delta t}\Phi \left(s{e}^{\sigma W_T}\right)\left({ϵ}_1-\mu \Delta t\right)\left({ϵ}_2-\mu \Delta t\right)\right].\end{array}$

This formula is divided into three parts, $G_2=G_2^1-2G_2^2+G_2^3$ , where $G_2^1\mathrm{,}{G}_2^2$ and $G_2^3$ are

$\begin{array}{c}{G}_2^1=\frac{1}{s^2{\sigma }^2\Delta t^2}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_1{ϵ}_2\right]\\ =\frac{1}{s^2{\sigma }^2\Delta t^2}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_i{ϵ}_{j\left(\ne i\right)}\right]\\ =\frac{1}{s^2{\sigma }^2{T}^2}\frac{N}{N-1}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right)\left\{{\left({\displaystyle \underset{i}{\sum }}{ϵ}_i\right)}^2-N\Delta t\right\}\right]\\ =\frac{1}{s^2{\sigma }^2{T}^2}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T^2-T\right)\right]+O\left(1/N\right)\end{array}$

$G_2^2=\frac{1}{s^2{\sigma }^2\Delta t^2}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_1\mu \Delta t\right]=\frac{\mu }{s^2{\sigma }^{2}T}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right)W_T\right]$

$G_2^3=\frac{1}{s^2{\sigma }^2\Delta t^2}E\left[e^{-rT}\Phi \left(s{e}^{\sigma W_T}\right){\mu }^2\Delta t^2\right]=\frac{{\mu }^2}{s^2{\sigma }^2}E\left[e^{-rN\Delta t}\Phi \left(s{e}^{\sigma W_T}\right)\right].$

These results are combined into the closed-form formulas for MS Gamma:

${\Gamma }_E^{MS}\left(s,0\right)=\frac{e^{-rT}}{s^2\sigma T}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left\{\frac{{\left(W_T-\mu T\right)}^2}{\sigma T}-\left(W_T-\mu T\right)-\frac{1}{\sigma }\right\}\right]+O\left(1/N\right).$

As discussed in the delta case, Walsh [22] yields,

$\begin{array}{c}{\Gamma }_E^{MS}\left(s,0\right)=\frac{e^{-rT}}{s^2\sigma T}E\left[\Phi \left(S_T\right)\left\{\frac{{\left(\log \left(S_T/s\right)-\mu \sigma T\right)}^2}{{\sigma }^{3}T}-\frac{\log \left(S_T/s\right)-\mu \sigma T}{\sigma }-\frac{1}{\sigma }\right\}\right]+O\left(1/N\right)\\ =\frac{e^{-rT}}{s^2\sigma T}E\left[\Phi \left(P_T\right)\left\{\frac{{\left(\log \left(P_T/s\right)-\mu \sigma T\right)}^2}{{\sigma }^{3}T}-\frac{\log \left(P_T/s\right)-\mu \sigma T}{\sigma }-\frac{1}{\sigma }\right\}\right]+O\left(1/\sqrt{N}\right)\\ =\frac{e^{-rT}}{s^2\sigma T}E\left[\Phi \left(P_T\right)\left\{\frac{Z_T^2}{\sigma T}-Z_T-\frac{1}{\sigma }\right\}\right]+O\left(1/\sqrt{N}\right)\\ ={\Gamma }_E^{BS}\left(s,0\right)+O\left(1/\sqrt{N}\right),\end{array}$

even if $\Phi (\cdot )$ is not smooth.

3. Vega

MS Vega for the binomial tree model is given by (16). The expectation $E\left[{\mathcal{V}}_1^i\right]$ is given by

$E\left[{\mathcal{V}}_1^i\right]=e^{-r\left(N-i\right)\Delta t}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}E\left[\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\Phi \left(s{e}^{\sigma W_T}\right)\right].$

Under the condition $i\ne N-1$ , the expectation $E\left[{\mathcal{V}}_2^i\right]$ is given by

$\begin{array}{c}E\left[{\mathcal{V}}_2^i\right]=E[e^{-r\Delta t}\frac{e^{-r\Delta t}}{s{e}^{\sigma \left(W_{i\Delta t}+{ϵ}_{i+1}\right)}\sigma \Delta t}\\ \times E\left[\mathcal{E}\left(s{e}^{\sigma W_{i\Delta t}+{ϵ}_{i+1}}{e}^{\sigma {ϵ}_{i+2}},\left(i+2\right)\Delta t;\sigma \right)\left({ϵ}_{i+2}-\mu \Delta t\right)|{\mathcal{F}}_{\left(i+1\right)\Delta t}\right]\times s{e}^{\sigma \left(W_{i\Delta t}+{ϵ}_{i+1}\right)}{ϵ}_{i+1}]\\ =\frac{e^{-r\left(N-i\right)\Delta t}}{\sigma \Delta t}\left\{E\left[\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_{i+2}{ϵ}_{i+1}\right]-E\left[\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_{i+1}\right]\mu \Delta t\right\}.\end{array}$

Formulas

$\begin{array}{c}E\left[\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_{i+2}{ϵ}_{i+1}\right]=E\left[\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_i{ϵ}_{j\left(\ne i\right)}\right]\\ =\frac{1}{N^2-N}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left\{{\left({\displaystyle \underset{i=1}{\overset{N}{\sum }}}{ϵ}_i\right)}^2-N\Delta t\right\}\right]\\ =\frac{1}{N^2-N}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T^2-T\right)\right]\end{array}$

$E\left[\Phi \left(s{e}^{\sigma W_T}\right){ϵ}_{i+1}\right]=\frac{1}{N}E\left[\Phi \left(s{e}^{\sigma W_T}\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}}{ϵ}_i\right]=\frac{1}{N}E\left[\Phi \left(s{e}^{\sigma W_T}\right)W_T\right]$

lead to

$E\left[{\mathcal{V}}_2^i\right]=\frac{e^{-r\left(N-i\right)\Delta t}}{\sigma T}\left\{\frac{1}{N-1}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T^2-T\right)\right]-E\left[\Phi \left(s{e}^{\sigma W_T}\right)W_T\right]\mu \Delta t\right\}\left(i\ne N-1\right)$

If the pay-off function $\Phi (\cdot )$ is a smooth one, the expectation $e^{-r\left(N-1\right)\Delta t}E\left[{\mathcal{V}}_2^{N-1}\right]$ is calculated as

$\begin{array}{l}{e}^{-r\left(N-1\right)\Delta t}E\left[{\mathcal{V}}_2^{N-1}\right]\\ =e^{-rT}s{e}^{\sigma W_{\left(N-1\right)\Delta t}}\sqrt{\Delta t}\times \left[p{\Phi }^{\prime }\left(s{e}^{\sigma \left(W_{\left(N-1\right)\Delta t}+\sqrt{\Delta t}\right)}\right)e^{\sigma \sqrt{\Delta t}}-\left(1-p\right){\Phi }^{\prime }\left(s{e}^{\sigma \left(W_{\left(N-1\right)\Delta t}-\sqrt{\Delta t}\right)}\right)e^{-\sigma \sqrt{\Delta t}}\right]\\ =e^{-rT}s{e}^{\sigma W_{\left(N-1\right)\Delta t}}\Delta t\times {\frac{\partial }{\partial z}{\Phi }^{\prime }\left(s{e}^{\sigma \left(W_{\left(N-1\right)\Delta t}+z\right)}\right)e^{\sigma z}|}_{z=0}+O\left(\Delta t^{3/2}\right)\\ =e^{-rT}s\sigma e^{\sigma W_{\left(N-1\right)\Delta t}}\Delta t\times \left({\Phi }^{″}\left(s{e}^{\sigma W_{\left(N-1\right)\Delta t}}\right)s{e}^{\sigma W_{\left(N-1\right)\Delta t}}+{\Phi }^{\prime }\left(s{e}^{\sigma W_{\left(N-1\right)\Delta t}}\right)\right)+O\left(\Delta t^{3/2}\right)\\ =O\left(1/N\right),\end{array}$

If the pay-off function $\Phi \left(x\right)$ is not smooth, the relation

$e^{-r\left(N-1\right)\Delta t}E\left[{\mathcal{V}}_2^{N-1}\right]=0\left(=O\left(1/N\right)\right)\mathrm{,}$

still satisfied, because we formally assumed ${\mathcal{V}}_2^{N-1}=0$ . These results are combi- ned into

$\begin{array}{c}{\mathcal{V}}_E^{MS}\left(s,0;\sigma \right)={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}}{e}^{-rT}\left\{-\left(1+\frac{2r}{{\sigma }^2}\right)\right\}E\left[\frac{{ϵ}_{i+1}-\mu \Delta t}{2}\Phi \left(s{e}^{\sigma W_T}\right)\right]\\ +{\displaystyle \underset{i=0}{\overset{N-2}{\sum }}}\frac{e^{-rT}}{\sigma T}\left\{\frac{1}{N-1}E\left[\Phi \left(s{e}^{\sigma W_T}\right)\left(W_T^2-T\right)\right]-E\left[\Phi \left(s{e}^{\sigma W_T}\right)W_T\right]\mu \Delta t\right\}\\ =e^{-rT}E\left[\left(\frac{{\left(W_T-\mu T\right)}^2}{\sigma T}-\left(W_T-\mu T\right)-\frac{1}{\sigma }\right)\Phi \left(s{e}^{\sigma W_T}\right)\right]+O\left(\frac{1}{N}\right).\end{array}$

As discussed in delta case, the approximation

$\begin{array}{c}{\mathcal{V}}_E^{MS}\left(s,0;\sigma \right)=e^{-rT}E\left[\left(\frac{{\left(\log \left(S_T/s\right)-\mu \sigma T\right)}^2}{{\sigma }^{3}T}-\frac{\log \left(S_T/s\right)-\mu \sigma T}{\sigma }-\frac{1}{\sigma }\right)\Phi \left(S_T\right)\right]+O\left(\frac{1}{N}\right)\\ =e^{-rT}E\left[\left(\frac{{\left(\log \left(P_T/s\right)-\mu \sigma T\right)}^2}{{\sigma }^{3}T}-\frac{\log \left(P_T/s\right)-\mu \sigma T}{\sigma }-\frac{1}{\sigma }\right)\Phi \left(P_T\right)\right]+O\left(\frac{1}{\sqrt{N}}\right)\\ =e^{-rT}E\left[\left(\frac{Z_T^2}{\sigma T}-Z_T-\frac{1}{\sigma }\right)\Phi \left(P_T\right)\right]+O\left(\frac{1}{\sqrt{N}}\right)\\ ={\mathcal{V}}_E^{BS}\left(s,0;\sigma \right)+O\left(\frac{1}{\sqrt{N}}\right)\end{array}$

is valid, even if the pay-off function $\Phi (\cdot )$ is not smooth.

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

$E\left[\left(\frac{{ϵ}_i-\mu \Delta t}{\sigma }-\Delta t\right)\mathcal{E}\left(s{e}^{\sigma W_{i\Delta t}},\Delta t;r\right)\right]=e^{-r\left(N-i\right)\Delta t}E\left[\left(\frac{{ϵ}_i-\mu \Delta t}{\sigma }-\Delta t\right)\Phi \left(s{e}^{\sigma W_{N\Delta t}}\right)\right]$

leads the further calculation of rho. This formula is plugged into the Formula (18) and we have

${\rho }_E^{MS}\left(s,0;r\right)=e^{-rT}E\left[\left(\frac{W_T-\mu T}{\sigma }-T\right)\Phi \left(S_T\right)\right]=T\left(s{\Delta }_E^{MS}\left(s,0\right)-\mathcal{E}\left(s,0;r\right)\right).$

The last equality comes from the closed-form formula for MS delta given by (9) (or (22)). As is the discussions in the previous cases, the formula

${\rho }_E^{MS}\left(s,0;r\right)=T\left(s{\Delta }_E^{BS}\left(s,0\right)-E\left[e^{-rT}\Phi \left(P_T\right)\right]\right)+O\left(1/\sqrt{N}\right)={\rho }_E^{BS}\left(s,0\right)+O\left(1/\sqrt{N}\right)$

must be satisfied under certain conditions given by Walsh [22] .

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Hull, J. (2008) Options, Futures and Other Derivatives. 7th Edition, Prentice Hall, Upper Saddle River, New Jersey.
[2]	Fournié, E., Laszry, J., Lebuchoux, J., Lions, P. and Touzi, N. (1999) Applications of Malliavin Calculus to Monte Carlo Methods in Finance. Finance and Stochastics, 3, 391-412. https://doi.org/10.1007/s007800050068
[3]	Kohatsu-Higa, A. and Montero, M. (2003) Malliavin Calculus Applied to Finance. Physica A, 320, 548-570. https://doi.org/10.1016/S0378-4371(02)01531-5
[4]	Bernis, G., Gobet E. and Kohatsu-Higa, A. (2003) Monte Carlo Evaluation of Greeks for Multidimensional Barrier and Lookback Options. Mathematical Finance, 13, 99-113. https://doi.org/10.1111/1467-9965.00008
[5]	Gobet, E. (2004) Revisiting the Greeks for European and American Options. Proceedings of the International Symposium on Stochastic Processes and Mathematical Finance at Ritsumeikan University, Kusatsu, 5-9 March 2003, 53-71. https://doi.org/10.1142/9789812702852_0003
[6]	Bally, V., Caramellino, L. and Zanette, A. (2005) Pricing and Hedging American Options by Monte Carlo Methods Using a Malliavin Calculus Approach. Monte Carlo Methods and Applications, 11, 97-133. https://doi.org/10.1515/156939605777585944
[7]	Longstaff, F. and Schwartz, E.S. (2001) Valuing American Options by Simulation: A Simple Least Squares Approach. Review of Financial Studies, 14, 113-147. https://doi.org/10.1093/rfs/14.1.113
[8]	Muroi, Y. and Suda, S. (2013) Discrete Malliavin Calculus and Computations of Greeks in the Binomial Tree. European Journal of Operational Research, 231, 349-361. https://doi.org/10.1016/j.ejor.2013.05.038
[9]	Muroi, Y. and Suda, S. (2017) Computation of Greeks in the Jump-Diffusion Model using Discrete Malliavin Calculus. Mathematics and Computers in Simulation, 140, 69-93. https://doi.org/10.1016/j.matcom.2017.03.002
[10]	Cox, J.C., Ross, S.A. and Rubinstein, M. (1979) Option Pricing: A Simplified Approach. Journal of Financial Economics, 7, 229-263. https://doi.org/10.1016/0304-405X(79)90015-1
[11]	Leitz-Martini, M. (2000) A Discrete Clark-Ocone Formula. Maphysto Research Report, 29.
[12]	Privault, N. (2008) Stochastic Analysis of Bernoulli Processes. Probability Surveys, 5, 435-483. https://doi.org/10.1214/08-PS139
[13]	Privault, N. (2009) Stochastic Analysis in Discrete and Continuous Settings. Springer, Berlin. https://doi.org/10.1007/978-3-642-02380-4
[14]	Privault, N. and Schoutens, W. (2002) Discrete Chaotic Calculus and Covariance Identities. Stochatics and Stochastic Reports, 72, 289-315. https://doi.org/10.1080/10451120290019230
[15]	Suda, S. and Muroi, Y. (2015) Computation of Greeks Using Binomial Trees in a Jump-Diffusion Model. Journal of Economic Dynamics and Control, 51, 93-110. https://doi.org/10.1016/j.jedc.2014.09.032
[16]	Chung, S.L., Hung, W., Lee, H.H. and Shih, P.T. (2011) On the Rate of Convergence of Binomial Greeks. Journal of Futures Markets, 31, 562-597. https://doi.org/10.1002/fut.20484
[17]	Pelsser, A. and Vorst, T. (1994) The Binomial Model and the Greeks. Journal of Derivatives, 1, 45-49. https://doi.org/10.3905/jod.1994.407888
[18]	Chung, S.L. and Shackleton, M. (2002) The Binomial Black-Scholes Model and the Greeks. Journal of Futures Markets, 22, 143-153. https://doi.org/10.1002/fut.2211
[19]	Tian, Y. (1993) A Modified Lattice Approach to Option Pricing. Journal of Futures Markets, 13, 563-577. https://doi.org/10.1002/fut.3990130509
[20]	Leisen, D. and Reimer, M. (1996) Binomial Models for Option Valuation-Examining and Improving Convergence. Applied Mathematical Finance, 3, 319-346. https://doi.org/10.1080/13504869600000015
[21]	Heston, S. and Zhou, G. (2000) On the Rate of Convergence of Discrete-Time Contingent Claims. Mathematical Finance, 10, 53-75. https://doi.org/10.1111/1467-9965.00080
[22]	Walsh, J. (2003) The Rate of Convergence of the Binomial Tree Scheme. Finance and Stochastics, 7, 337-361. https://doi.org/10.1007/s007800200094