I try to compute the local volatility in python in both formula, i.e. in terms of call price surface and total variance surface.
However, I get 2 different values. What did I do wrong?
import numpy as np
from scipy.stats import norm
today=44935
date1=44956
t1 = (date1-today)/365
a1=0.00403326433586633
b1=0.0393291646363571
sigma1=0.103678184440224
rho1=-0.147025189597259
m1=0.0154107775237231
date2=44957
t2 = (date2-today)/365
a2=0.00423386330856259
b2=0.0402496625734392
sigma2=0.106395807290502
rho2=-0.150124616359447
m2=0.0161557212111055
def rawSVI(logMoneyness,a,b,sigma,rho,m):
return a + b(rho(logMoneyness-m) + ((logMoneyness-m)2 + sigma2)**0.5)
def bsCall(s,k,r,t,vol):
d1 = (np.log(s/k) + (r+0.5vol2)t)/vol/ t0.5
d2 = d1 - vol*t0.5
return snorm.cdf(d1) - knp.exp(-rt)norm.cdf(d2)
w = rawSVI(0,a1,b1,sigma1,rho1,m1)
vol = np.sqrt(w/t1)
s = 1
r = 0
k = 0.65
f = snp.exp(rt1)
y = np.log(k/f)
#method 1
dk = 0.01
vol = np.sqrt(rawSVI(np.log(k/f),a1,b1,sigma1,rho1,m1)/t1)
call = bsCall(s,k,r,t1,vol)
k_plus = k*(1+dk)
vol_plus = np.sqrt(rawSVI(np.log(k_plus/f),a1,b1,sigma1,rho1,m1)/t1)
call_plus = bsCall(s,k_plus,r,t1,vol_plus)
k_minus = k*(1-dk)
vol_minus = np.sqrt(rawSVI(np.log(k_minus/f),a1,b1,sigma1,rho1,m1)/t1)
call_minus = bsCall(s,k_minus,r,t1,vol_minus)
call_dkdk = (call_plus - 2*call + call_minus)/(k_plus - k) / (k - k_minus)
vol_tPlus = np.sqrt(rawSVI(np.log(k/f),a2,b2,sigma2,rho2,m2)/t2)
call_tPlus = bsCall(s,k,r,t2,vol_tPlus)
call_dt = (call_tPlus - call)/(t2-t1)
localVol_1 = np.sqrt(2call_dt/k*2/call_dkdk)
print(localVol_1) #0.9224266360763037
#method 2
dy = 0.01
w = rawSVI(y,a1,b1,sigma1,rho1,m1)
w_plus = rawSVI(y(1+dy),a1,b1,sigma1,rho1,m1)
w_minus = rawSVI(y(1-dy),a1,b1,sigma1,rho1,m1)
w_dy = (w_plus - w_minus)/(2ydy)
w_dydy = (w_plus - 2w + w_minus)/(ydy)**2
w_tPlus = rawSVI(y,a2,b2,sigma2,rho2,m2)
w_dt = (w_tPlus - w)/(t2-t1)
localVol_2 = np.sqrt(w_dt/(1-y/ww_dy+0.25(-0.25-1/w+y2/w2)w_dy2+0.5w_dydy))
print(localVol_2) #0.8991455521477574
