Bakshi, Kapadia, and Madan (2003) risk-neutral moment estimators: A Gram–Charlier density approach

This paper is a sequel to Aschakulporn and Zhang (J Futures Mark 42(3):365–388, 2022). The errors of the Bakshi et al. (Rev Financ Stud 16(1):101–143, 2003) risk-neutral moment estimators is studied using the Gram–Charlier density—with the skewness and excess kurtosis specified. To obtain skewness with (absolute) errors less than 10-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-3}$$\end{document}, the range of strikes (Kmin,Kmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{\min }, K_{\max }$$\end{document}) must contain at least 3/4 to 4/3 of the forward price and have a step size (ΔK\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta K$$\end{document}) of no more than 0.1% of the forward price. The range of strikes and step size corresponds to truncation and discretization errors, respectively. This is consistent to Aschakulporn and Zhang (2022) for non-volatile market periods.


Introduction
This paper examines the  (BKM) risk-neutral estimators with explicit skewness and kurtosis as both the input and benchmark and finds the condition of strike prices required to bound the errors of skewness. The BKM estimators 1 3 have been widely used in the literature and are used by both academics and practitioners; for example, the BKM estimator for skewness is used to calculate the Chicago Board Options Exchange (CBOE) skewness (SKEW) index, 1 the forward-looking tail-risk or crash risk indicator of the S&P 500 index. Although the errors of the BKM estimators have been studied by Aschakulporn and Zhang (2022) (AZ), the relationship between the model's parameters and skewness is complicated. This paper analyzes the BKM risk-neutral estimators using simulated option prices generated by a closed-form option-pricing formula with skewness and kurtosis (derived using the Gram-Charlier density)-in doing so, the estimated values can be compared directly with their true values.
A literature review regarding the use and importance of the BKM estimators is presented in AZ. Four papers have examined BKM errors-AZ, Ammann and Feser (2019), Liu and van der Heijden (2016), and Lee and Yang (2015)-with AZ being the key piece of the literature. Prior to AZ, the most advanced model used was the Bates (1996) model-a combination of Heston's (1993) stochastic volatility and Merton's (1976) jump. AZ use the Duffie et al. (2000) double-jump affine jumpdiffusion model and derive the skewness analytically to get the closed-form solution; this removes the approximation aspect of previous literature. However, due to the complexity of the affine jump-diffusion model, the entire range of their parameters cannot be examined exhaustively, and the relationship between skewness and the model parameters is also lengthy.
Following AZ, the errors of BKM risk-neutral moment estimators caused by truncation and discretization are tested using simulated option prices. The BKM estimators are derived in a setting that has options with a continuum of strikes over an infinite domain; this, however, does not match reality, as options data have discrete strikes over a finite domain. These discrete strikes cause discretization errors, and a finite domain causes truncation errors. Having complete control over options allows the errors of the BKM estimators to be analyzed without external influences in a controlled manner. Using simulated option prices, the BKM estimators can be applied and the estimates can be compared with the true moments calculated from the model originally used to create the options. Similar to AZ, the true moments need not be found from simulations or approximated. The true moments are instead found analytically. For AZ, their formula is quite lengthy; whereas for the Gram-Charlier density, the true moments are specified, not obtained, as the input parameters are the moments themselves. Following Zhang and Xiang (2008), option prices are generated using an option pricing formula with skewness and excess kurtosis (hereinafter, kurtosis). The Gram-Charlier density is used as it provides a way to introduce higher-order moments, namely skewness and kurtosis, to the standard normal distribution. This method has been used by many, including Jarrow and Rudd (1982), Corrado and Su (1996), Longstaff (1995), Backus et al. (1997), Kochard (1999), and Corrado (2007). Using AZ's framework, the effects of constant and linear extrapolation, and linear, cubic spline, and kernel interpolation on the BKM estimators are similarly tested. In addition, this paper extends the analysis of AZ to examine the volatility, cubic, and quartic contracts used to calculate the BKM estimators in depth. The error of the approximate convexity term used to derive the BKM estimators is also examined.
The remainder of this paper is organized as follows. Section 2 presents an option pricing model with skewness and kurtosis. Section 3 presents the method used to create options prices and the method used to quantify and analyze the errors and convergence of the BKM method. Section 4 describes the data. Section 5 provides the numerical results, and Sect. 6 concludes. The appendix gives the details of key derivations, including part A which briefly presents the BKM risk-neutral moment estimators.

Option pricing formula with skewness and kurtosis
The return of the Black-Scholes model is normally distributed. The return has no skewness or excess kurtosis. The Gram-Charlier series is a statistical method that is commonly used to introduce skewness and kurtosis (and higher-order moments) to the normal distribution. 2 Combining the two, a parsimonious option-pricing formula with skewness and kurtosis can be derived. The two additional parameters, the skewness and kurtosis, are the parameters of interest for this paper. As there are only two additional parameters, examining their effects on errors comprehensively is possible.

Option pricing using the Gram-Charlier density
This paper follows Zhang and Xiang (2008) to derive the option-pricing formula. They first specify the underlying stock price at maturity, S T , to be where F T t , , , c , and y is the forward price, standard deviation, time to maturity ( T − t ), convexity adjustment term, and a random variable, respectively. Using the Gram-Charlier density, if y has probability density where n(y) = 1 √ 2 e − y 2 2 , then y will have a mean of zero, variance of one, skewness equal to 1 and an excess kurtosis of 2 . The forward price, F T t , is related to the current stock price, S t , by F T t = S t e (r−q) , where r is the risk-free rate and q is the continuous dividend rate. The convexity adjustment term is required to keep this model (2) f (y) = n(y) − 1 3! d 3 n(y) dy 3 + 2 4! d 4 n(y) dy 4 , arbitrage-free by ensuring that the stock price satisfies the martingale condition ( F T t = E Q t S T ) in the risk-neutral world and was found to be Using the risk-neutral valuation formula of Harrison and Kreps (1979) and Harrison and Pliska (1981), Proposition 1, the European call option pricing formula with skewness and kurtosis can be derived.

Proposition 1 The price of a European call option with skewness and kurtosis is
where The European put option price can be found using the put-call parity: Remark 1 If the skewness, 1 , and excess kurtosis, 2 , are both zero, then the convexity adjustment term c will also be zero and Eq. (4) will be equivalent to Black and Scholes (1973).
Remark 2 There is a minor error in Zhang and Xiang (2008), which, once fixed, allows the final result to be expressed more concisely. In Appendix B of Zhang and Xiang (2008), the proof of proposition 2, the call price is given as: where (5) c t − p t = F T t e −r − Ke −r .
, respectively-using their notation. They also extend from using Hermite polynomials to logistic polynomials (another orthogonal polynomial) and show that orthogonal polynomials can be used to span option payoffs. However, Hermite polynomials tend to diverge unlike logistic polynomials.
Remark 8 An advantage of using the Gram-Charlier density is that it can be easily extended to have an arbitrary number of additional moments/cumulants while satisfying the martingale condition.
For the purposes of this paper, the Gram-Charlier density (Eq. (2)), including only the skewness and kurtosis term is used. Figure 1 shows the region of skewness and kurtosis, which ensures a valid density function. The derivation of this region and an alternate formulation of the density is shown in the appendix, part D.
This formulation (Eq. 4) is ideal for testing the BKM estimators. The stock price model satisfies the no-arbitrage principle, unlike Corrado and Su (1996). The martingale restriction is not approximated like Backus et al. (1997). The expression is clearer and more elegant compared to Kochard (1999), Schlögl (2013), and Ki et al. (2005). The relationship between the model's skewness and excess kurtosis and inputs are clear as they are the same, unlike Corrado (2007) and Chateau and Dufresne (2017). Gram-Charlier valid region. This figure shows the region (shaded area) in which pairs of skewness and kurtosis values will yield a valid probability density function using the Gram-Charlier density. The curve represents the skewness and kurtosis required to make the density function vanish. Point A corresponds to the origin and points B to O correspond to various values of x between −∞ and ∞ 1 3  risk-neutral moment…

Methodology
To analyze the errors of BKM risk-neutral moment estimators, a controlled environment with clear benchmarks and parameters would be ideal. Part A of the appendix presents a brief derivation of the BKM estimator. The BKM estimator requires the use of option prices which can be provided by simulated option prices generated based on the Gram-Charlier density. The Gram-Charlier density is ideal, as it provides a malleable framework with skewness and kurtosis that are explicitly specified as parameters. These parameters are used as direct benchmarks against the moments calculated by the BKM estimator.
The procedure to test the BKM estimator is to (1) generate simulated option prices (Sect. 3.1), (2) apply the BKM estimator (Part A of the appendix), and (3) compare the estimates with their true values (Sect. 3.2). In addition to this, various interpolation and extrapolation techniques could be tested by introducing these techniques before applying the BKM estimator. This procedure follows the error analysis of AZ with the addition of a more comprehensive and in-depth analysis. Instead of testing two data points in two different market periods, the entire range of parameters (2515 points) is tested using the Gram-Charlier approach. Each contract used to calculate the BKM estimators are examined individually and with a smaller step size in Sect. 5.4. The rare disaster concern index (RIX) from Gao et al. (2018) and simple variance swaps (SVIX) from Martin (2016) are also tested.

Simulated option prices using the Gram-Charlier density
The options prices generated by Eq. (4) will be for an underlying asset with a specified return distribution. These values can be generated with specific sets of strike prices to test various errors. As the skewness and kurtosis is bounded by the valid region of the Gram-Charlier density (shown in Fig. 1), 2515 different pairs of skewness and kurtosis values within the valid region were chosen for analysis. These points are shown in Fig. 2.
Using the Gram-Charlier approach, the BKM estimator errors can be found directly. Similarly, other estimators can be tested with this approach. The CBOE VIX, rare disaster concern index, RIX, from Gao et al. (2018) and simple variance swaps, SVIX, from Martin (2016) derived under the Gram-Charlier model is shown in the appendix, part E.

Errors from truncation and discretization
Following AZ, the truncation and discretization errors are the two sources of errors that were tested. The same definitions are used to allow for consistency. These are presented in part B of the appendix for convenience. For calculations, the maximum and minimum strike prices are rounded to the nearest dollar and the boundary 1 3 controlling factor, a, is varied between 0.05 and 0.95 in 0.01 intervals. 3 The step size, ΔK , has been chosen to vary from 1 to 50, in increments of 1. 4 For each combination of a and ΔK , of which there are 4550 (91 × 50), the mesh of points (pairs of skewness and kurtosis) within the valid Gram-Charlier region have been used to generate simulated option prices (as described in Sect. 3.1). From this, the discretized BKM method is applied to each set of options to find the riskneutral moments. The estimation error is defined as following AZ.
(6) Estimation Error ∶= Estimated Moment − True Moment, The error for each point in the valid Gram-Charlier region for a specific combination of a and ΔK is averaged to find the average error for the whole valid region. The maximum of the absolute value of the error over the valid region was also recorded.
The specification of each parameter is shown in Table 1, and the skewness-kurtosis points used in the valid Gram-Charlier region are shown in Fig. 2.
An example of the error calculation for skewness and kurtosis done for the points from Fig. 2 are shown in Fig. 3. From this, both the average of the errors and the maximum absolute error for each region are plotted at their respective boundary controlling factor and step size.

Errors from the approximated convexity adjustment term
The approximation made by BKM (Eq. A2) can be analyzed analytically by entering the Black and Scholes (1973)  . Clearly, the returns are not the same. This additional term corresponds to the convexity adjustment term c as calculated by the BKM method. This is found by comparing the expected log returns of the BKM method and Gram-Charlier density return: Using Eq. (8) in place of the BKM , the errors caused by BKM's approximation can be calculated.

Error reduction by using interpolation and extrapolation
Following AZ, various interpolation and extrapolation techniques are applied to test whether the errors of the BKM estimator can be reduced. Constant and linear extrapolation as well as linear, cubic spline, and Gaussian kernel interpolation were applied to the implied volatility curve. What differs from AZ is that interpolation is done so that the largest step size between each strike is no larger than 0.1% of (7) The details are shown in the appendix, part F. the forward price (with a minimum of $1), and the extrapolation is done so that the domain contains strike prices between 3/4 and 4/3 of the forward price. In AZ, interpolation was done to 0.05% of the forward price and the boundary controlling factor was 1/4-resulting in smaller step sizes over a larger domain.

Data
As this paper uses simulated option prices to test for estimation errors, the data required are minimal and can be set arbitrarily. However, to direct the results towards US markets, specifically the S&P 500, the risk-free rate, volatility, and time to maturity have been set to 2.4%, 6 0.20, and one month, respectively. These values partially reflect the current market conditions in April 2019. The current market volatility based on the CBOE VIX 7 is oddly low, so the standard deviation has been chosen to be 0.20. A one month time to maturity has been chosen as these tend to be the most liquid (shortest term contracts).

BKM estimation
The mean errors of the BKM standard deviation, skewness, and kurtosis estimators are shown in Table 2. The mean error is the mean of the estimation error (Eq. (6)) calculated over each valid Gram-Charlier region. These errors show that, as expected, the smaller the boundary controlling factor, a, the smaller the approximation error. Similarly, a smaller step size, ΔK , corresponds to a smaller discretization error. As these two approach zero, the discretization and truncation of the BKM method tend towards an integral over an infinite range. For the average of skewness errors to be less than 10 −3 , the boundary controlling factor and step size, as shown in Table 2, must be less than 0.75 and $20, respectively. A more robust way to limit errors to be below 10 −3 is done by using the maximum value of the absolute error rather than the mean. This is shown in Table 3. The constraint is stronger and is reflected in the allowable boundary controlling factor and step size. The boundary controlling factor and step size to ensure that the errors are bounded by 10 −3 must also be bounded by 0.75 and $2 , respectively. The step size is dependent on the forward price; to generalize this result, the step size must be less than 0.1% of the forward price. Comparing the errors of the second to fourth moments shows that the BKM estimator, for these values, is more accurate for lower-order moment estimators. These results are consistent with AZ for a non-volatile market period. During volatile market periods, AZ show that the errors of the BKM skewness estimator are accurate when the boundary controlling factor is 0.25 and the step size is 0.05% of the forward price ($1). A visual form of Table 2 is presented in Fig. 4. This figure shows that the errors caused by truncation of the integral causes less predictable behaviour, whereas changes in the step size do not seem to affect the error as unpredictably. Although the truncation error is less predictable, when the boundary controlling factor is below 0.75, the behaviour is stable. Sensitivity of the Average Error for Skewness. Using the Gram-Charlier density, simulated option prices with known standard deviation , skewness 1 , and kurtosis 2 were created for a forward price F T t , time to maturity , risk-free rate r, and standard deviation , have been arbitrarily set to $2000, one month, 2.4%, and 0.20, respectively. In total, 2515 different pairs of skewness and kurtosis (within the Gram-Charlier region) were used. The errors are shown for a varying degree of fineness of strike discretization ( ΔK ∈ [1, 50] ). The range of strikes used is also varied symmetrically by adjusting  boundary of 10 −3 is shown in red. The same has been done for standard deviation and is shown in Fig. 7. Figures 6 and 8 show the same errors but with respect to the boundary controlling factors for step sizes of 1, 10, and 25.    The skewness errors can easily be restricted below 10 −3 by fixing ΔK = $2 and a ≤ 0.75 . Kurtosis, however, has much larger errors.
The standard deviation, skewness, and kurtosis using the BKM method have been calculated for a skewness and kurtosis of −1 and 2.5, respectively. This has been done for a boundary controlling factor value of 0.05, 0.40, and 0.75 and step size of 1 to 50. This is shown in Fig. 9. Figure 10 shows the same errors with respect to the boundary controlling factors.
To examine the independence of the skewness estimator from kurtosis, the skewness is estimated when the skewness, the boundary controlling factor, and step size are set to −1 , 0.25, and 1, respectively, whilst adjusting the kurtosis within the valid Gram-Charlier region. The same has been done for kurtosis when it is set to 2.5. As neither the true skewness nor the true kurtosis has been chosen to be zero, the relative error for each estimator can, therefore, be calculated. The (absolute) relative error is defined as Relative Error = | | |   is the true value. The maximum(minimum) relative error for skewness and kurtosis were found to be 4.476 × 10 −04 (4.476 × 10 −04 ) and 1.272 × 10 −03 (1.245 × 10 −03 ) , respectively. These relative error curves are monotonically increasing and decreasing, as shown in Fig. 11. These values cannot be used to compare the sensitivity of skewness and kurtosis directly. The difference of the maximum and minimum values for skewness and kurtosis are 5.866 × 10 −08 and 2.663 × 10 −05 , respectively. From this, the effects of kurtosis on the estimation of skewness and vice versa are shown to be insignificant. Therefore, for this case, the estimation of skewness is independent of kurtosis and vice versa.

Convexity adjustment term errors
The errors of using the Black-Scholes return in place of BKM's approximation have been tabulated in Tables 4 and 5 and presented in Figs. 12 and 13.    Tables 4 and 5 have been combined with Tables 2 and 3 in Tables 6 and 7, respectively, to allow for ease of comparison.
In general, there is a slight decrease in both mean and absolute maximum errors for skewness when the boundary controlling factor is less than 0.75. For larger boundary controlling factors, the alternate convexity adjustment term increases the errors.  . The boundary controlling factor and step size are set to 0.25 and 1, respectively. When testing skewness, the true skewness is set to −1 . For kurtosis, the true value is set to 2.5 (Color figure online)

Table 4
Errors for alternative convexity adjustment term

Interpolation and extrapolation
The results of interpolation and extrapolation are shown in Tables 8 and 9, respectively. Table 8 shows that all interpolation techniques reduce the magnitude of skewness error for almost all points. There are a few interpolated points which have  1 3 errors with greater magnitudes than the non-interpolated case; however, this is due to the oscillating nature of the error. All three interpolation techniques yield much more stable results. Cubic spline performs the best, followed by linear interpolation, then Gaussian kernel regression. The linear interpolation option shows how

Table 6
Comparison of errors for the original and alternative convexity adjustment term   smoothing the implied volatility curve is more beneficial than interpolating with the option prices directly. As the trapezium rule is used for integration, this effectively is equivalent to linear interpolation of the option prices, which corresponds to the no interpolation or "None" case in Table 8. A reason for this is the shape of the implied volatility curve compared with the option price curve, with respect to the strike price or moneyness. Under the Black-Scholes model, the implied volatility curve is simply a constant line, whereas the price curve is much more complicated. Simply put, linear interpolation of the implied volatility curve is effectively interpolating the price curve non-linearly. Table 8 also shows that if cubic spline interpolation is used, the condition of having a minimum step size no greater than 0.1% of the forward price can be relaxed to 2.5% of the forward price. Table 9 shows that both constant and linear extrapolation, in general, decreases the errors. Constant extrapolation tends to yield more consistent results than linear extrapolation; however, linear extrapolation yields more accurate results. Linear extrapolation, although occasionally more accurate, can lead to negative implied volatilities, depending on the end points of the existing implied volatility curve. Constant extrapolation essentially assumes that the tails of the return are normally distributed, in the same way as Black-Scholes' model has constant volatility. If the BKM skewness estimator is used directly with no interpolation or extrapolation techniques on S&P 500 index options, based on Eq. (4), the errors are expected to exceed 10 −3 . However, this can be easily remediated by using cubic spline interpolation with linear interpolation. For more stable and reliable results, cubic spline interpolation paired with constant extrapolation may be more appropriate, depending on the breadth of strikes available. These results are also consistent with AZ-although, the conditions are more strict as the affine jumpdiffusion model is able to model more volatile markets.

In-depth error analysis
A closer inspection of errors shows that the current level of granularity of the step size ΔK = 1 is not sufficient to show the overall form of the relationship between the errors and step size. By reducing the step size to ΔK = 0.001 , this reveals a form which is closer to the true form. This is shown in Fig. 14.
This table shows the maximum values of the absolute errors of the standard deviation, skewness, and kurtosis estimators which were calculated from the valid Gram-Charlier region. The difference is that the expected return of BKM has been replaced by E Q t [R(t, )] = r − 1 2 2 + c . These values have been scaled. The options created have a forward price F T t , time to maturity , risk-free rate r, and standard deviation which have been arbitrarily set to $2000, one month, 2.4%, and 0.20, respectively. This has been done for various boundary controlling factor (a) values and step sizes ( ΔK ). The error for each a and ΔK is defined as Error ∶= max | Estimated Moment − True Moment | 1 3  risk-neutral moment… The error does not decrease exponentially with respect to the step size, but more so quadratically. The error envelope, positive (+) and negative (−) , was found using a least-squares approximation of a relatively simple parsimonious quadratic curve: The coefficients + 1 , + 2 , − 1 , and − 2 were found to be 2.984 × 10 −04 , 1.969 × 10 −05 , −1.639 × 10 −04 , and −8.959 × 10 −06 , respectively. Using these values, the asymmetric envelope can be decomposed into two components, a symmetric envelope and a trend term. The symmetric envelope has parameters s 1 and s 2 equal to 2.311 × 10 −04 and 1.433 × 10 −05 , respectively. The trend component has parameters t 1 and t 2 equal to 6.726 × 10 −05 and 5.366 × 10 −06 , respectively. The error is oscillating with an increasing period. To capture this behaviour, Eq. (10) seems to be suitable. The increasing period increases approximately quadratically with 1 and 2 equal to 6.980 × 10 −04 and 4.967 × 10 −04 , respectively. This model is also presented in Fig. 14. The methodology used to obtain these parameters is presented in the appendix, part G.
(10) 1 ΔK + 2 (ΔK) 2 . As the risk-neutral moments are composed of multiple contracts, which in itself produces errors, the errors of each contract with respect to the step size are shown in Fig. 15. 8 The errors of the volatility contract can be modelled the same way as the skewness estimator. The same error model is not appropriate for the other contracts. The coefficients for the volatility contract ( Fig. 14 Granularity of step size ΔK . The errors for skewness for a boundary controlling factor of a = 0.05 with varying step sizes between 1 and 50 shows that, in general, as the step size gets larger, so do the errors. When the fineness of the granularity is increased from ΔK = 1 (black line) to ΔK = 0.001 (red line), the shape of the relationship between the error and the step size is revealed. The forward price F T t , time to maturity , risk-free rate r, standard deviation , skewness 1 , and kurtosis 2 are equal to $2000, one month, 2.4%, 0.20, −1 , and 2.5, respectively. The error model and envelope ( E ± ) is given by Error Model = X s sin 2 X ΔK + X t and ΔK . The relative errors for the volatility, cubic, and quartic contracts for a boundary controlling factor of a = 0.05 with varying step sizes between 1 and 50 show that, in general, as the step size gets larger, so do the errors. The forward price F T t , time to maturity , risk-free rate r, standard deviation , skewness 1 , and kurtosis 2 are equal to $2000, one month, 2.4%, 0.20, −1 , and 2.5, respectively. The value of the volatility, cubic, and quartic contracts are 3.327 × 10 −03 , −1.884 × 10 −04 , and 6.071 × 10 −05 , respectively. The relative error is defined as Relative Error ∶= Estimated − True

Conclusion
The importance of the BKM method is clear. It is used by many practitioners and academics as an indicator of tail risk and also to calculate the CBOE SKEW. This paper demonstrates a method to create simulated option prices using the Gram-Charlier density. These simulated prices are created by specifying the mean, variance, skewness, and kurtosis. By doing so, the risk-neutral moment estimators can be compared with the true values. The results show that in order to estimate skewness within 10 −3 of the true value, the range of strikes ( K min , K max ) must contain at least 3/4 to 4/3 of the forward price and have a step size ( ΔK ) of no more than 0.1% of the forward price. Rather than using the absolute error as the measure, if the average error is used, then the step size restriction can be relaxed to 1% of the forward price. Under the same boundary controlling factor and step size specification, absolute errors of the kurtosis and standard deviation are bounded by 5 × 10 −3 and 10 −4 , respectively.
In general, markets are not able to satisfy the conditions required to calculate skewness with errors less than 10 −3 , hence, interpolation and extrapolation techniques are required. Cubic spline interpolation performs better than linear interpolation and Gaussian kernel regression. For extrapolation, there is a trade-off between stability and accuracy, with constant extrapolation being more stable than the more accurate linear extrapolation. Using both cubic spline and either constant or linear regressions will help reduce discretization and truncation errors of skewness.
AZ provide the prequel to this paper by examining the BKM risk-neutral skewness estimator under the affine jump-diffusion model. The results of this paper are consistent with AZ for market periods that are not too volatile. When volatile periods are included, AZ recommend that linear interpolation to 0.05% of the forward price and constant extrapolation to include strikes at least a quarter to quadruple the forward price ( a = 1∕4 ) should be performed on the implied volatility curve in order to bound the errors of the BKM skewness estimator to be less than 10 −3 . This condition is stricter than the a = 3∕4 and ΔK = 0.1%F T t found using the Gram-Charlier approach-this is expected. The Gram-Charlier density introduces skewness and excess kurtosis, whereas the Duffie et al. (2000) model has parameters, such as the instantaneous variance, v t , mean-reverting speed, , longterm mean, , diffusive volatility, , correlation between the stock and variance, , jump intensity, , mean jump size in variance, y , mean jump size in price, x + J y , variance of jump size in price, 2 x , and correlation between jumps, J , which can be combined to get the skewness. With the additional freedom given by the number of additional parameters and the model's specification, the affine jump-diffusion model can better model the market. However, due to the complexity of the Duffie et al. (2000) model, the errors of the BKM estimators were not analyzed as comprehensively or as in-depth by AZ.
Beyond AZ, this paper finds that the errors of skewness oscillate with respect to the step size. Increasing the granularity of the step size decreases the error approximately quadratically, not exponentially. However, this oscillation can results in H S t = 0 , H x S t = 0 , and Finally, the expected value of each contract is given by where n specifies the type of power contract and Q(K) corresponds to the OTM option with strike K. If there exists both put and call at the at-the-money point, then the average of the two is taken. n = {2, 3, 4} corresponds to the volatility (V), cubic (W), and quadratic (X) payoff contracts, respectively. With these payoff contracts, the skewness (and kurtosis) can be calculated using Eq. (A1).
The standardized skewness is given by Eq. (A1) when n = 3 . Expanding this, the BKM formula for the risk-neutral skewness can be found. Similarly, the standardized kurtosis can be found when n = 4 . BKM defines , V, W, and X as The (annualized) variance 2 is given by

Appendix B Truncation and discretization errors
The truncation and discretization errors are defined following AZ. For convenience, they are presented here. As data from the options market are not continuous, the integral calculation of Eq. (A6) must be solved numerically. Using the trapezium rule, the integral can be discretized to where and m is the number of strikes. This discretization introduces errors. This slightly differs from CBOE's discretization methods, as shown in the appendix, part C.

Truncation errors
= e r X − 4 e r W + 6 2 e r V − 4 3 + 4 e r V − 2 + 2 2 (A8) = e r X − 4 e r W + 6 2 e r V − 3 4 This proposition allows the variance swap to be calculated directly from the distribution of the returns, specifically, the moments and the time to maturity. For a stock price model with moments no higher than kurtosis, the VIX is given by A simpler form of Eq. (E30) can be obtained by approximating the log and square root, as shown here: By suppressing kurtosis, this formula shows that introducing negative skewness to returns will result in the VIX becoming smaller than the standard deviation. Originally, without kurtosis and skewness, the two were simply linked by a scaling factor.
The CBOE VIX is used as a benchmark against the VIX, calculated using the BKM method. This VIX can be calculated with The volatility index can also be calculated using the BKM method: As , the approximation of E Q t [R(t, )] , is required, this method introduces errors; however, it does remain model-free.

E.2 Rare disaster concern index
In this paper, the RIX is defined as to include both upside and downside jumps. Originally in Gao et al. (2018), the RIX only included downside jumps: The RIX is calculated from the difference of the variance of log returns ( ) and the variance swap ( ).

E.3 Simple variance swaps
From Martin (2016), the SVIX is given by where R T = S T S t and R t = e r . By setting q = 0 and using a more general case of Eq.

E.4 Errors of the additional variables
The VIX, RIX, and SVIX errors are shown in Table 10 with the errors of the BKM risk-neutral standard for comparison. From this table, the general trend of having a smaller error when the boundary controlling factor and step sizes are small holds; however, this trend is weaker for the RIX.

Appendix F BKM and Black-Scholes risk-neutral log returns
Using the standard Black-Scholes stock price model, 10 With the same stock price model: but with y following the extended Gram-Charlier density (to the mth additional cumulant, m ). The probability density function is given by where y now has a mean of 0, variance of 1, skewness of 1 , excess kurtosis of 2 , and so on ( i−2 ≈ i ).
c is found from making E t S T = F T t . Which gives c = (   0.94 − 6.23 − 6.22 − 6.22 − 6.19 − 6.20 − 6.12 − 6.12 − 6.11 − 6.10 − 5.21 0.95 − 6.84 − 6.83 − 6.82 − 6.82 − 6.84 − 6.74 − 6.74 − 6.25 − 6.69 − 6.69 This table shows the errors of the risk-neutral standard deviation, VIX, RIX, and SVIX estimators which were calculated for a single point within the valid Gram-Charlier region. The skewness ( 1 ) and kurtosis ( 2 ) is −1 and 2.5, respectively. These values have been scaled. The options created have a forward price F T t , time to maturity , risk-free rate r, and standard deviation which have been arbitrarily set to $2000, one month, 2.4%, and 0.20, respectively. This has been done for various boundary controlling factor (a) values and step sizes ( ΔK ). The error for each a and ΔK is defined as Error ∶= Estimated Moment − True Moment . The exact VIX, RIX, and SVIX values are 19.8136, 0.00074363, and 0.19479, respectively Funding Open Access funding enabled and organized by CAUL and its Member Institutions. Pakorn Aschakulporn appreciates being awarded the University of Otago Doctoral Scholarship. Jin E. Zhang has been supported by an establishment grant from the University of Otago. The authors have no financial or proprietary interests in any material discussed in this article.

Data availibility
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.