Black Scholes option pricing & GPs

• The purpose of this notebook is to demonstrate the derivation of the greeks in a Gaussian Process Regression model (GP), fitted to option price data.

• In this notebook, European option prices are generated from the Black-Scholes model.

• The notebook begins by building a GP call model, where the input is the underlying price.

• The delta is then derived and compared with the Black-Scholes (BS) delta.

• The exercise is repeated, but using the volatility as the input instead of the underlying price.

• The vega of the GP is then derived and compared with the BS vega (vega measures sensitivity to volatility. Vega is the derivative of the option value with respect to the volatility of the underlying asset).

from BlackScholes import bsformula

import numpy as np
import scipy as sp

from sklearn import gaussian_process
from sklearn.gaussian_process.kernels import RBF

import matplotlib.pyplot as plt

%matplotlib inline

Black-Scholes Model

First, set the model parameters

KC = 130    # Call strike
KP = 70     # Put strike
r = 0.002   # risk-free rate
sigma = 0.4 # implied volatility
T = 2.0     # Time to maturity
S0 = 100    # Underlying spot

lb = 0   # lower bound on domain
ub = 300 # upper bound on domain
training_number = 100 # Number of training samples
testing_number = 50   # Number of testing samples
sigma_n = 1e-8 # additive noise in GP

Define the call and put prices using the BS model

call = lambda x, y: bsformula(1, lb+(ub-lb)*x, KC, r, T, y, 0)[0]
put = lambda x, y: bsformula(-1, lb+(ub-lb)*x, KP, r, T, y, 0)[0]

Delta

Generate the training and testing data, where the input is the gridded underlying and the output are the option prices.

x_train = np.array(np.linspace(0.01, 1.2, training_number), dtype='float32').reshape(training_number, 1)
x_test = np.array(np.linspace(0.01, 1.0, testing_number), dtype='float32').reshape(testing_number, 1)

y_train = []
    
for idx in range(len(x_train)):
    y_train.append(call(x_train[idx], sigma))
y_train = np.array(y_train)

Fit the GP model to the generated data

sk_kernel = RBF(length_scale=1.0, length_scale_bounds=(0.01, 10000.0))
gp = gaussian_process.GaussianProcessRegressor(kernel=sk_kernel, n_restarts_optimizer=20)
gp.fit(x_train, y_train)

/usr/local/lib/python3.6/dist-packages/sklearn/gaussian_process/_gpr.py:504: ConvergenceWarning: lbfgs failed to converge (status=2):
ABNORMAL_TERMINATION_IN_LNSRCH.

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
  _check_optimize_result("lbfgs", opt_res)

GaussianProcessRegressor(kernel=RBF(length_scale=1), n_restarts_optimizer=20)

Get the model's predicted outputs for each of the test inputs

y_pred, sigma_hat = gp.predict(x_test, return_std=True)

/usr/local/lib/python3.6/dist-packages/sklearn/gaussian_process/_gpr.py:370: UserWarning: Predicted variances smaller than 0. Setting those variances to 0.
  warnings.warn("Predicted variances smaller than 0. "

Derive the GP delta

l = gp.kernel_.length_scale
rbf = gaussian_process.kernels.RBF(length_scale=l)
 
Kernel = rbf(x_train, x_train)
K_y = Kernel + np.eye(training_number) * sigma_n
L = sp.linalg.cho_factor(K_y)
# cho_factor returns the (factor, lower_flag) pair that cho_solve expects;
# transposing that pair builds a ragged array (a hard error under NumPy 2).
alpha_p = sp.linalg.cho_solve(L, y_train)
    
k_s = rbf(x_test, x_train)
k_s_prime = (x_train.T - x_test) * k_s / l**2
f_prime = np.dot(k_s_prime, alpha_p) / (ub - lb)

/usr/local/lib/python3.6/dist-packages/numpy/core/_asarray.py:83: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order)

Calculate the BS delta

delta = lambda x, y: bsformula(1, lb+(ub-lb)*x, KC, r, T, y, 0)[1]
delta(x_test, sigma) - f_prime

array([[ 4.49894761e-04],
       [-1.48833441e-04],
       [-6.93913435e-05],
       [ 4.65385803e-05],
       [-4.98361790e-06],
       [-6.27400717e-05],
       [-8.62552503e-05],
       [-2.51381502e-05],
       [-1.47299604e-05],
       [-3.82167371e-05],
       [-7.22612021e-05],
       [-1.09925525e-05],
       [-7.32271155e-05],
       [-3.83089605e-05],
       [-3.29955983e-05],
       [-4.32648978e-05],
       [-4.68437001e-05],
       [-4.04393282e-05],
       [-2.67814069e-05],
       [-1.76363382e-05],
       [-1.94383865e-05],
       [-2.70635246e-05],
       [-3.05092406e-05],
       [-2.35631835e-05],
       [-1.48750391e-05],
       [-9.23441482e-06],
       [-1.68036306e-05],
       [-3.23303050e-05],
       [-4.23449670e-05],
       [-3.88301894e-05],
       [-2.62582636e-05],
       [-1.85345285e-05],
       [-2.87462109e-05],
       [-5.09098433e-05],
       [-6.50349044e-05],
       [-5.08255886e-05],
       [-3.78843033e-06],
       [ 5.67892481e-05],
       [ 1.10524320e-04],
       [ 1.55639129e-04],
       [ 2.19518691e-04],
       [ 3.43617331e-04],
       [ 5.47301648e-04],
       [ 8.15623470e-04],
       [ 1.10704450e-03],
       [ 1.39580542e-03],
       [ 1.69917364e-03],
       [ 2.07908613e-03],
       [ 2.58281389e-03],
       [ 3.18387898e-03]])

Compare the GP delta with the BS delta

plt.figure(figsize = (10,6),facecolor='white', edgecolor='black')
plt.plot(lb+(ub-lb)*x_test, delta(x_test,sigma), color = 'black', label = 'Exact')
plt.plot(lb+(ub-lb)*x_test, f_prime, color = 'red', label = 'GP')
plt.grid(True)
plt.xlabel('S')
plt.ylabel('$\Delta$')
plt.legend(loc = 'best', prop={'size':10});

Show the error between the GP delta and the BS delta

plt.figure(figsize = (10,6),facecolor='white', edgecolor='black')
plt.plot(lb+(ub-lb)*x_test, delta(x_test,sigma) - f_prime, color = 'black', label = 'GP Error')
plt.grid(True)
plt.xlabel('S')
plt.ylabel('Error in $\Delta$')
plt.legend(loc = 'best', prop={'size':10});

Vega

• Generate the training and testing data, where the input is the gridded underlying and the output are the option prices.
• The inputs are again scaled to the unit domain.

x_train = np.array(np.linspace(0.01, 1.2, training_number), dtype='float32').reshape(training_number, 1)
x_test = np.array(np.linspace(0.01, 1.0, testing_number), dtype='float32').reshape(testing_number, 1)

y_train = []
for idx in range(len(x_train)):
    y_train.append(call((S0-lb)/(ub-lb), x_train[idx]))
y_train = np.array(y_train)

Fit the GP model to the generated data

sk_kernel = RBF(length_scale=1.0, length_scale_bounds=(0.01, 10000.0))
gp = gaussian_process.GaussianProcessRegressor(kernel=sk_kernel, n_restarts_optimizer=20)
gp.fit(x_train, y_train)

GaussianProcessRegressor(kernel=RBF(length_scale=1), n_restarts_optimizer=20)

Get the model's predicted outputs for each of the test inputs

y_pred, sigma_hat = gp.predict(x_test, return_std=True)

/usr/local/lib/python3.6/dist-packages/sklearn/gaussian_process/_gpr.py:370: UserWarning: Predicted variances smaller than 0. Setting those variances to 0.
  warnings.warn("Predicted variances smaller than 0. "

Derive the GP vega

l = gp.kernel_.length_scale
rbf = gaussian_process.kernels.RBF(length_scale=l)
 
Kernel= rbf(x_train, x_train)
K_y = Kernel + np.eye(training_number) * sigma_n
L = sp.linalg.cho_factor(K_y)
# cho_factor returns the (factor, lower_flag) pair that cho_solve expects;
# transposing that pair builds a ragged array (a hard error under NumPy 2).
alpha_p = sp.linalg.cho_solve(L, y_train)
    
k_s = rbf(x_test, x_train)

k_s_prime = np.zeros([len(x_test), len(x_train)])
for i in range(len(x_test)):
    for j in range(len(x_train)):
        # x_train/x_test rows are length-1 arrays; NumPy 2 will not assign
        # one into a scalar slot, so take the element explicitly.
        k_s_prime[i, j] = (1.0/l**2) * (x_train[j].item() - x_test[i].item()) * k_s[i, j]
        
f_prime = np.dot(k_s_prime, alpha_p)

/usr/local/lib/python3.6/dist-packages/numpy/core/_asarray.py:83: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order)

Calculate the BS vega

vega = lambda x, y: bsformula(1, lb + (ub-lb) * x, KC, r, T, y, 0)[2]
vega((S0-lb)/(ub-lb), x_test) - f_prime

array([[-0.02976497],
       [-0.04088537],
       [ 0.05650961],
       [-0.01862286],
       [ 0.07929726],
       [-0.04356045],
       [ 0.0141282 ],
       [ 0.0549902 ],
       [-0.02229726],
       [ 0.01961521],
       [ 0.02189037],
       [ 0.00671363],
       [ 0.00702096],
       [ 0.01115929],
       [ 0.00410982],
       [ 0.0071057 ],
       [ 0.01868945],
       [ 0.00107479],
       [-0.00254281],
       [-0.00371633],
       [-0.00267023],
       [-0.00319543],
       [-0.00088096],
       [-0.00581987],
       [-0.00521972],
       [-0.00583219],
       [-0.00462355],
       [-0.00472223],
       [-0.00372284],
       [-0.00567935],
       [-0.0025467 ],
       [-0.00209934],
       [-0.00235141],
       [ 0.0006609 ],
       [ 0.00159   ],
       [-0.00203874],
       [-0.00097373],
       [ 0.00207862],
       [-0.00063463],
       [-0.00356933],
       [ 0.00375744],
       [-0.00283033],
       [-0.01011117],
       [-0.01430314],
       [-0.02503132],
       [-0.04477911],
       [-0.06429683],
       [-0.12263231],
       [-0.13711686],
       [-0.16173512]])

Compare the GP vega with the BS vega

plt.figure(figsize = (10,6), facecolor='white', edgecolor='black')
plt.plot(x_test, vega((S0-lb)/(ub-lb), x_test), color = 'black', label = 'Exact')
plt.plot(x_test, f_prime, color = 'red', label = 'GP')
plt.grid(True)
plt.xlabel('$\\sigma$')
plt.ylabel('$\\nu$')
plt.legend(loc = 'best', prop={'size':10});

Plot the error between the GP vega and the BS vega

plt.figure(figsize = (10,6), facecolor='white', edgecolor='black')
plt.plot(x_test, vega((S0-lb)/(ub-lb), x_test)-f_prime, color = 'black', label = 'GP Error')
plt.grid(True)
plt.xlabel('$\\sigma$')
plt.ylabel('Error in $\\nu$')
plt.legend(loc = 'best', prop={'size':10});