import numpy as np
import scipy
import matplotlib.pyplot as plt
Question 1: Probability & Monte Carlo Simulation¶
(a)
Consider a company “XYZ capital” currently trading at $100$ (at time = $t_0$) on the National Stock Exchange of Vol-land and it moves up by $1 or moves down by the same amount with equal probability at each minute.
Let the starting value be denoted as $S_0$. For each $t$ (unit: minute), we define the random variable: $$ Z_t = \begin{cases} 1, \quad \text{w.p. }\frac12\\ -1, \quad \text{w.p. }\frac12 \end{cases} $$ which are mutually independent and identically distributed for each $t$.
Then let $S_t = S_0 + \sum_{k=1}^tZ_k$, we have: $$ \mathbb{E}[S_t] = S_0 + \sum_{k=1}^t\mathbb{E}[Z_k] = S_0 $$ since $\mathbb{E}[Z_k] = \frac12\cdot 1 - \frac12\cdot 1$ for each $k\ge 1$. Therefore, we should expect that the values to be the same for any amount of time $t$.
(b)
We are now interested in the probability of the stock price ending at a level $a$ before $b$. Let $T_a := \inf\{t>0: S_t=a\}$ and $T_b := \inf\{t>0: S_t=b\}$, we are interested in the probability of the event $\{T_a<T_b\}$.
Define $T = \inf\{t>0: S_t=a\text{ or } S_t=b\}$. Using the law of total expectations: $$ \mathbb{E}[S_{T}] = \mathbb{E}[\mathbb{E}[S_{T}|T=t]] = S_0 $$ since the $\mathbb{E}[S_t]=S_0$ holds for any $t$. On the other hand, there are only two cases, $T=T_a$ or $T=T_b$, we have: $$ \mathbb{E}[S_T] = S_0 = \mathbb{P}\{T_a<T_b\}\cdot a + (1-\mathbb{P}\{T_a<T_b\})\cdot b $$ which implies: $$ \mathbb{P}\{T_a<T_b\} = \frac{S_0-b}{a-b} $$
We now numerically approximate this probability value and investigate convergence rate.
def radamacher(n):
# defines a random variable with
# +/-1 outcomes of equal probability
s = np.random.rand(n)
return np.sign((s < 0.5)-0.5)
def simulate_one_path(x0, a, b):
# simulates one path and save
x = []
x.append(x0)
while x[-1] < a and x[-1] > b:
x.append(x[-1] + radamacher(1).item())
return x
def simulate_stock_price_stopping(x0, a, b, n):
# given an initial stock price and boundary values,
# output probability of reaching `a` as a proportion
# of paths that satisfy this criterion (the other is
# reaching `b`).
i = 0
# number of paths stopped at a
count = 0
while i < n:
x = simulate_one_path(x0, a, b)
if x[-1] == a:
count = count + 1
i = i + 1
return count / n
# fix seed
np.random.seed(123)
# plot a few paths
x0 = 100
a = 102
b = 96
plt.figure(1, figsize=(20, 8));
for i in range(10):
plt.plot(np.array(simulate_one_path(x0, a, b)), lw=1.5, alpha=0.6)
plt.xlabel('Time', fontsize=20)
plt.ylabel('Stock Price', fontsize=20);
plt.xticks(fontsize=15)
plt.yticks(fontsize=15);
# plot horizontal lines
plt.axhline(y=a, color='r', linestyle='--')
plt.axhline(y=b, color='r', linestyle='--');
plt.title("Example Paths", fontsize=20);
plt.show()
# fix seed
np.random.seed(123)
# investigate convergence rate to converge to theoretical value
theo = (100-96)/(102-96)
def simulate_stock_price_stopping_fixed(x0, a, b, n):
# another function such that we only terminate simulation
# until we have gathered a fixed number of paths
i = 0
# number of paths stopped at a
count = 0
while count < n:
if count % 10000 == 0:
print("simulating stopping = {}".format(count))
x = simulate_one_path(x0, a, b)
if x[-1] == a:
count = count + 1
i = i + 1
return count / i
# parameters
x0 = 100
a = 102
b = 96
# test a range of n
n = np.array([10, 100, 500, 1000, 5000, 10000, 50000, 100000])
all_probs = []
for i in range(len(n)):
all_probs.append(simulate_stock_price_stopping_fixed(x0, a, b, n[i]))
simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 0 simulating stopping = 10000 simulating stopping = 20000 simulating stopping = 20000 simulating stopping = 30000 simulating stopping = 30000 simulating stopping = 40000 simulating stopping = 0 simulating stopping = 0 simulating stopping = 10000 simulating stopping = 20000 simulating stopping = 30000 simulating stopping = 30000 simulating stopping = 40000 simulating stopping = 40000 simulating stopping = 50000 simulating stopping = 60000 simulating stopping = 60000 simulating stopping = 70000 simulating stopping = 70000 simulating stopping = 80000 simulating stopping = 90000
# plot convergence
plt.figure(2, figsize=(20, 8))
plt.plot(n, all_probs, "-o", markersize=10, color="black", lw=1.5);
plt.axhline(y=theo, color='r', linestyle='--')
plt.xlabel('Number of Paths', fontsize=20)
plt.ylabel('Probability of Stopping', fontsize=20);
plt.xticks(fontsize=15)
plt.yticks(fontsize=15);
We visualize convergence rate.
import matplotlib.pyplot as plt
import numpy as np
# calculate error
errors = np.abs(np.array(all_probs) - theo)
# plot error with respect to number of paths on log scale
plt.figure(3, figsize=(20, 8))
plt.plot(n, errors, "-o", markersize=10, color="black", lw=1.5);
plt.axhline(y=1e-2, color='r', linestyle='--')
plt.xlabel('Number of Paths', fontsize=20)
plt.ylabel('Error', fontsize=20);
plt.xticks(fontsize=15)
plt.yticks(fontsize=15);
plt.yscale('log')
plt.xscale('log');
Question 2: Pricing of a Call Option¶
import numpy as np
from scipy.stats import norm
def black_scholes_call(S, K, T, r, sigma):
"""
Calculates the price of a European call option using the Black-Scholes formula.
Args:
S: Current stock price.
K: Strike price.
T: Time to expiration (in years).
r: Risk-free interest rate.
sigma: Implied volatility.
Returns:
The price of the European call option.
"""
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
# test call option price formula
S = 200
K = 180
T = 1/12
r = 0.02
sigma = 0.15
print("Price of call option = {}".format(black_scholes_call(S, K, T, r, sigma)))
Price of call option = 20.317700916589644
The Black-Scholes formula is used to calculate the theoretical price of a European call option. The formula is given by:
$$ C(S, t) = S_0 \cdot N(d_1) - K \cdot e^{-r(T - t)} \cdot N(d_2) $$
Where:
- $C(S, t)$ is the price of the call option.
- $S_0$ is the current stock price.
- $K$ is the strike price of the option.
- $T - t$ is the time to expiration of the option (in years).
- $r$ is the risk-free interest rate.
- $\sigma$ is the volatility of the stock's returns (standard deviation of returns).
- $N(\cdot)$ is the cumulative distribution function (CDF) of the standard normal distribution.
- $d_1$ and $d_2$ are calculated as follows:
$d_1 = \frac{\ln\left(\frac{S_0}{K}\right) + \left(r + \frac{\sigma^2}{2}\right)(T - t)}{\sigma \sqrt{T - t}}$
$d_2 = d_1 - \sigma \sqrt{T - t}$
Question 3: Pricing call option using Monte Carlo¶
def step(St, r, sigma, dt):
# take one step of the geometric brownian motion
return St * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * np.random.normal())
def get_path(S, r, sigma, tgrid):
# get one path of geometric brownian motion and save all points
St = S
path = [St]
for t in tgrid[1:]:
St = step(St, r, sigma, tgrid[1]-tgrid[0])
path.append(St)
return np.array(path)
# generate a few example paths
S = 200
K = 180
T = 1/12
r = 0.02
sigma = 0.15
# discretize time
tgrid = np.linspace(0, T, 240)
plt.figure(4, figsize=(20, 8))
for i in range(10):
# generate one path
path = get_path(S, r, sigma, tgrid)
plt.plot(tgrid, path, lw=1.5, alpha=0.6)
plt.xlabel('Time', fontsize=20)
plt.ylabel('Stock Price (GBM)', fontsize=20);
plt.xticks(fontsize=15)
plt.yticks(fontsize=15);
def monte_carlo_price(S, K, tgrid, r, sigma, n):
# approximate call option price using n paths of GBM
all_prices = []
for i in range(n):
if i % 1000 == 0:
print("simulating price = {}".format(i))
# generate one path
path = get_path(S, r, sigma, tgrid)
# calculate call option price using simulated path
option_price = np.maximum(path[-1] - K, 0)
# save
all_prices.append(option_price)
return np.mean(all_prices)
monte_carlo_price(S, K, tgrid, r, sigma, 100)
simulating price = 0
20.301260293551753
# fix seed
np.random.seed(123)
# test a range of n
n = np.array([10, 100, 500, 1000, 5000, 10000])
all_prices = []
for i in range(len(n)):
all_prices.append(monte_carlo_price(S, K, tgrid, r, sigma, n[i]))
# plot convergence
plt.figure(5, figsize=(20, 8))
plt.plot(n, all_prices, "-o", markersize=10, color="black", lw=1.5);
plt.axhline(y=black_scholes_call(S, K, T, r, sigma), color='r', linestyle='--')
plt.xlabel('Number of Paths', fontsize=20)
plt.ylabel('Call Option Price', fontsize=20);
plt.xticks(fontsize=15)
plt.yticks(fontsize=15);
simulating price = 0 simulating price = 0 simulating price = 0 simulating price = 0 simulating price = 0 simulating price = 1000 simulating price = 2000 simulating price = 3000 simulating price = 4000 simulating price = 0 simulating price = 1000 simulating price = 2000 simulating price = 3000 simulating price = 4000 simulating price = 5000 simulating price = 6000 simulating price = 7000 simulating price = 8000 simulating price = 9000
