Setup¶
import numpy as np
import matplotlib.pyplot as plt
import math
Explanation¶
A common numerical method to generate sample paths of a standard (mean 0, variance t) Brownian Motion is to use the $Euler-Maruyama$ approximation:
$ B(t + ∆t) = B(t) +\sqrt[]{∆t}N(0, 1) $
where $N(0, 1)$ is sample draw from a standard Normal distribution, mean 0, variance 1. Similarly, to generate a Brownian Motion with drift parameter μ and diffusion coefficient σ, we can adapt the above:
$B(t+∆t) = B(t)+μ∆t+σ ∆tN(0,1)$
In this case, $B(t) ∼ N(μt, σ2t)$.
Examples¶
(1) Generate 1,000 paths of a Brownian Motion for 1y with 12 monthly periods, ∆t = 1/12, r = 5%,μ = 4.5%,σ = 10%
Generate the 1y monthly evolution of a random asset starting with $A(0)=100$, and whose return is modeled by above BM, i.e., $A(t) = A(0)e^{B(t)}$
Price the 105-strike call option: $ D(1y)E[max(0, A(1y) − 105)]$, where $D(1y) = \frac{1}{(1 + 5\%)}$. Repeat this 10 times to observe run-to-run variability, i.e, simulation noise.
(2) Antithetic Variance Reduction: Generate 500 sample paths for the Brownian Motion as above, and for each sample path, create its mirror-image: $B(t, ω∗) = −B(t, ω)$ to come up with 1000 total sample paths. Repeat Steps 2, 3 above. Does the simulation noise improve?
(3) Add a 110 Knock-out feature to above option: The call option is knocked out (expires worth- less right away) if asset exceeds 110 before expiration (1y). Does this reduce the cost of the option?
1¶
## assigning values
n = 12
t = 1
mu = 0.045
sigma = 0.10
callprice = np.zeros(10)
df = 1/(1.05)
def generateIncrements(t,n,mu,sigma,numPaths):
dt=t/n
# Generate increments
randomSamples=np.random.normal(0,1,(numPaths,n))
X=mu*dt+sigma*math.sqrt(dt)*randomSamples
return(X)
def generateBM(X):
# Generate Brownian Motion
numPaths=len(X)
n=len(X[0])
BM=np.zeros((numPaths,n+1))
for t in range(1,n+1):
BM[:,t]=BM[:,t-1]+X[:,t-1]
return BM
def generateAssetPrices(A0,BM,plotFlag): # Generate Asset Prices
A=A0*np.exp(BM)
if plotFlag:
numPaths=len(A)
plt.figure(figsize=(20,10))
for i in range(numPaths):
plt.plot(A[i])
return(A)
def calculateCallPrice(A,K,df):
finalValue=A[:,len(A[0])-1]
C0=df*np.average(np.maximum(0,finalValue-K))
return(C0)
for i in range(10):
Increment = generateIncrements(t,n,mu,sigma,numPaths=1000)
Brownian = generateBM(Increment)
AssetPrices = generateAssetPrices(100, Brownian, False)
callprice[i] = calculateCallPrice(AssetPrices, 105, df)
print("The call option price: ", callprice[i])
print("The average call option price is ", np.mean(callprice))
print('The standard deviation of call option prices is ', np.sqrt(np.var(callprice)))
generateAssetPrices(100, Brownian, True)
Variance Reduction¶
def generateAntiTheticIncrements(t,n,mu,sigma,numPaths):
dt=t/n
# Generate half of the randomSamples
randomSamples1=np.random.normal(0,1,(int(numPaths/2),n))
# Generate the other half as "mirror" image (in this case reflected around 0)
randomSamples2=-1*randomSamples1
# Put them together
randomSamples=np.vstack((randomSamples1,randomSamples2))
# Generate Increments
X=mu*dt+sigma*math.sqrt(dt)*randomSamples
return(X)
2¶
for i in range(10):
Increment = generateAntiTheticIncrements(t,n,mu,sigma,numPaths=1000)
Brownian = generateBM(Increment)
AssetPrices = generateAssetPrices(100, Brownian, False)
callprice[i] = calculateCallPrice(AssetPrices, 105, df)
print("The call option price: ", callprice[i])
print("The average call option price is ", np.mean(callprice))
print('The standard deviation of call option prices is ', np.sqrt(np.var(callprice)))
AssetPrices = generateAssetPrices(100, Brownian, True)
The standard deviation is slightly lower; therefore, the simulation noise improves slightly¶
3¶
def calculateKnockCallPrice(A,K,df,KnockOut):
Alive = np.ones(len(A))
for i in range(1000):
for j in range(1,13):
if A[i,j]>=KnockOut:
Alive[i] = 0
#Multiply the matrix element-wise
finalValue=A[:,len(A[0])-1]*Alive
C0=df*np.average(np.maximum(0,finalValue-K))
return(C0)
Knockout effect on NORMAL increments¶
for i in range(10):
Increment = generateIncrements(t,n,mu,sigma,numPaths=1000)
Brownian = generateBM(Increment)
AssetPrices = generateAssetPrices(100, Brownian, False)
callprice[i] = calculateKnockCallPrice(AssetPrices, 105, df, 110)
print("The call option price: ", callprice[i])
print("The average call option price is ", np.mean(callprice))
print('The standard deviation of call option prices is ', np.sqrt(np.var(callprice)))
Knockout effect on ANTITHETIC increments¶
for i in range(10):
Increment = generateAntiTheticIncrements(t,n,mu,sigma,numPaths=1000)
Brownian = generateBM(Increment)
AssetPrices = generateAssetPrices(100, Brownian, False)
callprice[i] = calculateKnockCallPrice(AssetPrices, 105, df, 110)
print("The call option price: ", callprice[i])
print("The average call option price is ", np.mean(callprice))
print('The standard deviation of call option prices is ', np.sqrt(np.var(callprice)))