By: Matt Brigida

In this post we’ll calculate the value of a European call option on a non-dividend paying stock using both Black-Scholes and Monte Carlo. Both valuation methods should give you the same answer, so this is a good first exercise in Monte Carlo methods. Moreover, it illustrates the point that Monte Carlo should be used only when an analytical solution is impractical.

Here are the parameters for the call option.

```
S <- 50
K <- 45
r <- 0.01
vol <- 0.2
T <- 0.5
```

We first value the call option by Black Scholes. The equations is:

\[ \begin{aligned} C = S_0N(d_1)-e^{-rT}KN(d_2) \end{aligned} \] \[ \begin{aligned} d_1 = \frac{ln\frac{S_0}{K} + \left(r + \frac{\sigma^2}{2}\right)T}{\sigma\sqrt{T}} \end{aligned} \] \[ \begin{aligned} d_2 = d_1 - \sigma\sqrt{T} \end{aligned} \]

where \(N()\) denotes the cumulative *standard* normal distribution function, and we assume today is time 0, that is \(t = 0\) so \(T-t = T\).

And the option value is:

```
C <- S*pnorm((1/(vol*sqrt(T)))*(log(S/K)+(r+vol*vol/2)*T))-K*exp(-r*T)*pnorm((1/(vol*sqrt(T)))*(log(S/K)+(r+vol*vol/2)*T)-vol*sqrt(T))
print(C)
```

`## [1] 6.056`

Now we’ll value the same call option by Monte Carlo. To do so we’ll need to calculate the probability density function (pdf) of \(S_T\). Because European options are not path dependent, we do not need to simulate the paths of the stock prices.

To find the pdf, assume the stock price follows the stochastic differential equation \(dS_t = rS_tdt + \sigma S_tdB_t\) where \(B_t\) is a Brownian motion process. Let \(Y(S_t, t) \equiv ln(S_t)\) then by Ito’s Lemma (informally written \(dY = \frac{\partial Y}{\partial S_t}dS_t + \frac{\partial Y}{\partial t}dt + \frac{1}{2}\frac{\partial^2 Y}{\partial S^2_t}\) ) we have:

\[ \begin{aligned} dY = \frac{1}{S_t}\left(rS_tdt + \sigma S_tdB_t\right) - \frac{1}{2 S^2_t}\left(rS_tdt + \sigma S_tdB_t\right)^2 \end{aligned} \] \[ \begin{aligned} \Rightarrow dY \equiv d ln(S_t) = \left(r - \frac{1}{2}\sigma^2\right)dt + \sigma dB_t \end{aligned} \]

Then integrating through from \(0\) to \(T\) affords:

\[ \begin{aligned} ln\left(\frac{S_T}{S_0}\right) = \left(r-\frac{1}{2}\sigma^2\right)t + \sigma B_t \end{aligned} \]

and finally taking the antilog we have:

\[ \begin{aligned} S_T = S_0 e^{\left(r-\frac{1}{2}\sigma^2\right)T + \sigma B_T} \end{aligned} \]

Now, note the the Brownian motion \(B_T\) is normally distributed with a mean of 0 and variance \(T\). So we take random draws from \(N(0,T)\) and plug into the above equation to get random draws from \(S_T\). We then take \(max(S_T - K, 0)\) for each draw, average this over all draws, and discount to time \(0\). The code is below.

```
z <- rnorm(1000000)
C <- exp(-r*T)*mean(ifelse(S*exp((r-.5*vol^2)*T+vol*sqrt(T)*z)>K, S*exp((r-.5*vol^2)*T+vol*sqrt(T)*z)-K, 0))
print(C)
```

`## [1] 6.061`

We can also recalculate the value by Monte Carlo many times, to get the density (and standard deviation) of the Monte Carlo estimator. This gives us an idea of the precision of the calculation. In the code below I repeat the Monte Carlo option valuation 200 times, and then print the mean, standard deviation, and plot the density.

```
call <- 0
for (i in 1:200){
z <- rnorm(1000000)
call[i] <- exp(-r*T)*mean(ifelse(S*exp((r-.5*vol^2)*T+vol*sqrt(T)*z)>K, S*exp((r-.5*vol^2)*T+vol*sqrt(T)*z)-K, 0))
}
```

`mean(call)`

`## [1] 6.056`

`sd(call)`

`## [1] 0.006231`

As you can see, on average the Monte Carlo estimate varies by about 0.6 cents.

`plot(density(call),main='Empirical Density of the Monte Carlo Values \n of a European Call', xlab='Call Premium ($)')`