In this tutorial we will investigate ways we can reduce the variance of results from a Monte Carlo simulation method when valuing financial derivatives. The mathematic notation and examples are from Les Clewlow and Chris Strickland’s book Implementing Derivatives Models.
Unfortunately, although a great method for approximating option values with complex payoffs or high dimensionality, in order to get an acceptably accurate estimate we must perform a large number of simulations M. Instead we can lean on Variance Reduction methods which work on the same principles as that of hedging an option position. i.e. the variability of a hedged option portfolio will have a smaller variance that that of it’s unhedged counterpart.
# Import dependencies import time import math import datetime import numpy as np import pandas as pd import scipy.stats as stats from py_vollib.black_scholes import black_scholes as bs
General Control Variate Equation
For J control variates we have:
\(\Large C_0\exp(rT) = C_T – \sum^J_{i=j}\beta_j cv_j + \eta\)
where
– \(\beta_j\) are factors to account for the “true” linear relationship between the option pay-off and the control variate \(cv_j\)
– \(\eta\) accounts for errors: – discrete rebalancing – approximations in hedge sensitivities (calc. delta / gamma) Option price as the sum of the linear relationships with J control variates \( \large C_T =\beta_0 + \sum^J_{i=j}\beta_j cv_j + \eta\)
where \(\beta_0 = C_0\exp(rT)\) is the forward price of the option
If we perform M simulations at discrete time intervals N we can regard the pay-offs and control variates as samples of the linear relationship with some noise. We can estimate the true relationship between control variates and option pay-offs with least-squares regression:
\(\beta = (X’X)^{-1}X’Y\)
We don’t want biased estimates of \(\beta_j\) so these should be precomputed by least-squares regression to establish the relationship between types of control variates and options first. These estaimates of \(\beta_j\) values can then be used for \(cv_j\) for pricing any option.
# initial derivative parameters S = 101.15 #stock price K = 98.01 #strike price vol = 0.0991 #volatility (%) r = 0.015 #risk-free rate (%) N = 20 #number of time steps M = 1000 #number of simulations market_value = 3.86 #market price of option T = ((datetime.date(2022,3,17)-datetime.date(2022,1,17)).days+1)/365 #time in years print(T)
Sample Estimate – no variance reduction methods
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
lnS = np.log(S)
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_lnSt = nudt + volsdt*Z
lnSt = lnS + np.cumsum(delta_lnSt, axis=0)
lnSt = np.concatenate( (np.full(shape=(1, M), fill_value=lnS), lnSt ) )
# Compute Expectation and SE
ST = np.exp(lnSt)
CT = np.maximum(0, ST - K)
C0_se = np.exp(-r*T)*np.sum(CT[-1])/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT[-1] - C0_se)**2) / (M-1) )
SE_se = sigma/np.sqrt(M)
mc_time_se = time.time() - start_time
print("Sample Estimate: Call value is ${0} with SE +/- {1}".format(np.round(C0_se,2),np.round(SE_se,2)))
print("Computation time is: ", round(mc_time_se,4))
Implementation of Antithetic Variate
To implement an antithetic variate we create a hypothetical asset which is perfectly negatively correlated with the original asset. Implementation is very simple, and if we consider the example of the European Call Option (as in last weeks video). Our simulated pay-offs are under the following \(S_t\) dynamics:
\(\large S_{t+\Delta t} = S_{t} \exp( \nu \Delta t + \sigma (z_{t+\Delta t}- z_t) )\)
Where \((z_{t+\Delta t}- z_t) \sim N(0,\Delta t) \sim \sqrt{\Delta t} N(0,1) \sim \sqrt{\Delta t} \epsilon_i\)
Contract Simulation
– \(\large C_{T,i} = max(0, S \exp( \nu \Delta T + \sigma \sqrt{T} (\epsilon_i) ) – K)\)
– \(\large \bar{C}_{T,i} = max(0, S \exp( \nu \Delta T + \sigma \sqrt{T} (-\epsilon_i) ) – K)\)
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
lnS = np.log(S)
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_lnSt1 = nudt + volsdt*Z
delta_lnSt2 = nudt - volsdt*Z
lnSt1 = lnS + np.cumsum(delta_lnSt1, axis=0)
lnSt2 = lnS + np.cumsum(delta_lnSt2, axis=0)
# Compute Expectation and SE
ST1 = np.exp(lnSt1)
ST2 = np.exp(lnSt2)
CT = 0.5 * ( np.maximum(0, ST1[-1] - K) + np.maximum(0, ST2[-1] - K) )
C0_av = np.exp(-r*T)*np.sum(CT)/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0_av)**2) / (M-1) )
SE_av = sigma/np.sqrt(M)
mc_time_av = time.time() - start_time
print("Call value is ${0} with SE +/- {1}".format(np.round(C0_av,2),np.round(SE_av,2)))
print("Computation time is: ", round(mc_time_av,4))
Implementation of Delta-based Control Variates
\(\large cv_1 = \sum^{N-1}_{i=0} \frac{\delta C_{t_i}}{\delta S}(S_{t_{i+1}} – {\mathbb E}[S_{t_i}])\exp{r(T-t_{i+1})}\)
\(\large C_{t_0}\exp{rT} = C_T + \beta_1 cv_1 + \eta\)
where with GBM dynamics:
– \({\mathbb E}[S_{t_{i+1}}] = S_{t_{i}} \exp (r \Delta t_i)\)
– \(\beta_1 = -1\) which is the appropriate value where we have exact delta for European Option
def delta_calc(r, S, K, T, sigma, type="c"):
"Calculate delta of an option"
d1 = (np.log(S/K) + (r + sigma**2/2)*T)/(sigma*np.sqrt(T))
try:
if type == "c":
delta_calc = stats.norm.cdf(d1, 0, 1)
elif type == "p":
delta_calc = -stats.norm.cdf(-d1, 0, 1)
return delta_calc
except:
print("Please confirm option type, either 'c' for Call or 'p' for Put!")
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
erdt = np.exp(r*dt)
cv = 0
beta1 = -1
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_St = nudt + volsdt*Z
ST = S*np.cumprod( np.exp(delta_St), axis=0)
ST = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST ) )
deltaSt = delta_calc(r, ST[:-1].T, K, np.linspace(T,dt,N), vol, "c").T
cv = np.cumsum(deltaSt*(ST[1:] - ST[:-1]*erdt), axis=0)
CT = np.maximum(0, ST[-1] - K) + beta1*cv[-1]
C0_dv = np.exp(-r*T)*np.sum(CT)/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0_dv)**2) / (M-1) )
sigma = np.std(np.exp(-r*T)*CT)
SE_dv = sigma/np.sqrt(M)
mc_time_dv = time.time() - start_time
print("Call value is ${0} with SE +/- {1}".format(np.round(C0_dv,2),np.round(SE_dv,3)))
print("Computation time is: ", round(mc_time_dv,4))
Gamma Based Control Variate
The control variate is a random variable whose expected value we know, which is correlated with the varaible we are trying to estimate.
In the same way as for \(cv_1\) we can create other control variates, which are equivalent to other hedges.
For example a gamma-based control variate (\(cv_2\)):
\(\large cv_2 = \sum^{N-1}_{i=0} \frac{\delta^2 C_{t_i}}{\delta S^2}((\Delta S_{t_{i+1}})^2 – {\mathbb E}[(\Delta S_{t_i})^2])\exp{r(T-t_{i+1})}\)
Where \({\mathbb E}[(\Delta S_{t_i})^2] = S_{t_i}^2 (\exp([2r+\sigma^2]\Delta t_i)-2\exp(r\Delta t_i)+1)\)
def gamma_calc(r, S, K, T, sigma):
"Calculate delta of an option"
d1 = (np.log(S/K) + (r + sigma**2/2)*T)/(sigma*np.sqrt(T))
try:
gamma_calc = stats.norm.pdf(d1, 0, 1)/(S*sigma*np.sqrt(T))
return gamma_calc
except:
print("Please confirm option type, either 'c' for Call or 'p' for Put!")
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
erdt = np.exp(r*dt)
ergamma = np.exp((2*r+vol**2)*dt) - 2*erdt + 1
beta2 = -0.5
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_St = nudt + volsdt*Z
ST = S*np.cumprod( np.exp(delta_St), axis=0)
ST = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST ) )
gammaSt = gamma_calc(r, ST[:-1].T, K, np.linspace(T,dt,N), vol).T
cv2 = np.cumsum(gammaSt*((ST[1:] - ST[:-1])**2 - ergamma*ST[:-1]**2), axis=0)
CT = np.maximum(0, ST[-1] - K) + beta2*cv2[-1]
C0_gv = np.exp(-r*T)*np.sum(CT)/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0_gv)**2) / (M-1) )
sigma = np.std(np.exp(-r*T)*CT)
SE_gv = sigma/np.sqrt(M)
mc_time_gv = time.time() - start_time
print("Call value is ${0} with SE +/- {1}".format(np.round(C0_gv,2),np.round(SE_gv,3)))
print("Computation time is: ", round(mc_time_gv,4))
Combined Antithetic and Delta Variates
\(C_T = 0.5( max(0, S_{1,t} – K) + max(0, S_{2,t} – K) + \beta_1 cv_1)\)
where \(cv_1\) is delta variate but now we have to account for antithetic variates – two perfectly negatively correlated underlyings.
\(cv_1 = 0.5 * \beta_1 * (cv_{11} + cv_{12})\)
where:
– \(cv_{11} = \Delta_{S_{1,t}}[S_{1,t_{i+1}} – S_{1,t_{i}} \exp (r \Delta t_i)]\)
– \(cv_{12} = \Delta_{S_{2,t}}[S_{2,t_{i+1}} – S_{2,t_{i}} \exp (r \Delta t_i)]\)
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
erdt = np.exp(r*dt)
beta1 = -1
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_St1 = nudt + volsdt*Z
delta_St2 = nudt - volsdt*Z
ST1 = S*np.cumprod( np.exp(delta_St1), axis=0)
ST2 = S*np.cumprod( np.exp(delta_St2), axis=0)
ST1 = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST1 ) )
ST2 = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST2 ) )
# Calculate delta for both sets of underlying stock prices
deltaSt1 = delta_calc(r, ST1[:-1].T, K, np.linspace(T,dt,N), vol, "c").T
deltaSt2 = delta_calc(r, ST2[:-1].T, K, np.linspace(T,dt,N), vol, "c").T
# Calculate two sets of delta control variates for negatively correlated assets
cv11 = np.cumsum(deltaSt1*(ST1[1:] - ST1[:-1]*erdt), axis=0)
cv12 = np.cumsum(deltaSt2*(ST2[1:] - ST2[:-1]*erdt), axis=0)
CT = 0.5 * ( np.maximum(0, ST1[-1] - K) + beta1*cv11[-1]
+ np.maximum(0, ST2[-1] - K) + beta1*cv12[-1] )
C0_adv = np.exp(-r*T)*np.sum(CT)/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0_adv)**2) / (M-1) )
sigma = np.std(np.exp(-r*T)*CT)
SE_adv = sigma/np.sqrt(M)
mc_time_adv = time.time() - start_time
print("Call value is ${0} with SE +/- {1}".format(np.round(C0_adv,2),np.round(SE_adv,3)))
print("Computation time is: ", round(mc_time_adv,4))
Combined Antithetic, Delta and Gamma Variates
\(C_T = 0.5( max(0, S_{1,t} – K) + max(0, S_{2,t} – K) + \beta_1 cv_1 + \beta_2 cv_2)\)
where \(cv_1\) is delta variate and \(cv_2\) is the gamma variate. When combined with antithetic technique you have to apply the following!
\(cv_1 = 0.5 * \beta_1 * (cv_{11} + cv_{12})\)
– \(cv_{11} = \Delta_{S_{1,t}}[S_{1,t_{i+1}} – S_{1,t_{i}} \exp (r \Delta t_i)]\)
– \(cv_{12} = \Delta_{S_{2,t}}[S_{2,t_{i+1}} – S_{2,t_{i}} \exp (r \Delta t_i)]\)
\(cv_2 = 0.5 * \beta_2 * (cv_{21} + cv_{22})\)
– \(cv_{21} = \gamma_{S_{1,t}}[(S_{1,t_{i+1}} – S_{1,t_i})^2 – S_{1,t_i}^2 (\exp([2r+\sigma^2]\Delta t_i)-2\exp(r\Delta t_i)+1)]\)
– \(cv_{22} = \gamma_{S_{2,t}}[(S_{2,t_{i+1}} – S_{2,t_i})^2 – S_{2,t_i}^2 (\exp([2r+\sigma^2]\Delta t_i)-2\exp(r\Delta t_i)+1)]\)
start_time = time.time()
#precompute constants
dt = T/N
nudt = (r - 0.5*vol**2)*dt
volsdt = vol*np.sqrt(dt)
erdt = np.exp(r*dt)
ergamma = np.exp((2*r+vol**2)*dt) - 2*erdt + 1
beta1 = -1
beta2 = -0.5
# Monte Carlo Method
Z = np.random.normal(size=(N, M))
delta_St1 = nudt + volsdt*Z
delta_St2 = nudt - volsdt*Z
ST1 = S*np.cumprod( np.exp(delta_St1), axis=0)
ST2 = S*np.cumprod( np.exp(delta_St2), axis=0)
ST1 = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST1 ) )
ST2 = np.concatenate( (np.full(shape=(1, M), fill_value=S), ST2 ) )
# Calculate delta for both sets of underlying stock prices
deltaSt1 = delta_calc(r, ST1[:-1].T, K, np.linspace(T,dt,N), vol, "c").T
deltaSt2 = delta_calc(r, ST2[:-1].T, K, np.linspace(T,dt,N), vol, "c").T
# Calculate gamma for both sets of underlying stock prices
gammaSt1 = gamma_calc(r, ST1[:-1].T, K, np.linspace(T,dt,N), vol).T
gammaSt2 = gamma_calc(r, ST2[:-1].T, K, np.linspace(T,dt,N), vol).T
# Calculate two sets of delta control variates for negatively correlated assets
cv11 = np.cumsum(deltaSt1*(ST1[1:] - ST1[:-1]*erdt), axis=0)
cv12 = np.cumsum(deltaSt2*(ST2[1:] - ST2[:-1]*erdt), axis=0)
# Calculate two sets of gamma control variates for negatively correlated assets
cv21 = np.cumsum(gammaSt1*((ST1[1:] - ST1[:-1])**2 - ergamma*ST1[:-1]**2), axis=0)
cv22 = np.cumsum(gammaSt2*((ST2[1:] - ST2[:-1])**2 - ergamma*ST2[:-1]**2), axis=0)
CT = 0.5 * ( np.maximum(0, ST1[-1] - K) + beta1*cv11[-1] + beta2*cv21[-1]
+ np.maximum(0, ST2[-1] - K) + beta1*cv12[-1] + beta2*cv22[-1])
C0_adgv = np.exp(-r*T)*np.sum(CT)/M
sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0_adgv)**2) / (M-1) )
sigma = np.std(np.exp(-r*T)*CT)
SE_adgv = sigma/np.sqrt(M)
mc_time_adgv = time.time() - start_time
print("Call value is ${0} with SE +/- {1}".format(np.round(C0_adgv,2),np.round(SE_adgv,3)))
print("Computation time is: ", round(mc_time_adgv,4))
Comparing Reduction Methods
params = [K,round(T,2),S,vol,r,N,M,round(bs('c', S, K, T, r, vol),2), market_value]
params_rd = [round(param,2) for param in params]
data = {'Contract Terms':['Strike Price', 'Time to Maturity', 'Asset Price', 'Volatility',
'Riskfree Rate', 'Number of Time Steps', 'Number of Simuations',
'Standard European Call Price', 'Market Price'],
'Parameters': params_rd}
# Creates pandas DataFrame.
df = pd.DataFrame(data) #, index =['position1', 'position2', 'position3', 'position4'])
df
Trade-off between Error vs Computation Time
Comparison table of standard errors and relative computation time for each reduction method or combination
se_variates = [SE_se, SE_av, SE_dv, SE_gv, SE_adv, SE_adgv]
se_rd = [round(se,4) for se in se_variates]
se_red = [round(SE_se/se,2) for se in se_variates]
comp_time = [mc_time_se, mc_time_av, mc_time_dv, mc_time_gv, mc_time_adv, mc_time_adgv]
rel_time = [round(mc_time/mc_time_se,2) for mc_time in comp_time]
data = {'Standard Error (SE)': se_rd,
'SE Reduction Multiple': se_red,
'Relative Computation Time': rel_time}
# Creates pandas DataFrame.
df = pd.DataFrame(data, index =['Simple estimate', 'with antithetic variate',
'with delta-based control variate', 'with gamma-based control variate', 'with antithetic and delta variates', 'with all combined variates'])
df