Please see the attached code for an implementation of a standard Black-Scholes model for options pricing and risk management. Includes functions for valuation of first, second, and third order Greeks.

Clone Algorithm

Loading...

There was an error loading this backtest.

Backtest from
to
with
initial capital

Cumulative performance:

Algorithm
Benchmark

Custom data:

Total Returns

--

Alpha

--

Beta

--

Sharpe

--

Sortino

--

Max Drawdown

--

Benchmark Returns

--

Volatility

--

Returns | 1 Month | 3 Month | 6 Month | 12 Month |

Alpha | 1 Month | 3 Month | 6 Month | 12 Month |

Beta | 1 Month | 3 Month | 6 Month | 12 Month |

Sharpe | 1 Month | 3 Month | 6 Month | 12 Month |

Sortino | 1 Month | 3 Month | 6 Month | 12 Month |

Volatility | 1 Month | 3 Month | 6 Month | 12 Month |

Max Drawdown | 1 Month | 3 Month | 6 Month | 12 Month |

import numpy as np from math import log, sqrt, erf, exp import scipy.stats # Library for valuation of European call and put options using the # standard Black-Scholes model # S is the current stock price # K is the strike price # r is the risk-free rate # sigma is the annualized volatility # T is the time to expiration # q is the continuous dividend yield # call is a flag indicating whether the option is a call (1) or put (0) def d1(S,K,r,sigma,T,q): num = log(S/K)+(r-q+0.5*sigma**2)*T den = sigma*sqrt(T) return num/den def d2(S,K,r,sigma,T,q): num = log(S/K)+(r-q-0.5*sigma**2)*T den = sigma*sqrt(T) return num/den # Cumulative distribution function for the standard normal distribution def phi(x): return (1.0 + erf(x / sqrt(2.0))) / 2.0 # Current value of the option given market parameters def V(S,K,r,sigma,T,q,call): if call == 1: return S*exp(-q*T)*phi(d1(S,K,r,sigma,T,q)) - K*exp(-r*T)*phi(d2(S,K,r,sigma,T,q)) else: return exp(-r*T)*K*phi(-d2(S,K,r,sigma,T,q)) - S*exp(-q*T)*phi(-d1(S,K,r,sigma,T,q)) # Delta is the instantaneous rate of change of the option value with respect # to the price of the underlying def delta(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) if call == 1: return exp(-q*T)*scipy.stats.norm(0, 1).cdf(d1_new) else: return -exp(-q*T)*scipy.stats.norm(0, 1).cdf(-d1_new) # Gamma is the instantaneous rate of change of delta with respect to the price # of the underlying def gamma(S,K,r,sigma,T,q,call): # Same for both call and put d1_new = d1(S,K,r,sigma,T,q) return exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)/(S*sigma*sqrt(T)) # Vega measures the instantaneous rate of change of the option value with # respect to volatility. def vega(S,K,r,sigma,T,q,call): # Same for both call and put d2_new = d2(S,K,r,sigma,T,q) return K*exp(-r*T)*scipy.stats.norm(0, 1).pdf(d2_new)*sqrt(T)/100 # Theta, or time decay, measures the instantaneous rate of change of the option # value with respect to the passage of time. def theta(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) d2_new = d2(S,K,r,sigma,T,q) first_term = exp(-q*T)*S*sigma*scipy.stats.norm(0, 1).pdf(d1_new)/(2*sqrt(T)) second_term = r*exp(-r*T)*K*phi(d2_new) third_term = q*S*exp(-q*T)*phi(d1_new) if call == 1: return (-first_term - second_term + third_term)/365 else: return (-first_term + second_term - third_term)/365 # Rho measures sensitivity to the interest rate. It is the instantaneous rate # of change of the option value with respect to the risk-free interest rate. def rho(S,K,r,sigma,T,q,call): d2_new = d2(S,K,r,sigma,T,q) if call == 1: return K*T*exp(-r*T)*phi(d2_new)/100 else: return -K*T*exp(-r*T)*phi(-d2_new)/100 # Higher order Greeks # Charm, or delta decay, measures the instantaneous rate of change of delta # with the passage of time def charm(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) d2_new = d2(S,K,r,sigma,T,q) first_term = q*exp(-q*T)*phi(d1_new) second_term = exp(-q*T)*scipy.stats.norm(0, 1).pdf(d1_new)*(2*(r-q)*T - d2_new*sigma*sqrt(T))/(2*T*sigma*sqrt(T)) if call == 1: return first_term - second_term else: return - first_term - second_term # Speed is the rate of change of gamma with respect to the price of the # underlying def speed(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) gamma_new = gamma(S,K,r,sigma,T,q,call) return (-gamma_new/S)*(d1_new/(sigma*sqrt(T))+1) # Zomma measures the rate of change of gamma with respect to changes in volatility def zomma(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) d2_new = d2(S,K,r,sigma,T,q) gamma_new = gamma(S,K,r,sigma,T,q,call) return gamma_new*(d1_new*d2_new - 1)/sigma # Vomma measures second-order exposure to volatility. Vomma is the instantaneous # rate of change of vega with respect to volatility def vomma(S,K,r,sigma,T,q,call): d1_new = d1(S,K,r,sigma,T,q) d2_new = d2(S,K,r,sigma,T,q) vega_new = vega(S,K,r,sigma,T,q,call) return vega_new*d1_new*d2_new/sigma def initialize(context): context.stock = sid(24) def handle_data(context, data): pass

This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.