. Let us compute the expec- tation $E\left[\mathcal{E}\left(s{e}^{\sigma {ϵ}_{1}},\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 \left(\cdot \right)$ as a smooth function. We can estimate the accuracy of MS delta as

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

[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}^{±\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=±\sqrt{\Delta t}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\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=±\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=±\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}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}±{\frac{{\partial }^{2}}{\partial {z}^{2}}\mathcal{E}\left(s{e}^{\sigma z},\Delta t\right)\sqrt{\Delta t}|}_{z=±\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,0\right)={\Delta }_{E}^{MS}\left(s,0\right)+O\left(\Delta t\right).$ (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,\Delta t\right)|}_{x=s{e}^{±\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}^{±\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}^{±\sigma \sqrt{\Delta t}},\Delta t\right)}{s{e}^{\mp \sigma \sqrt{\Delta t}}-s{e}^{±\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,0\right)\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}}.$

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}^{±\sigma \sqrt{\Delta t}},\Delta t\right)$ given by (5) and (13) yield two

different approximation formulas for delta1. 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,0\right)={\Delta }_{E}^{D}\left(s,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

1Actually, MS delta is deeply related to Hull’s delta. We have a relation, ${\Delta }_{E}^{MS}\left(s,0\right)={e}^{-r\Delta t}{\Delta }_{E}^{Hull}\left(s,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,0\right)={e}^{-r\Delta t}\left\{p\frac{\partial }{\partial x}\mathcal{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\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}},\Delta t\right){e}^{-\sigma \sqrt{\Delta t}}\right\}\\ ={e}^{-r\Delta t}\left\{p{\Delta }_{E}\left(s{e}^{\sigma \sqrt{\Delta t}},\Delta t\right){e}^{\sigma \sqrt{\Delta t}}+\left(1-p\right){\Delta }_{E}\left(s{e}^{-\sigma \sqrt{\Delta t}},\Delta t\right){e}^{-\sigma \sqrt{\Delta t}}\right\}.\end{array}$

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

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

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}\left\{{p\frac{\partial }{\partial z}\left({\Delta }_{E}\left(s{e}^{\sigma z},\Delta t\right){e}^{\sigma z}\right)|}_{z=\sqrt{\Delta t}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\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}}\right\}.\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=±\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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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 \left(\cdot \right)$ 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,0\right)={\Gamma }_{E}^{MS}\left(s,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},i\Delta t;\sigma \right)$ and ${\mathcal{V}}_{E}\left({s}_{i\Delta t},i\Delta t;\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},i\Delta t;\sigma \right)={e}^{-r\Delta t}\frac{\partial }{\partial \sigma }\left\{p\mathcal{E}\left({s}_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(1-p\right)\mathcal{E}\left({s}_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right)\right\}\\ ={\mathcal{V}}_{1}^{i}+{\mathcal{V}}_{2}^{i}+{\mathcal{V}}_{3}^{i},\end{array}$

where ${\mathcal{V}}_{1}^{i},{\mathcal{V}}_{2}^{i},{\mathcal{V}}_{3}^{i}$ are given by

$\begin{array}{l}{\mathcal{V}}_{1}^{i}={e}^{-r\Delta t}\left\{\frac{\partial p}{\partial \sigma }\mathcal{E}\left({s}_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right)\\ \text{ }\text{ }+\frac{\partial \left(1-p\right)}{\partial \sigma }\mathcal{E}\left({s}_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right)\right\}\end{array}$

$\begin{array}{c}{\mathcal{V}}_{2}^{i}={e}^{-r\Delta t}\left\{p\frac{\partial }{\partial x}\mathcal{E}\left({s}_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right){s}_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}}\sqrt{\Delta t}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(1-p\right)\frac{\partial }{\partial x}\mathcal{E}\left({s}_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right){s}_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}}\left(-\sqrt{\Delta t}\right)\right\}\\ ={e}^{-r\Delta t}E\left[{\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]\end{array}$

${\mathcal{V}}_{3}^{i}={e}^{-r\Delta t}E\left[{\mathcal{V}}_{E}\left({s}_{i\Delta t}{e}^{\sigma {ϵ}_{i+1}},\left(i+1\right)\Delta t;\sigma \right)\right].$

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}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}×\frac{\mathcal{E}\left({s}_{i\Delta t}{e}^{\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\sigma \right)\sqrt{\Delta t}+\mathcal{E}\left({s}_{i\Delta t}{e}^{-\sigma \sqrt{\Delta t}},\left(i+1\right)\Delta t;\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}},\left(i+1\right)\Delta t;\sigma \right)\right]+O\left(\Delta {t}^{3/2}\right).\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},N\Delta t;\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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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 \left(\cdot \right)$ is a smooth function. The accuracy of the MS vega is given by

${\mathcal{V}}_{E}\left(s,0;\sigma \right)={\mathcal{V}}_{E}^{MS}\left(s,0;\sigma \right)+O\left(\sqrt{\Delta t}\right).$

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 \left(\cdot \right)$ 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},i\Delta t;r\right)$ and ${\rho }_{E}\left({s}_{i\Delta t},i\Delta t;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},i\Delta t;r\right)={\rho }_{1}^{i}+{\rho }_{2}^{i}+{\rho }_{3}^{i},$

where ${\rho }_{1}^{i},{\rho }_{2}^{i},{\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 }\left\{\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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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)\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}},\left(i+1\right)\Delta t;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[\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)\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[\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)\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 \\ =\underset{i=1}{\overset{N}{\sum }}\text{ }\text{ }{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,0;r\right)={\rho }_{E}^{MS}\left(s,0;r\right)+O\left(\sqrt{\Delta t}\right).$

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 \left(\cdot \right)$ 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,i\Delta t\right)$ . The price of American options is given by the backward induction algorithm:

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

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

${\mathcal{S}}_{i}=\left\{x|x\in \left(0,\infty \right),\text{\hspace{0.17em}}\mathcal{A}\left(x,i\Delta t\right)=\Phi \left(x\right)\right\}$ and ${\mathcal{C}}_{i}=\left(0,\infty \right)\{\mathcal{S}}_{i}\text{\hspace{0.17em}}\left(i=0,1,\cdots ,N\right)$

The continuous region and stopping region for American options are defined by $\mathcal{C}={\cup }_{i=0}^{N}{\mathcal{C}}_{i}$ and $\mathcal{S}={\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,0\right)={\partial }_{\theta }\Phi \left(s\right)$

if the function $\Phi \left(\cdot \right)$ is differentiable at $s\in {\mathcal{S}}_{0}$ . (Second derivatives such as gamma is given by ${\partial }_{ss}\mathcal{A}\left(s,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 by2

${\partial }_{\theta }\mathcal{A}\left(s,0\right)=\frac{\partial }{\partial \theta }\left[{e}^{-r\Delta t}\left\{p\mathcal{A}\left(su,\Delta t\right)+\left(1-p\right)\mathcal{A}\left(sd,\Delta t\right)\right\}\right].$ (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,0\right)=\left\{\begin{array}{ll}{\partial }_{s}\Phi \left(s\right)\hfill & \text{ }\text{if}\text{\hspace{0.17em}}s\in {\mathcal{S}}_{0}\text{ }\hfill \\ \frac{{e}^{-r\Delta t}}{s\sigma \Delta t}E\left[\mathcal{A}\left(s{e}^{\sigma {ϵ}_{1}},\Delta t\right)\left({ϵ}_{1}-\mu \Delta t\right)\right]\hfill & \text{ }\text{if}\text{\hspace{0.17em}}s\in {\mathcal{C}}_{0}.\hfill \end{array}$

2Strictly 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,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,\text{ }\sigma =0.3,\text{ }r=0.05,\text{ }K=100,\text{ }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,0;\sigma \right)=\frac{\mathcal{A}\left(s,0;\sigma +\Delta \sigma \right)-\mathcal{A}\left(s,0;\sigma -\Delta \sigma \right)}{2\Delta \sigma }$

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

where $\Delta s$ , $\Delta \sigma$ , and $\Delta r$ are small parameters. We take $\Delta s=s×\Delta h,\Delta \sigma =\sigma \Delta h$ , and $\Delta r=r×\Delta h,\text{\hspace{0.17em}}\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.

3It 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 steps3. 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 \left(\cdot \right)$ is a smooth function. The accuracy of the MS vega is given by

${\mathcal{V}}_{E}\left(s,0;\sigma \right)={\mathcal{V}}_{E}^{MS}\left(s,0;\sigma \right)+O\left(\sqrt{\Delta t}\right).$

[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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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]\text{ }\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)=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},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)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+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)|{\mathcal{F}}_{i\Delta t}\right]+O\left(\Delta {t}^{3/2}\right).\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\left[{\mathcal{V}}_{1}^{1}\left(s{e}^{\sigma {W}_{\Delta t}},\Delta t;\sigma \right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\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)\right]\\ =\cdots \\ =\underset{i=0}{\overset{N-1}{\sum }}\text{ }\text{ }{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 \left(\cdot \right)$ 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 \left(\cdot \right)$ , 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 \left(\cdot \right)$ is piecewise ${C}^{\left(2\right)}$ , (ii) at each x, the function

$\varphi \left(\cdot \right)$ 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\left(\cdot \right)$ 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\left[\Phi \left(s{e}^{\sigma {W}_{T}}\right)\left({ϵ}_{1}-\mu \Delta t\right)\right)\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},\cdots ,{ϵ}_{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}\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{\mathrm{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 \left(\cdot \right)$ is smooth, proposition 2.1 in Heston and Zhou [21] leads to

$\begin{array}{c}{\Delta }_{E}^{MS}\left(s,0\right)=\frac{{e}^{-rT}}{s\sigma T}E\left[\Phi \left({P}_{T}\right)\frac{\mathrm{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}$

4The order of error term is $O\left(\frac{1}{\sqrt{N}}\right)$ , if the discontinuity for the pay-off function $\Phi \left(\cdot \right)$ 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,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 \left(\cdot \right)$ 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 shows4

$\begin{array}{c}{\Delta }_{E}^{MS}\left(s,0\right)=\frac{{e}^{-rT}}{s\sigma T}E\left[\Phi \left({P}_{T}\right)\frac{\mathrm{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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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}\underset{i=1}{\overset{N}{\sum }}\text{ }\text{ }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}-2{G}_{2}^{2}+{G}_{2}^{3}$ , where ${G}_{2}^{1},{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(\underset{i}{\sum }\text{ }\text{ }{ϵ}_{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(\mathrm{log}\left({S}_{T}/s\right)-\mu \sigma T\right)}^{2}}{{\sigma }^{3}T}-\frac{\mathrm{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(\mathrm{log}\left({P}_{T}/s\right)-\mu \sigma T\right)}^{2}}{{\sigma }^{3}T}-\frac{\mathrm{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 \left(\cdot \right)$ 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\left[{e}^{-r\Delta t}\frac{{e}^{-r\Delta t}}{s{e}^{\sigma \left({W}_{i\Delta t}+{ϵ}_{i+1}\right)}\sigma \Delta t}\\ ×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]×s{e}^{\sigma \left({W}_{i\Delta t}+{ϵ}_{i+1}\right)}{ϵ}_{i+1}\right]\\ =\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(\underset{i=1}{\overset{N}{\sum }}\text{ }\text{ }{ϵ}_{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)\underset{i=1}{\overset{N}{\sum }}\text{ }\text{ }{ϵ}_{i}\right]=\frac{1}{N}E\left[\Phi \left(s{e}^{\sigma {W}_{T}}\right){W}_{T}\right]$

$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\}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(i\ne N-1\right)$

If the pay-off function $\Phi \left(\cdot \right)$ 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}×\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×{\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×\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\text{\hspace{0.17em}}\left(=O\left(1/N\right)\right),$

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)=\underset{i=0}{\overset{N-1}{\sum }}\text{ }\text{ }{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]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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(\mathrm{log}\left({S}_{T}/s\right)-\mu \sigma T\right)}^{2}}{{\sigma }^{3}T}-\frac{\mathrm{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(\mathrm{log}\left({P}_{T}/s\right)-\mu \sigma T\right)}^{2}}{{\sigma }^{3}T}-\frac{\mathrm{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 \left(\cdot \right)$ 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.

 [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