importing modules, and some basic stuff
In [9]:
# always yielding a real result, even dividing two integers
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats as st
# plot inline
%matplotlib inline
function for a monte carlo path sample, based on Generalized geometric Brownian motion:
In [2]:
def mcpath(nsteps, S0, r, vol, T):
# generate nsteps random numbers
sample=pd.Series(np.random.standard_normal(nsteps))
# define time increment
dt=T/nsteps
# initialize vectors
initial = np.zeros(1)
ts=np.repeat(dt,nsteps)
# calculate s1 and s2
s1=(r-(vol**2)/2)*dt
s2=vol*np.sqrt(dt)
steps=s1*np.repeat(1,nsteps)+s2*sample
path=steps.cumsum()
Sj=S0*np.exp(path)
# add initial time t=0
S = np.concatenate(([S0],Sj))
t = np.concatenate(([0],ts.cumsum()))
return pd.Series(S,index=t)
An auxiliary function to generate npaths
In [3]:
def mcpaths(npaths, nsteps, S0, r, vol, t):
# generate npaths using mcpath
paths=[mcpath(nsteps, S0, r, vol, t) for j in range(npaths)]
return paths
Let´s test
In [209]:
npaths=10
nsteps=100
S0=100
r=0.1
vol=0.2
T=1
test_mcpaths=mcpaths(npaths,nsteps,S0,r,vol,T)
In [210]:
[test_mcpaths[i].plot() for i in range(npaths)]
Out[210]:
Defining the Black formula, according to:
Where:
$$F = S_0e^{\left(r-q\right)t}$$$$d_1 = \frac{1}{\sigma\sqrt{t}} \left[ \ln{\left( \frac{S_0}{K} \right) } + \left( r -q + \frac{\sigma^2}{2} \right) \left(t \right) \right]$$$$d_2 = d_1 - \sigma\sqrt{t}$$
In [211]:
def bsv(phi,S,K,r,q,vol,t):
if(t>0):
# calculate F
F=S*np.exp((r-q)*t)
# calculate d1 and d2
sigma_rt = vol*np.sqrt(t)
d1=(np.log(S/K)+(r-q+vol**2/2)*t)/(sigma_rt)
d2=d1-sigma_rt
# calculate N(d1) and N(d2)
Nd1=st.norm.cdf(phi*d1)
Nd2=st.norm.cdf(phi*d2)
# calculate delta and premium
delta=phi*np.exp(-q*t)*Nd1
premium=S*delta-phi*K*np.exp(-r*t)*Nd2
else:
delta=0
premium=max(phi*(S-K),0)
return [premium,delta]
Testing the function we can verify that results match with this Black-Scholes Calculator from Drexel University
In [212]:
call = 1
put = -1
S0=100
r=0.1
q=.02
vol=0.2
T=1
K=100
[bsv(call,S0,K,r,q,vol,T),bsv(put,S0,K,r,q,vol,T)]
Out[212]:
Function for calculating portfolio P&L and cash flow value given a path
In [213]:
def calcportfolio(path,phi,K,r,q,vol,T):
nstp=len(path)-1
# calculate t,S(t),premium(t),delta(t) using Black Scholes function
bsvpath=[[path.index[i],path.values[i]]+bsv(call,path.values[i],K,r,q,vol,1-path.index[i])\
for i in range(len(path))]
df = pd.DataFrame(bsvpath,columns=['time','spot','premium','delta'])
# generate cfwprem column to include cash flows regarding premium
df['cfwprem']=0
# the first cash flow is the price payed to enter the option
df.loc[0,'cfwprem']=-df['premium'][0]
# at the maturity:
# - all premium, if any, becomes cash flow regarding option settlement
# - there is no delta as option expires and we should withdraw the position on asset
df.loc[nstp,'cfwprem']=df['premium'][nstp]
df.loc[nstp,'premium']=0
df.loc[nstp,'delta']=0
# calculate time intervals dt
df['timechg']=df['time'].diff()
df.loc[0,'timechg']=0
# calculate changes in delta
df['dltchg']=df['delta'].diff()
df.loc[0,'dltchg']=0
# calculate changes in spot price
df['spotchg']=df['spot'].diff()
df.loc[0,'spotchg']=0
# cashflows for heding the portfolio buying/selling delta quantities of the asset
df['cfwspot']=0
df.loc[0,'cfwspot']=df['delta'][0]*df['spot'][0]
df.loc[1:,'cfwspot']=df['dltchg'][1:]*df['spot'][1:]
# dividend cashflows
df['cfwdivid']=0
df.loc[1:,'cfwdivid']=-((df['delta'][0:nstp]*df['spot'][0:nstp]).values)*(np.exp(q*df['timechg'][1:].values)-1)
# cashflows before interest
df['cfwprer']=df['cfwprem']+df['cfwspot']+df['cfwdivid']
# interest and consolidation of cashflows
df['balance']=0
df.loc[0,'balance']=df['cfwprer'][0]
for j in range(1,nstp+1):
df.loc[j,'balance']=df['balance'][j-1]*(np.exp(r*df['timechg'][j]))+df['cfwprer'][j]
# portfolio
df['portf']=df['premium']-df['delta']*df['spot']+df['balance']
# consolidated discount factor
return df
some function test, just to see the data frame´s face
In [214]:
call=1
put=-1
nsteps=10
S0=100
r=0.1
q=0
vol=0.2
T=1
K=100
mypath = mcpath(nsteps,S0,r,vol,T)
calcportfolio(mypath,call,K,r,q,vol,T)
Out[214]:
generating a histogram of portfolio P&L
In [215]:
npaths=500
pls=pd.Series([calcportfolio(mcpath(nsteps, S0, r, vol, T),call,K,r,q,vol,T)['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=30)
Out[215]:
Change the the Black&Scholes and portfolio function to receive implied and realized volatity
In [216]:
def bsv2(phi,S,K,r,q,vol_pricing,vol_delta,t):
if(t>0):
# calculate F
F=S*np.exp((r-q)*t)
# calculate d1 and d2
sigma_rt = vol_pricing*np.sqrt(t)
d1=(np.log(S/K)+(r-q+vol_pricing**2/2)*t)/(sigma_rt)
d2=d1-sigma_rt
# calculate N(d1) and N(d2)
Nd1=st.norm.cdf(phi*d1)
Nd2=st.norm.cdf(phi*d2)
# calculate premium
delta=phi*np.exp(-q*t)*Nd1
premium=S*delta-phi*K*np.exp(-r*t)*Nd2
# calculate delta
sigma_rtdelta=vol_delta*np.sqrt(t)
d1delta=(np.log(S/K)+(r-q+vol_delta**2/2)*t)/(sigma_rtdelta)
Nd1delta=st.norm.cdf(phi*d1delta)
delta=phi*np.exp(-q*t)*Nd1delta
else:
delta=0
premium=max(phi*(S-K),0)
return [premium,delta]
In [217]:
def calcportfolio2(path,phi,K,r,q,vol_pricing,vol_delta,T):
nstp=len(path)-1
# calculate t,S(t),premium(t),delta(t) using Black Scholes function
bsvpath=[[path.index[i],path.values[i]]+bsv2(call,path.values[i],K,r,q,vol_pricing,vol_delta,1-path.index[i])\
for i in range(len(path))]
df = pd.DataFrame(bsvpath,columns=['time','spot','premium','delta'])
# generate cfwprem column to include cash flows regarding premium
df['cfwprem']=0
# the first cash flow is the price payed to enter the option
df.loc[0,'cfwprem']=-df['premium'][0]
# at the maturity:
# - all premium, if any, becomes cash flow regarding option settlement
# - there is no delta as option expires and we should withdraw the position on asset
df.loc[nstp,'cfwprem']=df['premium'][nstp]
df.loc[nstp,'premium']=0
df.loc[nstp,'delta']=0
# calculate time intervals dt
df['timechg']=df['time'].diff()
df.loc[0,'timechg']=0
# calculate changes in delta
df['dltchg']=df['delta'].diff()
df.loc[0,'dltchg']=0
# calculate changes in spot price
df['spotchg']=df['spot'].diff()
df.loc[0,'spotchg']=0
# cashflows for heding the portfolio buying/selling delta quantities of the asset
df['cfwspot']=0
df.loc[0,'cfwspot']=df['delta'][0]*df['spot'][0]
df.loc[1:,'cfwspot']=df['dltchg'][1:]*df['spot'][1:]
# dividend cashflows
df['cfwdivid']=0
df.loc[1:,'cfwdivid']=-((df['delta'][0:nstp]*df['spot'][0:nstp]).values)*(np.exp(q*df['timechg'][1:].values)-1)
# cashflows before interest
df['cfwprer']=df['cfwprem']+df['cfwspot']+df['cfwdivid']
# interest and consolidation of cashflows
df['balance']=0
df.loc[0,'balance']=df['cfwprer'][0]
for j in range(1,nstp+1):
df.loc[j,'balance']=df['balance'][j-1]*(np.exp(r*df['timechg'][j]))+df['cfwprer'][j]
# portfolio
df['portf']=df['premium']-df['delta']*df['spot']+df['balance']
# consolidated discount factor
return df
Define parameters
In [218]:
call=1
put=-1
S0=100
r=0
q=0
vol_pricing=.2
vol_delta=.3
T=1
K=100
nsteps=500
npaths=500
Generate paths
In [219]:
paths = [mcpath(nsteps, S0, r, vol_delta, T) for i in range(npaths)]
P&L portfolio envolution
In [220]:
portfsVolReal=[calcportfolio2(paths[i],call,K,r,q,vol_pricing,vol_delta,T) for i in range(npaths)]
df=pd.DataFrame([pd.Series(portfsVolReal[i]['portf'].values,index=portfsVolReal[i]['time']) for i in range(npaths)])
df.transpose().plot(legend=None)
Out[220]:
Final P&L distribution
In [221]:
pls=pd.Series([portfsVolReal[i]['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=10)
Out[221]:
In [222]:
portfsVolImpl=[calcportfolio2(paths[i],call,K,r,q,vol_pricing,vol_pricing,T) for i in range(npaths)]
df=pd.DataFrame([pd.Series(portfsVolImpl[i]['portf'].values,index=portfsVolImpl[i]['time']) for i in range(npaths)])
df.transpose().plot(legend=None)
Out[222]:
Final P&L distribution
In [223]:
pls=pd.Series([portfsVolImpl[i]['portf'][nsteps] for i in range(npaths)])
pls.plot(kind='hist',bins=10)
Out[223]:
In [224]:
import scipy.optimize as opt
In [242]:
def bev(vol,path,phi,K,r,q,T):
plfinal = calcportfolio(path,phi,K,r,q,vol,T)['portf'][len(path)-1]
return plfinal
In [243]:
bevs = pd.Series([opt.newton(bev, 0.2, args=(paths[i],call,K,r,q,T)) for i in range(npaths)])
The Break Event Vol distribution
In [244]:
bevs.plot(kind='hist',bins=20)
Out[244]:
In [245]:
def calcvolreal(path,T):
return np.std(path.values[1:] / path.values[0:-1] - 1)*np.sqrt(len(path)/T)
In [246]:
vols_real = [calcvolreal(paths[i],T) for i in range(npaths)]
Scatter plot between Realized Vol and Break Event Vol
In [247]:
plt.plot(vols_real,bevs.values,'o')
plt.xlabel('rvol')
plt.ylabel('bev')
Out[247]:
Function that returns Strike value given a Delta value
In [268]:
def strike_from_delta(phi,S,delta,r,q,vol,t):
Nd1=delta/(phi*np.exp(-q*t))
d1=st.norm.ppf(phi*Nd1)
sigma_rt=vol*np.sqrt(t)
K=S/(np.exp((d1*sigma_rt)-(r-q+vol**2/2)*t))
return K
Calculate strike values for delta 25, delta 10, delta -10 and delta -25
In [274]:
strikes=[strike_from_delta(call,S0,0.1,r,q,vol,T),\
strike_from_delta(call,S0,0.25,r,q,vol,T),\
S0,\
strike_from_delta(put,S0,0.75,r,q,vol,T),\
strike_from_delta(put,S0,0.9,r,q,vol,T)]
strikes
Out[274]:
Calculate the break even smile for the first 10 paths
In [262]:
bevsmile = pd.DataFrame([pd.Series([opt.newton(bev, 0.3, args=(path,call,strike,r,q,T)) for strike in strikes], index=strikes) \
for path in paths[:20]])
Skew of break even vol
In [266]:
bevsmile.transpose().plot(legend=None)
Out[266]:
In [298]:
plbystrike = pd.DataFrame([pd.Series([calcportfolio(path,call,strike,r,q,vol_pricing,T)['portf'][len(path)-1] \
for strike in strikes], index=strikes) \
for path in paths[:5]])
In [299]:
plbystrike.plot()
Out[299]:
In [300]:
plbystrike
Out[300]:
In [21]:
import pandas.io as web
import pandas.io.data as webdata
from pandas.io.data import Options
import datetime
In [17]:
ticker="PETR4.SA"
start = datetime.datetime(2007, 1, 1)
end = datetime.datetime(2007,6,30)
bullQuotes = webdata.DataReader(ticker, 'yahoo', start, end)
start = datetime.datetime(2015, 7, 1)
end = datetime.datetime(2015,12,31)
bearQuotes = webdata.DataReader(ticker, 'yahoo', start, end)
In [19]:
plt.plot(bullQuotes['Close'])
Out[19]:
In [20]:
plt.plot(bearQuotes['Close'])
Out[20]: