

In [1]:

    
import os, sys
import numpy as np
import pandas as pd
import datetime as dt

from IPython.display import display
from math import log, exp, sqrt, pow, erf, pi
from scipy.stats import norm

import ezvis3d as v3d

BlackScholes functions

The outputs below are computed from the following inputs

Inputs

Input	Symbol	Variable
Spot in currency	$$S$$	S
Strike in the same unit as Spot	$$K$$	K
Maturity in years	$$T$$	T
Volatility in % e.g. 15%=0.15	$$\sigma$$	v
Risk free interest rate in % e.g. 1.5%=0.015	$$r$$	r
Continuous dividend rate in % e.g. 2%=0.02	$$q$$	q

Outputs

Option	Call	Put
$$Payoff$$	$$Max(0, S-K)$$	$$Max(0, K-S)$$
$$Value=V$$	$$Se^{-qT}N(d_1)-Ke^{-rT}N(d_2)$$	$$Ke^{-rT}N(-d_2)-Se^{-qT}N(-d_1)$$
$$\Delta=\frac{\partial V}{\partial S}$$	$$e^{-qT}N(d_1)$$	$$-e^{-qT}N(-d_1)$$
$$\Gamma=\frac{\partial \Delta}{\partial S}$$	$$e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}$$	$$e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}$$
$$\nu=\frac{\partial V}{\partial \sigma}$$	$$Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}$$	$$Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}$$
$$\Theta=-\frac{\partial V}{\partial T}$$	$$-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}-rKe^{-rT}N(d_2)+qSe^{-qT}N(d_1)$$	$$-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}+rKe^{-rT}N(-d_2)-qSe^{-qT}N(-d_1)$$
$$\rho=\frac{\partial V}{\partial r}$$	$$KTe^{-rT}N(d_2)$$	$$-KTe^{-rT}N(-d_2)$$
$$Voma=\frac{\partial\nu}{\partial \sigma}$$	$$Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}$$	$$Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}$$

with

$$\frac{\partial V}{\partial t} + rS\frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2}=rV$$$$d_1=\frac{\log(S/K) + (r-q +\sigma^2/2)T}{\sigma \sqrt{T}}$$$$d_2=\frac{\log(S/K) + (r-q -\sigma^2/2)T}{\sigma \sqrt{T}}=d_1 - \sigma \sqrt{T - t}$$$$N(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-\frac{t^{2}}{2}}dt$$$$N'(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$$



In [2]:

    
def N(x):
    return (1.0 + erf(x / sqrt(2.0))) / 2.0

def Nprime(x):
    return exp(-x*x / 2.0) / sqrt(2.0 * pi)

# BlackScholes formula and Greeks - Call option
def BS_Greeks_Call(S, K, T, v, r, q):
    sqrt_T = sqrt(T)
    d1 = (log(S/K) + (r -q + v**2) * T) / (v *sqrt_T)
    d2 = d1 - v * sqrt_T
    N_d1 = N(d1)
    N_prime_d1 = Nprime(d1)
    N_d2 = N(d2)
    N_prime_d2 = Nprime(d2)
    PV = exp(-r * T)
    PV_K = K * PV
    D = exp(-q * T)
    value = N_d1 * S * D - N_d2 * PV_K
    
    delta = D * N_d1
    gamma = D * N_prime_d1 / (S * v * sqrt_T)
    vega = S *D * N_prime_d1 * sqrt_T
    theta = -D * S * N_prime_d1 * v / (2 * sqrt_T) - r * PV_K * N_d2 + q * S * D * N_d1
    rho = PV_K * T * N_d2
    voma = S * D * N_prime_d1 * sqrt_T * d1 * d2 / v

    payoff = max(0, S - K)
    PV_payoff = PV * payoff
 
    return {'d1': d1,
            'd2': d2,
            'N_d1_call': N_d1,
            'N_d2_call': N_d2,
            'N_prime_d1': N_prime_d1,
            'N_prime_d2': N_prime_d2,
            'PV': PV,
            'PV_K': PV_K,
            'D': D,
            'value': value,
            'delta': delta,
            'gamma': gamma,
            'vega': vega,
            'theta': theta,
            'rho': rho,
            'voma': voma,
            'payoff': payoff,
            'PV_payoff': PV_payoff}

# S = 100.0
# K = 100.0
# T = 5.0
# vol = 0.25
# r = 0.02
# q = 0.03
# typ = +1
# BS_Greeks_Call(S, K, T, vol, r, q)



In [3]:

    
# BlackScholes formula and Greeks - Put option
def BS_Greeks_Put(S, K, T, v, r, q):
    sqrt_T = sqrt(T)
    d1 = (1 / (v *sqrt_T)) * (log(S/K) + (r -q + v*v) * T)
    d2 = d1 - v * sqrt_T
    N_minus_d1 = N(-d1)
    N_prime_d1 = Nprime(d1)
    N_minus_d2 = N(-d2)
    N_prime_d2 = Nprime(d2)
    PV = exp(-r * T)
    PV_K = K * PV
    D = exp(-q * T)
    value = N_minus_d2 * PV_K - N_minus_d1 * S * D
    
    delta = - D * N_minus_d1
    gamma = D * N_prime_d1 / (S * v * sqrt_T)
    vega = S * D * N_prime_d1 * sqrt_T
    theta = -D * S * N_prime_d1 * v / (2 * sqrt_T) + r * PV_K * N_minus_d2 - q * S * D * N_minus_d1
    rho = - PV_K * T * N_minus_d2
    voma = S * D * N_prime_d1 * sqrt_T * d1 *d2 / v

    payoff = max(0, K - S)
    PV_payoff = PV * payoff

    return {'d1': d1,
            'd2': d2,
            'N_d1_put': N_minus_d1,
            'N_d2_put': N_minus_d2,
            'N_prime_d1': N_prime_d1,
            'N_prime_d2': N_prime_d2,
            'PV': PV,
            'PV_K': PV_K,
            'D': D,
            'value': value,
            'delta': delta,
            'gamma': gamma,
            'vega': vega,
            'theta': theta,
            'rho': rho,
            'voma': voma,
            'payoff': payoff,
            'PV_payoff': PV_payoff}



# S = 100.0
# K = 100.0
# T = 5.0
# vol = 0.25
# r = 0.02
# q = 0.03
# typ = +1
# BS_Greeks_Put(S, K, T, vol, r, q)

Plotting

Compute one output (z) as a function of 2 variables (x, y) among S, K, T, v, r, q while the others are kept fixed
The output must be a key in the dictionary returned by functions BS_Greeks_Put() or BS_Greeks_Call()
Pass argument save=True to function plot() as standalone HTML doc under ./saved



In [4]:

    
# S=Spot, K=Strike, T=mat, v=vol, r=rate, q=div rate
S = 100           
K = 100
T = 3
v = 15 / 100.0
r = 2 / 100.0
q = 1 / 100.0

x_var = 'Spot' # --------------------- choose x axis name
y_var = 'Maturity' # ----------------- choose y axis name
z = 'value' # ------------------------ choose z axis

def func(x, y):
    S_loc, K_loc, T_loc, v_loc, r_loc, q_loc = S, K, T, v, r, q
    S_loc = x #-------------- choose x axis
    T_loc = y #-------------- choose y axis
    res = BS_Greeks_Call(S_loc, K_loc, T_loc, v_loc, r_loc, q_loc)
    return res[z]

x_min, x_max, x_num = 0.05, 200, 50 # ---------------- choose x axis boundaries
y_min, y_max, y_num = 0.05, 10, 50 # ----------------- choose y axis boundaries

x_rng = np.linspace(x_min, x_max, x_num)
y_rng = np.linspace(y_min, y_max, y_num)
li_data = [{'x': x, 'y': y, 'z': func(x, y)}
            for y in y_rng
            for x in x_rng]
df_data = pd.DataFrame(li_data)

g = v3d.Vis3d()
g.width = '700px'
g.height = '700px'
g.style = 'surface'
g.showPerspective = True
g.showGrid = True
g.showShadow = False
g.keepAspectRatio = False
g.verticalRatio = 0.8
g.xLabel = x_var
g.yLabel = y_var
g.zLabel = z
g.cameraPosition = {'horizontal' : 0.9,
                    'vertical': 0.5,
                    'distance': 1.8}

g.plot(df_data, center=True, save=False)









    Out[4]:



In [5]:

    
# S=Spot, K=Strike, T=mat, v=vol, r=rate, q=div rate
S = 100           
K = 100
T = 3
v = 15 / 100.0
r = 2 / 100.0
q = 1 / 100.0

x_var = 'Spot' # --------------------- choose x axis name
y_var = 'Maturity' # ----------------- choose y axis name

def func(x, y):
    S_loc, K_loc, T_loc, v_loc, r_loc, q_loc = S, K, T, v, r, q
    S_loc = x #-------------- choose x axis
    T_loc = y #-------------- choose y axis
    res = BS_Greeks_Call(S_loc, K_loc, T_loc, v_loc, r_loc, q_loc)
    return res

x_min, x_max, x_num = 0.05, 200, 50 # ----------------- choose x axis boundaries
y_min, y_max, y_num = 0.05, 10, 50 # ----------------- choose y axis boundaries

x_rng = np.linspace(x_min, x_max, x_num)
y_rng = np.linspace(y_min, y_max, y_num)
li_data = [{'x': x, 'y': y, 'z': func(x, y)}
            for y in y_rng
            for x in x_rng]

for z in ['d1', 'd2', 'N_d1_call', 'N_d2_call', 'PV_K', 'value', 'delta',
          'gamma', 'vega', 'theta', 'rho', 'voma', 'payoff', 'PV_payoff']:
    li_data_z = [{'x': e['x'], 'y': e['y'], 'z': e['z'][z]} for e in li_data]
    df_data = pd.DataFrame(li_data_z)

    g = v3d.Vis3d()
    g.width = '400px'
    g.height = '400px'
    g.style = 'surface'
    g.showPerspective = True
    g.showGrid = True
    g.showShadow = False
    g.keepAspectRatio = False
    g.verticalRatio = 0.8
    g.xLabel = x_var
    g.yLabel = y_var
    g.zLabel = z
    g.cameraPosition = {'horizontal' : 0.9,
                        'vertical': 0.4,
                        'distance': 1.8}

    print('{} as a function of {} and {}'.format(z, x_var, y_var))
    display(g.plot(df_data, center=True, save=False))









    



d1 as a function of Spot and Maturity






    





        
        






    



d2 as a function of Spot and Maturity






    





        
        






    



N_d1_call as a function of Spot and Maturity






    





        
        






    



N_d2_call as a function of Spot and Maturity






    





        
        






    



PV_K as a function of Spot and Maturity






    





        
        






    



value as a function of Spot and Maturity






    





        
        






    



delta as a function of Spot and Maturity






    





        
        






    



gamma as a function of Spot and Maturity






    





        
        






    



vega as a function of Spot and Maturity






    





        
        






    



theta as a function of Spot and Maturity






    





        
        






    



rho as a function of Spot and Maturity






    





        
        






    



voma as a function of Spot and Maturity






    





        
        






    



payoff as a function of Spot and Maturity






    





        
        






    



PV_payoff as a function of Spot and Maturity



In [ ]:



In [ ]:



In [ ]:

Annex: Create Markdown



In [6]:

    
from jinja2 import Environment

ref_template_1 = """
| **Option** | **Call** | **Put** |
| :--------: |:--------:|:-------:|
{% for x in output_data %}| $${{x[0]}}$$   | $${{x[1]}}$$  | $${{x[2]}}$$   |
{% endfor %}
"""

ref_data_1 = [
	[r'Payoff', r'Max(0, S-K)', r'Max(0, K-S)'],
	[r"Value=V", r"Se^{-qT}N(d_1)-Ke^{-rT}N(d_2)", r'Ke^{-rT}N(-d_2)-Se^{-qT}N(-d_1)'],
	[r"\Delta=\frac{\partial V}{\partial S}", r"e^{-qT}N(d_1)", r'-e^{-qT}N(-d_1)'],
	[r"\Gamma=\frac{\partial \Delta}{\partial S}", r"e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}", r"e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}"],
	[r"\nu=\frac{\partial V}{\partial \sigma}", r"Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}", r"Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}"],
	[r"\Theta=-\frac{\partial V}{\partial T}", r"-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}-rKe^{-rT}N(d_2)+qSe^{-qT}N(d_1)", r"-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}+rKe^{-rT}N(-d_2)-qSe^{-qT}N(-d_1)"],
	[r"\rho=\frac{\partial V}{\partial r}", r"KTe^{-rT}N(d_2)", r"-KTe^{-rT}N(-d_2)"],
	[r"Voma=\frac{\partial\nu}{\partial \sigma}", r"Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}", r"Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}"],
]

res = Environment().from_string(ref_template_1).render(output_data=ref_data_1)
print(res)









    



| **Option** | **Call** | **Put** |
| :--------: |:--------:|:-------:|
| $$Payoff$$   | $$Max(0, S-K)$$  | $$Max(0, K-S)$$   |
| $$Value=V$$   | $$Se^{-qT}N(d_1)-Ke^{-rT}N(d_2)$$  | $$Ke^{-rT}N(-d_2)-Se^{-qT}N(-d_1)$$   |
| $$\Delta=\frac{\partial V}{\partial S}$$   | $$e^{-qT}N(d_1)$$  | $$-e^{-qT}N(-d_1)$$   |
| $$\Gamma=\frac{\partial \Delta}{\partial S}$$   | $$e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}$$  | $$e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}$$   |
| $$\nu=\frac{\partial V}{\partial \sigma}$$   | $$Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}$$  | $$Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}$$   |
| $$\Theta=-\frac{\partial V}{\partial T}$$   | $$-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}-rKe^{-rT}N(d_2)+qSe^{-qT}N(d_1)$$  | $$-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}+rKe^{-rT}N(-d_2)-qSe^{-qT}N(-d_1)$$   |
| $$\rho=\frac{\partial V}{\partial r}$$   | $$KTe^{-rT}N(d_2)$$  | $$-KTe^{-rT}N(-d_2)$$   |
| $$Voma=\frac{\partial\nu}{\partial \sigma}$$   | $$Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}$$  | $$Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}$$   |



In [7]:

    
from jinja2 import Environment

ref_template_2 = """
{% for x in output_data %}
$${{x[0]}}$$
{% endfor %}
"""

ref_data_2 = [
	[r"\frac{\partial V}{\partial t} + rS\frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2}=rV"],
	[r"d_1=\frac{\log(S/K) + (r-q +\sigma^2/2)T}{\sigma \sqrt{T}}"],
	[r"d_2=\frac{\log(S/K) + (r-q -\sigma^2/2)T}{\sigma \sqrt{T}}=d_1 - \sigma \sqrt{T - t}"],
	[r"N(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-\frac{t^{2}}{2}}dt"],
	[r"N'(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}"]
]


res = Environment().from_string(ref_template_2).render(output_data=ref_data_2)
print(res)









    




$$\frac{\partial V}{\partial t} + rS\frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2S^2\frac{\partial^2 V}{\partial S^2}=rV$$

$$d_1=\frac{\log(S/K) + (r-q +\sigma^2/2)T}{\sigma \sqrt{T}}$$

$$d_2=\frac{\log(S/K) + (r-q -\sigma^2/2)T}{\sigma \sqrt{T}}=d_1 - \sigma \sqrt{T - t}$$

$$N(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-\frac{t^{2}}{2}}dt$$

$$N'(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$$