Does Bakshi, Kapadia and Madan (2003) VIX building approach underestimate volatility?

Question

From a paper that shortly addresses an alternative approach to VIX-like index building:

To test this approach, I've built a fake book of B&S options with constant volatility equal to $\sigma=20\%$, $0\%$ risk free rate and $0\%$ cost of carry. I focused my attention on $T=15$ days maturity (roughly $0.0411$ years on $365$ basis).

Starting from the assumption that Bakshi, Kapadia and Madan (2003) approach is correct, one should get as a result a VIX-like object whose value is equal to $20\%$ if he applied the formula to the book above.

From partial moments to VIX-like value with fake $\sigma=20\%$ options:

H1  =   0,001647398
H2  =   -4,6183E-06
H3  =   7,30375E-06
mu  =   -0,000823234
VAR =   0,040070196
VIX =   sqrt(VAR) = 20%

My VBA code so you can easily reproduce my results:

Public Function VAR_T(F As Double, K_C As Range, K_P As Range, C As Range, P As Range, rf As Double, t As Double) As Double

    ' F         Underlying forward value
    ' K_C       Range of Call options' strike prices
    ' K_P       Range of Put options' strike prices
    ' C         Range of Call options' prices (same extent of K_C is mandatory)
    ' P         Range of Put options' prices (same extent of K_P is mandatory)
    ' rf        Risk free rate (as instance of EUR, EONIA spot is suggested)
    ' t         Time to expiry in years on 365 basis

    Dim kH1 As Double
    Dim kmu_t As Double

    kH1 = H1_T(F, K_C, K_P, C, P)
    kmu_t = mu_T(F, K_C, K_P, C, P, rf, t)
    VAR_T = (Exp(rf * t) * kH1 - (kmu_t) ^ 2) / t

    ' Squared root of VAR_T is the underlying volatility on expiry date equal to today + t

End Function

Public Function mu_T(F As Double, K_C As Range, K_P As Range, C As Range, P As Range, rf As Double, t As Double) As Double

    Dim kH1 As Double
    Dim kH2 As Double
    Dim kH3 As Double

    kH1 = H1_T(F, K_C, K_P, C, P)
    kH2 = H2_T(F, K_C, K_P, C, P)
    kH3 = H3_T(F, K_C, K_P, C, P)

    mu_T = Exp(rf * t) - 1 - Exp(rf * t) / 2 * kH1 - Exp(rf * t) / 6 * kH2 - Exp(rf * t) / 24 * kH3

End Function

Public Function H1_T(F As Double, K_C As Range, K_P As Range, C As Range, P As Range) As Double

    Dim n_C As Integer
    Dim n_P As Integer
    Dim partial_moment_C()
    Dim partial_moment_P()

    n_C = K_C.Count
    n_P = K_P.Count

    ReDim partial_moment_C(n_C)
    ReDim partial_moment_P(n_P)

    For i = 1 To n_C
        partial_moment_C(i) = 2 * (1 - Log(K_C(i) / F)) / (K_C(i) ^ 2) * C(i)
    Next

    For i = 1 To n_P
        partial_moment_P(i) = 2 * (1 + Log(F / K_P(i))) / (K_P(i) ^ 2) * P(i)
    Next

    H1_T = Application.Sum(partial_moment_C) + Application.Sum(partial_moment_P)

End Function

Public Function H2_T(F As Double, K_C As Range, K_P As Range, C As Range, P As Range) As Double

    Dim n_C As Integer
    Dim n_P As Integer
    Dim partial_moment_C()
    Dim partial_moment_P()

    n_C = K_C.Count
    n_P = K_P.Count

    ReDim partial_moment_C(n_C)
    ReDim partial_moment_P(n_P)

    For i = 1 To n_C
        partial_moment_C(i) = (6 * Log(K_C(i) / F) - 3 * (Log(K_C(i) / F)) ^ 2) / (K_C(i) ^ 2) * C(i)
    Next

    For i = 1 To n_P
        partial_moment_P(i) = (6 * Log(F / K_P(i)) + 3 * (Log(F / K_P(i))) ^ 2) / (K_P(i) ^ 2) * P(i)
    Next

    H2_T = Application.Sum(partial_moment_C) - Application.Sum(partial_moment_P)

End Function

Public Function H3_T(F As Double, K_C As Range, K_P As Range, C As Range, P As Range) As Double

    Dim n_C As Integer
    Dim n_P As Integer
    Dim partial_moment_C()
    Dim partial_moment_P()

    n_C = K_C.Count
    n_P = K_P.Count

    ReDim partial_moment_C(n_C)
    ReDim partial_moment_P(n_P)

    For i = 1 To n_C
        partial_moment_C(i) = (12 * (Log(K_C(i) / F)) ^ 2 - 4 * (Log(K_C(i) / F)) ^ 3) / (K_C(i) ^ 2) * C(i)
    Next

    For i = 1 To n_P
        partial_moment_P(i) = (12 * (Log(F / K_P(i))) ^ 2 + 4 * (Log(F / K_P(i))) ^ 3) / (K_P(i) ^ 2) * P(i)
    Next

    H3_T = Application.Sum(partial_moment_C) + Application.Sum(partial_moment_P)

End Function

If you use the functions above to calculate the VIX-like value of the fake options, you will get $20\%$ as expected; but, if you build fake options using a huge implied volatility, something like $\sigma=80\%$, my code returns an implied model free volatility of... $62\%$!

Possibilities:

1. I am wrong with my understanding of the aforementioned formulas;
2. I am wrong with the implementation of the aforementioned formulas in VBA code;
3. Bakshi, Kapadia and Madan (2003) approach underestimate volatility, which sounds unlikely... what am I not taking into account?

Answer 1

In your implementation, you are approximating continuous integrals over the strike domain by Riemann sums. This introduces an error.

More specifically, for a fixed time to expiry $\tau$, you're integrating (scaled) OTM price curves (see equations $(3)-(4)-(5)$). As volatility increases, these curves will stretch over wider and wider spatial domains.

When using a fixed strike range for approximating the integrals, the truncation error will therefore mechanically increase as you will less and less be able to capture what happens in the `tails'.