-
Notifications
You must be signed in to change notification settings - Fork 2
/
Heston Options Pricer.py
166 lines (132 loc) · 5.63 KB
/
Heston Options Pricer.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
import numpy as np
import matplotlib.pyplot as plt
from py_vollib_vectorized import vectorized_implied_volatility as impvol
#Parameters for generating the correlated Brownian Motions
steps=252 #252 trading days in a year
rho=-0.65
mu=np.array([0, 0])
sig=np.array([[1, rho],[rho, 1]])
#Parameters for simulating the volatility process
nu0=0.2 #Initial volatility
theta=0.2 #Long-term mean of the volatility
kappa=3 #Mean reversion rate
ksi=0.5 #Volatility of the volatility process
T=1 #Time until expiry of the option. Here we take T=1 year, and each step corresponds to 1 trading day
dt=T/steps #Split into equal time steps
#Parameters for simulating the underlying stock process
S0=100.0 #Initial stock price
r=0.001 #Risk-free rate
M=1000 #Number of simulations
def heston_model(steps, r, rho, S0, nu0, theta, kappa, ksi, T, M):
dt=T/steps #Split into equal time steps
#Parameters for generating the correlated Brownian Motions
mu=np.array([0, 0])
sig=np.array([[1, rho],[rho, 1]])
#M pairs of Brownian Motion increments drawn from the bivariate normal distribution
bminc=np.random.multivariate_normal(mu, sig, (steps, M))
Wnuinc=np.squeeze(bminc[:,:,0]) #Increments of the BM in the volatility process.
WSinc=np.squeeze(bminc[:,:,1]) #Increments of the BM in the stock price process.
#Plot a pair of correlated BMs
plt.figure(figsize=(10, 6))
plt.plot(Wnuinc[:,0].cumsum(axis=0), label=r"$W^\nu_t$", linestyle='-', color='blue')
plt.plot(WSinc[:,0].cumsum(axis=0), label=r"$W^S_t$", linestyle='-', color='orange')
plt.xlabel("Time Steps")
plt.ylabel("Brownian Motion")
plt.title(rf"Correlated Brownian Motion, $\rho$={rho}")
plt.legend(prop={'size': 15})
plt.grid(True)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
#Simulating the volatility process & the stock price process
St=np.full(shape=(steps,M), fill_value=S0)
nut=np.full(shape=(steps,M), fill_value=nu0)
for i in range(1,steps):
nut[i]=np.abs(nut[i-1] + kappa*(theta-nut[i-1])*dt + ksi*np.sqrt(nut[i-1]*dt)*Wnuinc[i-1,:])
St[i]=St[i-1]+r*St[i-1]*dt+np.sqrt(nut[i-1]*dt)*St[i-1]*WSinc[i-1,:]
return St, nut
t=np.linspace(0, T, steps)
St, nut = heston_model(steps, r, rho, S0, nu0, theta, kappa, ksi, T, M)
#Plot the volatility process along with the long-term mean
plt.figure(figsize=(10, 6))
plt.plot(t, nut)
plt.axhline(y=theta, color='black', linestyle="--", label=r"Long-term mean, $\theta$")
plt.legend()
plt.xlabel("Time")
plt.ylabel(r"Volatility $\nu_t$")
plt.title("Heston Volatility Paths")
plt.grid(True)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
#Plot the stock process
plt.figure(figsize=(10, 6))
plt.plot(t, St, label="Stock Price")
plt.xlabel("Time")
plt.ylabel(r"Price $S(t)$")
plt.title("Heston Underlying Paths")
plt.grid(True)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
ST=St[-1,:] #An array of the underlying price at maturity for each path
K=np.arange(50,200,2) #An array of strikes
#Option price under the risk-neutral measure is the expectation value
#of the payoff function, discounted to today
callprices=np.array([np.mean(np.exp(-r*T)*np.maximum(ST-k,0)) for k in K])
putprices=np.array([np.mean(np.exp(-r*T)*np.maximum(k-ST,0)) for k in K])
#Calculate implied volatilities under Black-Scholes & plot them as a function of strike price ,
#to show that the Heston model captures volatility smiles/smirks/skews.
callimpvols=impvol(callprices, S0, K, T, r, flag="c", q=0, return_as="numpy")
putimpvols=impvol(putprices, S0, K, T, r, flag="p", q=0, return_as="numpy")
plt.figure(figsize=(10, 6))
plt.plot(K, callimpvols, label="Calls")
plt.plot(K, putimpvols, label="Puts")
plt.title("Implied Volatility as a Function of Strike Price")
plt.legend()
plt.xlabel("Strike Price")
plt.ylabel("Implied Volatility")
plt.grid(True)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
#Compute implied volatilities as a function of time to maturity
Ts=np.arange(1,3,1/12)
volsurfcall=np.array([impvol(callprices, S0, K, t, r, flag="c", q=0, return_as="numpy") for t in Ts])
volsurfput=np.array([impvol(putprices, S0, K, t, r, flag="p", q=0, return_as="numpy") for t in Ts])
#Volatility surface for calls
X, Y = np.meshgrid(Ts, K, indexing="ij")
Z=volsurfcall
fig = plt.figure(figsize=(8,8), dpi=150)
ax = fig.add_subplot(111, projection='3d')
volsurface = ax.plot_surface(X, Y, Z, cmap='viridis', antialiased=True)
ax.set_xlabel('Maturity (Years)')
ax.set_ylabel('Strike Price')
ax.zaxis.set_rotate_label(False)
ax.set_zlabel('Implied Volatility', rotation=90)
ax.xaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.yaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.zaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.view_init(elev=20, azim=65)
plt.title("Heston Volatility Surface, Calls", y=1)
plt.show()
#Volatility surface for puts
X, Y = np.meshgrid(Ts, K, indexing="ij")
Z=volsurfput
fig = plt.figure(figsize=(8,8), dpi=150)
ax = fig.add_subplot(111, projection='3d')
volsurface = ax.plot_surface(X, Y, Z, cmap='viridis', antialiased=True)
ax.set_xlabel('Maturity (Years)')
ax.set_ylabel('Strike Price')
ax.zaxis.set_rotate_label(False)
ax.set_zlabel('Implied Volatility', rotation=90)
ax.xaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.yaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.zaxis._axinfo['grid'].update(color='grey', linestyle='dotted')
ax.view_init(elev=20, azim=65)
plt.title("Heston Volatility Surface, Puts", y=1)
plt.show()