Deriving Local Volatility from an Implied Volatility Surface

# Creating the local volatility surface from the IV surface


# Creating the local volatility surface from the IV surface












I have been using the dupire equation:

$$ \sigma_{LV} (K,T) = \frac{\sigma_{i m p}^2+2 \sigma_{i m p} T\left(\frac{\partial \sigma_{i m p}}{\partial T}+(r-q) K \frac{\partial \sigma_{i m p}}{\partial K}\right)}{\left(1 + d_1K\sqrt{T}\frac{\partial \sigma_{i m p}}{\partial K}\right)^2 + \sigma_{i m p}TK^2\left(\frac{\partial^2 \sigma_{i m p}}{\partial \sigma_{K}^2} - d_1\sqrt{T}\left(\frac{\partial \sigma_{i m p}}{\partial K}\right)^2\right)} $$

I have used linear interpolation for the IV data, but I sometimes get negative denominator values, which shouldn't happen. My data is: https://jmp.sh/s/dxUT9M9wE6j2R00Ak5Yg

(You will have to change the csv file name from `data - .csv` to `data.csv`)

which comes from this page: https://financetrainingcourse.com/education/2014/05/implied-and-local-volatility-surfaces-in-excel-final-steps/

This figures I get are

Which is very clearly different to what the website claims it is:

I have been stuck on this for a couple of days and am not sure where I have gone wrong. The interpolation technique for the IV should work since it's just linear interpolation...

My python code is:

```
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import interpolate

def implied_to_local_vol(implied_interp,eps,T,K,S,r,q):
    
    # The diffrence amount for numerical differentiation 
    up = 1+eps
    down = 1-eps
    
    sigma = implied_interp(T,K)
   
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    
    # 1st and 2nd derivatives using central difference
    iv_T = (implied_interp(T*up,K) -  implied_interp(T*down,K)) / ((up-down)*T)
    iv_K = (implied_interp(T,K*up) -  implied_interp(T,K*down)) / ((up-down)*K)
    iv_K2 = (implied_interp(T,K*up) - 2 * sigma +  implied_interp(T,K*down)) / (((up-down)*K)**2)

    numerator =  sigma**2 + 2*sigma *T * (iv_T + (r-q)*K*iv_K)

    denominator = (1 + d1*K*np.sqrt(T)*iv_K)**2 + sigma*T*(K**2)*( iv_K2 - d1*np.sqrt(T)*(iv_K**2))

    local_vol = np.sqrt(numerator/denominator)
    
    return local_vol

data = pd.read_csv('tools/data.csv',index_col=0)

# Importing data and cleaning. 
# T is shape(9,) and K is shape(13,)
T = (data.columns).values
T = T.astype(float)
K = data.index
K = K.astype(float)

imp_vol = data.to_numpy() # is shape (13,9) with strikes being the rows

# Turning strikes and expiries into a (13,9) meshgrid.
Tii,Kii = np.meshgrid(T,K)

"""
Stacking the meshgrids on top of each other
so the strikes and expiries become (117,) (because 13*9=117)
"""
Tii =Tii.flatten()
Kii = Kii.flatten()
imp_vol = imp_vol.flatten()

# Interpolaing and creating an interpolate object
interpolated_iv = interpolate.Rbf(Tii, Kii, imp_vol / 100, function='linear')

# Creating new data for the interpolated object that has 100 data points
ti = np.linspace(np.amin(T),np.amax(T),100)
ki = np.linspace(np.amin(K),np.amax(K),100)

# Constants
eps=0.01
r=0.01
S=15.7
q=0

local_vol = np.empty((100,100))
for j in range(100):
    for i in range(100):
        # Local vol is (K,T)
        local_vol[j,i] = implied_to_local_vol(interpolated_iv,eps,ti[i],ki[j],S,r,q) * 100

# Turning the data back into the meshgrid so we can plot 
ti,ki = np.meshgrid(ti,ki)
Tii,Kii = np.meshgrid(T,K)
imp_vol = imp_vol.reshape(13,9)

fig = plt.figure(figsize =(14, 9))
ax = plt.axes(projection ='3d')
ax.plot_surface(Tii, Kii, imp_vol)
ax.set_xlabel('T')
ax.set_ylabel('K')
ax.set_zlim([0,40])
ax.set_title('implied vol')
plt.show()

ax2 = plt.axes(projection ='3d')
ax2.plot_surface(ti, ki, local_vol)
ax2.set_xlabel('T')
ax2.set_ylabel('K')
ax2.set_zlim([0,40])
ax2.set_title('local vol')
plt.show()
```