simple example simulation of delta hedging with python

Here we will present simple python code of delta hedging example of a call option .
it’s a minimal example with zero interest rates , no dividends. On day 1 we sell 10 near ATM call options and start delta hedging i.e. buying/selling stock so that change in stock price neutralizes change in options value.The portfolio is then $$\pi= delta*Stock-callOption+cashAccount$$

to simulate stock prices we will use log-normal dynamics $$ \frac{dS_t}{S_t} =\sigma_t dW_t $$

We will also simulate implied volatility as log-normal $$ \frac{d\sigma_t}{\sigma_t}=\xi dB_t $$
each day of simulation we will store in DataFrame df , so it will be easy to print and plot with pandas library.

to run the python code you will need pandas library installed in your distribution

apart from usual NPV of hedging portfolio and values of call option and greeks we will also display P&L prediction which is $$ PNL \approx \frac{\Gamma}{2} (dS)^2 + \theta dt $$

delta hedging simulation

first 10 days of delta hedging:

       spot   vol  shares      cash  option   npv    vega  gamma  theta  pnlPredict   pnl  error
0    100.00 30.00    0.00      0.00    0.00  0.00    0.00   0.00   0.00        0.00   nan    nan
1    102.89 29.83    6.35   -401.20  252.38  0.00 -772.52  -0.06   0.12       -0.14  0.00   0.00
2    103.65 29.89    6.40   -406.13  257.56 -0.36 -774.36  -0.06   0.12        0.10 -0.36  -0.46
3     98.00 29.34    6.02   -368.79  218.14  2.91 -754.75  -0.07   0.11       -0.96  3.27   4.23
4     96.83 29.61    5.95   -362.24  213.06  0.93 -748.61  -0.07   0.11        0.06 -1.98  -2.04
5     95.76 29.35    5.86   -353.97  204.68  2.96 -743.77  -0.07   0.11        0.07  2.03   1.96
6     95.40 28.75    5.81   -348.68  198.03  7.53 -742.86  -0.07   0.11        0.10  4.57   4.46
7     94.12 28.71    5.71   -339.71  190.26  7.85 -735.94  -0.07   0.11        0.05  0.32   0.28
8     94.28 27.96    5.69   -336.97  185.56 13.47 -737.79  -0.07   0.10        0.10  5.62   5.52
9     95.42 27.56    5.75   -343.04  189.00 16.52 -744.17  -0.07   0.10        0.05  3.05   3.00
10    96.53 28.02    5.85   -353.01  198.72 13.16 -748.41  -0.07   0.11        0.06 -3.36  -3.42

P&L profile

running this code for number of days 10,100,1000,10000 we can see P&L profile start to look narrower and more centered at 0 with number of days , in the limit of infinite days (equivalent to infinite frequency of hedging of Black Scholes assumptions) it would look like a dirac function.

pnl profile delta hedging 10 100 1000 10000 days

Why is it skewed to the left? as Gamma is always negative and high near ATM, if there’s jump down, PNL is proportional to gamma * jump size $$ PNL \approx \frac{\Gamma S^2}{2}((\frac{dS}{S})^2-\sigma^2)$$ and delta makes a high decrease, so we would need to sell many shares cheap. In case there’s a jump up, we would still need to buy more stock (as delta of call option goes up), so still would loose money in this case.That makes the daily PNL profile asymmetric.

Theta

setting Spot prices constant in the code we can see the Theta effect on the hedging strategy:
here’s the example of running the python code with constant Spot:

theta OTM delta hedging constant spot

volaility dynamics

here we will run the code first on constant volaility , and then on stochastic implied volaility:
we can see that pnl has much wider distribution in the latter case. this is due to the fact that black sholes assumption does not hold anymore (we would need to hedge with other options to get better pnl profile)

constant vs stochastic volatility pnl delta hedging

import numpy as np
import pandas as pd
pd.set_option('display.width', 320)
pd.set_option('display.max_rows', 100)
pd.options.display.float_format = '{:,.2f}'.format
from scipy.stats import norm
import matplotlib.pyplot as plt

def BlackScholes(tau, S, K, sigma):
    d1=np.log(S/K)/sigma/np.sqrt(tau)+0.5*sigma*np.sqrt(tau)
    d2=d1-sigma*np.sqrt(tau)
    npv=(S*norm.cdf(d1)-K*norm.cdf(d2))
    delta=norm.cdf(d1)
    gamma=norm.pdf(d1)/(S*sigma*np.sqrt(tau))
    vega=S*norm.pdf(d1)*np.sqrt(tau)
    theta=-.5*S*norm.pdf(d1)*sigma/np.sqrt(tau)
    return {'npv':npv,'delta':delta,'gamma':gamma,'vega':vega,'theta':theta}

class Call(object):
    def __init__(self,start,T,K,N):
        self.T=T
        self.K=K
        self.start=start  #day to sell   option
        self.N=N

    def calc(self,today,vol,S):
        if today<self.start:
            return {'delta':0,'npv':0,'vega':0,'gamma':0,'theta':0,'intrinsic':0}
        if today>self.T:
            return {'delta':0,'npv':0,'vega':0,'gamma':0,'theta':0,'intrinsic':0}
        if today==self.T:
            return {'delta':0,'npv':0,'vega':0,'gamma':0,'theta':0,'intrinsic':self.N*max(0,S-self.K)}
        tau=(self.T-today)/250.
        call=BlackScholes(tau, S, self.K, vol)
        return {'delta':self.N*call['delta'],'npv':self.N*call['npv'],'vega':self.N*call['vega'],'gamma':self.N*call['gamma'],'theta':self.N*call['theta'],'intrinsic':self.N*max(0,S-self.K)}

def deltahedge(Ndays,Sdynamics="S*=(1.0+vol*np.sqrt(dt)*np.random.randn())",volDynamics="vol=.30"):
    S=90.0
    Strike=100.0
    vol=.30
    columns=('spot','vol','shares','cash','option','npv','vega','gamma','theta','pnlPredict')
    df = pd.DataFrame([[S,vol,0,0,0,0,0,0,0,0]],columns=columns)
    dt=1/250.
    cash=0
    dayToSellCall=1
    maturityCall=Ndays-1
    call=Call(dayToSellCall,maturityCall,Strike,-10)# sell one call on dayToSellCall day
    for day in xrange(1,Ndays+1):
        exec Sdynamics
        exec volDynamics

        if day==dayToSellCall: #sell call
          callValue=call.calc(day,vol,S)
          cash-=callValue['npv']

        #delta hedge
        callValue=call.calc(day,vol,S)
        delta=callValue['delta']
        currentNumberShares=df.iloc[day-1].shares
        sharesBuy=-currentNumberShares-delta
        cash-=sharesBuy*S
        if day==maturityCall:
            cash+=call.calc(day,vol,S)['intrinsic'] #settle call

        gamma=callValue['gamma']
        theta=callValue['theta']
        dS=S-df.iloc[day-1].spot
        pnlPredict=0.5*gamma*dS*dS+theta*dt
        dfnew=pd.DataFrame([[S,vol,-delta,cash,-callValue['npv'],cash+callValue['npv']-delta*S,callValue['vega'],gamma,theta/250.,pnlPredict]],columns=columns)
        df=df.append(dfnew,ignore_index=True)

    df['pnl'] = df['npv'] - df['npv'].shift(1)
    df['vol']=100.0*df['vol']
    df['error']=df['pnl']-df['pnlPredict']
    df.set_value(dayToSellCall, "error", 0)
    #df.loc[:,['vol','spot']].plot(title='Spot and implied Volatility')
    df.loc[:,['npv','spot','option']].plot(title='-Call+delta*S+cash vs Spot {0} {1}'.format(Sdynamics,volDynamics))
    df.loc[:,['theta']].plot(title='Theta {0} {1}'.format(Sdynamics,volDynamics))
    df.loc[:,['pnl']].hist(bins=50)
    #df.loc[:,['error']].hist(bins=50)
    print df.loc[:,['pnl']].describe()
    #print df

deltahedge(1000)#constant vol
deltahedge(1000,volDynamics="vol*=(1.0+0.5*np.sqrt(dt)*np.random.randn())")#stochastic vol
deltahedge(1000,Sdynamics="S=90") #consant OTM
deltahedge(10000,Sdynamics="S=100") #consant ATM
deltahedge(10000,Sdynamics="S=110")#consant ITM
deltahedge(10000,Sdynamics="S+=S*dt") #growing stock
plt.show()