Calibration Heston Model via Python

Semi-Closed-Form of Pricing Formula

import numpy as np

import scipy.integrate as spi

def Heston_Call_Value_Int(kappa, theta, sigma, rho, v0, r, T, s0, K):
a = s0*Heston_P_Value(kappa, theta, sigma, rho, v0, r, T, s0, K, 1)
b = K*np.exp(-r*T)*Heston_P_Value(kappa, theta, sigma, rho, v0, r, T, s0, K, 2)
#print (a,b)
return a-b

def Heston_P_Value(kappa, theta, sigma, rho, v0, r, T, s0, K, type):
ifun = lambda phi: Int_Function_1(phi, kappa, theta, sigma, rho, v0, r, T, s0, K, type)
return 0.5 + (1/np.pi)*spi.quad(ifun, 0, 100)[0]

def Int_Function_1(phi, kappa, theta, sigma, rho, v0, r, T, s0, K, type):
temp = (np.exp(-1*1j*phi*np.log(K))*Int_Function_2(phi, kappa, theta, sigma, rho, v0, r, T, s0, K, type))
return temp

def Int_Function_2(phi, kappa, theta, sigma, rho, v0, r, T, s0, K, type):
if type==1:
u = 0.5
b = kappa - rho*sigma
else:
u = -0.5
b = kappa
a = kappa*theta
x = np.log(s0)
d = np.sqrt((rho*sigma*phi*1j-b)**2 - sigma**2*(2*u*phi*1j-phi**2))
g = (b-rho*sigma*phi*1j+d)/(b-rho*sigma*phi*1j-d)
D = r*phi*1j*T + (a/sigma**2)*((b-rho*sigma*phi*1j+d)*T - 2*np.log((1-g*np.exp(d*T))/(1-g)))
E = ((b-rho*sigma*phi*1j+d)/sigma**2)*(1-np.exp(d*T))/(1-g*np.exp(d*T))

return np.exp(D + E*v0 + 1j*phi*x)

Define object functions, minimize pricing errors
import numpy as np
from A_Heston_Pricing_Integral import *

def Error_Function(p0, vK, price, v0, s0, r, T):
kappa, theta, sigma, rho = p0
price_e = np.zeros(len(vK), np.float)

for i in range(len(vK)):
price_e[i] = Heston_Call_Value_Int(kappa, theta, sigma, rho, v0, r, T, s0, vK[i])
RMSE = np.sqrt(np.sum(price-price_e)**2/len(price))
print np.vstack((price_e,price)).T
return RMSE

Use optimization method to locate parameters

import numpy as np
import scipy as sp
import scipy.optimize as spopt
import time as tm
from E_Error_Function_Definition import *

def Heston_Option_Values(vK, kappa, theta, sigma, rho, v0, r, T, s0):
HOV = np.zeros(len(vK), dtype=np.float)
for i in range(len(vK)):
HOV[i] = Heston_Call_Value_Int(kappa, theta, sigma, rho, v0, r, T, s0, vK[i])
return HOV

# generate theoretical option prices,
# for a vector of strike price
# for known parameters
kappa = 4.0
theta = 0.04
sigma = 0.1
rho = -0.5
v0 = 0.04
r = 0.01
T = 0.5
s0 = 1.0
K = np.arange(1.2, 0.795, -0.05)
price = Heston_Option_Values(K, kappa, theta, sigma, rho, v0, r, T, s0)

t1 = tm.time()

p0 = np.array([3.0, 0.05, 0.15,-0.7]) # initial guess for parameters
opt = spopt.fmin(Error_Function, p0, maxiter=75, maxfun=75, args=(K, price, v0, s0, r, T))
t2 = tm.time()

print "Duration in seconds", t2 - t1
print opt