BKM risk neutral moments in python

Question

I am trying to compute the BKM implied moments (Bakshi, Kapadia and Madan 2003) in python by following this paper:

Neumann, Skiadopoulos: Predictable Dynamics in Higher Order Risk-Neutral Moments: Evidence from the S&P 500 Options, Journal of Financial and Quantitative Analysis (JFQA), 2013, p947-977 link

which in page 7 and 8, describe the integrals for variance, skewness and kurtosis by equations (8), (9) and (10).

In page 8 part B. Empirical implementation, the authors go on about their methodology of extracting the implied moments. What I don't understand is:

a) Why they have to interpolate across the implied volatilities as a function of delta (by creating an artificial 'delta grid'), and just to convert these deltas back to the respective strikes afterwards via BS? They say the integrals require 'a continuum of OTM call and put options across strikes', and that we only observe discrete strike prices, hence the reason for using 'delta space'. Why does interpolating in delta space make it any more 'continuous'?

b) From page 9 of the paper: Then, we compute the constant maturity moments [equations (5), (6), (7)] by evaluating the integrals in formulae (8), (9), and (10) using trapezoidal approximation. How do you evaluate those definite integrals in python? I am stuck staring at those 3 integrals not knowing how to continue in python.

Example (part of eq (8)):

Answer 1

I think the reason to interpolate using a delta grid is twofold.

First, note that the delta of an option (or lets talk about a given strike) gives you some information about how OTM/ITM the option is. Using strikes without checking the deltas could imply that for some maturity you could choose a nice strikes-grid to sample, but however for another set of variable choices that grid would not be so useful or representative, therefore leading to numerical errors in your code.

This leads to the second reason, which is related to the first one. From the coding point of view, implementing a sampling routine can be made much more robust using a delta-grid as compared to a strike-based choice. Think about it this way: no matter what are your variable choices, the order of magnitude of the underlying price (which could range from pennies to thousands of dollars for a stock, or could be 1.10-1.20 for some FX, while being 500 for another, etc). If you implement a routine using deltas, then the interpolation mechanism/routine would not care about this issue, as it's defined in terms of delta no matter the underlying, which is actually a quantity that means something (as I said, it gives you an idea on how much an option is ITM/OTM). Therefore, your code does not need you to hard-code a solution on how to interpolate depending on the underlying level, etc.

Regarding question b), I suggest you look for some library that allows you to compute those integrals. Unfortunately, I am not as familiarized with those as to suggest one. I think scipy should have some routines, maybe you could take a look at them. Edit: I did not see Kermittfrog comment on this. Thant links points directly to it, so I recommend you go through it.