Difference of polynomial interpolation for volatility smile

Question

I am using 5 volatility points to build a volatility smile : put 10D, put 25D, ATMF, call 25D and call 10D. I have thus 5 pairs of data : (Delta, Vol) let's say for example (10;5.75) ; (25; 5.50) ; (50;5.25) ; (75;5.60) ; (90;5.70).

I am doing a simple lagrange interpolation as below :

x = np.array([10, 25, 50, 75, 90]) #LIST OF DELTA
y = np.array([])

for index, elt in enumerate(vols[currency][tenor]):
    y=np.append(y,vols[currency][tenor][elt])

interpolation = lagrange(x, y) #TO KEEP
xpol = np.linspace(x[0], x[4], 90)
ypol_lag = interpolation(xpol)
plt.scatter(x, y, marker='s', c='r')
plt.plot(xpol, ypol_lag, "b")
plt.show()

Which returns a nice smile, very realistic given the Market data I can observe in Bloomberg :

Then, I am transforming each delta in strike given the formula provided by Wystup, 2010 :

Therefore I have all my pairs (Strike, vol), let's say for example : (10;1.1650) ; (25; 1.1710) ; (50;1.1800) ; (75;1.1840) ; (90;1.1950) (non real numbers).

If I plot these five numbers and interpolate them using the same function as I did for the deltas :

x = np.array([K10p, K25p, K50, K25c, K10c])
y = np.array([])

for index, elt in enumerate(vols[currency][tenor]):
    y=np.append(y,vols[currency][tenor][elt])

interpolation = lagrange(x, y)

xpol = np.linspace(x[0], x[4], 90)
ypol_lag = interpolation(xpol)
plt.scatter(x, y, marker='s', c='r')
plt.plot(xpol, ypol_lag, "r")
plt.show()

Surprisingly, the interpolation curb is different and does not look like a smile :

I have no idea how to get a similar interpolation with the strikes...

Thank you.

Answer 1

This is not really an answer but it's too long for a comment.

The lagrange / cublic spline interpolation is very sensitive to the input data, given slightly different input data it can produce vastly different interpolation. It also affects not only neighbouring points but also different points on the curve. Consider this example.

x = np.array([10, 25, 50, 75, 90]) #LIST OF DELTA
y = np.array([5.60, 5.40, 5.25, 5.28, 5.35])

interpolation = lagrange(x, y)
xpol = np.linspace(x[0], x[4], 90)
ypol_lag = interpolation(xpol)
plt.scatter(x, y, marker='s', c='r')
plt.plot(xpol, ypol_lag, "b")
plt.show()

Now let's change only the 2nd point:

x = np.array([10, 35, 50, 75, 90]) #change from 25 to 35
y = np.array([5.60, 5.40, 5.25, 5.28, 5.35])

the result:

We have changed only 2nd point but the whole curve changed (also interpolation between 4th and 5th point changed!). The problem with cublic splines is that it finds the parameters that solve linear system of equations, and the parameters that are found affect whole curve. Cubic spline is "unstable" in that manner. When you are transforming each delta in strike, you are changing the input in non linear fashion which results in change in the whole curve.

If you want to use cubic spline interpolation anyway, the way around this problem is to firstly interpolate the data before transformation (so on deltas), then apply transformation on interpolated data. You should get similiar shape before and after transformation. To be more precise, you have 5 deltas, interpolate the smile for example for 20 deltas. Then transform each point of 20 in strike by the formula you used. Then interpolate on 20 strikes and plot the result.

EDIT: I have tried this method above for lagrange interpolation but I got numerically unstable output. According to numpy documentation this function is unstable when we have 20 data points.

Warning: This implementation is numerically unstable. Do not expect to be able to use more than about 20 points even if they are chosen optimally.

Therefore I switched to CubicSpline function which produces the same looking graphs as lagrange, but is numerically stable.

When you run the code below:

x = np.array([10, 25, 50, 75, 90])
y = np.array([5.60, 5.40, 5.25, 5.28, 5.35])

# We first interpolate data on deltas to create 20 points
interpolation = CubicSpline(x, y)
xpol = np.linspace(x[0], x[4], 20)
ypol_lag_org = interpolation(xpol)

# We then transform delta in strike by the equation for all points interpolated, we now have 20 points
K_list = []
for a, b in zip(xpol,ypol_lag_org):
    K = np.exp(norm.ppf(a/100)*(b/100)+0.5*((b/100)**2))
    K_list.append(K)
    
K_list = np.array(K_list) # this is vector of 20 numbers

# And then we interpolate the rest of the points and produce a graph
interpolation = CubicSpline(K_list, ypol_lag_org)
xpol = np.linspace(K_list[0], K_list[-1], 90)
ypol_lag = interpolation(xpol)
plt.scatter(K_list, ypol_lag_org, marker='s', c='r')
plt.plot(xpol, ypol_lag, "b")
plt.show()

you will get something that looks like smile but is in (strike,vol) space. The additional points enforce the stability of the curve.