#!/usr/bin/env python # coding: utf-8 # 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) # 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)) # 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) # 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)