The value of an option is the risk-neutral expectation of its discounted payoff. **Monte Carlo** methods approximate the expectation operator by an average of a large number of simulated values. For instance, suppose a random variable, $x$ is distributed as $\mathcal{N}(\mu, \sigma^2)$. If we take $M$ draws from this distribution and average them we get an estimate for the expected value:
$$\mathbb{E}[x] \approx \frac{1}{M}\sum_{j=1}^{M} x_j $$

This works well if we are willing to make an assumption about the distribution. You may recall that we use bootstrapping methods when we are not willing to make such an assumption (through resampling).

As a first step to using **Monte Carlo** methods to price options we want to be able to simulate the path of the underlying security, $\{S_t\}_{t=0}^T$. Suppose this process takes the form:
$$dS_t = (r - \delta) S_t dt + \sigma S_t dz_t $$
Where $dz_t$ is a Brownian Motion increment. The above process is a **Geometric Brownian Motion**. In this case, the solution to the above stochastic differential euqation is:

$$S_T = S_0 e^{(r - \delta - \frac{1}{2}\sigma^2)T + \sigma z_T}$$

$$\mathbb{E}[S_T] = S_0 e^{(r - \delta - \frac{1}{2}\sigma^2)T} $$

However, for general cases without a known solution it is useful to be able to simulate the above stochastic process. First, we use **Ito's Lemma** for the process for $x = f(S) = \ln S$.

$$ dx_t = \frac{\partial f(S)}{\partial s}dS + \frac{1}{2}\frac{\partial^2 f(S)}{\partial S^2}dS^2 = (r - \delta - \frac{1}{2}\sigma^2)dt + \sigma dz_t = \nu dt + \sigma dz_t$$

This has been expressed in terms of the infinitesimals, $dz,dt,dx$. We change this to discrete time notation:

$$\Delta x = \nu \Delta t + \sigma \Delta z$$

Which we can express as:

$$x_{t+\Delta t} = x_t + \nu \Delta t + \sigma (z_{t+\Delta t} - z_t) $$

Where $z_{t+\Delta t} - z_t$ is distributed $\mathcal{N}(0, \Delta t)$. Thus, we can draw a sequence of random normals and simulate 100 sequences as follows:

In [43]:

```
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as ss
%matplotlib inline
def sample_path(N = 365, T = 1.0, sig = 0.2, r = .06, delta = .02, S = 100.0):
dt = 1.0 * T / N
x = np.zeros(N+1)
s = np.zeros(N+1)
eps = sig* np.sqrt(dt) * ss.norm.rvs(size = N, loc = 0, scale = 1)
nu = r - delta - 0.5 * sig ** 2
x[0] = np.log(S)
for t in range(N):
x[t+1] = x[t] + nu * dt + eps[t]
s = np.exp(x)
return x, s
```

In [44]:

```
sumS = 0.0
M = 100
for j in range(M):
X,S = sample_path()
sumS = sumS + S[-1]
plt.plot(S)
plt.xlabel('Time')
plt.ylabel('Stock Price')
plt.title('Simulated Sample Paths of a Geometric Brownian Motion')
plt.show()
meanST = sumS / M
expST = 100 * np.exp((.06 - .02 - .5*.2**2)*1)
print 'The mean of the terminal Values and the Expected Value'
print meanST, expST
```

We see that the simulations tend to their expected value, which is promising. As well, this gives us a sense of what the stochastic process looks like.

Monte Carlo Methods have lots to recommend them for European option pricing problems. We can easily price options with complex underlying stochastic processes by simulating the sample paths and averaging the discounted payoffs. We will go through the example of a European Call Option. For the assignment problem I will have you guys simulate a different stochastic process than the Geometric Brownian Motion.

In the case of a GBM we know that at time $T$ the stock price is distributed as $\mathcal{N}(r - \delta - \frac{1}{2}\sigma^2)T, \sigma^2 T)$. This makes simulation easier and we use this fact in the function below. For general stochastic processes we will use the *SamplePath* function to find the time $T$ stock price.

- For $m = 1, 2, \dots, M$ find the time $T$ stock price, $S_T$.
- For each $S_T^m$ find the option payoff: $\max\{0, S_T - K\}$.
- Average the discounted payoffs: $\frac{1}{M}\sum_{m=1}^{M}e^{-rT}C_T^m$

We may also want to find a measure of how accurate the simulation is. We use the **numerical standard error** for this. We define the numerical standard error as:

$$NSE = \frac{SD(C_T)}{\sqrt{M}} $$

So long as the variance of the call price simulations is finite we can make the error due to sampling variation arbitrarily small by increasing $M$. However, it may take a very large $M$ (i.e. $M > 1,000,000$) to make the NSE very small. This will be apparent in the example below. For this reason **Clemlow and Strickland** discuss various variance reduction techniques. We will not get heavily into these issues, but you should be aware.

We are now ready to define a function which prices an option via Monte Carlo methods.

In [70]:

```
def SamplePath(N = 10, T = 1.0, r = .06, div = .03, sig = .2, S = 100.):
'''
Simulates the sample path of the process:
dS = nu S dt + sigma S dz
With parameters:
N - numer of time intervals
T - time until exercise date
r - risk free interest rate
div - the dividend payout rate
sig - the standard deviation or diffusion term
S - initial security price
Returns:
S_T - the value of the stock at time T.
'''
dt = 1.0 * T / N
nu = r - div - 0.5 * sig ** 2
x = np.log(S)
for t in range(N):
eps = sig * np.sqrt(dt) * ss.norm.rvs(size = 1, loc = 0, scale = 1)
x = x + nu * dt + eps
s = np.exp(x)
return s
def EuroCall(K = 100.0, T = 1.0, S = 100., sig = .2, r = .06, div = .03, N = 10, M = 100, GBM = True):
'''
Uses Monte Carlo to price a European Call Option. There are two options for the underlying stochastic process:
- GBM == True implies that we can draw from a lognormal distribution.
- GBM == False implies we have to use the SamplePath function to find the time T stock price.
Parameters:
N - numer of time intervals
T - time until exercise date
r - risk free interest rate
div - the dividend payout rate
sig - the standard deviation or diffusion term
S - initial security price
K - the strike price
M - number of Monte Carlo simulations
Output:
0 - The Call price
1 - The Numerical Standard Deviation
'''
dt = (1.0 * T) / N
nudt = (r - div - .5 * sig ** 2 ) * dt
sigdt = sig * np.sqrt(dt)
lnS = np.log(S)
sumCT = 0.
sumCT2 = 0.
nuT = (r - div - 0.5 * sig ** 2) * T
sigT = sig * np.sqrt(T)
if GBM:
for j in range(M):
lnSt = ss.norm.rvs(size = 1, loc = nuT, scale = sigT)
ST = np.exp(lnSt + np.log(S))
CT = max(0, ST - K)
sumCT = sumCT + CT
sumCT2 = sumCT2 + CT * CT
else:
for j in range(M):
ST = SamplePath(N, T, r, div, sig, S)
CT = max(0, ST - K)
sumCT = sumCT + CT
sumCT2 = sumCT2 + CT * CT
callprice = np.exp(-r * T) * sumCT / M
SD = np.sqrt( (sumCT2 - sumCT2 / M) * np.exp(-2*r*T) / (M-1))
SE = SD / np.sqrt(M)
return callprice, SE
Cprice, Sterr = EuroCall()
print 'The Call Price is : ', Cprice
print 'The Numerical Standard Error is: ', Sterr
```

You will notice that I took the time to introduce what is called a *docstring* to the two functions above. This is good practice for functions that you will reuse often as you may forget how they work. Notice that when you call the python help function it will return the *docstring*. It is good practice to properly comment your functions.

In [71]:

```
help(SamplePath)
help(EuroCall)
```

Valuing American style options via Monte Carlo methods is more challenging due to early exercise. If there are no dividends then we can use the fact that it is not optimal to exercise a call early to simplify the problem (for call options only!). However, the beauty of the Monte Carlo approach is that we can flexibly account for dividends in the stochastic process for the stock price. This introduces the possibility of early exercise, which is problematic. The chapter on Monte Carlo Methods in Clemlow and Strickland is silent on how to handle American Options. I found a couple papers by students that summarize approaches to the problem that may be useful as starting points for those interested.

- http://www.math.wustl.edu/~feres/Math350Fall2012/Projects/mathproj07.pdf
- http://www2.math.uu.se/research/pub/Jia1.pdf

You may find this type of report at a more accessible level. If you really wanted to get into this material then the references cited within could be useful.

The assignment problem is to:

Write a function which simulates the sample path for the following stochastic process:

$$dS = \theta(r - \delta - S)dt + \sigma dz $$

Where $dz$ is a

**Brownian Motion**increment.Plot 100 sample paths of the process and compare the results to the

**Geometric Brownian Motion**presented in the lecture notes. What is a key difference between the two processes?*Hint: This is an Ornstein-Uhlenbeck Process and it has a nice Wikipedia article! Think about what happens to the variance as $T\to \infty$.*Write a function that uses the sample path function to generate 200 realizations of the process and return the Monte Carlo estimate of the Call Price as well as the numerical standard deviation of the estimate.

**The following are bonus questions.**Now write a new sample path function that simulates the following

**Cox-Ingersoll-Ross**process, which is used to model interest rates:$$dS = \theta(r - \delta - S)dt + \sigma\sqrt{S} dz $$

Simulate its sample paths and compare it to the previous two stochastic processes. It should be easy to amend your call pricing function to price a security with this underlying process. Do so.

Write a loop that uses the call pricing function and any of the above stochastic processes to calculate the price and numerical standard deviation of a European Call for $M = \{200, 300, 500\}$. Plot the estimated price and standard deviation. What do you notice?