import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from numpy.random import standard_normal, seed
import scipy.stats as stats
from scipy.stats import norm
import sys
sys.path.append("..")
import grading
import datetime
import time
import bspline
import bspline.splinelab as splinelab
### ONLY FOR GRADING. DO NOT EDIT ###
submissions=dict()
assignment_key="J_L65CoiEeiwfQ53m1Mlug"
all_parts=["9jLRK","YoMns","Wc3NN","fcl3r"]
### ONLY FOR GRADING. DO NOT EDIT ###
COURSERA_TOKEN ='UYFOfJYBW8iN3RfJ'# the key provided to the Student under his/her email on submission page
COURSERA_EMAIL = 'bai_liping@outlook.com'# the email
# The Black-Scholes prices
def bs_put(t, S0, K, r, sigma, T):
d1 = (np.log(S0/K) + (r + 1/2 * sigma**2) * (T-t)) / sigma / np.sqrt(T-t)
d2 = (np.log(S0/K) + (r - 1/2 * sigma**2) * (T-t)) / sigma / np.sqrt(T-t)
price = K * np.exp(-r * (T-t)) * norm.cdf(-d2) - S0 * norm.cdf(-d1)
return price
def bs_call(t, S0, K, r, sigma, T):
d1 = (np.log(S0/K) + (r + 1/2 * sigma**2) * (T-t)) / sigma / np.sqrt(T-t)
d2 = (np.log(S0/K) + (r - 1/2 * sigma**2) * (T-t)) / sigma / np.sqrt(T-t)
price = S0 * norm.cdf(d1) - K * np.exp(-r * (T-t)) * norm.cdf(d2)
return price
def d1(S0, K, r, sigma, T):
return (np.log(S0/K) + (r + sigma**2 / 2) * T)/(sigma * np.sqrt(T))
def d2(S0, K, r, sigma, T):
return (np.log(S0 / K) + (r - sigma**2 / 2) * T) / (sigma * np.sqrt(T))
Simulate \(N_{MC}\) stock price sample paths with \(T\) steps by the classical Black-Sholes formula.
\[dS_t=\mu S_tdt+\sigma S_tdW_t\quad\quad S_{t+1}=S_te^{\left(\mu-\frac{1}{2}\sigma^2\right)\Delta t+\sigma\sqrt{\Delta t}Z}\]where \(Z\) is a standard normal random variable.
MC paths are simulated by GeneratePaths() method of DiscreteBlackScholes class.
Part 1
Class DiscreteBlackScholes implements the above calculations with class variables to math symbols mapping of:
\[\Delta S_t=S_{t+1} - e^{-r\Delta t} S_t\space \quad t=T-1,...,0\]Instructions: Some portions of code in DiscreteBlackScholes have bee taken out. You are to implement the missing portions of code in DiscreteBlackScholes class.
\[\Pi_t=e^{-r\Delta t}\left[\Pi_{t+1}-u_t \Delta S_t\right]\quad t=T-1,...,0\]-
implement DiscreteBlackScholes.function_A_vec() method \(A_{nm}^{\left(t\right)}=\sum_{k=1}^{N_{MC}}{\Phi_n\left(X_t^k\right)\Phi_m\left(X_t^k\right)\left(\Delta\hat{S}_t^k\right)^2}\quad\quad\)
- implement DiscreteBlackScholes.function_B_vec() method \(B_n^{\left(t\right)}=\sum_{k=1}^{N_{MC}}{\Phi_n\left(X_t^k\right)\left[\hat\Pi_{t+1}^k\Delta\hat{S}_t^k+\frac{1}{2\gamma\lambda}\Delta S_t^k\right]}\)
- implement DiscreteBlackScholes.gen_paths() method using the following relation: \(S_{t+1}=S_te^{\left(\mu-\frac{1}{2}\sigma^2\right)\Delta t+\sigma\sqrt{\Delta t}Z}\) where \(Z \sim N(0,1)\)
- implement parts of DiscreteBlackScholes.roll_backward()
- DiscreteBlackScholes.bVals corresponds to \(B_t\) and is computed as \(B_t = e^{-r\Delta t}\left[B_{t+1} + (u_{t+1} - u_t)S_{t+1}\right]\quad t=T-1,...,0\)
DiscreteBlackScholes.opt_hedge corresponds to \(\phi_t\) and is computed as \(\phi_t=\mathbf A_t^{-1}\mathbf B_t\)
class DiscreteBlackScholes:
"""
Class implementing discrete Black Scholes
DiscreteBlackScholes is class for pricing and hedging under
the real-world measure for a one-dimensional Black-Scholes setting
"""
def __init__(self,
s0,
strike,
vol,
T,
r,
mu,
numSteps,
numPaths):
"""
:param s0: initial price of the underlying
:param strike: option strike
:param vol: volatility
:param T: time to maturity, in years
:param r: risk-free rate,
:param mu: real drift, asset drift
:param numSteps: number of time steps
:param numPaths: number of Monte Carlo paths
"""
self.s0 = s0
self.strike = strike
self.vol = vol
self.T = T
self.r = r
self.mu = mu
self.numSteps = numSteps
self.numPaths = numPaths
self.dt = self.T / self.numSteps # time step
self.gamma = np.exp(-r * self.dt) # discount factor for one time step, i.e. gamma in the QLBS paper
self.sVals = np.zeros((self.numPaths, self.numSteps + 1), 'float') # matrix of stock values
# initialize half of the paths with stock price values ranging from 0.5 to 1.5 of s0
# the other half of the paths start with s0
half_paths = int(numPaths / 2)
if False:
# Grau (2010) "Applications of Least-Squares Regressions to Pricing and Hedging of Financial Derivatives"
self.sVals[:, 0] = (np.hstack((np.linspace(0.5 * s0, 1.5 * s0, half_paths),
s0 * np.ones(half_paths, 'float')))).T
self.sVals[:, 0] = s0 * np.ones(numPaths, 'float')
self.optionVals = np.zeros((self.numPaths, self.numSteps + 1), 'float') # matrix of option values
self.intrinsicVals = np.zeros((self.numPaths, self.numSteps + 1), 'float')
self.bVals = np.zeros((self.numPaths, self.numSteps + 1), 'float') # matrix of cash position values
self.opt_hedge = np.zeros((self.numPaths, self.numSteps + 1),
'float') # matrix of optimal hedges calculated from cross-sectional information F_t
self.X = None
self.data = None # matrix of features, i.e. self.X as sum of basis functions
self.delta_S_hat = None
# coef = 1.0/(2 * gamma * risk_lambda)
# override it by zero to have pure risk hedge
self.coef = 0.
def gen_paths(self):
"""
A simplest path generator
"""
np.random.seed(42)
# Spline basis of order p on knots k
### START CODE HERE ### (≈ 3-4 lines of code)
# self.sVals = your code goes here ...
# for-loop or while loop is allowed heres
Z = np.random.normal(0, 1, size=(self.numSteps + 1, self.numPaths)).T
for t in range(0, self.numSteps):
self.sVals[:, t + 1] = self.sVals[:, t] * np.exp((self.mu - 0.5 * self.vol**2) * self.dt + (self.vol * np.sqrt(self.dt) * Z[:, t + 1]))
print(self.sVals)
### END CODE HERE ###
# like in QLBS
delta_S = self.sVals[:, 1:] - np.exp(self.r * self.dt) * self.sVals[:, :self.numSteps]
self.delta_S_hat = np.apply_along_axis(lambda x: x - np.mean(x), axis=0, arr=delta_S)
# state variable
# delta_t here is due to their conventions
self.X = - (self.mu - 0.5 * self.vol ** 2) * np.arange(self.numSteps + 1) * self.dt + np.log(self.sVals)
X_min = np.min(np.min(self.X))
X_max = np.max(np.max(self.X))
print('X.shape = ', self.X.shape)
print('X_min, X_max = ', X_min, X_max)
p = 4 # order of spline (as-is; 3 = cubic, 4: B-spline?)
ncolloc = 12
tau = np.linspace(X_min, X_max, ncolloc) # These are the sites to which we would like to interpolate
# k is a knot vector that adds endpoints repeats as appropriate for a spline of order p
# To get meaningful results, one should have ncolloc >= p+1
k = splinelab.aptknt(tau, p)
basis = bspline.Bspline(k, p)
num_basis = ncolloc # len(k) #
self.data = np.zeros((self.numSteps + 1, self.numPaths, num_basis))
print('num_basis = ', num_basis)
print('dim self.data = ', self.data.shape)
# fill it, expand function in finite dimensional space
# in neural network the basis is the neural network itself
t_0 = time.time()
for ix in np.arange(self.numSteps + 1):
x = self.X[:, ix]
self.data[ix, :, :] = np.array([basis(el) for el in x])
t_end = time.time()
print('\nTime Cost of basis expansion:', t_end - t_0, 'seconds')
def function_A_vec(self, t, reg_param=1e-3):
"""
function_A_vec - compute the matrix A_{nm} from Eq. (52) (with a regularization!)
Eq. (52) in QLBS Q-Learner in the Black-Scholes-Merton article
Arguments:
t - time index, a scalar, an index into time axis of data_mat
reg_param - a scalar, regularization parameter
Return:
- np.array, i.e. matrix A_{nm} of dimension num_basis x num_basis
"""
X_mat = self.data[t, :, :]
num_basis_funcs = X_mat.shape[1]
this_dS = self.delta_S_hat[:, t]
hat_dS2 = (this_dS ** 2).reshape(-1, 1)
A_mat = np.dot(X_mat.T, X_mat * hat_dS2) + reg_param * np.eye(num_basis_funcs)
return A_mat
def function_B_vec(self, t, Pi_hat):
"""
function_B_vec - compute vector B_{n} from Eq. (52) QLBS Q-Learner in the Black-Scholes-Merton article
Arguments:
t - time index, a scalar, an index into time axis of delta_S_hat
Pi_hat - pandas.DataFrame of dimension N_MC x T of portfolio values
Return:
B_vec - np.array() of dimension num_basis x 1
"""
tmp = Pi_hat * self.delta_S_hat[:, t] + self.coef * (np.exp((self.mu - self.r) * self.dt)) * self.sVals[:, t]
X_mat = self.data[t, :, :] # matrix of dimension N_MC x num_basis
B_vec = np.dot(X_mat.T, tmp)
return B_vec
def seed_intrinsic(self, strike=None, cp='P'):
"""
initilaize option value and intrinsic value for each node
"""
if strike is not None:
self.strike = strike
if cp == 'P':
# payoff function at maturity T: max(K - S(T),0) for all paths
self.optionVals = np.maximum(self.strike - self.sVals[:, -1], 0).copy()
# payoff function for all paths, at all time slices
self.intrinsicVals = np.maximum(self.strike - self.sVals, 0).copy()
elif cp == 'C':
# payoff function at maturity T: max(S(T) -K,0) for all paths
self.optionVals = np.maximum(self.sVals[:, -1] - self.strike, 0).copy()
# payoff function for all paths, at all time slices
self.intrinsicVals = np.maximum(self.sVals - self.strike, 0).copy()
else:
raise Exception('Invalid parameter: %s'% cp)
self.bVals[:, -1] = self.intrinsicVals[:, -1]
def roll_backward(self):
"""
Roll the price and optimal hedge back in time starting from maturity
"""
for t in range(self.numSteps - 1, -1, -1):
# determine the expected portfolio value at the next time node
piNext = self.bVals[:, t+1] + self.opt_hedge[:, t+1] * self.sVals[:, t+1]
pi_hat = piNext - np.mean(piNext)
A_mat = self.function_A_vec(t)
B_vec = self.function_B_vec(t, pi_hat)
phi = np.dot(np.linalg.inv(A_mat), B_vec)
self.opt_hedge[:, t] = np.dot(self.data[t, :, :], phi)
### START CODE HERE ### (≈ 1-2 lines of code)
# implement code to update self.bVals
# self.bVals[:,t] = your code goes here ....
self.bVals[:,t] = np.exp(-self.r * self.dt) * (self.bVals[:, t+1] + (self.opt_hedge[:, t+1] - self.opt_hedge[:, t]) * self.sVals[:, t+1])
### END CODE HERE ###
# calculate the initial portfolio value
initPortfolioVal = self.bVals[:, 0] + self.opt_hedge[:, 0] * self.sVals[:, 0]
# use only the second half of the paths generated with paths starting from S0
optionVal = np.mean(initPortfolioVal)
optionValVar = np.std(initPortfolioVal)
delta = np.mean(self.opt_hedge[:, 0])
return optionVal, delta, optionValVar
np.random.seed(42)
strike_k = 95
test_vol = 0.2
test_mu = 0.03
dt = 0.01
rfr = 0.05
num_paths = 100
num_periods = 252
hMC = DiscreteBlackScholes(100, strike_k, test_vol, 1., rfr, test_mu, num_periods, num_paths)
hMC.gen_paths()
t = hMC.numSteps - 1
piNext = hMC.bVals[:, t+1] + 0.1 * hMC.sVals[:, t+1]
pi_hat = piNext - np.mean(piNext)
A_mat = hMC.function_A_vec(t)
B_vec = hMC.function_B_vec(t, pi_hat)
phi = np.dot(np.linalg.inv(A_mat), B_vec)
opt_hedge = np.dot(hMC.data[t, :, :], phi)
# plot the results
fig = plt.figure(figsize=(12,4))
ax1 = fig.add_subplot(121)
ax1.scatter(hMC.sVals[:,t], pi_hat)
ax1.set_title(r'Expected $\Pi_0$ vs. $S_t$')
ax1.set_xlabel(r'$S_t$')
ax1.set_ylabel(r'$\Pi_0$')
[[ 100. 98.23650359 98.6842395 ..., 111.52820537
111.93345414 111.50088104]
[ 100. 99.47538589 100.18466561 ..., 69.58859259
69.36721589 68.46903615]
[ 100. 99.57310236 100.94511135 ..., 110.66761375
110.53260244 110.37282496]
...,
[ 100. 100.19783913 100.59050962 ..., 151.7887043 151.63565543
152.14692905]
[ 100. 100.07733423 101.11151453 ..., 103.08321744
101.41095506 101.46651123]
[ 100. 98.57422289 99.36322314 ..., 91.31429149
91.06798685 92.50219743]]
X.shape = (100, 253)
X_min, X_max = 4.10743882917 5.16553756345
num_basis = 12
dim self.data = (253, 100, 12)
Time Cost of basis expansion: 4.927450180053711 seconds
<matplotlib.text.Text at 0x7f91ad760f98>
### GRADED PART (DO NOT EDIT) ###
part_1 = list(pi_hat)
try:
part1 = " ".join(map(repr, part_1))
except TypeError:
part1 = repr(part_1)
submissions[all_parts[0]]=part1
grading.submit(COURSERA_EMAIL, COURSERA_TOKEN, assignment_key,all_parts[:1],all_parts,submissions)
pi_hat
### GRADED PART (DO NOT EDIT) ###
Submission successful, please check on the coursera grader page for the status
array([ 0.81274895, -3.49043554, 0.69994334, 1.61239986, -0.25153316,
-3.19082265, 0.8848621 , -2.0380868 , 0.45033564, 3.74872863,
-0.6568227 , 1.74148929, 0.94314331, -4.19716113, 1.72135256,
-0.66188482, 6.95675041, -2.20512677, -0.14942482, 0.30067272,
3.33419402, 0.68536713, 1.65097153, 2.69898611, 1.22528159,
1.47188744, -2.48129898, -0.37360224, 0.81064666, -1.05269459,
0.02476551, -1.88267258, 0.11748169, -0.9038195 , 0.69753811,
-0.54805029, 1.97594593, -0.44331403, 0.62134931, -1.86191032,
-3.21226413, 2.24508097, -2.23451292, -0.13488281, 3.64364848,
-0.11270281, -1.15582237, -3.30169455, 1.74454841, -1.10425448,
2.10192819, 1.80570507, -1.68587001, -1.42113397, -2.70292006,
0.79454199, -2.05396827, 3.13973887, -1.08786662, 0.42347686,
1.32787012, 0.55924965, -3.54140814, -3.70258632, 2.14853641,
1.11495458, 3.69639676, 0.62864736, -2.62282995, -0.05315552,
1.05789698, 1.8023196 , -3.35217374, -2.30436466, -2.68609519,
0.95284884, -1.35963013, -0.56273408, -0.08311276, 0.79044269,
0.46247485, -1.04921463, -2.18122285, 1.82920128, 1.05635272,
0.90161346, -1.93870347, -0.37549305, -1.96383274, 1.9772888 ,
-1.37386984, 0.95230068, 0.88842589, -1.42214528, -2.60256696,
-1.53509699, 4.47491253, 4.87735375, -0.19068803, -1.08711941])
# input parameters
s0 = 100.0
strike = 100.0
r = 0.05
mu = 0.07 # 0.05
vol = 0.4
T = 1.0
# Simulation Parameters
numPaths = 50000 # number of Monte Carlo trials
numSteps = 6
# create the class object
hMC = DiscreteBlackScholes(s0, strike, vol, T, r, mu, numSteps, numPaths)
# calculation
hMC.gen_paths()
hMC.seed_intrinsic()
option_val, delta, option_val_variance = hMC.roll_backward()
bs_call_value = bs_put(0, s0, K=strike, r=r, sigma=vol, T=T)
print('Option value = ', option_val)
print('Option value variance = ', option_val_variance)
print('Option delta = ', delta)
print('BS value', bs_call_value)
[[ 100. 101.44740793 119.84140463 ..., 192.78653975 210.7076386
167.37134738]
[ 100. 98.79378416 81.67103247 ..., 78.75163254
104.69106128 114.29766651]
[ 100. 116.62110943 127.89787986 ..., 85.9631909 79.72061217
78.03372489]
...,
[ 100. 106.73222875 103.49782882 ..., 108.30352919
96.76512324 114.08668191]
[ 100. 96.45073828 98.70345177 ..., 89.5899346 75.07626471
91.91332688]
[ 100. 101.81014094 115.21893111 ..., 68.72837469
64.71929858 65.04500528]]
X.shape = (50000, 7)
X_min, X_max = 2.96880459823 6.37164911461
num_basis = 12
dim self.data = (7, 50000, 12)
Time Cost of basis expansion: 59.888269901275635 seconds
Option value = 13.1083498881
Option value variance = 5.17079676287
Option delta = -0.356133671254
BS value 13.1458939003
### GRADED PART (DO NOT EDIT) ###
part2 = str(option_val)
submissions[all_parts[1]]=part2
grading.submit(COURSERA_EMAIL, COURSERA_TOKEN, assignment_key,all_parts[:2],all_parts,submissions)
option_val
### GRADED PART (DO NOT EDIT) ###
Submission successful, please check on the coursera grader page for the status
13.108349888127721
strikes = np.linspace(85, 110, 6)
results = [None] * len(strikes)
bs_prices = np.zeros(len(strikes))
bs_deltas = np.zeros(len(strikes))
numPaths = 50000
hMC = DiscreteBlackScholes(s0, strike, vol, T, r, mu, numSteps, numPaths)
hMC.gen_paths()
for ix, k_strike in enumerate(strikes):
hMC.seed_intrinsic(k_strike)
results[ix] = hMC.roll_backward()
bs_prices[ix] = bs_put(0, s0, K=k_strike, r=r, sigma=vol, T=T)
bs_deltas[ix] = norm.cdf(d1(s0, K=k_strike, r=r, sigma=vol, T=T)) - 1
bs_prices
[[ 100. 101.44740793 119.84140463 ..., 192.78653975 210.7076386
167.37134738]
[ 100. 98.79378416 81.67103247 ..., 78.75163254
104.69106128 114.29766651]
[ 100. 116.62110943 127.89787986 ..., 85.9631909 79.72061217
78.03372489]
...,
[ 100. 106.73222875 103.49782882 ..., 108.30352919
96.76512324 114.08668191]
[ 100. 96.45073828 98.70345177 ..., 89.5899346 75.07626471
91.91332688]
[ 100. 101.81014094 115.21893111 ..., 68.72837469
64.71929858 65.04500528]]
X.shape = (50000, 7)
X_min, X_max = 2.96880459823 6.37164911461
num_basis = 12
dim self.data = (7, 50000, 12)
Time Cost of basis expansion: 60.33323049545288 seconds
array([ 6.70326307, 8.59543726, 10.74614496, 13.1458939 ,
15.78197485, 18.63949388])
mc_prices = np.array([x[0] for x in results])
mc_deltas = np.array([x[1] for x in results])
price_variances = np.array([x[-1] for x in results])
prices_diff = mc_prices - bs_prices
deltas_diff = mc_deltas - bs_deltas
# price_variances
### GRADED PART (DO NOT EDIT) ###
part_3 = list(prices_diff)
try:
part3 = " ".join(map(repr, part_3))
except TypeError:
part3 = repr(part_3)
submissions[all_parts[2]]=part3
grading.submit(COURSERA_EMAIL, COURSERA_TOKEN, assignment_key,all_parts[:3],all_parts,submissions)
prices_diff
### GRADED PART (DO NOT EDIT) ###
Submission successful, please check on the coursera grader page for the status
array([-0.03641516, -0.0403414 , -0.039966 , -0.03754401, -0.03240018,
-0.02997068])
### GRADED PART (DO NOT EDIT) ###
part_4 = list(deltas_diff)
try:
part4 = " ".join(map(repr, part_4))
except TypeError:
part4= repr(part_4)
submissions[all_parts[3]]=part4
grading.submit(COURSERA_EMAIL, COURSERA_TOKEN, assignment_key,all_parts[:4],all_parts,submissions)
deltas_diff
### GRADED PART (DO NOT EDIT) ###
Submission successful, please check on the coursera grader page for the status
array([ 0.01279812, 0.01416024, 0.01532707, 0.01645686, 0.01715368,
0.01780667])